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 |
|
24 |
|
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 |
|
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 |
|
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 |
|
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 |
|
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 |
|
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 |
|
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 |
|
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 |
|
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 |
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 |
|
821 |
angleIn.getline(inLine,999,'\n'); |
822 |
|
823 |
while (!angleIn.eof()) { |
824 |
|
825 |
token = strtok(inLine,delimit); |
826 |
|
827 |
strcpy(inValue,token); |
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 |
|