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