1 |
+ |
// Thermodynamic integration is not multiprocessor friendly right now |
2 |
+ |
|
3 |
+ |
|
4 |
+ |
|
5 |
|
#include <iostream> |
6 |
+ |
|
7 |
|
#include <stdlib.h> |
8 |
+ |
|
9 |
|
#include <cstdio> |
10 |
+ |
|
11 |
|
#include <fstream> |
12 |
+ |
|
13 |
|
#include <iomanip> |
14 |
+ |
|
15 |
|
#include <string> |
16 |
+ |
|
17 |
|
#include <cstring> |
18 |
+ |
|
19 |
|
#include <math.h> |
20 |
|
|
21 |
+ |
|
22 |
+ |
|
23 |
|
using namespace std; |
24 |
|
|
12 |
– |
#include "Restraints.hpp" |
13 |
– |
#include "SimInfo.hpp" |
14 |
– |
#include "simError.h" |
25 |
|
|
26 |
+ |
|
27 |
+ |
#include "restraints/Restraints.hpp" |
28 |
+ |
|
29 |
+ |
#include "brains/SimInfo.hpp" |
30 |
+ |
|
31 |
+ |
#include "utils/simError.h" |
32 |
+ |
|
33 |
+ |
|
34 |
+ |
|
35 |
|
#define PI 3.14159265359 |
36 |
+ |
|
37 |
|
#define TWO_PI 6.28318530718 |
38 |
|
|
39 |
+ |
|
40 |
+ |
|
41 |
|
Restraints::Restraints(double lambdaVal, double lambdaExp){ |
42 |
+ |
|
43 |
|
lambdaValue = lambdaVal; |
44 |
+ |
|
45 |
|
lambdaK = lambdaExp; |
46 |
|
|
47 |
+ |
std::vector<double> resConsts; |
48 |
+ |
|
49 |
|
const char *jolt = " \t\n;,"; |
50 |
|
|
51 |
+ |
|
52 |
+ |
|
53 |
|
#ifdef IS_MPI |
54 |
+ |
|
55 |
|
if(worldRank == 0 ){ |
56 |
+ |
|
57 |
|
#endif // is_mpi |
58 |
|
|
59 |
+ |
|
60 |
+ |
|
61 |
|
strcpy(springName, "HarmSpringConsts.txt"); |
62 |
+ |
|
63 |
|
|
64 |
+ |
|
65 |
|
ifstream springs(springName); |
66 |
+ |
|
67 |
|
|
68 |
+ |
|
69 |
|
if (!springs) { |
70 |
+ |
|
71 |
|
sprintf(painCave.errMsg, |
72 |
< |
"In Restraints: Unable to open HarmSpringConsts.txt for reading.\n" |
73 |
< |
"\tDefault spring constants will be loaded. If you want to specify\n" |
74 |
< |
"\tspring constants, include a three line HarmSpringConsts.txt file\n" |
75 |
< |
"\tin the current directory.\n"); |
72 |
> |
|
73 |
> |
"Unable to open HarmSpringConsts.txt for reading, so the\n" |
74 |
> |
|
75 |
> |
"\tdefault spring constants will be loaded. If you want\n" |
76 |
> |
|
77 |
> |
"\tto specify spring constants, include a three line\n" |
78 |
> |
|
79 |
> |
"\tHarmSpringConsts.txt file in the execution directory.\n"); |
80 |
> |
|
81 |
|
painCave.severity = OOPSE_WARNING; |
82 |
+ |
|
83 |
|
painCave.isFatal = 0; |
84 |
+ |
|
85 |
|
simError(); |
86 |
+ |
|
87 |
|
|
88 |
+ |
|
89 |
|
// load default spring constants |
90 |
+ |
|
91 |
|
kDist = 6; // spring constant in units of kcal/(mol*ang^2) |
92 |
+ |
|
93 |
|
kTheta = 7.5; // in units of kcal/mol |
94 |
+ |
|
95 |
|
kOmega = 13.5; // in units of kcal/mol |
96 |
+ |
|
97 |
|
} else { |
98 |
+ |
|
99 |
|
|
100 |
+ |
|
101 |
|
springs.getline(inLine,999,'\n'); |
102 |
+ |
|
103 |
+ |
// the file is blank! |
104 |
+ |
|
105 |
+ |
if (springs.eof()){ |
106 |
+ |
|
107 |
+ |
sprintf(painCave.errMsg, |
108 |
+ |
|
109 |
+ |
"HarmSpringConsts.txt file is not valid.\n" |
110 |
+ |
|
111 |
+ |
"\tThe file should contain four rows, the last three containing\n" |
112 |
+ |
|
113 |
+ |
"\ta label and the spring constant value. They should be listed\n" |
114 |
+ |
|
115 |
+ |
"\tin the following order: kDist (positional restrant), kTheta\n" |
116 |
+ |
|
117 |
+ |
"\t(rot. restraint: deflection of z-axis), and kOmega (rot.\n" |
118 |
+ |
|
119 |
+ |
"\trestraint: rotation about the z-axis).\n"); |
120 |
+ |
|
121 |
+ |
painCave.severity = OOPSE_ERROR; |
122 |
+ |
|
123 |
+ |
painCave.isFatal = 1; |
124 |
+ |
|
125 |
+ |
simError(); |
126 |
+ |
|
127 |
+ |
} |
128 |
+ |
|
129 |
+ |
// read in spring constants and check to make sure it is a valid file |
130 |
+ |
|
131 |
|
springs.getline(inLine,999,'\n'); |
132 |
< |
token = strtok(inLine,jolt); |
133 |
< |
token = strtok(NULL,jolt); |
134 |
< |
strcpy(inValue,token); |
135 |
< |
kDist = (atof(inValue)); |
136 |
< |
springs.getline(inLine,999,'\n'); |
137 |
< |
token = strtok(inLine,jolt); |
138 |
< |
token = strtok(NULL,jolt); |
139 |
< |
strcpy(inValue,token); |
140 |
< |
kTheta = (atof(inValue)); |
141 |
< |
springs.getline(inLine,999,'\n'); |
142 |
< |
token = strtok(inLine,jolt); |
143 |
< |
token = strtok(NULL,jolt); |
144 |
< |
strcpy(inValue,token); |
145 |
< |
kOmega = (atof(inValue)); |
146 |
< |
springs.close(); |
132 |
> |
|
133 |
> |
while (!springs.eof()){ |
134 |
> |
|
135 |
> |
if (NULL != inLine){ |
136 |
> |
|
137 |
> |
token = strtok(inLine,jolt); |
138 |
> |
|
139 |
> |
token = strtok(NULL,jolt); |
140 |
> |
|
141 |
> |
if (NULL != token){ |
142 |
> |
|
143 |
> |
strcpy(inValue,token); |
144 |
> |
|
145 |
> |
resConsts.push_back(atof(inValue)); |
146 |
> |
|
147 |
> |
} |
148 |
> |
|
149 |
> |
} |
150 |
> |
|
151 |
> |
springs.getline(inLine,999,'\n'); |
152 |
> |
|
153 |
> |
} |
154 |
> |
|
155 |
> |
if (resConsts.size() == 3){ |
156 |
> |
|
157 |
> |
kDist = resConsts[0]; |
158 |
> |
|
159 |
> |
kTheta = resConsts[1]; |
160 |
> |
|
161 |
> |
kOmega = resConsts[2]; |
162 |
> |
|
163 |
> |
} |
164 |
> |
|
165 |
> |
else { |
166 |
> |
|
167 |
> |
sprintf(painCave.errMsg, |
168 |
> |
|
169 |
> |
"HarmSpringConsts.txt file is not valid.\n" |
170 |
> |
|
171 |
> |
"\tThe file should contain four rows, the last three containing\n" |
172 |
> |
|
173 |
> |
"\ta label and the spring constant value. They should be listed\n" |
174 |
> |
|
175 |
> |
"\tin the following order: kDist (positional restrant), kTheta\n" |
176 |
> |
|
177 |
> |
"\t(rot. restraint: deflection of z-axis), and kOmega (rot.\n" |
178 |
> |
|
179 |
> |
"\trestraint: rotation about the z-axis).\n"); |
180 |
> |
|
181 |
> |
painCave.severity = OOPSE_ERROR; |
182 |
> |
|
183 |
> |
painCave.isFatal = 1; |
184 |
> |
|
185 |
> |
simError(); |
186 |
> |
|
187 |
> |
} |
188 |
> |
|
189 |
|
} |
190 |
+ |
|
191 |
|
#ifdef IS_MPI |
192 |
+ |
|
193 |
|
} |
194 |
+ |
|
195 |
|
|
196 |
+ |
|
197 |
|
MPI_Bcast(&kDist, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); |
198 |
+ |
|
199 |
|
MPI_Bcast(&kTheta, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); |
200 |
+ |
|
201 |
|
MPI_Bcast(&kOmega, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); |
202 |
+ |
|
203 |
|
|
204 |
+ |
|
205 |
|
sprintf( checkPointMsg, |
206 |
+ |
|
207 |
|
"Sucessfully opened and read spring file.\n"); |
208 |
+ |
|
209 |
|
MPIcheckPoint(); |
210 |
|
|
211 |
+ |
|
212 |
+ |
|
213 |
|
#endif // is_mpi |
214 |
+ |
|
215 |
|
|
216 |
+ |
|
217 |
|
sprintf(painCave.errMsg, |
218 |
+ |
|
219 |
|
"The spring constants for thermodynamic integration are:\n" |
220 |
+ |
|
221 |
|
"\tkDist = %lf\n" |
222 |
+ |
|
223 |
|
"\tkTheta = %lf\n" |
224 |
+ |
|
225 |
|
"\tkOmega = %lf\n", kDist, kTheta, kOmega); |
226 |
+ |
|
227 |
|
painCave.severity = OOPSE_INFO; |
228 |
+ |
|
229 |
|
painCave.isFatal = 0; |
230 |
+ |
|
231 |
|
simError(); |
232 |
+ |
|
233 |
|
} |
234 |
|
|
235 |
+ |
|
236 |
+ |
|
237 |
|
Restraints::~Restraints(){ |
238 |
+ |
|
239 |
|
} |
240 |
|
|
241 |
+ |
|
242 |
+ |
|
243 |
|
void Restraints::Calc_rVal(double position[3], int currentMol){ |
244 |
+ |
|
245 |
|
delRx = position[0] - cofmPosX[currentMol]; |
246 |
+ |
|
247 |
|
delRy = position[1] - cofmPosY[currentMol]; |
248 |
+ |
|
249 |
|
delRz = position[2] - cofmPosZ[currentMol]; |
250 |
|
|
251 |
+ |
|
252 |
+ |
|
253 |
|
return; |
254 |
+ |
|
255 |
|
} |
256 |
|
|
257 |
+ |
|
258 |
+ |
|
259 |
|
void Restraints::Calc_body_thetaVal(double matrix[3][3], int currentMol){ |
260 |
+ |
|
261 |
|
ub0x = matrix[0][0]*uX0[currentMol] + matrix[0][1]*uY0[currentMol] |
262 |
+ |
|
263 |
|
+ matrix[0][2]*uZ0[currentMol]; |
264 |
+ |
|
265 |
|
ub0y = matrix[1][0]*uX0[currentMol] + matrix[1][1]*uY0[currentMol] |
266 |
+ |
|
267 |
|
+ matrix[1][2]*uZ0[currentMol]; |
268 |
+ |
|
269 |
|
ub0z = matrix[2][0]*uX0[currentMol] + matrix[2][1]*uY0[currentMol] |
270 |
+ |
|
271 |
|
+ matrix[2][2]*uZ0[currentMol]; |
272 |
|
|
273 |
+ |
|
274 |
+ |
|
275 |
|
normalize = sqrt(ub0x*ub0x + ub0y*ub0y + ub0z*ub0z); |
276 |
+ |
|
277 |
|
ub0x = ub0x/normalize; |
278 |
+ |
|
279 |
|
ub0y = ub0y/normalize; |
280 |
+ |
|
281 |
|
ub0z = ub0z/normalize; |
282 |
|
|
283 |
+ |
|
284 |
+ |
|
285 |
|
// Theta is the dot product of the reference and new z-axes |
286 |
+ |
|
287 |
|
theta = acos(ub0z); |
288 |
|
|
289 |
+ |
|
290 |
+ |
|
291 |
|
return; |
292 |
+ |
|
293 |
|
} |
294 |
|
|
295 |
+ |
|
296 |
+ |
|
297 |
|
void Restraints::Calc_body_omegaVal(double matrix[3][3], double zAngle){ |
298 |
+ |
|
299 |
|
double zRotator[3][3]; |
300 |
+ |
|
301 |
|
double tempOmega; |
302 |
+ |
|
303 |
|
double wholeTwoPis; |
304 |
+ |
|
305 |
|
// Use the omega accumulated from the rotation propagation |
306 |
+ |
|
307 |
|
omega = zAngle; |
308 |
|
|
309 |
+ |
|
310 |
+ |
|
311 |
|
// translate the omega into a range between -PI and PI |
312 |
+ |
|
313 |
|
if (omega < -PI){ |
314 |
+ |
|
315 |
|
tempOmega = omega / -TWO_PI; |
316 |
+ |
|
317 |
|
wholeTwoPis = floor(tempOmega); |
318 |
+ |
|
319 |
|
tempOmega = omega + TWO_PI*wholeTwoPis; |
320 |
+ |
|
321 |
|
if (tempOmega < -PI) |
322 |
+ |
|
323 |
|
omega = tempOmega + TWO_PI; |
324 |
+ |
|
325 |
|
else |
326 |
+ |
|
327 |
|
omega = tempOmega; |
328 |
+ |
|
329 |
|
} |
330 |
+ |
|
331 |
|
if (omega > PI){ |
332 |
+ |
|
333 |
|
tempOmega = omega / TWO_PI; |
334 |
+ |
|
335 |
|
wholeTwoPis = floor(tempOmega); |
336 |
+ |
|
337 |
|
tempOmega = omega - TWO_PI*wholeTwoPis; |
338 |
+ |
|
339 |
|
if (tempOmega > PI) |
340 |
+ |
|
341 |
|
omega = tempOmega - TWO_PI; |
342 |
+ |
|
343 |
|
else |
344 |
+ |
|
345 |
|
omega = tempOmega; |
346 |
+ |
|
347 |
|
} |
348 |
|
|
349 |
+ |
|
350 |
+ |
|
351 |
|
vb0x = sin(omega); |
352 |
+ |
|
353 |
|
vb0y = cos(omega); |
354 |
+ |
|
355 |
|
vb0z = 0.0; |
356 |
|
|
357 |
+ |
|
358 |
+ |
|
359 |
|
normalize = sqrt(vb0x*vb0x + vb0y*vb0y + vb0z*vb0z); |
360 |
+ |
|
361 |
|
vb0x = vb0x/normalize; |
362 |
+ |
|
363 |
|
vb0y = vb0y/normalize; |
364 |
+ |
|
365 |
|
vb0z = vb0z/normalize; |
366 |
|
|
367 |
+ |
|
368 |
+ |
|
369 |
|
return; |
370 |
+ |
|
371 |
|
} |
372 |
|
|
373 |
+ |
|
374 |
+ |
|
375 |
|
double Restraints::Calc_Restraint_Forces(vector<StuntDouble*> vecParticles){ |
376 |
+ |
|
377 |
|
double pos[3]; |
378 |
+ |
|
379 |
|
double A[3][3]; |
380 |
+ |
|
381 |
|
double tolerance; |
382 |
+ |
|
383 |
|
double tempPotent; |
384 |
+ |
|
385 |
|
double factor; |
386 |
+ |
|
387 |
|
double spaceTrq[3]; |
388 |
+ |
|
389 |
|
double omegaPass; |
390 |
|
|
391 |
+ |
|
392 |
+ |
|
393 |
|
tolerance = 5.72957795131e-7; |
394 |
|
|
395 |
+ |
|
396 |
+ |
|
397 |
|
harmPotent = 0.0; // zero out the global harmonic potential variable |
398 |
|
|
399 |
+ |
|
400 |
+ |
|
401 |
|
factor = 1 - pow(lambdaValue, lambdaK); |
402 |
|
|
403 |
+ |
|
404 |
+ |
|
405 |
|
for (i=0; i<vecParticles.size(); i++){ |
406 |
+ |
|
407 |
|
if (vecParticles[i]->isDirectional()){ |
408 |
< |
vecParticles[i]->getPos(pos); |
408 |
> |
|
409 |
> |
pos = vecParticles[i]->getPos(); |
410 |
> |
|
411 |
|
vecParticles[i]->getA(A); |
412 |
+ |
|
413 |
|
Calc_rVal( pos, i ); |
414 |
+ |
|
415 |
|
Calc_body_thetaVal( A, i ); |
416 |
+ |
|
417 |
|
omegaPass = vecParticles[i]->getZangle(); |
418 |
+ |
|
419 |
|
Calc_body_omegaVal( A, omegaPass ); |
420 |
|
|
183 |
– |
if (omega > PI || omega < -PI) |
184 |
– |
cout << "oops... " << omega << "\n"; |
421 |
|
|
422 |
+ |
|
423 |
|
// first we calculate the derivatives |
424 |
+ |
|
425 |
|
dVdrx = -kDist*delRx; |
426 |
+ |
|
427 |
|
dVdry = -kDist*delRy; |
428 |
+ |
|
429 |
|
dVdrz = -kDist*delRz; |
430 |
|
|
431 |
+ |
|
432 |
+ |
|
433 |
|
// uTx... and vTx... are the body-fixed z and y unit vectors |
434 |
+ |
|
435 |
|
uTx = 0.0; |
436 |
+ |
|
437 |
|
uTy = 0.0; |
438 |
+ |
|
439 |
|
uTz = 1.0; |
440 |
+ |
|
441 |
|
vTx = 0.0; |
442 |
+ |
|
443 |
|
vTy = 1.0; |
444 |
+ |
|
445 |
|
vTz = 0.0; |
446 |
|
|
447 |
+ |
|
448 |
+ |
|
449 |
|
dVdux = 0; |
450 |
+ |
|
451 |
|
dVduy = 0; |
452 |
+ |
|
453 |
|
dVduz = 0; |
454 |
+ |
|
455 |
|
dVdvx = 0; |
456 |
+ |
|
457 |
|
dVdvy = 0; |
458 |
+ |
|
459 |
|
dVdvz = 0; |
460 |
|
|
461 |
+ |
|
462 |
+ |
|
463 |
|
if (fabs(theta) > tolerance) { |
464 |
+ |
|
465 |
|
dVdux = -(kTheta*theta/sin(theta))*ub0x; |
466 |
+ |
|
467 |
|
dVduy = -(kTheta*theta/sin(theta))*ub0y; |
468 |
+ |
|
469 |
|
dVduz = -(kTheta*theta/sin(theta))*ub0z; |
470 |
+ |
|
471 |
|
} |
472 |
|
|
473 |
+ |
|
474 |
+ |
|
475 |
|
if (fabs(omega) > tolerance) { |
476 |
+ |
|
477 |
|
dVdvx = -(kOmega*omega/sin(omega))*vb0x; |
478 |
+ |
|
479 |
|
dVdvy = -(kOmega*omega/sin(omega))*vb0y; |
480 |
+ |
|
481 |
|
dVdvz = -(kOmega*omega/sin(omega))*vb0z; |
482 |
+ |
|
483 |
|
} |
484 |
|
|
485 |
+ |
|
486 |
+ |
|
487 |
|
// next we calculate the restraint forces and torques |
488 |
+ |
|
489 |
|
restraintFrc[0] = dVdrx; |
490 |
+ |
|
491 |
|
restraintFrc[1] = dVdry; |
492 |
+ |
|
493 |
|
restraintFrc[2] = dVdrz; |
494 |
+ |
|
495 |
|
tempPotent = 0.5*kDist*(delRx*delRx + delRy*delRy + delRz*delRz); |
496 |
|
|
497 |
+ |
|
498 |
+ |
|
499 |
|
restraintTrq[0] = 0.0; |
500 |
+ |
|
501 |
|
restraintTrq[1] = 0.0; |
502 |
+ |
|
503 |
|
restraintTrq[2] = 0.0; |
504 |
|
|
505 |
+ |
|
506 |
+ |
|
507 |
|
if (fabs(omega) > tolerance) { |
508 |
+ |
|
509 |
|
restraintTrq[0] += 0.0; |
510 |
+ |
|
511 |
|
restraintTrq[1] += 0.0; |
512 |
+ |
|
513 |
|
restraintTrq[2] += vTy*dVdvx; |
514 |
+ |
|
515 |
|
tempPotent += 0.5*(kOmega*omega*omega); |
516 |
+ |
|
517 |
|
} |
518 |
+ |
|
519 |
|
if (fabs(theta) > tolerance) { |
520 |
+ |
|
521 |
|
restraintTrq[0] += (uTz*dVduy); |
522 |
+ |
|
523 |
|
restraintTrq[1] += -(uTz*dVdux); |
524 |
+ |
|
525 |
|
restraintTrq[2] += 0.0; |
526 |
+ |
|
527 |
|
tempPotent += 0.5*(kTheta*theta*theta); |
528 |
+ |
|
529 |
|
} |
530 |
|
|
531 |
+ |
|
532 |
+ |
|
533 |
|
for (j = 0; j < 3; j++) { |
534 |
+ |
|
535 |
|
restraintFrc[j] *= factor; |
536 |
+ |
|
537 |
|
restraintTrq[j] *= factor; |
538 |
+ |
|
539 |
|
} |
540 |
|
|
541 |
+ |
|
542 |
+ |
|
543 |
|
harmPotent += tempPotent; |
544 |
|
|
545 |
+ |
|
546 |
+ |
|
547 |
|
// now we need to convert from body-fixed torques to space-fixed torques |
548 |
+ |
|
549 |
|
spaceTrq[0] = A[0][0]*restraintTrq[0] + A[1][0]*restraintTrq[1] |
550 |
+ |
|
551 |
|
+ A[2][0]*restraintTrq[2]; |
552 |
+ |
|
553 |
|
spaceTrq[1] = A[0][1]*restraintTrq[0] + A[1][1]*restraintTrq[1] |
554 |
+ |
|
555 |
|
+ A[2][1]*restraintTrq[2]; |
556 |
+ |
|
557 |
|
spaceTrq[2] = A[0][2]*restraintTrq[0] + A[1][2]*restraintTrq[1] |
558 |
+ |
|
559 |
|
+ A[2][2]*restraintTrq[2]; |
560 |
|
|
561 |
+ |
|
562 |
+ |
|
563 |
|
// now it's time to pass these temporary forces and torques |
564 |
+ |
|
565 |
|
// to the total forces and torques |
566 |
+ |
|
567 |
|
vecParticles[i]->addFrc(restraintFrc); |
568 |
+ |
|
569 |
|
vecParticles[i]->addTrq(spaceTrq); |
570 |
+ |
|
571 |
|
} |
572 |
+ |
|
573 |
|
} |
574 |
|
|
575 |
+ |
|
576 |
+ |
|
577 |
|
// and we can return the appropriately scaled potential energy |
578 |
+ |
|
579 |
|
tempPotent = harmPotent * factor; |
580 |
+ |
|
581 |
|
return tempPotent; |
582 |
+ |
|
583 |
|
} |
584 |
|
|
585 |
+ |
|
586 |
+ |
|
587 |
|
void Restraints::Store_Init_Info(vector<StuntDouble*> vecParticles){ |
588 |
+ |
|
589 |
+ |
int idealSize; |
590 |
+ |
|
591 |
|
double pos[3]; |
592 |
+ |
|
593 |
|
double A[3][3]; |
594 |
+ |
|
595 |
|
double RfromQ[3][3]; |
596 |
+ |
|
597 |
|
double quat0, quat1, quat2, quat3; |
598 |
+ |
|
599 |
|
double dot; |
600 |
< |
// char *token; |
601 |
< |
// char fileName[200]; |
602 |
< |
// char angleName[200]; |
277 |
< |
// char inLine[1000]; |
278 |
< |
// char inValue[200]; |
600 |
> |
|
601 |
> |
std::vector<double> tempZangs; |
602 |
> |
|
603 |
|
const char *delimit = " \t\n;,"; |
604 |
|
|
605 |
+ |
|
606 |
+ |
|
607 |
|
//open the idealCrystal.in file and zAngle.ang file |
608 |
+ |
|
609 |
|
strcpy(fileName, "idealCrystal.in"); |
610 |
+ |
|
611 |
|
strcpy(angleName, "zAngle.ang"); |
612 |
+ |
|
613 |
|
|
614 |
+ |
|
615 |
|
ifstream crystalIn(fileName); |
616 |
+ |
|
617 |
|
ifstream angleIn(angleName); |
618 |
|
|
619 |
+ |
|
620 |
+ |
|
621 |
+ |
// check to see if these files are present in the execution directory |
622 |
+ |
|
623 |
|
if (!crystalIn) { |
624 |
+ |
|
625 |
|
sprintf(painCave.errMsg, |
626 |
+ |
|
627 |
|
"Restraints Error: Unable to open idealCrystal.in for reading.\n" |
628 |
< |
"\tMake sure a reference crystal file is in the current directory.\n"); |
628 |
> |
|
629 |
> |
"\tMake sure a ref. crystal file is in the working directory.\n"); |
630 |
> |
|
631 |
> |
painCave.severity = OOPSE_ERROR; |
632 |
> |
|
633 |
|
painCave.isFatal = 1; |
634 |
+ |
|
635 |
|
simError(); |
636 |
< |
|
295 |
< |
return; |
636 |
> |
|
637 |
|
} |
638 |
|
|
639 |
+ |
|
640 |
+ |
|
641 |
+ |
// it's not fatal to lack a zAngle.ang file, it just means you're starting |
642 |
+ |
|
643 |
+ |
// from the ideal crystal state |
644 |
+ |
|
645 |
|
if (!angleIn) { |
646 |
+ |
|
647 |
|
sprintf(painCave.errMsg, |
648 |
+ |
|
649 |
|
"Restraints Warning: The lack of a zAngle.ang file is mildly\n" |
650 |
+ |
|
651 |
|
"\tunsettling... This means the simulation is starting from the\n" |
652 |
+ |
|
653 |
|
"\tidealCrystal.in reference configuration, so the omega values\n" |
654 |
< |
"\twill all be set to zero. If this is not the case, you should\n" |
655 |
< |
"\tquestion your results.\n"); |
654 |
> |
|
655 |
> |
"\twill all be set to zero. If this is not the case, the energy\n" |
656 |
> |
|
657 |
> |
"\tcalculations will be wrong.\n"); |
658 |
> |
|
659 |
> |
painCave.severity = OOPSE_WARNING; |
660 |
> |
|
661 |
|
painCave.isFatal = 0; |
662 |
+ |
|
663 |
|
simError(); |
664 |
+ |
|
665 |
|
} |
666 |
|
|
667 |
+ |
|
668 |
+ |
|
669 |
|
// A rather specific reader for OOPSE .eor files... |
670 |
+ |
|
671 |
|
// Let's read in the perfect crystal file |
672 |
+ |
|
673 |
|
crystalIn.getline(inLine,999,'\n'); |
674 |
+ |
|
675 |
+ |
// check to see if the crystal file is the same length as starting config. |
676 |
+ |
|
677 |
+ |
token = strtok(inLine,delimit); |
678 |
+ |
|
679 |
+ |
strcpy(inValue,token); |
680 |
+ |
|
681 |
+ |
idealSize = atoi(inValue); |
682 |
+ |
|
683 |
+ |
if (idealSize != vecParticles.size()) { |
684 |
+ |
|
685 |
+ |
sprintf(painCave.errMsg, |
686 |
+ |
|
687 |
+ |
"Restraints Error: Reference crystal file is not valid.\n" |
688 |
+ |
|
689 |
+ |
"\tMake sure the idealCrystal.in file is the same size as the\n" |
690 |
+ |
|
691 |
+ |
"\tstarting configuration. Using an incompatable crystal will\n" |
692 |
+ |
|
693 |
+ |
"\tlead to energy calculation failures.\n"); |
694 |
+ |
|
695 |
+ |
painCave.severity = OOPSE_ERROR; |
696 |
+ |
|
697 |
+ |
painCave.isFatal = 1; |
698 |
+ |
|
699 |
+ |
simError(); |
700 |
+ |
|
701 |
+ |
} |
702 |
+ |
|
703 |
+ |
// else, the file is okay... let's continue |
704 |
+ |
|
705 |
|
crystalIn.getline(inLine,999,'\n'); |
706 |
+ |
|
707 |
|
|
708 |
+ |
|
709 |
|
for (i=0; i<vecParticles.size(); i++) { |
710 |
+ |
|
711 |
|
crystalIn.getline(inLine,999,'\n'); |
712 |
+ |
|
713 |
|
token = strtok(inLine,delimit); |
714 |
+ |
|
715 |
|
token = strtok(NULL,delimit); |
716 |
+ |
|
717 |
|
strcpy(inValue,token); |
718 |
+ |
|
719 |
|
cofmPosX.push_back(atof(inValue)); |
720 |
+ |
|
721 |
|
token = strtok(NULL,delimit); |
722 |
+ |
|
723 |
|
strcpy(inValue,token); |
724 |
+ |
|
725 |
|
cofmPosY.push_back(atof(inValue)); |
726 |
+ |
|
727 |
|
token = strtok(NULL,delimit); |
728 |
+ |
|
729 |
|
strcpy(inValue,token); |
730 |
+ |
|
731 |
|
cofmPosZ.push_back(atof(inValue)); |
732 |
+ |
|
733 |
|
token = strtok(NULL,delimit); |
734 |
+ |
|
735 |
|
token = strtok(NULL,delimit); |
736 |
+ |
|
737 |
|
token = strtok(NULL,delimit); |
738 |
+ |
|
739 |
|
token = strtok(NULL,delimit); |
740 |
+ |
|
741 |
|
strcpy(inValue,token); |
742 |
+ |
|
743 |
|
quat0 = atof(inValue); |
744 |
+ |
|
745 |
|
token = strtok(NULL,delimit); |
746 |
+ |
|
747 |
|
strcpy(inValue,token); |
748 |
+ |
|
749 |
|
quat1 = atof(inValue); |
750 |
+ |
|
751 |
|
token = strtok(NULL,delimit); |
752 |
+ |
|
753 |
|
strcpy(inValue,token); |
754 |
+ |
|
755 |
|
quat2 = atof(inValue); |
756 |
+ |
|
757 |
|
token = strtok(NULL,delimit); |
758 |
+ |
|
759 |
|
strcpy(inValue,token); |
760 |
+ |
|
761 |
|
quat3 = atof(inValue); |
762 |
|
|
763 |
+ |
|
764 |
+ |
|
765 |
|
// now build the rotation matrix and find the unit vectors |
766 |
+ |
|
767 |
|
RfromQ[0][0] = quat0*quat0 + quat1*quat1 - quat2*quat2 - quat3*quat3; |
768 |
+ |
|
769 |
|
RfromQ[0][1] = 2*(quat1*quat2 + quat0*quat3); |
770 |
+ |
|
771 |
|
RfromQ[0][2] = 2*(quat1*quat3 - quat0*quat2); |
772 |
+ |
|
773 |
|
RfromQ[1][0] = 2*(quat1*quat2 - quat0*quat3); |
774 |
+ |
|
775 |
|
RfromQ[1][1] = quat0*quat0 - quat1*quat1 + quat2*quat2 - quat3*quat3; |
776 |
+ |
|
777 |
|
RfromQ[1][2] = 2*(quat2*quat3 + quat0*quat1); |
778 |
+ |
|
779 |
|
RfromQ[2][0] = 2*(quat1*quat3 + quat0*quat2); |
780 |
+ |
|
781 |
|
RfromQ[2][1] = 2*(quat2*quat3 - quat0*quat1); |
782 |
+ |
|
783 |
|
RfromQ[2][2] = quat0*quat0 - quat1*quat1 - quat2*quat2 + quat3*quat3; |
784 |
+ |
|
785 |
|
|
786 |
+ |
|
787 |
|
normalize = sqrt(RfromQ[2][0]*RfromQ[2][0] + RfromQ[2][1]*RfromQ[2][1] |
788 |
+ |
|
789 |
|
+ RfromQ[2][2]*RfromQ[2][2]); |
790 |
+ |
|
791 |
|
uX0.push_back(RfromQ[2][0]/normalize); |
792 |
+ |
|
793 |
|
uY0.push_back(RfromQ[2][1]/normalize); |
794 |
+ |
|
795 |
|
uZ0.push_back(RfromQ[2][2]/normalize); |
796 |
|
|
797 |
+ |
|
798 |
+ |
|
799 |
|
normalize = sqrt(RfromQ[1][0]*RfromQ[1][0] + RfromQ[1][1]*RfromQ[1][1] |
800 |
+ |
|
801 |
|
+ RfromQ[1][2]*RfromQ[1][2]); |
802 |
+ |
|
803 |
|
vX0.push_back(RfromQ[1][0]/normalize); |
804 |
+ |
|
805 |
|
vY0.push_back(RfromQ[1][1]/normalize); |
806 |
+ |
|
807 |
|
vZ0.push_back(RfromQ[1][2]/normalize); |
808 |
+ |
|
809 |
|
} |
810 |
|
|
811 |
< |
// now we can read in the zAngle.ang file |
811 |
> |
crystalIn.close(); |
812 |
> |
|
813 |
> |
|
814 |
> |
|
815 |
> |
// now we read in the zAngle.ang file |
816 |
> |
|
817 |
|
if (angleIn){ |
818 |
+ |
|
819 |
|
angleIn.getline(inLine,999,'\n'); |
820 |
< |
for (i=0; i<vecParticles.size(); i++) { |
821 |
< |
angleIn.getline(inLine,999,'\n'); |
820 |
> |
|
821 |
> |
angleIn.getline(inLine,999,'\n'); |
822 |
> |
|
823 |
> |
while (!angleIn.eof()) { |
824 |
> |
|
825 |
|
token = strtok(inLine,delimit); |
826 |
+ |
|
827 |
|
strcpy(inValue,token); |
828 |
< |
vecParticles[i]->setZangle(atof(inValue)); |
828 |
> |
|
829 |
> |
tempZangs.push_back(atof(inValue)); |
830 |
> |
|
831 |
> |
angleIn.getline(inLine,999,'\n'); |
832 |
> |
|
833 |
|
} |
834 |
+ |
|
835 |
+ |
|
836 |
+ |
|
837 |
+ |
// test to make sure the zAngle.ang file is the proper length |
838 |
+ |
|
839 |
+ |
if (tempZangs.size() == vecParticles.size()) |
840 |
+ |
|
841 |
+ |
for (i=0; i<vecParticles.size(); i++) |
842 |
+ |
|
843 |
+ |
vecParticles[i]->setZangle(tempZangs[i]); |
844 |
+ |
|
845 |
+ |
else { |
846 |
+ |
|
847 |
+ |
sprintf(painCave.errMsg, |
848 |
+ |
|
849 |
+ |
"Restraints Error: the supplied zAngle file is not valid.\n" |
850 |
+ |
|
851 |
+ |
"\tMake sure the zAngle.ang file matches with the initial\n" |
852 |
+ |
|
853 |
+ |
"\tconfiguration (i.e. they're the same length). Using the wrong\n" |
854 |
+ |
|
855 |
+ |
"\tzAngle file will lead to errors in the energy calculations.\n"); |
856 |
+ |
|
857 |
+ |
painCave.severity = OOPSE_ERROR; |
858 |
+ |
|
859 |
+ |
painCave.isFatal = 1; |
860 |
+ |
|
861 |
+ |
simError(); |
862 |
+ |
|
863 |
+ |
} |
864 |
+ |
|
865 |
|
} |
866 |
|
|
867 |
+ |
angleIn.close(); |
868 |
+ |
|
869 |
+ |
|
870 |
+ |
|
871 |
|
return; |
872 |
+ |
|
873 |
|
} |
874 |
|
|
875 |
+ |
|
876 |
+ |
|
877 |
|
void Restraints::Write_zAngle_File(vector<StuntDouble*> vecParticles){ |
878 |
|
|
879 |
+ |
|
880 |
+ |
|
881 |
|
char zOutName[200]; |
882 |
|
|
883 |
+ |
|
884 |
+ |
|
885 |
|
strcpy(zOutName,"zAngle.ang"); |
886 |
|
|
887 |
+ |
|
888 |
+ |
|
889 |
|
ofstream angleOut(zOutName); |
890 |
+ |
|
891 |
|
angleOut << "This file contains the omega values for the .eor file\n"; |
892 |
+ |
|
893 |
|
for (i=0; i<vecParticles.size(); i++) { |
894 |
+ |
|
895 |
|
angleOut << vecParticles[i]->getZangle() << "\n"; |
896 |
+ |
|
897 |
|
} |
898 |
+ |
|
899 |
|
return; |
900 |
+ |
|
901 |
|
} |
902 |
|
|
903 |
+ |
|
904 |
+ |
|
905 |
|
double Restraints::getVharm(){ |
906 |
+ |
|
907 |
|
return harmPotent; |
908 |
+ |
|
909 |
|
} |
910 |
|
|
911 |
+ |
|
912 |
+ |
|