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 |
} |