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