/* program for making initial coordinate and velocity */
/* by K.Saitoh  2002.5.15 */
/* compile: gcc with -lm option */ 

#include	<stdio.h>
#include	<math.h>
#include	<stdlib.h>	/* due to necessity of RAND_MAX */
#define	NM	10000		/* maximum array treatable */
#define KB	1.38066e-23	/* Boltzmann constant [J/K]*/
#define PI	3.14159265358979323e0
#define	SIGMA	2.0			/* standard deviation */
#define	Mass	(63.546e-03/6.022e23)	/* mass of a copper atom */

main()
{
	int	nx, ny, nz;
	double	a0;		/* lattice constant */
	int	i, j, k, natom;
	double	x[NM], y[NM], z[NM];
	double	vx[NM], vy[NM], vz[NM];
	double	x1, y1, z1;
	double	temp;
	
	fprintf(stderr, "(x periodicity)?\n", nx);
	scanf("%d", &nx);
	fprintf(stderr, "(y periodicity)?\n", ny);
	scanf("%d", &ny);
	fprintf(stderr, "(z periodicity)?\n", nz);
	scanf("%d", &nz);
	fprintf(stderr, "(lattice constant[m])?\n", a0);
	scanf("%lf", &a0);
	fprintf(stderr, "(temperature[K])?\n", temp);
	scanf("%lf", &temp);
	
/*** making initial coordinate */
	natom = 0;
	for (k = -nz; k < nz; k++) {
		z1 = (double)k * a0;
		for (j = -ny; j < ny; j++) {
			y1 = (double)j * a0;
			for (i = -nx; i < nx; i++) {
				x1 = (double)i * a0;
				x[natom] = x1;
				y[natom] = y1;
				z[natom] = z1;
				natom++;
			}
		}
		for (j = -ny; j < ny; j++) {
			y1 = (double)(j + 0.5) * a0;
			for (i = -nx; i < nx; i++) {
				x1 = (double)(i + 0.5) * a0;
				x[natom] = x1;
				y[natom] = y1;
				z[natom] = z1;
				natom++;
			}
		}
	}
	for (k = -nz; k < nz; k++) {
		z1 = (double)(k + 0.5) * a0;
		for (j = -ny; j < ny; j++) {
			y1 = (double)j * a0;
			for (i = -nx; i < nx; i++) {
				x1 = (double)(i + 0.5) * a0;
				x[natom] = x1;
				y[natom] = y1;
				z[natom] = z1;
				natom++;
			}
		}
		for (j = -ny; j < ny; j++) {
			y1 = (double)(j + 0.5) * a0;
			for (i = -nx; i < nx; i++) {
				x1 = (double)i * a0;
				x[natom] = x1;
				y[natom] = y1;
				z[natom] = z1;
				natom++;
			}
		}
	}
	fprintf(stderr, "(number of atoms)%d\n", natom);
/*** setting velocity using temperature ***/
	maxwell(natom, temp, &vx, &vy, &vz);

/*** data output ***/
	for (i = 0; i < natom; i++) 
		printf("%10d %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e\n", 
			i + 1, x[i], y[i], z[i], vx[i], vy[i], vz[i]);
}


/*** subroutine for initial velocity */
maxwell(natom, temp, vx1, vy1, vz1)
int	natom;
double	vx1[], vy1[], vz1[], temp;
{
	int	i, n;
	double	a, b, vmax, vmin, vx, vy, vz, ek, temp1;
	double	x1, x2, x3, x4, w;
	double	urand();

	fprintf(stderr, "(set temperature[K]) %12.5e\n", temp);

	b = 0.5 * Mass / (KB * temp);
	a = (double)natom * pow(b / PI, 1.5);

	vmax = SIGMA * sqrt(KB * temp / Mass);	/* maximum velocity */
	vmin = -vmax;				/* maximum velocity */

	n = 0; ek = 0.0;
	for (i = 0; i < natom; i++) {
		do {
			x1 = urand();	/* generation of random number (x)*/
			x2 = urand();	/* generation of random number (y)*/
			x3 = urand();	/* generation of random number (z)*/
			x4 = urand();	/* generation of random number */
			vx = vmin + (vmax - vmin) * x1;
			vy = vmin + (vmax - vmin) * x2;
			vz = vmin + (vmax - vmin) * x3;
			w = a * exp(-b * (vx * vx + vy * vy + vz * vz));
		} while (x4 * a > w);
		vx1[i] = vx;
		vy1[i] = vy;
		vz1[i] = vz;
		n++;
		ek += 0.5 * Mass * (vx * vx + vy * vy + vz * vz);
					/* ek : total kinetic energy */
	}
	fprintf(stderr, "(kinetic energy[J]) %12.5e\n", ek);
	temp1 = 2.0 * ek / (3.0 * natom * KB);		/* real temperature */
	fprintf(stderr, "(real temperature[K]) %12.5e\n", temp1);
}

/*** function for generating uniform random number ***/
double	urand()
{
	double	x;
	x = (double)rand() / (double)RAND_MAX;
	return(x);
}