OpenMD 3.0
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
Integrator.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 "integrators/Integrator.hpp"
46
47#include <memory>
48#include <utility>
49
50#include "brains/Snapshot.hpp"
51#include "flucq/FluctuatingChargeDamped.hpp"
52#include "flucq/FluctuatingChargeLangevin.hpp"
53#include "flucq/FluctuatingChargeNVE.hpp"
54#include "flucq/FluctuatingChargeNVT.hpp"
55#include "integrators/DLM.hpp"
56#include "rnemd/MethodFactory.hpp"
57#include "rnemd/RNEMD.hpp"
58#include "rnemd/SPFForceManager.hpp"
59#include "utils/CI_String.hpp"
60#include "utils/CaseConversion.hpp"
61#include "utils/simError.h"
62
63namespace OpenMD {
64
65 Integrator::Integrator(SimInfo* info) :
66 info_(info), forceMan_(NULL), rotAlgo_(NULL), flucQ_(NULL), rattle_(NULL),
67 velocitizer_(nullptr), needPotential(false), needVirial(false),
68 needReset(false), needVelocityScaling(false), useRNEMD(false),
69 dumpWriter(NULL), statWriter(NULL), thermo(info_),
70 snap(info_->getSnapshotManager()->getCurrentSnapshot()) {
71 simParams = info->getSimParams();
72
73 if (simParams->haveDt()) {
74 dt = simParams->getDt();
75 dt2 = 0.5 * dt;
76 } else {
77 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
78 "Integrator Error: dt is not set\n");
79 painCave.isFatal = 1;
80 simError();
81 }
82
83 if (simParams->haveRunTime()) {
84 runTime = simParams->getRunTime();
85 } else {
86 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
87 "Integrator Error: runTime is not set\n");
88 painCave.isFatal = 1;
89 simError();
90 }
91
92 // set the status, sample, and thermal kick times
93 if (simParams->haveSampleTime()) {
94 sampleTime = simParams->getSampleTime();
95 statusTime = sampleTime;
96 } else {
97 sampleTime = simParams->getRunTime();
98 statusTime = sampleTime;
99 }
100
101 if (simParams->haveStatusTime()) {
102 statusTime = simParams->getStatusTime();
103 }
104
105 if (simParams->haveThermalTime()) {
106 thermalTime = simParams->getThermalTime();
107 } else {
108 thermalTime = simParams->getRunTime();
109 }
110
111 if (!simParams->getUseInitalTime()) { snap->setTime(0.0); }
112
113 if (simParams->haveResetTime()) {
114 needReset = true;
115 resetTime = simParams->getResetTime();
116 }
117
118 // check for the temperature set flag (velocity scaling)
119 needVelocityScaling = false;
120 if (simParams->haveTempSet()) {
121 needVelocityScaling = simParams->getTempSet();
122 }
123
124 if (needVelocityScaling) {
125 if (simParams->haveTargetTemp()) {
126 targetScalingTemp = simParams->getTargetTemp();
127 } else {
128 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
129 "Integrator Error: Target Temperature must be set to turn on "
130 "tempSet\n");
131 painCave.isFatal = 1;
132 simError();
133 }
134 }
135
136 // Create a default a velocitizer: If the subclass wants to use
137 // a different velocitizer, use setVelocitizer
138 velocitizer_ = std::make_unique<Velocitizer>(info);
139
140 if (simParams->getRNEMDParameters()->haveUseRNEMD()) {
141 useRNEMD = simParams->getRNEMDParameters()->getUseRNEMD();
142
143 if (useRNEMD) {
144 // ForceManager will only be changed if SPF-RNEMD is enabled
145 if (Utils::traits_cast<Utils::ci_char_traits>(
146 simParams->getRNEMDParameters()->getMethod()) == "SPF") {
147 if (toUpperCopy(simParams->getEnsemble()) != "SPF") {
148 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
149 "The SPF RNEMD method requires the SPF Ensemble.\n");
150 painCave.isFatal = 1;
151 simError();
152 }
153
154 forceMan_ = new RNEMD::SPFForceManager(info);
155 }
156
157 RNEMD::MethodFactory rnemdMethod {
158 simParams->getRNEMDParameters()->getMethod()};
159 rnemd_ = rnemdMethod.create(info, forceMan_);
160
161 if (simParams->getRNEMDParameters()->haveExchangeTime()) {
162 RNEMD_exchangeTime =
163 simParams->getRNEMDParameters()->getExchangeTime();
164
165 // check to make sure exchange time is a multiple of dt;
166 RealType newET = ceil(RNEMD_exchangeTime / dt) * dt;
167 if (fabs(newET - RNEMD_exchangeTime) > 1e-6) {
168 RNEMD_exchangeTime = newET;
169 }
170 }
171 }
172 }
173
174 if (forceMan_ == NULL) forceMan_ = new ForceManager(info);
175
176 rotAlgo_ = new DLM();
177 rattle_ = new Rattle(info);
178
179 if (simParams->getFluctuatingChargeParameters()->havePropagator()) {
180 std::string prop = toUpperCopy(
181 simParams->getFluctuatingChargeParameters()->getPropagator());
182 if (prop.compare("NVT") == 0) {
183 flucQ_ = new FluctuatingChargeNVT(info);
184 } else if (prop.compare("NVE") == 0) {
185 flucQ_ = new FluctuatingChargeNVE(info);
186 } else if (prop.compare("LANGEVIN") == 0) {
187 flucQ_ = new FluctuatingChargeLangevin(info);
188 } else if (prop.compare("DAMPED") == 0) {
189 flucQ_ = new FluctuatingChargeDamped(info);
190 } else {
191 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
192 "Integrator Error: Unknown Fluctuating Charge propagator (%s) "
193 "requested\n",
194 simParams->getFluctuatingChargeParameters()
195 ->getPropagator()
196 .c_str());
197 painCave.isFatal = 1;
198 }
199 }
200
201 flucQ_->setForceManager(forceMan_);
202 }
203
205 delete forceMan_;
206 delete flucQ_;
207 delete rotAlgo_;
208 delete rattle_;
209 delete dumpWriter;
210 delete stats;
211 delete statWriter;
212 }
213
214 void Integrator::updateSizes() {
215 doUpdateSizes();
216 flucQ_->updateSizes();
217 }
218
219 void Integrator::setVelocitizer(std::unique_ptr<Velocitizer> velocitizer) {
220 velocitizer_ = std::move(velocitizer);
221 }
222
223 void Integrator::setFluctuatingChargePropagator(
225 if (prop != flucQ_ && flucQ_ != NULL) { delete flucQ_; }
226 flucQ_ = prop;
227 if (forceMan_ != NULL) { flucQ_->setForceManager(forceMan_); }
228 }
229
230 void Integrator::setRotationAlgorithm(RotationAlgorithm* algo) {
231 if (algo != rotAlgo_ && rotAlgo_ != NULL) { delete rotAlgo_; }
232
233 rotAlgo_ = algo;
234 }
235
236 void Integrator::setRNEMD(std::unique_ptr<RNEMD::RNEMD> rnemd) {
237 rnemd_ = std::move(rnemd);
238 }
239
240 void Integrator::integrate() {
241 initialize();
242
243 while (snap->getTime() <= runTime) {
244 preStep();
245 step();
246 postStep();
247 }
248
249 finalize();
250 }
251
252 void Integrator::saveConservedQuantity() {
253 snap->setConservedQuantity(calcConservedQuantity());
254 }
255
256 void Integrator::initialize() {
257 forceMan_->initialize();
258
259 // remove center of mass drift velocity (in case we passed in a
260 // configuration that was drifting)
261 velocitizer_->removeComDrift();
262
263 // find the initial fluctuating charges.
264 flucQ_->initialize();
265 // initialize the forces before the first step
266 calcForce();
267 // execute the constraint algorithm to make sure that the system is
268 // constrained at the very beginning
269 if (info_->getNGlobalConstraints() > 0) {
270 rattle_->constraintA();
271 calcForce();
272 rattle_->constraintB();
273 // copy the current snapshot to previous snapshot
274 info_->getSnapshotManager()->advance();
275 }
276
277 if (needVelocityScaling) { velocitizer_->randomize(targetScalingTemp); }
278
279 dumpWriter = createDumpWriter();
280 statWriter = createStatWriter();
281 dumpWriter->writeDumpAndEor();
282
283 progressBar = std::make_unique<ProgressBar>();
284
285 // save statistics, before writeStat, we must save statistics
286 saveConservedQuantity();
287 stats->collectStats();
288
289 if (useRNEMD) rnemd_->getStarted();
290
291 statWriter->writeStat();
292
293 currSample = sampleTime + snap->getTime();
294 currStatus = statusTime + snap->getTime();
295 currThermal = thermalTime + snap->getTime();
296 if (needReset) { currReset = resetTime + snap->getTime(); }
297 if (useRNEMD) { currRNEMD = RNEMD_exchangeTime + snap->getTime(); }
298 needPotential = false;
299 needVirial = false;
300 }
301
302 void Integrator::preStep() {
303 RealType difference = snap->getTime() + dt - currStatus;
304
305 if (difference > 0 || fabs(difference) <= OpenMD::epsilon) {
306 needPotential = true;
307 needVirial = true;
308 }
309 }
310
311 void Integrator::calcForce() {
312 forceMan_->calcForces();
313 flucQ_->applyConstraints();
314 }
315
316 void Integrator::postStep() {
317 RealType difference;
318
319 saveConservedQuantity();
320
321 if (needVelocityScaling) {
322 difference = snap->getTime() - currThermal;
323
324 if (difference > 0 || fabs(difference) <= OpenMD::epsilon) {
325 velocitizer_->randomize(targetScalingTemp);
326 currThermal += thermalTime;
327 }
328 }
329
330 if (useRNEMD) {
331 difference = snap->getTime() - currRNEMD;
332
333 if (difference > 0 || fabs(difference) <= OpenMD::epsilon) {
334 rnemd_->doRNEMD();
335 currRNEMD += RNEMD_exchangeTime;
336 }
337 rnemd_->collectData();
338 }
339
340 difference = snap->getTime() - currSample;
341
342 if (difference > 0 || fabs(difference) <= OpenMD::epsilon) {
343 dumpWriter->writeDumpAndEor();
344 currSample += sampleTime;
345 }
346
347 difference = snap->getTime() - currStatus;
348
349 if (difference > 0 || fabs(difference) <= OpenMD::epsilon) {
350 stats->collectStats();
351
352 if (useRNEMD) { rnemd_->writeOutputFile(); }
353
354 statWriter->writeStat();
355
356 progressBar->setStatus(snap->getTime(), runTime);
357 progressBar->update();
358
359 needPotential = false;
360 needVirial = false;
361 currStatus += statusTime;
362 }
363
364 difference = snap->getTime() - currReset;
365
366 if (needReset && (difference > 0 || fabs(difference) <= OpenMD::epsilon)) {
367 resetIntegrator();
368 currReset += resetTime;
369 }
370 // save snapshot
371 info_->getSnapshotManager()->advance();
372
373 // increase time
374 snap->increaseTime(dt);
375 }
376
377 void Integrator::finalize() {
378 dumpWriter->writeEor();
379 if (useRNEMD) { rnemd_->writeOutputFile(); }
380 progressBar->setStatus(runTime, runTime);
381 progressBar->update();
382
383 statWriter->writeStatReport();
384 }
385
386 DumpWriter* Integrator::createDumpWriter() { return new DumpWriter(info_); }
387
388 StatWriter* Integrator::createStatWriter() {
389 stats = new Stats(info_);
390 statWriter = new StatWriter(info_->getStatFileName(), stats);
391 statWriter->setReportFileName(info_->getReportFileName());
392
393 return statWriter;
394 }
395} // namespace OpenMD
abstract class for propagating fluctuating charge variables
ForceManager is responsible for calculating both the short range (bonded) interactions and long range...
virtual ~Integrator()
Default Destructor.
virtual void step()=0
Computes an integration step from t to t+dt.
Velocity Verlet Constraint Algorithm.
Definition Rattle.hpp:58
abstract class for rotation
One of the heavy-weight classes of OpenMD, SimInfo maintains objects and variables relating to the cu...
Definition SimInfo.hpp:93
SnapshotManager * getSnapshotManager()
Returns the snapshot manager.
Definition SimInfo.hpp:248
A configurable Statistics Writer.
This basic Periodic Table class was originally taken from the data.cpp file in OpenBabel.