#include	<stdio.h>
#include	<stdlib.h>
#include	<math.h>
/* (2013new) malloc variables are used, instead of using stack. */
#define	NM	300		/* Maximum number of atoms */
#define	PM	((NM * NM) / 2)	/* Maximum number of pairs */
/* */
#define	KB	1.38e-23		/* Boltzmann Constant */
#define	EPS	(100.0 * KB)		/* energy unit */
#define	SIGMA	2.0e-10			/* length unit*/
#define	TDELTA	5.0e-15			/* increment of time step */
#define	WM	(63.546e-03/6.022e23)	/* mass of atom (copper) */
#define	TSTEP	1000			/* the number of time step */
main()
{
	int	i, j, kstep, natom, npair, *ip, *jp;
	double	*x, *y, *z, *vx, *vy, *vz;
	double	*rx, *ry, *rz, *r, rs;
	double	*phi, *dphi, dphix, dphiy, dphiz;
	double	*fx, *fy, *fz;
	double	*x1, *y1, *z1, *vx1, *vy1, *vz1;
	double	ep, ek;
	double	time;
	double	tau, tstep;
	FILE	*fp;


/* normalization unit for time*/
	tau = SIGMA * sqrt(WM / EPS);
	tstep = TDELTA / tau;

	fp = fopen("energy.dat", "w");

	natom = 256;

/* malloc variables (i.e. dynamically allocating memory) */
	x = (double *)malloc( sizeof(double) * natom);
	y = (double *)malloc( sizeof(double) * natom);
	z = (double *)malloc( sizeof(double) * natom);
	vx = (double *)malloc( sizeof(double) * natom);
	vy = (double *)malloc( sizeof(double) * natom);
	vz = (double *)malloc( sizeof(double) * natom);
	fx = (double *)malloc( sizeof(double) * natom);
	fy = (double *)malloc( sizeof(double) * natom);
	fz = (double *)malloc( sizeof(double) * natom);
	x1 = (double *)malloc( sizeof(double) * natom);
	y1 = (double *)malloc( sizeof(double) * natom);
	z1 = (double *)malloc( sizeof(double) * natom);
	vx1 = (double *)malloc( sizeof(double) * natom);
	vy1 = (double *)malloc( sizeof(double) * natom);
	vz1 = (double *)malloc( sizeof(double) * natom);
	rx = (double *)malloc( sizeof(double) * (natom * natom / 2));
	ry = (double *)malloc( sizeof(double) * (natom * natom / 2));
	rz = (double *)malloc( sizeof(double) * (natom * natom / 2));
	r = (double *)malloc( sizeof(double) * (natom * natom / 2));
	phi = (double *)malloc( sizeof(double) * (natom * natom / 2));
	dphi = (double *)malloc( sizeof(double) * (natom * natom / 2));
	ip = (int *)malloc( sizeof(int) * (natom * natom / 2));
	jp = (int *)malloc( sizeof(int) * (natom * natom / 2));
/*	*/

	for (i = 0; i < natom; i++) {
		scanf("%*d%lf%lf%lf%lf%lf%lf", 
			&x[i], &y[i], &z[i], &vx[i], &vy[i], &vz[i]);
		x[i] /= SIGMA;
		y[i] /= SIGMA;
		z[i] /= SIGMA;
		vx[i] *= tau / SIGMA;
		vy[i] *= tau / SIGMA;
		vz[i] *= tau / SIGMA;
	}

	time = 0.0;	/* start time */
	for (kstep = 0; kstep < TSTEP; kstep++) {
	/* one step begins */
		fprintf(stderr, "(step)%5d ", kstep);

		npair = 0;
		for (i = 0; i < natom - 1; i++) {
			for (j = i + 1; j < natom; j++) {
				rx[npair] = x[i] - x[j];
				ry[npair] = y[i] - y[j];
				rz[npair] = z[i] - z[j];
				r[npair] = sqrt(rx[npair] * rx[npair] 
						+ ry[npair] * ry[npair]
						+ rz[npair] * rz[npair]);
				ip[npair] = i;
				jp[npair] = j;
				npair = npair + 1;
			}
		}
		fprintf(stderr, "(pair)%5d\n", npair);

		for (i = 0; i < natom; i++) {
			/* initialization of force vectors */
			fx[i] = 0.0; fy[i] = 0.0; fz[i] = 0.0;
		}
		ep = 0.0; ek = 0.0;	/* initialization of energies */
		for (i = 0; i < npair; i++) {
			rs = 1.0 / r[i];
			phi[i] = 4.0 * (pow(rs, 12.0) - pow(rs, 6.0));
			dphi[i] = -24.0 * (2.0 * pow(rs, 13.0) 
							- pow(rs, 7.0));
			dphix = -dphi[i] * rx[i] / r[i];
			dphiy = -dphi[i] * ry[i] / r[i];
			dphiz = -dphi[i] * rz[i] / r[i];
			fx[ip[i]] += dphix;
			fx[jp[i]] -= dphix;
			fy[ip[i]] += dphiy;
			fy[jp[i]] -= dphiy;
			fz[ip[i]] += dphiz;
			fz[jp[i]] -= dphiz;
			ep += phi[i];	/* summation of pot. energy */
		}
		for (i = 0; i < natom; i++) {
			ek += (vx[i] * vx[i] + vy[i] * vy[i] + vz[i] * vz[i]);
		}
		ek *= 0.5 * WM * SIGMA * SIGMA / tau / tau;
		ep *= EPS;

		/* output of energies */
		fprintf(fp, "%12.5e %12.5e %12.5e\n", time, ep, ek);

		for (i = 0; i < natom; i++) {
			vx1[i] = vx[i] + fx[i] * tstep;
			vy1[i] = vy[i] + fy[i] * tstep;
			vz1[i] = vz[i] + fz[i] * tstep;
			x1[i] = x[i] + vx1[i] * tstep;
			y1[i] = y[i] + vy1[i] * tstep;
			z1[i] = z[i] + vz1[i] * tstep;
		}
		time += TDELTA;
		for (i = 0; i < natom; i++) {
			printf("%2d %15.8e %15.8e %15.8e %15.8e %15.8e %15.8e\n", 
				i, x1[i] * SIGMA, y1[i] * SIGMA, z1[i] * SIGMA, 
				vx1[i] * SIGMA / tau, vy1[i] * SIGMA / tau, 
				vz1[i] * SIGMA / tau);
		}
		for (i = 0; i < natom; i++) {
			x[i] = x1[i]; y[i] = y1[i]; z[i] = z1[i];
			vx[i] = vx1[i]; vy[i] = vy1[i]; vz[i] = vz1[i];
		}
	}
	/* one step ends */
	fclose(fp);

/* freeing for malloc */
	free(x); free(y); free(z);
	free(vx); free(vy); free(vz);
	free(fx); free(fy); free(fz);
	free(x1); free(y1); free(z1);
	free(vx1); free(vy1); free(vz1);
	free(rx); free(ry); free(rz); free(r);
	free(phi); free(dphi);
	free(ip); free(jp);

}