OpenMD 3.2
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
Swap.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/Swap.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 SwapMethod::SwapMethod(SimInfo* info, ForceManager* forceMan) :
83 RNEMD {info, forceMan} {
84 rnemdMethodLabel_ = "Swap";
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 hasCorrectFlux = hasKineticFlux;
97 break;
98 case rnemdPx:
99 case rnemdPy:
100 case rnemdPz:
101 hasCorrectFlux = hasMomentumFlux;
102 break;
103 default:
104 methodFluxMismatch = true;
105 break;
106 }
107
108 if (methodFluxMismatch) {
109 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
110 "RNEMD: The current method,\n"
111 "\t\tSwap\n"
112 "\tcannot be used with the current flux type, %s\n",
113 rnemdFluxTypeLabel_.c_str());
114 painCave.isFatal = 1;
115 painCave.severity = OPENMD_ERROR;
116 simError();
117 }
118
119 if (!hasCorrectFlux) {
120 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
121 "RNEMD: The current method, Swap, and flux type, %s,\n"
122 "\tdid not have the correct flux value specified. Options\n"
123 "\tinclude: kineticFlux and momentumFlux.\n",
124 rnemdFluxTypeLabel_.c_str());
125 painCave.isFatal = 1;
126 painCave.severity = OPENMD_ERROR;
127 simError();
128 }
129
130 if (hasKineticFlux) {
131 setKineticFlux(rnemdParams->getKineticFlux());
132 } else {
133 setKineticFlux(0.0);
134 }
135
136 if (hasMomentumFlux) {
137 RealType momentumFlux = rnemdParams->getMomentumFlux();
138 std::vector<RealType> momentumFluxVector(3);
139
140 switch (rnemdFluxType_) {
141 case rnemdPx:
142 momentumFluxVector[0] = momentumFlux;
143 break;
144 case rnemdPy:
145 momentumFluxVector[1] = momentumFlux;
146 break;
147 case rnemdPz:
148 momentumFluxVector[2] = momentumFlux;
149 break;
150 default:
151 break;
152 }
153
154 setMomentumFluxVector(momentumFluxVector);
155 }
156 }
157
158 void SwapMethod::doRNEMDImpl(SelectionManager& smanA,
159 SelectionManager& smanB) {
160 if (!doRNEMD_) return;
161 int selei;
162 int selej;
163
164 StuntDouble* sd;
165
166 RealType min_val(0.0);
167 int min_found = 0;
168 StuntDouble* min_sd = NULL;
169
170 RealType max_val(0.0);
171 int max_found = 0;
172 StuntDouble* max_sd = NULL;
173
174 for (sd = smanA.beginSelected(selei); sd != NULL;
175 sd = smanA.nextSelected(selei)) {
176 Vector3d pos = sd->getPos();
177
178 // wrap the stuntdouble's position back into the box:
179
180 if (usePeriodicBoundaryConditions_) currentSnap_->wrapVector(pos);
181
182 RealType mass = sd->getMass();
183 Vector3d vel = sd->getVel();
184 RealType value(0.0);
185
186 switch (rnemdFluxType_) {
187 case rnemdKE:
188
189 value = mass * vel.lengthSquare();
190
191 if (sd->isDirectional()) {
192 Vector3d angMom = sd->getJ();
193 Mat3x3d I = sd->getI();
194
195 if (sd->isLinear()) {
196 int i = sd->linearAxis();
197 int j = (i + 1) % 3;
198 int k = (i + 2) % 3;
199 value += angMom[j] * angMom[j] / I(j, j) +
200 angMom[k] * angMom[k] / I(k, k);
201 } else {
202 value += angMom[0] * angMom[0] / I(0, 0) +
203 angMom[1] * angMom[1] / I(1, 1) +
204 angMom[2] * angMom[2] / I(2, 2);
205 }
206 } // angular momenta exchange enabled
207 value *= 0.5;
208 break;
209 case rnemdPx:
210 value = mass * vel[0];
211 break;
212 case rnemdPy:
213 value = mass * vel[1];
214 break;
215 case rnemdPz:
216 value = mass * vel[2];
217 break;
218 default:
219 break;
220 }
221 if (!max_found) {
222 max_val = value;
223 max_sd = sd;
224 max_found = 1;
225 } else {
226 if (max_val < value) {
227 max_val = value;
228 max_sd = sd;
229 }
230 }
231 }
232
233 for (sd = smanB.beginSelected(selej); sd != NULL;
234 sd = smanB.nextSelected(selej)) {
235 Vector3d pos = sd->getPos();
236
237 // wrap the stuntdouble's position back into the box:
238
239 if (usePeriodicBoundaryConditions_) currentSnap_->wrapVector(pos);
240
241 RealType mass = sd->getMass();
242 Vector3d vel = sd->getVel();
243 RealType value(0.0);
244
245 switch (rnemdFluxType_) {
246 case rnemdKE:
247
248 value = mass * vel.lengthSquare();
249
250 if (sd->isDirectional()) {
251 Vector3d angMom = sd->getJ();
252 Mat3x3d I = sd->getI();
253
254 if (sd->isLinear()) {
255 int i = sd->linearAxis();
256 int j = (i + 1) % 3;
257 int k = (i + 2) % 3;
258 value += angMom[j] * angMom[j] / I(j, j) +
259 angMom[k] * angMom[k] / I(k, k);
260 } else {
261 value += angMom[0] * angMom[0] / I(0, 0) +
262 angMom[1] * angMom[1] / I(1, 1) +
263 angMom[2] * angMom[2] / I(2, 2);
264 }
265 } // angular momenta exchange enabled
266 value *= 0.5;
267 break;
268 case rnemdPx:
269 value = mass * vel[0];
270 break;
271 case rnemdPy:
272 value = mass * vel[1];
273 break;
274 case rnemdPz:
275 value = mass * vel[2];
276 break;
277 default:
278 break;
279 }
280
281 if (!min_found) {
282 min_val = value;
283 min_sd = sd;
284 min_found = 1;
285 } else {
286 if (min_val > value) {
287 min_val = value;
288 min_sd = sd;
289 }
290 }
291 }
292
293#ifdef IS_MPI
294 int worldRank;
295 MPI_Comm_rank(MPI_COMM_WORLD, &worldRank);
296
297 int my_min_found = min_found;
298 int my_max_found = max_found;
299
300 // Even if we didn't find a minimum, did someone else?
301 MPI_Allreduce(&my_min_found, &min_found, 1, MPI_INT, MPI_LOR,
302 MPI_COMM_WORLD);
303 // Even if we didn't find a maximum, did someone else?
304 MPI_Allreduce(&my_max_found, &max_found, 1, MPI_INT, MPI_LOR,
305 MPI_COMM_WORLD);
306#endif
307
308 if (max_found && min_found) {
309#ifdef IS_MPI
310 struct {
311 RealType val;
312 int rank;
313 } max_vals, min_vals;
314
315 if (my_min_found) {
316 min_vals.val = min_val;
317 } else {
318 min_vals.val = HONKING_LARGE_VALUE;
319 }
320 min_vals.rank = worldRank;
321
322 // Who had the minimum?
323 MPI_Allreduce(&min_vals, &min_vals, 1, MPI_REALTYPE_INT, MPI_MINLOC,
324 MPI_COMM_WORLD);
325 min_val = min_vals.val;
326
327 if (my_max_found) {
328 max_vals.val = max_val;
329 } else {
330 max_vals.val = -HONKING_LARGE_VALUE;
331 }
332 max_vals.rank = worldRank;
333
334 // Who had the maximum?
335 MPI_Allreduce(&max_vals, &max_vals, 1, MPI_REALTYPE_INT, MPI_MAXLOC,
336 MPI_COMM_WORLD);
337 max_val = max_vals.val;
338#endif
339
340 if (min_val < max_val) {
341#ifdef IS_MPI
342 if (max_vals.rank == worldRank && min_vals.rank == worldRank) {
343 // I have both maximum and minimum, so proceed like a single
344 // processor version:
345#endif
346
347 Vector3d min_vel = min_sd->getVel();
348 Vector3d max_vel = max_sd->getVel();
349 RealType temp_vel;
350
351 switch (rnemdFluxType_) {
352 case rnemdKE:
353 min_sd->setVel(max_vel);
354 max_sd->setVel(min_vel);
355 if (min_sd->isDirectional() && max_sd->isDirectional()) {
356 Vector3d min_angMom = min_sd->getJ();
357 Vector3d max_angMom = max_sd->getJ();
358 min_sd->setJ(max_angMom);
359 max_sd->setJ(min_angMom);
360 } // angular momenta exchange enabled
361 // assumes same rigid body identity
362 break;
363 case rnemdPx:
364 temp_vel = min_vel.x();
365 min_vel.x() = max_vel.x();
366 max_vel.x() = temp_vel;
367 min_sd->setVel(min_vel);
368 max_sd->setVel(max_vel);
369 break;
370 case rnemdPy:
371 temp_vel = min_vel.y();
372 min_vel.y() = max_vel.y();
373 max_vel.y() = temp_vel;
374 min_sd->setVel(min_vel);
375 max_sd->setVel(max_vel);
376 break;
377 case rnemdPz:
378 temp_vel = min_vel.z();
379 min_vel.z() = max_vel.z();
380 max_vel.z() = temp_vel;
381 min_sd->setVel(min_vel);
382 max_sd->setVel(max_vel);
383 break;
384 default:
385 break;
386 }
387
388#ifdef IS_MPI
389 // the rest of the cases only apply in parallel simulations:
390 } else if (max_vals.rank == worldRank) {
391 // I had the max, but not the minimum
392
393 Vector3d min_vel;
394 Vector3d max_vel = max_sd->getVel();
395 MPI_Status status;
396
397 // point-to-point swap of the velocity vector
398 MPI_Sendrecv(max_vel.getArrayPointer(), 3, MPI_REALTYPE,
399 min_vals.rank, 0, min_vel.getArrayPointer(), 3,
400 MPI_REALTYPE, min_vals.rank, 0, MPI_COMM_WORLD, &status);
401
402 switch (rnemdFluxType_) {
403 case rnemdKE:
404 max_sd->setVel(min_vel);
405 // angular momenta exchange enabled
406 if (max_sd->isDirectional()) {
407 Vector3d min_angMom;
408 Vector3d max_angMom = max_sd->getJ();
409
410 // point-to-point swap of the angular momentum vector
411 MPI_Sendrecv(max_angMom.getArrayPointer(), 3, MPI_REALTYPE,
412 min_vals.rank, 1, min_angMom.getArrayPointer(), 3,
413 MPI_REALTYPE, min_vals.rank, 1, MPI_COMM_WORLD,
414 &status);
415
416 max_sd->setJ(min_angMom);
417 }
418 break;
419 case rnemdPx:
420 max_vel.x() = min_vel.x();
421 max_sd->setVel(max_vel);
422 break;
423 case rnemdPy:
424 max_vel.y() = min_vel.y();
425 max_sd->setVel(max_vel);
426 break;
427 case rnemdPz:
428 max_vel.z() = min_vel.z();
429 max_sd->setVel(max_vel);
430 break;
431 default:
432 break;
433 }
434 } else if (min_vals.rank == worldRank) {
435 // I had the minimum but not the maximum:
436
437 Vector3d max_vel;
438 Vector3d min_vel = min_sd->getVel();
439 MPI_Status status;
440
441 // point-to-point swap of the velocity vector
442 MPI_Sendrecv(min_vel.getArrayPointer(), 3, MPI_REALTYPE,
443 max_vals.rank, 0, max_vel.getArrayPointer(), 3,
444 MPI_REALTYPE, max_vals.rank, 0, MPI_COMM_WORLD, &status);
445
446 switch (rnemdFluxType_) {
447 case rnemdKE:
448 min_sd->setVel(max_vel);
449 // angular momenta exchange enabled
450 if (min_sd->isDirectional()) {
451 Vector3d min_angMom = min_sd->getJ();
452 Vector3d max_angMom;
453
454 // point-to-point swap of the angular momentum vector
455 MPI_Sendrecv(min_angMom.getArrayPointer(), 3, MPI_REALTYPE,
456 max_vals.rank, 1, max_angMom.getArrayPointer(), 3,
457 MPI_REALTYPE, max_vals.rank, 1, MPI_COMM_WORLD,
458 &status);
459
460 min_sd->setJ(max_angMom);
461 }
462 break;
463 case rnemdPx:
464 min_vel.x() = max_vel.x();
465 min_sd->setVel(min_vel);
466 break;
467 case rnemdPy:
468 min_vel.y() = max_vel.y();
469 min_sd->setVel(min_vel);
470 break;
471 case rnemdPz:
472 min_vel.z() = max_vel.z();
473 min_sd->setVel(min_vel);
474 break;
475 default:
476 break;
477 }
478 }
479#endif
480
481 switch (rnemdFluxType_) {
482 case rnemdKE:
483 kineticExchange_ += max_val - min_val;
484 break;
485 case rnemdPx:
486 momentumExchange_.x() += max_val - min_val;
487 break;
488 case rnemdPy:
489 momentumExchange_.y() += max_val - min_val;
490 break;
491 case rnemdPz:
492 momentumExchange_.z() += max_val - min_val;
493 break;
494 default:
495 break;
496 }
497 } else {
498 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
499 "RNEMD::doSwap exchange NOT performed "
500 "because min_val > max_val\n");
501 painCave.isFatal = 0;
502 painCave.severity = OPENMD_INFO;
503 simError();
504 failTrialCount_++;
505 }
506 } else {
507 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
508 "Swap exchange NOT performed because selected object\n"
509 "\twas not present in at least one of the two slabs.\n");
510 painCave.isFatal = 0;
511 painCave.severity = OPENMD_INFO;
512 simError();
513 failTrialCount_++;
514 }
515 }
516} // namespace OpenMD::RNEMD