OpenMD 3.0
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
NIVS.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/NIVS.hpp"
46
47#include <algorithm>
48#include <cmath>
49#include <map>
50#include <set>
51#include <sstream>
52#include <string>
53#include <vector>
54
55#ifdef IS_MPI
56#include <mpi.h>
57#endif
58
59#include "brains/ForceManager.hpp"
60#include "brains/Thermo.hpp"
61#include "io/Globals.hpp"
62#include "math/ConvexHull.hpp"
63#include "math/Polynomial.hpp"
64#include "math/SquareMatrix3.hpp"
65#include "math/Vector.hpp"
66#include "math/Vector3.hpp"
67#include "primitives/Molecule.hpp"
68#include "primitives/StuntDouble.hpp"
69#include "rnemd/RNEMD.hpp"
70#include "rnemd/RNEMDParameters.hpp"
71#include "types/FixedChargeAdapter.hpp"
72#include "types/FluctuatingChargeAdapter.hpp"
73#include "utils/Constants.hpp"
74
75#define HONKING_LARGE_VALUE 1.0e10
76
77namespace OpenMD::RNEMD {
78
79 NIVSMethod::NIVSMethod(SimInfo* info, ForceManager* forceMan) :
80 RNEMD {info, forceMan} {
81 rnemdMethodLabel_ = "NIVS";
82
83 RNEMDParameters* rnemdParams = info->getSimParams()->getRNEMDParameters();
84
85 bool hasKineticFlux = rnemdParams->haveKineticFlux();
86 bool hasMomentumFlux = rnemdParams->haveMomentumFlux();
87
88 bool methodFluxMismatch = false;
89 bool hasCorrectFlux = false;
90
91 switch (rnemdFluxType_) {
92 case rnemdKE:
93 case rnemdRotKE:
94 case rnemdFullKE:
95 hasCorrectFlux = hasKineticFlux;
96 break;
97 case rnemdPx:
98 case rnemdPy:
99 case rnemdPz:
100 hasCorrectFlux = hasMomentumFlux;
101 break;
102 case rnemdKePx:
103 case rnemdKePy:
104 hasCorrectFlux = hasMomentumFlux && hasKineticFlux;
105 break;
106 default:
107 methodFluxMismatch = true;
108 break;
109 }
110
111 if (methodFluxMismatch) {
112 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
113 "RNEMD: The current method,\n"
114 "\t\tNIVS\n"
115 "\tcannot be used with the current flux type, %s\n",
116 rnemdFluxTypeLabel_.c_str());
117 painCave.isFatal = 1;
118 painCave.severity = OPENMD_ERROR;
119 simError();
120 }
121
122 if (!hasCorrectFlux) {
123 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
124 "RNEMD: The current method, NIVS, and flux type, %s,\n"
125 "\tdid not have the correct flux value specified. Options\n"
126 "\tinclude: kineticFlux and momentumFlux.\n",
127 rnemdFluxTypeLabel_.c_str());
128 painCave.isFatal = 1;
129 painCave.severity = OPENMD_ERROR;
130 simError();
131 }
132
133 if (hasKineticFlux) {
134 setKineticFlux(rnemdParams->getKineticFlux());
135 } else {
136 setKineticFlux(0.0);
137 }
138
139 if (hasMomentumFlux) {
140 RealType momentumFlux = rnemdParams->getMomentumFlux();
141 std::vector<RealType> momentumFluxVector(3);
142
143 switch (rnemdFluxType_) {
144 case rnemdPx:
145 momentumFluxVector[0] = momentumFlux;
146 break;
147 case rnemdPy:
148 momentumFluxVector[1] = momentumFlux;
149 break;
150 case rnemdPz:
151 momentumFluxVector[2] = momentumFlux;
152 break;
153 default:
154 break;
155 }
156
157 setMomentumFluxVector(momentumFluxVector);
158 }
159 }
160
161 void NIVSMethod::doRNEMDImpl(SelectionManager& smanA,
162 SelectionManager& smanB) {
163 if (!doRNEMD_) return;
164 int selei;
165 int selej;
166
167 StuntDouble* sd;
168
169 std::vector<StuntDouble*> hotBin, coldBin;
170
171 RealType Phx = 0.0;
172 RealType Phy = 0.0;
173 RealType Phz = 0.0;
174 RealType Khx = 0.0;
175 RealType Khy = 0.0;
176 RealType Khz = 0.0;
177 RealType Khw = 0.0;
178 RealType Pcx = 0.0;
179 RealType Pcy = 0.0;
180 RealType Pcz = 0.0;
181 RealType Kcx = 0.0;
182 RealType Kcy = 0.0;
183 RealType Kcz = 0.0;
184 RealType Kcw = 0.0;
185
186 for (sd = smanA.beginSelected(selei); sd != NULL;
187 sd = smanA.nextSelected(selei)) {
188 Vector3d pos = sd->getPos();
189
190 // wrap the stuntdouble's position back into the box:
191
192 if (usePeriodicBoundaryConditions_) currentSnap_->wrapVector(pos);
193
194 RealType mass = sd->getMass();
195 Vector3d vel = sd->getVel();
196
197 hotBin.push_back(sd);
198 Phx += mass * vel.x();
199 Phy += mass * vel.y();
200 Phz += mass * vel.z();
201 Khx += mass * vel.x() * vel.x();
202 Khy += mass * vel.y() * vel.y();
203 Khz += mass * vel.z() * vel.z();
204 if (sd->isDirectional()) {
205 Vector3d angMom = sd->getJ();
206 Mat3x3d I = sd->getI();
207 if (sd->isLinear()) {
208 int i = sd->linearAxis();
209 int j = (i + 1) % 3;
210 int k = (i + 2) % 3;
211 Khw +=
212 angMom[j] * angMom[j] / I(j, j) + angMom[k] * angMom[k] / I(k, k);
213 } else {
214 Khw += angMom[0] * angMom[0] / I(0, 0) +
215 angMom[1] * angMom[1] / I(1, 1) +
216 angMom[2] * angMom[2] / I(2, 2);
217 }
218 }
219 }
220 for (sd = smanB.beginSelected(selej); sd != NULL;
221 sd = smanB.nextSelected(selej)) {
222 Vector3d pos = sd->getPos();
223
224 // wrap the stuntdouble's position back into the box:
225
226 if (usePeriodicBoundaryConditions_) currentSnap_->wrapVector(pos);
227
228 RealType mass = sd->getMass();
229 Vector3d vel = sd->getVel();
230
231 coldBin.push_back(sd);
232 Pcx += mass * vel.x();
233 Pcy += mass * vel.y();
234 Pcz += mass * vel.z();
235 Kcx += mass * vel.x() * vel.x();
236 Kcy += mass * vel.y() * vel.y();
237 Kcz += mass * vel.z() * vel.z();
238 if (sd->isDirectional()) {
239 Vector3d angMom = sd->getJ();
240 Mat3x3d I = sd->getI();
241 if (sd->isLinear()) {
242 int i = sd->linearAxis();
243 int j = (i + 1) % 3;
244 int k = (i + 2) % 3;
245 Kcw +=
246 angMom[j] * angMom[j] / I(j, j) + angMom[k] * angMom[k] / I(k, k);
247 } else {
248 Kcw += angMom[0] * angMom[0] / I(0, 0) +
249 angMom[1] * angMom[1] / I(1, 1) +
250 angMom[2] * angMom[2] / I(2, 2);
251 }
252 }
253 }
254
255 Khx *= 0.5;
256 Khy *= 0.5;
257 Khz *= 0.5;
258 Khw *= 0.5;
259 Kcx *= 0.5;
260 Kcy *= 0.5;
261 Kcz *= 0.5;
262 Kcw *= 0.5;
263
264#ifdef IS_MPI
265 MPI_Allreduce(MPI_IN_PLACE, &Phx, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
266 MPI_Allreduce(MPI_IN_PLACE, &Phy, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
267 MPI_Allreduce(MPI_IN_PLACE, &Phz, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
268 MPI_Allreduce(MPI_IN_PLACE, &Pcx, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
269 MPI_Allreduce(MPI_IN_PLACE, &Pcy, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
270 MPI_Allreduce(MPI_IN_PLACE, &Pcz, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
271
272 MPI_Allreduce(MPI_IN_PLACE, &Khx, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
273 MPI_Allreduce(MPI_IN_PLACE, &Khy, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
274 MPI_Allreduce(MPI_IN_PLACE, &Khz, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
275 MPI_Allreduce(MPI_IN_PLACE, &Khw, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
276
277 MPI_Allreduce(MPI_IN_PLACE, &Kcx, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
278 MPI_Allreduce(MPI_IN_PLACE, &Kcy, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
279 MPI_Allreduce(MPI_IN_PLACE, &Kcz, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
280 MPI_Allreduce(MPI_IN_PLACE, &Kcw, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
281#endif
282
283 // solve coldBin coeff's first
284 RealType px = Pcx / Phx;
285 RealType py = Pcy / Phy;
286 RealType pz = Pcz / Phz;
287 RealType c(0.0), x(0.0), y(0.0), z(0.0);
288 bool successfulScale = false;
289 if ((rnemdFluxType_ == rnemdFullKE) || (rnemdFluxType_ == rnemdRotKE)) {
290 // may need sanity check Khw & Kcw > 0
291
292 if (rnemdFluxType_ == rnemdFullKE) {
293 c = 1.0 - kineticTarget_ / (Kcx + Kcy + Kcz + Kcw);
294 } else {
295 c = 1.0 - kineticTarget_ / Kcw;
296 }
297
298 if ((c > 0.81) && (c < 1.21)) { // restrict scaling coefficients
299 c = sqrt(c);
300
301 RealType w = 0.0;
302 if (rnemdFluxType_ == rnemdFullKE) {
303 x = 1.0 + px * (1.0 - c);
304 y = 1.0 + py * (1.0 - c);
305 z = 1.0 + pz * (1.0 - c);
306 /* more complicated way
307 w = 1.0 + (Kcw - Kcw * c * c - (c * c * (Kcx + Kcy + Kcz
308 + Khx * px * px + Khy * py * py + Khz * pz * pz)
309 - 2.0 * c * (Khx * px * (1.0 + px) + Khy * py * (1.0 + py)
310 + Khz * pz * (1.0 + pz)) + Khx * px * (2.0 + px)
311 + Khy * py * (2.0 + py) + Khz * pz * (2.0 + pz)
312 - Kcx - Kcy - Kcz)) / Khw; the following is simpler
313 */
314 if ((std::fabs(x - 1.0) < 0.1) && (std::fabs(y - 1.0) < 0.1) &&
315 (std::fabs(z - 1.0) < 0.1)) {
316 w = 1.0 + (kineticTarget_ + Khx * (1.0 - x * x) +
317 Khy * (1.0 - y * y) + Khz * (1.0 - z * z)) /
318 Khw;
319 } // no need to calculate w if x, y or z is out of range
320 } else {
321 w = 1.0 + kineticTarget_ / Khw;
322 }
323 if ((w > 0.81) && (w < 1.21)) { // restrict scaling coefficients
324 // if w is in the right range, so should be x, y, z.
325 std::vector<StuntDouble*>::iterator sdi;
326 Vector3d vel;
327 for (sdi = coldBin.begin(); sdi != coldBin.end(); ++sdi) {
328 if (rnemdFluxType_ == rnemdFullKE) {
329 vel = (*sdi)->getVel() * c;
330 (*sdi)->setVel(vel);
331 }
332 if ((*sdi)->isDirectional()) {
333 Vector3d angMom = (*sdi)->getJ() * c;
334 (*sdi)->setJ(angMom);
335 }
336 }
337 w = sqrt(w);
338 for (sdi = hotBin.begin(); sdi != hotBin.end(); ++sdi) {
339 if (rnemdFluxType_ == rnemdFullKE) {
340 vel = (*sdi)->getVel();
341 vel.x() *= x;
342 vel.y() *= y;
343 vel.z() *= z;
344 (*sdi)->setVel(vel);
345 }
346 if ((*sdi)->isDirectional()) {
347 Vector3d angMom = (*sdi)->getJ() * w;
348 (*sdi)->setJ(angMom);
349 }
350 }
351 successfulScale = true;
352 kineticExchange_ += kineticTarget_;
353 }
354 }
355 } else {
356 RealType a000(0.0), a110(0.0), c0(0.0);
357 RealType a001(0.0), a111(0.0), b01(0.0), b11(0.0), c1(0.0);
358 switch (rnemdFluxType_) {
359 case rnemdKE:
360 /* used hotBin coeff's & only scale x & y dimensions
361 RealType px = Phx / Pcx;
362 RealType py = Phy / Pcy;
363 a110 = Khy;
364 c0 = - Khx - Khy - kineticTarget_;
365 a000 = Khx;
366 a111 = Kcy * py * py;
367 b11 = -2.0 * Kcy * py * (1.0 + py);
368 c1 = Kcy * py * (2.0 + py) + Kcx * px * ( 2.0 + px) + kineticTarget_;
369 b01 = -2.0 * Kcx * px * (1.0 + px);
370 a001 = Kcx * px * px;
371 */
372 // scale all three dimensions, let c_x = c_y
373 a000 = Kcx + Kcy;
374 a110 = Kcz;
375 c0 = kineticTarget_ - Kcx - Kcy - Kcz;
376 a001 = Khx * px * px + Khy * py * py;
377 a111 = Khz * pz * pz;
378 b01 = -2.0 * (Khx * px * (1.0 + px) + Khy * py * (1.0 + py));
379 b11 = -2.0 * Khz * pz * (1.0 + pz);
380 c1 = Khx * px * (2.0 + px) + Khy * py * (2.0 + py) +
381 Khz * pz * (2.0 + pz) - kineticTarget_;
382 break;
383 case rnemdPx:
384 c = 1 - momentumTarget_.x() / Pcx;
385 a000 = Kcy;
386 a110 = Kcz;
387 c0 = Kcx * c * c - Kcx - Kcy - Kcz;
388 a001 = py * py * Khy;
389 a111 = pz * pz * Khz;
390 b01 = -2.0 * Khy * py * (1.0 + py);
391 b11 = -2.0 * Khz * pz * (1.0 + pz);
392 c1 = Khy * py * (2.0 + py) + Khz * pz * (2.0 + pz) +
393 Khx * (fastpow(c * px - px - 1.0, 2) - 1.0);
394 break;
395 case rnemdPy:
396 c = 1 - momentumTarget_.y() / Pcy;
397 a000 = Kcx;
398 a110 = Kcz;
399 c0 = Kcy * c * c - Kcx - Kcy - Kcz;
400 a001 = px * px * Khx;
401 a111 = pz * pz * Khz;
402 b01 = -2.0 * Khx * px * (1.0 + px);
403 b11 = -2.0 * Khz * pz * (1.0 + pz);
404 c1 = Khx * px * (2.0 + px) + Khz * pz * (2.0 + pz) +
405 Khy * (fastpow(c * py - py - 1.0, 2) - 1.0);
406 break;
407 case rnemdPz: // we don't really do this, do we?
408 c = 1 - momentumTarget_.z() / Pcz;
409 a000 = Kcx;
410 a110 = Kcy;
411 c0 = Kcz * c * c - Kcx - Kcy - Kcz;
412 a001 = px * px * Khx;
413 a111 = py * py * Khy;
414 b01 = -2.0 * Khx * px * (1.0 + px);
415 b11 = -2.0 * Khy * py * (1.0 + py);
416 c1 = Khx * px * (2.0 + px) + Khy * py * (2.0 + py) +
417 Khz * (fastpow(c * pz - pz - 1.0, 2) - 1.0);
418 break;
419 default:
420 break;
421 }
422
423 RealType v1 = a000 * a111 - a001 * a110;
424 RealType v2 = a000 * b01;
425 RealType v3 = a000 * b11;
426 RealType v4 = a000 * c1 - a001 * c0;
427 RealType v8 = a110 * b01;
428 RealType v10 = -b01 * c0;
429
430 RealType u0 = v2 * v10 - v4 * v4;
431 RealType u1 = -2.0 * v3 * v4;
432 RealType u2 = -v2 * v8 - v3 * v3 - 2.0 * v1 * v4;
433 RealType u3 = -2.0 * v1 * v3;
434 RealType u4 = -v1 * v1;
435 // rescale coefficients
436 RealType maxAbs = fabs(u0);
437 if (maxAbs < fabs(u1)) maxAbs = fabs(u1);
438 if (maxAbs < fabs(u2)) maxAbs = fabs(u2);
439 if (maxAbs < fabs(u3)) maxAbs = fabs(u3);
440 if (maxAbs < fabs(u4)) maxAbs = fabs(u4);
441 u0 /= maxAbs;
442 u1 /= maxAbs;
443 u2 /= maxAbs;
444 u3 /= maxAbs;
445 u4 /= maxAbs;
446 // max_element(start, end) is also available.
447 Polynomial<RealType> poly; // same as DoublePolynomial poly;
448 poly.setCoefficient(4, u4);
449 poly.setCoefficient(3, u3);
450 poly.setCoefficient(2, u2);
451 poly.setCoefficient(1, u1);
452 poly.setCoefficient(0, u0);
453 vector<RealType> realRoots = poly.FindRealRoots();
454
455 vector<RealType>::iterator ri;
456 RealType r1, r2, alpha0;
457 vector<pair<RealType, RealType>> rps;
458 for (ri = realRoots.begin(); ri != realRoots.end(); ++ri) {
459 r2 = *ri;
460 // Check to see if FindRealRoots() gave the right answer:
461 if (fabs(u0 + r2 * (u1 + r2 * (u2 + r2 * (u3 + r2 * u4)))) > 1e-6) {
462 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
463 "RNEMD Warning: polynomial solve seems to have an error!");
464 painCave.isFatal = 0;
465 simError();
466 failRootCount_++;
467 }
468 // Might not be useful w/o rescaling coefficients
469 alpha0 = -c0 - a110 * r2 * r2;
470 if (alpha0 >= 0.0) {
471 r1 = sqrt(alpha0 / a000);
472 if (fabs(c1 + r1 * (b01 + r1 * a001) + r2 * (b11 + r2 * a111)) <
473 1e-6) {
474 rps.push_back(make_pair(r1, r2));
475 }
476 if (r1 > 1e-6) { // r1 non-negative
477 r1 = -r1;
478 if (fabs(c1 + r1 * (b01 + r1 * a001) + r2 * (b11 + r2 * a111)) <
479 1e-6) {
480 rps.push_back(make_pair(r1, r2));
481 }
482 }
483 }
484 }
485 // Consider combining together the part for solving for the pair
486 // w/ the searching for the best solution part so that we don't
487 // need the pairs vector:
488 if (!rps.empty()) {
489 RealType smallestDiff = HONKING_LARGE_VALUE;
490 RealType diff(0.0);
491 std::pair<RealType, RealType> bestPair = std::make_pair(1.0, 1.0);
492 std::vector<std::pair<RealType, RealType>>::iterator rpi;
493 for (rpi = rps.begin(); rpi != rps.end(); ++rpi) {
494 r1 = (*rpi).first;
495 r2 = (*rpi).second;
496 switch (rnemdFluxType_) {
497 case rnemdKE:
498 diff = fastpow(1.0 - r1, 2) + fastpow(1.0 - r2, 2) +
499 fastpow(r1 * r1 / r2 / r2 - Kcz / Kcx, 2) +
500 fastpow(r1 * r1 / r2 / r2 - Kcz / Kcy, 2);
501 break;
502 case rnemdPx:
503 diff = fastpow(1.0 - r1, 2) + fastpow(1.0 - r2, 2) +
504 fastpow(r1 * r1 / r2 / r2 - Kcz / Kcy, 2);
505 break;
506 case rnemdPy:
507 diff = fastpow(1.0 - r1, 2) + fastpow(1.0 - r2, 2) +
508 fastpow(r1 * r1 / r2 / r2 - Kcz / Kcx, 2);
509 break;
510 case rnemdPz:
511 diff = fastpow(1.0 - r1, 2) + fastpow(1.0 - r2, 2) +
512 fastpow(r1 * r1 / r2 / r2 - Kcy / Kcx, 2);
513 default:
514 break;
515 }
516 if (diff < smallestDiff) {
517 smallestDiff = diff;
518 bestPair = *rpi;
519 }
520 }
521#ifdef IS_MPI
522 if (worldRank == 0) {
523#endif
524 // snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
525 // "RNEMD: roots r1= %lf\tr2 = %lf\n",
526 // bestPair.first, bestPair.second);
527 // painCave.isFatal = 0;
528 // painCave.severity = OPENMD_INFO;
529 // simError();
530#ifdef IS_MPI
531 }
532#endif
533
534 switch (rnemdFluxType_) {
535 case rnemdKE:
536 x = bestPair.first;
537 y = bestPair.first;
538 z = bestPair.second;
539 break;
540 case rnemdPx:
541 x = c;
542 y = bestPair.first;
543 z = bestPair.second;
544 break;
545 case rnemdPy:
546 x = bestPair.first;
547 y = c;
548 z = bestPair.second;
549 break;
550 case rnemdPz:
551 x = bestPair.first;
552 y = bestPair.second;
553 z = c;
554 break;
555 default:
556 break;
557 }
558 vector<StuntDouble*>::iterator sdi;
559 Vector3d vel;
560 for (sdi = coldBin.begin(); sdi != coldBin.end(); ++sdi) {
561 vel = (*sdi)->getVel();
562 vel.x() *= x;
563 vel.y() *= y;
564 vel.z() *= z;
565 (*sdi)->setVel(vel);
566 }
567 // convert to hotBin coefficient
568 x = 1.0 + px * (1.0 - x);
569 y = 1.0 + py * (1.0 - y);
570 z = 1.0 + pz * (1.0 - z);
571 for (sdi = hotBin.begin(); sdi != hotBin.end(); ++sdi) {
572 vel = (*sdi)->getVel();
573 vel.x() *= x;
574 vel.y() *= y;
575 vel.z() *= z;
576 (*sdi)->setVel(vel);
577 }
578 successfulScale = true;
579 switch (rnemdFluxType_) {
580 case rnemdKE:
581 kineticExchange_ += kineticTarget_;
582 break;
583 case rnemdPx:
584 case rnemdPy:
585 case rnemdPz:
586 momentumExchange_ += momentumTarget_;
587 break;
588 default:
589 break;
590 }
591 }
592 }
593 if (successfulScale != true) {
594 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
595 "NIVS exchange NOT performed - roots that solve\n"
596 "\tthe constraint equations may not exist or there may be\n"
597 "\tno selected objects in one or both slabs.\n");
598 painCave.isFatal = 0;
599 painCave.severity = OPENMD_INFO;
600 simError();
601 failTrialCount_++;
602 }
603 }
604} // namespace OpenMD::RNEMD
OpenMD::ForceManager
ForceManager is responsible for calculating both the short range (bonded) interactions and long range...
Definition ForceManager.hpp:80
OpenMD::SimInfo
One of the heavy-weight classes of OpenMD, SimInfo maintains objects and variables relating to the cu...
Definition SimInfo.hpp:93
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::getJ
Vector3d getJ()
Returns the current angular momentum of this stuntDouble (body -fixed).
Definition StuntDouble.hpp:374
OpenMD::StuntDouble::isDirectional
bool isDirectional()
Tests if this stuntDouble is a directional one.
Definition StuntDouble.hpp:159
OpenMD::Vector3::z
Real & z()
Returns reference of the third element of Vector3.
Definition Vector3.hpp:120
OpenMD::Vector3::x
Real & x()
Returns reference of the first element of Vector3.
Definition Vector3.hpp:96
OpenMD::Vector3::y
Real & y()
Returns reference of the second element of Vector3.
Definition Vector3.hpp:108