ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/branches/mmeineke/madProps/chrisProps.cpp
Revision: 38
Committed: Fri Jul 19 01:37:38 2002 UTC (23 years, 5 months ago) by mmeineke
File size: 27754 byte(s)
Log Message: 
A half polished version of props

File Contents

# User Rev Content
1 mmeineke 38 // Bloated Analysis Program
2     // Shamefully written by Chris Fennell
3     // Updated: 4/17/02
4     // Rewritten: 5/15/02
5    
6     #include <iostream>
7     #include <fstream>
8     #include <cmath>
9     #include <cstdio>
10     #include <cstdlib>
11     #include <string>
12     #include <cstring>
13     #include <vector>
14    
15     using namespace std;
16    
17     class Props {
18     public:
19     Props(string fileName, const char *splitter, vector<int> &pT, int delR);
20     ~Props();
21     void DoRMSDCorrelation();
22     void DoVelCorrelation();
23     void DoGofR();
24     void DoDipoleCorr();
25     void DoOrientCorr();
26     void DoMagFileCorr();
27     void DoCosineCorr();
28     private:
29     string filer;
30     const char *delimiter,*file;
31     char *token;
32     char inLine[1000],inValue[200], allList[200], velList[200];
33     char smallList[200], largeList[200],gofrList[200],sepFiler[200];
34     char atomNeon[3],atomX[3],atomO[3],atomH[3],atomSSD[4];
35     int nAtoms, i, j, k, l, simTime, simStep, counter, coorCount, count;
36     int blackHole, m, start, end, atomer, countSm, countBg, fineGrain;
37     int bin, option, outTyper, counter1, counter2, counter3, counter4, counter5;
38     // int *pAtomIdent;
39     vector<int> pAtomIdent;
40     double result, rX1, rY1, rZ1, rX2, rY2, rZ2, vX1, vY1, vZ1, vX2, vY2, vZ2;
41     double uX1, uY1, uZ1, uX2, uY2, uZ2, velDot, velDotNaught, dotDivide;
42     double mSdist, sum, average, xsum, ysum, xxsum, xysum, boxLength;
43     double densityStar, rijDist, rXij, rYij, rZij, deltaR, constant;
44     double rLower, rUpper, nIdeal, rho, rhoA, rhoB, rhoO, rhoH, rhoX;
45     double constantA, constantB, constantX, constantO, constantH, nIdealA;
46     double nIdealB, nIdealO, nIdealH, valueOuts, earlyTime, lateTime;
47     double earlyAvg, lateAvg, totalAvg, timeDiffer, crossTime;
48     double boxx, boxy, boxz, dipoleXo, dipoleYo, dipoleZo;
49     double dipoleXt, dipoleYt, dipoleZt, dipODotdipO, dipTDotdipO;
50     };
51    
52     Props::Props(string fileName, const char *splitter, vector<int> &pT, int delR){
53     filer = fileName;
54     file = filer.c_str();
55     delimiter = splitter;
56     pAtomIdent = pT;
57     fineGrain = delR;
58    
59    
60     strcpy(atomNeon,"Ne");
61     strcpy(atomX,"X");
62     strcpy(atomO,"O");
63     strcpy(atomH,"H");
64     strcpy(atomSSD,"SSD");
65    
66     ifstream inputer(file);
67    
68     inputer.getline(inLine,999,'\n');
69     token = strtok(inLine, delimiter);
70     strcpy(inValue,token);
71     nAtoms = atoi(inValue);
72    
73     counter = 0; counter1 = 0; counter2 = 0; counter3 = 0; counter4 = 0;
74     counter5 = 0;
75    
76     // Identify the particles we're dealing with
77     inputer.getline(inLine,999,'\n');
78     for (i = 0; i < nAtoms; i++) {
79     inputer.getline(inLine,999,'\n');
80     token = strtok(inLine, delimiter);
81     strcpy(inValue,token);
82     if (strcmp(atomNeon,inValue) == 0) {
83     pAtomIdent.push_back(0);
84     counter++;
85     }
86     else if (strcmp(atomX,inValue) == 0) {
87     pAtomIdent.push_back(2);
88     counter2++;
89     }
90     else if (strcmp(atomO,inValue) == 0) {
91     pAtomIdent.push_back(3);
92     counter3++;
93     }
94     else if (strcmp(atomH,inValue) == 0) {
95     pAtomIdent.push_back(4);
96     counter4++;
97     }
98     else if (strcmp(atomSSD, inValue) == 0) {
99     pAtomIdent.push_back(5);
100     counter5++;
101     }
102     else {
103     pAtomIdent.push_back(1);
104     counter1++;
105     }
106     }
107    
108     inputer.clear();
109     }
110    
111     Props::~Props() {
112     }
113    
114    
115     void Props::DoRMSDCorrelation() {
116     strcpy(allList,file);
117     strcat(allList,"rcorr");
118    
119     counter = nAtoms;
120    
121     if (counter5 > 0)
122     counter = nAtoms;
123     else {
124     nAtoms = nAtoms/4; //A quick hack for SSD water (4 atoms in each SSD)
125     counter = nAtoms;
126     }
127    
128     vector<double> time;
129     vector<double> pxSpot;
130     vector<double> pySpot;
131     vector<double> pzSpot;
132     vector<double> avg;
133    
134     cout << "Reading in Trajectory...\n";
135    
136     // Fill the vectors with the needed information from the .dump or .xyz file
137     ifstream inputer(file);
138     coorCount = 0;
139    
140     if (counter5 > 0) {
141     inputer.getline(inLine,999,'\n');
142     while (!inputer.eof()) {
143     inputer.getline(inLine,999,'\n');
144     token = strtok(inLine, delimiter);
145     strcpy(inValue,token);
146     time.push_back(atof(inValue));
147     for (i = 0; i < nAtoms; i++) {
148     inputer.getline(inLine,999,'\n');
149     token = strtok(inLine, delimiter);
150     token = strtok(NULL, delimiter);
151     strcpy(inValue,token);
152     pxSpot.push_back(atof(inValue));
153     token = strtok(NULL, delimiter);
154     strcpy(inValue,token);
155     pySpot.push_back(atof(inValue));
156     token = strtok(NULL, delimiter);
157     strcpy(inValue,token);
158     pzSpot.push_back(atof(inValue));
159     }
160     coorCount++;
161     inputer.getline(inLine,999,'\n');
162     }
163     }
164     else {
165     inputer.getline(inLine,999,'\n');
166     while (!inputer.eof()) {
167     inputer.getline(inLine,999,'\n');
168     token = strtok(inLine, delimiter);
169     strcpy(inValue,token);
170     time.push_back(atof(inValue));
171     for (i = 0; i < nAtoms; i++) {
172     inputer.getline(inLine,999,'\n');
173     token = strtok(inLine, delimiter);
174     token = strtok(NULL, delimiter);
175     strcpy(inValue,token);
176     pxSpot.push_back(atof(inValue));
177     token = strtok(NULL, delimiter);
178     strcpy(inValue,token);
179     pySpot.push_back(atof(inValue));
180     token = strtok(NULL, delimiter);
181     strcpy(inValue,token);
182     pzSpot.push_back(atof(inValue));
183     inputer.getline(inLine,999,'\n');
184     inputer.getline(inLine,999,'\n');
185     inputer.getline(inLine,999,'\n');
186     }
187     coorCount++;
188     inputer.getline(inLine,999,'\n');
189     }
190     }
191     cout << "Performing RMSD Calculations...\n";
192    
193     // The main calculation loop...
194     for (j = 0; j < coorCount-1; j++) {
195     sum = 0;
196     blackHole = 0;
197     for (i = 0; (i + j) < coorCount; i++) {
198     for (k = 0; k < counter; k++) {
199     rX1 = pxSpot[k + (i * counter)];
200     rY1 = pySpot[k + (i * counter)];
201     rZ1 = pzSpot[k + (i * counter)];
202     rX2 = pxSpot[k + (i * counter) + (j * counter)];
203     rY2 = pySpot[k + (i * counter) + (j * counter)];
204     rZ2 = pzSpot[k + (i * counter) + (j * counter)];
205    
206     mSdist = (pow((rX2-rX1),2) + pow((rY2-rY1),2)
207     + pow((rZ2-rZ1),2));
208    
209     sum += mSdist;
210     }
211     blackHole++;
212     }
213     average = sum / counter;
214     average /= blackHole;
215     avg.push_back(average);
216     }
217    
218     // Write out the results...
219     ofstream deltaOut(allList);
220     deltaOut << "#dt\tRMSD\n";
221     for (i = 0; i < coorCount-1; i++)
222     deltaOut << time[i]-time[0] << "\t" << avg[i] << "\n";
223    
224     inputer.clear();
225     }
226    
227    
228     void Props::DoVelCorrelation() {
229     strcpy(velList,file);
230     strcat(velList,"vcorr");
231    
232     counter = nAtoms;
233    
234     vector<double> time;
235     vector<double> pxVel;
236     vector<double> pyVel;
237     vector<double> pzVel;
238     vector<double> avg;
239    
240     ifstream inputer(file);
241    
242     if (counter5 == 0) {
243     cout << "Note: VCorr should only be run on a .dump file\n";
244     return;
245     }
246    
247     //Load all velocities into the appropriate vectors declared above.
248     cout << "Reading in Trajectory...\n";
249     coorCount = 0;
250     inputer.getline(inLine,999,'\n');
251     while (!inputer.eof()) {
252     inputer.getline(inLine,999,'\n');
253     token = strtok(inLine, delimiter);
254     strcpy(inValue,token);
255     time.push_back(atof(inValue));
256     for (i = 0; i < nAtoms; i++) {
257     inputer.getline(inLine,999,'\n');
258     token = strtok(inLine, delimiter);
259     token = strtok(NULL, delimiter);
260     token = strtok(NULL, delimiter);
261     token = strtok(NULL, delimiter);
262     token = strtok(NULL, delimiter);
263     strcpy(inValue,token);
264     pxVel.push_back(atof(inValue));
265     token = strtok(NULL, delimiter);
266     strcpy(inValue,token);
267     pyVel.push_back(atof(inValue));
268     token = strtok(NULL, delimiter);
269     strcpy(inValue,token);
270     pzVel.push_back(atof(inValue));
271     }
272     coorCount++;
273     inputer.getline(inLine,999,'\n');
274     }
275    
276     // Do Velocity Correlation Calculations...
277     cout << "Performing VCorr Calculations...\n";
278     for (j = 0; j < coorCount-1; j++) {
279     sum = 0;
280     blackHole = 0;
281     for (i = 0; (i + j) < coorCount-1; i++) {
282     for (k = 0; k < nAtoms; k++) {
283     vX1 = pxVel[(k + (i * nAtoms))];
284     vY1 = pyVel[(k + (i * nAtoms))];
285     vZ1 = pzVel[(k + (i * nAtoms))];
286     vX2 = pxVel[k + (i * nAtoms) + (j * nAtoms)];
287     vY2 = pyVel[k + (i * nAtoms) + (j * nAtoms)];
288     vZ2 = pzVel[k + (i * nAtoms) + (j * nAtoms)];
289     velDot = vX1*vX2 + vY1*vY2 + vZ1*vZ2;
290     velDotNaught = vX1*vX1 + vY1*vY1 + vZ1*vZ1;
291     dotDivide = velDot / velDotNaught;
292    
293     sum += dotDivide;
294     }
295     blackHole++;
296     }
297     average = sum / nAtoms;
298     average /= blackHole;
299     avg.push_back(average);
300     }
301    
302     // Write out the results...
303     ofstream deltaOut(velList);
304     deltaOut << "#dt\tVCorr\n";
305     for (i=0; i < coorCount-1; i++)
306     deltaOut << (time[i]-time[0]) << "\t" << avg[i] << "\n";
307     inputer.clear();
308     }
309    
310    
311     void Props::DoGofR() {
312     strcpy(gofrList,file);
313     strcat(gofrList,"gofr");
314    
315     count = nAtoms/4;
316    
317     vector<double> pxSpot;
318     vector<double> pySpot;
319     vector<double> pzSpot;
320     // vector<int> pAtomIdent;
321     int *histogram = new int[fineGrain];
322     int *OOhistogram = new int[fineGrain];
323     int *OHhistogram = new int[fineGrain];
324     int *HHhistogram = new int[fineGrain];
325     double *prValues= new double[fineGrain];
326     double *OOGofR = new double[fineGrain];
327     double *HHGofR = new double[fineGrain];
328     double *OHGofR = new double[fineGrain];
329     double *GofR = new double[fineGrain];
330    
331     ifstream inputer(file);
332    
333     //Load all positions into the appropriate vectors declared above.
334     coorCount = 0;
335     cout << "Reading in Trajectory...\n";
336     inputer.getline(inLine,999,'\n');
337     inputer.getline(inLine,999,'\n');
338     token = strtok(inLine,delimiter);
339     token = strtok(NULL,delimiter);
340     strcpy(inValue,token);
341     boxx = atof(inValue);
342     token = strtok(NULL,delimiter);
343     strcpy(inValue,token);
344     boxy = atof(inValue);
345     token = strtok(NULL,delimiter);
346     strcpy(inValue,token);
347     boxz = atof(inValue);
348     while (!inputer.eof()) {
349     for (i = 0; i < nAtoms; i++) {
350     inputer.getline(inLine,999,'\n');
351     token = strtok(inLine, delimiter);
352     token = strtok(NULL, delimiter);
353     strcpy(inValue,token);
354     pxSpot.push_back(atof(inValue));
355     token = strtok(NULL, delimiter);
356     strcpy(inValue,token);
357     pySpot.push_back(atof(inValue));
358     token = strtok(NULL, delimiter);
359     strcpy(inValue,token);
360     pzSpot.push_back(atof(inValue));
361     }
362     coorCount++;
363     inputer.getline(inLine,999,'\n');
364     inputer.getline(inLine,999,'\n');
365     }
366    
367     boxLength = boxy;
368     deltaR = (boxLength/2)/fineGrain;
369     rho = (nAtoms / pow(boxLength, 3));
370     constant = (4.0 * M_PI * rho) / 3.0;
371     rhoO = ((nAtoms/4) / pow(boxLength, 3));
372     rhoH = ((nAtoms/2) / pow(boxLength, 3));
373     rhoX = ((nAtoms) / pow(boxLength, 3));
374     constantX = (4.0 * M_PI * rhoX) / 3.0;
375     constantO = (4.0 * M_PI * rhoO) / 3.0;
376     constantH = (4.0 * M_PI * rhoH) / 3.0;
377    
378     //Zero out all of the g(r)'s and histogram bins
379     for (i = 0; i < fineGrain; i++) {
380     prValues[i] = 0;
381     OOGofR[i] = 0;
382     HHGofR[i] = 0;
383     OHGofR[i] = 0;
384     GofR[i] = 0;
385     histogram[i] = 0;
386     OOhistogram[i] = 0;
387     HHhistogram[i] = 0;
388     OHhistogram[i] = 0;
389     }
390    
391     cout << "Performing G(r) Calculations...\n";
392    
393     // The main calculation loops...
394     for (j = 0; j < coorCount; j++) {
395     for (i = 0; i < nAtoms-1; i++) {
396     rX1 = pxSpot[(i + (j * nAtoms))];
397     rY1 = pySpot[(i + (j * nAtoms))];
398     rZ1 = pzSpot[(i + (j * nAtoms))];
399     for (k = 1; (k + i) < nAtoms; k++) {
400     rX2 = pxSpot[(k + i + (j * nAtoms))];
401     rY2 = pySpot[(k + i + (j * nAtoms))];
402     rZ2 = pzSpot[(k + i + (j * nAtoms))];
403    
404     rXij = (rX1-rX2);
405     rYij = (rY1-rY2);
406     rZij = (rZ1-rZ2);
407    
408     //Do minimum image filter...
409     rXij -= boxLength*copysign(1.0,rXij)*floor(fabs(rXij/boxLength)+0.5);
410     rYij -= boxLength*copysign(1.0,rYij)*floor(fabs(rYij/boxLength)+0.5);
411     rZij -= boxLength*copysign(1.0,rZij)*floor(fabs(rZij/boxLength)+0.5);
412     rijDist = sqrt(rXij*rXij + rYij*rYij + rZij*rZij);
413    
414     //Add to appropriate histogram bin...
415     bin = (int)(rijDist / deltaR);
416     if (bin < fineGrain) {
417     if (pAtomIdent[i] == 3 && pAtomIdent[(k+i)] == 3) {
418     OOhistogram[bin] = OOhistogram[bin]+2;
419     }
420     else if(pAtomIdent[i] == 4 && pAtomIdent[(k+i)] == 4) {
421     HHhistogram[bin] = HHhistogram[bin]+2;
422     }
423     else if(pAtomIdent[i] == 3 && pAtomIdent[(k+i)] == 4) {
424     OHhistogram[bin] = OHhistogram[bin]+2;
425     }
426     else if(pAtomIdent[i] == 4 && pAtomIdent[(k+i)] == 3) {
427     OHhistogram[bin] = OHhistogram[bin]+2;
428     }
429     else if (counter5 > 0) {
430     histogram[bin] = histogram[bin]+2;
431     }
432     }
433     }
434     }
435    
436     // Calculate the G(r) values...
437     for (i = 0; i < fineGrain; i++) {
438     rLower = i * deltaR;
439     rUpper = rLower + deltaR;
440     nIdeal = constant * (pow(rUpper, 3) - pow(rLower, 3));
441     nIdealO = constantO * (pow(rUpper, 3) - pow(rLower, 3));
442     nIdealH = constantH * (pow(rUpper, 3) - pow(rLower, 3));
443    
444     GofR[i] = histogram[i] / (coorCount * nAtoms * nIdeal);
445     OOGofR[i] = OOhistogram[i] / (coorCount * (nAtoms/4) * nIdealO);
446     HHGofR[i] = HHhistogram[i] / (coorCount * (nAtoms/2) * nIdealH);
447     //Check if this one works...
448     OHGofR[i] = OHhistogram[i] / (coorCount * (nAtoms/4) * nIdealO * 4);
449    
450     prValues[i] = rLower + (deltaR/2);
451     }
452     }
453    
454     // Write out the results...
455     strcpy(sepFiler, gofrList);
456     ofstream deltaOut(gofrList);
457     if (counter5 > 0) {
458     deltaOut << "#rDist\tSSDG(r)\n";
459     for (i=0; i<fineGrain; i++)
460     deltaOut << prValues[i] << "\t" << GofR[i] << "\n";
461     }
462     else {
463     deltaOut << "#rDist\tOOG(r)\tHHG(r)\tOHG(r)\n";
464     for (i=0; i < fineGrain; i++)
465     deltaOut << prValues[i] << "\t" << OOGofR[i] << " \t"
466     << HHGofR[i] << " \t" << OHGofR[i] << "\n";
467     }
468     inputer.clear();
469     delete[] OOhistogram; delete[] OHhistogram; delete[] HHhistogram;
470     delete[] prValues; delete[] GofR; delete[] OOGofR;
471     delete[] OHGofR; delete[] HHGofR;
472     }
473    
474    
475     void Props::DoDipoleCorr() {
476     strcpy(allList,file);
477     strcat(allList,"dcorr");
478    
479     counter = nAtoms;
480    
481     if (counter5 > 0) {
482     cout << "Not yet working for this file type. Use an .xyz file.\n";
483     return;
484     }
485     else {
486     nAtoms = nAtoms/4; //A quick hack for SSD water (4 atoms in each SSD)
487     counter = nAtoms;
488     }
489    
490     vector<double> time;
491     vector<double> pUnitX;
492     vector<double> pUnitY;
493     vector<double> pUnitZ;
494     double dot;
495     double zeroval;
496     int nloops;
497    
498     ifstream inputer(file);
499    
500     // Load all positions and unit vectors into the appropriate
501     // vectors declared above.
502     coorCount = 0;
503     cout << "Reading in Trajectory...\n";
504     inputer.getline(inLine,999,'\n');
505     while (!inputer.eof()) {
506     inputer.getline(inLine,999,'\n');
507     token = strtok(inLine, delimiter);
508     strcpy(inValue,token);
509     time.push_back(atof(inValue));
510     for (i = 0; i < nAtoms; i++) {
511     inputer.getline(inLine,999,'\n');
512     token = strtok(inLine, delimiter);
513     token = strtok(NULL, delimiter);
514     token = strtok(NULL, delimiter);
515     token = strtok(NULL, delimiter);
516     token = strtok(NULL, delimiter);
517     strcpy(inValue,token);
518     pUnitX.push_back(atof(inValue));
519     token = strtok(NULL, delimiter);
520     strcpy(inValue,token);
521     pUnitY.push_back(atof(inValue));
522     token = strtok(NULL, delimiter);
523     strcpy(inValue,token);
524     pUnitZ.push_back(atof(inValue));
525     inputer.getline(inLine,999,'\n');
526     inputer.getline(inLine,999,'\n');
527     inputer.getline(inLine,999,'\n');
528     }
529     coorCount++;
530     inputer.getline(inLine,999,'\n');
531     }
532    
533     double *mx = new double[coorCount];
534     double *my = new double[coorCount];
535     double *mz = new double[coorCount];
536     double *avg = new double[coorCount];
537    
538     cout << "Performing DCorr Calculations...\n";
539    
540     // Compute total magnetizations
541     for (j = 0; j < coorCount; j++) {
542     mx[j] = 0;
543     my[j] = 0;
544     mz[j] = 0;
545     for (k = 0; k < counter; k++) {
546     mx[j] += pUnitX[k + (j*counter)];
547     my[j] += pUnitY[k + (j*counter)];
548     mz[j] += pUnitZ[k + (j*counter)];
549     }
550     }
551    
552     // Perform time correlations
553     for (j = 0; j < coorCount-1; j++) {
554     nloops = 0;
555     dot = 0.0;
556    
557     for (i = 0 ; (i+j)<coorCount-1; i++) {
558    
559     dot += mx[i]*mx[i+j] + my[i]*my[i+j] + mz[i]*mz[i+j];
560     nloops++;
561     }
562    
563     avg[j] = dot / (double)nloops;
564    
565     if (j == 0) {
566     zeroval = avg[j];
567     }
568    
569     avg[j] = avg[j] / zeroval;
570     }
571    
572     // Write out the results...
573     ofstream deltaOut(allList);
574     deltaOut << "#dt\tDipole Corr\n";
575     for (i = 0; i < time.size()-1; i++)
576     deltaOut << (time[i]-time[0]) << "\t" << avg[i] << "\n";
577    
578     inputer.clear();
579     delete[] mx; delete[] my; delete[] mz; delete[] avg;
580     }
581    
582     void Props::DoOrientCorr() {
583     strcpy(allList,file);
584     strcat(allList,"mucorr");
585    
586     counter = nAtoms;
587    
588     if (counter5 > 0) {
589     cout << "Not yet working for this file type. Use an .xyz file.\n";
590     return;
591     }
592     else {
593     nAtoms = nAtoms/4; //A quick hack for SSD water (4 atoms in each SSD)
594     counter = nAtoms;
595     }
596    
597     vector<double> time;
598     vector<double> pUnitX;
599     vector<double> pUnitY;
600     vector<double> pUnitZ;
601     double dipDot;
602     double dipDotNaught;
603     double legendreNum;
604    
605     ifstream inputer(file);
606    
607     // Load all positions and unit vectors into the appropriate
608     // vectors declared above.
609     coorCount = 0;
610     cout << "Reading in Trajectory...\n";
611     inputer.getline(inLine,999,'\n');
612     while (!inputer.eof()) {
613     inputer.getline(inLine,999,'\n');
614     token = strtok(inLine, delimiter);
615     strcpy(inValue,token);
616     time.push_back(atof(inValue));
617     for (i = 0; i < counter; i++) {
618     inputer.getline(inLine,999,'\n');
619     token = strtok(inLine, delimiter);
620     token = strtok(NULL, delimiter);
621     token = strtok(NULL, delimiter);
622     token = strtok(NULL, delimiter);
623     token = strtok(NULL, delimiter);
624     strcpy(inValue,token);
625     pUnitX.push_back(atof(inValue));
626     token = strtok(NULL, delimiter);
627     strcpy(inValue,token);
628     pUnitY.push_back(atof(inValue));
629     token = strtok(NULL, delimiter);
630     strcpy(inValue,token);
631     pUnitZ.push_back(atof(inValue));
632     inputer.getline(inLine,999,'\n'); // Jump over O and H atom lines
633     inputer.getline(inLine,999,'\n'); // "
634     inputer.getline(inLine,999,'\n'); // "
635     }
636     coorCount++;
637     inputer.getline(inLine,999,'\n');
638     }
639    
640     double *pC1 = new double[coorCount];
641     double *pC2 = new double[coorCount];
642     double *dipDotSum = new double[coorCount];
643     double *dipDotSum2 = new double[coorCount];
644    
645     cout << "Performing muCorr Calculations...\n";
646     // The main calculation loops...
647     for (j = 0; j < coorCount-1; j++) {
648     dipDotSum[j] = 0;
649     dipDotSum2[j] = 0;
650     blackHole = 0;
651     for (i = 0; (i + j) < coorCount-1; i++) {
652     for (k = 0; k < counter; k++) {
653     uX1 = pUnitX[k + (i * counter)];
654     uY1 = pUnitY[k + (i * counter)];
655     uZ1 = pUnitZ[k + (i * counter)];
656     uX2 = pUnitX[k + (i * counter) + (j * counter)];
657     uY2 = pUnitY[k + (i * counter) + (j * counter)];
658     uZ2 = pUnitZ[k + (i * counter) + (j * counter)];
659    
660     dipDot = uX1*uX2 + uY1*uY2 + uZ1*uZ2;
661     legendreNum = 1.5 * dipDot * dipDot - 0.5;
662    
663     dipDotSum[j] += dipDot;
664     dipDotSum2[j] += legendreNum;
665     }
666     blackHole++;
667     }
668    
669     dipDotSum[j] = dipDotSum[j] / (blackHole*counter);
670     dipDotSum2[j] = dipDotSum2[j] / (blackHole*counter);
671     }
672    
673     for (j = 0; j < coorCount-1; j++) {
674     pC1[j] = dipDotSum[j] / dipDotSum[0];
675     pC2[j] = dipDotSum2[j] / dipDotSum2[0];
676     }
677    
678    
679     ofstream deltaOut(allList);
680     deltaOut << "#dt\tC1(t)\tC2(t)\n";
681     for (i = 0; i < coorCount-1; i++)
682     deltaOut << (time[i]-time[0]) << "\t" << pC1[i] << "\t" << pC2[i] << "\n";
683     inputer.clear();
684    
685     delete[] pC1; delete[] pC2;
686     }
687    
688    
689     void Props::DoMagFileCorr() {
690     strcpy(allList,file);
691     strcat(allList,"magcorr");
692    
693     vector<double> time;
694     vector<double> pUnitX;
695     vector<double> pUnitY;
696     vector<double> pUnitZ;
697     double dot;
698     double zeroval;
699     int nloops;
700    
701     ifstream inputer(file);
702    
703     // Load all positions and unit vectors into the appropriate
704     // vectors declared above.
705     coorCount = 0;
706     cout << "Reading in Trajectory...\n";
707     inputer.getline(inLine,999,'\n');
708     while (!inputer.eof()) {
709     token = strtok(inLine, delimiter);
710     strcpy(inValue,token);
711     time.push_back(atof(inValue));
712     token = strtok(NULL, delimiter);
713     strcpy(inValue,token);
714     pUnitX.push_back(atof(inValue));
715     token = strtok(NULL, delimiter);
716     strcpy(inValue,token);
717     pUnitY.push_back(atof(inValue));
718     token = strtok(NULL, delimiter);
719     strcpy(inValue,token);
720     pUnitZ.push_back(atof(inValue));
721     coorCount++;
722     inputer.getline(inLine,999,'\n');
723     }
724    
725     double *avg = new double[time.size()];
726    
727     cout << "Performing magCorr Calculations...\n";
728    
729     // Perform time correlations
730     for (j = 0; j < coorCount-1; j++) {
731     nloops = 0;
732     dot = 0.0;
733    
734     for (i = 0 ; (i+j)<coorCount-1; i++) {
735     dot += pUnitX[i]*pUnitX[i+j] + pUnitY[i]*pUnitY[i+j]
736     + pUnitZ[i]*pUnitZ[i+j];
737     nloops++;
738     }
739    
740     avg[j] = dot / (double)nloops;
741    
742     if (j == 0) {
743     zeroval = avg[j];
744     }
745     avg[j] = avg[j] / zeroval;
746     }
747    
748     // Write out the results...
749     ofstream deltaOut(allList);
750     deltaOut << "#dt\tMag Corr\n";
751     for (i = 0; i < time.size()-1; i++)
752     deltaOut << (time[i]-time[0]) << "\t" << avg[i] << "\n";
753    
754     inputer.clear();
755     delete[] avg;
756     }
757    
758     void Props::DoCosineCorr() {
759     strcpy(gofrList,file);
760     strcat(gofrList,"coscorr");
761    
762     counter = nAtoms/4;
763    
764     vector<double> pxSpot;
765     vector<double> pySpot;
766     vector<double> pzSpot;
767     vector<double> pUnitX;
768     vector<double> pUnitY;
769     vector<double> pUnitZ;
770     double dot;
771     int *histogram = new int[fineGrain];
772     double *prValues= new double[fineGrain];
773     double *avgCos = new double[fineGrain];
774    
775     if (counter5 > 0) {
776     cout << "File type currently unsupported. Please use a .xyz file.\n";
777     return;
778     }
779    
780     ifstream inputer(file);
781    
782     //Load all positions into the appropriate vectors declared above.
783     coorCount = 0;
784     cout << "Reading in Trajectory...\n";
785     inputer.getline(inLine,999,'\n');
786     inputer.getline(inLine,999,'\n');
787     token = strtok(inLine,delimiter);
788     token = strtok(NULL,delimiter);
789     strcpy(inValue,token);
790     boxx = atof(inValue);
791     token = strtok(NULL,delimiter);
792     strcpy(inValue,token);
793     boxy = atof(inValue);
794     token = strtok(NULL,delimiter);
795     strcpy(inValue,token);
796     boxz = atof(inValue);
797     while (!inputer.eof()) {
798     for (i = 0; i < counter; i++) {
799     inputer.getline(inLine,999,'\n');
800     token = strtok(inLine, delimiter);
801     token = strtok(NULL, delimiter);
802     strcpy(inValue,token);
803     pxSpot.push_back(atof(inValue));
804     token = strtok(NULL, delimiter);
805     strcpy(inValue,token);
806     pySpot.push_back(atof(inValue));
807     token = strtok(NULL, delimiter);
808     strcpy(inValue,token);
809     pzSpot.push_back(atof(inValue));
810     token = strtok(NULL, delimiter);
811     strcpy(inValue,token);
812     pUnitX.push_back(atof(inValue));
813     token = strtok(NULL, delimiter);
814     strcpy(inValue,token);
815     pUnitY.push_back(atof(inValue));
816     token = strtok(NULL, delimiter);
817     strcpy(inValue,token);
818     pUnitZ.push_back(atof(inValue));
819     inputer.getline(inLine,999,'\n');
820     inputer.getline(inLine,999,'\n');
821     inputer.getline(inLine,999,'\n');
822     }
823     coorCount++;
824     inputer.getline(inLine,999,'\n');
825     inputer.getline(inLine,999,'\n');
826     }
827    
828     boxLength = boxy;
829     deltaR = (boxLength/2)/fineGrain;
830    
831     //Zero out all of the g(r)'s and histogram bins
832     for (i = 0; i < fineGrain; i++) {
833     prValues[i] = 0;
834     histogram[i] = 0;
835     avgCos[i] = 0;
836     }
837    
838     cout << "Performing cosCorr Calculations...\n";
839     // Main sorting loops...
840     for (j = 0; j < coorCount; j++) {
841     for (i = 0; i < counter-1; i++) {
842     rX1 = pxSpot[(i + (j * counter))];
843     rY1 = pySpot[(i + (j * counter))];
844     rZ1 = pzSpot[(i + (j * counter))];
845     uX1 = pUnitX[(i + (j * counter))];
846     uY1 = pUnitY[(i + (j * counter))];
847     uZ1 = pUnitZ[(i + (j * counter))];
848     for (k = 1; (k + i) < counter; k++) {
849     rX2 = pxSpot[(k + i + (j * counter))];
850     rY2 = pySpot[(k + i + (j * counter))];
851     rZ2 = pzSpot[(k + i + (j * counter))];
852    
853     rXij = (rX1-rX2);
854     rYij = (rY1-rY2);
855     rZij = (rZ1-rZ2);
856    
857     //Do minimum image filter...
858     rXij -= boxLength*copysign(1.0,rXij)*floor(fabs(rXij/boxLength)+0.5);
859     rYij -= boxLength*copysign(1.0,rYij)*floor(fabs(rYij/boxLength)+0.5);
860     rZij -= boxLength*copysign(1.0,rZij)*floor(fabs(rZij/boxLength)+0.5);
861     rijDist = sqrt(rXij*rXij + rYij*rYij + rZij*rZij);
862    
863     //Now do the average cosine b/t dipoles
864     uX2 = pUnitX[(k + i + (j * counter))];
865     uY2 = pUnitY[(k + i + (j * counter))];
866     uZ2 = pUnitZ[(k + i + (j * counter))];
867    
868     dot = uX1*uX2 + uY1*uY2 + uZ1*uZ2;
869    
870     //Add to appropriate histogram bin...
871     bin = (int)(rijDist / deltaR);
872     if (bin < fineGrain) {
873     histogram[bin] = histogram[bin] + 1;
874     avgCos[bin] = avgCos[bin] + dot;
875     }
876     }
877     }
878     }
879    
880     for (i = 0; i < fineGrain; i++) {
881     if (histogram[i] != 0)
882     avgCos[i] = avgCos[i] / histogram[i];
883     }
884    
885     for (i = 0; i < fineGrain; i++) {
886     rLower = i * deltaR;
887     prValues[i] = rLower + (deltaR/2);
888     }
889    
890     // Write out the results...
891     strcpy(sepFiler, gofrList);
892     ofstream deltaOut(gofrList);
893     deltaOut << "#Radius\tcosCorr\n";
894     for (i=0; i<fineGrain; i++)
895     deltaOut << prValues[i] << "\t" << avgCos[i] << "\n";
896    
897     inputer.clear();
898     delete[] prValues; delete[] histogram; delete[] avgCos;
899     }
900    
901    
902    
903     int main(int argc, char *argv[]) {
904     char choice[200], inplace[200], atomChar[10], inLine[1000], inValue[200];
905     char temp[4];
906     char fracString[4] = "";
907     char *token;
908     const char *delimit = " \t\n";
909     int selection,nParticles,i,grain;
910     vector<int> pAtomType;
911     string strungName;
912    
913     // Command line input format test
914     if (argc != 3) {
915     cout
916     << "Usage: " << argv[0] << " <file name> <correlation choice>\n"
917     << "\n"
918     << "Correlation choice options:\n"
919     << "\n"
920     << " rcorr - Root Mean-Squared Displacement\n"
921     << " vcorr - Velocity Autocorrelation Function\n"
922     << " gofr - Pair Correlation Function\n"
923     << " dcorr - Collective Dipole Correlation Function\n"
924     << " mucorr - Self Dipolar Orientational Correlation Function (C1 & C2)\n"
925     << " magfile - Collective Dipole Correlation Function\n"
926     << " coscorr - Distance-Dependent Average Cosine Correlation\n";
927     return 0;
928     }
929    
930     ifstream prayer(argv[1]);
931    
932     // Make sure the file exists
933     if (!prayer) {
934     cout << "Unable to open " << argv[1] << " for reading.\n";
935     return 0;
936     }
937    
938     prayer.close();
939    
940     // Build a filename string
941     strungName = argv[1];
942     grain = 500; //detail of g(r) curve (# of bins) - increase for smoother curve
943    
944     // Assign the correlation to perform
945     if (!strcmp(argv[2], "rcorr"))
946     selection = 1;
947     else if (!strcmp(argv[2], "vcorr"))
948     selection = 2;
949     else if (!strcmp(argv[2], "gofr"))
950     selection = 3;
951     else if (!strcmp(argv[2], "dcorr"))
952     selection = 4;
953     else if (!strcmp(argv[2], "mucorr"))
954     selection = 5;
955     else if (!strcmp(argv[2], "magfile"))
956     selection = 6;
957     else if (!strcmp(argv[2], "coscorr"))
958     selection = 7;
959     else {
960     cout
961     << "Not a valid correlation. Choose one of the following:\n"
962     << "\n"
963     << " rcorr - Root Mean-Squared Displacement\n"
964     << " vcorr - Velocity Autocorrelation Function\n"
965     << " gofr - Pair Correlation Function\n"
966     << " dcorr - Collective Dipole Correlation Function\n"
967     << " mucorr - Self Dipolar Orientational Correlation Function (C1 & C2)\n"
968     << " magfile - Collective Dipole Correlation Function\n"
969     << " coscorr - Distance-Dependent Average Cosine Correlation\n";
970     return 0;
971     }
972    
973     // Call up correlation calculations
974    
975     Props analyzer(strungName,delimit,pAtomType,grain);
976    
977     switch( selection ){
978    
979     case 1:
980     analyzer.DoRMSDCorrelation();
981     break;
982    
983     case 2:
984     analyzer.DoVelCorrelation();
985     break;
986    
987     case 3:
988     analyzer.DoGofR();
989     break;
990    
991     case 4:
992     analyzer.DoDipoleCorr();
993     break;
994    
995     case 5:
996     analyzer.DoOrientCorr();
997     break;
998    
999     case 6:
1000     analyzer.DoMagFileCorr();
1001     break;
1002    
1003     case 7:
1004     analyzer.DoCosineCorr();
1005     break;
1006    
1007     default:
1008     cout << "Error in understanding the selection\n";
1009     }
1010    
1011     return 0;
1012     }