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