1 |
/* |
2 |
* Copyright (c) 2009 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. Redistributions of source code must retain the above copyright |
10 |
* notice, this list of conditions and the following disclaimer. |
11 |
* |
12 |
* 2. Redistributions in binary form must reproduce the above copyright |
13 |
* notice, this list of conditions and the following disclaimer in the |
14 |
* documentation and/or other materials provided with the |
15 |
* distribution. |
16 |
* |
17 |
* This software is provided "AS IS," without a warranty of any |
18 |
* kind. All express or implied conditions, representations and |
19 |
* warranties, including any implied warranty of merchantability, |
20 |
* fitness for a particular purpose or non-infringement, are hereby |
21 |
* excluded. The University of Notre Dame and its licensors shall not |
22 |
* be liable for any damages suffered by licensee as a result of |
23 |
* using, modifying or distributing the software or its |
24 |
* derivatives. In no event will the University of Notre Dame or its |
25 |
* licensors be liable for any lost revenue, profit or data, or for |
26 |
* direct, indirect, special, consequential, incidental or punitive |
27 |
* damages, however caused and regardless of the theory of liability, |
28 |
* arising out of the use of or inability to use software, even if the |
29 |
* University of Notre Dame has been advised of the possibility of |
30 |
* such damages. |
31 |
* |
32 |
* SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your |
33 |
* research, please cite the appropriate papers when you publish your |
34 |
* work. Good starting points are: |
35 |
* |
36 |
* [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005). |
37 |
* [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006). |
38 |
* [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008). |
39 |
* [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). |
40 |
* [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). |
41 |
*/ |
42 |
|
43 |
#ifdef IS_MPI |
44 |
#include <mpi.h> |
45 |
#endif |
46 |
|
47 |
#include "config.h" |
48 |
#include <cmath> |
49 |
|
50 |
#include "restraints/RestraintForceManager.hpp" |
51 |
#include "restraints/MolecularRestraint.hpp" |
52 |
#include "restraints/ObjectRestraint.hpp" |
53 |
#include "io/RestReader.hpp" |
54 |
#include "utils/simError.h" |
55 |
#include "utils/PhysicalConstants.hpp" |
56 |
#include "utils/StringUtils.hpp" |
57 |
#include "selection/SelectionEvaluator.hpp" |
58 |
#include "selection/SelectionManager.hpp" |
59 |
|
60 |
namespace OpenMD { |
61 |
|
62 |
RestraintForceManager::RestraintForceManager(SimInfo* info): ForceManager(info) { |
63 |
|
64 |
// order of affairs: |
65 |
// |
66 |
// 1) create restraints from the restraintStamps found in the MD |
67 |
// file. |
68 |
// |
69 |
// 2) Create RestraintReader to parse the input files for the ideal |
70 |
// structures. This reader will set reference structures, and will |
71 |
// calculate molecular centers of mass, etc. |
72 |
// |
73 |
// 3) sit around and wait for calcForces to be called. When it comes, |
74 |
// call the normal force manager calcForces, then loop through the |
75 |
// restrained objects and do their restraint forces. |
76 |
|
77 |
currSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot(); |
78 |
Globals* simParam = info_->getSimParams(); |
79 |
|
80 |
if (simParam->haveStatusTime()){ |
81 |
restTime_ = simParam->getStatusTime(); |
82 |
} else { |
83 |
sprintf(painCave.errMsg, |
84 |
"Restraint warning: If you use restraints without setting\n" |
85 |
"\tstatusTime, no restraint data will be written to the rest\n" |
86 |
"\tfile.\n"); |
87 |
painCave.isFatal = 0; |
88 |
simError(); |
89 |
restTime_ = simParam->getRunTime(); |
90 |
} |
91 |
|
92 |
int nRestraintStamps = simParam->getNRestraintStamps(); |
93 |
std::vector<RestraintStamp*> stamp = simParam->getRestraintStamps(); |
94 |
|
95 |
std::vector<int> stuntDoubleIndex; |
96 |
|
97 |
for (int i = 0; i < nRestraintStamps; i++){ |
98 |
|
99 |
std::string myType = toUpperCopy(stamp[i]->getType()); |
100 |
|
101 |
if (myType.compare("MOLECULAR")==0){ |
102 |
|
103 |
int molIndex; |
104 |
Vector3d refCom; |
105 |
|
106 |
if (!stamp[i]->haveMolIndex()) { |
107 |
sprintf(painCave.errMsg, |
108 |
"Restraint Error: A molecular restraint was specified\n" |
109 |
"\twithout providing a value for molIndex.\n"); |
110 |
painCave.isFatal = 1; |
111 |
simError(); |
112 |
} else { |
113 |
molIndex = stamp[i]->getMolIndex(); |
114 |
} |
115 |
|
116 |
if (molIndex < 0) { |
117 |
sprintf(painCave.errMsg, |
118 |
"Restraint Error: A molecular restraint was specified\n" |
119 |
"\twith a molIndex that was less than 0\n"); |
120 |
painCave.isFatal = 1; |
121 |
simError(); |
122 |
} |
123 |
if (molIndex >= info_->getNGlobalMolecules()) { |
124 |
sprintf(painCave.errMsg, |
125 |
"Restraint Error: A molecular restraint was specified with\n" |
126 |
"\ta molIndex that was greater than the total number of molecules\n"); |
127 |
painCave.isFatal = 1; |
128 |
simError(); |
129 |
} |
130 |
|
131 |
Molecule* mol = info_->getMoleculeByGlobalIndex(molIndex); |
132 |
|
133 |
if (mol == NULL) { |
134 |
#ifdef IS_MPI |
135 |
// getMoleculeByGlobalIndex returns a NULL in parallel if |
136 |
// this proc doesn't have the molecule. Do a quick check to |
137 |
// make sure another processor is supposed to have it. |
138 |
|
139 |
int myrank; |
140 |
MPI_Comm_rank( MPI_COMM_WORLD, &myrank); |
141 |
|
142 |
if (info_->getMolToProc(molIndex) == myrank) { |
143 |
|
144 |
// If we were supposed to have it but got a null, then freak out. |
145 |
#endif |
146 |
|
147 |
sprintf(painCave.errMsg, |
148 |
"Restraint Error: A molecular restraint was specified, but\n" |
149 |
"\tno molecule was found with global index %d.\n", |
150 |
molIndex); |
151 |
painCave.isFatal = 1; |
152 |
simError(); |
153 |
|
154 |
#ifdef IS_MPI |
155 |
} |
156 |
#endif |
157 |
} |
158 |
|
159 |
|
160 |
#ifdef IS_MPI |
161 |
// only handle this molecular restraint if this processor owns the |
162 |
// molecule |
163 |
int myrank; |
164 |
MPI_Comm_rank( MPI_COMM_WORLD, &myrank); |
165 |
if (info_->getMolToProc(molIndex) == myrank) { |
166 |
|
167 |
#endif |
168 |
|
169 |
MolecularRestraint* rest = new MolecularRestraint(); |
170 |
|
171 |
std::string restPre("mol_"); |
172 |
std::stringstream restName; |
173 |
restName << restPre << molIndex; |
174 |
rest->setRestraintName(restName.str()); |
175 |
|
176 |
if (stamp[i]->haveDisplacementSpringConstant()) { |
177 |
rest->setDisplacementForceConstant(stamp[i]->getDisplacementSpringConstant()); |
178 |
} |
179 |
if (stamp[i]->haveTwistSpringConstant()) { |
180 |
rest->setTwistForceConstant(stamp[i]->getTwistSpringConstant()); |
181 |
} |
182 |
if (stamp[i]->haveSwingXSpringConstant()) { |
183 |
rest->setSwingXForceConstant(stamp[i]->getSwingXSpringConstant()); |
184 |
} |
185 |
if (stamp[i]->haveSwingYSpringConstant()) { |
186 |
rest->setSwingYForceConstant(stamp[i]->getSwingYSpringConstant()); |
187 |
} |
188 |
if (stamp[i]->haveRestrainedTwistAngle()) { |
189 |
rest->setRestrainedTwistAngle(stamp[i]->getRestrainedTwistAngle() * M_PI/180.0); |
190 |
} |
191 |
if (stamp[i]->haveRestrainedSwingYAngle()) { |
192 |
rest->setRestrainedSwingYAngle(stamp[i]->getRestrainedSwingYAngle() * M_PI/180.0); |
193 |
} |
194 |
if (stamp[i]->haveRestrainedSwingXAngle()) { |
195 |
rest->setRestrainedSwingXAngle(stamp[i]->getRestrainedSwingXAngle() * M_PI/180.0); |
196 |
} |
197 |
if (stamp[i]->havePrint()) { |
198 |
rest->setPrintRestraint(stamp[i]->getPrint()); |
199 |
} |
200 |
|
201 |
restraints_.push_back(rest); |
202 |
mol->addProperty(new RestraintData("Restraint", rest)); |
203 |
restrainedMols_.push_back(mol); |
204 |
#ifdef IS_MPI |
205 |
} |
206 |
#endif |
207 |
} else if (myType.compare("OBJECT") == 0) { |
208 |
|
209 |
std::string objectSelection; |
210 |
|
211 |
if (!stamp[i]->haveObjectSelection()) { |
212 |
sprintf(painCave.errMsg, |
213 |
"Restraint Error: An object restraint was specified\n" |
214 |
"\twithout providing a selection script in the\n" |
215 |
"\tobjectSelection variable.\n"); |
216 |
painCave.isFatal = 1; |
217 |
simError(); |
218 |
} else { |
219 |
objectSelection = stamp[i]->getObjectSelection(); |
220 |
} |
221 |
|
222 |
SelectionEvaluator evaluator(info); |
223 |
SelectionManager seleMan(info); |
224 |
|
225 |
evaluator.loadScriptString(objectSelection); |
226 |
seleMan.setSelectionSet(evaluator.evaluate()); |
227 |
int selectionCount = seleMan.getSelectionCount(); |
228 |
|
229 |
sprintf(painCave.errMsg, |
230 |
"Restraint Info: The specified restraint objectSelection,\n" |
231 |
"\t\t%s\n" |
232 |
"\twill result in %d integrable objects being\n" |
233 |
"\trestrained.\n", objectSelection.c_str(), selectionCount); |
234 |
painCave.severity = OPENMD_INFO; |
235 |
painCave.isFatal = 0; |
236 |
simError(); |
237 |
|
238 |
int selei; |
239 |
StuntDouble* sd; |
240 |
|
241 |
for (sd = seleMan.beginSelected(selei); sd != NULL; |
242 |
sd = seleMan.nextSelected(selei)) { |
243 |
stuntDoubleIndex.push_back(sd->getGlobalIntegrableObjectIndex()); |
244 |
|
245 |
ObjectRestraint* rest = new ObjectRestraint(); |
246 |
|
247 |
if (stamp[i]->haveDisplacementSpringConstant()) { |
248 |
rest->setDisplacementForceConstant(stamp[i]->getDisplacementSpringConstant()); |
249 |
} |
250 |
if (stamp[i]->haveTwistSpringConstant()) { |
251 |
rest->setTwistForceConstant(stamp[i]->getTwistSpringConstant()); |
252 |
} |
253 |
if (stamp[i]->haveSwingXSpringConstant()) { |
254 |
rest->setSwingXForceConstant(stamp[i]->getSwingXSpringConstant()); |
255 |
} |
256 |
if (stamp[i]->haveSwingYSpringConstant()) { |
257 |
rest->setSwingYForceConstant(stamp[i]->getSwingYSpringConstant()); |
258 |
} |
259 |
if (stamp[i]->haveRestrainedTwistAngle()) { |
260 |
rest->setRestrainedTwistAngle(stamp[i]->getRestrainedTwistAngle()); |
261 |
} |
262 |
if (stamp[i]->haveRestrainedSwingXAngle()) { |
263 |
rest->setRestrainedSwingXAngle(stamp[i]->getRestrainedSwingXAngle()); |
264 |
} |
265 |
if (stamp[i]->haveRestrainedSwingYAngle()) { |
266 |
rest->setRestrainedSwingYAngle(stamp[i]->getRestrainedSwingYAngle()); |
267 |
} |
268 |
if (stamp[i]->havePrint()) { |
269 |
rest->setPrintRestraint(stamp[i]->getPrint()); |
270 |
} |
271 |
|
272 |
restraints_.push_back(rest); |
273 |
sd->addProperty(new RestraintData("Restraint", rest)); |
274 |
restrainedObjs_.push_back(sd); |
275 |
} |
276 |
|
277 |
} |
278 |
} |
279 |
|
280 |
// ThermodynamicIntegration subclasses RestraintForceManager, and there |
281 |
// are times when it won't use restraints at all, so only open the |
282 |
// restraint file if we are actually using restraints: |
283 |
|
284 |
if (simParam->getUseRestraints()) { |
285 |
std::string refFile = simParam->getRestraint_file(); |
286 |
RestReader* rr = new RestReader(info, refFile, stuntDoubleIndex); |
287 |
rr->readReferenceStructure(); |
288 |
} |
289 |
|
290 |
restOutput_ = getPrefix(info_->getFinalConfigFileName()) + ".rest"; |
291 |
restOut = new RestWriter(info_, restOutput_.c_str(), restraints_); |
292 |
if(!restOut){ |
293 |
sprintf(painCave.errMsg, "Restraint error: Failed to create RestWriter\n"); |
294 |
painCave.isFatal = 1; |
295 |
simError(); |
296 |
} |
297 |
|
298 |
// todo: figure out the scale factor. Right now, just scale them all to 1 |
299 |
std::vector<Restraint*>::const_iterator resti; |
300 |
for(resti=restraints_.begin(); resti != restraints_.end(); ++resti){ |
301 |
(*resti)->setScaleFactor(1.0); |
302 |
} |
303 |
} |
304 |
|
305 |
RestraintForceManager::~RestraintForceManager(){ |
306 |
if (restOut) |
307 |
delete restOut; |
308 |
} |
309 |
|
310 |
void RestraintForceManager::init() { |
311 |
currRestTime_ = currSnapshot_->getTime(); |
312 |
} |
313 |
|
314 |
void RestraintForceManager::calcForces(){ |
315 |
|
316 |
ForceManager::calcForces(); |
317 |
RealType restPot_local, restPot; |
318 |
|
319 |
restPot_local = doRestraints(1.0); |
320 |
|
321 |
#ifdef IS_MPI |
322 |
MPI_Allreduce(&restPot_local, &restPot, 1, |
323 |
MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); |
324 |
#else |
325 |
restPot = restPot_local; |
326 |
#endif |
327 |
currSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot(); |
328 |
RealType pot = currSnapshot_->getLongRangePotential(); |
329 |
pot += restPot; |
330 |
currSnapshot_->setLongRangePotential(pot); |
331 |
currSnapshot_->setRestraintPotential(restPot); |
332 |
|
333 |
//write out forces and current positions of restrained molecules |
334 |
if (currSnapshot_->getTime() >= currRestTime_){ |
335 |
restOut->writeRest(restInfo_); |
336 |
currRestTime_ += restTime_; |
337 |
} |
338 |
} |
339 |
|
340 |
RealType RestraintForceManager::doRestraints(RealType scalingFactor){ |
341 |
std::vector<Molecule*>::const_iterator rm; |
342 |
GenericData* data; |
343 |
Molecule::IntegrableObjectIterator ioi; |
344 |
MolecularRestraint* mRest; |
345 |
StuntDouble* sd; |
346 |
|
347 |
std::vector<StuntDouble*>::const_iterator ro; |
348 |
ObjectRestraint* oRest; |
349 |
|
350 |
std::map<int, Restraint::RealPair> restInfo; |
351 |
|
352 |
unscaledPotential_ = 0.0; |
353 |
|
354 |
restInfo_.clear(); |
355 |
|
356 |
for(rm=restrainedMols_.begin(); rm != restrainedMols_.end(); ++rm){ |
357 |
|
358 |
// make sure this molecule (*rm) has a generic data for restraints: |
359 |
data = (*rm)->getPropertyByName("Restraint"); |
360 |
if (data != NULL) { |
361 |
// make sure we can reinterpret the generic data as restraint data: |
362 |
RestraintData* restData= dynamic_cast<RestraintData*>(data); |
363 |
if (restData != NULL) { |
364 |
// make sure we can reinterpet the restraint data as a pointer to |
365 |
// an MolecularRestraint: |
366 |
mRest = dynamic_cast<MolecularRestraint*>(restData->getData()); |
367 |
if (mRest == NULL) { |
368 |
sprintf( painCave.errMsg, |
369 |
"Can not cast RestraintData to MolecularRestraint\n"); |
370 |
painCave.severity = OPENMD_ERROR; |
371 |
painCave.isFatal = 1; |
372 |
simError(); |
373 |
} |
374 |
} else { |
375 |
sprintf( painCave.errMsg, |
376 |
"Can not cast GenericData to RestraintData\n"); |
377 |
painCave.severity = OPENMD_ERROR; |
378 |
painCave.isFatal = 1; |
379 |
simError(); |
380 |
} |
381 |
} else { |
382 |
sprintf( painCave.errMsg, "Can not find Restraint for RestrainedObject\n"); |
383 |
painCave.severity = OPENMD_ERROR; |
384 |
painCave.isFatal = 1; |
385 |
simError(); |
386 |
} |
387 |
|
388 |
// phew. At this point, we should have the pointer to the |
389 |
// correct MolecularRestraint in the variable mRest. |
390 |
|
391 |
Vector3d molCom = (*rm)->getCom(); |
392 |
|
393 |
std::vector<Vector3d> struc; |
394 |
std::vector<Vector3d> forces; |
395 |
|
396 |
|
397 |
for(sd = (*rm)->beginIntegrableObject(ioi); sd != NULL; |
398 |
sd = (*rm)->nextIntegrableObject(ioi)) { |
399 |
struc.push_back(sd->getPos()); |
400 |
} |
401 |
|
402 |
mRest->setScaleFactor(scalingFactor); |
403 |
mRest->calcForce(struc, molCom); |
404 |
forces = mRest->getRestraintForces(); |
405 |
int index = 0; |
406 |
|
407 |
for(sd = (*rm)->beginIntegrableObject(ioi); sd != NULL; |
408 |
sd = (*rm)->nextIntegrableObject(ioi)) { |
409 |
sd->addFrc(forces[index]); |
410 |
struc.push_back(sd->getPos()); |
411 |
index++; |
412 |
} |
413 |
|
414 |
unscaledPotential_ += mRest->getUnscaledPotential(); |
415 |
|
416 |
restInfo = mRest->getRestraintInfo(); |
417 |
|
418 |
// only collect data on restraints that we're going to print: |
419 |
if (mRest->getPrintRestraint()) |
420 |
restInfo_.push_back(restInfo); |
421 |
} |
422 |
|
423 |
for(ro=restrainedObjs_.begin(); ro != restrainedObjs_.end(); ++ro){ |
424 |
// make sure this object (*ro) has a generic data for restraints: |
425 |
data = (*ro)->getPropertyByName("Restraint"); |
426 |
if (data != NULL) { |
427 |
// make sure we can reinterpret the generic data as restraint data: |
428 |
RestraintData* restData= dynamic_cast<RestraintData*>(data); |
429 |
if (restData != NULL) { |
430 |
// make sure we can reinterpet the restraint data as a pointer to |
431 |
// an ObjectRestraint: |
432 |
oRest = dynamic_cast<ObjectRestraint*>(restData->getData()); |
433 |
if (oRest == NULL) { |
434 |
sprintf( painCave.errMsg, |
435 |
"Can not cast RestraintData to ObjectRestraint\n"); |
436 |
painCave.severity = OPENMD_ERROR; |
437 |
painCave.isFatal = 1; |
438 |
simError(); |
439 |
} |
440 |
} else { |
441 |
sprintf( painCave.errMsg, |
442 |
"Can not cast GenericData to RestraintData\n"); |
443 |
painCave.severity = OPENMD_ERROR; |
444 |
painCave.isFatal = 1; |
445 |
simError(); |
446 |
} |
447 |
} else { |
448 |
sprintf( painCave.errMsg, "Can not find Restraint for RestrainedObject\n"); |
449 |
painCave.severity = OPENMD_ERROR; |
450 |
painCave.isFatal = 1; |
451 |
simError(); |
452 |
} |
453 |
|
454 |
// phew. At this point, we should have the pointer to the |
455 |
// correct Object restraint in the variable oRest. |
456 |
oRest->setScaleFactor(scalingFactor); |
457 |
|
458 |
Vector3d pos = (*ro)->getPos(); |
459 |
|
460 |
if ( (*ro)->isDirectional() ) { |
461 |
|
462 |
// directional objects may have orientational restraints as well |
463 |
// as positional, so get the rotation matrix first: |
464 |
|
465 |
RotMat3x3d A = (*ro)->getA(); |
466 |
oRest->calcForce(pos, A); |
467 |
(*ro)->addFrc(oRest->getRestraintForce()); |
468 |
(*ro)->addTrq(oRest->getRestraintTorque()); |
469 |
} else { |
470 |
|
471 |
// plain vanilla positional restraints: |
472 |
|
473 |
oRest->calcForce(pos); |
474 |
(*ro)->addFrc(oRest->getRestraintForce()); |
475 |
} |
476 |
|
477 |
unscaledPotential_ += oRest->getUnscaledPotential(); |
478 |
|
479 |
restInfo = oRest->getRestraintInfo(); |
480 |
|
481 |
// only collect data on restraints that we're going to print: |
482 |
if (oRest->getPrintRestraint()) |
483 |
restInfo_.push_back(restInfo); |
484 |
} |
485 |
|
486 |
return unscaledPotential_ * scalingFactor; |
487 |
} |
488 |
} |