OpenMD 3.1
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
RNEMD.cpp
1/*
2 * Copyright (c) 2004-present, The University of Notre Dame. All rights
3 * reserved.
4 *
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions are met:
7 *
8 * 1. Redistributions of source code must retain the above copyright notice,
9 * this list of conditions and the following disclaimer.
10 *
11 * 2. Redistributions in binary form must reproduce the above copyright notice,
12 * this list of conditions and the following disclaimer in the documentation
13 * and/or other materials provided with the distribution.
14 *
15 * 3. Neither the name of the copyright holder nor the names of its
16 * contributors may be used to endorse or promote products derived from
17 * this software without specific prior written permission.
18 *
19 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
20 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
21 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
22 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
23 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
24 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
25 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
26 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
27 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
28 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
29 * POSSIBILITY OF SUCH DAMAGE.
30 *
31 * SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your
32 * research, please cite the appropriate papers when you publish your
33 * work. Good starting points are:
34 *
35 * [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
36 * [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
37 * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).
38 * [4] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
39 * [5] Kuang & Gezelter, Mol. Phys., 110, 691-701 (2012).
40 * [6] Lamichhane, Gezelter & Newman, J. Chem. Phys. 141, 134109 (2014).
41 * [7] Lamichhane, Newman & Gezelter, J. Chem. Phys. 141, 134110 (2014).
42 * [8] Bhattarai, Newman & Gezelter, Phys. Rev. B 99, 094106 (2019).
43 */
44
45#include "rnemd/RNEMD.hpp"
46
47#include <algorithm>
48#include <cmath>
49#include <map>
50#include <memory>
51#include <set>
52#include <sstream>
53#include <string>
54#include <typeinfo>
55#include <utility>
56#include <vector>
57
58#ifdef IS_MPI
59#include <mpi.h>
60#endif
61
62#include "brains/ForceManager.hpp"
63#include "brains/Thermo.hpp"
64#include "io/Globals.hpp"
65#include "math/ConvexHull.hpp"
66#include "math/Polynomial.hpp"
67#include "math/SquareMatrix3.hpp"
68#include "math/Vector.hpp"
69#include "math/Vector3.hpp"
70#include "primitives/Molecule.hpp"
71#include "primitives/StuntDouble.hpp"
72#include "rnemd/RNEMDParameters.hpp"
73#include "types/FixedChargeAdapter.hpp"
74#include "types/FluctuatingChargeAdapter.hpp"
75#include "utils/Accumulator.hpp"
76#include "utils/AccumulatorView.hpp"
77#include "utils/Constants.hpp"
78
79using namespace OpenMD::Utils;
80
81namespace OpenMD::RNEMD {
82
83 RNEMD::RNEMD(SimInfo* info, ForceManager*) :
84 info_(info), commonA_(info_), commonB_(info_), evaluator_(info_),
85 seleMan_(info_), evaluatorA_(info_), evaluatorB_(info_), seleManA_(info_),
86 seleManB_(info_), outputEvaluator_(info_), outputSeleMan_(info_) {
87 trialCount_ = 0;
88 failTrialCount_ = 0;
89 failRootCount_ = 0;
90
91 Globals* simParams = info->getSimParams();
92 RNEMDParameters* rnemdParams = simParams->getRNEMDParameters();
93
94 usePeriodicBoundaryConditions_ =
95 simParams->getUsePeriodicBoundaryConditions();
96
97 doRNEMD_ = rnemdParams->getUseRNEMD();
98 if (!doRNEMD_) return;
99
100 // Determine Flux Type
101 std::map<std::string, RNEMDFluxType> stringToFluxType;
102
103 stringToFluxType["KE"] = rnemdKE;
104 stringToFluxType["Px"] = rnemdPx;
105 stringToFluxType["Py"] = rnemdPy;
106 stringToFluxType["Pz"] = rnemdPz;
107 stringToFluxType["Pvector"] = rnemdPvector;
108 stringToFluxType["Lx"] = rnemdLx;
109 stringToFluxType["Ly"] = rnemdLy;
110 stringToFluxType["Lz"] = rnemdLz;
111 stringToFluxType["Lvector"] = rnemdLvector;
112 stringToFluxType["Particle"] = rnemdParticle;
113 stringToFluxType["Particle+KE"] = rnemdParticleKE;
114 stringToFluxType["KE+Px"] = rnemdKePx;
115 stringToFluxType["KE+Py"] = rnemdKePy;
116 stringToFluxType["KE+Pvector"] = rnemdKePvector;
117 stringToFluxType["KE+Lx"] = rnemdKeLx;
118 stringToFluxType["KE+Ly"] = rnemdKeLy;
119 stringToFluxType["KE+Lz"] = rnemdKeLz;
120 stringToFluxType["KE+Lvector"] = rnemdKeLvector;
121
122 if (rnemdParams->haveFluxType()) {
123 rnemdFluxTypeLabel_ = rnemdParams->getFluxType();
124 rnemdFluxType_ = stringToFluxType.find(rnemdFluxTypeLabel_)->second;
125 } else {
126 std::string allowedFluxTypes;
127 int currentLineLength = 0;
128
129 for (std::map<std::string, RNEMDFluxType>::iterator fluxStrIter =
130 stringToFluxType.begin();
131 fluxStrIter != stringToFluxType.end(); ++fluxStrIter) {
132 allowedFluxTypes += fluxStrIter->first + ", ";
133 currentLineLength += fluxStrIter->first.length() + 2;
134
135 if (currentLineLength >= 50) {
136 allowedFluxTypes += "\n\t\t";
137 currentLineLength = 0;
138 }
139 }
140
141 allowedFluxTypes.erase(allowedFluxTypes.length() - 2, 2);
142
143 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
144 "RNEMD: No fluxType was set in the omd file. This parameter\n"
145 "\tmust be set to use RNEMD, and can take any of these values:\n"
146 "\t\t%s.\n",
147 allowedFluxTypes.c_str());
148 painCave.isFatal = 1;
149 painCave.severity = OPENMD_ERROR;
150 simError();
151 }
152
153 // Determine Privileged Axis
154 const std::string privAxis = rnemdParams->getPrivilegedAxis();
155
156 if (privAxis == "x") {
157 rnemdAxisLabel_ = "x";
158 rnemdPrivilegedAxis_ = rnemdX;
159 } else if (privAxis == "y") {
160 rnemdAxisLabel_ = "y";
161 rnemdPrivilegedAxis_ = rnemdY;
162 } else {
163 rnemdAxisLabel_ = "z";
164 rnemdPrivilegedAxis_ = rnemdZ;
165 }
166
167 runTime_ = simParams->getRunTime();
168 statusTime_ = simParams->getStatusTime();
169
170 rnemdObjectSelection_ = rnemdParams->getObjectSelection();
171
172 bool hasSlabWidth = rnemdParams->haveSlabWidth();
173 bool hasSlabACenter = rnemdParams->haveSlabACenter();
174 bool hasSlabBCenter = rnemdParams->haveSlabBCenter();
175 bool hasSphereARadius = rnemdParams->haveSphereARadius();
176 bool hasSphereBRadius = rnemdParams->haveSphereBRadius();
177
178 hasSelectionA_ = rnemdParams->haveSelectionA();
179 hasSelectionB_ = rnemdParams->haveSelectionB();
180
181 hasDividingArea_ = rnemdParams->haveDividingArea();
182 dividingArea_ = rnemdParams->getDividingArea();
183
184 bool hasCoordinateOrigin = rnemdParams->haveCoordinateOrigin();
185 bool hasOutputFileName = rnemdParams->haveOutputFileName();
186 bool hasOutputFields = rnemdParams->haveOutputFields();
187 bool hasOutputSelection = rnemdParams->haveOutputSelection();
188
189 if (hasOutputSelection) {
190 outputSelection_ = rnemdParams->getOutputSelection();
191 } else {
192 outputSelection_ = rnemdObjectSelection_;
193 }
194
195 if (hasCoordinateOrigin) {
196 std::vector<RealType> co = rnemdParams->getCoordinateOrigin();
197 if (co.size() != 3) {
198 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
199 "RNEMD: Incorrect number of parameters specified for "
200 "coordinateOrigin.\n"
201 "\tthere should be 3 parameters, but %zu were specified.\n",
202 co.size());
203 painCave.isFatal = 1;
204 simError();
205 }
206 coordinateOrigin_.x() = co[0];
207 coordinateOrigin_.y() = co[1];
208 coordinateOrigin_.z() = co[2];
209 } else {
210 coordinateOrigin_ = V3Zero;
211 }
212
213 outputEvaluator_.loadScriptString(outputSelection_);
214 outputSeleMan_.setSelectionSet(outputEvaluator_.evaluate());
215
216 SelectionManager tempOutputSeleMan =
217 outputSeleMan_.replaceRigidBodiesWithAtoms();
218
219 AtomTypeSet osTypes = tempOutputSeleMan.getSelectedAtomTypes();
220 std::copy(osTypes.begin(), osTypes.end(), std::back_inserter(outputTypes_));
221
222 nBins_ = rnemdParams->getOutputBins();
223 binWidth_ = rnemdParams->getOutputBinWidth();
224
225 // Pre-load the OutputData
226 data_.resize(RNEMD::ENDINDEX);
227
228 OutputData z;
229 z.units = "Angstroms";
230 z.title = rnemdAxisLabel_;
231 for (unsigned int i = 0; i < nBins_; i++)
232 z.accumulator.push_back(
233 std::make_unique<AccumulatorView<RealAccumulator>>());
234 data_[Z] = std::move(z);
235 outputMap_["Z"] = Z;
236
237 OutputData r;
238 r.units = "Angstroms";
239 r.title = "R";
240 for (unsigned int i = 0; i < nBins_; i++)
241 r.accumulator.push_back(
242 std::make_unique<AccumulatorView<RealAccumulator>>());
243 data_[R] = std::move(r);
244 outputMap_["R"] = R;
245
246 OutputData temperature;
247 temperature.units = "K";
248 temperature.title = "Temperature";
249 for (unsigned int i = 0; i < nBins_; i++)
250 temperature.accumulator.push_back(
251 std::make_unique<AccumulatorView<RealAccumulator>>());
252 data_[TEMPERATURE] = std::move(temperature);
253 outputMap_["TEMPERATURE"] = TEMPERATURE;
254
255 OutputData velocity;
256 velocity.units = "angstroms/fs";
257 velocity.title = "Velocity";
258 for (unsigned int i = 0; i < nBins_; i++)
259 velocity.accumulator.push_back(
260 std::make_unique<AccumulatorView<Vector3dAccumulator>>());
261 data_[VELOCITY] = std::move(velocity);
262 outputMap_["VELOCITY"] = VELOCITY;
263
264 OutputData angularVelocity;
265 angularVelocity.units = "angstroms^2/fs";
266 angularVelocity.title = "AngularVelocity";
267 for (unsigned int i = 0; i < nBins_; i++)
268 angularVelocity.accumulator.push_back(
269 std::make_unique<AccumulatorView<Vector3dAccumulator>>());
270 data_[ANGULARVELOCITY] = std::move(angularVelocity);
271 outputMap_["ANGULARVELOCITY"] = ANGULARVELOCITY;
272
273 OutputData density;
274 density.units = "g cm^-3";
275 density.title = "Density";
276 for (unsigned int i = 0; i < nBins_; i++)
277 density.accumulator.push_back(
278 std::make_unique<AccumulatorView<RealAccumulator>>());
279 data_[DENSITY] = std::move(density);
280 outputMap_["DENSITY"] = DENSITY;
281
282 OutputData activity;
283 activity.units = "unitless";
284 activity.title = "Activity";
285 for (unsigned int i = 0; i < nBins_; i++)
286 activity.accumulator.push_back(
287 std::make_unique<AccumulatorView<StdVectorAccumulator>>());
288 data_[ACTIVITY] = std::move(activity);
289 outputMap_["ACTIVITY"] = ACTIVITY;
290
291 OutputData eField;
292 eField.units = "kcal/mol/angstroms/e";
293 eField.title = "Electric Field";
294 for (unsigned int i = 0; i < nBins_; i++)
295 eField.accumulator.push_back(
296 std::make_unique<AccumulatorView<Vector3dAccumulator>>());
297 data_[ELECTRICFIELD] = std::move(eField);
298 outputMap_["ELECTRICFIELD"] = ELECTRICFIELD;
299
300 OutputData ePot;
301 ePot.units = "kcal/mol/e";
302 ePot.title = "Electrostatic Potential";
303 for (unsigned int i = 0; i < nBins_; i++)
304 ePot.accumulator.push_back(
305 std::make_unique<AccumulatorView<RealAccumulator>>());
306 data_[ELECTROSTATICPOTENTIAL] = std::move(ePot);
307 outputMap_["ELECTROSTATICPOTENTIAL"] = ELECTROSTATICPOTENTIAL;
308
309 if (hasOutputFields) {
310 parseOutputFileFormat(rnemdParams->getOutputFields());
311 } else {
312 if (usePeriodicBoundaryConditions_)
313 outputMask_.set(Z);
314 else
315 outputMask_.set(R);
316 switch (rnemdFluxType_) {
317 case rnemdKE:
318 case rnemdRotKE:
319 case rnemdFullKE:
320 outputMask_.set(TEMPERATURE);
321 break;
322 case rnemdPx:
323 case rnemdPy:
324 outputMask_.set(VELOCITY);
325 break;
326 case rnemdPz:
327 case rnemdPvector:
328 outputMask_.set(VELOCITY);
329 outputMask_.set(DENSITY);
330 break;
331 case rnemdLx:
332 case rnemdLy:
333 case rnemdLz:
334 case rnemdLvector:
335 outputMask_.set(ANGULARVELOCITY);
336 break;
337 case rnemdKeLx:
338 case rnemdKeLy:
339 case rnemdKeLz:
340 case rnemdKeLvector:
341 outputMask_.set(TEMPERATURE);
342 outputMask_.set(ANGULARVELOCITY);
343 break;
344 case rnemdKePx:
345 case rnemdKePy:
346 outputMask_.set(TEMPERATURE);
347 outputMask_.set(VELOCITY);
348 break;
349 case rnemdKePvector:
350 outputMask_.set(TEMPERATURE);
351 outputMask_.set(VELOCITY);
352 outputMask_.set(DENSITY);
353 break;
354 default:
355 break;
356 }
357 }
358
359 if (hasOutputFileName) {
360 rnemdFileName_ = rnemdParams->getOutputFileName();
361 } else {
362 rnemdFileName_ = getPrefix(info->getFinalConfigFileName()) + ".rnemd";
363 }
364
365 // Exchange time should not be less than time step and should be a
366 // multiple of dt
367 exchangeTime_ = rnemdParams->getExchangeTime();
368 RealType dt = simParams->getDt();
369 RealType newET = std::ceil(exchangeTime_ / dt) * dt;
370
371 if (std::fabs(newET - exchangeTime_) > 1e-6) {
372 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
373 "RNEMD: The exchangeTime was reset to %lf,\n"
374 "\t\twhich is a multiple of dt, %lf.\n",
375 newET, dt);
376 painCave.isFatal = 0;
377 painCave.severity = OPENMD_WARNING;
378 simError();
379 exchangeTime_ = newET;
380 }
381
382 currentSnap_ = info->getSnapshotManager()->getCurrentSnapshot();
383 hmat_ = currentSnap_->getHmat();
384
385 // Set up the slab selection logic
386 std::ostringstream selectionAstream;
387 std::ostringstream selectionBstream;
388
389 if (hasSelectionA_) {
390 selectionA_ = rnemdParams->getSelectionA();
391 } else {
392 if (usePeriodicBoundaryConditions_) {
393 slabWidth_ =
394 hasSlabWidth ?
395 rnemdParams->getSlabWidth() :
396 hmat_(rnemdPrivilegedAxis_, rnemdPrivilegedAxis_) / 10.0;
397
398 slabACenter_ = hasSlabACenter ? rnemdParams->getSlabACenter() : 0.0;
399
400 selectionA_ = this->setSelection(slabACenter_);
401 } else {
402 if (hasSphereARadius)
403 sphereARadius_ = rnemdParams->getSphereARadius();
404 else {
405 // use an initial guess to the size of the inner slab to be 1/10 the
406 // radius of an approximately spherical hull:
407 Thermo thermo(info);
408 RealType hVol = thermo.getHullVolume();
409 sphereARadius_ =
410 0.1 * pow((3.0 * hVol / (4.0 * Constants::PI)), 1.0 / 3.0);
411 }
412 selectionAstream << "select r < " << sphereARadius_;
413 selectionA_ = selectionAstream.str();
414 }
415 }
416
417 if (hasSelectionB_) {
418 selectionB_ = rnemdParams->getSelectionB();
419 } else {
420 if (usePeriodicBoundaryConditions_) {
421 slabWidth_ =
422 hasSlabWidth ?
423 rnemdParams->getSlabWidth() :
424 hmat_(rnemdPrivilegedAxis_, rnemdPrivilegedAxis_) / 10.0;
425
426 slabBCenter_ =
427 hasSlabBCenter ?
428 rnemdParams->getSlabBCenter() :
429 hmat_(rnemdPrivilegedAxis_, rnemdPrivilegedAxis_) / 2.0;
430
431 selectionB_ = this->setSelection(slabBCenter_);
432 } else {
433 if (hasSphereBRadius) {
434 sphereBRadius_ = rnemdParams->getSphereBRadius();
435 selectionBstream << "select r > " << sphereBRadius_;
436 selectionB_ = selectionBstream.str();
437 } else {
438 selectionB_ = "select hull";
439 hasSelectionB_ = true;
440 }
441 }
442 }
443
444 // Static object evaluators
445 evaluator_.loadScriptString(rnemdObjectSelection_);
446 if (!evaluator_.isDynamic())
447 seleMan_.setSelectionSet(evaluator_.evaluate());
448
449 evaluatorA_.loadScriptString(selectionA_);
450 if (!evaluatorA_.isDynamic())
451 seleManA_.setSelectionSet(evaluatorA_.evaluate());
452
453 evaluatorB_.loadScriptString(selectionB_);
454 if (!evaluatorB_.isDynamic())
455 seleManB_.setSelectionSet(evaluatorB_.evaluate());
456
457 // Do some sanity checking
458 int selectionCount =
459 seleMan_.removeAtomsInRigidBodies().getSelectionCount();
460 int nIntegrable = info->getNGlobalIntegrableObjects();
461
462 if (selectionCount > nIntegrable) {
463 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
464 "RNEMD: The current objectSelection,\n"
465 "\t\t%s\n"
466 "\thas resulted in %d selected objects. However,\n"
467 "\tthe total number of integrable objects in the system\n"
468 "\tis only %d. This is almost certainly not what you want\n"
469 "\tto do.\n",
470 rnemdObjectSelection_.c_str(), selectionCount, nIntegrable);
471 painCave.isFatal = 0;
472 painCave.severity = OPENMD_WARNING;
473 simError();
474 }
475 }
476
477 RNEMD::~RNEMD() {
478 if (!doRNEMD_) return;
479#ifdef IS_MPI
480 if (worldRank == 0) {
481#endif
482 writeOutputFile();
483
484 rnemdFile_.close();
485#ifdef IS_MPI
486 }
487#endif
488 }
489
490 void RNEMD::getStarted() {
491 if (!doRNEMD_) return;
492 collectData();
493 writeOutputFile();
494 }
495
496 void RNEMD::doRNEMD() {
497 if (!doRNEMD_) return;
498 trialCount_++;
499 hmat_ = currentSnap_->getHmat();
500
501 // dynamic object evaluators:
502 evaluator_.loadScriptString(rnemdObjectSelection_);
503 if (evaluator_.isDynamic()) seleMan_.setSelectionSet(evaluator_.evaluate());
504
505 evaluatorA_.loadScriptString(selectionA_);
506 if (evaluatorA_.isDynamic())
507 seleManA_.setSelectionSet(evaluatorA_.evaluate());
508
509 evaluatorB_.loadScriptString(selectionB_);
510 if (evaluatorB_.isDynamic())
511 seleManB_.setSelectionSet(evaluatorB_.evaluate());
512
513 commonA_ = seleManA_ & seleMan_;
514 commonB_ = seleManB_ & seleMan_;
515
516 auto reducedCommonA = commonA_.removeAtomsInRigidBodies();
517 auto reducedCommonB = commonB_.removeAtomsInRigidBodies();
518
519 // Target exchange quantities (in each exchange) = flux * dividingArea *
520 // dt flux = target flux dividingArea = smallest dividing surface between
521 // the two regions dt = exchange time interval
522
523 RealType area = getDefaultDividingArea();
524
525 kineticTarget_ = kineticFlux_ * exchangeTime_ * area;
526 momentumTarget_ = momentumFluxVector_ * exchangeTime_ * area;
527 angularMomentumTarget_ = angularMomentumFluxVector_ * exchangeTime_ * area;
528 particleTarget_ = particleFlux_ * exchangeTime_ * area;
529
530 if (std::fabs(particleTarget_) > 1.0) {
531 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
532 "RNEMD: The current particleFlux,\n"
533 "\t\t%f\n"
534 "\thas resulted in a target particle exchange of %f.\n"
535 "\tThis is equivalent to moving more than one particle\n"
536 "\tduring each exchange. Please reduce your particleFlux.\n",
537 particleFlux_, particleTarget_);
538 painCave.isFatal = 1;
539 painCave.severity = OPENMD_ERROR;
540 simError();
541 }
542
543 if (rnemdFluxType_ == rnemdParticle || rnemdFluxType_ == rnemdParticleKE) {
544 SelectionManager tempCommonA = seleManA_ & outputSeleMan_;
545 SelectionManager tempCommonB = seleManB_ & outputSeleMan_;
546
547 auto reducedTempCommonA = tempCommonA.removeAtomsInRigidBodies();
548 auto reducedTempCommonB = tempCommonB.removeAtomsInRigidBodies();
549
550 this->doRNEMDImpl(reducedTempCommonA, reducedTempCommonB);
551 } else {
552 this->doRNEMDImpl(reducedCommonA, reducedCommonB);
553 }
554 }
555
556 void RNEMD::collectData() {
557 if (!doRNEMD_) return;
558 currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
559 hmat_ = currentSnap_->getHmat();
560
561 // collectData() can be called more frequently than the doRNEMD(), so use
562 // the computed area from the last exchange time:
563 RealType area = getDefaultDividingArea();
564 areaAccumulator_.add(area);
565
566 Vector3d u = angularMomentumFluxVector_;
567 u.normalize();
568
569 // throw an error if isDynamic instead?
570 if (outputEvaluator_.isDynamic()) {
571 outputSeleMan_.setSelectionSet(outputEvaluator_.evaluate());
572 }
573
574 int binNo {};
575 int typeIndex(-1);
576 RealType mass {};
577 Vector3d vel {};
578 Vector3d rPos {};
579 RealType KE {};
580 Vector3d L {};
581 Mat3x3d I {};
582 RealType r2 {};
583 Vector3d eField {};
584 int DOF {};
585
586 vector<RealType> binMass(nBins_, 0.0);
587 vector<Vector3d> binP(nBins_, V3Zero);
588 vector<RealType> binOmega(nBins_, 0.0);
589 vector<Vector3d> binL(nBins_, V3Zero);
590 vector<Mat3x3d> binI(nBins_);
591 vector<RealType> binKE(nBins_, 0.0);
592 vector<Vector3d> binEField(nBins_, V3Zero);
593 vector<int> binDOF(nBins_, 0);
594 vector<int> binCount(nBins_, 0);
595 vector<int> binEFieldCount(nBins_, 0);
596 vector<vector<int>> binTypeCounts;
597
598 if (outputMask_[ACTIVITY]) {
599 binTypeCounts.resize(nBins_);
600 for (unsigned int i = 0; i < nBins_; i++) {
601 binTypeCounts[i].resize(outputTypes_.size(), 0);
602 }
603 }
604
605 SimInfo::MoleculeIterator miter;
606 std::vector<StuntDouble*>::iterator iiter;
607 std::vector<AtomType*>::iterator at;
608 Molecule* mol;
609 StuntDouble* sd;
610 AtomType* atype;
611 ConstraintPair* consPair;
612 Molecule::ConstraintPairIterator cpi;
613
614 for (mol = info_->beginMolecule(miter); mol != NULL;
615 mol = info_->nextMolecule(miter)) {
616 for (sd = mol->beginIntegrableObject(iiter); sd != NULL;
617 sd = mol->nextIntegrableObject(iiter)) {
618 if (outputSeleMan_.isSelected(sd)) {
619 Vector3d pos = sd->getPos();
620 binNo = getBin(pos);
621
622 mass = sd->getMass();
623 vel = sd->getVel();
624 rPos = sd->getPos() - coordinateOrigin_;
625 KE = 0.5 * mass * vel.lengthSquare();
626 DOF = 3;
627
628 if (sd->isDirectional()) {
629 Vector3d angMom = sd->getJ();
630 Mat3x3d Ia = sd->getI();
631 if (sd->isLinear()) {
632 int i = sd->linearAxis();
633 int j = (i + 1) % 3;
634 int k = (i + 2) % 3;
635 KE += 0.5 * (angMom[j] * angMom[j] / Ia(j, j) +
636 angMom[k] * angMom[k] / Ia(k, k));
637 DOF += 2;
638 } else {
639 KE += 0.5 * (angMom[0] * angMom[0] / Ia(0, 0) +
640 angMom[1] * angMom[1] / Ia(1, 1) +
641 angMom[2] * angMom[2] / Ia(2, 2));
642 DOF += 3;
643 }
644 }
645
646 L = mass * cross(rPos, vel);
647 I = outProduct(rPos, rPos) * mass;
648 r2 = rPos.lengthSquare();
649 I(0, 0) += mass * r2;
650 I(1, 1) += mass * r2;
651 I(2, 2) += mass * r2;
652
653 if (outputMask_[ACTIVITY]) {
654 typeIndex = -1;
655 if (sd->isRigidBody()) {
656 int atomBinNo;
657 RigidBody* rb = static_cast<RigidBody*>(sd);
658 std::vector<Atom*>::iterator ai;
659 Atom* atom;
660 for (atom = rb->beginAtom(ai); atom != NULL;
661 atom = rb->nextAtom(ai)) {
662 atomBinNo = getBin(atom->getPos());
663
664 atype = static_cast<Atom*>(atom)->getAtomType();
665 at = std::find(outputTypes_.begin(), outputTypes_.end(), atype);
666 if (at != outputTypes_.end()) {
667 typeIndex = std::distance(outputTypes_.begin(), at);
668 }
669
670 if (atomBinNo >= 0 && atomBinNo < int(nBins_)) {
671 if (typeIndex != -1) binTypeCounts[atomBinNo][typeIndex]++;
672 }
673 }
674 } else if (sd->isAtom()) {
675 atype = static_cast<Atom*>(sd)->getAtomType();
676 at = std::find(outputTypes_.begin(), outputTypes_.end(), atype);
677 if (at != outputTypes_.end()) {
678 typeIndex = std::distance(outputTypes_.begin(), at);
679 }
680 }
681 }
682
683 if (binNo >= 0 && binNo < int(nBins_)) {
684 binCount[binNo]++;
685 binMass[binNo] += mass;
686 binP[binNo] += mass * vel;
687 binKE[binNo] += KE;
688 binI[binNo] += I;
689 binL[binNo] += L;
690 binDOF[binNo] += DOF;
691
692 if (outputMask_[ACTIVITY] && typeIndex != -1)
693 binTypeCounts[binNo][typeIndex]++;
694 }
695 }
696
697 // Calculate the electric field (kcal/mol/e/Angstrom) for all atoms in
698 // the box
699 if (outputMask_[ELECTRICFIELD]) {
700 int atomBinNo;
701 if (sd->isRigidBody()) {
702 RigidBody* rb = static_cast<RigidBody*>(sd);
703 std::vector<Atom*>::iterator ai;
704 Atom* atom;
705 for (atom = rb->beginAtom(ai); atom != NULL;
706 atom = rb->nextAtom(ai)) {
707 atomBinNo = getBin(atom->getPos());
708 eField = atom->getElectricField();
709
710 binEFieldCount[atomBinNo]++;
711 binEField[atomBinNo] += eField;
712 }
713 } else {
714 eField = sd->getElectricField();
715 atomBinNo = getBin(sd->getPos());
716
717 binEFieldCount[atomBinNo]++;
718 binEField[atomBinNo] += eField;
719 }
720 }
721 }
722
723 // we need to subtract out degrees of freedom from constraints
724 // belonging in this bin:
725 if (outputSeleMan_.isSelected(mol)) {
726 for (consPair = mol->beginConstraintPair(cpi); consPair != NULL;
727 consPair = mol->nextConstraintPair(cpi)) {
728 Vector3d posA = consPair->getConsElem1()->getPos();
729 Vector3d posB = consPair->getConsElem2()->getPos();
730
731 if (usePeriodicBoundaryConditions_) {
732 currentSnap_->wrapVector(posA);
733 currentSnap_->wrapVector(posB);
734 }
735
736 Vector3d coc = 0.5 * (posA + posB);
737 int binCons = getBin(coc);
738 binDOF[binCons] -= 1;
739 }
740 }
741 }
742
743#ifdef IS_MPI
744 for (unsigned int i = 0; i < nBins_; i++) {
745 MPI_Allreduce(MPI_IN_PLACE, &binCount[i], 1, MPI_INT, MPI_SUM,
746 MPI_COMM_WORLD);
747 MPI_Allreduce(MPI_IN_PLACE, &binMass[i], 1, MPI_REALTYPE, MPI_SUM,
748 MPI_COMM_WORLD);
749 MPI_Allreduce(MPI_IN_PLACE, binP[i].getArrayPointer(), 3, MPI_REALTYPE,
750 MPI_SUM, MPI_COMM_WORLD);
751 MPI_Allreduce(MPI_IN_PLACE, binL[i].getArrayPointer(), 3, MPI_REALTYPE,
752 MPI_SUM, MPI_COMM_WORLD);
753 MPI_Allreduce(MPI_IN_PLACE, binI[i].getArrayPointer(), 9, MPI_REALTYPE,
754 MPI_SUM, MPI_COMM_WORLD);
755 MPI_Allreduce(MPI_IN_PLACE, &binKE[i], 1, MPI_REALTYPE, MPI_SUM,
756 MPI_COMM_WORLD);
757 MPI_Allreduce(MPI_IN_PLACE, &binDOF[i], 1, MPI_INT, MPI_SUM,
758 MPI_COMM_WORLD);
759 MPI_Allreduce(MPI_IN_PLACE, &binEFieldCount[i], 1, MPI_INT, MPI_SUM,
760 MPI_COMM_WORLD);
761
762 if (outputMask_[ELECTRICFIELD]) {
763 MPI_Allreduce(MPI_IN_PLACE, binEField[i].getArrayPointer(), 3,
764 MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
765 }
766 if (outputMask_[ACTIVITY]) {
767 MPI_Allreduce(MPI_IN_PLACE, &binTypeCounts[i][0], outputTypes_.size(),
768 MPI_INT, MPI_SUM, MPI_COMM_WORLD);
769 }
770 }
771#endif
772
773 Vector3d omega;
774 RealType z, r, temp, binVolume, den(0.0), dz(0.0);
775 std::vector<RealType> nden(outputTypes_.size(), 0.0);
776 RealType boxVolume = currentSnap_->getVolume();
777 RealType ePot(0.0);
778
779 for (unsigned int i = 0; i < nBins_; i++) {
780 if (usePeriodicBoundaryConditions_) {
781 z = (((RealType)i + 0.5) / (RealType)nBins_) *
782 hmat_(rnemdPrivilegedAxis_, rnemdPrivilegedAxis_);
783 data_[Z].accumulator[i]->add(z);
784
785 binVolume = boxVolume / nBins_;
786 dz = hmat_(rnemdPrivilegedAxis_, rnemdPrivilegedAxis_) /
787 (RealType)nBins_;
788 } else {
789 r = (((RealType)i + 0.5) * binWidth_);
790 data_[R].accumulator[i]->add(r);
791
792 RealType rinner = (RealType)i * binWidth_;
793 RealType router = (RealType)(i + 1) * binWidth_;
794 binVolume =
795 (4.0 * Constants::PI * (pow(router, 3) - pow(rinner, 3))) / 3.0;
796 }
797
798 // The calculations of the following properties are done regardless
799 // of whether or not the selected species are present in the bin
800 if (outputMask_[ELECTRICFIELD] && binEFieldCount[i] > 0) {
801 eField = binEField[i] / RealType(binEFieldCount[i]);
802 data_[ELECTRICFIELD].accumulator[i]->add(eField);
803 }
804
805 if (outputMask_[ELECTROSTATICPOTENTIAL]) {
806 if (usePeriodicBoundaryConditions_ && binEFieldCount[i] > 0) {
807 ePot += eField[rnemdPrivilegedAxis_] * dz;
808 data_[ELECTROSTATICPOTENTIAL].accumulator[i]->add(ePot);
809 }
810 }
811
812 // For the following properties, zero should be added if the selected
813 // species is not present in the bin
814 if (outputMask_[DENSITY]) {
815 den = binMass[i] * Constants::densityConvert / binVolume;
816 data_[DENSITY].accumulator[i]->add(den);
817 }
818
819 if (outputMask_[ACTIVITY]) {
820 for (unsigned int j = 0; j < outputTypes_.size(); j++) {
821 nden[j] = (binTypeCounts[i][j] / binVolume) *
822 Constants::concentrationConvert;
823 }
824 data_[ACTIVITY].accumulator[i]->add(nden);
825 }
826
827 if (binCount[i] > 0) {
828 // The calculations of the following properties are meaningless if
829 // the selected species is not found in the bin
830 if (outputMask_[VELOCITY]) {
831 vel = binP[i] / binMass[i];
832 data_[VELOCITY].accumulator[i]->add(vel);
833 }
834
835 if (outputMask_[ANGULARVELOCITY]) {
836 omega = binI[i].inverse() * binL[i];
837 data_[ANGULARVELOCITY].accumulator[i]->add(omega);
838 }
839
840 if (outputMask_[TEMPERATURE]) {
841 if (binDOF[i] > 0) {
842 temp = 2.0 * binKE[i] /
843 (binDOF[i] * Constants::kb * Constants::energyConvert);
844 data_[TEMPERATURE].accumulator[i]->add(temp);
845 } else {
846 std::cerr << "No degrees of freedom in this bin?\n";
847 }
848 }
849 }
850 }
851
852 hasData_ = true;
853 }
854
855 void RNEMD::writeOutputFile() {
856 if (!doRNEMD_) return;
857 if (!hasData_) return;
858
859#ifdef IS_MPI
860 // If we're the primary node, should we print out the results
861 int worldRank;
862 MPI_Comm_rank(MPI_COMM_WORLD, &worldRank);
863
864 if (worldRank == 0) {
865#endif
866 rnemdFile_.open(rnemdFileName_.c_str(), std::ios::out | std::ios::trunc);
867
868 if (!rnemdFile_) {
869 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
870 "Could not open \"%s\" for RNEMD output.\n",
871 rnemdFileName_.c_str());
872 painCave.isFatal = 1;
873 simError();
874 }
875
876 RealType time = currentSnap_->getTime();
877 RealType avgArea = areaAccumulator_.getAverage();
878
879 RealType Jz(0.0);
880 Vector3d JzP(V3Zero);
881 Vector3d JzL(V3Zero);
882 RealType Jpart(0.0);
883
884 if (time >= info_->getSimParams()->getDt()) {
885 Jz = kineticExchange_ / (time * avgArea) / Constants::energyConvert;
886 JzP = momentumExchange_ / (time * avgArea);
887 JzL = angularMomentumExchange_ / (time * avgArea);
888 Jpart = particleExchange_ / (time * avgArea);
889 }
890
891 rnemdFile_ << "#######################################################\n";
892 rnemdFile_ << "# RNEMD {\n";
893 rnemdFile_ << "# exchangeMethod = \"" << rnemdMethodLabel_ << "\";\n";
894 rnemdFile_ << "# fluxType = \"" << rnemdFluxTypeLabel_ << "\";\n";
895
896 if (usePeriodicBoundaryConditions_)
897 rnemdFile_ << "# privilegedAxis = " << rnemdAxisLabel_ << ";\n";
898
899 rnemdFile_ << "# exchangeTime = " << exchangeTime_ << ";\n";
900 rnemdFile_ << "# objectSelection = \"" << rnemdObjectSelection_
901 << "\";\n";
902 rnemdFile_ << "# selectionA = \"" << selectionA_ << "\";\n";
903 rnemdFile_ << "# selectionB = \"" << selectionB_ << "\";\n";
904 rnemdFile_ << "# outputSelection = \"" << outputSelection_ << "\";\n";
905 rnemdFile_ << "# }\n";
906 rnemdFile_ << "#######################################################\n";
907 rnemdFile_ << "# RNEMD report:\n";
908 rnemdFile_ << "# running time = " << time << " fs\n";
909 rnemdFile_ << "# Target flux:\n";
910 rnemdFile_ << "# kinetic = "
911 << kineticFlux_ / Constants::energyConvert
912 << " (kcal/mol/A^2/fs)\n";
913 rnemdFile_ << "# momentum = " << momentumFluxVector_
914 << " (amu/A/fs^2)\n";
915 rnemdFile_ << "# angular momentum = " << angularMomentumFluxVector_
916 << " (amu/A^2/fs^2)\n";
917 rnemdFile_ << "# particle = " << particleFlux_
918 << " (particles/A^2/fs)\n";
919
920 rnemdFile_ << "# Target one-time exchanges:\n";
921 rnemdFile_ << "# kinetic = "
922 << kineticTarget_ / Constants::energyConvert
923 << " (kcal/mol)\n";
924 rnemdFile_ << "# momentum = " << momentumTarget_
925 << " (amu*A/fs)\n";
926 rnemdFile_ << "# angular momentum = " << angularMomentumTarget_
927 << " (amu*A^2/fs)\n";
928 rnemdFile_ << "# particle = " << particleTarget_
929 << " (particles)\n";
930 rnemdFile_ << "# Actual exchange totals:\n";
931 rnemdFile_ << "# kinetic = "
932 << kineticExchange_ / Constants::energyConvert
933 << " (kcal/mol)\n";
934 rnemdFile_ << "# momentum = " << momentumExchange_
935 << " (amu*A/fs)\n";
936 rnemdFile_ << "# angular momentum = " << angularMomentumExchange_
937 << " (amu*A^2/fs)\n";
938 rnemdFile_ << "# particles = " << particleExchange_
939 << " (particles)\n";
940
941 rnemdFile_ << "# Actual flux:\n";
942 rnemdFile_ << "# kinetic = " << Jz << " (kcal/mol/A^2/fs)\n";
943 rnemdFile_ << "# momentum = " << JzP << " (amu/A/fs^2)\n";
944 rnemdFile_ << "# angular momentum = " << JzL << " (amu/A^2/fs^2)\n";
945 rnemdFile_ << "# particle = " << Jpart
946 << " (particles/A^2/fs)\n";
947 rnemdFile_ << "# Exchange statistics:\n";
948 rnemdFile_ << "# attempted = " << trialCount_ << "\n";
949 rnemdFile_ << "# failed = " << failTrialCount_ << "\n";
950 if (rnemdMethodLabel_ == "NIVS") {
951 rnemdFile_ << "# NIVS root-check errors = " << failRootCount_ << "\n";
952 }
953 rnemdFile_ << "#######################################################\n";
954
955 // write title
956 rnemdFile_ << "#";
957 for (unsigned int i = 0; i < outputMask_.size(); ++i) {
958 if (outputMask_[i]) {
959 rnemdFile_ << "\t" << data_[i].title << "(" << data_[i].units << ")";
960
961 // add some extra tabs for column alignment
962 if (data_[i].accumulator[0]->getType() ==
963 std::type_index(typeid(Vector3d))) {
964 rnemdFile_ << "\t\t";
965 }
966
967 if (data_[i].accumulator[0]->getType() ==
968 std::type_index(typeid(std::vector<RealType>))) {
969 rnemdFile_ << "(";
970 for (unsigned int type = 0; type < outputTypes_.size(); type++) {
971 rnemdFile_ << outputTypes_[type]->getName() << "\t";
972 }
973 rnemdFile_ << ")\t";
974 }
975 }
976 }
977
978 rnemdFile_ << std::endl;
979
980 std::vector<int> nonEmptyAccumulators(nBins_);
981 int numberOfAccumulators {};
982
983 for (unsigned int i = 0; i < outputMask_.size(); ++i) {
984 if (outputMask_[i]) {
985 for (unsigned int bin = 0; bin < nBins_; bin++) {
986 nonEmptyAccumulators[bin] +=
987 static_cast<int>(data_[i].accumulator[bin]->getCount() != 0);
988 }
989
990 numberOfAccumulators++;
991 }
992 }
993
994 rnemdFile_.precision(8);
995
996 for (unsigned int bin = 0; bin < nBins_; bin++) {
997 if (nonEmptyAccumulators[bin] == numberOfAccumulators) {
998 for (unsigned int i = 0; i < outputMask_.size(); ++i) {
999 if (outputMask_[i]) {
1000 std::string message =
1001 "RNEMD detected a numerical error writing: " +
1002 data_[i].title + " for bin " + std::to_string(bin);
1003
1004 data_[i].accumulator[bin]->writeData(rnemdFile_, message);
1005 }
1006 }
1007
1008 rnemdFile_ << std::endl;
1009 }
1010 }
1011
1012 rnemdFile_ << "#######################################################\n";
1013 rnemdFile_ << "# 95% confidence intervals in those quantities follow:\n";
1014 rnemdFile_ << "#######################################################\n";
1015
1016 for (unsigned int bin = 0; bin < nBins_; bin++) {
1017 if (nonEmptyAccumulators[bin] == numberOfAccumulators) {
1018 rnemdFile_ << "#";
1019
1020 for (unsigned int i = 0; i < outputMask_.size(); ++i) {
1021 if (outputMask_[i]) {
1022 std::string message =
1023 "RNEMD detected a numerical error writing: " +
1024 data_[i].title + " std. dev. for bin " + std::to_string(bin);
1025
1026 data_[i].accumulator[bin]->writeErrorBars(rnemdFile_, message);
1027 }
1028 }
1029
1030 rnemdFile_ << std::endl;
1031 }
1032 }
1033
1034 rnemdFile_.flush();
1035 rnemdFile_.close();
1036#ifdef IS_MPI
1037 }
1038#endif
1039 }
1040
1041 void RNEMD::setKineticFlux(RealType kineticFlux) {
1042 // convert the kcal / mol / Angstroms^2 / fs values in the md file
1043 // into amu / fs^3:
1044 kineticFlux_ = kineticFlux * Constants::energyConvert;
1045 }
1046
1047 void RNEMD::setParticleFlux(RealType particleFlux) {
1048 particleFlux_ = particleFlux;
1049 }
1050
1051 void RNEMD::setMomentumFluxVector(
1052 const std::vector<RealType>& momentumFluxVector) {
1053 if (momentumFluxVector.size() != 3) {
1054 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
1055 "RNEMD: Incorrect number of parameters specified for "
1056 "momentumFluxVector.\n"
1057 "\tthere should be 3 parameters, but %zu were specified.\n",
1058 momentumFluxVector.size());
1059 painCave.isFatal = 1;
1060 simError();
1061 }
1062
1063 momentumFluxVector_.x() = momentumFluxVector[0];
1064 momentumFluxVector_.y() = momentumFluxVector[1];
1065 momentumFluxVector_.z() = momentumFluxVector[2];
1066 }
1067
1068 void RNEMD::setAngularMomentumFluxVector(
1069 const std::vector<RealType>& angularMomentumFluxVector) {
1070 if (angularMomentumFluxVector.size() != 3) {
1071 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
1072 "RNEMD: Incorrect number of parameters specified for "
1073 "angularMomentumFluxVector.\n"
1074 "\tthere should be 3 parameters, but %zu were specified.\n",
1075 angularMomentumFluxVector.size());
1076 painCave.isFatal = 1;
1077 simError();
1078 }
1079
1080 angularMomentumFluxVector_.x() = angularMomentumFluxVector[0];
1081 angularMomentumFluxVector_.y() = angularMomentumFluxVector[1];
1082 angularMomentumFluxVector_.z() = angularMomentumFluxVector[2];
1083 }
1084
1085 void RNEMD::parseOutputFileFormat(const std::string& format) {
1086 if (!doRNEMD_) return;
1087 StringTokenizer tokenizer(format, " ,;|\t\n\r");
1088
1089 while (tokenizer.hasMoreTokens()) {
1090 std::string token(tokenizer.nextToken());
1091 toUpper(token);
1092 OutputMapType::iterator i = outputMap_.find(token);
1093 if (i != outputMap_.end()) {
1094 outputMask_.set(i->second);
1095 } else {
1096 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
1097 "RNEMD::parseOutputFileFormat: %s is not a recognized\n"
1098 "\toutputFileFormat keyword.\n",
1099 token.c_str());
1100 painCave.isFatal = 0;
1101 painCave.severity = OPENMD_ERROR;
1102 simError();
1103 }
1104 }
1105 }
1106
1107 std::string RNEMD::setSelection(RealType& slabCenter) {
1108 bool printSlabCenterWarning {false};
1109
1110 Vector3d tempSlabCenter {V3Zero};
1111 tempSlabCenter[rnemdPrivilegedAxis_] = slabCenter;
1112
1113 RealType hmat_2 = hmat_(rnemdPrivilegedAxis_, rnemdPrivilegedAxis_) / 2.0;
1114
1115 if (slabCenter > hmat_2) {
1116 currentSnap_->wrapVector(tempSlabCenter);
1117 printSlabCenterWarning = true;
1118 } else if (slabCenter < -hmat_2) {
1119 currentSnap_->wrapVector(tempSlabCenter);
1120 printSlabCenterWarning = true;
1121 }
1122
1123 if (printSlabCenterWarning) {
1124 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
1125 "The given slab center was set to %0.2f. In the wrapped "
1126 "coordinates\n"
1127 "\t[-Hmat/2, +Hmat/2], this has been remapped to %0.2f.\n",
1128 slabCenter, tempSlabCenter[rnemdPrivilegedAxis_]);
1129 painCave.isFatal = 0;
1130 painCave.severity = OPENMD_WARNING;
1131 simError();
1132
1133 slabCenter = tempSlabCenter[rnemdPrivilegedAxis_];
1134 }
1135
1136 Vector3d leftSlab {V3Zero};
1137 const RealType& leftSlabBoundary = leftSlab[rnemdPrivilegedAxis_];
1138 leftSlab[rnemdPrivilegedAxis_] = slabCenter - 0.5 * slabWidth_;
1139 currentSnap_->wrapVector(leftSlab);
1140
1141 Vector3d rightSlab {V3Zero};
1142 const RealType& rightSlabBoundary = rightSlab[rnemdPrivilegedAxis_];
1143 rightSlab[rnemdPrivilegedAxis_] = slabCenter + 0.5 * slabWidth_;
1144 currentSnap_->wrapVector(rightSlab);
1145
1146 std::ostringstream selectionStream;
1147
1148 selectionStream << "select wrapped" << rnemdAxisLabel_
1149 << " >= " << leftSlabBoundary;
1150
1151 if (leftSlabBoundary > rightSlabBoundary)
1152 selectionStream << " || wrapped" << rnemdAxisLabel_ << " < "
1153 << rightSlabBoundary;
1154 else
1155 selectionStream << " && wrapped" << rnemdAxisLabel_ << " < "
1156 << rightSlabBoundary;
1157
1158 return selectionStream.str();
1159 }
1160
1161 RealType RNEMD::getDefaultDividingArea() {
1162 if (hasDividingArea_) return dividingArea_;
1163
1164 Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot();
1165
1166 if (hasSelectionA_) {
1167 if (evaluatorA_.hasSurfaceArea()) {
1168 areaA_ = evaluatorA_.getSurfaceArea();
1169 volumeA_ = evaluatorA_.getVolume();
1170 } else {
1171 int isd;
1172 StuntDouble* sd;
1173 std::vector<StuntDouble*> aSites;
1174 seleManA_.setSelectionSet(evaluatorA_.evaluate());
1175 for (sd = seleManA_.beginSelected(isd); sd != NULL;
1176 sd = seleManA_.nextSelected(isd)) {
1177 aSites.push_back(sd);
1178 }
1179#if defined(HAVE_QHULL)
1180 ConvexHull* surfaceMeshA = new ConvexHull();
1181 surfaceMeshA->computeHull(aSites);
1182 areaA_ = surfaceMeshA->getArea();
1183 volumeA_ = surfaceMeshA->getVolume();
1184 delete surfaceMeshA;
1185#else
1186 snprintf(
1187 painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
1188 "RNEMD::getDividingArea : Hull calculation is not possible\n"
1189 "\twithout libqhull. Please rebuild OpenMD with qhull enabled.");
1190 painCave.severity = OPENMD_ERROR;
1191 painCave.isFatal = 1;
1192 simError();
1193#endif
1194 }
1195 } else {
1196 if (usePeriodicBoundaryConditions_) {
1197 // in periodic boundaries, the surface area is twice the
1198 // area of the current box, normal to the privileged axis:
1199 switch (rnemdPrivilegedAxis_) {
1200 case rnemdX:
1201 areaA_ = 2.0 * snap->getYZarea();
1202 break;
1203 case rnemdY:
1204 areaA_ = 2.0 * snap->getXZarea();
1205 break;
1206 case rnemdZ:
1207 default:
1208 areaA_ = 2.0 * snap->getXYarea();
1209 }
1210
1211 volumeA_ = areaA_ * slabWidth_;
1212 } else {
1213 // in non-periodic simulations, without explicitly setting
1214 // selections, the sphere radius sets the surface area of the
1215 // dividing surface:
1216 areaA_ = 4.0 * Constants::PI * std::pow(sphereARadius_, 2);
1217 volumeA_ = 4.0 * Constants::PI * std::pow(sphereARadius_, 3) / 3.0;
1218 }
1219 }
1220
1221 if (hasSelectionB_) {
1222 if (evaluatorB_.hasSurfaceArea()) {
1223 areaB_ = evaluatorB_.getSurfaceArea();
1224 volumeB_ = evaluatorB_.getVolume();
1225 } else {
1226 int isd;
1227 StuntDouble* sd;
1228 std::vector<StuntDouble*> bSites;
1229 seleManB_.setSelectionSet(evaluatorB_.evaluate());
1230 for (sd = seleManB_.beginSelected(isd); sd != NULL;
1231 sd = seleManB_.nextSelected(isd)) {
1232 bSites.push_back(sd);
1233 }
1234
1235#if defined(HAVE_QHULL)
1236 ConvexHull* surfaceMeshB = new ConvexHull();
1237 surfaceMeshB->computeHull(bSites);
1238 areaB_ = surfaceMeshB->getArea();
1239 volumeB_ = surfaceMeshB->getVolume();
1240 delete surfaceMeshB;
1241#else
1242 snprintf(
1243 painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
1244 "RNEMD::getDividingArea : Hull calculation is not possible\n"
1245 "\twithout libqhull. Please rebuild OpenMD with qhull enabled.");
1246 painCave.severity = OPENMD_ERROR;
1247 painCave.isFatal = 1;
1248 simError();
1249#endif
1250 }
1251 } else {
1252 if (usePeriodicBoundaryConditions_) {
1253 // in periodic boundaries, the surface area is twice the
1254 // area of the current box, normal to the privileged axis:
1255 switch (rnemdPrivilegedAxis_) {
1256 case rnemdX:
1257 areaB_ = 2.0 * snap->getYZarea();
1258 break;
1259 case rnemdY:
1260 areaB_ = 2.0 * snap->getXZarea();
1261 break;
1262 case rnemdZ:
1263 default:
1264 areaB_ = 2.0 * snap->getXYarea();
1265 }
1266
1267 volumeB_ = areaB_ * slabWidth_;
1268 } else {
1269 // in non-periodic simulations, without explicitly setting
1270 // selections, but if a sphereBradius has been set, just use that:
1271 areaB_ = 4.0 * Constants::PI * pow(sphereBRadius_, 2);
1272 Thermo thermo(info_);
1273 RealType hVol = thermo.getHullVolume();
1274 volumeB_ = hVol - 4.0 * Constants::PI * pow(sphereBRadius_, 3) / 3.0;
1275 }
1276 }
1277
1278 dividingArea_ = min(areaA_, areaB_);
1279 hasDividingArea_ = true;
1280 return dividingArea_;
1281 }
1282
1283 int RNEMD::getBin(Vector3d pos) {
1284 if (usePeriodicBoundaryConditions_) {
1285 currentSnap_->wrapVector(pos);
1286
1287 return int(nBins_ *
1288 (pos[rnemdPrivilegedAxis_] /
1289 hmat_(rnemdPrivilegedAxis_, rnemdPrivilegedAxis_) +
1290 0.5)) %
1291 nBins_;
1292 } else {
1293 Vector3d rPos = pos - coordinateOrigin_;
1294 return int(rPos.length() / binWidth_);
1295 }
1296 }
1297} // namespace OpenMD::RNEMD
OpenMD::AtomType
AtomType is what OpenMD looks to for unchanging data about an atom.
Definition AtomType.hpp:66
OpenMD::ConstraintElem::getPos
Vector3d getPos()
Returns the current position of this stuntdouble.
Definition ConstraintElem.hpp:177
OpenMD::ConstraintPair::getConsElem1
ConstraintElem * getConsElem1()
Return the first constraint elemet.
Definition ConstraintPair.hpp:69
OpenMD::ConstraintPair::getConsElem2
ConstraintElem * getConsElem2()
Retunr the second constraint element.
Definition ConstraintPair.hpp:72
OpenMD::ForceManager
ForceManager is responsible for calculating both the short range (bonded) interactions and long range...
Definition ForceManager.hpp:80
OpenMD::RigidBody
Definition RigidBody.hpp:62
OpenMD::SelectionManager::getSelectedAtomTypes
AtomTypeSet getSelectedAtomTypes()
getSelectedAtomTypes
Definition SelectionManager.cpp:547
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::getNGlobalIntegrableObjects
int getNGlobalIntegrableObjects()
Returns the total number of integrable objects (total number of rigid bodies plus the total number of...
Definition SimInfo.hpp:139
OpenMD::SimInfo::getSnapshotManager
SnapshotManager * getSnapshotManager()
Returns the snapshot manager.
Definition SimInfo.hpp:248
OpenMD::Snapshot
The Snapshot class is a repository storing dynamic data during a Simulation.
Definition Snapshot.hpp:147
OpenMD::Snapshot::getHmat
Mat3x3d getHmat()
Returns the H-Matrix.
Definition Snapshot.cpp:214
OpenMD::SnapshotManager::getCurrentSnapshot
Snapshot * getCurrentSnapshot()
Returns the pointer of current snapshot.
Definition SnapshotManager.hpp:89
OpenMD::StringTokenizer
The string tokenizer class allows an application to break a string into tokens The set of delimiters ...
Definition StringTokenizer.hpp:71
OpenMD::StuntDouble
"Don't move, or you're dead! Stand up! Captain, we've got them!"
Definition StuntDouble.hpp:91
OpenMD::StuntDouble::getVel
Vector3d getVel()
Returns the current velocity of this stuntDouble.
Definition StuntDouble.hpp:256
OpenMD::StuntDouble::linearAxis
int linearAxis()
Returns the linear axis of the rigidbody, atom and directional atom will always return -1.
Definition StuntDouble.hpp:1503
OpenMD::StuntDouble::getMass
RealType getMass()
Returns the mass of this stuntDouble.
Definition StuntDouble.hpp:1506
OpenMD::StuntDouble::getI
virtual Mat3x3d getI()=0
Returns the inertia tensor of this stuntDouble.
OpenMD::StuntDouble::isLinear
bool isLinear()
Tests the if this stuntDouble is a linear rigidbody.
Definition StuntDouble.hpp:1493
OpenMD::StuntDouble::getPos
Vector3d getPos()
Returns the current position of this stuntDouble.
Definition StuntDouble.hpp:196
OpenMD::StuntDouble::getElectricField
Vector3d getElectricField()
Returns the current electric field of this stuntDouble.
Definition StuntDouble.hpp:1200
OpenMD::StuntDouble::isRigidBody
bool isRigidBody()
Tests if this stuntDouble is a rigid body.
Definition StuntDouble.hpp:152
OpenMD::StuntDouble::getJ
Vector3d getJ()
Returns the current angular momentum of this stuntDouble (body -fixed).
Definition StuntDouble.hpp:374
OpenMD::StuntDouble::isAtom
bool isAtom()
Tests if this stuntDouble is an atom.
Definition StuntDouble.hpp:139
OpenMD::StuntDouble::isDirectional
bool isDirectional()
Tests if this stuntDouble is a directional one.
Definition StuntDouble.hpp:159
OpenMD::Vector::length
Real length()
Returns the length of this vector.
Definition Vector.hpp:393
OpenMD::Vector::normalize
void normalize()
Normalizes this vector in place.
Definition Vector.hpp:402
OpenMD::Vector::lengthSquare
Real lengthSquare()
Returns the squared length of this vector.
Definition Vector.hpp:399
OpenMD::Vector::add
void add(const Vector< Real, Dim > &v1)
Sets the value of this vector to the sum of itself and v1 (*this += v1).
Definition Vector.hpp:218
OpenMD::cross
Vector3< Real > cross(const Vector3< Real > &v1, const Vector3< Real > &v2)
Returns the cross product of two Vectors.
Definition Vector3.hpp:136
OpenMD::getPrefix
std::string getPrefix(const std::string &str)
Definition StringUtils.cpp:169