1. #include <stdio.h>  
2. #include <math.h>  
3. #include <stdlib.h>  
4. #include <time.h>  
5. #include <string.h>  
6. #define nDim 3  
7. #define sigma 3.4e-10  
8. #define epsilon 1.65e-21  
9. #define masskg 6.6e-26  
10.   
11. double random1();  
12. double maxwell();  
13. double gauss(double std, double mean );  
14. void velscaler(double *velocity, double T, int n);  
15. void initialization(double *positions, double *velocities, int n, double L, double T);  
16. int neighborhood(double *position, int *celllist, int n, double L, double rcut);  
17. void LJ_Forces( double *force, double *potential, double *position, int n, int Nmax, int *celllist, double L, double rcut, double rmax, double *px);  
18. void andersen(int switch1, int n, double *force, double *velocity, double *positions,double *tracer, double dt, double T, double L, int *blinker );  
19. void properties(double *position, double *tracer, double *velocities, double *forces, double *potential, int n, double *U, double *KE, double *energy, double *avevel, double *Temp, double *px, double *pressure, double v, double rcut, double rho, double L, double *Ds, double timecount);  
20. void position_printer(FILE * fp, double * positions, int n, int *blinker, double L, int flag);  
21. void velocity_printer(FILE *fp, double * velocity, int n );  
22.   
23. main() {  
24.     int n, i, j, ncel, *link, *celllist, *neighborlist, nMax, neighborcounter, cycles, track, *blinker;  
25.     double Rho, T, dt, dt_squared, L, vol, mass, ndouble, *positions, *velocities, rcut, *force;  
26.     double rcel, *potential, kb, kT, *acceleration, Lcell, rmax, timecounter, error, v, U, totale;  
27.     double *energies, KE, *avevel, Temp, *forcematrix, pressure, *pix, *tracer, Ds, tau, *forcestore;  
28.     double squarevel[3], aveke, avepe, avetotal, aveT, aveP, aveDs;   
29.     double andersentime,veltime, codetime,productiontime, inittime;  
30.     char buffer[50], vyes[20];  
31.     char yes[20];  
32.     FILE * fp;  
33.     FILE * tracker;  
34.     FILE * melter;  
35.     FILE * velocityfp;  
36.     FILE * OUT;  
37.     clock_t codebegin, codeend, initbegin, initend, scalebegin, scaleend, andersenbegin, andersenend;   
38.     clock_t productionbegin, productionend;  
39.     codebegin = clock();  
40.     srandom(time(NULL));  
41.     printf("How many cycles?\n");  
42.     scanf("%d",&cycles);  
43.     //cycles = 3000;  
44.     printf("How many atoms (perfect cube)?\n");  
45.     scanf("%d",&n);  
46.     //n = 3375;  
47.     rcut = 2.5; //2.5sigma  
48.     rmax = rcut + 1;   
49.     printf("Rho (Lennard-Jones units)?\n");  
50.     scanf("%lf",&Rho);  
51.     //Rho = 0.90399;   
52.     printf("Temp (Lennard-Jones units)?\n");  
53.     scanf("%lf",&T);  
54.     //T = 1.17145;  
55.     mass = 1;  
56.     tau = sigma*sqrt(masskg/epsilon);  
57.     dt = 1e-15/tau;  
58.     L = pow((double)n / Rho, 0.333333);  
59.     vol = n/Rho;  
60.     //printf("%lf %lf\n", L, vol);  
61.     positions = (double*)malloc(nDim*n*sizeof(double));  
62.     velocities = (double*)malloc(nDim*n*sizeof(double));   
63.     force = (double*)malloc(nDim*n*sizeof(double));  
64.     potential = (double*)malloc(n*sizeof(double));  
65.     neighborlist = (int*)malloc(n*n*sizeof(int));  
66.     avevel = (double*)malloc(nDim*sizeof(double));  
67.     pix = (double*)malloc(nDim*sizeof(double));  
68.     tracer = (double*)malloc(nDim*n*sizeof(double));  
69.     blinker = (int*)malloc(nDim*n*sizeof(int));  
70.     memset(positions,0.0,nDim*n*sizeof(double));  
71.     memset(velocities,0.0,nDim*n*sizeof(double));  
72.     memset(force,0.0,nDim*n*sizeof(double));  
73.     memset(potential,0.0,n*sizeof(double));  
74.     memset(neighborlist,0,n*n*sizeof(int));  
75.     memset(avevel,0.0,3*sizeof(double));  
76.     memset(tracer,0.0,nDim*n*sizeof(double));  
77.     memset(tracer,0.0,nDim*n*sizeof(double));  
78.     memset(blinker,0,nDim*n*sizeof(int));  
79.     OUT = fopen("output","w");  
80.     printf("Would you like to apply an Andersen Thermostat? [Y/N]\n");  
81.     scanf("%s",yes);  
82.     printf("Would you like to Velocity Scale? [Y/N]\n");  
83.     scanf("%s",vyes);  
84.     printf("What particle would you like to track?\n");  
85.     scanf("%d",&track);  
86.     initbegin = clock();  
87.     fp = fopen("Initial.xyz","w");  
88.     velocityfp = fopen("Initialvel","w");  
89.     initialization(positions,velocities,n,L,T);  
90.     velocity_printer(velocityfp,velocities,n);  
91.     position_printer(fp,positions,n,blinker,L,0);  
92.     fclose(fp);  
93.     fclose(velocityfp);  
94.     melter = fopen("Melting.xyz","w");  
95.     position_printer(melter,positions,n, blinker, L,1);  
96.     nMax = neighborhood(positions,neighborlist,n,L,rmax);   
97.     if (nMax == 0 ) {  
98.         return 0;  
99.     }  
100.     LJ_Forces(force,potential,positions,n,nMax,neighborlist,L,rcut,rmax, pix);  
101.     initend = clock();  
102.     inittime = (double)(initend - initbegin) / CLOCKS_PER_SEC;  
103.     timecounter = 0.0;  
104.     neighborcounter = 0;  
105.     if (strcmp(vyes,"Y") == 0) {   
106.     scalebegin = clock();  
107.     printf("\n\nMelting at constant Temp\n\n");  
108.     fprintf(OUT,"\n\nMelting at constant Temp\n\n");  
109.         for(i=0;i<cycles;i++) {  
110.             neighborcounter += 1;  
111.             andersen(1,n,force,velocities,positions,tracer, dt, T, L, blinker);  
112.             LJ_Forces(force,potential,positions,n,nMax,neighborlist,L,rcut,rmax, pix);  
113.             andersen(2,n,force,velocities,positions,tracer,dt, T, L, blinker);  
114.             if (neighborcounter > 10) {  
115.                 nMax = neighborhood(positions,neighborlist,n,L,rmax);  
116.                 neighborcounter = 0;  
117.                 velscaler(velocities,T,n);  
118.                 position_printer(melter,positions,n, blinker, L,1);  
119.             }  
120.             timecounter += dt;  
121.             properties(positions,tracer,velocities,force,potential,n,&U,&KE,&totale, avevel, &Temp,pix, &pressure, vol, rcut, Rho, L, &Ds, timecounter);  
122.             printf("%lf %lf %lf %lf %lf %lf %lf\n", timecounter, KE, U, totale, pressure, Ds, Temp);  
123.             fprintf(OUT,"%lf %lf %lf %lf %lf %lf %lf\n", timecounter, KE, U, totale, pressure, Ds, Temp);  
124.         }   
125.         velocityfp = fopen("Meltedvel","w");  
126.         velocity_printer(velocityfp,velocities,n);  
127.         fclose(velocityfp);  
128.         sprintf(buffer,"%lf.xyz",dt*i);  
129.         fp = fopen("melted.xyz","w");  
130.         position_printer(fp,positions,n,blinker,L,0);  
131.         fclose(fp);  
132.     scaleend = clock();  
133.     veltime = (double)(scaleend - scalebegin) / CLOCKS_PER_SEC;  
134.     }  
135.     if( strcmp(yes,"Y") == 0 ) {  
136.     andersenbegin = clock();  
137.         printf("\n\nAndersen Thermostat\n\n");  
138.         fprintf(OUT,"\n\nAndersen Thermostat\n\n");  
139.         velocityfp = fopen("Andersenvel","w");  
140.         neighborcounter = 0;   
141.         for(i=0;i<cycles;i++) {  
142.             neighborcounter += 1;  
143.             andersen(1,n,force,velocities,positions,tracer, dt, T, L, blinker);  
144.             LJ_Forces(force,potential,positions,n,nMax,neighborlist,L,rcut,rmax, pix);  
145.             andersen(3,n,force,velocities,positions,tracer,dt, T, L, blinker);  
146.             if ((i % 10) == 0) {  
147.                 nMax = neighborhood(positions,neighborlist,n,L,rmax);  
148.                 neighborcounter = 0;  
149.                 position_printer(melter,positions,n, blinker, L,1);  
150.             }  
151.             timecounter += dt;  
152.             properties(positions,tracer,velocities,force,potential,n,&U,&KE,&totale, avevel, &Temp,pix, &pressure, vol, rcut, Rho, L, &Ds, timecounter);  
153.             printf("%lf %lf %lf %lf %lf %lf %lf\n", timecounter, KE, U, totale, pressure, Ds, Temp);  
154.             fprintf(OUT,"%lf %lf %lf %lf %lf %lf %lf\n", timecounter, KE, U, totale, pressure, Ds, Temp);  
155.               
156.         }   
157.         velocity_printer(velocityfp,velocities,n);  
158.         fclose(velocityfp);  
159.     andersenend = clock();  
160.     andersentime = (double)(andersenend - andersenbegin) / CLOCKS_PER_SEC;   
161.     }  
162.       
163.     productionbegin = clock();  
164.     memset(tracer,0.0,nDim*n*sizeof(double));  
165.     cycles = 10000;  
166.     tracker = fopen("tracking","w");  
167.     fprintf(tracker,"0 0 0\n");  
168.     neighborcounter = 0;  
169.     aveke=avepe=avetotal=aveT=aveP=aveDs=0.0;  
170.     printf("\n\nNVE Production\n\n");  
171.     fprintf(OUT,"\n\nNVE Production\n\n");  
172.     velocityfp = fopen("FinalVel","w");  
173.     for(i=0;i<cycles;i++) {  
174.         timecounter += dt;  
175.         neighborcounter += 1;  
176.         andersen(1,n,force,velocities,positions,tracer,dt, T, L, blinker);  
177.         fprintf(tracker,"%lf %lf %lf\n", *(tracer + track), *(tracer + n + track), *(tracer + 2*n + track));  
178.         LJ_Forces(force,potential,positions,n,nMax,neighborlist,L,rcut,rmax,pix);  
179.         andersen(2,n,force,velocities,positions,tracer,dt, T, L, blinker);  
180.         if ((i % 10) == 0) {  
181.             nMax = neighborhood(positions,neighborlist,n,L,rmax);  
182.             neighborcounter = 0;  
183.             position_printer(melter,positions,n, blinker, L,1);  
184.         }  
185.         timecounter += dt;  
186.         properties(positions,tracer,velocities,force,potential,n,&U,&KE,&totale, avevel, &Temp,pix, &pressure, vol, rcut, Rho, L, &Ds, timecounter);  
187.         printf("%lf %lf %lf %lf %lf %lf %lf\n", timecounter, KE, U, totale, pressure, Ds, Temp);  
188.         fprintf(OUT,"%lf %lf %lf %lf %lf %lf %lf\n", timecounter, KE, U, totale, pressure, Ds, Temp);  
189.         aveke += KE;  
190.         avepe += U;  
191.         avetotal += totale;  
192.         aveP += pressure;  
193.         aveT += Temp;  
194.         aveDs += Ds;  
195.     }   
196.     productionend = clock();  
197.     productiontime = (double)(productionend - productionbegin) / CLOCKS_PER_SEC;  
198.     velocity_printer(velocityfp,velocities,n);  
199.     fclose(velocityfp);  
200.     fclose(tracker);  
201.     fclose(melter);  
202.     aveke /= cycles;  
203.     avepe /= cycles;  
204.     avetotal /=cycles;  
205.     aveP /= cycles;  
206.     aveT /= cycles;  
207.     aveDs /= cycles;  
208.       
209.     Ds = Ds / (cycles*dt);  
210.       
211.     printf("<ke>\t<u>\t<e>\tD\t<P>\t<T>\t\n");  
212.     printf("%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t\n",aveke,avepe,avetotal,Ds,aveP,aveT);  
213.     fprintf(OUT,"<ke>\t<u>\t<e>\tD\t<P>\t<T>\t\n");  
214.     fprintf(OUT,"%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t\n",aveke,avepe,avetotal,Ds,aveP,aveT);    
215.       
216.     fclose(OUT);  
217.   
218.     free(positions);  
219.     free(velocities);  
220.     free(force);  
221.     free(potential);  
222.     free(neighborlist);  
223.     free(avevel);   
224.     free(pix);  
225.     free(tracer);   
226.     free(blinker);   
227.     codeend = clock();  
228.     codetime = (double)(codeend - codebegin) / CLOCKS_PER_SEC;   
229.       
230.     printf("Timing results\n");  
231.     if (strcmp(vyes,"Y") == 0) {   
232.         printf("Whole Code\tInitialization\tVelscale\tProduction\n");  
233.         printf("%lf\t%lf\t%lf\t%lf\n",codetime,inittime,veltime,productiontime);  
234.       
235.     } else if( strcmp(yes,"Y") == 0 ) {  
236.         printf("Whole Code\tInitialization\tAndersen\tProduction\n");  
237.         printf("%lf\t%lf\t%lf\t%lf\n",codetime,inittime,andersentime,productiontime);  
238.     }  
239.       
240. }  
241.   
242. double random1() {  
243.     return (double)random()/RAND_MAX;  
244. }  
245.   
246. double maxwell() { //johnk's random number (from a maxwellian distribution) generator  
247. //referenced from: http://www.stat.rice.edu/~dobelman/textfiles/DistributionsHandbook.pdf  
248.     double x1, x2, w, y1, y2;  
249.     do {  
250.         x1 = 2.0 * random1() - 1.0;  
251.         x2 = 2.0 * random1() - 1.0;  
252.         w = x1 * x1 + x2 * x2;  
253.     } while ( w >= 1.0 );  
254.     w = sqrt( (-2.0 * log10( w ) ) / w );  
255.     y1 = x1 * w;  
256.     y2 = x2 * w;  
257.     return y2;  
258.   
259. }  
260.   
261. double gauss(double std, double mean ) {  
262. //Used by Andersen (properly scaled velocities)  
263.     double rand;   
264.     rand = random1();  
265.       
266.     rand = mean + std*rand;  
267.     return rand;  
268. }  
269.   
270. void velscaler(double *velocity, double T, int n) {  
271. //Scales the velocity to proper temp, T.  
272.     int i,j;  
273.     double square_vel[nDim], velscale;  
274.     for (i = 0; i<nDim; i++) {  
275.         square_vel[i] = 0.0;  
276.     }  
277.       
278.     for (i=0; i<n; i++) {  
279.         for (j=0;j<nDim;j++) {  
280.             square_vel[j] +=  *(velocity + j*n+ i) * *(velocity +j*n + i);  
281.         }  
282.     }  
283.     velscale = sqrt(3*T/((square_vel[0]+square_vel[1]+square_vel[2])/n));  
284.     for (i=0;i<n;i++) {  
285.         for (j=0;j<nDim;j++) {  
286.             *(velocity + j*n + i) *= velscale;  
287.         }  
288.     }  
289.   
290. }  
291.   
292. void initialization(double *positions, double *velocities, int n, double L, double T) {  
293.     int N, i, j, k, flag;  
294.     double del, ndouble, per_dimension, *sum_velocity, *squared_velocity, velscale, rcutoff;  
295.     double r, tester;  
296.     per_dimension = pow(n,1/3.0);  
297.     rcutoff = 0.212; //Argon diameter divided by sigma  
298.     del = L/per_dimension;  
299.     //printf("%d %lf %lf\n", n, per_dimension, del);  
300.     N = 0;  
301.     flag = 0;  
302.     sum_velocity = malloc(3*sizeof(double));  
303.     squared_velocity = malloc(3*sizeof(double));  
304.     for(i=0;i<nDim;i++) {  
305.         *(sum_velocity + i) = 0.0;  
306.         *(squared_velocity + i) = 0.0;  
307.     }  
308.   
309.     for (i=0;i< per_dimension; i++) {  
310.         for(j=0;j< per_dimension; j++) {  
311.             for(k=0; k< per_dimension;k++) {  
312.                 *(positions + 0*n + N) = i * del;  
313.                 *(positions + 1*n + N) = j * del;  
314.                 *(positions + 2*n + N) = k * del;  
315.                 N++;  
316.             }  
317.         }  
318.     }  
319.       
320.       
321.     for(i=0; i<n; i++) {  
322.         *(velocities + i) = maxwell();  
323.         *(velocities + 1*n + i) = maxwell();  
324.         *(velocities + 2*n + i) = maxwell();  
325.         //below for calculating center of mass velocity  
326.         *(sum_velocity) += *(velocities + i)/n;   
327.         *(sum_velocity + 1) += *(velocities + 1*n + i)/n;  
328.         *(sum_velocity + 2) += *(velocities + 2*n + i)/n;  
329.     }  
330.     //correcting for CM drift  
331.     for (i=0;i<n;i++) {  
332.         *(velocities + i) -= *(sum_velocity);  
333.         *(velocities + 1*n + i) -= *(sum_velocity + 1);  
334.         *(velocities + 2*n + i) -= *(sum_velocity + 2);  
335.     }  
336.     velscaler(velocities,T,n);  
337. }  
338.   
339. int neighborhood(double *position, int *celllist, int n, double L, double rcut) {  
340. //Neighborhood list generation  
341.     int i,j,k,xcell,ycell,zcell, Nmax;  
342.     double x,y,z, halfL, distance[nDim], r;  
343.     Nmax = 0;   
344.     for (i = 0; i<n; i++) {  
345.         for(j=i+1;j<n; j++) {  
346.             r = 0.0;  
347.             if(i == j) { break;}  
348.             distance[0] = *(position + i) - *(position + j);  
349.             distance[1] = *(position + n + i) - *(position + n + j);  
350.             distance[2] = *(position + 2*n + i) - *(position + 2*n + j);  
351.             //periodic boundary condition:  
352.             halfL = L/2.0;  
353.             for (k=0;k<nDim;k++) {  
354.                 if (distance[k] > halfL) {   
355.                     distance[k] = distance[k] - L;   
356.                 } else if( distance[k] < -halfL ) {  
357.                     distance[k] = distance[k] + L;  
358.                 }  
359.                 r = r + distance[k]*distance[k];  
360.             }  
361.             r = sqrt(r);   
362.             // generates a Nx2 matrix with the list of atoms   
363.             //each row gives i,j corresponding to atoms near eachother  
364.             if (r <= rcut) {  
365.                 *(celllist + Nmax*2 + 1) = i;  
366.                 *(celllist + Nmax*2 + 2) = j;  
367.                 Nmax = Nmax + 1;   
368.             }  
369.         }  
370.     }  
371.     if (Nmax == 0 ) {  
372.         printf("No near atoms\n");  
373.         return 0;   
374.     }  
375.     return Nmax;  
376. }  
377.   
378. void LJ_Forces( double *force, double *potential, double *position, int n, int Nmax, int *celllist, double L, double rcut, double rmax, double *px){  
379.     int i, j, k, x, y, z, num;   
380.     double distance[nDim], halfL,ir6, U, F, r, r6, ir, r2, numerator, den, rmax2, rcut2, virial, tester;   
381.     double ucut;  
382.     ucut = 4*(1/(rcut*rcut*rcut*rcut*rcut*rcut))*((1/(rcut*rcut*rcut*rcut*rcut*rcut)) - 1);  
383.     virial = 0.0;  
384.     memset(potential,0.0,n*sizeof(double));  
385.     memset(force,0.0,nDim*n*sizeof(double));  
386.     memset(px,0.0,nDim*sizeof(double));  
387.   
388.     virial = 0.0;  
389.     halfL = L/2.0;  
390.     for ( num = 0; num<Nmax; num++) {  
391.             r = 0.0;  
392.             i = *(celllist + num*2 + 1);  
393.             j = *(celllist + num*2 + 2);  
394.             if( i != j) {    
395.                 distance[0] = *(position + i) - *(position + j);  
396.                 distance[1] = *(position + n + i) - *(position + n + j);  
397.                 distance[2] = *(position + 2*n + i) - *(position + 2*n + j);  
398.                 for (k=0;k<nDim;k++) {  
399.                     if (distance[k] > halfL) {   
400.                         distance[k] = distance[k] - L;   
401.                     }else if( distance[k] < -halfL ) {  
402.                         distance[k] = distance[k] + L;  
403.                     }  
404.                 }  
405.                 r2 = distance[0]*distance[0]+distance[1]*distance[1]+distance[2]*distance[2];  
406.                 r = sqrt(r2);   
407.                 //printf("%lf\n",r);  
408.                 if (r <= rcut) {  
409.                     ir = 1.0/r;  
410.                     r6 = r2*r2*r2;  
411.                     ir6 = 1.0/r6;   
412.                     U = 4.0*ir6*(ir6-1.);   
413.                     F = 48.0*(ir6*ir6-0.5*ir6);  
414.                     //printf("%lf %lf %lf\n", r, U, F);  
415.                     *(potential + i) += U;  
416.                     for (k =0; k<nDim; k++) {  
417.                         *(force + k*n + i) += F*distance[k]*ir*ir;  
418.                         *(force + k*n + j) -= F*distance[k]*ir*ir; //equal and opposite forces  
419.                     }  
420.                     virial += F;  
421.                   
422.                 }else if(r>rcut && r< rmax) {  
423.                     r6 = r*r*r*r*r*r;  
424.                     r2 = r*r;  
425.                     rmax2 = rmax*rmax;  
426.                     rcut2 = rcut*rcut;  
427.                     numerator = 24.0*(rmax*rmax-r2)*(rmax2*rmax2*(r6-2.0)-rmax2*(r6-2.0)*(3.0*rcut2-r2)+rcut2*r2*(r6-4.0)+2.0*r2*r2);  
428.                     den = r6*r6*r*(rmax2-rcut2)*(rmax2-rcut2)*(rmax2-rcut2);  
429.                     ir6 = 1./r6;   
430.                     U = 4.0*ir6*(ir6-1.0)*(rmax*rmax-r*r)*(rmax*rmax-r*r)*(rmax*rmax+2.0*r*r-3.0*rcut*rcut)/((rmax*rmax-rcut*rcut)*(rmax*rmax-rcut*rcut)*(rmax*rmax-rcut*rcut));   
431.                     *(potential + i) += U;  
432.                     F = -numerator/den;  
433.                     for (k =0; k<nDim; k++) {  
434.                         *(force + k*n + i) += F*distance[k]/r;  
435.                         *(force + k*n + j) -= F*distance[k]/r; //equal and opposite forces  
436.                     }  
437.                 }  
438.               
439.             }         
440.     }  
441.     *(px) = virial;  
442. }  
443.   
444. void andersen(int switch1, int n, double *force, double *velocity, double *positions,double *tracer, double dt, double T, double L, int *blinker ) {  
445. //Andersen thermostat  
446.     int i,j;  
447.     double Temp, rand;   
448.     double nu = 0.2; //collision with heat bath (needs to be a controllable variable)  
449.     if (switch1 == 1) {  
450.         for (i=0;i<n;i++) {  
451.             for (j = 0; j<nDim; j++) {  
452.                 *(positions + j*n + i ) += *(velocity +j*n + i)*dt + 0.5 * *(force +j*n + i)*dt*dt;  
453.                 *(velocity + j*n + i ) += 0.5 * (*(force + j*n + i)) * dt;  
454.                 *(tracer + j*n + i) += *(velocity +j*n + i)*dt + 0.5 * *(force +j*n + i)*dt*dt;  
455.                 if( *(positions + j*n + i) < 0.0) {   
456.                     *(positions + j*n + i) = *(positions + j*n + i) + L;   
457.                     *(blinker + j*n + i) -= 1;  
458.                 }  
459.                 if( *(positions + j*n + i) > L ) {   
460.                     *(positions + j*n + i) = *(positions + j*n + i) - L;   
461.                     *(blinker + j*n + i) += 1;  
462.                 }     
463.             }  
464.         }  
465.     } else if(switch1 == 3 || switch1 == 2 ) {  
466.         Temp = 0.0; //initialize temp  
467.         for (i = 0;i<n;i++) {   
468.             for (j=0;j<nDim;j++) {  
469.             *(velocity + j*n + i ) = *(velocity + j*n + i ) + 0.5 * *(force + j*n + i) * dt;  
470.             }  
471.         }  
472.         if(switch1 == 3) {  
473.             for(i=0;i<n;i++) {  
474.                 for(j=0;j<nDim;j++) {  
475.                     rand = random1();  
476.                     if(rand < nu*dt) {  
477.                         //printf("%lf hit\n",rand);  
478.                         *(velocity + j*n + i) = gauss(sqrt(T),0);   
479.                     }  
480.                 }  
481.             }  
482.         }  
483.     }   
484. }  
485.   
486. void properties(double *position, double *tracer, double *velocities, double *forces, double *potential, int n, double *U, double *KE, double *energy, double *avevel, double *Temp, double *px, double *pressure, double v, double rcut, double rho, double L, double *Ds, double timecount) {  
487.     int i,j,k;  
488.     double fx, fy, fz, radius, radiust, radiusz, virial;  
489.     double distance[3], F, radius6, radius12, ucor, pcor, rcut9, rcut3, halfL;  
490.     rcut3 = rcut*rcut*rcut;  
491.     rcut9 = rcut3*rcut3*rcut3;  
492.     ucor = 8.0*3.1415926*rho*(1.0/(9.0*rcut9)-1.0/(3.0*rcut3));  
493.     pcor = 16.0*3.1415926*rho*(2.0/(9.0*rcut9)-1.0/(3.0*rcut3));  
494.     halfL = L/2.0;  
495.     *U = 0.0;  
496.     *KE = 0.0;   
497.     *(avevel) = 0.0;  
498.     *(avevel + 1 ) = 0.0;  
499.     *(avevel + 2) = 0.0;  
500.     radiust = 0.0;  
501.     *Ds = 0.0;  
502.     for (j=0;j<n;j++) {  
503.         *(avevel) += *(velocities + j);  
504.         *(avevel + 1 ) += *(velocities + n + j);  
505.         *(avevel + 2) += *(velocities + 2*n + j);  
506.         *U += *(potential + j);  
507.         *Ds += *(tracer + j) * *(tracer + j) + *(tracer + n + j) * *(tracer + n +j) + *(tracer + 2*n + j) * *(tracer + 2*n +j);  
508.           
509.     }  
510.     *Ds = *Ds / (2*nDim*n) ;  
511.     // Below is just a sanity check to make sure that we have no "Flying Ice"  
512.     *(avevel) /= n;  
513.     *(avevel + 1 ) /= n;  
514.     *(avevel + 2) /= n;  
515.     //Above should all be identically 0.0, andersen thermostat may negate this  
516.     for(j = 0; j<n; j++) {  
517.         *KE += (*(velocities + j) - *(avevel)) * (*(velocities + j) - *(avevel)) + (*(velocities + n + j)- *(avevel + 1)) * (*(velocities + n + j)- *(avevel + 1)) + (*(velocities + 2*n + j)-*(avevel + 2)) * (*(velocities + 2*n + j)-*(avevel + 2));  
518.     }  
519.     *U /= n;  
520.     *U += ucor;  
521.     *KE = *KE/(2.0 * n);  
522.     *Temp = (2.0/3.0) * *KE ;   
523.     *energy = *KE + *U;  
524.     *pressure = 0.0;  
525.     *pressure = *px / 3.0;  
526.     *pressure += n * *Temp;  
527.     *pressure /= v;  
528.     *pressure += pcor;   
529. }  
530.   
531. void position_printer(FILE * fp, double * positions, int n, int *blinker, double L, int flag) {  
532.     int i;  
533.     double length;  
534.     length = sigma*10.0/(1.0e-9); //converting distance to angstrom (for animation)  
535.     fprintf(fp,"%d\n\n",n);  
536.     if (flag == 1) { //animation  
537.         for (i=0;i<n;i++) {  
538.             fprintf(fp,"Ar %lf %lf %lf\n", *(positions + i)*length + *(blinker + i)*L*length, *(positions + n + i)*length + *(blinker + n + i)*L*length, *(positions + 2*n + i)*length + *(blinker + 2*n + i)*L*length);  
539.         }  
540.     } else if(flag == 0) {  
541.         for (i=0;i<n;i++) {  
542.             fprintf(fp,"Ar %lf %lf %lf\n", *(positions + i)*length, *(positions + n + i)*length, *(positions + 2*n + i)*length);  
543.         }  
544.     }  
545. }  
546.   
547. void velocity_printer(FILE *fp, double * velocity, int n ) {  
548.     int i;  
549.     double Vx, Vy, Vz;  
550.     double V;  
551.     for(i=0;i<n;i++) {  
552.         Vx = *(velocity + i) * *(velocity + i);  
553.         Vy = *(velocity + n + i) * *(velocity + n + i);  
554.         Vz = *(velocity + 2*n + i) * *(velocity + 2*n + i);  
555.         V = sqrt(Vx + Vy + Vz);  
556.         fprintf(fp,"%lf\n",V);  
557.     }  
558. }