1 |
// Thermodynamic integration is not multiprocessor friendly right now |
2 |
|
3 |
#include <iostream> |
4 |
#include <stdlib.h> |
5 |
#include <cstdio> |
6 |
#include <fstream> |
7 |
#include <iomanip> |
8 |
#include <string> |
9 |
#include <cstring> |
10 |
#include <math.h> |
11 |
|
12 |
using namespace std; |
13 |
|
14 |
#include "restraints/Restraints.hpp" |
15 |
#include "brains/SimInfo.hpp" |
16 |
#include "utils/simError.h" |
17 |
#include "io/basic_ifstrstream.hpp" |
18 |
|
19 |
#define PI 3.14159265359 |
20 |
#define TWO_PI 6.28318530718 |
21 |
|
22 |
Restraints::Restraints(double lambdaVal, double lambdaExp){ |
23 |
lambdaValue = lambdaVal; |
24 |
lambdaK = lambdaExp; |
25 |
vector<double> resConsts; |
26 |
const char *jolt = " \t\n;,"; |
27 |
|
28 |
#ifdef IS_MPI |
29 |
if(worldRank == 0 ){ |
30 |
#endif // is_mpi |
31 |
|
32 |
strcpy(springName, "HarmSpringConsts.txt"); |
33 |
|
34 |
ifstream springs(springName); |
35 |
|
36 |
if (!springs) { |
37 |
sprintf(painCave.errMsg, |
38 |
"Unable to open HarmSpringConsts.txt for reading, so the\n" |
39 |
"\tdefault spring constants will be loaded. If you want\n" |
40 |
"\tto specify spring constants, include a three line\n" |
41 |
"\tHarmSpringConsts.txt file in the execution directory.\n"); |
42 |
painCave.severity = OOPSE_WARNING; |
43 |
painCave.isFatal = 0; |
44 |
simError(); |
45 |
|
46 |
// load default spring constants |
47 |
kDist = 6; // spring constant in units of kcal/(mol*ang^2) |
48 |
kTheta = 7.5; // in units of kcal/mol |
49 |
kOmega = 13.5; // in units of kcal/mol |
50 |
} else { |
51 |
|
52 |
springs.getline(inLine,999,'\n'); |
53 |
// the file is blank! |
54 |
if (springs.eof()){ |
55 |
sprintf(painCave.errMsg, |
56 |
"HarmSpringConsts.txt file is not valid.\n" |
57 |
"\tThe file should contain four rows, the last three containing\n" |
58 |
"\ta label and the spring constant value. They should be listed\n" |
59 |
"\tin the following order: kDist (positional restrant), kTheta\n" |
60 |
"\t(rot. restraint: deflection of z-axis), and kOmega (rot.\n" |
61 |
"\trestraint: rotation about the z-axis).\n"); |
62 |
painCave.severity = OOPSE_ERROR; |
63 |
painCave.isFatal = 1; |
64 |
simError(); |
65 |
} |
66 |
// read in spring constants and check to make sure it is a valid file |
67 |
springs.getline(inLine,999,'\n'); |
68 |
while (!springs.eof()){ |
69 |
if (NULL != inLine){ |
70 |
token = strtok(inLine,jolt); |
71 |
token = strtok(NULL,jolt); |
72 |
if (NULL != token){ |
73 |
strcpy(inValue,token); |
74 |
resConsts.push_back(atof(inValue)); |
75 |
} |
76 |
} |
77 |
springs.getline(inLine,999,'\n'); |
78 |
} |
79 |
if (resConsts.size() == 3){ |
80 |
kDist = resConsts[0]; |
81 |
kTheta = resConsts[1]; |
82 |
kOmega = resConsts[2]; |
83 |
} |
84 |
else { |
85 |
sprintf(painCave.errMsg, |
86 |
"HarmSpringConsts.txt file is not valid.\n" |
87 |
"\tThe file should contain four rows, the last three containing\n" |
88 |
"\ta label and the spring constant value. They should be listed\n" |
89 |
"\tin the following order: kDist (positional restrant), kTheta\n" |
90 |
"\t(rot. restraint: deflection of z-axis), and kOmega (rot.\n" |
91 |
"\trestraint: rotation about the z-axis).\n"); |
92 |
painCave.severity = OOPSE_ERROR; |
93 |
painCave.isFatal = 1; |
94 |
simError(); |
95 |
} |
96 |
} |
97 |
#ifdef IS_MPI |
98 |
} |
99 |
|
100 |
MPI_Bcast(&kDist, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); |
101 |
MPI_Bcast(&kTheta, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); |
102 |
MPI_Bcast(&kOmega, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); |
103 |
|
104 |
sprintf( checkPointMsg, |
105 |
"Sucessfully opened and read spring file.\n"); |
106 |
MPIcheckPoint(); |
107 |
|
108 |
#endif // is_mpi |
109 |
|
110 |
sprintf(painCave.errMsg, |
111 |
"The spring constants for thermodynamic integration are:\n" |
112 |
"\tkDist = %lf\n" |
113 |
"\tkTheta = %lf\n" |
114 |
"\tkOmega = %lf\n", kDist, kTheta, kOmega); |
115 |
painCave.severity = OOPSE_INFO; |
116 |
painCave.isFatal = 0; |
117 |
simError(); |
118 |
} |
119 |
|
120 |
Restraints::~Restraints(){ |
121 |
} |
122 |
|
123 |
void Restraints::Calc_rVal(double position[3], int currentMol){ |
124 |
delRx = position[0] - cofmPosX[currentMol]; |
125 |
delRy = position[1] - cofmPosY[currentMol]; |
126 |
delRz = position[2] - cofmPosZ[currentMol]; |
127 |
|
128 |
return; |
129 |
} |
130 |
|
131 |
void Restraints::Calc_body_thetaVal(double matrix[3][3], int currentMol){ |
132 |
ub0x = matrix[0][0]*uX0[currentMol] + matrix[0][1]*uY0[currentMol] |
133 |
+ matrix[0][2]*uZ0[currentMol]; |
134 |
ub0y = matrix[1][0]*uX0[currentMol] + matrix[1][1]*uY0[currentMol] |
135 |
+ matrix[1][2]*uZ0[currentMol]; |
136 |
ub0z = matrix[2][0]*uX0[currentMol] + matrix[2][1]*uY0[currentMol] |
137 |
+ matrix[2][2]*uZ0[currentMol]; |
138 |
|
139 |
normalize = sqrt(ub0x*ub0x + ub0y*ub0y + ub0z*ub0z); |
140 |
ub0x = ub0x/normalize; |
141 |
ub0y = ub0y/normalize; |
142 |
ub0z = ub0z/normalize; |
143 |
|
144 |
// Theta is the dot product of the reference and new z-axes |
145 |
theta = acos(ub0z); |
146 |
|
147 |
return; |
148 |
} |
149 |
|
150 |
void Restraints::Calc_body_omegaVal(double matrix[3][3], double zAngle){ |
151 |
double zRotator[3][3]; |
152 |
double tempOmega; |
153 |
double wholeTwoPis; |
154 |
// Use the omega accumulated from the rotation propagation |
155 |
omega = zAngle; |
156 |
|
157 |
// translate the omega into a range between -PI and PI |
158 |
if (omega < -PI){ |
159 |
tempOmega = omega / -TWO_PI; |
160 |
wholeTwoPis = floor(tempOmega); |
161 |
tempOmega = omega + TWO_PI*wholeTwoPis; |
162 |
if (tempOmega < -PI) |
163 |
omega = tempOmega + TWO_PI; |
164 |
else |
165 |
omega = tempOmega; |
166 |
} |
167 |
if (omega > PI){ |
168 |
tempOmega = omega / TWO_PI; |
169 |
wholeTwoPis = floor(tempOmega); |
170 |
tempOmega = omega - TWO_PI*wholeTwoPis; |
171 |
if (tempOmega > PI) |
172 |
omega = tempOmega - TWO_PI; |
173 |
else |
174 |
omega = tempOmega; |
175 |
} |
176 |
|
177 |
vb0x = sin(omega); |
178 |
vb0y = cos(omega); |
179 |
vb0z = 0.0; |
180 |
|
181 |
normalize = sqrt(vb0x*vb0x + vb0y*vb0y + vb0z*vb0z); |
182 |
vb0x = vb0x/normalize; |
183 |
vb0y = vb0y/normalize; |
184 |
vb0z = vb0z/normalize; |
185 |
|
186 |
return; |
187 |
} |
188 |
|
189 |
double Restraints::Calc_Restraint_Forces(vector<StuntDouble*> vecParticles){ |
190 |
double pos[3]; |
191 |
double A[3][3]; |
192 |
double tolerance; |
193 |
double tempPotent; |
194 |
double factor; |
195 |
double spaceTrq[3]; |
196 |
double omegaPass; |
197 |
|
198 |
tolerance = 5.72957795131e-7; |
199 |
|
200 |
harmPotent = 0.0; // zero out the global harmonic potential variable |
201 |
|
202 |
factor = 1 - pow(lambdaValue, lambdaK); |
203 |
|
204 |
for (i=0; i<vecParticles.size(); i++){ |
205 |
if (vecParticles[i]->isDirectional()){ |
206 |
vecParticles[i]->getPos(pos); |
207 |
vecParticles[i]->getA(A); |
208 |
Calc_rVal( pos, i ); |
209 |
Calc_body_thetaVal( A, i ); |
210 |
omegaPass = vecParticles[i]->getZangle(); |
211 |
Calc_body_omegaVal( A, omegaPass ); |
212 |
|
213 |
// first we calculate the derivatives |
214 |
dVdrx = -kDist*delRx; |
215 |
dVdry = -kDist*delRy; |
216 |
dVdrz = -kDist*delRz; |
217 |
|
218 |
// uTx... and vTx... are the body-fixed z and y unit vectors |
219 |
uTx = 0.0; |
220 |
uTy = 0.0; |
221 |
uTz = 1.0; |
222 |
vTx = 0.0; |
223 |
vTy = 1.0; |
224 |
vTz = 0.0; |
225 |
|
226 |
dVdux = 0; |
227 |
dVduy = 0; |
228 |
dVduz = 0; |
229 |
dVdvx = 0; |
230 |
dVdvy = 0; |
231 |
dVdvz = 0; |
232 |
|
233 |
if (fabs(theta) > tolerance) { |
234 |
dVdux = -(kTheta*theta/sin(theta))*ub0x; |
235 |
dVduy = -(kTheta*theta/sin(theta))*ub0y; |
236 |
dVduz = -(kTheta*theta/sin(theta))*ub0z; |
237 |
} |
238 |
|
239 |
if (fabs(omega) > tolerance) { |
240 |
dVdvx = -(kOmega*omega/sin(omega))*vb0x; |
241 |
dVdvy = -(kOmega*omega/sin(omega))*vb0y; |
242 |
dVdvz = -(kOmega*omega/sin(omega))*vb0z; |
243 |
} |
244 |
|
245 |
// next we calculate the restraint forces and torques |
246 |
restraintFrc[0] = dVdrx; |
247 |
restraintFrc[1] = dVdry; |
248 |
restraintFrc[2] = dVdrz; |
249 |
tempPotent = 0.5*kDist*(delRx*delRx + delRy*delRy + delRz*delRz); |
250 |
|
251 |
restraintTrq[0] = 0.0; |
252 |
restraintTrq[1] = 0.0; |
253 |
restraintTrq[2] = 0.0; |
254 |
|
255 |
if (fabs(omega) > tolerance) { |
256 |
restraintTrq[0] += 0.0; |
257 |
restraintTrq[1] += 0.0; |
258 |
restraintTrq[2] += vTy*dVdvx; |
259 |
tempPotent += 0.5*(kOmega*omega*omega); |
260 |
} |
261 |
if (fabs(theta) > tolerance) { |
262 |
restraintTrq[0] += (uTz*dVduy); |
263 |
restraintTrq[1] += -(uTz*dVdux); |
264 |
restraintTrq[2] += 0.0; |
265 |
tempPotent += 0.5*(kTheta*theta*theta); |
266 |
} |
267 |
|
268 |
for (j = 0; j < 3; j++) { |
269 |
restraintFrc[j] *= factor; |
270 |
restraintTrq[j] *= factor; |
271 |
} |
272 |
|
273 |
harmPotent += tempPotent; |
274 |
|
275 |
// now we need to convert from body-fixed torques to space-fixed torques |
276 |
spaceTrq[0] = A[0][0]*restraintTrq[0] + A[1][0]*restraintTrq[1] |
277 |
+ A[2][0]*restraintTrq[2]; |
278 |
spaceTrq[1] = A[0][1]*restraintTrq[0] + A[1][1]*restraintTrq[1] |
279 |
+ A[2][1]*restraintTrq[2]; |
280 |
spaceTrq[2] = A[0][2]*restraintTrq[0] + A[1][2]*restraintTrq[1] |
281 |
+ A[2][2]*restraintTrq[2]; |
282 |
|
283 |
// now it's time to pass these temporary forces and torques |
284 |
// to the total forces and torques |
285 |
vecParticles[i]->addFrc(restraintFrc); |
286 |
vecParticles[i]->addTrq(spaceTrq); |
287 |
} |
288 |
} |
289 |
|
290 |
// and we can return the appropriately scaled potential energy |
291 |
tempPotent = harmPotent * factor; |
292 |
return tempPotent; |
293 |
} |
294 |
|
295 |
void Restraints::Store_Init_Info(vector<StuntDouble*> vecParticles){ |
296 |
int idealSize; |
297 |
double pos[3]; |
298 |
double A[3][3]; |
299 |
double RfromQ[3][3]; |
300 |
double quat0, quat1, quat2, quat3; |
301 |
double dot; |
302 |
vector<double> tempZangs; |
303 |
const char *delimit = " \t\n;,"; |
304 |
|
305 |
//open the idealCrystal.in file and zAngle.ang file |
306 |
strcpy(fileName, "idealCrystal.in"); |
307 |
strcpy(angleName, "zAngle.ang"); |
308 |
|
309 |
ifstrstream crystalIn(fileName); |
310 |
ifstrstream angleIn(angleName); |
311 |
|
312 |
// check to see if these files are present in the execution directory |
313 |
if (!crystalIn) { |
314 |
sprintf(painCave.errMsg, |
315 |
"Restraints Error: Unable to open idealCrystal.in for reading.\n" |
316 |
"\tMake sure a ref. crystal file is in the working directory.\n"); |
317 |
painCave.severity = OOPSE_ERROR; |
318 |
painCave.isFatal = 1; |
319 |
simError(); |
320 |
} |
321 |
|
322 |
// it's not fatal to lack a zAngle.ang file, it just means you're starting |
323 |
// from the ideal crystal state |
324 |
if (!angleIn) { |
325 |
sprintf(painCave.errMsg, |
326 |
"Restraints Warning: The lack of a zAngle.ang file is mildly\n" |
327 |
"\tunsettling... This means the simulation is starting from the\n" |
328 |
"\tidealCrystal.in reference configuration, so the omega values\n" |
329 |
"\twill all be set to zero. If this is not the case, the energy\n" |
330 |
"\tcalculations will be wrong.\n"); |
331 |
painCave.severity = OOPSE_WARNING; |
332 |
painCave.isFatal = 0; |
333 |
simError(); |
334 |
} |
335 |
|
336 |
// A rather specific reader for OOPSE .eor files... |
337 |
// Let's read in the perfect crystal file |
338 |
crystalIn.getline(inLine,999,'\n'); |
339 |
// check to see if the crystal file is the same length as starting config. |
340 |
token = strtok(inLine,delimit); |
341 |
strcpy(inValue,token); |
342 |
idealSize = atoi(inValue); |
343 |
if (idealSize != vecParticles.size()) { |
344 |
sprintf(painCave.errMsg, |
345 |
"Restraints Error: Reference crystal file is not valid.\n" |
346 |
"\tMake sure the idealCrystal.in file is the same size as the\n" |
347 |
"\tstarting configuration. Using an incompatable crystal will\n" |
348 |
"\tlead to energy calculation failures.\n"); |
349 |
painCave.severity = OOPSE_ERROR; |
350 |
painCave.isFatal = 1; |
351 |
simError(); |
352 |
} |
353 |
// else, the file is okay... let's continue |
354 |
crystalIn.getline(inLine,999,'\n'); |
355 |
|
356 |
for (i=0; i<vecParticles.size(); i++) { |
357 |
crystalIn.getline(inLine,999,'\n'); |
358 |
token = strtok(inLine,delimit); |
359 |
token = strtok(NULL,delimit); |
360 |
strcpy(inValue,token); |
361 |
cofmPosX.push_back(atof(inValue)); |
362 |
token = strtok(NULL,delimit); |
363 |
strcpy(inValue,token); |
364 |
cofmPosY.push_back(atof(inValue)); |
365 |
token = strtok(NULL,delimit); |
366 |
strcpy(inValue,token); |
367 |
cofmPosZ.push_back(atof(inValue)); |
368 |
token = strtok(NULL,delimit); |
369 |
token = strtok(NULL,delimit); |
370 |
token = strtok(NULL,delimit); |
371 |
token = strtok(NULL,delimit); |
372 |
strcpy(inValue,token); |
373 |
quat0 = atof(inValue); |
374 |
token = strtok(NULL,delimit); |
375 |
strcpy(inValue,token); |
376 |
quat1 = atof(inValue); |
377 |
token = strtok(NULL,delimit); |
378 |
strcpy(inValue,token); |
379 |
quat2 = atof(inValue); |
380 |
token = strtok(NULL,delimit); |
381 |
strcpy(inValue,token); |
382 |
quat3 = atof(inValue); |
383 |
|
384 |
// now build the rotation matrix and find the unit vectors |
385 |
RfromQ[0][0] = quat0*quat0 + quat1*quat1 - quat2*quat2 - quat3*quat3; |
386 |
RfromQ[0][1] = 2*(quat1*quat2 + quat0*quat3); |
387 |
RfromQ[0][2] = 2*(quat1*quat3 - quat0*quat2); |
388 |
RfromQ[1][0] = 2*(quat1*quat2 - quat0*quat3); |
389 |
RfromQ[1][1] = quat0*quat0 - quat1*quat1 + quat2*quat2 - quat3*quat3; |
390 |
RfromQ[1][2] = 2*(quat2*quat3 + quat0*quat1); |
391 |
RfromQ[2][0] = 2*(quat1*quat3 + quat0*quat2); |
392 |
RfromQ[2][1] = 2*(quat2*quat3 - quat0*quat1); |
393 |
RfromQ[2][2] = quat0*quat0 - quat1*quat1 - quat2*quat2 + quat3*quat3; |
394 |
|
395 |
normalize = sqrt(RfromQ[2][0]*RfromQ[2][0] + RfromQ[2][1]*RfromQ[2][1] |
396 |
+ RfromQ[2][2]*RfromQ[2][2]); |
397 |
uX0.push_back(RfromQ[2][0]/normalize); |
398 |
uY0.push_back(RfromQ[2][1]/normalize); |
399 |
uZ0.push_back(RfromQ[2][2]/normalize); |
400 |
|
401 |
normalize = sqrt(RfromQ[1][0]*RfromQ[1][0] + RfromQ[1][1]*RfromQ[1][1] |
402 |
+ RfromQ[1][2]*RfromQ[1][2]); |
403 |
vX0.push_back(RfromQ[1][0]/normalize); |
404 |
vY0.push_back(RfromQ[1][1]/normalize); |
405 |
vZ0.push_back(RfromQ[1][2]/normalize); |
406 |
} |
407 |
crystalIn.close(); |
408 |
|
409 |
// now we read in the zAngle.ang file |
410 |
if (angleIn){ |
411 |
angleIn.getline(inLine,999,'\n'); |
412 |
angleIn.getline(inLine,999,'\n'); |
413 |
while (!angleIn.eof()) { |
414 |
token = strtok(inLine,delimit); |
415 |
strcpy(inValue,token); |
416 |
tempZangs.push_back(atof(inValue)); |
417 |
angleIn.getline(inLine,999,'\n'); |
418 |
} |
419 |
|
420 |
// test to make sure the zAngle.ang file is the proper length |
421 |
if (tempZangs.size() == vecParticles.size()) |
422 |
for (i=0; i<vecParticles.size(); i++) |
423 |
vecParticles[i]->setZangle(tempZangs[i]); |
424 |
else { |
425 |
sprintf(painCave.errMsg, |
426 |
"Restraints Error: the supplied zAngle file is not valid.\n" |
427 |
"\tMake sure the zAngle.ang file matches with the initial\n" |
428 |
"\tconfiguration (i.e. they're the same length). Using the wrong\n" |
429 |
"\tzAngle file will lead to errors in the energy calculations.\n"); |
430 |
painCave.severity = OOPSE_ERROR; |
431 |
painCave.isFatal = 1; |
432 |
simError(); |
433 |
} |
434 |
} |
435 |
angleIn.close(); |
436 |
|
437 |
return; |
438 |
} |
439 |
|
440 |
void Restraints::Write_zAngle_File(vector<StuntDouble*> vecParticles){ |
441 |
|
442 |
char zOutName[200]; |
443 |
|
444 |
strcpy(zOutName,"zAngle.ang"); |
445 |
|
446 |
ofstream angleOut(zOutName); |
447 |
angleOut << "This file contains the omega values for the .eor file\n"; |
448 |
for (i=0; i<vecParticles.size(); i++) { |
449 |
angleOut << vecParticles[i]->getZangle() << "\n"; |
450 |
} |
451 |
return; |
452 |
} |
453 |
|
454 |
double Restraints::getVharm(){ |
455 |
return harmPotent; |
456 |
} |
457 |
|