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