ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/restraints/RestraintForceManager.cpp
Revision: 2020
Committed: Mon Sep 22 19:18:35 2014 UTC (10 years, 7 months ago) by gezelter
File size: 17897 byte(s)
Log Message:
Fixes for restraints, renaming of UniformField

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 #ifdef IS_MPI
230 MPI_Allreduce(MPI_IN_PLACE, &selectionCount, 1, MPI_INT, MPI_SUM,
231 MPI_COMM_WORLD);
232 #endif
233
234 sprintf(painCave.errMsg,
235 "Restraint Info: The specified restraint objectSelection,\n"
236 "\t\t%s\n"
237 "\twill result in %d integrable objects being\n"
238 "\trestrained.\n", objectSelection.c_str(), selectionCount);
239 painCave.severity = OPENMD_INFO;
240 painCave.isFatal = 0;
241 simError();
242
243 int selei;
244 StuntDouble* sd;
245
246 for (sd = seleMan.beginSelected(selei); sd != NULL;
247 sd = seleMan.nextSelected(selei)) {
248 stuntDoubleIndex.push_back(sd->getGlobalIntegrableObjectIndex());
249
250 ObjectRestraint* rest = new ObjectRestraint();
251
252 if (stamp[i]->haveDisplacementSpringConstant()) {
253 rest->setDisplacementForceConstant(stamp[i]->getDisplacementSpringConstant());
254 }
255 if (stamp[i]->haveTwistSpringConstant()) {
256 rest->setTwistForceConstant(stamp[i]->getTwistSpringConstant());
257 }
258 if (stamp[i]->haveSwingXSpringConstant()) {
259 rest->setSwingXForceConstant(stamp[i]->getSwingXSpringConstant());
260 }
261 if (stamp[i]->haveSwingYSpringConstant()) {
262 rest->setSwingYForceConstant(stamp[i]->getSwingYSpringConstant());
263 }
264 if (stamp[i]->haveRestrainedTwistAngle()) {
265 rest->setRestrainedTwistAngle(stamp[i]->getRestrainedTwistAngle());
266 }
267 if (stamp[i]->haveRestrainedSwingXAngle()) {
268 rest->setRestrainedSwingXAngle(stamp[i]->getRestrainedSwingXAngle());
269 }
270 if (stamp[i]->haveRestrainedSwingYAngle()) {
271 rest->setRestrainedSwingYAngle(stamp[i]->getRestrainedSwingYAngle());
272 }
273 if (stamp[i]->havePrint()) {
274 rest->setPrintRestraint(stamp[i]->getPrint());
275 }
276
277 restraints_.push_back(rest);
278 sd->addProperty(new RestraintData("Restraint", rest));
279 restrainedObjs_.push_back(sd);
280 }
281
282 }
283 }
284
285 // ThermodynamicIntegration subclasses RestraintForceManager, and there
286 // are times when it won't use restraints at all, so only open the
287 // restraint file if we are actually using restraints:
288
289 if (simParam->getUseRestraints()) {
290 std::string refFile = simParam->getRestraint_file();
291 RestReader* rr = new RestReader(info, refFile, stuntDoubleIndex);
292 rr->readReferenceStructure();
293 }
294
295 restOutput_ = getPrefix(info_->getFinalConfigFileName()) + ".rest";
296 restOut = new RestWriter(info_, restOutput_.c_str(), restraints_);
297 if(!restOut){
298 sprintf(painCave.errMsg, "Restraint error: Failed to create RestWriter\n");
299 painCave.isFatal = 1;
300 simError();
301 }
302
303 // todo: figure out the scale factor. Right now, just scale them all to 1
304 std::vector<Restraint*>::const_iterator resti;
305 for(resti=restraints_.begin(); resti != restraints_.end(); ++resti){
306 (*resti)->setScaleFactor(1.0);
307 }
308 }
309
310 RestraintForceManager::~RestraintForceManager(){
311 if (restOut)
312 delete restOut;
313 }
314
315 void RestraintForceManager::init() {
316 currRestTime_ = currSnapshot_->getTime();
317 }
318
319 void RestraintForceManager::calcForces(){
320
321 ForceManager::calcForces();
322 RealType restPot(0.0);
323
324 restPot = doRestraints(1.0);
325
326 #ifdef IS_MPI
327 MPI_Allreduce(MPI_IN_PLACE, &restPot, 1, MPI_REALTYPE, MPI_SUM,
328 MPI_COMM_WORLD);
329 #endif
330
331 currSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
332 RealType pot = currSnapshot_->getLongRangePotential();
333 pot += restPot;
334 currSnapshot_->setLongRangePotential(pot);
335 currSnapshot_->setRestraintPotential(restPot);
336
337 //write out forces and current positions of restrained molecules
338 if (currSnapshot_->getTime() >= currRestTime_){
339 restOut->writeRest(restInfo_);
340 currRestTime_ += restTime_;
341 }
342 }
343
344 RealType RestraintForceManager::doRestraints(RealType scalingFactor){
345 std::vector<Molecule*>::const_iterator rm;
346 GenericData* data;
347 Molecule::IntegrableObjectIterator ioi;
348 MolecularRestraint* mRest;
349 StuntDouble* sd;
350
351 std::vector<StuntDouble*>::const_iterator ro;
352 ObjectRestraint* oRest;
353
354 std::map<int, Restraint::RealPair> restInfo;
355
356 unscaledPotential_ = 0.0;
357
358 restInfo_.clear();
359
360 for(rm=restrainedMols_.begin(); rm != restrainedMols_.end(); ++rm){
361
362 // make sure this molecule (*rm) has a generic data for restraints:
363 data = (*rm)->getPropertyByName("Restraint");
364 if (data != NULL) {
365 // make sure we can reinterpret the generic data as restraint data:
366 RestraintData* restData= dynamic_cast<RestraintData*>(data);
367 if (restData != NULL) {
368 // make sure we can reinterpet the restraint data as a pointer to
369 // an MolecularRestraint:
370 mRest = dynamic_cast<MolecularRestraint*>(restData->getData());
371 if (mRest == NULL) {
372 sprintf( painCave.errMsg,
373 "Can not cast RestraintData to MolecularRestraint\n");
374 painCave.severity = OPENMD_ERROR;
375 painCave.isFatal = 1;
376 simError();
377 }
378 } else {
379 sprintf( painCave.errMsg,
380 "Can not cast GenericData to RestraintData\n");
381 painCave.severity = OPENMD_ERROR;
382 painCave.isFatal = 1;
383 simError();
384 }
385 } else {
386 sprintf( painCave.errMsg, "Can not find Restraint for RestrainedObject\n");
387 painCave.severity = OPENMD_ERROR;
388 painCave.isFatal = 1;
389 simError();
390 }
391
392 // phew. At this point, we should have the pointer to the
393 // correct MolecularRestraint in the variable mRest.
394
395 Vector3d molCom = (*rm)->getCom();
396
397 std::vector<Vector3d> struc;
398 std::vector<Vector3d> forces;
399
400
401 for(sd = (*rm)->beginIntegrableObject(ioi); sd != NULL;
402 sd = (*rm)->nextIntegrableObject(ioi)) {
403 struc.push_back(sd->getPos());
404 }
405
406 mRest->setScaleFactor(scalingFactor);
407 mRest->calcForce(struc, molCom);
408 forces = mRest->getRestraintForces();
409 int index = 0;
410
411 for(sd = (*rm)->beginIntegrableObject(ioi); sd != NULL;
412 sd = (*rm)->nextIntegrableObject(ioi)) {
413 sd->addFrc(forces[index]);
414 struc.push_back(sd->getPos());
415 index++;
416 }
417
418 unscaledPotential_ += mRest->getUnscaledPotential();
419
420 // only collect data on restraints that we're going to print:
421 if (mRest->getPrintRestraint())
422 restInfo = mRest->getRestraintInfo();
423 restInfo_.push_back(restInfo);
424 }
425
426 for(ro=restrainedObjs_.begin(); ro != restrainedObjs_.end(); ++ro){
427 // make sure this object (*ro) has a generic data for restraints:
428 data = (*ro)->getPropertyByName("Restraint");
429 if (data != NULL) {
430 // make sure we can reinterpret the generic data as restraint data:
431 RestraintData* restData= dynamic_cast<RestraintData*>(data);
432 if (restData != NULL) {
433 // make sure we can reinterpet the restraint data as a pointer to
434 // an ObjectRestraint:
435 oRest = dynamic_cast<ObjectRestraint*>(restData->getData());
436 if (oRest == NULL) {
437 sprintf( painCave.errMsg,
438 "Can not cast RestraintData to ObjectRestraint\n");
439 painCave.severity = OPENMD_ERROR;
440 painCave.isFatal = 1;
441 simError();
442 }
443 } else {
444 sprintf( painCave.errMsg,
445 "Can not cast GenericData to RestraintData\n");
446 painCave.severity = OPENMD_ERROR;
447 painCave.isFatal = 1;
448 simError();
449 }
450 } else {
451 sprintf( painCave.errMsg, "Can not find Restraint for RestrainedObject\n");
452 painCave.severity = OPENMD_ERROR;
453 painCave.isFatal = 1;
454 simError();
455 }
456
457 // phew. At this point, we should have the pointer to the
458 // correct Object restraint in the variable oRest.
459 oRest->setScaleFactor(scalingFactor);
460
461 Vector3d pos = (*ro)->getPos();
462
463 if ( (*ro)->isDirectional() ) {
464
465 // directional objects may have orientational restraints as well
466 // as positional, so get the rotation matrix first:
467
468 RotMat3x3d A = (*ro)->getA();
469 oRest->calcForce(pos, A);
470 (*ro)->addFrc(oRest->getRestraintForce());
471 (*ro)->addTrq(oRest->getRestraintTorque());
472 } else {
473
474 // plain vanilla positional restraints:
475
476 oRest->calcForce(pos);
477 (*ro)->addFrc(oRest->getRestraintForce());
478 }
479
480 unscaledPotential_ += oRest->getUnscaledPotential();
481
482 // only collect data on restraints that we're going to print:
483 if (oRest->getPrintRestraint())
484 restInfo = oRest->getRestraintInfo();
485 restInfo_.push_back(restInfo);
486 }
487
488 return unscaledPotential_ * scalingFactor;
489 }
490 }

Properties

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