#include <stdio.h>
	#include <stdlib.h>
	#include <time.h>
	#include <math.h>
	#include "atoms.h"
	#include "header.h"
	#define XMUL 2
	#define YMUL 1
	#define ZMUL 4
	#define N 100
	#define TEMP 350.0
	const double sqrt2 = 1.4142135623730950488016887242097;
	const double sqrt3 = 1.7320508075688772935274463415059;
	const double kB = 1.380658E-23;
	const double RMo2 = (double)RAND_MAX / 2.0;
	int main()
	{
	Headerinfo tempheader;
	Atom* allatoms;
	int currentatom;
	int numatoms;
	FILE* file;
	int i,j,k,l,n;
	double currentT;
	double vMul;
	double netmx,netmy,netmz;
	double a;
	file = fopen("infile.atm", "wb");
	if(file == NULL)
	{
	printf("Failed to open infile.atm for writing\n\n");
	exit(6);
	}
	numatoms = (2+6*N)*XMUL*YMUL*ZMUL;
	allatoms = (Atom*)malloc(numatoms*sizeof(Atom));
	if(allatoms == NULL)
	{
	printf("unable to allocate memory for allatoms\n\n");
	fclose(file);
	exit(5);
	}
	/*srand((unsigned int)time(NULL));*/
	srand(2197);
	a = (elements[(enum AtomTypes)H].radius + elements[(enum AtomTypes)C].radius) * 2.0 / sqrt3;
	netmx = netmy = netmz = 0.0;
	currentT = 0.0;
	tempheader.iterations = 1000000;
	tempheader.mincutoff = 10.0 * a;
	tempheader.dt = 1.0E-15;
	tempheader.bigX = a + 6.4 * (double)XMUL * a + a*(double)N;
	tempheader.bigY = a + 8.4 * (double)YMUL * a + a*(double)N;
	tempheader.bigZ = a + 7.8 * (double)ZMUL * a + a*(double)(N/4);
	tempheader.T = TEMP;
	tempheader.desiredT = TEMP;
	tempheader.numatoms = numatoms;
	currentatom = 0;
	for(i=0; i<XMUL; i++)
	{
	for(j=0; j<YMUL; j++)
	{
	for(k=0; k<ZMUL; k++)
	{
	allatoms[currentatom].type = C;
	allatoms[currentatom].x = 2.0 * a + 6.4 * (double)i * a;
	allatoms[currentatom].y = 2.0 * a + 8.4 * (double)j * a;
	allatoms[currentatom].z = 2.0 * a + 7.8 * (double)k * a + 0.25 * a * (double)(N/4);
	allatoms[currentatom].vx = ((double)rand()-RMo2)/RMo2;
	allatoms[currentatom].vy = ((double)rand()-RMo2)/RMo2;
	allatoms[currentatom].vz = ((double)rand()-RMo2)/RMo2;
	allatoms[currentatom].ax = 0.0;
	allatoms[currentatom].ay = 0.0;
	allatoms[currentatom].az = 0.0;
	allatoms[currentatom].ndx = 0.0;
	allatoms[currentatom].ndy = 0.0;
	allatoms[currentatom].ndz = 0.0;
	allatoms[currentatom].potential = 0.0;
	allatoms[currentatom].newpotential = 0.0;
	for(l=0; l<NUMBONDS; l++) allatoms[currentatom].bonds[l] = 0;
	allatoms[currentatom].bonds[0] = 1;
	allatoms[currentatom].bonds[1] = 2;
	allatoms[currentatom].bonds[2] = 3;
	allatoms[currentatom].bonds[3] = 4;
	currentT += 0.5 * elements[allatoms[currentatom].type].mass * (allatoms[currentatom].vx*allatoms[currentatom].vx + allatoms[currentatom].vy*allatoms[currentatom].vy + allatoms[currentatom].vz*allatoms[currentatom].vz);
	netmx += elements[allatoms[currentatom].type].mass * allatoms[currentatom].vx;
	netmy += elements[allatoms[currentatom].type].mass * allatoms[currentatom].vy;
	netmz += elements[allatoms[currentatom].type].mass * allatoms[currentatom].vz;
	currentatom++;
	allatoms[currentatom].type = H;
	allatoms[currentatom].x = 1.5 * a + 6.4 * (double)i * a;
	allatoms[currentatom].y = 1.5 * a + 8.4 * (double)j * a;
	allatoms[currentatom].z = 1.5 * a + 7.8 * (double)k * a + 0.25 * a * (double)(N/4);
	allatoms[currentatom].vx = ((double)rand()-RMo2)/RMo2;
	allatoms[currentatom].vy = ((double)rand()-RMo2)/RMo2;
	allatoms[currentatom].vz = ((double)rand()-RMo2)/RMo2;
	allatoms[currentatom].ax = 0.0;
	allatoms[currentatom].ay = 0.0;
	allatoms[currentatom].az = 0.0;
	allatoms[currentatom].ndx = 0.0;
	allatoms[currentatom].ndy = 0.0;
	allatoms[currentatom].ndz = 0.0;
	allatoms[currentatom].potential = 0.0;
	allatoms[currentatom].newpotential = 0.0;
	for(l=0; l<NUMBONDS; l++) allatoms[currentatom].bonds[l] = 0;
	allatoms[currentatom].bonds[0] = -1;
	currentT += 0.5 * elements[allatoms[currentatom].type].mass * (allatoms[currentatom].vx*allatoms[currentatom].vx + allatoms[currentatom].vy*allatoms[currentatom].vy + allatoms[currentatom].vz*allatoms[currentatom].vz);
	netmx += elements[allatoms[currentatom].type].mass * allatoms[currentatom].vx;
	netmy += elements[allatoms[currentatom].type].mass * allatoms[currentatom].vy;
	netmz += elements[allatoms[currentatom].type].mass * allatoms[currentatom].vz;
	currentatom++;
	allatoms[currentatom].type = H;
	allatoms[currentatom].x = 1.5 * a + 6.4 * (double)i * a;
	allatoms[currentatom].y = 2.5 * a + 8.4 * (double)j * a;
	allatoms[currentatom].z = 2.5 * a + 7.8 * (double)k * a + 0.25 * a * (double)(N/4);
	allatoms[currentatom].vx = ((double)rand()-RMo2)/RMo2;
	allatoms[currentatom].vy = ((double)rand()-RMo2)/RMo2;
	allatoms[currentatom].vz = ((double)rand()-RMo2)/RMo2;
	allatoms[currentatom].ax = 0.0;
	allatoms[currentatom].ay = 0.0;
	allatoms[currentatom].az = 0.0;
	allatoms[currentatom].ndx = 0.0;
	allatoms[currentatom].ndy = 0.0;
	allatoms[currentatom].ndz = 0.0;
	allatoms[currentatom].potential = 0.0;
	allatoms[currentatom].newpotential = 0.0;
	for(l=0; l<NUMBONDS; l++) allatoms[currentatom].bonds[l] = 0;
	allatoms[currentatom].bonds[0] = -2;
	currentT += 0.5 * elements[allatoms[currentatom].type].mass * (allatoms[currentatom].vx*allatoms[currentatom].vx + allatoms[currentatom].vy*allatoms[currentatom].vy + allatoms[currentatom].vz*allatoms[currentatom].vz);
	netmx += elements[allatoms[currentatom].type].mass * allatoms[currentatom].vx;
	netmy += elements[allatoms[currentatom].type].mass * allatoms[currentatom].vy;
	netmz += elements[allatoms[currentatom].type].mass * allatoms[currentatom].vz;
	currentatom++;
	allatoms[currentatom].type = H;
	allatoms[currentatom].x = 2.5 * a + 6.4 * (double)i * a;
	allatoms[currentatom].y = 1.5 * a + 8.4 * (double)j * a;
	allatoms[currentatom].z = 2.5 * a + 7.8 * (double)k * a + 0.25 * a * (double)(N/4);
	allatoms[currentatom].vx = ((double)rand()-RMo2)/RMo2;
	allatoms[currentatom].vy = ((double)rand()-RMo2)/RMo2;
	allatoms[currentatom].vz = ((double)rand()-RMo2)/RMo2;
	allatoms[currentatom].ax = 0.0;
	allatoms[currentatom].ay = 0.0;
	allatoms[currentatom].az = 0.0;
	allatoms[currentatom].ndx = 0.0;
	allatoms[currentatom].ndy = 0.0;
	allatoms[currentatom].ndz = 0.0;
	allatoms[currentatom].potential = 0.0;
	allatoms[currentatom].newpotential = 0.0;
	for(l=0; l<NUMBONDS; l++) allatoms[currentatom].bonds[l] = 0;
	allatoms[currentatom].bonds[0] = -3;
	currentT += 0.5 * elements[allatoms[currentatom].type].mass * (allatoms[currentatom].vx*allatoms[currentatom].vx + allatoms[currentatom].vy*allatoms[currentatom].vy + allatoms[currentatom].vz*allatoms[currentatom].vz);
	netmx += elements[allatoms[currentatom].type].mass * allatoms[currentatom].vx;
	netmy += elements[allatoms[currentatom].type].mass * allatoms[currentatom].vy;
	netmz += elements[allatoms[currentatom].type].mass * allatoms[currentatom].vz;
	currentatom++;
	for(n=0; n<N-1; n++)
	{
	allatoms[currentatom].type = C;
	allatoms[currentatom].x = 2.5 * a + 6.4 * (double)i * a + a*n;
	allatoms[currentatom].y = 2.5 * a + 8.4 * (double)j * a + a*n;
	allatoms[currentatom].z = 1.5 * a + 7.8 * (double)k * a + 0.25 * a * (double)(N/4) + 0.00*a*n;
	allatoms[currentatom].vx = ((double)rand()-RMo2)/RMo2;
	allatoms[currentatom].vy = ((double)rand()-RMo2)/RMo2;
	allatoms[currentatom].vz = ((double)rand()-RMo2)/RMo2;
	allatoms[currentatom].ax = 0.0;
	allatoms[currentatom].ay = 0.0;
	allatoms[currentatom].az = 0.0;
	allatoms[currentatom].ndx = 0.0;
	allatoms[currentatom].ndy = 0.0;
	allatoms[currentatom].ndz = 0.0;
	allatoms[currentatom].potential = 0.0;
	allatoms[currentatom].newpotential = 0.0;
	for(l=0; l<NUMBONDS; l++) allatoms[currentatom].bonds[l] = 0;
	allatoms[currentatom].bonds[0] = -4;
	allatoms[currentatom].bonds[1] = 1;
	allatoms[currentatom].bonds[2] = 2;
	allatoms[currentatom].bonds[3] = 3;
	currentT += 0.5 * elements[allatoms[currentatom].type].mass * (allatoms[currentatom].vx*allatoms[currentatom].vx + allatoms[currentatom].vy*allatoms[currentatom].vy + allatoms[currentatom].vz*allatoms[currentatom].vz);
	netmx += elements[allatoms[currentatom].type].mass * allatoms[currentatom].vx;
	netmy += elements[allatoms[currentatom].type].mass * allatoms[currentatom].vy;
	netmz += elements[allatoms[currentatom].type].mass * allatoms[currentatom].vz;
	currentatom++;
	allatoms[currentatom].type = H;
	allatoms[currentatom].x = 2.0 * a + 6.4 * (double)i * a + a*n;
	allatoms[currentatom].y = 3.0 * a + 8.4 * (double)j * a + a*n;
	allatoms[currentatom].z = 1.0 * a + 7.8 * (double)k * a + 0.25 * a * (double)(N/4) + 0.00*a*n;
	allatoms[currentatom].vx = ((double)rand()-RMo2)/RMo2;
	allatoms[currentatom].vy = ((double)rand()-RMo2)/RMo2;
	allatoms[currentatom].vz = ((double)rand()-RMo2)/RMo2;
	allatoms[currentatom].ax = 0.0;
	allatoms[currentatom].ay = 0.0;
	allatoms[currentatom].az = 0.0;
	allatoms[currentatom].ndx = 0.0;
	allatoms[currentatom].ndy = 0.0;
	allatoms[currentatom].ndz = 0.0;
	allatoms[currentatom].potential = 0.0;
	allatoms[currentatom].newpotential = 0.0;
	for(l=0; l<NUMBONDS; l++) allatoms[currentatom].bonds[l] = 0;
	allatoms[currentatom].bonds[0] = -1;
	currentT += 0.5 * elements[allatoms[currentatom].type].mass * (allatoms[currentatom].vx*allatoms[currentatom].vx + allatoms[currentatom].vy*allatoms[currentatom].vy + allatoms[currentatom].vz*allatoms[currentatom].vz);
	netmx += elements[allatoms[currentatom].type].mass * allatoms[currentatom].vx;
	netmy += elements[allatoms[currentatom].type].mass * allatoms[currentatom].vy;
	netmz += elements[allatoms[currentatom].type].mass * allatoms[currentatom].vz;
	currentatom++;
	allatoms[currentatom].type = H;
	allatoms[currentatom].x = 3.0 * a + 6.4 * (double)i * a + a*n;
	allatoms[currentatom].y = 2.0 * a + 8.4 * (double)j * a + a*n;
	allatoms[currentatom].z = 1.0 * a + 7.8 * (double)k * a + 0.25 * a * (double)(N/4) + 0.00*a*n;
	allatoms[currentatom].vx = ((double)rand()-RMo2)/RMo2;
	allatoms[currentatom].vy = ((double)rand()-RMo2)/RMo2;
	allatoms[currentatom].vz = ((double)rand()-RMo2)/RMo2;
	allatoms[currentatom].ax = 0.0;
	allatoms[currentatom].ay = 0.0;
	allatoms[currentatom].az = 0.0;
	allatoms[currentatom].ndx = 0.0;
	allatoms[currentatom].ndy = 0.0;
	allatoms[currentatom].ndz = 0.0;
	allatoms[currentatom].potential = 0.0;
	allatoms[currentatom].newpotential = 0.0;
	for(l=0; l<NUMBONDS; l++) allatoms[currentatom].bonds[l] = 0;
	allatoms[currentatom].bonds[0] = -2;
	currentT += 0.5 * elements[allatoms[currentatom].type].mass * (allatoms[currentatom].vx*allatoms[currentatom].vx + allatoms[currentatom].vy*allatoms[currentatom].vy + allatoms[currentatom].vz*allatoms[currentatom].vz);
	netmx += elements[allatoms[currentatom].type].mass * allatoms[currentatom].vx;
	netmy += elements[allatoms[currentatom].type].mass * allatoms[currentatom].vy;
	netmz += elements[allatoms[currentatom].type].mass * allatoms[currentatom].vz;
	currentatom++;
	allatoms[currentatom].type = C;
	allatoms[currentatom].x = 2.0 * a + 6.4 * (double)i * a + a*(n+1);
	allatoms[currentatom].y = 2.0 * a + 8.4 * (double)j * a + a*(n+1);
	allatoms[currentatom].z = 2.0 * a + 7.8 * (double)k * a + 0.25 * a * (double)(N/4) + 0.00*a*n;
	allatoms[currentatom].vx = ((double)rand()-RMo2)/RMo2;
	allatoms[currentatom].vy = ((double)rand()-RMo2)/RMo2;
	allatoms[currentatom].vz = ((double)rand()-RMo2)/RMo2;
	allatoms[currentatom].ax = 0.0;
	allatoms[currentatom].ay = 0.0;
	allatoms[currentatom].az = 0.0;
	allatoms[currentatom].ndx = 0.0;
	allatoms[currentatom].ndy = 0.0;
	allatoms[currentatom].ndz = 0.0;
	allatoms[currentatom].potential = 0.0;
	allatoms[currentatom].newpotential = 0.0;
	for(l=0; l<NUMBONDS; l++) allatoms[currentatom].bonds[l] = 0;
	allatoms[currentatom].bonds[0] = -3;
	allatoms[currentatom].bonds[1] = 1;
	allatoms[currentatom].bonds[2] = 2;
	allatoms[currentatom].bonds[3] = 3;
	currentT += 0.5 * elements[allatoms[currentatom].type].mass * (allatoms[currentatom].vx*allatoms[currentatom].vx + allatoms[currentatom].vy*allatoms[currentatom].vy + allatoms[currentatom].vz*allatoms[currentatom].vz);
	netmx += elements[allatoms[currentatom].type].mass * allatoms[currentatom].vx;
	netmy += elements[allatoms[currentatom].type].mass * allatoms[currentatom].vy;
	netmz += elements[allatoms[currentatom].type].mass * allatoms[currentatom].vz;
	currentatom++;
	allatoms[currentatom].type = H;
	allatoms[currentatom].x = 1.5 * a + 6.4 * (double)i * a + a*(n+1);
	allatoms[currentatom].y = 2.5 * a + 8.4 * (double)j * a + a*(n+1);
	allatoms[currentatom].z = 2.5 * a + 7.8 * (double)k * a + 0.25 * a * (double)(N/4) + 0.00*a*n;
	allatoms[currentatom].vx = ((double)rand()-RMo2)/RMo2;
	allatoms[currentatom].vy = ((double)rand()-RMo2)/RMo2;
	allatoms[currentatom].vz = ((double)rand()-RMo2)/RMo2;
	allatoms[currentatom].ax = 0.0;
	allatoms[currentatom].ay = 0.0;
	allatoms[currentatom].az = 0.0;
	allatoms[currentatom].ndx = 0.0;
	allatoms[currentatom].ndy = 0.0;
	allatoms[currentatom].ndz = 0.0;
	allatoms[currentatom].potential = 0.0;
	allatoms[currentatom].newpotential = 0.0;
	for(l=0; l<NUMBONDS; l++) allatoms[currentatom].bonds[l] = 0;
	allatoms[currentatom].bonds[0] = -1;
	currentT += 0.5 * elements[allatoms[currentatom].type].mass * (allatoms[currentatom].vx*allatoms[currentatom].vx + allatoms[currentatom].vy*allatoms[currentatom].vy + allatoms[currentatom].vz*allatoms[currentatom].vz);
	netmx += elements[allatoms[currentatom].type].mass * allatoms[currentatom].vx;
	netmy += elements[allatoms[currentatom].type].mass * allatoms[currentatom].vy;
	netmz += elements[allatoms[currentatom].type].mass * allatoms[currentatom].vz;
	currentatom++;
	allatoms[currentatom].type = H;
	allatoms[currentatom].x = 2.5 * a + 6.4 * (double)i * a + a*(n+1);
	allatoms[currentatom].y = 1.5 * a + 8.4 * (double)j * a + a*(n+1);
	allatoms[currentatom].z = 2.5 * a + 7.8 * (double)k * a + 0.25 * a * (double)(N/4) + 0.00*a*n;
	allatoms[currentatom].vx = ((double)rand()-RMo2)/RMo2;
	allatoms[currentatom].vy = ((double)rand()-RMo2)/RMo2;
	allatoms[currentatom].vz = ((double)rand()-RMo2)/RMo2;
	allatoms[currentatom].ax = 0.0;
	allatoms[currentatom].ay = 0.0;
	allatoms[currentatom].az = 0.0;
	allatoms[currentatom].ndx = 0.0;
	allatoms[currentatom].ndy = 0.0;
	allatoms[currentatom].ndz = 0.0;
	allatoms[currentatom].potential = 0.0;
	allatoms[currentatom].newpotential = 0.0;
	for(l=0; l<NUMBONDS; l++) allatoms[currentatom].bonds[l] = 0;
	allatoms[currentatom].bonds[0] = -2;
	currentT += 0.5 * elements[allatoms[currentatom].type].mass * (allatoms[currentatom].vx*allatoms[currentatom].vx + allatoms[currentatom].vy*allatoms[currentatom].vy + allatoms[currentatom].vz*allatoms[currentatom].vz);
	netmx += elements[allatoms[currentatom].type].mass * allatoms[currentatom].vx;
	netmy += elements[allatoms[currentatom].type].mass * allatoms[currentatom].vy;
	netmz += elements[allatoms[currentatom].type].mass * allatoms[currentatom].vz;
	currentatom++;
	}
	allatoms[currentatom].type = C;
	allatoms[currentatom].x = 2.5 * a + 6.4 * (double)i * a + a*(N-1);
	allatoms[currentatom].y = 2.5 * a + 8.4 * (double)j * a + a*(N-1);
	allatoms[currentatom].z = 1.5 * a + 7.8 * (double)k * a + 0.25 * a * (double)(N/4) + 0.00*a*(N-1);
	allatoms[currentatom].vx = ((double)rand()-RMo2)/RMo2;
	allatoms[currentatom].vy = ((double)rand()-RMo2)/RMo2;
	allatoms[currentatom].vz = ((double)rand()-RMo2)/RMo2;
	allatoms[currentatom].ax = 0.0;
	allatoms[currentatom].ay = 0.0;
	allatoms[currentatom].az = 0.0;
	allatoms[currentatom].ndx = 0.0;
	allatoms[currentatom].ndy = 0.0;
	allatoms[currentatom].ndz = 0.0;
	allatoms[currentatom].potential = 0.0;
	allatoms[currentatom].newpotential = 0.0;
	for(l=0; l<NUMBONDS; l++) allatoms[currentatom].bonds[l] = 0;
	allatoms[currentatom].bonds[0] = -4;
	allatoms[currentatom].bonds[1] = 1;
	allatoms[currentatom].bonds[2] = 2;
	allatoms[currentatom].bonds[3] = 3;
	currentT += 0.5 * elements[allatoms[currentatom].type].mass * (allatoms[currentatom].vx*allatoms[currentatom].vx + allatoms[currentatom].vy*allatoms[currentatom].vy + allatoms[currentatom].vz*allatoms[currentatom].vz);
	netmx += elements[allatoms[currentatom].type].mass * allatoms[currentatom].vx;
	netmy += elements[allatoms[currentatom].type].mass * allatoms[currentatom].vy;
	netmz += elements[allatoms[currentatom].type].mass * allatoms[currentatom].vz;
	currentatom++;
	allatoms[currentatom].type = H;
	allatoms[currentatom].x = 2.0 * a + 6.4 * (double)i * a + a*(N-1);
	allatoms[currentatom].y = 3.0 * a + 8.4 * (double)j * a + a*(N-1);
	allatoms[currentatom].z = 1.0 * a + 7.8 * (double)k * a + 0.25 * a * (double)(N/4) + 0.00*a*(N-1);
	allatoms[currentatom].vx = ((double)rand()-RMo2)/RMo2;
	allatoms[currentatom].vy = ((double)rand()-RMo2)/RMo2;
	allatoms[currentatom].vz = ((double)rand()-RMo2)/RMo2;
	allatoms[currentatom].ax = 0.0;
	allatoms[currentatom].ay = 0.0;
	allatoms[currentatom].az = 0.0;
	allatoms[currentatom].ndx = 0.0;
	allatoms[currentatom].ndy = 0.0;
	allatoms[currentatom].ndz = 0.0;
	allatoms[currentatom].potential = 0.0;
	allatoms[currentatom].newpotential = 0.0;
	for(l=0; l<NUMBONDS; l++) allatoms[currentatom].bonds[l] = 0;
	allatoms[currentatom].bonds[0] = -1;
	currentT += 0.5 * elements[allatoms[currentatom].type].mass * (allatoms[currentatom].vx*allatoms[currentatom].vx + allatoms[currentatom].vy*allatoms[currentatom].vy + allatoms[currentatom].vz*allatoms[currentatom].vz);
	netmx += elements[allatoms[currentatom].type].mass * allatoms[currentatom].vx;
	netmy += elements[allatoms[currentatom].type].mass * allatoms[currentatom].vy;
	netmz += elements[allatoms[currentatom].type].mass * allatoms[currentatom].vz;
	currentatom++;
	allatoms[currentatom].type = H;
	allatoms[currentatom].x = 3.0 * a + 6.4 * (double)i * a + a*(N-1);
	allatoms[currentatom].y = 2.0 * a + 8.4 * (double)j * a + a*(N-1);
	allatoms[currentatom].z = 1.0 * a + 7.8 * (double)k * a + 0.25 * a * (double)(N/4) + 0.00*a*(N-1);
	allatoms[currentatom].vx = ((double)rand()-RMo2)/RMo2;
	allatoms[currentatom].vy = ((double)rand()-RMo2)/RMo2;
	allatoms[currentatom].vz = ((double)rand()-RMo2)/RMo2;
	allatoms[currentatom].ax = 0.0;
	allatoms[currentatom].ay = 0.0;
	allatoms[currentatom].az = 0.0;
	allatoms[currentatom].ndx = 0.0;
	allatoms[currentatom].ndy = 0.0;
	allatoms[currentatom].ndz = 0.0;
	allatoms[currentatom].potential = 0.0;
	allatoms[currentatom].newpotential = 0.0;
	for(l=0; l<NUMBONDS; l++) allatoms[currentatom].bonds[l] = 0;
	allatoms[currentatom].bonds[0] = -2;
	currentT += 0.5 * elements[allatoms[currentatom].type].mass * (allatoms[currentatom].vx*allatoms[currentatom].vx + allatoms[currentatom].vy*allatoms[currentatom].vy + allatoms[currentatom].vz*allatoms[currentatom].vz);
	netmx += elements[allatoms[currentatom].type].mass * allatoms[currentatom].vx;
	netmy += elements[allatoms[currentatom].type].mass * allatoms[currentatom].vy;
	netmz += elements[allatoms[currentatom].type].mass * allatoms[currentatom].vz;
	currentatom++;
	allatoms[currentatom].type = H;
	allatoms[currentatom].x = 3.0 * a + 6.4 * (double)i * a + a*(N-1);
	allatoms[currentatom].y = 3.0 * a + 8.4 * (double)j * a + a*(N-1);
	allatoms[currentatom].z = 2.0 * a + 7.8 * (double)k * a + 0.25 * a * (double)(N/4) + 0.00*a*(N-1);
	allatoms[currentatom].vx = ((double)rand()-RMo2)/RMo2;
	allatoms[currentatom].vy = ((double)rand()-RMo2)/RMo2;
	allatoms[currentatom].vz = ((double)rand()-RMo2)/RMo2;
	allatoms[currentatom].ax = 0.0;
	allatoms[currentatom].ay = 0.0;
	allatoms[currentatom].az = 0.0;
	allatoms[currentatom].ndx = 0.0;
	allatoms[currentatom].ndy = 0.0;
	allatoms[currentatom].ndz = 0.0;
	allatoms[currentatom].potential = 0.0;
	allatoms[currentatom].newpotential = 0.0;
	for(l=0; l<NUMBONDS; l++) allatoms[currentatom].bonds[l] = 0;
	allatoms[currentatom].bonds[0] = -3;
	currentT += 0.5 * elements[allatoms[currentatom].type].mass * (allatoms[currentatom].vx*allatoms[currentatom].vx + allatoms[currentatom].vy*allatoms[currentatom].vy + allatoms[currentatom].vz*allatoms[currentatom].vz);
	netmx += elements[allatoms[currentatom].type].mass * allatoms[currentatom].vx;
	netmy += elements[allatoms[currentatom].type].mass * allatoms[currentatom].vy;
	netmz += elements[allatoms[currentatom].type].mass * allatoms[currentatom].vz;
	currentatom++;
	}
	}
	}
	currentT /= (double)numatoms;
	currentT /= 1.5 * kB;
	netmx /= (double)numatoms;
	netmy /= (double)numatoms;
	netmz /= (double)numatoms;
	vMul = sqrt(TEMP/currentT);
	for(i=0; i<numatoms; i++)
	{
	while(allatoms[i].x > tempheader.bigX) allatoms[i].x -= tempheader.bigX;
	while(allatoms[i].x < 0.0) allatoms[i].x += tempheader.bigX;
	while(allatoms[i].y > tempheader.bigY) allatoms[i].y -= tempheader.bigY;
	while(allatoms[i].y < 0.0) allatoms[i].y += tempheader.bigY;
	while(allatoms[i].z > tempheader.bigZ) allatoms[i].z -= tempheader.bigZ;
	while(allatoms[i].z < 0.0) allatoms[i].z += tempheader.bigZ;
	allatoms[i].vx = (allatoms[i].vx - netmx/elements[allatoms[i].type].mass) * vMul;
	allatoms[i].vy = (allatoms[i].vy - netmy/elements[allatoms[i].type].mass) * vMul;
	allatoms[i].vz = (allatoms[i].vz - netmz/elements[allatoms[i].type].mass) * vMul;
	}
	fwrite(&tempheader, sizeof(Headerinfo), 1, file);
	fwrite(allatoms, sizeof(Atom), numatoms, file);
	free(allatoms);
	fclose(file);
	return 1;
	}