1 |
/* |
2 |
* Copyright (c) 2005 The University of Notre Dame. All Rights Reserved. |
3 |
* |
4 |
* The University of Notre Dame grants you ("Licensee") a |
5 |
* non-exclusive, royalty free, license to use, modify and |
6 |
* redistribute this software in source and binary code form, provided |
7 |
* that the following conditions are met: |
8 |
* |
9 |
* 1. Acknowledgement of the program authors must be made in any |
10 |
* publication of scientific results based in part on use of the |
11 |
* program. An acceptable form of acknowledgement is citation of |
12 |
* the article in which the program was described (Matthew |
13 |
* A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher |
14 |
* J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented |
15 |
* Parallel Simulation Engine for Molecular Dynamics," |
16 |
* J. Comput. Chem. 26, pp. 252-271 (2005)) |
17 |
* |
18 |
* 2. Redistributions of source code must retain the above copyright |
19 |
* notice, this list of conditions and the following disclaimer. |
20 |
* |
21 |
* 3. Redistributions in binary form must reproduce the above copyright |
22 |
* notice, this list of conditions and the following disclaimer in the |
23 |
* documentation and/or other materials provided with the |
24 |
* distribution. |
25 |
* |
26 |
* This software is provided "AS IS," without a warranty of any |
27 |
* kind. All express or implied conditions, representations and |
28 |
* warranties, including any implied warranty of merchantability, |
29 |
* fitness for a particular purpose or non-infringement, are hereby |
30 |
* excluded. The University of Notre Dame and its licensors shall not |
31 |
* be liable for any damages suffered by licensee as a result of |
32 |
* using, modifying or distributing the software or its |
33 |
* derivatives. In no event will the University of Notre Dame or its |
34 |
* licensors be liable for any lost revenue, profit or data, or for |
35 |
* direct, indirect, special, consequential, incidental or punitive |
36 |
* damages, however caused and regardless of the theory of liability, |
37 |
* arising out of the use of or inability to use software, even if the |
38 |
* University of Notre Dame has been advised of the possibility of |
39 |
* such damages. |
40 |
*/ |
41 |
|
42 |
#include <stdlib.h> |
43 |
#include <math.h> |
44 |
|
45 |
using namespace std; |
46 |
|
47 |
#include "restraints/Restraints.hpp" |
48 |
#include "primitives/Molecule.hpp" |
49 |
#include "utils/simError.h" |
50 |
|
51 |
#define PI 3.14159265359 |
52 |
#define TWO_PI 6.28318530718 |
53 |
|
54 |
namespace oopse { |
55 |
|
56 |
Restraints::Restraints(SimInfo* info, double lambdaVal, double lambdaExp){ |
57 |
info_ = info; |
58 |
Globals* simParam = info_->getSimParams(); |
59 |
|
60 |
lambdaValue = lambdaVal; |
61 |
lambdaK = lambdaExp; |
62 |
|
63 |
if (simParam->getUseSolidThermInt()) { |
64 |
if (simParam->haveDistSpringConst()) { |
65 |
kDist = simParam->getDistSpringConst(); |
66 |
} |
67 |
else{ |
68 |
kDist = 6.0; |
69 |
sprintf(painCave.errMsg, |
70 |
"ThermoIntegration Warning: the spring constant for the\n" |
71 |
"\ttranslational restraint was not specified. OOPSE will use\n" |
72 |
"\ta default value of %f. To set it to something else, use\n" |
73 |
"\tthe thermIntDistSpringConst variable.\n", |
74 |
kDist); |
75 |
painCave.isFatal = 0; |
76 |
simError(); |
77 |
} |
78 |
if (simParam->haveThetaSpringConst()) { |
79 |
kTheta = simParam->getThetaSpringConst(); |
80 |
} |
81 |
else{ |
82 |
kTheta = 7.5; |
83 |
sprintf(painCave.errMsg, |
84 |
"ThermoIntegration Warning: the spring constant for the\n" |
85 |
"\tdeflection orientational restraint was not specified.\n" |
86 |
"\tOOPSE will use a default value of %f. To set it to\n" |
87 |
"\tsomething else, use the thermIntThetaSpringConst variable.\n", |
88 |
kTheta); |
89 |
painCave.isFatal = 0; |
90 |
simError(); |
91 |
} |
92 |
if (simParam->haveOmegaSpringConst()) { |
93 |
kOmega = simParam->getOmegaSpringConst(); |
94 |
} |
95 |
else{ |
96 |
kOmega = 13.5; |
97 |
sprintf(painCave.errMsg, |
98 |
"ThermoIntegration Warning: the spring constant for the\n" |
99 |
"\tspin orientational restraint was not specified. OOPSE\n" |
100 |
"\twill use a default value of %f. To set it to something\n" |
101 |
"\telse, use the thermIntOmegaSpringConst variable.\n", |
102 |
kOmega); |
103 |
painCave.isFatal = 0; |
104 |
simError(); |
105 |
} |
106 |
} |
107 |
|
108 |
// build a RestReader and read in important information |
109 |
|
110 |
restRead_ = new RestReader(info_); |
111 |
restRead_->readIdealCrystal(); |
112 |
restRead_->readZangle(); |
113 |
|
114 |
delete restRead_; |
115 |
restRead_ = NULL; |
116 |
|
117 |
} |
118 |
|
119 |
Restraints::~Restraints(){ |
120 |
} |
121 |
|
122 |
void Restraints::Calc_rVal(Vector3d &position, double refPosition[3]){ |
123 |
delRx = position.x() - refPosition[0]; |
124 |
delRy = position.y() - refPosition[1]; |
125 |
delRz = position.z() - refPosition[2]; |
126 |
|
127 |
return; |
128 |
} |
129 |
|
130 |
void Restraints::Calc_body_thetaVal(RotMat3x3d &matrix, double refUnit[3]){ |
131 |
ub0x = matrix(0,0)*refUnit[0] + matrix(0,1)*refUnit[1] |
132 |
+ matrix(0,2)*refUnit[2]; |
133 |
ub0y = matrix(1,0)*refUnit[0] + matrix(1,1)*refUnit[1] |
134 |
+ matrix(1,2)*refUnit[2]; |
135 |
ub0z = matrix(2,0)*refUnit[0] + matrix(2,1)*refUnit[1] |
136 |
+ matrix(2,2)*refUnit[2]; |
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 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(){ |
189 |
SimInfo::MoleculeIterator mi; |
190 |
Molecule* mol; |
191 |
Molecule::IntegrableObjectIterator ii; |
192 |
StuntDouble* integrableObject; |
193 |
Vector3d pos; |
194 |
RotMat3x3d A; |
195 |
double refPos[3]; |
196 |
double refVec[3]; |
197 |
double tolerance; |
198 |
double tempPotent; |
199 |
double factor; |
200 |
double spaceTrq[3]; |
201 |
double omegaPass; |
202 |
GenericData* data; |
203 |
DoubleGenericData* doubleData; |
204 |
|
205 |
tolerance = 5.72957795131e-7; |
206 |
|
207 |
harmPotent = 0.0; // zero out the global harmonic potential variable |
208 |
|
209 |
factor = 1 - pow(lambdaValue, lambdaK); |
210 |
|
211 |
for (mol = info_->beginMolecule(mi); mol != NULL; |
212 |
mol = info_->nextMolecule(mi)) { |
213 |
for (integrableObject = mol->beginIntegrableObject(ii); |
214 |
integrableObject != NULL; |
215 |
integrableObject = mol->nextIntegrableObject(ii)) { |
216 |
|
217 |
// obtain the current and reference positions |
218 |
pos = integrableObject->getPos(); |
219 |
|
220 |
data = integrableObject->getPropertyByName("refPosX"); |
221 |
if (data){ |
222 |
doubleData = dynamic_cast<DoubleGenericData*>(data); |
223 |
if (!doubleData){ |
224 |
cerr << "Can't obtain refPosX from StuntDouble\n"; |
225 |
return 0.0; |
226 |
} |
227 |
else refPos[0] = doubleData->getData(); |
228 |
} |
229 |
data = integrableObject->getPropertyByName("refPosY"); |
230 |
if (data){ |
231 |
doubleData = dynamic_cast<DoubleGenericData*>(data); |
232 |
if (!doubleData){ |
233 |
cerr << "Can't obtain refPosY from StuntDouble\n"; |
234 |
return 0.0; |
235 |
} |
236 |
else refPos[1] = doubleData->getData(); |
237 |
} |
238 |
data = integrableObject->getPropertyByName("refPosZ"); |
239 |
if (data){ |
240 |
doubleData = dynamic_cast<DoubleGenericData*>(data); |
241 |
if (!doubleData){ |
242 |
cerr << "Can't obtain refPosZ from StuntDouble\n"; |
243 |
return 0.0; |
244 |
} |
245 |
else refPos[2] = doubleData->getData(); |
246 |
} |
247 |
|
248 |
// calculate the displacement |
249 |
Calc_rVal( pos, refPos ); |
250 |
|
251 |
// calculate the derivatives |
252 |
dVdrx = -kDist*delRx; |
253 |
dVdry = -kDist*delRy; |
254 |
dVdrz = -kDist*delRz; |
255 |
|
256 |
// next we calculate the restraint forces |
257 |
restraintFrc[0] = dVdrx; |
258 |
restraintFrc[1] = dVdry; |
259 |
restraintFrc[2] = dVdrz; |
260 |
tempPotent = 0.5*kDist*(delRx*delRx + delRy*delRy + delRz*delRz); |
261 |
|
262 |
// apply the lambda scaling factor to the forces |
263 |
for (j = 0; j < 3; j++) restraintFrc[j] *= factor; |
264 |
|
265 |
// and add the temporary force to the total force |
266 |
integrableObject->addFrc(restraintFrc); |
267 |
|
268 |
// if the particle is directional, we accumulate the rot. restraints |
269 |
if (integrableObject->isDirectional()){ |
270 |
|
271 |
// get the current rotation matrix and reference vector |
272 |
A = integrableObject->getA(); |
273 |
|
274 |
data = integrableObject->getPropertyByName("refVectorX"); |
275 |
if (data){ |
276 |
doubleData = dynamic_cast<DoubleGenericData*>(data); |
277 |
if (!doubleData){ |
278 |
cerr << "Can't obtain refVectorX from StuntDouble\n"; |
279 |
return 0.0; |
280 |
} |
281 |
else refVec[0] = doubleData->getData(); |
282 |
} |
283 |
data = integrableObject->getPropertyByName("refVectorY"); |
284 |
if (data){ |
285 |
doubleData = dynamic_cast<DoubleGenericData*>(data); |
286 |
if (!doubleData){ |
287 |
cerr << "Can't obtain refVectorY from StuntDouble\n"; |
288 |
return 0.0; |
289 |
} |
290 |
else refVec[1] = doubleData->getData(); |
291 |
} |
292 |
data = integrableObject->getPropertyByName("refVectorZ"); |
293 |
if (data){ |
294 |
doubleData = dynamic_cast<DoubleGenericData*>(data); |
295 |
if (!doubleData){ |
296 |
cerr << "Can't obtain refVectorZ from StuntDouble\n"; |
297 |
return 0.0; |
298 |
} |
299 |
else refVec[2] = doubleData->getData(); |
300 |
} |
301 |
|
302 |
// calculate the theta and omega displacements |
303 |
Calc_body_thetaVal( A, refVec ); |
304 |
omegaPass = integrableObject->getZangle(); |
305 |
Calc_body_omegaVal( omegaPass ); |
306 |
|
307 |
// uTx... and vTx... are the body-fixed z and y unit vectors |
308 |
uTx = 0.0; |
309 |
uTy = 0.0; |
310 |
uTz = 1.0; |
311 |
vTx = 0.0; |
312 |
vTy = 1.0; |
313 |
vTz = 0.0; |
314 |
|
315 |
dVdux = 0.0; |
316 |
dVduy = 0.0; |
317 |
dVduz = 0.0; |
318 |
dVdvx = 0.0; |
319 |
dVdvy = 0.0; |
320 |
dVdvz = 0.0; |
321 |
|
322 |
if (fabs(theta) > tolerance) { |
323 |
dVdux = -(kTheta*theta/sin(theta))*ub0x; |
324 |
dVduy = -(kTheta*theta/sin(theta))*ub0y; |
325 |
dVduz = -(kTheta*theta/sin(theta))*ub0z; |
326 |
} |
327 |
|
328 |
if (fabs(omega) > tolerance) { |
329 |
dVdvx = -(kOmega*omega/sin(omega))*vb0x; |
330 |
dVdvy = -(kOmega*omega/sin(omega))*vb0y; |
331 |
dVdvz = -(kOmega*omega/sin(omega))*vb0z; |
332 |
} |
333 |
|
334 |
// next we calculate the restraint torques |
335 |
restraintTrq[0] = 0.0; |
336 |
restraintTrq[1] = 0.0; |
337 |
restraintTrq[2] = 0.0; |
338 |
|
339 |
if (fabs(omega) > tolerance) { |
340 |
restraintTrq[0] += 0.0; |
341 |
restraintTrq[1] += 0.0; |
342 |
restraintTrq[2] += vTy*dVdvx; |
343 |
tempPotent += 0.5*(kOmega*omega*omega); |
344 |
} |
345 |
if (fabs(theta) > tolerance) { |
346 |
restraintTrq[0] += (uTz*dVduy); |
347 |
restraintTrq[1] += -(uTz*dVdux); |
348 |
restraintTrq[2] += 0.0; |
349 |
tempPotent += 0.5*(kTheta*theta*theta); |
350 |
} |
351 |
|
352 |
// apply the lambda scaling factor to these torques |
353 |
for (j = 0; j < 3; j++) restraintTrq[j] *= factor; |
354 |
|
355 |
// now we need to convert from body-fixed to space-fixed torques |
356 |
spaceTrq[0] = A(0,0)*restraintTrq[0] + A(1,0)*restraintTrq[1] |
357 |
+ A(2,0)*restraintTrq[2]; |
358 |
spaceTrq[1] = A(0,1)*restraintTrq[0] + A(1,1)*restraintTrq[1] |
359 |
+ A(2,1)*restraintTrq[2]; |
360 |
spaceTrq[2] = A(0,2)*restraintTrq[0] + A(1,2)*restraintTrq[1] |
361 |
+ A(2,2)*restraintTrq[2]; |
362 |
|
363 |
// now pass this temporary torque vector to the total torque |
364 |
integrableObject->addTrq(spaceTrq); |
365 |
} |
366 |
|
367 |
// update the total harmonic potential with this object's contribution |
368 |
harmPotent += tempPotent; |
369 |
} |
370 |
|
371 |
} |
372 |
|
373 |
// we can finish by returning the appropriately scaled potential energy |
374 |
tempPotent = harmPotent * factor; |
375 |
|
376 |
return tempPotent; |
377 |
|
378 |
} |
379 |
|
380 |
}// end namespace oopse |