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