ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/madProps/chrisProps.cpp
Revision: 39
Committed: Fri Jul 19 01:37:38 2002 UTC (23 years, 3 months ago) by mmeineke
File size: 27754 byte(s)
Log Message: 
This commit was generated by cvs2svn to compensate for changes in r38, which
included commits to RCS files with non-trunk default branches.

File Contents

# Content
1 // 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 }