ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/restraints/RestraintForceManager.cpp
Revision: 2024
Committed: Thu Oct 16 19:13:51 2014 UTC (10 years, 6 months ago) by gezelter
File size: 17771 byte(s)
Log Message:
Added Radial and Z-projected velocity autocorrelation functions
Started to add SequentialProps program
Mucking about with angular restraint potentials

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 currSnapshot_->setRestraintPotential(restPot);
333
334 //write out forces and current positions of restrained molecules
335 if (currSnapshot_->getTime() >= currRestTime_){
336 restOut->writeRest(restInfo_);
337 currRestTime_ += restTime_;
338 }
339 }
340
341 RealType RestraintForceManager::doRestraints(RealType scalingFactor){
342 std::vector<Molecule*>::const_iterator rm;
343 GenericData* data;
344 Molecule::IntegrableObjectIterator ioi;
345 MolecularRestraint* mRest;
346 StuntDouble* sd;
347
348 std::vector<StuntDouble*>::const_iterator ro;
349 ObjectRestraint* oRest;
350
351 std::map<int, Restraint::RealPair> restInfo;
352
353 unscaledPotential_ = 0.0;
354
355 restInfo_.clear();
356
357 for(rm=restrainedMols_.begin(); rm != restrainedMols_.end(); ++rm){
358
359 // make sure this molecule (*rm) has a generic data for restraints:
360 data = (*rm)->getPropertyByName("Restraint");
361 if (data != NULL) {
362 // make sure we can reinterpret the generic data as restraint data:
363 RestraintData* restData= dynamic_cast<RestraintData*>(data);
364 if (restData != NULL) {
365 // make sure we can reinterpet the restraint data as a pointer to
366 // an MolecularRestraint:
367 mRest = dynamic_cast<MolecularRestraint*>(restData->getData());
368 if (mRest == NULL) {
369 sprintf( painCave.errMsg,
370 "Can not cast RestraintData to MolecularRestraint\n");
371 painCave.severity = OPENMD_ERROR;
372 painCave.isFatal = 1;
373 simError();
374 }
375 } else {
376 sprintf( painCave.errMsg,
377 "Can not cast GenericData to RestraintData\n");
378 painCave.severity = OPENMD_ERROR;
379 painCave.isFatal = 1;
380 simError();
381 }
382 } else {
383 sprintf( painCave.errMsg, "Can not find Restraint for RestrainedObject\n");
384 painCave.severity = OPENMD_ERROR;
385 painCave.isFatal = 1;
386 simError();
387 }
388
389 // phew. At this point, we should have the pointer to the
390 // correct MolecularRestraint in the variable mRest.
391
392 Vector3d molCom = (*rm)->getCom();
393
394 std::vector<Vector3d> struc;
395 std::vector<Vector3d> forces;
396
397
398 for(sd = (*rm)->beginIntegrableObject(ioi); sd != NULL;
399 sd = (*rm)->nextIntegrableObject(ioi)) {
400 struc.push_back(sd->getPos());
401 }
402
403 mRest->setScaleFactor(scalingFactor);
404 mRest->calcForce(struc, molCom);
405 forces = mRest->getRestraintForces();
406 int index = 0;
407
408 for(sd = (*rm)->beginIntegrableObject(ioi); sd != NULL;
409 sd = (*rm)->nextIntegrableObject(ioi)) {
410 sd->addFrc(forces[index]);
411 struc.push_back(sd->getPos());
412 index++;
413 }
414
415 unscaledPotential_ += mRest->getUnscaledPotential();
416
417 // only collect data on restraints that we're going to print:
418 if (mRest->getPrintRestraint())
419 restInfo = mRest->getRestraintInfo();
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 // only collect data on restraints that we're going to print:
480 if (oRest->getPrintRestraint())
481 restInfo = oRest->getRestraintInfo();
482 restInfo_.push_back(restInfo);
483 }
484
485 return unscaledPotential_ * scalingFactor;
486 }
487 }

Properties

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