ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/restraints/RestraintForceManager.cpp
Revision: 1973
Committed: Thu Mar 6 19:34:22 2014 UTC (11 years, 1 month ago) by gezelter
File size: 17787 byte(s)
Log Message:
Some updates for restraints and fluctuating charges in Ewald

File Contents

# Content
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 }

Properties

Name Value
svn:executable *
svn:keywords Author Id Revision Date