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