OpenMD 3.0
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
Shake.cpp
1/*
2 * Copyright (c) 2004-present, The University of Notre Dame. All rights
3 * reserved.
4 *
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions are met:
7 *
8 * 1. Redistributions of source code must retain the above copyright notice,
9 * this list of conditions and the following disclaimer.
10 *
11 * 2. Redistributions in binary form must reproduce the above copyright notice,
12 * this list of conditions and the following disclaimer in the documentation
13 * and/or other materials provided with the distribution.
14 *
15 * 3. Neither the name of the copyright holder nor the names of its
16 * contributors may be used to endorse or promote products derived from
17 * this software without specific prior written permission.
18 *
19 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
20 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
21 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
22 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
23 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
24 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
25 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
26 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
27 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
28 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
29 * POSSIBILITY OF SUCH DAMAGE.
30 *
31 * SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your
32 * research, please cite the appropriate papers when you publish your
33 * work. Good starting points are:
34 *
35 * [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
36 * [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
37 * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).
38 * [4] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
39 * [5] Kuang & Gezelter, Mol. Phys., 110, 691-701 (2012).
40 * [6] Lamichhane, Gezelter & Newman, J. Chem. Phys. 141, 134109 (2014).
41 * [7] Lamichhane, Newman & Gezelter, J. Chem. Phys. 141, 134110 (2014).
42 * [8] Bhattarai, Newman & Gezelter, Phys. Rev. B 99, 094106 (2019).
43 */
44
45#include "constraints/Shake.hpp"
46
47#ifdef IS_MPI
48#include <mpi.h>
49#endif
50
52#include "utils/simError.h"
53
54namespace OpenMD {
55
56 Shake::Shake(SimInfo* info) :
57 info_(info), maxConsIteration_(10), consTolerance_(1.0e-6),
58 doShake_(false), currConstraintTime_(0.0) {
59 if (info_->getNGlobalConstraints() > 0) doShake_ = true;
60
61 if (!doShake_) return;
62
63 Globals* simParams = info_->getSimParams();
64
65 currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
66 if (simParams->haveConstraintTime()) {
67 constraintTime_ = simParams->getConstraintTime();
68 } else {
69 constraintTime_ = simParams->getStatusTime();
70 }
71
72 constraintOutputFile_ =
73 getPrefix(info_->getFinalConfigFileName()) + ".constraintForces";
74
75 // create ConstraintWriter
76 constraintWriter_ =
77 new ConstraintWriter(info_, constraintOutputFile_.c_str());
78
79 if (!constraintWriter_) {
80 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
81 "Failed to create ConstraintWriter\n");
82 painCave.isFatal = 1;
83 simError();
84 }
85 }
86
87 void Shake::constraintR() {
88 if (!doShake_) return;
89 doConstraint(&Shake::constraintPairR);
90 }
91 void Shake::constraintF() {
92 if (!doShake_) return;
93 doConstraint(&Shake::constraintPairF);
94
95 if (currentSnapshot_->getTime() >= currConstraintTime_) {
96 Molecule* mol;
97 SimInfo::MoleculeIterator mi;
98 ConstraintPair* consPair;
99 Molecule::ConstraintPairIterator cpi;
100 std::list<ConstraintPair*> constraints;
101 for (mol = info_->beginMolecule(mi); mol != NULL;
102 mol = info_->nextMolecule(mi)) {
103 for (consPair = mol->beginConstraintPair(cpi); consPair != NULL;
104 consPair = mol->nextConstraintPair(cpi)) {
105 constraints.push_back(consPair);
106 }
107 }
108
109 constraintWriter_->writeConstraintForces(constraints);
110 currConstraintTime_ += constraintTime_;
111 }
112 }
113
114 void Shake::doConstraint(ConstraintPairFuncPtr func) {
115 if (!doShake_) return;
116
117 Molecule* mol;
118 SimInfo::MoleculeIterator mi;
119 ConstraintElem* consElem;
120 Molecule::ConstraintElemIterator cei;
121 ConstraintPair* consPair;
122 Molecule::ConstraintPairIterator cpi;
123
124 for (mol = info_->beginMolecule(mi); mol != NULL;
125 mol = info_->nextMolecule(mi)) {
126 for (consElem = mol->beginConstraintElem(cei); consElem != NULL;
127 consElem = mol->nextConstraintElem(cei)) {
128 consElem->setMoved(true);
129 consElem->setMoving(false);
130 }
131 }
132
133 // main loop of constraint algorithm
134 int done = 0;
135 int iteration = 0;
136 while (!done && iteration < maxConsIteration_) {
137 done = 1;
138
139 // loop over every constraint pair
140
141 for (mol = info_->beginMolecule(mi); mol != NULL;
142 mol = info_->nextMolecule(mi)) {
143 for (consPair = mol->beginConstraintPair(cpi); consPair != NULL;
144 consPair = mol->nextConstraintPair(cpi)) {
145 // dispatch constraint algorithm
146 if (consPair->isMoved()) {
147 int exeStatus = (this->*func)(consPair);
148
149 switch (exeStatus) {
150 case consFail:
151 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
152 "Constraint failure in Shake::constrainA, "
153 "Constraint Fail\n");
154 painCave.isFatal = 1;
155 simError();
156
157 break;
158 case consSuccess:
159 // constrain the pair by moving two elements
160 done = 0;
161 consPair->getConsElem1()->setMoving(true);
162 consPair->getConsElem2()->setMoving(true);
163 break;
164 case consAlready:
165 // current pair is already constrained, do not need to
166 // move the elements
167 break;
168 default:
169 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
170 "ConstraintAlgorithm::doConstraint() "
171 "Error: unrecognized status");
172 painCave.isFatal = 1;
173 simError();
174 break;
175 }
176 }
177 }
178 } // end for(iter->first())
179
180#ifdef IS_MPI
181 MPI_Allreduce(MPI_IN_PLACE, &done, 1, MPI_INT, MPI_LAND, MPI_COMM_WORLD);
182#endif
183
184 errorCheckPoint();
185
186 for (mol = info_->beginMolecule(mi); mol != NULL;
187 mol = info_->nextMolecule(mi)) {
188 for (consElem = mol->beginConstraintElem(cei); consElem != NULL;
189 consElem = mol->nextConstraintElem(cei)) {
190 consElem->setMoved(consElem->getMoving());
191 consElem->setMoving(false);
192 }
193 }
194
195 iteration++;
196 } // end while
197
198 if (!done) {
199 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
200 "Constraint failure in Shake::constrainA, "
201 "too many iterations: %d\n",
202 iteration);
203 painCave.isFatal = 1;
204 simError();
205 }
206
207 errorCheckPoint();
208 }
209
210 /**
211 * remove constraint force along the bond direction
212 */
213 int Shake::constraintPairR(ConstraintPair* consPair) {
214 ConstraintElem* consElem1 = consPair->getConsElem1();
215 ConstraintElem* consElem2 = consPair->getConsElem2();
216
217 Vector3d posA = consElem1->getPos();
218 Vector3d posB = consElem2->getPos();
219
220 Vector3d pab = posA - posB;
221
222 // periodic boundary condition
223
224 currentSnapshot_->wrapVector(pab);
225
226 RealType pabsq = pab.lengthSquare();
227
228 RealType rabsq = consPair->getConsDistSquare();
229 RealType diffsq = rabsq - pabsq;
230
231 // the original rattle code from alan tidesley
232 if (fabs(diffsq) > (consTolerance_ * rabsq * 2)) {
233 Vector3d oldPosA = consElem1->getPrevPos();
234 Vector3d oldPosB = consElem2->getPrevPos();
235
236 Vector3d rab = oldPosA - oldPosB;
237
238 currentSnapshot_->wrapVector(rab);
239
240 RealType rpab = dot(rab, pab);
241 RealType rpabsq = rpab * rpab;
242
243 if (rpabsq < (rabsq * -diffsq)) { return consFail; }
244
245 RealType rma = 1.0 / consElem1->getMass();
246 RealType rmb = 1.0 / consElem2->getMass();
247
248 RealType gab = diffsq / (2.0 * (rma + rmb) * rpab);
249
250 Vector3d delta = rab * gab;
251
252 // set atom1's position
253 posA += rma * delta;
254 consElem1->setPos(posA);
255
256 // set atom2's position
257 posB -= rmb * delta;
258 consElem2->setPos(posB);
259
260 // report the constraint force back to the constraint pair:
261 consPair->setConstraintForce(gab);
262 return consSuccess;
263 } else
264 return consAlready;
265 }
266
267 /**
268 * remove constraint force along the bond direction
269 */
270 int Shake::constraintPairF(ConstraintPair* consPair) {
271 ConstraintElem* consElem1 = consPair->getConsElem1();
272 ConstraintElem* consElem2 = consPair->getConsElem2();
273
274 Vector3d posA = consElem1->getPos();
275 Vector3d posB = consElem2->getPos();
276
277 Vector3d rab = posA - posB;
278
279 currentSnapshot_->wrapVector(rab);
280
281 Vector3d frcA = consElem1->getFrc();
282 Vector3d frcB = consElem2->getFrc();
283
284 RealType rma = 1.0 / consElem1->getMass();
285 RealType rmb = 1.0 / consElem2->getMass();
286
287 Vector3d fpab = frcA * rma - frcB * rmb;
288
289 RealType gab = fpab.lengthSquare();
290 if (gab < 1.0) gab = 1.0;
291
292 RealType rabsq = rab.lengthSquare();
293 RealType rfab = dot(rab, fpab);
294
295 if (fabs(rfab) > sqrt(rabsq * gab) * consTolerance_) {
296 gab = -rfab / (rabsq * (rma + rmb));
297
298 frcA += rab * gab;
299 frcB -= rab * gab;
300
301 consElem1->setFrc(frcA);
302 consElem2->setFrc(frcB);
303
304 // report the constraint force back to the constraint pair:
305 consPair->setConstraintForce(gab);
306 return consSuccess;
307 } else
308 return consAlready;
309 }
310} // namespace OpenMD
Real lengthSquare()
Returns the squared length of this vector.
Definition Vector.hpp:399
This basic Periodic Table class was originally taken from the data.cpp file in OpenBabel.
Real dot(const DynamicVector< Real > &v1, const DynamicVector< Real > &v2)
Returns the dot product of two DynamicVectors.
std::string getPrefix(const std::string &str)