OpenMD 3.2
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
ForceManager.cpp
Go to the documentation of this file.
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 following paper when you publish your work:
33 *
34 * [1] Drisko et al., J. Open Source Softw. 9, 7004 (2024).
35 *
36 * Good starting points for code and simulation methodology are:
37 *
38 * [2] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
39 * [3] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
40 * [4] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).
41 * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
42 * [6] Kuang & Gezelter, Mol. Phys., 110, 691-701 (2012).
43 * [7] Lamichhane, Gezelter & Newman, J. Chem. Phys. 141, 134109 (2014).
44 * [8] Bhattarai, Newman & Gezelter, Phys. Rev. B 99, 094106 (2019).
45 * [9] Drisko & Gezelter, J. Chem. Theory Comput. 20, 4986-4997 (2024).
46 */
47
48/**
49 * @file ForceManager.cpp
50 * @author tlin
51 * @date 11/09/2004
52 * @version 1.0
53 */
54
55#include "brains/ForceManager.hpp"
56
57#define __OPENMD_C
58
59#include <cstdio>
60#include <iomanip>
61#include <iostream>
62
63#include "constraints/ZconstraintForceModifier.hpp"
64#include "integrators/LDForceModifier.hpp"
65#include "integrators/LHDForceModifier.hpp"
66#include "integrators/LangevinHullForceModifier.hpp"
67#include "nonbonded/NonBondedInteraction.hpp"
68#include "parallel/ForceMatrixDecomposition.hpp"
69#include "perturbations/Light.hpp"
70#include "perturbations/LightParameters.hpp"
71#include "perturbations/MagneticField.hpp"
72#include "perturbations/UniformField.hpp"
73#include "perturbations/UniformGradient.hpp"
74#include "primitives/Bend.hpp"
75#include "primitives/Bond.hpp"
76#include "primitives/Inversion.hpp"
77#include "primitives/Molecule.hpp"
78#include "primitives/Torsion.hpp"
79#include "restraints/RestraintForceModifier.hpp"
80#include "restraints/ThermoIntegrationForceModifier.hpp"
81#include "utils/MemoryUtils.hpp"
82#include "utils/simError.h"
83
84using namespace std;
85namespace OpenMD {
86
87 ForceManager::ForceManager(SimInfo* info) :
88 initialized_(false), info_(info), switcher_(NULL), seleMan_(info),
89 evaluator_(info) {
90 forceField_ = info_->getForceField();
91 interactionMan_ = new InteractionManager();
92 fDecomp_ = new ForceMatrixDecomposition(info_, interactionMan_);
93 thermo = new Thermo(info_);
94 }
95
96 ForceManager::~ForceManager() {
97 Utils::deletePointers(forceModifiers_);
98
99 delete switcher_;
100 delete interactionMan_;
101 delete fDecomp_;
102 delete thermo;
103 }
104
105 /**
106 * setupCutoffs
107 *
108 * Sets the values of cutoffRadius, switchingRadius, and cutoffMethod
109 *
110 * cutoffRadius : realType
111 * If the cutoffRadius was explicitly set, use that value.
112 * If the cutoffRadius was not explicitly set:
113 * Are there electrostatic atoms? Use 12.0 Angstroms.
114 * No electrostatic atoms? Poll the atom types present in the
115 * simulation for suggested cutoff values (e.g. 2.5 * sigma).
116 * Use the maximum suggested value that was found.
117 *
118 * cutoffMethod : (one of HARD, SWITCHED, SHIFTED_FORCE, TAYLOR_SHIFTED,
119 * SHIFTED_POTENTIAL, or EWALD_FULL)
120 * If cutoffMethod was explicitly set, use that choice.
121 * If cutoffMethod was not explicitly set, use SHIFTED_FORCE
122 *
123 * switchingRadius : realType
124 * If the cutoffMethod was set to SWITCHED:
125 * If the switchingRadius was explicitly set, use that value
126 * (but do a sanity check first).
127 * If the switchingRadius was not explicitly set: use 0.85 *
128 * cutoffRadius_
129 * If the cutoffMethod was not set to SWITCHED:
130 * Set switchingRadius equal to cutoffRadius for safety.
131 */
132 void ForceManager::setupCutoffs() {
133 Globals* simParams_ = info_->getSimParams();
134 int mdFileVersion;
135 rCut_ = 0.0; // Needs a value for a later max() call;
136
137 if (simParams_->haveMDfileVersion())
138 mdFileVersion = simParams_->getMDfileVersion();
139 else
140 mdFileVersion = 0;
141
142 // We need the list of simulated atom types to figure out cutoffs
143 // as well as long range corrections.
144
145 AtomTypeSet::iterator i;
146 AtomTypeSet atomTypes_;
147 atomTypes_ = info_->getSimulatedAtomTypes();
148
149 if (simParams_->haveCutoffRadius()) {
150 rCut_ = simParams_->getCutoffRadius();
151 } else {
152 if (info_->usesElectrostaticAtoms()) {
153 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
154 "ForceManager::setupCutoffs: No value was set for the "
155 "cutoffRadius.\n"
156 "\tOpenMD will use a default value of 12.0 angstroms "
157 "for the cutoffRadius.\n");
158 painCave.isFatal = 0;
159 painCave.severity = OPENMD_INFO;
160 simError();
161 rCut_ = 12.0;
162 } else {
163 RealType thisCut;
164 for (i = atomTypes_.begin(); i != atomTypes_.end(); ++i) {
165 thisCut = interactionMan_->getSuggestedCutoffRadius((*i));
166 rCut_ = max(thisCut, rCut_);
167 }
168 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
169 "ForceManager::setupCutoffs: No value was set for the "
170 "cutoffRadius.\n"
171 "\tOpenMD will use %lf angstroms.\n",
172 rCut_);
173 painCave.isFatal = 0;
174 painCave.severity = OPENMD_INFO;
175 simError();
176 }
177 }
178
179 fDecomp_->setCutoffRadius(rCut_);
180 interactionMan_->setCutoffRadius(rCut_);
181 rCutSq_ = rCut_ * rCut_;
182
183 map<string, CutoffMethod> stringToCutoffMethod;
184 stringToCutoffMethod["HARD"] = HARD;
185 stringToCutoffMethod["SWITCHED"] = SWITCHED;
186 stringToCutoffMethod["SHIFTED_POTENTIAL"] = SHIFTED_POTENTIAL;
187 stringToCutoffMethod["SHIFTED_FORCE"] = SHIFTED_FORCE;
188 stringToCutoffMethod["TAYLOR_SHIFTED"] = TAYLOR_SHIFTED;
189 stringToCutoffMethod["EWALD_FULL"] = EWALD_FULL;
190
191 if (simParams_->haveCutoffMethod()) {
192 string cutMeth = toUpperCopy(simParams_->getCutoffMethod());
193 map<string, CutoffMethod>::iterator i;
194 i = stringToCutoffMethod.find(cutMeth);
195 if (i == stringToCutoffMethod.end()) {
196 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
197 "ForceManager::setupCutoffs: Could not find chosen "
198 "cutoffMethod %s\n"
199 "\tShould be one of: "
200 "HARD, SWITCHED, SHIFTED_POTENTIAL, TAYLOR_SHIFTED,\n"
201 "\tSHIFTED_FORCE, or EWALD_FULL\n",
202 cutMeth.c_str());
203 painCave.isFatal = 1;
204 painCave.severity = OPENMD_ERROR;
205 simError();
206 } else {
207 cutoffMethod_ = i->second;
208 }
209 } else {
210 if (mdFileVersion > 1) {
211 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
212 "ForceManager::setupCutoffs: No value was set "
213 "for the cutoffMethod.\n"
214 "\tOpenMD will use SHIFTED_FORCE.\n");
215 painCave.isFatal = 0;
216 painCave.severity = OPENMD_INFO;
217 simError();
218 cutoffMethod_ = SHIFTED_FORCE;
219 } else {
220 // handle the case where the old file version was in play
221 // (there should be no cutoffMethod, so we have to deduce it
222 // from other data).
223
224 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
225 "ForceManager::setupCutoffs : DEPRECATED FILE FORMAT!\n"
226 "\tOpenMD found a file which does not set a cutoffMethod.\n"
227 "\tOpenMD will attempt to deduce a cutoffMethod using the\n"
228 "\tbehavior of the older (version 1) code. To remove this\n"
229 "\twarning, add an explicit cutoffMethod and change the top\n"
230 "\tof the file so that it begins with <OpenMD version=2>\n");
231 painCave.isFatal = 0;
232 painCave.severity = OPENMD_WARNING;
233 simError();
234
235 // The old file version tethered the shifting behavior to the
236 // electrostaticSummationMethod keyword.
237
238 if (simParams_->haveElectrostaticSummationMethod()) {
239 string myMethod = simParams_->getElectrostaticSummationMethod();
240 toUpper(myMethod);
241
242 if (myMethod == "SHIFTED_POTENTIAL") {
243 cutoffMethod_ = SHIFTED_POTENTIAL;
244 } else if (myMethod == "SHIFTED_FORCE") {
245 cutoffMethod_ = SHIFTED_FORCE;
246 } else if (myMethod == "TAYLOR_SHIFTED") {
247 cutoffMethod_ = TAYLOR_SHIFTED;
248 } else if (myMethod == "EWALD_FULL") {
249 cutoffMethod_ = EWALD_FULL;
250 useSurfaceTerm_ = true;
251 }
252
253 if (simParams_->haveSwitchingRadius())
254 rSwitch_ = simParams_->getSwitchingRadius();
255
256 if (myMethod == "SHIFTED_POTENTIAL" || myMethod == "SHIFTED_FORCE" ||
257 myMethod == "TAYLOR_SHIFTED" || myMethod == "EWALD_FULL") {
258 if (simParams_->haveSwitchingRadius()) {
259 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
260 "ForceManager::setupCutoffs : DEPRECATED ERROR MESSAGE\n"
261 "\tA value was set for the switchingRadius\n"
262 "\teven though the electrostaticSummationMethod was\n"
263 "\tset to %s\n",
264 myMethod.c_str());
265 painCave.severity = OPENMD_WARNING;
266 painCave.isFatal = 1;
267 simError();
268 }
269 }
270 if (abs(rCut_ - rSwitch_) < 0.0001) {
271 if (cutoffMethod_ == SHIFTED_FORCE) {
272 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
273 "ForceManager::setupCutoffs : DEPRECATED BEHAVIOR\n"
274 "\tcutoffRadius and switchingRadius are set to the\n"
275 "\tsame value. OpenMD will use shifted force\n"
276 "\tpotentials instead of switching functions.\n");
277 painCave.isFatal = 0;
278 painCave.severity = OPENMD_WARNING;
279 simError();
280 } else {
281 cutoffMethod_ = SHIFTED_POTENTIAL;
282 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
283 "ForceManager::setupCutoffs : DEPRECATED BEHAVIOR\n"
284 "\tcutoffRadius and switchingRadius are set to the\n"
285 "\tsame value. OpenMD will use shifted potentials\n"
286 "\tinstead of switching functions.\n");
287 painCave.isFatal = 0;
288 painCave.severity = OPENMD_WARNING;
289 simError();
290 }
291 }
292 }
293 }
294 }
295
296 // create the switching function object:
297
298 switcher_ = new SwitchingFunction();
299
300 if (cutoffMethod_ == SWITCHED) {
301 if (simParams_->haveSwitchingRadius()) {
302 rSwitch_ = simParams_->getSwitchingRadius();
303 if (rSwitch_ > rCut_) {
304 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
305 "ForceManager::setupCutoffs: switchingRadius (%f) is larger "
306 "than the cutoffRadius(%f)\n",
307 rSwitch_, rCut_);
308 painCave.isFatal = 1;
309 painCave.severity = OPENMD_ERROR;
310 simError();
311 }
312 } else {
313 rSwitch_ = 0.85 * rCut_;
314 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
315 "ForceManager::setupCutoffs: No value was set for the "
316 "switchingRadius.\n"
317 "\tOpenMD will use a default value of 85 percent of the "
318 "cutoffRadius.\n"
319 "\tswitchingRadius = %f. for this simulation\n",
320 rSwitch_);
321 painCave.isFatal = 0;
322 painCave.severity = OPENMD_WARNING;
323 simError();
324 }
325 } else {
326 if (mdFileVersion > 1) {
327 // throw an error if we define a switching radius and don't need one.
328 // older file versions should not do this.
329 if (simParams_->haveSwitchingRadius()) {
330 map<string, CutoffMethod>::const_iterator it;
331 string theMeth;
332 for (it = stringToCutoffMethod.begin();
333 it != stringToCutoffMethod.end(); ++it) {
334 if (it->second == cutoffMethod_) {
335 theMeth = it->first;
336 break;
337 }
338 }
339 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
340 "ForceManager::setupCutoffs: the cutoffMethod (%s)\n"
341 "\tis not set to SWITCHED, so switchingRadius value\n"
342 "\twill be ignored for this simulation\n",
343 theMeth.c_str());
344 painCave.isFatal = 0;
345 painCave.severity = OPENMD_WARNING;
346 simError();
347 }
348 }
349 rSwitch_ = rCut_;
350 }
351
352 // Default to cubic switching function.
353 sft_ = cubic;
354 if (simParams_->haveSwitchingFunctionType()) {
355 string funcType = simParams_->getSwitchingFunctionType();
356 toUpper(funcType);
357 if (funcType == "CUBIC") {
358 sft_ = cubic;
359 } else {
360 if (funcType == "FIFTH_ORDER_POLYNOMIAL") {
361 sft_ = fifth_order_poly;
362 } else {
363 // throw error
364 snprintf(
365 painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
366 "ForceManager::setupSwitching : Unknown switchingFunctionType. "
367 "(Input "
368 "file specified %s .)\n"
369 "\tswitchingFunctionType must be one of: "
370 "\"cubic\" or \"fifth_order_polynomial\".",
371 funcType.c_str());
372 painCave.isFatal = 1;
373 painCave.severity = OPENMD_ERROR;
374 simError();
375 }
376 }
377 }
378 switcher_->setSwitchType(sft_);
379 switcher_->setSwitch(rSwitch_, rCut_);
380 }
381
382 void ForceManager::initialize() {
383 if (!info_->isTopologyDone()) {
384 info_->update();
385 interactionMan_->setSimInfo(info_);
386 interactionMan_->initialize();
387
388 useSurfaceTerm_ = false; // default, can be set by Ewald or directly
389 useSlabGeometry_ = false;
390
391 //! We want to delay the cutoffs until after the interaction
392 //! manager has set up the atom-atom interactions so that we can
393 //! query them for suggested cutoff values
394 setupCutoffs();
395
396 info_->prepareTopology();
397
398 doParticlePot_ = info_->getSimParams()->getOutputParticlePotential();
399 doHeatFlux_ = info_->getSimParams()->getPrintHeatFlux();
400 if (doHeatFlux_) doParticlePot_ = true;
401
402 doElectricField_ = info_->getSimParams()->getOutputElectricField();
403 doElectricField_ |=
404 info_->getSimParams()->getRNEMDParameters()->requiresElectricField();
405 doSitePotential_ = info_->getSimParams()->getOutputSitePotential();
406 if (info_->getSimParams()->haveUseSurfaceTerm() &&
407 info_->getSimParams()->getUseSurfaceTerm()) {
408 if (info_->usesElectrostaticAtoms()) {
409 useSurfaceTerm_ = info_->getSimParams()->getUseSurfaceTerm();
410 if (info_->getSimParams()->haveUseSlabGeometry()) {
411 useSlabGeometry_ = info_->getSimParams()->getUseSlabGeometry();
412
413 string axis =
414 toUpperCopy(info_->getSimParams()->getPrivilegedAxis());
415 if (axis.compare("X") == 0)
416 axis_ = 0;
417 else if (axis.compare("Y") == 0)
418 axis_ = 1;
419 else
420 axis_ = 2;
421 }
422 } else {
423 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
424 "ForceManager::initialize : useSurfaceTerm was set true\n"
425 "\tbut no electrostatic atoms are present. OpenMD will\n"
426 "\tignore this setting.\n");
427 painCave.isFatal = 0;
428 painCave.severity = OPENMD_WARNING;
429 simError();
430 }
431 }
432 }
433
434 ForceFieldOptions& fopts = forceField_->getForceFieldOptions();
435
436 //! Force fields can set options on how to scale van der Waals and
437 //! electrostatic interactions for atoms connected via bonds, bends
438 //! and torsions in this case the topological distance between
439 //! atoms is:
440 //! 0 = topologically unconnected
441 //! 1 = bonded together
442 //! 2 = connected via a bend
443 //! 3 = connected via a torsion
444
445 vdwScale_.resize(4);
446 fill(vdwScale_.begin(), vdwScale_.end(), 0.0);
447
448 electrostaticScale_.resize(4);
449 fill(electrostaticScale_.begin(), electrostaticScale_.end(), 0.0);
450
451 vdwScale_[0] = 1.0;
452 vdwScale_[1] = fopts.getvdw12scale();
453 vdwScale_[2] = fopts.getvdw13scale();
454 vdwScale_[3] = fopts.getvdw14scale();
455
456 electrostaticScale_[0] = 1.0;
457 electrostaticScale_[1] = fopts.getelectrostatic12scale();
458 electrostaticScale_[2] = fopts.getelectrostatic13scale();
459 electrostaticScale_[3] = fopts.getelectrostatic14scale();
460
461 // Initialize the perturbations
462 if (info_->getSimParams()->haveUniformField()) {
463 UniformField* eField = new UniformField(info_);
464 forceModifiers_.push_back(eField);
465 }
466
467 if (info_->getSimParams()->haveMagneticField()) {
468 MagneticField* mField = new MagneticField(info_);
469 forceModifiers_.push_back(mField);
470 }
471
472 if (info_->getSimParams()->haveUniformGradientStrength() ||
473 info_->getSimParams()->haveUniformGradientDirection1() ||
474 info_->getSimParams()->haveUniformGradientDirection2()) {
475 UniformGradient* eGrad = new UniformGradient(info_);
476 forceModifiers_.push_back(eGrad);
477 }
478
479 if (info_->getSimParams()->getLightParameters()->getUseLight()) {
480 Perturbations::Light* light = new Perturbations::Light(info_);
481 forceModifiers_.push_back(light);
482 }
483
484 // Initialize the force modifiers (order matters)
485 if (info_->getSimParams()->getUseThermodynamicIntegration()) {
486 ThermoIntegrationForceModifier* thermoInt =
487 new ThermoIntegrationForceModifier(info_);
488 forceModifiers_.push_back(thermoInt);
489 } else if (info_->getSimParams()->getUseRestraints()) {
490 RestraintForceModifier* restraint = new RestraintForceModifier(info_);
491 forceModifiers_.push_back(restraint);
492 }
493
494 std::string ensembleParam =
495 toUpperCopy(info_->getSimParams()->getEnsemble());
496
497 if (ensembleParam == "LHULL" || ensembleParam == "LANGEVINHULL" ||
498 ensembleParam == "SMIPD") {
499#if defined(HAVE_QHULL)
500 LangevinHullForceModifier* langevinHullFM =
501 new LangevinHullForceModifier(info_);
502 forceModifiers_.push_back(langevinHullFM);
503#else
504 snprintf(
505 painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
506 "ForceManager: Cannot use the LangevinHull ensembles without qhull.\n"
507 "\tPlease rebuild OpenMD with qhull enabled.");
508 painCave.severity = OPENMD_ERROR;
509 painCave.isFatal = 1;
510 simError();
511#endif
512 } else if (ensembleParam == "LANGEVINDYNAMICS" || ensembleParam == "LD") {
513 LDForceModifier* langevinDynamicsFM = new LDForceModifier(info_);
514 forceModifiers_.push_back(langevinDynamicsFM);
515 } else if (ensembleParam == "RPYDYNAMICS" || ensembleParam == "RPY"
516 || ensembleParam =="LHD" || ensembleParam =="LANGEVINHYDRODYNAMICS" ) {
517 LHDForceModifier* lhdFM = new LHDForceModifier(info_);
518 forceModifiers_.push_back(lhdFM);
519 }
520
521 if (info_->getSimParams()->getNZconsStamps() > 0) {
522 info_->setNZconstraint(info_->getSimParams()->getNZconsStamps());
523 ZConstraintForceModifier* zCons = new ZConstraintForceModifier(info_);
524 forceModifiers_.push_back(zCons);
525 }
526
527 usePeriodicBoundaryConditions_ =
528 info_->getSimParams()->getUsePeriodicBoundaryConditions();
529
530 fDecomp_->distributeInitialData();
531
532 doPotentialSelection_ = false;
533 if (info_->getSimParams()->havePotentialSelection()) {
534 doPotentialSelection_ = true;
535 selectionScript_ = info_->getSimParams()->getPotentialSelection();
536 evaluator_.loadScriptString(selectionScript_);
537 if (!evaluator_.isDynamic()) {
538 seleMan_.setSelectionSet(evaluator_.evaluate());
539 }
540 }
541
542 initialized_ = true;
543 }
544
545 void ForceManager::calcForces() {
546 if (!initialized_) initialize();
547
548 preCalculation();
549 shortRangeInteractions();
550 if (info_->getSimParams()->getSkipPairLoop() == false)
551 longRangeInteractions();
552 postCalculation();
553 }
554
555 void ForceManager::preCalculation() {
556 SimInfo::MoleculeIterator mi;
557 Molecule* mol;
558 Molecule::AtomIterator ai;
559 Atom* atom;
560 Molecule::RigidBodyIterator rbIter;
561 RigidBody* rb;
562 Molecule::CutoffGroupIterator ci;
563 CutoffGroup* cg;
564
565 // forces and potentials are zeroed here, before any are
566 // accumulated.
567
568 Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot();
569
570 fDecomp_->setSnapshot(snap);
571
572 snap->setBondPotential(0.0);
573 snap->setBendPotential(0.0);
574 snap->setTorsionPotential(0.0);
575 snap->setInversionPotential(0.0);
576
577 potVec zeroPot(0.0);
578 snap->setLongRangePotentials(zeroPot);
579
580 snap->setExcludedPotentials(zeroPot);
581 if (doPotentialSelection_) snap->setSelectionPotentials(zeroPot);
582
583 snap->setSelfPotentials(0.0);
584 snap->setRestraintPotential(0.0);
585 snap->setRawPotential(0.0);
586
587 for (mol = info_->beginMolecule(mi); mol != NULL;
588 mol = info_->nextMolecule(mi)) {
589 for (atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
590 atom->zeroForcesAndTorques();
591 }
592
593 // change the positions of atoms which belong to the rigidbodies
594 for (rb = mol->beginRigidBody(rbIter); rb != NULL;
595 rb = mol->nextRigidBody(rbIter)) {
596 rb->zeroForcesAndTorques();
597 }
598
599 if (info_->getNGlobalCutoffGroups() != info_->getNGlobalAtoms()) {
600 for (cg = mol->beginCutoffGroup(ci); cg != NULL;
601 cg = mol->nextCutoffGroup(ci)) {
602 // calculate the center of mass of cutoff group
603 cg->updateCOM();
604 }
605 }
606 }
607
608 // Zero out the virial tensor
609 virialTensor *= 0.0;
610 // Zero out the heatFlux
611 fDecomp_->setHeatFlux(Vector3d(0.0));
612
613 if (doPotentialSelection_) {
614 if (evaluator_.isDynamic()) {
615 seleMan_.setSelectionSet(evaluator_.evaluate());
616 }
617 }
618 }
619
620 void ForceManager::shortRangeInteractions() {
621 Molecule* mol;
622 RigidBody* rb;
623 Bond* bond;
624 Bend* bend;
625 Torsion* torsion;
626 Inversion* inversion;
627 SimInfo::MoleculeIterator mi;
628 Molecule::RigidBodyIterator rbIter;
629 Molecule::BondIterator bondIter;
630
631 Molecule::BendIterator bendIter;
632 Molecule::TorsionIterator torsionIter;
633 Molecule::InversionIterator inversionIter;
634 RealType bondPotential = 0.0;
635 RealType bendPotential = 0.0;
636 RealType torsionPotential = 0.0;
637 RealType inversionPotential = 0.0;
638 potVec selectionPotential(0.0);
639
640 // calculate short range interactions
641 for (mol = info_->beginMolecule(mi); mol != NULL;
642 mol = info_->nextMolecule(mi)) {
643 // change the positions of atoms which belong to the rigidbodies
644 for (rb = mol->beginRigidBody(rbIter); rb != NULL;
645 rb = mol->nextRigidBody(rbIter)) {
646 rb->updateAtoms();
647 }
648
649 for (bond = mol->beginBond(bondIter); bond != NULL;
650 bond = mol->nextBond(bondIter)) {
651 bond->calcForce(doParticlePot_);
652 bondPotential += bond->getPotential();
653 if (doPotentialSelection_) {
654 if (seleMan_.isSelected(bond->getAtomA()) ||
655 seleMan_.isSelected(bond->getAtomB())) {
656 selectionPotential[BONDED_FAMILY] += bond->getPotential();
657 }
658 }
659 }
660
661 for (bend = mol->beginBend(bendIter); bend != NULL;
662 bend = mol->nextBend(bendIter)) {
663 RealType angle;
664 bend->calcForce(angle, doParticlePot_);
665 RealType currBendPot = bend->getPotential();
666
667 bendPotential += currBendPot;
668 map<Bend*, BendDataSet>::iterator i = bendDataSets.find(bend);
669 if (i == bendDataSets.end()) {
670 BendDataSet dataSet;
671 dataSet.prev.angle = dataSet.curr.angle = angle;
672 dataSet.prev.potential = dataSet.curr.potential = currBendPot;
673 dataSet.deltaV = 0.0;
674 bendDataSets.insert(
675 map<Bend*, BendDataSet>::value_type(bend, dataSet));
676 } else {
677 i->second.prev.angle = i->second.curr.angle;
678 i->second.prev.potential = i->second.curr.potential;
679 i->second.curr.angle = angle;
680 i->second.curr.potential = currBendPot;
681 i->second.deltaV =
682 fabs(i->second.curr.potential - i->second.prev.potential);
683 }
684 if (doPotentialSelection_) {
685 if (seleMan_.isSelected(bend->getAtomA()) ||
686 seleMan_.isSelected(bend->getAtomB()) ||
687 seleMan_.isSelected(bend->getAtomC())) {
688 selectionPotential[BONDED_FAMILY] += currBendPot;
689 }
690 }
691 }
692
693 for (torsion = mol->beginTorsion(torsionIter); torsion != NULL;
694 torsion = mol->nextTorsion(torsionIter)) {
695 RealType angle;
696 torsion->calcForce(angle, doParticlePot_);
697 RealType currTorsionPot = torsion->getPotential();
698 torsionPotential += currTorsionPot;
699 map<Torsion*, TorsionDataSet>::iterator i =
700 torsionDataSets.find(torsion);
701 if (i == torsionDataSets.end()) {
702 TorsionDataSet dataSet;
703 dataSet.prev.angle = dataSet.curr.angle = angle;
704 dataSet.prev.potential = dataSet.curr.potential = currTorsionPot;
705 dataSet.deltaV = 0.0;
706 torsionDataSets.insert(
707 map<Torsion*, TorsionDataSet>::value_type(torsion, dataSet));
708 } else {
709 i->second.prev.angle = i->second.curr.angle;
710 i->second.prev.potential = i->second.curr.potential;
711 i->second.curr.angle = angle;
712 i->second.curr.potential = currTorsionPot;
713 i->second.deltaV =
714 fabs(i->second.curr.potential - i->second.prev.potential);
715 }
716 if (doPotentialSelection_) {
717 if (seleMan_.isSelected(torsion->getAtomA()) ||
718 seleMan_.isSelected(torsion->getAtomB()) ||
719 seleMan_.isSelected(torsion->getAtomC()) ||
720 seleMan_.isSelected(torsion->getAtomD())) {
721 selectionPotential[BONDED_FAMILY] += currTorsionPot;
722 }
723 }
724 }
725
726 for (inversion = mol->beginInversion(inversionIter); inversion != NULL;
727 inversion = mol->nextInversion(inversionIter)) {
728 RealType angle;
729 inversion->calcForce(angle, doParticlePot_);
730 RealType currInversionPot = inversion->getPotential();
731 inversionPotential += currInversionPot;
732 map<Inversion*, InversionDataSet>::iterator i =
733 inversionDataSets.find(inversion);
734 if (i == inversionDataSets.end()) {
735 InversionDataSet dataSet;
736 dataSet.prev.angle = dataSet.curr.angle = angle;
737 dataSet.prev.potential = dataSet.curr.potential = currInversionPot;
738 dataSet.deltaV = 0.0;
739 inversionDataSets.insert(
740 map<Inversion*, InversionDataSet>::value_type(inversion,
741 dataSet));
742 } else {
743 i->second.prev.angle = i->second.curr.angle;
744 i->second.prev.potential = i->second.curr.potential;
745 i->second.curr.angle = angle;
746 i->second.curr.potential = currInversionPot;
747 i->second.deltaV =
748 fabs(i->second.curr.potential - i->second.prev.potential);
749 }
750 if (doPotentialSelection_) {
751 if (seleMan_.isSelected(inversion->getAtomA()) ||
752 seleMan_.isSelected(inversion->getAtomB()) ||
753 seleMan_.isSelected(inversion->getAtomC()) ||
754 seleMan_.isSelected(inversion->getAtomD())) {
755 selectionPotential[BONDED_FAMILY] += currInversionPot;
756 }
757 }
758 }
759 }
760
761#ifdef IS_MPI
762 // Collect from all nodes. This should eventually be moved into a
763 // SystemDecomposition, but this is a better place than in
764 // Thermo to do the collection.
765
766 MPI_Allreduce(MPI_IN_PLACE, &bondPotential, 1, MPI_REALTYPE, MPI_SUM,
767 MPI_COMM_WORLD);
768 MPI_Allreduce(MPI_IN_PLACE, &bendPotential, 1, MPI_REALTYPE, MPI_SUM,
769 MPI_COMM_WORLD);
770 MPI_Allreduce(MPI_IN_PLACE, &torsionPotential, 1, MPI_REALTYPE, MPI_SUM,
771 MPI_COMM_WORLD);
772 MPI_Allreduce(MPI_IN_PLACE, &inversionPotential, 1, MPI_REALTYPE, MPI_SUM,
773 MPI_COMM_WORLD);
774 MPI_Allreduce(MPI_IN_PLACE, &selectionPotential[BONDED_FAMILY], 1,
775 MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
776#endif
777
778 Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
779 fDecomp_->setSnapshot(curSnapshot);
780
781 curSnapshot->setBondPotential(bondPotential);
782 curSnapshot->setBendPotential(bendPotential);
783 curSnapshot->setTorsionPotential(torsionPotential);
784 curSnapshot->setInversionPotential(inversionPotential);
785 curSnapshot->setSelectionPotentials(selectionPotential);
786
787 // RealType shortRangePotential = bondPotential + bendPotential +
788 // torsionPotential + inversionPotential;
789
790 // curSnapshot->setShortRangePotential(shortRangePotential);
791 }
792
793 void ForceManager::longRangeInteractions() {
794 Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
795 fDecomp_->setSnapshot(curSnapshot);
796
797 DataStorage* config = &(curSnapshot->atomData);
798 DataStorage* cgConfig = &(curSnapshot->cgData);
799
800 // calculate the center of mass of cutoff group
801
802 SimInfo::MoleculeIterator mi;
803 Molecule* mol;
804 Molecule::CutoffGroupIterator ci;
805 CutoffGroup* cg;
806
807 if (info_->getNCutoffGroups() != info_->getNAtoms()) {
808 for (mol = info_->beginMolecule(mi); mol != NULL;
809 mol = info_->nextMolecule(mi)) {
810 for (cg = mol->beginCutoffGroup(ci); cg != NULL;
811 cg = mol->nextCutoffGroup(ci)) {
812 cg->updateCOM();
813 }
814 }
815 } else {
816 // center of mass of the group is the same as position of the atom
817 // if cutoff group does not exist
818 cgConfig->position = config->position;
819 cgConfig->velocity = config->velocity;
820 }
821
822 fDecomp_->zeroWorkArrays();
823 fDecomp_->distributeData();
824
825 int cg1, cg2, atom1, atom2;
826 Vector3d d_grp, dag, gvel2, vel2;
827 RealType rgrpsq, rgrp;
828 RealType vij(0.0);
829 Vector3d fij, fg, f1;
830 bool in_switching_region;
831 RealType dswdr, swderiv;
832 vector<int> atomListColumn, atomListRow;
833
834 potVec longRangePotential(0.0);
835 potVec selfPotential(0.0);
836 RealType reciprocalPotential(0.0);
837 RealType surfacePotential(0.0);
838 potVec selectionPotential(0.0);
839
840 RealType mf;
841 bool newAtom1;
842 int gid1, gid2;
843
844 vector<int>::iterator ia, jb;
845
846 int loopStart, loopEnd;
847
848 idat.rcut = rCut_;
849 idat.shiftedPot = (cutoffMethod_ == SHIFTED_POTENTIAL) ? true : false;
850 idat.shiftedForce =
851 (cutoffMethod_ == SHIFTED_FORCE || cutoffMethod_ == TAYLOR_SHIFTED) ?
852 true :
853 false;
854 idat.doParticlePot = doParticlePot_;
855 idat.doElectricField = doElectricField_;
856 idat.doSitePotential = doSitePotential_;
857 sdat.doParticlePot = doParticlePot_;
858
859 loopEnd = PAIR_LOOP;
860 if (info_->requiresPrepair()) {
861 loopStart = PREPAIR_LOOP;
862 } else {
863 loopStart = PAIR_LOOP;
864 }
865 for (int iLoop = loopStart; iLoop <= loopEnd; iLoop++) {
866 if (iLoop == loopStart) {
867 bool update_nlist = fDecomp_->checkNeighborList(savedPositions_);
868
869 if (update_nlist) {
870 if (!usePeriodicBoundaryConditions_)
871 Mat3x3d bbox = thermo->getBoundingBox();
872 fDecomp_->buildNeighborList(neighborList_, point_, savedPositions_);
873 }
874 }
875
876 for (cg1 = 0; cg1 < int(point_.size()) - 1; cg1++) {
877 atomListRow = fDecomp_->getAtomsInGroupRow(cg1);
878 newAtom1 = true;
879
880 for (int m2 = point_[cg1]; m2 < point_[cg1 + 1]; m2++) {
881 cg2 = neighborList_[m2];
882
883 d_grp = fDecomp_->getIntergroupVector(cg1, cg2);
884
885 // already wrapped in the getIntergroupVector call:
886 // curSnapshot->wrapVector(d_grp);
887 rgrpsq = d_grp.lengthSquare();
888
889 if (rgrpsq < rCutSq_) {
890 if (iLoop == PAIR_LOOP) {
891 vij = 0.0;
892 fij.zero();
893 idat.eField1.zero();
894 idat.eField2.zero();
895 idat.sPot1 = 0.0;
896 idat.sPot2 = 0.0;
897 }
898
899 in_switching_region =
900 switcher_->getSwitch(rgrpsq, idat.sw, dswdr, rgrp);
901
902 atomListColumn = fDecomp_->getAtomsInGroupColumn(cg2);
903
904 if (doHeatFlux_) gvel2 = fDecomp_->getGroupVelocityColumn(cg2);
905
906 for (ia = atomListRow.begin(); ia != atomListRow.end(); ++ia) {
907 atom1 = (*ia);
908 if (atomListRow.size() > 1) newAtom1 = true;
909
910 if (doPotentialSelection_) {
911 gid1 = fDecomp_->getGlobalIDRow(atom1);
912 idat.isSelected = seleMan_.isGlobalIDSelected(gid1);
913 }
914
915 for (jb = atomListColumn.begin(); jb != atomListColumn.end();
916 ++jb) {
917 atom2 = (*jb);
918
919 if (doPotentialSelection_) {
920 gid2 = fDecomp_->getGlobalIDCol(atom2);
921 idat.isSelected |= seleMan_.isGlobalIDSelected(gid2);
922 }
923
924 if (!fDecomp_->skipAtomPair(atom1, atom2, cg1, cg2)) {
925 idat.vpair = 0.0;
926 idat.pot = 0.0;
927 idat.excludedPot = 0.0;
928 idat.selePot = 0.0;
929 idat.f1.zero();
930 idat.dVdFQ1 = 0.0;
931 idat.dVdFQ2 = 0.0;
932
933 fDecomp_->fillInteractionData(idat, atom1, atom2, newAtom1);
934 newAtom1 = false;
935
936 idat.topoDist =
937 fDecomp_->getTopologicalDistance(atom1, atom2);
938 idat.vdwMult = vdwScale_[idat.topoDist];
939 idat.electroMult = electrostaticScale_[idat.topoDist];
940
941 if (atomListRow.size() == 1 && atomListColumn.size() == 1) {
942 idat.d = d_grp;
943 idat.r2 = rgrpsq;
944 if (doHeatFlux_) vel2 = gvel2;
945 } else {
946 idat.d = fDecomp_->getInteratomicVector(atom1, atom2);
947 curSnapshot->wrapVector(idat.d);
948 idat.r2 = idat.d.lengthSquare();
949 if (doHeatFlux_)
950 vel2 = fDecomp_->getAtomVelocityColumn(atom2);
951 }
952
953 idat.rij = sqrt(idat.r2);
954
955 if (iLoop == PREPAIR_LOOP) {
956 interactionMan_->doPrePair(idat);
957 fDecomp_->unpackPrePairData(idat, atom1, atom2);
958 } else {
959 interactionMan_->doPair(idat);
960 fDecomp_->unpackInteractionData(idat, atom1, atom2);
961 vij += idat.vpair;
962 fij += idat.f1;
963 virialTensor -= outProduct(idat.d, idat.f1);
964 if (doHeatFlux_)
965 fDecomp_->addToHeatFlux(idat.d * dot(idat.f1, vel2));
966 }
967 }
968 }
969 }
970
971 if (iLoop == PAIR_LOOP) {
972 if (in_switching_region) {
973 swderiv = vij * dswdr / rgrp;
974 fg = swderiv * d_grp;
975 fij += fg;
976
977 if (atomListRow.size() == 1 && atomListColumn.size() == 1) {
978 if (!fDecomp_->skipAtomPair(atomListRow[0], atomListColumn[0],
979 cg1, cg2)) {
980 virialTensor -= outProduct(idat.d, fg);
981 if (doHeatFlux_)
982 fDecomp_->addToHeatFlux(idat.d * dot(fg, vel2));
983 }
984 }
985
986 for (ia = atomListRow.begin(); ia != atomListRow.end(); ++ia) {
987 atom1 = (*ia);
988 mf = fDecomp_->getMassFactorRow(atom1);
989 // fg is the force on atom ia due to cutoff group's
990 // presence in switching region
991 fg = swderiv * d_grp * mf;
992 fDecomp_->addForceToAtomRow(atom1, fg);
993 if (atomListRow.size() > 1) {
994 if (info_->usesAtomicVirial()) {
995 // find the distance between the atom
996 // and the center of the cutoff group:
997 dag = fDecomp_->getAtomToGroupVectorRow(atom1, cg1);
998 virialTensor -= outProduct(dag, fg);
999 if (doHeatFlux_)
1000 fDecomp_->addToHeatFlux(dag * dot(fg, vel2));
1001 }
1002 }
1003 }
1004 for (jb = atomListColumn.begin(); jb != atomListColumn.end();
1005 ++jb) {
1006 atom2 = (*jb);
1007 mf = fDecomp_->getMassFactorColumn(atom2);
1008 // fg is the force on atom jb due to cutoff group's
1009 // presence in switching region
1010 fg = -swderiv * d_grp * mf;
1011 fDecomp_->addForceToAtomColumn(atom2, fg);
1012
1013 if (atomListColumn.size() > 1) {
1014 if (info_->usesAtomicVirial()) {
1015 // find the distance between the atom
1016 // and the center of the cutoff group:
1017 dag = fDecomp_->getAtomToGroupVectorColumn(atom2, cg2);
1018 virialTensor -= outProduct(dag, fg);
1019 if (doHeatFlux_)
1020 fDecomp_->addToHeatFlux(dag * dot(fg, vel2));
1021 }
1022 }
1023 }
1024 }
1025 // if (!info_->usesAtomicVirial()) {
1026 // virialTensor -= outProduct(d_grp, fij);
1027 // if (doHeatFlux_)
1028 // fDecomp_->addToHeatFlux( d_grp * dot(fij, vel2));
1029 //}
1030 }
1031 }
1032 }
1033 }
1034
1035 if (iLoop == PREPAIR_LOOP) {
1036 if (info_->requiresPrepair()) {
1037 fDecomp_->collectIntermediateData();
1038
1039 for (unsigned int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) {
1040 if (doPotentialSelection_) {
1041 gid1 = fDecomp_->getGlobalID(atom1);
1042 sdat.isSelected = seleMan_.isGlobalIDSelected(gid1);
1043 }
1044 fDecomp_->fillPreForceData(sdat, atom1);
1045 interactionMan_->doPreForce(sdat);
1046 fDecomp_->unpackPreForceData(sdat, atom1);
1047 }
1048
1049 fDecomp_->distributeIntermediateData();
1050 }
1051 }
1052 }
1053
1054 // collects pairwise information
1055 fDecomp_->collectData();
1056 if (cutoffMethod_ == EWALD_FULL) {
1057 interactionMan_->doReciprocalSpaceSum(reciprocalPotential);
1058 curSnapshot->setReciprocalPotential(reciprocalPotential);
1059 }
1060
1061 if (useSurfaceTerm_) {
1062 interactionMan_->doSurfaceTerm(useSlabGeometry_, axis_, surfacePotential);
1063 curSnapshot->setSurfacePotential(surfacePotential);
1064 }
1065
1066 if (info_->requiresSelfCorrection()) {
1067 for (unsigned int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) {
1068 if (doPotentialSelection_) {
1069 gid1 = fDecomp_->getGlobalID(atom1);
1070 sdat.isSelected = seleMan_.isGlobalIDSelected(gid1);
1071 }
1072 fDecomp_->fillSelfData(sdat, atom1);
1073 interactionMan_->doSelfCorrection(sdat);
1074 fDecomp_->unpackSelfData(sdat, atom1);
1075 }
1076 }
1077
1078 // collects single-atom information
1079 fDecomp_->collectSelfData();
1080
1081 longRangePotential = fDecomp_->getPairwisePotential();
1082 curSnapshot->setLongRangePotentials(longRangePotential);
1083
1084 selfPotential = fDecomp_->getSelfPotential();
1085 curSnapshot->setSelfPotentials(selfPotential);
1086
1087 curSnapshot->setExcludedPotentials(fDecomp_->getExcludedSelfPotential() +
1088 fDecomp_->getExcludedPotential());
1089
1090 if (doPotentialSelection_) {
1091 selectionPotential = curSnapshot->getSelectionPotentials();
1092 selectionPotential += fDecomp_->getSelectedSelfPotential();
1093 selectionPotential += fDecomp_->getSelectedPotential();
1094 curSnapshot->setSelectionPotentials(selectionPotential);
1095 }
1096 }
1097
1098 void ForceManager::postCalculation() {
1099 for (auto& forceModifier : forceModifiers_)
1100 forceModifier->modifyForces();
1101
1102 // Modify the rigid bodies in response to the applied force modifications
1103 SimInfo::MoleculeIterator mi;
1104 Molecule* mol;
1105 Molecule::RigidBodyIterator rbIter;
1106 RigidBody* rb;
1107 Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
1108
1109 // Collect the atomic forces onto rigid bodies
1110 for (mol = info_->beginMolecule(mi); mol != NULL;
1111 mol = info_->nextMolecule(mi)) {
1112 for (rb = mol->beginRigidBody(rbIter); rb != NULL;
1113 rb = mol->nextRigidBody(rbIter)) {
1114 Mat3x3d rbTau = rb->calcForcesAndTorquesAndVirial();
1115 virialTensor += rbTau;
1116 }
1117 }
1118
1119#ifdef IS_MPI
1120 MPI_Allreduce(MPI_IN_PLACE, virialTensor.getArrayPointer(), 9, MPI_REALTYPE,
1121 MPI_SUM, MPI_COMM_WORLD);
1122#endif
1123 curSnapshot->setVirialTensor(virialTensor);
1124
1125 /*
1126 if (info_->getSimParams()->getUseLongRangeCorrections()) {
1127 RealType vol = curSnapshot->getVolume();
1128 RealType Elrc(0.0);
1129 RealType Wlrc(0.0);
1130
1131 AtomTypeSet::iterator i;
1132 AtomTypeSet::iterator j;
1133
1134 RealType n_i, n_j;
1135 RealType rho_i, rho_j;
1136 pair<RealType, RealType> LRI;
1137
1138 for (i = atomTypes_.begin(); i != atomTypes_.end(); ++i) {
1139 n_i = RealType(info_->getGlobalCountOfType(*i));
1140 rho_i = n_i / vol;
1141 for (j = atomTypes_.begin(); j != atomTypes_.end(); ++j) {
1142 n_j = RealType(info_->getGlobalCountOfType(*j));
1143 rho_j = n_j / vol;
1144
1145 LRI = interactionMan_->getLongRangeIntegrals( (*i), (*j) );
1146
1147 Elrc += n_i * rho_j * LRI.first;
1148 Wlrc -= rho_i * rho_j * LRI.second;
1149 }
1150 }
1151 Elrc *= 2.0 * Constants::PI;
1152 Wlrc *= 2.0 * Constants::PI;
1153
1154 RealType lrp = curSnapshot->getLongRangePotential();
1155 curSnapshot->setLongRangePotential(lrp + Elrc);
1156 virialTensor += Wlrc * SquareMatrix3<RealType>::identity();
1157 curSnapshot->setVirialTensor(virialTensor);
1158 }
1159 */
1160
1161 // Sanity check on PE - don't wait for StatWriter or DumpWriter.
1162 RealType pe = curSnapshot->getPotentialEnergy();
1163 if (std::isinf(pe) || std::isnan(pe)) {
1164 snprintf(
1165 painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
1166 "ForceManager detected a numerical error in the potential energy.\n"
1167 "\tStopping the simulation.");
1168 painCave.isFatal = 1;
1169 simError();
1170 }
1171 Vector3d f;
1172 Vector3d t;
1173
1174 Molecule::IntegrableObjectIterator ii;
1175 StuntDouble* sd;
1176
1177 for (mol = info_->beginMolecule(mi); mol != NULL;
1178 mol = info_->nextMolecule(mi)) {
1179 for (sd = mol->beginIntegrableObject(ii); sd != NULL;
1180 sd = mol->nextIntegrableObject(ii)) {
1181 f = sd->getFrc();
1182 if (std::isinf(f[0]) || std::isnan(f[0]) || std::isinf(f[1]) ||
1183 std::isnan(f[1]) || std::isinf(f[2]) || std::isnan(f[2])) {
1184 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
1185 "ForceManager detected a numerical error in the force"
1186 " for object %d",
1187 sd->getGlobalIndex());
1188 painCave.isFatal = 1;
1189 simError();
1190 }
1191 if (sd->isDirectional()) {
1192 t = sd->getTrq();
1193 if (std::isinf(t[0]) || std::isnan(t[0]) || std::isinf(t[1]) ||
1194 std::isnan(t[1]) || std::isinf(t[2]) || std::isnan(t[2])) {
1195 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
1196 "ForceManager detected a numerical error in the torque"
1197 " for object %d",
1198 sd->getGlobalIndex());
1199 painCave.isFatal = 1;
1200 simError();
1201 }
1202 }
1203 }
1204 }
1205 errorCheckPoint();
1206 }
1207
1208 void ForceManager::calcSelectedForces(Molecule* mol1, Molecule* mol2) {
1209 if (!initialized_) initialize();
1210
1211 selectedPreCalculation(mol1, mol2);
1212 selectedShortRangeInteractions(mol1, mol2);
1213 selectedLongRangeInteractions(mol1, mol2);
1214 selectedPostCalculation(mol1, mol2);
1215 }
1216
1217 void ForceManager::selectedPreCalculation(Molecule* mol1, Molecule* mol2) {
1218 SimInfo::MoleculeIterator mi;
1219 Molecule::AtomIterator ai;
1220 Atom* atom;
1221 Molecule::RigidBodyIterator rbIter;
1222 RigidBody* rb;
1223 Molecule::CutoffGroupIterator ci;
1224 CutoffGroup* cg;
1225
1226 // forces and potentials are zeroed here, before any are
1227 // accumulated.
1228
1229 Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot();
1230
1231 snap->setBondPotential(0.0);
1232 snap->setBendPotential(0.0);
1233 snap->setTorsionPotential(0.0);
1234 snap->setInversionPotential(0.0);
1235
1236 potVec zeroPot(0.0);
1237 snap->setLongRangePotentials(zeroPot);
1238 snap->setExcludedPotentials(zeroPot);
1239 if (doPotentialSelection_) snap->setSelectionPotentials(zeroPot);
1240
1241 snap->setRestraintPotential(0.0);
1242 snap->setRawPotential(0.0);
1243
1244 // First we zero out for mol1
1245 for (atom = mol1->beginAtom(ai); atom != NULL; atom = mol1->nextAtom(ai)) {
1246 atom->zeroForcesAndTorques();
1247 }
1248 // change the positions of atoms which belong to the rigidbodies
1249 for (rb = mol1->beginRigidBody(rbIter); rb != NULL;
1250 rb = mol1->nextRigidBody(rbIter)) {
1251 rb->zeroForcesAndTorques();
1252 }
1253 if (info_->getNGlobalCutoffGroups() != info_->getNGlobalAtoms()) {
1254 for (cg = mol1->beginCutoffGroup(ci); cg != NULL;
1255 cg = mol1->nextCutoffGroup(ci)) {
1256 // calculate the center of mass of cutoff group
1257 cg->updateCOM();
1258 }
1259 }
1260
1261 // Next we zero out for mol2
1262 for (atom = mol2->beginAtom(ai); atom != NULL; atom = mol2->nextAtom(ai)) {
1263 atom->zeroForcesAndTorques();
1264 }
1265 // change the positions of atoms which belong to the rigidbodies
1266 for (rb = mol2->beginRigidBody(rbIter); rb != NULL;
1267 rb = mol2->nextRigidBody(rbIter)) {
1268 rb->zeroForcesAndTorques();
1269 }
1270 if (info_->getNGlobalCutoffGroups() != info_->getNGlobalAtoms()) {
1271 for (cg = mol2->beginCutoffGroup(ci); cg != NULL;
1272 cg = mol2->nextCutoffGroup(ci)) {
1273 // calculate the center of mass of cutoff group
1274 cg->updateCOM();
1275 }
1276 }
1277
1278 // Zero out the virial tensor
1279 virialTensor *= 0.0;
1280 // Zero out the heatFlux
1281 fDecomp_->setHeatFlux(Vector3d(0.0));
1282 }
1283
1284 void ForceManager::selectedShortRangeInteractions(Molecule* mol1,
1285 Molecule* mol2) {
1286 RigidBody* rb;
1287 Bond* bond;
1288 Bend* bend;
1289 Torsion* torsion;
1290 Inversion* inversion;
1291 SimInfo::MoleculeIterator mi;
1292 Molecule::RigidBodyIterator rbIter;
1293 Molecule::BondIterator bondIter;
1294 Molecule::BendIterator bendIter;
1295 Molecule::TorsionIterator torsionIter;
1296 Molecule::InversionIterator inversionIter;
1297 RealType bondPotential = 0.0;
1298 RealType bendPotential = 0.0;
1299 RealType torsionPotential = 0.0;
1300 RealType inversionPotential = 0.0;
1301 potVec selectionPotential(0.0);
1302
1303 // First compute for mol1
1304 for (rb = mol1->beginRigidBody(rbIter); rb != NULL;
1305 rb = mol1->nextRigidBody(rbIter)) {
1306 rb->updateAtoms();
1307 }
1308
1309 for (bond = mol1->beginBond(bondIter); bond != NULL;
1310 bond = mol1->nextBond(bondIter)) {
1311 bond->calcForce(doParticlePot_);
1312 bondPotential += bond->getPotential();
1313 if (doPotentialSelection_) {
1314 if (seleMan_.isSelected(bond->getAtomA()) ||
1315 seleMan_.isSelected(bond->getAtomB())) {
1316 selectionPotential[BONDED_FAMILY] += bond->getPotential();
1317 }
1318 }
1319 }
1320
1321 for (bend = mol1->beginBend(bendIter); bend != NULL;
1322 bend = mol1->nextBend(bendIter)) {
1323 RealType angle;
1324 bend->calcForce(angle, doParticlePot_);
1325 RealType currBendPot = bend->getPotential();
1326
1327 bendPotential += bend->getPotential();
1328 map<Bend*, BendDataSet>::iterator i = bendDataSets.find(bend);
1329 if (i == bendDataSets.end()) {
1330 BendDataSet dataSet;
1331 dataSet.prev.angle = dataSet.curr.angle = angle;
1332 dataSet.prev.potential = dataSet.curr.potential = currBendPot;
1333 dataSet.deltaV = 0.0;
1334 bendDataSets.insert(map<Bend*, BendDataSet>::value_type(bend, dataSet));
1335 } else {
1336 i->second.prev.angle = i->second.curr.angle;
1337 i->second.prev.potential = i->second.curr.potential;
1338 i->second.curr.angle = angle;
1339 i->second.curr.potential = currBendPot;
1340 i->second.deltaV =
1341 fabs(i->second.curr.potential - i->second.prev.potential);
1342 }
1343 if (doPotentialSelection_) {
1344 if (seleMan_.isSelected(bend->getAtomA()) ||
1345 seleMan_.isSelected(bend->getAtomB()) ||
1346 seleMan_.isSelected(bend->getAtomC())) {
1347 selectionPotential[BONDED_FAMILY] += bend->getPotential();
1348 }
1349 }
1350 }
1351
1352 for (torsion = mol1->beginTorsion(torsionIter); torsion != NULL;
1353 torsion = mol1->nextTorsion(torsionIter)) {
1354 RealType angle;
1355 torsion->calcForce(angle, doParticlePot_);
1356 RealType currTorsionPot = torsion->getPotential();
1357 torsionPotential += torsion->getPotential();
1358 map<Torsion*, TorsionDataSet>::iterator i = torsionDataSets.find(torsion);
1359 if (i == torsionDataSets.end()) {
1360 TorsionDataSet dataSet;
1361 dataSet.prev.angle = dataSet.curr.angle = angle;
1362 dataSet.prev.potential = dataSet.curr.potential = currTorsionPot;
1363 dataSet.deltaV = 0.0;
1364 torsionDataSets.insert(
1365 map<Torsion*, TorsionDataSet>::value_type(torsion, dataSet));
1366 } else {
1367 i->second.prev.angle = i->second.curr.angle;
1368 i->second.prev.potential = i->second.curr.potential;
1369 i->second.curr.angle = angle;
1370 i->second.curr.potential = currTorsionPot;
1371 i->second.deltaV =
1372 fabs(i->second.curr.potential - i->second.prev.potential);
1373 }
1374 if (doPotentialSelection_) {
1375 if (seleMan_.isSelected(torsion->getAtomA()) ||
1376 seleMan_.isSelected(torsion->getAtomB()) ||
1377 seleMan_.isSelected(torsion->getAtomC()) ||
1378 seleMan_.isSelected(torsion->getAtomD())) {
1379 selectionPotential[BONDED_FAMILY] += torsion->getPotential();
1380 }
1381 }
1382 }
1383
1384 for (inversion = mol1->beginInversion(inversionIter); inversion != NULL;
1385 inversion = mol1->nextInversion(inversionIter)) {
1386 RealType angle;
1387 inversion->calcForce(angle, doParticlePot_);
1388 RealType currInversionPot = inversion->getPotential();
1389 inversionPotential += inversion->getPotential();
1390 map<Inversion*, InversionDataSet>::iterator i =
1391 inversionDataSets.find(inversion);
1392 if (i == inversionDataSets.end()) {
1393 InversionDataSet dataSet;
1394 dataSet.prev.angle = dataSet.curr.angle = angle;
1395 dataSet.prev.potential = dataSet.curr.potential = currInversionPot;
1396 dataSet.deltaV = 0.0;
1397 inversionDataSets.insert(
1398 map<Inversion*, InversionDataSet>::value_type(inversion, dataSet));
1399 } else {
1400 i->second.prev.angle = i->second.curr.angle;
1401 i->second.prev.potential = i->second.curr.potential;
1402 i->second.curr.angle = angle;
1403 i->second.curr.potential = currInversionPot;
1404 i->second.deltaV =
1405 fabs(i->second.curr.potential - i->second.prev.potential);
1406 }
1407 if (doPotentialSelection_) {
1408 if (seleMan_.isSelected(inversion->getAtomA()) ||
1409 seleMan_.isSelected(inversion->getAtomB()) ||
1410 seleMan_.isSelected(inversion->getAtomC()) ||
1411 seleMan_.isSelected(inversion->getAtomD())) {
1412 selectionPotential[BONDED_FAMILY] += inversion->getPotential();
1413 }
1414 }
1415 }
1416
1417 // Next compute for mol2
1418 for (rb = mol2->beginRigidBody(rbIter); rb != NULL;
1419 rb = mol2->nextRigidBody(rbIter)) {
1420 rb->updateAtoms();
1421 }
1422
1423 for (bond = mol2->beginBond(bondIter); bond != NULL;
1424 bond = mol2->nextBond(bondIter)) {
1425 bond->calcForce(doParticlePot_);
1426 bondPotential += bond->getPotential();
1427 if (doPotentialSelection_) {
1428 if (seleMan_.isSelected(bond->getAtomA()) ||
1429 seleMan_.isSelected(bond->getAtomB())) {
1430 selectionPotential[BONDED_FAMILY] += bond->getPotential();
1431 }
1432 }
1433 }
1434
1435 for (bend = mol2->beginBend(bendIter); bend != NULL;
1436 bend = mol2->nextBend(bendIter)) {
1437 RealType angle;
1438 bend->calcForce(angle, doParticlePot_);
1439 RealType currBendPot = bend->getPotential();
1440
1441 bendPotential += bend->getPotential();
1442 map<Bend*, BendDataSet>::iterator i = bendDataSets.find(bend);
1443 if (i == bendDataSets.end()) {
1444 BendDataSet dataSet;
1445 dataSet.prev.angle = dataSet.curr.angle = angle;
1446 dataSet.prev.potential = dataSet.curr.potential = currBendPot;
1447 dataSet.deltaV = 0.0;
1448 bendDataSets.insert(map<Bend*, BendDataSet>::value_type(bend, dataSet));
1449 } else {
1450 i->second.prev.angle = i->second.curr.angle;
1451 i->second.prev.potential = i->second.curr.potential;
1452 i->second.curr.angle = angle;
1453 i->second.curr.potential = currBendPot;
1454 i->second.deltaV =
1455 fabs(i->second.curr.potential - i->second.prev.potential);
1456 }
1457 if (doPotentialSelection_) {
1458 if (seleMan_.isSelected(bend->getAtomA()) ||
1459 seleMan_.isSelected(bend->getAtomB()) ||
1460 seleMan_.isSelected(bend->getAtomC())) {
1461 selectionPotential[BONDED_FAMILY] += bend->getPotential();
1462 }
1463 }
1464 }
1465
1466 for (torsion = mol2->beginTorsion(torsionIter); torsion != NULL;
1467 torsion = mol2->nextTorsion(torsionIter)) {
1468 RealType angle;
1469 torsion->calcForce(angle, doParticlePot_);
1470 RealType currTorsionPot = torsion->getPotential();
1471 torsionPotential += torsion->getPotential();
1472 map<Torsion*, TorsionDataSet>::iterator i = torsionDataSets.find(torsion);
1473 if (i == torsionDataSets.end()) {
1474 TorsionDataSet dataSet;
1475 dataSet.prev.angle = dataSet.curr.angle = angle;
1476 dataSet.prev.potential = dataSet.curr.potential = currTorsionPot;
1477 dataSet.deltaV = 0.0;
1478 torsionDataSets.insert(
1479 map<Torsion*, TorsionDataSet>::value_type(torsion, dataSet));
1480 } else {
1481 i->second.prev.angle = i->second.curr.angle;
1482 i->second.prev.potential = i->second.curr.potential;
1483 i->second.curr.angle = angle;
1484 i->second.curr.potential = currTorsionPot;
1485 i->second.deltaV =
1486 fabs(i->second.curr.potential - i->second.prev.potential);
1487 }
1488 if (doPotentialSelection_) {
1489 if (seleMan_.isSelected(torsion->getAtomA()) ||
1490 seleMan_.isSelected(torsion->getAtomB()) ||
1491 seleMan_.isSelected(torsion->getAtomC()) ||
1492 seleMan_.isSelected(torsion->getAtomD())) {
1493 selectionPotential[BONDED_FAMILY] += torsion->getPotential();
1494 }
1495 }
1496 }
1497
1498 for (inversion = mol2->beginInversion(inversionIter); inversion != NULL;
1499 inversion = mol2->nextInversion(inversionIter)) {
1500 RealType angle;
1501 inversion->calcForce(angle, doParticlePot_);
1502 RealType currInversionPot = inversion->getPotential();
1503 inversionPotential += inversion->getPotential();
1504 map<Inversion*, InversionDataSet>::iterator i =
1505 inversionDataSets.find(inversion);
1506 if (i == inversionDataSets.end()) {
1507 InversionDataSet dataSet;
1508 dataSet.prev.angle = dataSet.curr.angle = angle;
1509 dataSet.prev.potential = dataSet.curr.potential = currInversionPot;
1510 dataSet.deltaV = 0.0;
1511 inversionDataSets.insert(
1512 map<Inversion*, InversionDataSet>::value_type(inversion, dataSet));
1513 } else {
1514 i->second.prev.angle = i->second.curr.angle;
1515 i->second.prev.potential = i->second.curr.potential;
1516 i->second.curr.angle = angle;
1517 i->second.curr.potential = currInversionPot;
1518 i->second.deltaV =
1519 fabs(i->second.curr.potential - i->second.prev.potential);
1520 }
1521 if (doPotentialSelection_) {
1522 if (seleMan_.isSelected(inversion->getAtomA()) ||
1523 seleMan_.isSelected(inversion->getAtomB()) ||
1524 seleMan_.isSelected(inversion->getAtomC()) ||
1525 seleMan_.isSelected(inversion->getAtomD())) {
1526 selectionPotential[BONDED_FAMILY] += inversion->getPotential();
1527 }
1528 }
1529 }
1530
1531#ifdef IS_MPI
1532 // Collect from all nodes. This should eventually be moved into a
1533 // SystemDecomposition, but this is a better place than in
1534 // Thermo to do the collection.
1535
1536 MPI_Allreduce(MPI_IN_PLACE, &bondPotential, 1, MPI_REALTYPE, MPI_SUM,
1537 MPI_COMM_WORLD);
1538 MPI_Allreduce(MPI_IN_PLACE, &bendPotential, 1, MPI_REALTYPE, MPI_SUM,
1539 MPI_COMM_WORLD);
1540 MPI_Allreduce(MPI_IN_PLACE, &torsionPotential, 1, MPI_REALTYPE, MPI_SUM,
1541 MPI_COMM_WORLD);
1542 MPI_Allreduce(MPI_IN_PLACE, &inversionPotential, 1, MPI_REALTYPE, MPI_SUM,
1543 MPI_COMM_WORLD);
1544 MPI_Allreduce(MPI_IN_PLACE, &selectionPotential[BONDED_FAMILY], 1,
1545 MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1546#endif
1547
1548 Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
1549
1550 curSnapshot->setBondPotential(bondPotential);
1551 curSnapshot->setBendPotential(bendPotential);
1552 curSnapshot->setTorsionPotential(torsionPotential);
1553 curSnapshot->setInversionPotential(inversionPotential);
1554 curSnapshot->setSelectionPotentials(selectionPotential);
1555
1556 // RealType shortRangePotential = bondPotential + bendPotential +
1557 // torsionPotential + inversionPotential;
1558
1559 // curSnapshot->setShortRangePotential(shortRangePotential);
1560 }
1561
1562 void ForceManager::selectedLongRangeInteractions(Molecule* mol1,
1563 Molecule* mol2) {
1564 Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
1565 DataStorage* config = &(curSnapshot->atomData);
1566 DataStorage* cgConfig = &(curSnapshot->cgData);
1567
1568 // calculate the center of mass of cutoff group
1569
1570 SimInfo::MoleculeIterator mi;
1571 Molecule::CutoffGroupIterator ci;
1572 CutoffGroup* cg;
1573
1574 if (info_->getNCutoffGroups() != info_->getNAtoms()) {
1575 for (cg = mol1->beginCutoffGroup(ci); cg != NULL;
1576 cg = mol1->nextCutoffGroup(ci)) {
1577 cg->updateCOM();
1578 }
1579 for (cg = mol2->beginCutoffGroup(ci); cg != NULL;
1580 cg = mol2->nextCutoffGroup(ci)) {
1581 cg->updateCOM();
1582 }
1583 } else {
1584 // center of mass of the group is the same as position of the atom
1585 // if cutoff group does not exist
1586 cgConfig->position = config->position;
1587 cgConfig->velocity = config->velocity;
1588 }
1589
1590 fDecomp_->zeroWorkArrays();
1591 fDecomp_->distributeData();
1592
1593 int cg1, cg2, atom1, atom2;
1594 Vector3d d_grp, dag, gvel2, vel2;
1595 RealType rgrpsq, rgrp;
1596 RealType vij(0.0);
1597 Vector3d fij, fg;
1598 bool in_switching_region;
1599 RealType dswdr, swderiv;
1600 vector<int> atomListColumn, atomListRow;
1601
1602 potVec longRangePotential(0.0);
1603 potVec selfPotential(0.0);
1604 potVec selectionPotential(0.0);
1605 RealType reciprocalPotential(0.0);
1606 RealType surfacePotential(0.0);
1607
1608 RealType mf;
1609 bool newAtom1;
1610 int gid1, gid2;
1611
1612 vector<int>::iterator ia, jb;
1613
1614 int loopStart, loopEnd;
1615
1616 idat.rcut = rCut_;
1617 idat.shiftedPot = (cutoffMethod_ == SHIFTED_POTENTIAL) ? true : false;
1618 idat.shiftedForce =
1619 (cutoffMethod_ == SHIFTED_FORCE || cutoffMethod_ == TAYLOR_SHIFTED) ?
1620 true :
1621 false;
1622 idat.doParticlePot = doParticlePot_;
1623 idat.doElectricField = doElectricField_;
1624 idat.doSitePotential = doSitePotential_;
1625 sdat.doParticlePot = doParticlePot_;
1626
1627 loopEnd = PAIR_LOOP;
1628 if (info_->requiresPrepair()) {
1629 loopStart = PREPAIR_LOOP;
1630 } else {
1631 loopStart = PAIR_LOOP;
1632 }
1633 for (int iLoop = loopStart; iLoop <= loopEnd; iLoop++) {
1634 if (iLoop == loopStart) {
1635 bool update_nlist = fDecomp_->checkNeighborList(savedPositions_);
1636
1637 if (update_nlist) {
1638 if (!usePeriodicBoundaryConditions_)
1639 Mat3x3d bbox = thermo->getBoundingBox();
1640 fDecomp_->buildNeighborList(neighborList_, point_, savedPositions_);
1641 }
1642 }
1643
1644 for (cg1 = 0; cg1 < int(point_.size()) - 1; cg1++) {
1645 atomListRow = fDecomp_->getAtomsInGroupRow(cg1);
1646 newAtom1 = true;
1647
1648 for (int m2 = point_[cg1]; m2 < point_[cg1 + 1]; m2++) {
1649 cg2 = neighborList_[m2];
1650
1651 d_grp = fDecomp_->getIntergroupVector(cg1, cg2);
1652
1653 // already wrapped in the getIntergroupVector call:
1654 // curSnapshot->wrapVector(d_grp);
1655 rgrpsq = d_grp.lengthSquare();
1656
1657 if (rgrpsq < rCutSq_) {
1658 if (iLoop == PAIR_LOOP) {
1659 vij = 0.0;
1660 fij.zero();
1661 idat.eField1.zero();
1662 idat.eField2.zero();
1663 idat.sPot1 = 0.0;
1664 idat.sPot2 = 0.0;
1665 }
1666
1667 in_switching_region =
1668 switcher_->getSwitch(rgrpsq, idat.sw, dswdr, rgrp);
1669
1670 atomListColumn = fDecomp_->getAtomsInGroupColumn(cg2);
1671
1672 if (doHeatFlux_) gvel2 = fDecomp_->getGroupVelocityColumn(cg2);
1673
1674 for (ia = atomListRow.begin(); ia != atomListRow.end(); ++ia) {
1675 if (atomListRow.size() > 1) newAtom1 = true;
1676 atom1 = (*ia);
1677
1678 if (doPotentialSelection_) {
1679 gid1 = fDecomp_->getGlobalIDRow(atom1);
1680 idat.isSelected = seleMan_.isGlobalIDSelected(gid1);
1681 }
1682
1683 for (jb = atomListColumn.begin(); jb != atomListColumn.end();
1684 ++jb) {
1685 atom2 = (*jb);
1686
1687 if (doPotentialSelection_) {
1688 gid2 = fDecomp_->getGlobalIDCol(atom2);
1689 idat.isSelected |= seleMan_.isGlobalIDSelected(gid2);
1690 }
1691
1692 if (!fDecomp_->skipAtomPair(atom1, atom2, cg1, cg2)) {
1693 idat.vpair = 0.0;
1694 idat.pot = 0.0;
1695 idat.excludedPot = 0.0;
1696 idat.selePot = 0.0;
1697 idat.f1.zero();
1698 idat.dVdFQ1 = 0.0;
1699 idat.dVdFQ2 = 0.0;
1700
1701 fDecomp_->fillInteractionData(idat, atom1, atom2, newAtom1);
1702 newAtom1 = false;
1703
1704 idat.topoDist =
1705 fDecomp_->getTopologicalDistance(atom1, atom2);
1706 idat.vdwMult = vdwScale_[idat.topoDist];
1707 idat.electroMult = electrostaticScale_[idat.topoDist];
1708
1709 if (atomListRow.size() == 1 && atomListColumn.size() == 1) {
1710 idat.d = d_grp;
1711 idat.r2 = rgrpsq;
1712 if (doHeatFlux_) vel2 = gvel2;
1713 } else {
1714 idat.d = fDecomp_->getInteratomicVector(atom1, atom2);
1715 curSnapshot->wrapVector(idat.d);
1716 idat.r2 = idat.d.lengthSquare();
1717 if (doHeatFlux_)
1718 vel2 = fDecomp_->getAtomVelocityColumn(atom2);
1719 }
1720
1721 idat.rij = sqrt(idat.r2);
1722
1723 if (iLoop == PREPAIR_LOOP) {
1724 interactionMan_->doPrePair(idat);
1725 fDecomp_->unpackPrePairData(idat, atom1, atom2);
1726 } else {
1727 interactionMan_->doPair(idat);
1728 fDecomp_->unpackInteractionData(idat, atom1, atom2);
1729 vij += idat.vpair;
1730 fij += idat.f1;
1731 virialTensor -= outProduct(idat.d, idat.f1);
1732 if (doHeatFlux_)
1733 fDecomp_->addToHeatFlux(idat.d * dot(idat.f1, vel2));
1734 }
1735 }
1736 }
1737 }
1738
1739 if (iLoop == PAIR_LOOP) {
1740 if (in_switching_region) {
1741 swderiv = vij * dswdr / rgrp;
1742 fg = swderiv * d_grp;
1743 fij += fg;
1744
1745 if (atomListRow.size() == 1 && atomListColumn.size() == 1) {
1746 if (!fDecomp_->skipAtomPair(atomListRow[0], atomListColumn[0],
1747 cg1, cg2)) {
1748 virialTensor -= outProduct(idat.d, fg);
1749 if (doHeatFlux_)
1750 fDecomp_->addToHeatFlux(idat.d * dot(fg, vel2));
1751 }
1752 }
1753
1754 for (ia = atomListRow.begin(); ia != atomListRow.end(); ++ia) {
1755 atom1 = (*ia);
1756 mf = fDecomp_->getMassFactorRow(atom1);
1757 // fg is the force on atom ia due to cutoff group's
1758 // presence in switching region
1759 fg = swderiv * d_grp * mf;
1760 fDecomp_->addForceToAtomRow(atom1, fg);
1761 if (atomListRow.size() > 1) {
1762 if (info_->usesAtomicVirial()) {
1763 // find the distance between the atom
1764 // and the center of the cutoff group:
1765 dag = fDecomp_->getAtomToGroupVectorRow(atom1, cg1);
1766 virialTensor -= outProduct(dag, fg);
1767 if (doHeatFlux_)
1768 fDecomp_->addToHeatFlux(dag * dot(fg, vel2));
1769 }
1770 }
1771 }
1772 for (jb = atomListColumn.begin(); jb != atomListColumn.end();
1773 ++jb) {
1774 atom2 = (*jb);
1775 mf = fDecomp_->getMassFactorColumn(atom2);
1776 // fg is the force on atom jb due to cutoff group's
1777 // presence in switching region
1778 fg = -swderiv * d_grp * mf;
1779 fDecomp_->addForceToAtomColumn(atom2, fg);
1780
1781 if (atomListColumn.size() > 1) {
1782 if (info_->usesAtomicVirial()) {
1783 // find the distance between the atom
1784 // and the center of the cutoff group:
1785 dag = fDecomp_->getAtomToGroupVectorColumn(atom2, cg2);
1786 virialTensor -= outProduct(dag, fg);
1787 if (doHeatFlux_)
1788 fDecomp_->addToHeatFlux(dag * dot(fg, vel2));
1789 }
1790 }
1791 }
1792 }
1793 // if (!info_->usesAtomicVirial()) {
1794 // virialTensor -= outProduct(d_grp, fij);
1795 // if (doHeatFlux_)
1796 // fDecomp_->addToHeatFlux( d_grp * dot(fij, vel2));
1797 //}
1798 }
1799 }
1800 }
1801 }
1802
1803 if (iLoop == PREPAIR_LOOP) {
1804 if (info_->requiresPrepair()) {
1805 fDecomp_->collectIntermediateData();
1806
1807 for (unsigned int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) {
1808 if (doPotentialSelection_) {
1809 gid1 = fDecomp_->getGlobalID(atom1);
1810 sdat.isSelected = seleMan_.isGlobalIDSelected(gid1);
1811 }
1812
1813 fDecomp_->fillPreForceData(sdat, atom1);
1814 interactionMan_->doPreForce(sdat);
1815 fDecomp_->unpackPreForceData(sdat, atom1);
1816 }
1817
1818 fDecomp_->distributeIntermediateData();
1819 }
1820 }
1821 }
1822
1823 // collects pairwise information
1824 fDecomp_->collectData();
1825 if (cutoffMethod_ == EWALD_FULL) {
1826 interactionMan_->doReciprocalSpaceSum(reciprocalPotential);
1827 curSnapshot->setReciprocalPotential(reciprocalPotential);
1828 }
1829
1830 if (useSurfaceTerm_) {
1831 interactionMan_->doSurfaceTerm(useSlabGeometry_, axis_, surfacePotential);
1832 curSnapshot->setSurfacePotential(surfacePotential);
1833 }
1834
1835 if (info_->requiresSelfCorrection()) {
1836 for (unsigned int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) {
1837 if (doPotentialSelection_) {
1838 gid1 = fDecomp_->getGlobalID(atom1);
1839 sdat.isSelected = seleMan_.isGlobalIDSelected(gid1);
1840 }
1841
1842 fDecomp_->fillSelfData(sdat, atom1);
1843 interactionMan_->doSelfCorrection(sdat);
1844 fDecomp_->unpackSelfData(sdat, atom1);
1845 }
1846 }
1847
1848 // collects single-atom information
1849 fDecomp_->collectSelfData();
1850
1851 longRangePotential = fDecomp_->getPairwisePotential();
1852 curSnapshot->setLongRangePotentials(longRangePotential);
1853
1854 selfPotential = fDecomp_->getSelfPotential();
1855 curSnapshot->setSelfPotentials(selfPotential);
1856
1857 curSnapshot->setExcludedPotentials(fDecomp_->getExcludedSelfPotential() +
1858 fDecomp_->getExcludedPotential());
1859
1860 if (doPotentialSelection_) {
1861 selectionPotential = curSnapshot->getSelectionPotentials();
1862 selectionPotential += fDecomp_->getSelectedSelfPotential();
1863 selectionPotential += fDecomp_->getSelectedPotential();
1864 curSnapshot->setSelectionPotentials(selectionPotential);
1865 }
1866 }
1867
1868 void ForceManager::selectedPostCalculation(Molecule* mol1, Molecule* mol2) {
1869 for (auto& forceModifier : forceModifiers_)
1870 forceModifier->modifyForces();
1871
1872 // Modify the rigid bodies in response to the applied force
1873 // modifications
1874 SimInfo::MoleculeIterator mi;
1875 Molecule::RigidBodyIterator rbIter;
1876 RigidBody* rb;
1877 Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
1878
1879 // Collect the atomic forces onto rigid bodies
1880 for (rb = mol1->beginRigidBody(rbIter); rb != NULL;
1881 rb = mol1->nextRigidBody(rbIter)) {
1882 Mat3x3d rbTau = rb->calcForcesAndTorquesAndVirial();
1883 virialTensor += rbTau;
1884 }
1885
1886 for (rb = mol2->beginRigidBody(rbIter); rb != NULL;
1887 rb = mol2->nextRigidBody(rbIter)) {
1888 Mat3x3d rbTau = rb->calcForcesAndTorquesAndVirial();
1889 virialTensor += rbTau;
1890 }
1891
1892#ifdef IS_MPI
1893 MPI_Allreduce(MPI_IN_PLACE, virialTensor.getArrayPointer(), 9, MPI_REALTYPE,
1894 MPI_SUM, MPI_COMM_WORLD);
1895#endif
1896 curSnapshot->setVirialTensor(virialTensor);
1897
1898 /*
1899 if (info_->getSimParams()->getUseLongRangeCorrections()) {
1900
1901 RealType vol = curSnapshot->getVolume();
1902 RealType Elrc(0.0);
1903 RealType Wlrc(0.0);
1904
1905 AtomTypeSet::iterator i;
1906 AtomTypeSet::iterator j;
1907
1908 RealType n_i, n_j;
1909 RealType rho_i, rho_j;
1910 pair<RealType, RealType> LRI;
1911
1912 for (i = atomTypes_.begin(); i != atomTypes_.end(); ++i) {
1913 n_i = RealType(info_->getGlobalCountOfType(*i));
1914 rho_i = n_i / vol;
1915 for (j = atomTypes_.begin(); j != atomTypes_.end(); ++j) {
1916 n_j = RealType(info_->getGlobalCountOfType(*j));
1917 rho_j = n_j / vol;
1918
1919 LRI = interactionMan_->getLongRangeIntegrals( (*i), (*j) );
1920
1921 Elrc += n_i * rho_j * LRI.first;
1922 Wlrc -= rho_i * rho_j * LRI.second;
1923 }
1924 }
1925 Elrc *= 2.0 * Constants::PI;
1926 Wlrc *= 2.0 * Constants::PI;
1927
1928 RealType lrp = curSnapshot->getLongRangePotential();
1929 curSnapshot->setLongRangePotential(lrp + Elrc);
1930 virialTensor += Wlrc * SquareMatrix3<RealType>::identity();
1931 curSnapshot->setVirialTensor(virialTensor);
1932 }
1933 */
1934 }
1935} // namespace OpenMD
LHDForceModifier.hpp
Langevin force modifier with intramolecular RPY hydrodynamic interactions for flexible bead molecules...
Light.hpp
Rotating Electric Field perturbation.
MagneticField.hpp
Uniform Magnetic Field perturbation.
UniformField.hpp
Uniform Electric Field perturbation.
UniformGradient.hpp
Uniform Electric Field Gradient perturbation.
OpenMD::ForceManager::setupCutoffs
virtual void setupCutoffs()
setupCutoffs
Definition ForceManager.cpp:132
OpenMD::ForceManager::sft_
SwitchingFunctionType sft_
Type of switching function in use.
Definition ForceManager.hpp:122
OpenMD::ForceManager::rCut_
RealType rCut_
cutoff radius for non-bonded interactions
Definition ForceManager.hpp:123
OpenMD::ForceManager::cutoffMethod_
CutoffMethod cutoffMethod_
Cutoff Method for most non-bonded interactions.
Definition ForceManager.hpp:127
OpenMD::ForceManager::rSwitch_
RealType rSwitch_
inner radius of switching function
Definition ForceManager.hpp:125
OpenMD::LDForceModifier
Force modifier for Langevin Dynamics applying friction and random forces as well as torques.
Definition LDForceModifier.hpp:81
OpenMD::LangevinHullForceModifier
Force modifier for NPT Langevin Hull Dynamics applying friction and random forces as well as torques.
Definition LangevinHullForceModifier.hpp:74
OpenMD::MagneticField
Applies a uniform (vector) magnetic field to the system.
Definition MagneticField.hpp:84
OpenMD::Perturbations::Light
Applies an EM field representing light to the system.
Definition Light.hpp:107
OpenMD::SimInfo
One of the heavy-weight classes of OpenMD, SimInfo maintains objects and variables relating to the cu...
Definition SimInfo.hpp:96
OpenMD::UniformField
Applies a uniform (vector) electric field to the system.
Definition UniformField.hpp:86
OpenMD::UniformGradient
Applies a uniform electric field gradient to the system.
Definition UniformGradient.hpp:125
OpenMD
This basic Periodic Table class was originally taken from the data.cpp file in OpenBabel.
Definition ActionCorrFunc.cpp:63
OpenMD::dot
Real dot(const DynamicVector< Real > &v1, const DynamicVector< Real > &v2)
Returns the dot product of two DynamicVectors.
Definition DynamicVector.hpp:452
OpenMD::BONDED_FAMILY
@ BONDED_FAMILY
directly bonded 1-2, 1-3, or 1-4 interactions
Definition NonBondedInteraction.hpp:75