OpenMD 3.0
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 %lu 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 AtomTypeSet osTypes = outputSeleMan_.getSelectedAtomTypes();
216 std::copy(osTypes.begin(), osTypes.end(), std::back_inserter(outputTypes_));
217
218 nBins_ = rnemdParams->getOutputBins();
219 binWidth_ = rnemdParams->getOutputBinWidth();
220
221 // Pre-load the OutputData
222 data_.resize(RNEMD::ENDINDEX);
223
224 OutputData z;
225 z.units = "Angstroms";
226 z.title = rnemdAxisLabel_;
227 for (unsigned int i = 0; i < nBins_; i++)
228 z.accumulator.push_back(
229 std::make_unique<AccumulatorView<RealAccumulator>>());
230 data_[Z] = std::move(z);
231 outputMap_["Z"] = Z;
232
233 OutputData r;
234 r.units = "Angstroms";
235 r.title = "R";
236 for (unsigned int i = 0; i < nBins_; i++)
237 r.accumulator.push_back(
238 std::make_unique<AccumulatorView<RealAccumulator>>());
239 data_[R] = std::move(r);
240 outputMap_["R"] = R;
241
242 OutputData temperature;
243 temperature.units = "K";
244 temperature.title = "Temperature";
245 for (unsigned int i = 0; i < nBins_; i++)
246 temperature.accumulator.push_back(
247 std::make_unique<AccumulatorView<RealAccumulator>>());
248 data_[TEMPERATURE] = std::move(temperature);
249 outputMap_["TEMPERATURE"] = TEMPERATURE;
250
251 OutputData velocity;
252 velocity.units = "angstroms/fs";
253 velocity.title = "Velocity";
254 for (unsigned int i = 0; i < nBins_; i++)
255 velocity.accumulator.push_back(
256 std::make_unique<AccumulatorView<Vector3dAccumulator>>());
257 data_[VELOCITY] = std::move(velocity);
258 outputMap_["VELOCITY"] = VELOCITY;
259
260 OutputData angularVelocity;
261 angularVelocity.units = "angstroms^2/fs";
262 angularVelocity.title = "AngularVelocity";
263 for (unsigned int i = 0; i < nBins_; i++)
264 angularVelocity.accumulator.push_back(
265 std::make_unique<AccumulatorView<Vector3dAccumulator>>());
266 data_[ANGULARVELOCITY] = std::move(angularVelocity);
267 outputMap_["ANGULARVELOCITY"] = ANGULARVELOCITY;
268
269 OutputData density;
270 density.units = "g cm^-3";
271 density.title = "Density";
272 for (unsigned int i = 0; i < nBins_; i++)
273 density.accumulator.push_back(
274 std::make_unique<AccumulatorView<RealAccumulator>>());
275 data_[DENSITY] = std::move(density);
276 outputMap_["DENSITY"] = DENSITY;
277
278 OutputData activity;
279 activity.units = "unitless";
280 activity.title = "Activity";
281 for (unsigned int i = 0; i < nBins_; i++)
282 activity.accumulator.push_back(
283 std::make_unique<AccumulatorView<StdVectorAccumulator>>());
284 data_[ACTIVITY] = std::move(activity);
285 outputMap_["ACTIVITY"] = ACTIVITY;
286
287 OutputData eField;
288 eField.units = "kcal/mol/angstroms/e";
289 eField.title = "Electric Field";
290 for (unsigned int i = 0; i < nBins_; i++)
291 eField.accumulator.push_back(
292 std::make_unique<AccumulatorView<Vector3dAccumulator>>());
293 data_[ELECTRICFIELD] = std::move(eField);
294 outputMap_["ELECTRICFIELD"] = ELECTRICFIELD;
295
296 OutputData ePot;
297 ePot.units = "kcal/mol/e";
298 ePot.title = "Electrostatic Potential";
299 for (unsigned int i = 0; i < nBins_; i++)
300 ePot.accumulator.push_back(
301 std::make_unique<AccumulatorView<RealAccumulator>>());
302 data_[ELECTROSTATICPOTENTIAL] = std::move(ePot);
303 outputMap_["ELECTROSTATICPOTENTIAL"] = ELECTROSTATICPOTENTIAL;
304
305 if (hasOutputFields) {
306 parseOutputFileFormat(rnemdParams->getOutputFields());
307 } else {
308 if (usePeriodicBoundaryConditions_)
309 outputMask_.set(Z);
310 else
311 outputMask_.set(R);
312 switch (rnemdFluxType_) {
313 case rnemdKE:
314 case rnemdRotKE:
315 case rnemdFullKE:
316 outputMask_.set(TEMPERATURE);
317 break;
318 case rnemdPx:
319 case rnemdPy:
320 outputMask_.set(VELOCITY);
321 break;
322 case rnemdPz:
323 case rnemdPvector:
324 outputMask_.set(VELOCITY);
325 outputMask_.set(DENSITY);
326 break;
327 case rnemdLx:
328 case rnemdLy:
329 case rnemdLz:
330 case rnemdLvector:
331 outputMask_.set(ANGULARVELOCITY);
332 break;
333 case rnemdKeLx:
334 case rnemdKeLy:
335 case rnemdKeLz:
336 case rnemdKeLvector:
337 outputMask_.set(TEMPERATURE);
338 outputMask_.set(ANGULARVELOCITY);
339 break;
340 case rnemdKePx:
341 case rnemdKePy:
342 outputMask_.set(TEMPERATURE);
343 outputMask_.set(VELOCITY);
344 break;
345 case rnemdKePvector:
346 outputMask_.set(TEMPERATURE);
347 outputMask_.set(VELOCITY);
348 outputMask_.set(DENSITY);
349 break;
350 default:
351 break;
352 }
353 }
354
355 if (hasOutputFileName) {
356 rnemdFileName_ = rnemdParams->getOutputFileName();
357 } else {
358 rnemdFileName_ = getPrefix(info->getFinalConfigFileName()) + ".rnemd";
359 }
360
361 // Exchange time should not be less than time step and should be a
362 // multiple of dt
363 exchangeTime_ = rnemdParams->getExchangeTime();
364 RealType dt = simParams->getDt();
365 RealType newET = std::ceil(exchangeTime_ / dt) * dt;
366
367 if (std::fabs(newET - exchangeTime_) > 1e-6) {
368 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
369 "RNEMD: The exchangeTime was reset to %lf,\n"
370 "\t\twhich is a multiple of dt, %lf.\n",
371 newET, dt);
372 painCave.isFatal = 0;
373 painCave.severity = OPENMD_WARNING;
374 simError();
375 exchangeTime_ = newET;
376 }
377
378 currentSnap_ = info->getSnapshotManager()->getCurrentSnapshot();
379 hmat_ = currentSnap_->getHmat();
380
381 // Set up the slab selection logic
382 std::ostringstream selectionAstream;
383 std::ostringstream selectionBstream;
384
385 if (hasSelectionA_) {
386 selectionA_ = rnemdParams->getSelectionA();
387 } else {
388 if (usePeriodicBoundaryConditions_) {
389 slabWidth_ =
390 hasSlabWidth ?
391 rnemdParams->getSlabWidth() :
392 hmat_(rnemdPrivilegedAxis_, rnemdPrivilegedAxis_) / 10.0;
393
394 slabACenter_ = hasSlabACenter ? rnemdParams->getSlabACenter() : 0.0;
395
396 selectionA_ = this->setSelection(slabACenter_);
397 } else {
398 if (hasSphereARadius)
399 sphereARadius_ = rnemdParams->getSphereARadius();
400 else {
401 // use an initial guess to the size of the inner slab to be 1/10 the
402 // radius of an approximately spherical hull:
403 Thermo thermo(info);
404 RealType hVol = thermo.getHullVolume();
405 sphereARadius_ =
406 0.1 * pow((3.0 * hVol / (4.0 * Constants::PI)), 1.0 / 3.0);
407 }
408 selectionAstream << "select r < " << sphereARadius_;
409 selectionA_ = selectionAstream.str();
410 }
411 }
412
413 if (hasSelectionB_) {
414 selectionB_ = rnemdParams->getSelectionB();
415 } else {
416 if (usePeriodicBoundaryConditions_) {
417 slabWidth_ =
418 hasSlabWidth ?
419 rnemdParams->getSlabWidth() :
420 hmat_(rnemdPrivilegedAxis_, rnemdPrivilegedAxis_) / 10.0;
421
422 slabBCenter_ =
423 hasSlabBCenter ?
424 rnemdParams->getSlabBCenter() :
425 hmat_(rnemdPrivilegedAxis_, rnemdPrivilegedAxis_) / 2.0;
426
427 selectionB_ = this->setSelection(slabBCenter_);
428 } else {
429 if (hasSphereBRadius) {
430 sphereBRadius_ = rnemdParams->getSphereBRadius();
431 selectionBstream << "select r > " << sphereBRadius_;
432 selectionB_ = selectionBstream.str();
433 } else {
434 selectionB_ = "select hull";
435 hasSelectionB_ = true;
436 }
437 }
438 }
439
440 // Static object evaluators
441 evaluator_.loadScriptString(rnemdObjectSelection_);
442 if (!evaluator_.isDynamic())
443 seleMan_.setSelectionSet(evaluator_.evaluate());
444
445 evaluatorA_.loadScriptString(selectionA_);
446 if (!evaluatorA_.isDynamic())
447 seleManA_.setSelectionSet(evaluatorA_.evaluate());
448
449 evaluatorB_.loadScriptString(selectionB_);
450 if (!evaluatorB_.isDynamic())
451 seleManB_.setSelectionSet(evaluatorB_.evaluate());
452
453 // Do some sanity checking
454 int selectionCount = seleMan_.getSelectionCount();
455 int nIntegrable = info->getNGlobalIntegrableObjects();
456
457 if (selectionCount > nIntegrable) {
458 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
459 "RNEMD: The current objectSelection,\n"
460 "\t\t%s\n"
461 "\thas resulted in %d selected objects. However,\n"
462 "\tthe total number of integrable objects in the system\n"
463 "\tis only %d. This is almost certainly not what you want\n"
464 "\tto do. A likely cause of this is forgetting the _RB_0\n"
465 "\tselector in the selection script!\n",
466 rnemdObjectSelection_.c_str(), selectionCount, nIntegrable);
467 painCave.isFatal = 0;
468 painCave.severity = OPENMD_WARNING;
469 simError();
470 }
471 }
472
473 RNEMD::~RNEMD() {
474 if (!doRNEMD_) return;
475#ifdef IS_MPI
476 if (worldRank == 0) {
477#endif
478 writeOutputFile();
479
480 rnemdFile_.close();
481#ifdef IS_MPI
482 }
483#endif
484 }
485
486 void RNEMD::getStarted() {
487 if (!doRNEMD_) return;
488 collectData();
489 writeOutputFile();
490 }
491
492 void RNEMD::doRNEMD() {
493 if (!doRNEMD_) return;
494 trialCount_++;
495 hmat_ = currentSnap_->getHmat();
496
497 // dynamic object evaluators:
498 evaluator_.loadScriptString(rnemdObjectSelection_);
499 if (evaluator_.isDynamic()) seleMan_.setSelectionSet(evaluator_.evaluate());
500
501 evaluatorA_.loadScriptString(selectionA_);
502 if (evaluatorA_.isDynamic())
503 seleManA_.setSelectionSet(evaluatorA_.evaluate());
504
505 evaluatorB_.loadScriptString(selectionB_);
506 if (evaluatorB_.isDynamic())
507 seleManB_.setSelectionSet(evaluatorB_.evaluate());
508
509 commonA_ = seleManA_ & seleMan_;
510 commonB_ = seleManB_ & seleMan_;
511
512 // Target exchange quantities (in each exchange) = flux * dividingArea *
513 // dt flux = target flux dividingArea = smallest dividing surface between
514 // the two regions dt = exchange time interval
515
516 RealType area = getDefaultDividingArea();
517
518 kineticTarget_ = kineticFlux_ * exchangeTime_ * area;
519 momentumTarget_ = momentumFluxVector_ * exchangeTime_ * area;
520 angularMomentumTarget_ = angularMomentumFluxVector_ * exchangeTime_ * area;
521 particleTarget_ = particleFlux_ * exchangeTime_ * area;
522
523 if (std::fabs(particleTarget_) > 1.0) {
524 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
525 "RNEMD: The current particleFlux,\n"
526 "\t\t%f\n"
527 "\thas resulted in a target particle exchange of %f.\n"
528 "\tThis is equivalent to moving more than one particle\n"
529 "\tduring each exchange. Please reduce your particleFlux.\n",
530 particleFlux_, particleTarget_);
531 painCave.isFatal = 1;
532 painCave.severity = OPENMD_ERROR;
533 simError();
534 }
535
536 if (rnemdFluxType_ == rnemdParticle || rnemdFluxType_ == rnemdParticleKE) {
537 SelectionManager tempCommonA = seleManA_ & outputSeleMan_;
538 SelectionManager tempCommonB = seleManB_ & outputSeleMan_;
539
540 this->doRNEMDImpl(tempCommonA, tempCommonB);
541 } else {
542 this->doRNEMDImpl(commonA_, commonB_);
543 }
544 }
545
546 void RNEMD::collectData() {
547 if (!doRNEMD_) return;
548 currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
549 hmat_ = currentSnap_->getHmat();
550
551 // collectData() can be called more frequently than the doRNEMD(), so use
552 // the computed area from the last exchange time:
553 RealType area = getDefaultDividingArea();
554 areaAccumulator_.add(area);
555
556 Vector3d u = angularMomentumFluxVector_;
557 u.normalize();
558
559 // throw an error if isDynamic instead?
560 if (outputEvaluator_.isDynamic()) {
561 outputSeleMan_.setSelectionSet(outputEvaluator_.evaluate());
562 }
563
564 int binNo {};
565 int typeIndex(-1);
566 RealType mass {};
567 Vector3d vel {};
568 Vector3d rPos {};
569 RealType KE {};
570 Vector3d L {};
571 Mat3x3d I {};
572 RealType r2 {};
573 Vector3d eField {};
574 int DOF {};
575
576 vector<RealType> binMass(nBins_, 0.0);
577 vector<Vector3d> binP(nBins_, V3Zero);
578 vector<RealType> binOmega(nBins_, 0.0);
579 vector<Vector3d> binL(nBins_, V3Zero);
580 vector<Mat3x3d> binI(nBins_);
581 vector<RealType> binKE(nBins_, 0.0);
582 vector<Vector3d> binEField(nBins_, V3Zero);
583 vector<int> binDOF(nBins_, 0);
584 vector<int> binCount(nBins_, 0);
585 vector<int> binEFieldCount(nBins_, 0);
586 vector<vector<int>> binTypeCounts;
587
588 if (outputMask_[ACTIVITY]) {
589 binTypeCounts.resize(nBins_);
590 for (unsigned int i = 0; i < nBins_; i++) {
591 binTypeCounts[i].resize(outputTypes_.size(), 0);
592 }
593 }
594
595 SimInfo::MoleculeIterator miter;
596 std::vector<StuntDouble*>::iterator iiter;
597 std::vector<AtomType*>::iterator at;
598 Molecule* mol;
599 StuntDouble* sd;
600 AtomType* atype;
601 ConstraintPair* consPair;
602 Molecule::ConstraintPairIterator cpi;
603
604 for (mol = info_->beginMolecule(miter); mol != NULL;
605 mol = info_->nextMolecule(miter)) {
606 for (sd = mol->beginIntegrableObject(iiter); sd != NULL;
607 sd = mol->nextIntegrableObject(iiter)) {
608 if (outputSeleMan_.isSelected(sd)) {
609 Vector3d pos = sd->getPos();
610 binNo = getBin(pos);
611
612 mass = sd->getMass();
613 vel = sd->getVel();
614 rPos = sd->getPos() - coordinateOrigin_;
615 KE = 0.5 * mass * vel.lengthSquare();
616 DOF = 3;
617
618 if (sd->isDirectional()) {
619 Vector3d angMom = sd->getJ();
620 Mat3x3d Ia = sd->getI();
621 if (sd->isLinear()) {
622 int i = sd->linearAxis();
623 int j = (i + 1) % 3;
624 int k = (i + 2) % 3;
625 KE += 0.5 * (angMom[j] * angMom[j] / Ia(j, j) +
626 angMom[k] * angMom[k] / Ia(k, k));
627 DOF += 2;
628 } else {
629 KE += 0.5 * (angMom[0] * angMom[0] / Ia(0, 0) +
630 angMom[1] * angMom[1] / Ia(1, 1) +
631 angMom[2] * angMom[2] / Ia(2, 2));
632 DOF += 3;
633 }
634 }
635
636 L = mass * cross(rPos, vel);
637 I = outProduct(rPos, rPos) * mass;
638 r2 = rPos.lengthSquare();
639 I(0, 0) += mass * r2;
640 I(1, 1) += mass * r2;
641 I(2, 2) += mass * r2;
642
643 if (outputMask_[ACTIVITY]) {
644 typeIndex = -1;
645 if (sd->isAtom()) {
646 atype = static_cast<Atom*>(sd)->getAtomType();
647 at = std::find(outputTypes_.begin(), outputTypes_.end(), atype);
648 if (at != outputTypes_.end()) {
649 typeIndex = std::distance(outputTypes_.begin(), at);
650 }
651 }
652 }
653
654 if (binNo >= 0 && binNo < int(nBins_)) {
655 binCount[binNo]++;
656 binMass[binNo] += mass;
657 binP[binNo] += mass * vel;
658 binKE[binNo] += KE;
659 binI[binNo] += I;
660 binL[binNo] += L;
661 binDOF[binNo] += DOF;
662
663 if (outputMask_[ACTIVITY] && typeIndex != -1)
664 binTypeCounts[binNo][typeIndex]++;
665 }
666 }
667
668 // Calculate the electric field (kcal/mol/e/Angstrom) for all atoms in
669 // the box
670 if (outputMask_[ELECTRICFIELD]) {
671 int atomBinNo;
672 if (sd->isRigidBody()) {
673 RigidBody* rb = static_cast<RigidBody*>(sd);
674 std::vector<Atom*>::iterator ai;
675 Atom* atom;
676 for (atom = rb->beginAtom(ai); atom != NULL;
677 atom = rb->nextAtom(ai)) {
678 atomBinNo = getBin(atom->getPos());
679 eField = atom->getElectricField();
680
681 binEFieldCount[atomBinNo]++;
682 binEField[atomBinNo] += eField;
683 }
684 } else {
685 eField = sd->getElectricField();
686 atomBinNo = getBin(sd->getPos());
687
688 binEFieldCount[atomBinNo]++;
689 binEField[atomBinNo] += eField;
690 }
691 }
692 }
693
694 // we need to subtract out degrees of freedom from constraints
695 // belonging in this bin:
696 if (outputSeleMan_.isSelected(mol)) {
697 for (consPair = mol->beginConstraintPair(cpi); consPair != NULL;
698 consPair = mol->nextConstraintPair(cpi)) {
699 Vector3d posA = consPair->getConsElem1()->getPos();
700 Vector3d posB = consPair->getConsElem2()->getPos();
701
702 if (usePeriodicBoundaryConditions_) {
703 currentSnap_->wrapVector(posA);
704 currentSnap_->wrapVector(posB);
705 }
706
707 Vector3d coc = 0.5 * (posA + posB);
708 int binCons = getBin(coc);
709 binDOF[binCons] -= 1;
710 }
711 }
712 }
713
714#ifdef IS_MPI
715 for (unsigned int i = 0; i < nBins_; i++) {
716 MPI_Allreduce(MPI_IN_PLACE, &binCount[i], 1, MPI_INT, MPI_SUM,
717 MPI_COMM_WORLD);
718 MPI_Allreduce(MPI_IN_PLACE, &binMass[i], 1, MPI_REALTYPE, MPI_SUM,
719 MPI_COMM_WORLD);
720 MPI_Allreduce(MPI_IN_PLACE, binP[i].getArrayPointer(), 3, MPI_REALTYPE,
721 MPI_SUM, MPI_COMM_WORLD);
722 MPI_Allreduce(MPI_IN_PLACE, binL[i].getArrayPointer(), 3, MPI_REALTYPE,
723 MPI_SUM, MPI_COMM_WORLD);
724 MPI_Allreduce(MPI_IN_PLACE, binI[i].getArrayPointer(), 9, MPI_REALTYPE,
725 MPI_SUM, MPI_COMM_WORLD);
726 MPI_Allreduce(MPI_IN_PLACE, &binKE[i], 1, MPI_REALTYPE, MPI_SUM,
727 MPI_COMM_WORLD);
728 MPI_Allreduce(MPI_IN_PLACE, &binDOF[i], 1, MPI_INT, MPI_SUM,
729 MPI_COMM_WORLD);
730 MPI_Allreduce(MPI_IN_PLACE, &binEFieldCount[i], 1, MPI_INT, MPI_SUM,
731 MPI_COMM_WORLD);
732
733 if (outputMask_[ELECTRICFIELD]) {
734 MPI_Allreduce(MPI_IN_PLACE, binEField[i].getArrayPointer(), 3,
735 MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
736 }
737 if (outputMask_[ACTIVITY]) {
738 MPI_Allreduce(MPI_IN_PLACE, &binTypeCounts[i][0], outputTypes_.size(),
739 MPI_INT, MPI_SUM, MPI_COMM_WORLD);
740 }
741 }
742#endif
743
744 Vector3d omega;
745 RealType z, r, temp, binVolume, den(0.0), dz(0.0);
746 std::vector<RealType> nden(outputTypes_.size(), 0.0);
747 RealType boxVolume = currentSnap_->getVolume();
748 RealType ePot(0.0);
749
750 for (unsigned int i = 0; i < nBins_; i++) {
751 if (usePeriodicBoundaryConditions_) {
752 z = (((RealType)i + 0.5) / (RealType)nBins_) *
753 hmat_(rnemdPrivilegedAxis_, rnemdPrivilegedAxis_);
754 data_[Z].accumulator[i]->add(z);
755
756 binVolume = boxVolume / nBins_;
757 dz = hmat_(rnemdPrivilegedAxis_, rnemdPrivilegedAxis_) /
758 (RealType)nBins_;
759 } else {
760 r = (((RealType)i + 0.5) * binWidth_);
761 data_[R].accumulator[i]->add(r);
762
763 RealType rinner = (RealType)i * binWidth_;
764 RealType router = (RealType)(i + 1) * binWidth_;
765 binVolume =
766 (4.0 * Constants::PI * (pow(router, 3) - pow(rinner, 3))) / 3.0;
767 }
768
769 // The calculations of the following properties are done regardless
770 // of whether or not the selected species are present in the bin
771 if (outputMask_[ELECTRICFIELD] && binEFieldCount[i] > 0) {
772 eField = binEField[i] / RealType(binEFieldCount[i]);
773 data_[ELECTRICFIELD].accumulator[i]->add(eField);
774 }
775
776 if (outputMask_[ELECTROSTATICPOTENTIAL]) {
777 if (usePeriodicBoundaryConditions_ && binEFieldCount[i] > 0) {
778 ePot += eField[rnemdPrivilegedAxis_] * dz;
779 data_[ELECTROSTATICPOTENTIAL].accumulator[i]->add(ePot);
780 }
781 }
782
783 // For the following properties, zero should be added if the selected
784 // species is not present in the bin
785 if (outputMask_[DENSITY]) {
786 den = binMass[i] * Constants::densityConvert / binVolume;
787 data_[DENSITY].accumulator[i]->add(den);
788 }
789
790 if (outputMask_[ACTIVITY]) {
791 for (unsigned int j = 0; j < outputTypes_.size(); j++) {
792 nden[j] = (binTypeCounts[i][j] / binVolume) *
793 Constants::concentrationConvert;
794 }
795 data_[ACTIVITY].accumulator[i]->add(nden);
796 }
797
798 if (binCount[i] > 0) {
799 // The calculations of the following properties are meaningless if
800 // the selected species is not found in the bin
801 if (outputMask_[VELOCITY]) {
802 vel = binP[i] / binMass[i];
803 data_[VELOCITY].accumulator[i]->add(vel);
804 }
805
806 if (outputMask_[ANGULARVELOCITY]) {
807 omega = binI[i].inverse() * binL[i];
808 data_[ANGULARVELOCITY].accumulator[i]->add(omega);
809 }
810
811 if (outputMask_[TEMPERATURE]) {
812 if (binDOF[i] > 0) {
813 temp = 2.0 * binKE[i] /
814 (binDOF[i] * Constants::kb * Constants::energyConvert);
815 data_[TEMPERATURE].accumulator[i]->add(temp);
816 } else {
817 std::cerr << "No degrees of freedom in this bin?\n";
818 }
819 }
820 }
821 }
822
823 hasData_ = true;
824 }
825
826 void RNEMD::writeOutputFile() {
827 if (!doRNEMD_) return;
828 if (!hasData_) return;
829
830#ifdef IS_MPI
831 // If we're the primary node, should we print out the results
832 int worldRank;
833 MPI_Comm_rank(MPI_COMM_WORLD, &worldRank);
834
835 if (worldRank == 0) {
836#endif
837 rnemdFile_.open(rnemdFileName_.c_str(), std::ios::out | std::ios::trunc);
838
839 if (!rnemdFile_) {
840 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
841 "Could not open \"%s\" for RNEMD output.\n",
842 rnemdFileName_.c_str());
843 painCave.isFatal = 1;
844 simError();
845 }
846
847 RealType time = currentSnap_->getTime();
848 RealType avgArea = areaAccumulator_.getAverage();
849
850 RealType Jz(0.0);
851 Vector3d JzP(V3Zero);
852 Vector3d JzL(V3Zero);
853 RealType Jpart(0.0);
854
855 if (time >= info_->getSimParams()->getDt()) {
856 Jz = kineticExchange_ / (time * avgArea) / Constants::energyConvert;
857 JzP = momentumExchange_ / (time * avgArea);
858 JzL = angularMomentumExchange_ / (time * avgArea);
859 Jpart = particleExchange_ / (time * avgArea);
860 }
861
862 rnemdFile_ << "#######################################################\n";
863 rnemdFile_ << "# RNEMD {\n";
864 rnemdFile_ << "# exchangeMethod = \"" << rnemdMethodLabel_ << "\";\n";
865 rnemdFile_ << "# fluxType = \"" << rnemdFluxTypeLabel_ << "\";\n";
866
867 if (usePeriodicBoundaryConditions_)
868 rnemdFile_ << "# privilegedAxis = " << rnemdAxisLabel_ << ";\n";
869
870 rnemdFile_ << "# exchangeTime = " << exchangeTime_ << ";\n";
871 rnemdFile_ << "# objectSelection = \"" << rnemdObjectSelection_
872 << "\";\n";
873 rnemdFile_ << "# selectionA = \"" << selectionA_ << "\";\n";
874 rnemdFile_ << "# selectionB = \"" << selectionB_ << "\";\n";
875 rnemdFile_ << "# outputSelection = \"" << outputSelection_ << "\";\n";
876 rnemdFile_ << "# }\n";
877 rnemdFile_ << "#######################################################\n";
878 rnemdFile_ << "# RNEMD report:\n";
879 rnemdFile_ << "# running time = " << time << " fs\n";
880 rnemdFile_ << "# Target flux:\n";
881 rnemdFile_ << "# kinetic = "
882 << kineticFlux_ / Constants::energyConvert
883 << " (kcal/mol/A^2/fs)\n";
884 rnemdFile_ << "# momentum = " << momentumFluxVector_
885 << " (amu/A/fs^2)\n";
886 rnemdFile_ << "# angular momentum = " << angularMomentumFluxVector_
887 << " (amu/A^2/fs^2)\n";
888 rnemdFile_ << "# particle = " << particleFlux_
889 << " (particles/A^2/fs)\n";
890
891 rnemdFile_ << "# Target one-time exchanges:\n";
892 rnemdFile_ << "# kinetic = "
893 << kineticTarget_ / Constants::energyConvert
894 << " (kcal/mol)\n";
895 rnemdFile_ << "# momentum = " << momentumTarget_
896 << " (amu*A/fs)\n";
897 rnemdFile_ << "# angular momentum = " << angularMomentumTarget_
898 << " (amu*A^2/fs)\n";
899 rnemdFile_ << "# particle = " << particleTarget_
900 << " (particles)\n";
901 rnemdFile_ << "# Actual exchange totals:\n";
902 rnemdFile_ << "# kinetic = "
903 << kineticExchange_ / Constants::energyConvert
904 << " (kcal/mol)\n";
905 rnemdFile_ << "# momentum = " << momentumExchange_
906 << " (amu*A/fs)\n";
907 rnemdFile_ << "# angular momentum = " << angularMomentumExchange_
908 << " (amu*A^2/fs)\n";
909 rnemdFile_ << "# particles = " << particleExchange_
910 << " (particles)\n";
911
912 rnemdFile_ << "# Actual flux:\n";
913 rnemdFile_ << "# kinetic = " << Jz << " (kcal/mol/A^2/fs)\n";
914 rnemdFile_ << "# momentum = " << JzP << " (amu/A/fs^2)\n";
915 rnemdFile_ << "# angular momentum = " << JzL << " (amu/A^2/fs^2)\n";
916 rnemdFile_ << "# particle = " << Jpart
917 << " (particles/A^2/fs)\n";
918 rnemdFile_ << "# Exchange statistics:\n";
919 rnemdFile_ << "# attempted = " << trialCount_ << "\n";
920 rnemdFile_ << "# failed = " << failTrialCount_ << "\n";
921 if (rnemdMethodLabel_ == "NIVS") {
922 rnemdFile_ << "# NIVS root-check errors = " << failRootCount_ << "\n";
923 }
924 rnemdFile_ << "#######################################################\n";
925
926 // write title
927 rnemdFile_ << "#";
928 for (unsigned int i = 0; i < outputMask_.size(); ++i) {
929 if (outputMask_[i]) {
930 rnemdFile_ << "\t" << data_[i].title << "(" << data_[i].units << ")";
931
932 // add some extra tabs for column alignment
933 if (data_[i].accumulator[0]->getType() ==
934 std::type_index(typeid(Vector3d))) {
935 rnemdFile_ << "\t\t";
936 }
937
938 if (data_[i].accumulator[0]->getType() ==
939 std::type_index(typeid(std::vector<RealType>))) {
940 rnemdFile_ << "(";
941 for (unsigned int type = 0; type < outputTypes_.size(); type++) {
942 rnemdFile_ << outputTypes_[type]->getName() << "\t";
943 }
944 rnemdFile_ << ")\t";
945 }
946 }
947 }
948
949 rnemdFile_ << std::endl;
950
951 std::vector<int> nonEmptyAccumulators(nBins_);
952 int numberOfAccumulators {};
953
954 for (unsigned int i = 0; i < outputMask_.size(); ++i) {
955 if (outputMask_[i]) {
956 for (unsigned int bin = 0; bin < nBins_; bin++) {
957 nonEmptyAccumulators[bin] +=
958 static_cast<int>(data_[i].accumulator[bin]->getCount() != 0);
959 }
960
961 numberOfAccumulators++;
962 }
963 }
964
965 rnemdFile_.precision(8);
966
967 for (unsigned int bin = 0; bin < nBins_; bin++) {
968 if (nonEmptyAccumulators[bin] == numberOfAccumulators) {
969 for (unsigned int i = 0; i < outputMask_.size(); ++i) {
970 if (outputMask_[i]) {
971 std::string message =
972 "RNEMD detected a numerical error writing: " +
973 data_[i].title + " for bin " + std::to_string(bin);
974
975 data_[i].accumulator[bin]->writeData(rnemdFile_, message);
976 }
977 }
978
979 rnemdFile_ << std::endl;
980 }
981 }
982
983 rnemdFile_ << "#######################################################\n";
984 rnemdFile_ << "# 95% confidence intervals in those quantities follow:\n";
985 rnemdFile_ << "#######################################################\n";
986
987 for (unsigned int bin = 0; bin < nBins_; bin++) {
988 if (nonEmptyAccumulators[bin] == numberOfAccumulators) {
989 rnemdFile_ << "#";
990
991 for (unsigned int i = 0; i < outputMask_.size(); ++i) {
992 if (outputMask_[i]) {
993 std::string message =
994 "RNEMD detected a numerical error writing: " +
995 data_[i].title + " std. dev. for bin " + std::to_string(bin);
996
997 data_[i].accumulator[bin]->writeErrorBars(rnemdFile_, message);
998 }
999 }
1000
1001 rnemdFile_ << std::endl;
1002 }
1003 }
1004
1005 rnemdFile_.flush();
1006 rnemdFile_.close();
1007#ifdef IS_MPI
1008 }
1009#endif
1010 }
1011
1012 void RNEMD::setKineticFlux(RealType kineticFlux) {
1013 // convert the kcal / mol / Angstroms^2 / fs values in the md file
1014 // into amu / fs^3:
1015 kineticFlux_ = kineticFlux * Constants::energyConvert;
1016 }
1017
1018 void RNEMD::setParticleFlux(RealType particleFlux) {
1019 particleFlux_ = particleFlux;
1020 }
1021
1022 void RNEMD::setMomentumFluxVector(
1023 const std::vector<RealType>& momentumFluxVector) {
1024 if (momentumFluxVector.size() != 3) {
1025 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
1026 "RNEMD: Incorrect number of parameters specified for "
1027 "momentumFluxVector.\n"
1028 "\tthere should be 3 parameters, but %lu were specified.\n",
1029 momentumFluxVector.size());
1030 painCave.isFatal = 1;
1031 simError();
1032 }
1033
1034 momentumFluxVector_.x() = momentumFluxVector[0];
1035 momentumFluxVector_.y() = momentumFluxVector[1];
1036 momentumFluxVector_.z() = momentumFluxVector[2];
1037 }
1038
1039 void RNEMD::setAngularMomentumFluxVector(
1040 const std::vector<RealType>& angularMomentumFluxVector) {
1041 if (angularMomentumFluxVector.size() != 3) {
1042 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
1043 "RNEMD: Incorrect number of parameters specified for "
1044 "angularMomentumFluxVector.\n"
1045 "\tthere should be 3 parameters, but %lu were specified.\n",
1046 angularMomentumFluxVector.size());
1047 painCave.isFatal = 1;
1048 simError();
1049 }
1050
1051 angularMomentumFluxVector_.x() = angularMomentumFluxVector[0];
1052 angularMomentumFluxVector_.y() = angularMomentumFluxVector[1];
1053 angularMomentumFluxVector_.z() = angularMomentumFluxVector[2];
1054 }
1055
1056 void RNEMD::parseOutputFileFormat(const std::string& format) {
1057 if (!doRNEMD_) return;
1058 StringTokenizer tokenizer(format, " ,;|\t\n\r");
1059
1060 while (tokenizer.hasMoreTokens()) {
1061 std::string token(tokenizer.nextToken());
1062 toUpper(token);
1063 OutputMapType::iterator i = outputMap_.find(token);
1064 if (i != outputMap_.end()) {
1065 outputMask_.set(i->second);
1066 } else {
1067 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
1068 "RNEMD::parseOutputFileFormat: %s is not a recognized\n"
1069 "\toutputFileFormat keyword.\n",
1070 token.c_str());
1071 painCave.isFatal = 0;
1072 painCave.severity = OPENMD_ERROR;
1073 simError();
1074 }
1075 }
1076 }
1077
1078 std::string RNEMD::setSelection(RealType& slabCenter) {
1079 bool printSlabCenterWarning {false};
1080
1081 Vector3d tempSlabCenter {V3Zero};
1082 tempSlabCenter[rnemdPrivilegedAxis_] = slabCenter;
1083
1084 RealType hmat_2 = hmat_(rnemdPrivilegedAxis_, rnemdPrivilegedAxis_) / 2.0;
1085
1086 if (slabCenter > hmat_2) {
1087 currentSnap_->wrapVector(tempSlabCenter);
1088 printSlabCenterWarning = true;
1089 } else if (slabCenter < -hmat_2) {
1090 currentSnap_->wrapVector(tempSlabCenter);
1091 printSlabCenterWarning = true;
1092 }
1093
1094 if (printSlabCenterWarning) {
1095 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
1096 "The given slab center was set to %0.2f. In the wrapped "
1097 "coordinates\n"
1098 "\t[-Hmat/2, +Hmat/2], this has been remapped to %0.2f.\n",
1099 slabCenter, tempSlabCenter[rnemdPrivilegedAxis_]);
1100 painCave.isFatal = 0;
1101 painCave.severity = OPENMD_WARNING;
1102 simError();
1103
1104 slabCenter = tempSlabCenter[rnemdPrivilegedAxis_];
1105 }
1106
1107 Vector3d leftSlab {V3Zero};
1108 const RealType& leftSlabBoundary = leftSlab[rnemdPrivilegedAxis_];
1109 leftSlab[rnemdPrivilegedAxis_] = slabCenter - 0.5 * slabWidth_;
1110 currentSnap_->wrapVector(leftSlab);
1111
1112 Vector3d rightSlab {V3Zero};
1113 const RealType& rightSlabBoundary = rightSlab[rnemdPrivilegedAxis_];
1114 rightSlab[rnemdPrivilegedAxis_] = slabCenter + 0.5 * slabWidth_;
1115 currentSnap_->wrapVector(rightSlab);
1116
1117 std::ostringstream selectionStream;
1118
1119 selectionStream << "select wrapped" << rnemdAxisLabel_
1120 << " >= " << leftSlabBoundary;
1121
1122 if (leftSlabBoundary > rightSlabBoundary)
1123 selectionStream << " || wrapped" << rnemdAxisLabel_ << " < "
1124 << rightSlabBoundary;
1125 else
1126 selectionStream << " && wrapped" << rnemdAxisLabel_ << " < "
1127 << rightSlabBoundary;
1128
1129 return selectionStream.str();
1130 }
1131
1132 RealType RNEMD::getDefaultDividingArea() {
1133 if (hasDividingArea_) return dividingArea_;
1134
1135 Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot();
1136
1137 if (hasSelectionA_) {
1138 if (evaluatorA_.hasSurfaceArea()) {
1139 areaA_ = evaluatorA_.getSurfaceArea();
1140 volumeA_ = evaluatorA_.getVolume();
1141 } else {
1142 int isd;
1143 StuntDouble* sd;
1144 std::vector<StuntDouble*> aSites;
1145 seleManA_.setSelectionSet(evaluatorA_.evaluate());
1146 for (sd = seleManA_.beginSelected(isd); sd != NULL;
1147 sd = seleManA_.nextSelected(isd)) {
1148 aSites.push_back(sd);
1149 }
1150#if defined(HAVE_QHULL)
1151 ConvexHull* surfaceMeshA = new ConvexHull();
1152 surfaceMeshA->computeHull(aSites);
1153 areaA_ = surfaceMeshA->getArea();
1154 volumeA_ = surfaceMeshA->getVolume();
1155 delete surfaceMeshA;
1156#else
1157 snprintf(
1158 painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
1159 "RNEMD::getDividingArea : Hull calculation is not possible\n"
1160 "\twithout libqhull. Please rebuild OpenMD with qhull enabled.");
1161 painCave.severity = OPENMD_ERROR;
1162 painCave.isFatal = 1;
1163 simError();
1164#endif
1165 }
1166 } else {
1167 if (usePeriodicBoundaryConditions_) {
1168 // in periodic boundaries, the surface area is twice the
1169 // area of the current box, normal to the privileged axis:
1170 switch (rnemdPrivilegedAxis_) {
1171 case rnemdX:
1172 areaA_ = 2.0 * snap->getYZarea();
1173 break;
1174 case rnemdY:
1175 areaA_ = 2.0 * snap->getXZarea();
1176 break;
1177 case rnemdZ:
1178 default:
1179 areaA_ = 2.0 * snap->getXYarea();
1180 }
1181
1182 volumeA_ = areaA_ * slabWidth_;
1183 } else {
1184 // in non-periodic simulations, without explicitly setting
1185 // selections, the sphere radius sets the surface area of the
1186 // dividing surface:
1187 areaA_ = 4.0 * Constants::PI * std::pow(sphereARadius_, 2);
1188 volumeA_ = 4.0 * Constants::PI * std::pow(sphereARadius_, 3) / 3.0;
1189 }
1190 }
1191
1192 if (hasSelectionB_) {
1193 if (evaluatorB_.hasSurfaceArea()) {
1194 areaB_ = evaluatorB_.getSurfaceArea();
1195 volumeB_ = evaluatorB_.getVolume();
1196 } else {
1197 int isd;
1198 StuntDouble* sd;
1199 std::vector<StuntDouble*> bSites;
1200 seleManB_.setSelectionSet(evaluatorB_.evaluate());
1201 for (sd = seleManB_.beginSelected(isd); sd != NULL;
1202 sd = seleManB_.nextSelected(isd)) {
1203 bSites.push_back(sd);
1204 }
1205
1206#if defined(HAVE_QHULL)
1207 ConvexHull* surfaceMeshB = new ConvexHull();
1208 surfaceMeshB->computeHull(bSites);
1209 areaB_ = surfaceMeshB->getArea();
1210 volumeB_ = surfaceMeshB->getVolume();
1211 delete surfaceMeshB;
1212#else
1213 snprintf(
1214 painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
1215 "RNEMD::getDividingArea : Hull calculation is not possible\n"
1216 "\twithout libqhull. Please rebuild OpenMD with qhull enabled.");
1217 painCave.severity = OPENMD_ERROR;
1218 painCave.isFatal = 1;
1219 simError();
1220#endif
1221 }
1222 } else {
1223 if (usePeriodicBoundaryConditions_) {
1224 // in periodic boundaries, the surface area is twice the
1225 // area of the current box, normal to the privileged axis:
1226 switch (rnemdPrivilegedAxis_) {
1227 case rnemdX:
1228 areaB_ = 2.0 * snap->getYZarea();
1229 break;
1230 case rnemdY:
1231 areaB_ = 2.0 * snap->getXZarea();
1232 break;
1233 case rnemdZ:
1234 default:
1235 areaB_ = 2.0 * snap->getXYarea();
1236 }
1237
1238 volumeB_ = areaB_ * slabWidth_;
1239 } else {
1240 // in non-periodic simulations, without explicitly setting
1241 // selections, but if a sphereBradius has been set, just use that:
1242 areaB_ = 4.0 * Constants::PI * pow(sphereBRadius_, 2);
1243 Thermo thermo(info_);
1244 RealType hVol = thermo.getHullVolume();
1245 volumeB_ = hVol - 4.0 * Constants::PI * pow(sphereBRadius_, 3) / 3.0;
1246 }
1247 }
1248
1249 dividingArea_ = min(areaA_, areaB_);
1250 hasDividingArea_ = true;
1251 return dividingArea_;
1252 }
1253
1254 int RNEMD::getBin(Vector3d pos) {
1255 if (usePeriodicBoundaryConditions_) {
1256 currentSnap_->wrapVector(pos);
1257
1258 return int(nBins_ *
1259 (pos[rnemdPrivilegedAxis_] /
1260 hmat_(rnemdPrivilegedAxis_, rnemdPrivilegedAxis_) +
1261 0.5)) %
1262 nBins_;
1263 } else {
1264 Vector3d rPos = pos - coordinateOrigin_;
1265 return int(rPos.length() / binWidth_);
1266 }
1267 }
1268} // namespace OpenMD::RNEMD
OpenMD::AtomType
AtomType is what OpenMD looks to for unchanging data about an atom.
Definition AtomType.hpp:66
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::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::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::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