ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/rnemd/RNEMD.cpp
Revision: 1969
Committed: Wed Feb 26 14:14:50 2014 UTC (11 years, 2 months ago) by gezelter
File size: 76508 byte(s)
Log Message: 
Fixes to deal with deprecation of MPI C++ bindings.  We've reverted back to the
C calls.

File Contents

# Content
1 /*
2 * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3 *
4 * The University of Notre Dame grants you ("Licensee") a
5 * non-exclusive, royalty free, license to use, modify and
6 * redistribute this software in source and binary code form, provided
7 * that the following conditions are met:
8 *
9 * 1. Redistributions of source code must retain the above copyright
10 * notice, this list of conditions and the following disclaimer.
11 *
12 * 2. Redistributions in binary form must reproduce the above copyright
13 * notice, this list of conditions and the following disclaimer in the
14 * documentation and/or other materials provided with the
15 * distribution.
16 *
17 * This software is provided "AS IS," without a warranty of any
18 * kind. All express or implied conditions, representations and
19 * warranties, including any implied warranty of merchantability,
20 * fitness for a particular purpose or non-infringement, are hereby
21 * excluded. The University of Notre Dame and its licensors shall not
22 * be liable for any damages suffered by licensee as a result of
23 * using, modifying or distributing the software or its
24 * derivatives. In no event will the University of Notre Dame or its
25 * licensors be liable for any lost revenue, profit or data, or for
26 * direct, indirect, special, consequential, incidental or punitive
27 * damages, however caused and regardless of the theory of liability,
28 * arising out of the use of or inability to use software, even if the
29 * University of Notre Dame has been advised of the possibility of
30 * such damages.
31 *
32 * SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your
33 * research, please cite the appropriate papers when you publish your
34 * work. Good starting points are:
35 *
36 * [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
37 * [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
38 * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).
39 * [4] Vardeman & Gezelter, in progress (2009).
40 */
41 #ifdef IS_MPI
42 #include <mpi.h>
43 #endif
44
45 #include <cmath>
46 #include <sstream>
47 #include <string>
48
49 #include "rnemd/RNEMD.hpp"
50 #include "math/Vector3.hpp"
51 #include "math/Vector.hpp"
52 #include "math/SquareMatrix3.hpp"
53 #include "math/Polynomial.hpp"
54 #include "primitives/Molecule.hpp"
55 #include "primitives/StuntDouble.hpp"
56 #include "utils/PhysicalConstants.hpp"
57 #include "utils/Tuple.hpp"
58 #include "brains/Thermo.hpp"
59 #include "math/ConvexHull.hpp"
60
61 #ifdef _MSC_VER
62 #define isnan(x) _isnan((x))
63 #define isinf(x) (!_finite(x) && !_isnan(x))
64 #endif
65
66 #define HONKING_LARGE_VALUE 1.0e10
67
68 using namespace std;
69 namespace OpenMD {
70
71 RNEMD::RNEMD(SimInfo* info) : info_(info), evaluator_(info), seleMan_(info),
72 evaluatorA_(info), seleManA_(info),
73 commonA_(info), evaluatorB_(info),
74 seleManB_(info), commonB_(info),
75 hasData_(false), hasDividingArea_(false),
76 usePeriodicBoundaryConditions_(info->getSimParams()->getUsePeriodicBoundaryConditions()) {
77
78 trialCount_ = 0;
79 failTrialCount_ = 0;
80 failRootCount_ = 0;
81
82 Globals* simParams = info->getSimParams();
83 RNEMDParameters* rnemdParams = simParams->getRNEMDParameters();
84
85 doRNEMD_ = rnemdParams->getUseRNEMD();
86 if (!doRNEMD_) return;
87
88 stringToMethod_["Swap"] = rnemdSwap;
89 stringToMethod_["NIVS"] = rnemdNIVS;
90 stringToMethod_["VSS"] = rnemdVSS;
91
92 stringToFluxType_["KE"] = rnemdKE;
93 stringToFluxType_["Px"] = rnemdPx;
94 stringToFluxType_["Py"] = rnemdPy;
95 stringToFluxType_["Pz"] = rnemdPz;
96 stringToFluxType_["Pvector"] = rnemdPvector;
97 stringToFluxType_["Lx"] = rnemdLx;
98 stringToFluxType_["Ly"] = rnemdLy;
99 stringToFluxType_["Lz"] = rnemdLz;
100 stringToFluxType_["Lvector"] = rnemdLvector;
101 stringToFluxType_["KE+Px"] = rnemdKePx;
102 stringToFluxType_["KE+Py"] = rnemdKePy;
103 stringToFluxType_["KE+Pvector"] = rnemdKePvector;
104 stringToFluxType_["KE+Lx"] = rnemdKeLx;
105 stringToFluxType_["KE+Ly"] = rnemdKeLy;
106 stringToFluxType_["KE+Lz"] = rnemdKeLz;
107 stringToFluxType_["KE+Lvector"] = rnemdKeLvector;
108
109 runTime_ = simParams->getRunTime();
110 statusTime_ = simParams->getStatusTime();
111
112 const string methStr = rnemdParams->getMethod();
113 bool hasFluxType = rnemdParams->haveFluxType();
114
115 rnemdObjectSelection_ = rnemdParams->getObjectSelection();
116
117 string fluxStr;
118 if (hasFluxType) {
119 fluxStr = rnemdParams->getFluxType();
120 } else {
121 sprintf(painCave.errMsg,
122 "RNEMD: No fluxType was set in the md file. This parameter,\n"
123 "\twhich must be one of the following values:\n"
124 "\tKE, Px, Py, Pz, Pvector, Lx, Ly, Lz, Lvector,\n"
125 "\tKE+Px, KE+Py, KE+Pvector, KE+Lx, KE+Ly, KE+Lz, KE+Lvector\n"
126 "\tmust be set to use RNEMD\n");
127 painCave.isFatal = 1;
128 painCave.severity = OPENMD_ERROR;
129 simError();
130 }
131
132 bool hasKineticFlux = rnemdParams->haveKineticFlux();
133 bool hasMomentumFlux = rnemdParams->haveMomentumFlux();
134 bool hasMomentumFluxVector = rnemdParams->haveMomentumFluxVector();
135 bool hasAngularMomentumFlux = rnemdParams->haveAngularMomentumFlux();
136 bool hasAngularMomentumFluxVector = rnemdParams->haveAngularMomentumFluxVector();
137 hasSelectionA_ = rnemdParams->haveSelectionA();
138 hasSelectionB_ = rnemdParams->haveSelectionB();
139 bool hasSlabWidth = rnemdParams->haveSlabWidth();
140 bool hasSlabACenter = rnemdParams->haveSlabACenter();
141 bool hasSlabBCenter = rnemdParams->haveSlabBCenter();
142 bool hasSphereARadius = rnemdParams->haveSphereARadius();
143 hasSphereBRadius_ = rnemdParams->haveSphereBRadius();
144 bool hasCoordinateOrigin = rnemdParams->haveCoordinateOrigin();
145 bool hasOutputFileName = rnemdParams->haveOutputFileName();
146 bool hasOutputFields = rnemdParams->haveOutputFields();
147
148 map<string, RNEMDMethod>::iterator i;
149 i = stringToMethod_.find(methStr);
150 if (i != stringToMethod_.end())
151 rnemdMethod_ = i->second;
152 else {
153 sprintf(painCave.errMsg,
154 "RNEMD: The current method,\n"
155 "\t\t%s is not one of the recognized\n"
156 "\texchange methods: Swap, NIVS, or VSS\n",
157 methStr.c_str());
158 painCave.isFatal = 1;
159 painCave.severity = OPENMD_ERROR;
160 simError();
161 }
162
163 map<string, RNEMDFluxType>::iterator j;
164 j = stringToFluxType_.find(fluxStr);
165 if (j != stringToFluxType_.end())
166 rnemdFluxType_ = j->second;
167 else {
168 sprintf(painCave.errMsg,
169 "RNEMD: The current fluxType,\n"
170 "\t\t%s\n"
171 "\tis not one of the recognized flux types.\n",
172 fluxStr.c_str());
173 painCave.isFatal = 1;
174 painCave.severity = OPENMD_ERROR;
175 simError();
176 }
177
178 bool methodFluxMismatch = false;
179 bool hasCorrectFlux = false;
180 switch(rnemdMethod_) {
181 case rnemdSwap:
182 switch (rnemdFluxType_) {
183 case rnemdKE:
184 hasCorrectFlux = hasKineticFlux;
185 break;
186 case rnemdPx:
187 case rnemdPy:
188 case rnemdPz:
189 hasCorrectFlux = hasMomentumFlux;
190 break;
191 default :
192 methodFluxMismatch = true;
193 break;
194 }
195 break;
196 case rnemdNIVS:
197 switch (rnemdFluxType_) {
198 case rnemdKE:
199 case rnemdRotKE:
200 case rnemdFullKE:
201 hasCorrectFlux = hasKineticFlux;
202 break;
203 case rnemdPx:
204 case rnemdPy:
205 case rnemdPz:
206 hasCorrectFlux = hasMomentumFlux;
207 break;
208 case rnemdKePx:
209 case rnemdKePy:
210 hasCorrectFlux = hasMomentumFlux && hasKineticFlux;
211 break;
212 default:
213 methodFluxMismatch = true;
214 break;
215 }
216 break;
217 case rnemdVSS:
218 switch (rnemdFluxType_) {
219 case rnemdKE:
220 case rnemdRotKE:
221 case rnemdFullKE:
222 hasCorrectFlux = hasKineticFlux;
223 break;
224 case rnemdPx:
225 case rnemdPy:
226 case rnemdPz:
227 hasCorrectFlux = hasMomentumFlux;
228 break;
229 case rnemdLx:
230 case rnemdLy:
231 case rnemdLz:
232 hasCorrectFlux = hasAngularMomentumFlux;
233 break;
234 case rnemdPvector:
235 hasCorrectFlux = hasMomentumFluxVector;
236 break;
237 case rnemdLvector:
238 hasCorrectFlux = hasAngularMomentumFluxVector;
239 break;
240 case rnemdKePx:
241 case rnemdKePy:
242 hasCorrectFlux = hasMomentumFlux && hasKineticFlux;
243 break;
244 case rnemdKeLx:
245 case rnemdKeLy:
246 case rnemdKeLz:
247 hasCorrectFlux = hasAngularMomentumFlux && hasKineticFlux;
248 break;
249 case rnemdKePvector:
250 hasCorrectFlux = hasMomentumFluxVector && hasKineticFlux;
251 break;
252 case rnemdKeLvector:
253 hasCorrectFlux = hasAngularMomentumFluxVector && hasKineticFlux;
254 break;
255 default:
256 methodFluxMismatch = true;
257 break;
258 }
259 default:
260 break;
261 }
262
263 if (methodFluxMismatch) {
264 sprintf(painCave.errMsg,
265 "RNEMD: The current method,\n"
266 "\t\t%s\n"
267 "\tcannot be used with the current flux type, %s\n",
268 methStr.c_str(), fluxStr.c_str());
269 painCave.isFatal = 1;
270 painCave.severity = OPENMD_ERROR;
271 simError();
272 }
273 if (!hasCorrectFlux) {
274 sprintf(painCave.errMsg,
275 "RNEMD: The current method, %s, and flux type, %s,\n"
276 "\tdid not have the correct flux value specified. Options\n"
277 "\tinclude: kineticFlux, momentumFlux, angularMomentumFlux,\n"
278 "\tmomentumFluxVector, and angularMomentumFluxVector.\n",
279 methStr.c_str(), fluxStr.c_str());
280 painCave.isFatal = 1;
281 painCave.severity = OPENMD_ERROR;
282 simError();
283 }
284
285 if (hasKineticFlux) {
286 // convert the kcal / mol / Angstroms^2 / fs values in the md file
287 // into amu / fs^3:
288 kineticFlux_ = rnemdParams->getKineticFlux()
289 * PhysicalConstants::energyConvert;
290 } else {
291 kineticFlux_ = 0.0;
292 }
293 if (hasMomentumFluxVector) {
294 momentumFluxVector_ = rnemdParams->getMomentumFluxVector();
295 } else {
296 momentumFluxVector_ = V3Zero;
297 if (hasMomentumFlux) {
298 RealType momentumFlux = rnemdParams->getMomentumFlux();
299 switch (rnemdFluxType_) {
300 case rnemdPx:
301 momentumFluxVector_.x() = momentumFlux;
302 break;
303 case rnemdPy:
304 momentumFluxVector_.y() = momentumFlux;
305 break;
306 case rnemdPz:
307 momentumFluxVector_.z() = momentumFlux;
308 break;
309 case rnemdKePx:
310 momentumFluxVector_.x() = momentumFlux;
311 break;
312 case rnemdKePy:
313 momentumFluxVector_.y() = momentumFlux;
314 break;
315 default:
316 break;
317 }
318 }
319 if (hasAngularMomentumFluxVector) {
320 angularMomentumFluxVector_ = rnemdParams->getAngularMomentumFluxVector();
321 } else {
322 angularMomentumFluxVector_ = V3Zero;
323 if (hasAngularMomentumFlux) {
324 RealType angularMomentumFlux = rnemdParams->getAngularMomentumFlux();
325 switch (rnemdFluxType_) {
326 case rnemdLx:
327 angularMomentumFluxVector_.x() = angularMomentumFlux;
328 break;
329 case rnemdLy:
330 angularMomentumFluxVector_.y() = angularMomentumFlux;
331 break;
332 case rnemdLz:
333 angularMomentumFluxVector_.z() = angularMomentumFlux;
334 break;
335 case rnemdKeLx:
336 angularMomentumFluxVector_.x() = angularMomentumFlux;
337 break;
338 case rnemdKeLy:
339 angularMomentumFluxVector_.y() = angularMomentumFlux;
340 break;
341 case rnemdKeLz:
342 angularMomentumFluxVector_.z() = angularMomentumFlux;
343 break;
344 default:
345 break;
346 }
347 }
348 }
349
350 if (hasCoordinateOrigin) {
351 coordinateOrigin_ = rnemdParams->getCoordinateOrigin();
352 } else {
353 coordinateOrigin_ = V3Zero;
354 }
355
356 // do some sanity checking
357
358 int selectionCount = seleMan_.getSelectionCount();
359
360 int nIntegrable = info->getNGlobalIntegrableObjects();
361
362 if (selectionCount > nIntegrable) {
363 sprintf(painCave.errMsg,
364 "RNEMD: The current objectSelection,\n"
365 "\t\t%s\n"
366 "\thas resulted in %d selected objects. However,\n"
367 "\tthe total number of integrable objects in the system\n"
368 "\tis only %d. This is almost certainly not what you want\n"
369 "\tto do. A likely cause of this is forgetting the _RB_0\n"
370 "\tselector in the selection script!\n",
371 rnemdObjectSelection_.c_str(),
372 selectionCount, nIntegrable);
373 painCave.isFatal = 0;
374 painCave.severity = OPENMD_WARNING;
375 simError();
376 }
377
378 areaAccumulator_ = new Accumulator();
379
380 nBins_ = rnemdParams->getOutputBins();
381 binWidth_ = rnemdParams->getOutputBinWidth();
382
383 data_.resize(RNEMD::ENDINDEX);
384 OutputData z;
385 z.units = "Angstroms";
386 z.title = "Z";
387 z.dataType = "RealType";
388 z.accumulator.reserve(nBins_);
389 for (int i = 0; i < nBins_; i++)
390 z.accumulator.push_back( new Accumulator() );
391 data_[Z] = z;
392 outputMap_["Z"] = Z;
393
394 OutputData r;
395 r.units = "Angstroms";
396 r.title = "R";
397 r.dataType = "RealType";
398 r.accumulator.reserve(nBins_);
399 for (int i = 0; i < nBins_; i++)
400 r.accumulator.push_back( new Accumulator() );
401 data_[R] = r;
402 outputMap_["R"] = R;
403
404 OutputData temperature;
405 temperature.units = "K";
406 temperature.title = "Temperature";
407 temperature.dataType = "RealType";
408 temperature.accumulator.reserve(nBins_);
409 for (int i = 0; i < nBins_; i++)
410 temperature.accumulator.push_back( new Accumulator() );
411 data_[TEMPERATURE] = temperature;
412 outputMap_["TEMPERATURE"] = TEMPERATURE;
413
414 OutputData velocity;
415 velocity.units = "angstroms/fs";
416 velocity.title = "Velocity";
417 velocity.dataType = "Vector3d";
418 velocity.accumulator.reserve(nBins_);
419 for (int i = 0; i < nBins_; i++)
420 velocity.accumulator.push_back( new VectorAccumulator() );
421 data_[VELOCITY] = velocity;
422 outputMap_["VELOCITY"] = VELOCITY;
423
424 OutputData angularVelocity;
425 angularVelocity.units = "angstroms^2/fs";
426 angularVelocity.title = "AngularVelocity";
427 angularVelocity.dataType = "Vector3d";
428 angularVelocity.accumulator.reserve(nBins_);
429 for (int i = 0; i < nBins_; i++)
430 angularVelocity.accumulator.push_back( new VectorAccumulator() );
431 data_[ANGULARVELOCITY] = angularVelocity;
432 outputMap_["ANGULARVELOCITY"] = ANGULARVELOCITY;
433
434 OutputData density;
435 density.units = "g cm^-3";
436 density.title = "Density";
437 density.dataType = "RealType";
438 density.accumulator.reserve(nBins_);
439 for (int i = 0; i < nBins_; i++)
440 density.accumulator.push_back( new Accumulator() );
441 data_[DENSITY] = density;
442 outputMap_["DENSITY"] = DENSITY;
443
444 if (hasOutputFields) {
445 parseOutputFileFormat(rnemdParams->getOutputFields());
446 } else {
447 if (usePeriodicBoundaryConditions_)
448 outputMask_.set(Z);
449 else
450 outputMask_.set(R);
451 switch (rnemdFluxType_) {
452 case rnemdKE:
453 case rnemdRotKE:
454 case rnemdFullKE:
455 outputMask_.set(TEMPERATURE);
456 break;
457 case rnemdPx:
458 case rnemdPy:
459 outputMask_.set(VELOCITY);
460 break;
461 case rnemdPz:
462 case rnemdPvector:
463 outputMask_.set(VELOCITY);
464 outputMask_.set(DENSITY);
465 break;
466 case rnemdLx:
467 case rnemdLy:
468 case rnemdLz:
469 case rnemdLvector:
470 outputMask_.set(ANGULARVELOCITY);
471 break;
472 case rnemdKeLx:
473 case rnemdKeLy:
474 case rnemdKeLz:
475 case rnemdKeLvector:
476 outputMask_.set(TEMPERATURE);
477 outputMask_.set(ANGULARVELOCITY);
478 break;
479 case rnemdKePx:
480 case rnemdKePy:
481 outputMask_.set(TEMPERATURE);
482 outputMask_.set(VELOCITY);
483 break;
484 case rnemdKePvector:
485 outputMask_.set(TEMPERATURE);
486 outputMask_.set(VELOCITY);
487 outputMask_.set(DENSITY);
488 break;
489 default:
490 break;
491 }
492 }
493
494 if (hasOutputFileName) {
495 rnemdFileName_ = rnemdParams->getOutputFileName();
496 } else {
497 rnemdFileName_ = getPrefix(info->getFinalConfigFileName()) + ".rnemd";
498 }
499
500 exchangeTime_ = rnemdParams->getExchangeTime();
501
502 Snapshot* currentSnap_ = info->getSnapshotManager()->getCurrentSnapshot();
503 // total exchange sums are zeroed out at the beginning:
504
505 kineticExchange_ = 0.0;
506 momentumExchange_ = V3Zero;
507 angularMomentumExchange_ = V3Zero;
508
509 std::ostringstream selectionAstream;
510 std::ostringstream selectionBstream;
511
512 if (hasSelectionA_) {
513 selectionA_ = rnemdParams->getSelectionA();
514 } else {
515 if (usePeriodicBoundaryConditions_) {
516 Mat3x3d hmat = currentSnap_->getHmat();
517
518 if (hasSlabWidth)
519 slabWidth_ = rnemdParams->getSlabWidth();
520 else
521 slabWidth_ = hmat(2,2) / 10.0;
522
523 if (hasSlabACenter)
524 slabACenter_ = rnemdParams->getSlabACenter();
525 else
526 slabACenter_ = 0.0;
527
528 selectionAstream << "select wrappedz > "
529 << slabACenter_ - 0.5*slabWidth_
530 << " && wrappedz < "
531 << slabACenter_ + 0.5*slabWidth_;
532 selectionA_ = selectionAstream.str();
533 } else {
534 if (hasSphereARadius)
535 sphereARadius_ = rnemdParams->getSphereARadius();
536 else {
537 // use an initial guess to the size of the inner slab to be 1/10 the
538 // radius of an approximately spherical hull:
539 Thermo thermo(info);
540 RealType hVol = thermo.getHullVolume();
541 sphereARadius_ = 0.1 * pow((3.0 * hVol / (4.0 * M_PI)), 1.0/3.0);
542 }
543 selectionAstream << "select r < " << sphereARadius_;
544 selectionA_ = selectionAstream.str();
545 }
546 }
547
548 if (hasSelectionB_) {
549 selectionB_ = rnemdParams->getSelectionB();
550
551 } else {
552 if (usePeriodicBoundaryConditions_) {
553 Mat3x3d hmat = currentSnap_->getHmat();
554
555 if (hasSlabWidth)
556 slabWidth_ = rnemdParams->getSlabWidth();
557 else
558 slabWidth_ = hmat(2,2) / 10.0;
559
560 if (hasSlabBCenter)
561 slabBCenter_ = rnemdParams->getSlabBCenter();
562 else
563 slabBCenter_ = hmat(2,2) / 2.0;
564
565 selectionBstream << "select wrappedz > "
566 << slabBCenter_ - 0.5*slabWidth_
567 << " && wrappedz < "
568 << slabBCenter_ + 0.5*slabWidth_;
569 selectionB_ = selectionBstream.str();
570 } else {
571 if (hasSphereBRadius_) {
572 sphereBRadius_ = rnemdParams->getSphereBRadius();
573 selectionBstream << "select r > " << sphereBRadius_;
574 selectionB_ = selectionBstream.str();
575 } else {
576 selectionB_ = "select hull";
577 BisHull_ = true;
578 hasSelectionB_ = true;
579 }
580 }
581 }
582 }
583
584 // object evaluator:
585 evaluator_.loadScriptString(rnemdObjectSelection_);
586 seleMan_.setSelectionSet(evaluator_.evaluate());
587 evaluatorA_.loadScriptString(selectionA_);
588 evaluatorB_.loadScriptString(selectionB_);
589 seleManA_.setSelectionSet(evaluatorA_.evaluate());
590 seleManB_.setSelectionSet(evaluatorB_.evaluate());
591 commonA_ = seleManA_ & seleMan_;
592 commonB_ = seleManB_ & seleMan_;
593 }
594
595
596 RNEMD::~RNEMD() {
597 if (!doRNEMD_) return;
598 #ifdef IS_MPI
599 if (worldRank == 0) {
600 #endif
601
602 writeOutputFile();
603
604 rnemdFile_.close();
605
606 #ifdef IS_MPI
607 }
608 #endif
609
610 // delete all of the objects we created:
611 delete areaAccumulator_;
612 data_.clear();
613 }
614
615 void RNEMD::doSwap(SelectionManager& smanA, SelectionManager& smanB) {
616 if (!doRNEMD_) return;
617 int selei;
618 int selej;
619
620 Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
621 Mat3x3d hmat = currentSnap_->getHmat();
622
623 StuntDouble* sd;
624
625 RealType min_val;
626 bool min_found = false;
627 StuntDouble* min_sd;
628
629 RealType max_val;
630 bool max_found = false;
631 StuntDouble* max_sd;
632
633 for (sd = seleManA_.beginSelected(selei); sd != NULL;
634 sd = seleManA_.nextSelected(selei)) {
635
636 Vector3d pos = sd->getPos();
637
638 // wrap the stuntdouble's position back into the box:
639
640 if (usePeriodicBoundaryConditions_)
641 currentSnap_->wrapVector(pos);
642
643 RealType mass = sd->getMass();
644 Vector3d vel = sd->getVel();
645 RealType value;
646
647 switch(rnemdFluxType_) {
648 case rnemdKE :
649
650 value = mass * vel.lengthSquare();
651
652 if (sd->isDirectional()) {
653 Vector3d angMom = sd->getJ();
654 Mat3x3d I = sd->getI();
655
656 if (sd->isLinear()) {
657 int i = sd->linearAxis();
658 int j = (i + 1) % 3;
659 int k = (i + 2) % 3;
660 value += angMom[j] * angMom[j] / I(j, j) +
661 angMom[k] * angMom[k] / I(k, k);
662 } else {
663 value += angMom[0]*angMom[0]/I(0, 0)
664 + angMom[1]*angMom[1]/I(1, 1)
665 + angMom[2]*angMom[2]/I(2, 2);
666 }
667 } //angular momenta exchange enabled
668 value *= 0.5;
669 break;
670 case rnemdPx :
671 value = mass * vel[0];
672 break;
673 case rnemdPy :
674 value = mass * vel[1];
675 break;
676 case rnemdPz :
677 value = mass * vel[2];
678 break;
679 default :
680 break;
681 }
682 if (!max_found) {
683 max_val = value;
684 max_sd = sd;
685 max_found = true;
686 } else {
687 if (max_val < value) {
688 max_val = value;
689 max_sd = sd;
690 }
691 }
692 }
693
694 for (sd = seleManB_.beginSelected(selej); sd != NULL;
695 sd = seleManB_.nextSelected(selej)) {
696
697 Vector3d pos = sd->getPos();
698
699 // wrap the stuntdouble's position back into the box:
700
701 if (usePeriodicBoundaryConditions_)
702 currentSnap_->wrapVector(pos);
703
704 RealType mass = sd->getMass();
705 Vector3d vel = sd->getVel();
706 RealType value;
707
708 switch(rnemdFluxType_) {
709 case rnemdKE :
710
711 value = mass * vel.lengthSquare();
712
713 if (sd->isDirectional()) {
714 Vector3d angMom = sd->getJ();
715 Mat3x3d I = sd->getI();
716
717 if (sd->isLinear()) {
718 int i = sd->linearAxis();
719 int j = (i + 1) % 3;
720 int k = (i + 2) % 3;
721 value += angMom[j] * angMom[j] / I(j, j) +
722 angMom[k] * angMom[k] / I(k, k);
723 } else {
724 value += angMom[0]*angMom[0]/I(0, 0)
725 + angMom[1]*angMom[1]/I(1, 1)
726 + angMom[2]*angMom[2]/I(2, 2);
727 }
728 } //angular momenta exchange enabled
729 value *= 0.5;
730 break;
731 case rnemdPx :
732 value = mass * vel[0];
733 break;
734 case rnemdPy :
735 value = mass * vel[1];
736 break;
737 case rnemdPz :
738 value = mass * vel[2];
739 break;
740 default :
741 break;
742 }
743
744 if (!min_found) {
745 min_val = value;
746 min_sd = sd;
747 min_found = true;
748 } else {
749 if (min_val > value) {
750 min_val = value;
751 min_sd = sd;
752 }
753 }
754 }
755
756 #ifdef IS_MPI
757 int worldRank;
758 MPI_Comm_rank( MPI_COMM_WORLD, &worldRank);
759
760 int my_min_found = min_found;
761 int my_max_found = max_found;
762
763 // Even if we didn't find a minimum, did someone else?
764 MPI_Allreduce(&my_min_found, &min_found, 1, MPI_INT, MPI_LOR,
765 MPI_COMM_WORLD);
766 // Even if we didn't find a maximum, did someone else?
767 MPI_Allreduce(&my_max_found, &max_found, 1, MPI_INT, MPI_LOR,
768 MPI_COMM_WORLD);
769 #endif
770
771 if (max_found && min_found) {
772
773 #ifdef IS_MPI
774 struct {
775 RealType val;
776 int rank;
777 } max_vals, min_vals;
778
779 if (my_min_found) {
780 min_vals.val = min_val;
781 } else {
782 min_vals.val = HONKING_LARGE_VALUE;
783 }
784 min_vals.rank = worldRank;
785
786 // Who had the minimum?
787 MPI_Allreduce(&min_vals, &min_vals,
788 1, MPI_REALTYPE_INT, MPI_MINLOC, MPI_COMM_WORLD);
789 min_val = min_vals.val;
790
791 if (my_max_found) {
792 max_vals.val = max_val;
793 } else {
794 max_vals.val = -HONKING_LARGE_VALUE;
795 }
796 max_vals.rank = worldRank;
797
798 // Who had the maximum?
799 MPI_Allreduce(&max_vals, &max_vals,
800 1, MPI_REALTYPE_INT, MPI_MAXLOC, MPI_COMM_WORLD);
801 max_val = max_vals.val;
802 #endif
803
804 if (min_val < max_val) {
805
806 #ifdef IS_MPI
807 if (max_vals.rank == worldRank && min_vals.rank == worldRank) {
808 // I have both maximum and minimum, so proceed like a single
809 // processor version:
810 #endif
811
812 Vector3d min_vel = min_sd->getVel();
813 Vector3d max_vel = max_sd->getVel();
814 RealType temp_vel;
815
816 switch(rnemdFluxType_) {
817 case rnemdKE :
818 min_sd->setVel(max_vel);
819 max_sd->setVel(min_vel);
820 if (min_sd->isDirectional() && max_sd->isDirectional()) {
821 Vector3d min_angMom = min_sd->getJ();
822 Vector3d max_angMom = max_sd->getJ();
823 min_sd->setJ(max_angMom);
824 max_sd->setJ(min_angMom);
825 }//angular momenta exchange enabled
826 //assumes same rigid body identity
827 break;
828 case rnemdPx :
829 temp_vel = min_vel.x();
830 min_vel.x() = max_vel.x();
831 max_vel.x() = temp_vel;
832 min_sd->setVel(min_vel);
833 max_sd->setVel(max_vel);
834 break;
835 case rnemdPy :
836 temp_vel = min_vel.y();
837 min_vel.y() = max_vel.y();
838 max_vel.y() = temp_vel;
839 min_sd->setVel(min_vel);
840 max_sd->setVel(max_vel);
841 break;
842 case rnemdPz :
843 temp_vel = min_vel.z();
844 min_vel.z() = max_vel.z();
845 max_vel.z() = temp_vel;
846 min_sd->setVel(min_vel);
847 max_sd->setVel(max_vel);
848 break;
849 default :
850 break;
851 }
852
853 #ifdef IS_MPI
854 // the rest of the cases only apply in parallel simulations:
855 } else if (max_vals.rank == worldRank) {
856 // I had the max, but not the minimum
857
858 Vector3d min_vel;
859 Vector3d max_vel = max_sd->getVel();
860 MPI_Status* status;
861
862 // point-to-point swap of the velocity vector
863 MPI_Sendrecv(max_vel.getArrayPointer(), 3, MPI_REALTYPE,
864 min_vals.rank, 0,
865 min_vel.getArrayPointer(), 3, MPI_REALTYPE,
866 min_vals.rank, 0, MPI_COMM_WORLD, status);
867
868 switch(rnemdFluxType_) {
869 case rnemdKE :
870 max_sd->setVel(min_vel);
871 //angular momenta exchange enabled
872 if (max_sd->isDirectional()) {
873 Vector3d min_angMom;
874 Vector3d max_angMom = max_sd->getJ();
875
876 // point-to-point swap of the angular momentum vector
877 MPI_Sendrecv(max_angMom.getArrayPointer(), 3,
878 MPI_REALTYPE, min_vals.rank, 1,
879 min_angMom.getArrayPointer(), 3,
880 MPI_REALTYPE, min_vals.rank, 1,
881 MPI_COMM_WORLD, status);
882
883 max_sd->setJ(min_angMom);
884 }
885 break;
886 case rnemdPx :
887 max_vel.x() = min_vel.x();
888 max_sd->setVel(max_vel);
889 break;
890 case rnemdPy :
891 max_vel.y() = min_vel.y();
892 max_sd->setVel(max_vel);
893 break;
894 case rnemdPz :
895 max_vel.z() = min_vel.z();
896 max_sd->setVel(max_vel);
897 break;
898 default :
899 break;
900 }
901 } else if (min_vals.rank == worldRank) {
902 // I had the minimum but not the maximum:
903
904 Vector3d max_vel;
905 Vector3d min_vel = min_sd->getVel();
906 MPI_Status* status;
907
908 // point-to-point swap of the velocity vector
909 MPI_Sendrecv(min_vel.getArrayPointer(), 3, MPI_REALTYPE,
910 max_vals.rank, 0,
911 max_vel.getArrayPointer(), 3, MPI_REALTYPE,
912 max_vals.rank, 0, MPI_COMM_WORLD, status);
913
914 switch(rnemdFluxType_) {
915 case rnemdKE :
916 min_sd->setVel(max_vel);
917 //angular momenta exchange enabled
918 if (min_sd->isDirectional()) {
919 Vector3d min_angMom = min_sd->getJ();
920 Vector3d max_angMom;
921
922 // point-to-point swap of the angular momentum vector
923 MPI_Sendrecv(min_angMom.getArrayPointer(), 3,
924 MPI_REALTYPE, max_vals.rank, 1,
925 max_angMom.getArrayPointer(), 3,
926 MPI_REALTYPE, max_vals.rank, 1,
927 MPI_COMM_WORLD, status);
928
929 min_sd->setJ(max_angMom);
930 }
931 break;
932 case rnemdPx :
933 min_vel.x() = max_vel.x();
934 min_sd->setVel(min_vel);
935 break;
936 case rnemdPy :
937 min_vel.y() = max_vel.y();
938 min_sd->setVel(min_vel);
939 break;
940 case rnemdPz :
941 min_vel.z() = max_vel.z();
942 min_sd->setVel(min_vel);
943 break;
944 default :
945 break;
946 }
947 }
948 #endif
949
950 switch(rnemdFluxType_) {
951 case rnemdKE:
952 kineticExchange_ += max_val - min_val;
953 break;
954 case rnemdPx:
955 momentumExchange_.x() += max_val - min_val;
956 break;
957 case rnemdPy:
958 momentumExchange_.y() += max_val - min_val;
959 break;
960 case rnemdPz:
961 momentumExchange_.z() += max_val - min_val;
962 break;
963 default:
964 break;
965 }
966 } else {
967 sprintf(painCave.errMsg,
968 "RNEMD::doSwap exchange NOT performed because min_val > max_val\n");
969 painCave.isFatal = 0;
970 painCave.severity = OPENMD_INFO;
971 simError();
972 failTrialCount_++;
973 }
974 } else {
975 sprintf(painCave.errMsg,
976 "RNEMD::doSwap exchange NOT performed because selected object\n"
977 "\twas not present in at least one of the two slabs.\n");
978 painCave.isFatal = 0;
979 painCave.severity = OPENMD_INFO;
980 simError();
981 failTrialCount_++;
982 }
983 }
984
985 void RNEMD::doNIVS(SelectionManager& smanA, SelectionManager& smanB) {
986 if (!doRNEMD_) return;
987 int selei;
988 int selej;
989
990 Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
991 RealType time = currentSnap_->getTime();
992 Mat3x3d hmat = currentSnap_->getHmat();
993
994 StuntDouble* sd;
995
996 vector<StuntDouble*> hotBin, coldBin;
997
998 RealType Phx = 0.0;
999 RealType Phy = 0.0;
1000 RealType Phz = 0.0;
1001 RealType Khx = 0.0;
1002 RealType Khy = 0.0;
1003 RealType Khz = 0.0;
1004 RealType Khw = 0.0;
1005 RealType Pcx = 0.0;
1006 RealType Pcy = 0.0;
1007 RealType Pcz = 0.0;
1008 RealType Kcx = 0.0;
1009 RealType Kcy = 0.0;
1010 RealType Kcz = 0.0;
1011 RealType Kcw = 0.0;
1012
1013 for (sd = smanA.beginSelected(selei); sd != NULL;
1014 sd = smanA.nextSelected(selei)) {
1015
1016 Vector3d pos = sd->getPos();
1017
1018 // wrap the stuntdouble's position back into the box:
1019
1020 if (usePeriodicBoundaryConditions_)
1021 currentSnap_->wrapVector(pos);
1022
1023
1024 RealType mass = sd->getMass();
1025 Vector3d vel = sd->getVel();
1026
1027 hotBin.push_back(sd);
1028 Phx += mass * vel.x();
1029 Phy += mass * vel.y();
1030 Phz += mass * vel.z();
1031 Khx += mass * vel.x() * vel.x();
1032 Khy += mass * vel.y() * vel.y();
1033 Khz += mass * vel.z() * vel.z();
1034 if (sd->isDirectional()) {
1035 Vector3d angMom = sd->getJ();
1036 Mat3x3d I = sd->getI();
1037 if (sd->isLinear()) {
1038 int i = sd->linearAxis();
1039 int j = (i + 1) % 3;
1040 int k = (i + 2) % 3;
1041 Khw += angMom[j] * angMom[j] / I(j, j) +
1042 angMom[k] * angMom[k] / I(k, k);
1043 } else {
1044 Khw += angMom[0]*angMom[0]/I(0, 0)
1045 + angMom[1]*angMom[1]/I(1, 1)
1046 + angMom[2]*angMom[2]/I(2, 2);
1047 }
1048 }
1049 }
1050 for (sd = smanB.beginSelected(selej); sd != NULL;
1051 sd = smanB.nextSelected(selej)) {
1052 Vector3d pos = sd->getPos();
1053
1054 // wrap the stuntdouble's position back into the box:
1055
1056 if (usePeriodicBoundaryConditions_)
1057 currentSnap_->wrapVector(pos);
1058
1059 RealType mass = sd->getMass();
1060 Vector3d vel = sd->getVel();
1061
1062 coldBin.push_back(sd);
1063 Pcx += mass * vel.x();
1064 Pcy += mass * vel.y();
1065 Pcz += mass * vel.z();
1066 Kcx += mass * vel.x() * vel.x();
1067 Kcy += mass * vel.y() * vel.y();
1068 Kcz += mass * vel.z() * vel.z();
1069 if (sd->isDirectional()) {
1070 Vector3d angMom = sd->getJ();
1071 Mat3x3d I = sd->getI();
1072 if (sd->isLinear()) {
1073 int i = sd->linearAxis();
1074 int j = (i + 1) % 3;
1075 int k = (i + 2) % 3;
1076 Kcw += angMom[j] * angMom[j] / I(j, j) +
1077 angMom[k] * angMom[k] / I(k, k);
1078 } else {
1079 Kcw += angMom[0]*angMom[0]/I(0, 0)
1080 + angMom[1]*angMom[1]/I(1, 1)
1081 + angMom[2]*angMom[2]/I(2, 2);
1082 }
1083 }
1084 }
1085
1086 Khx *= 0.5;
1087 Khy *= 0.5;
1088 Khz *= 0.5;
1089 Khw *= 0.5;
1090 Kcx *= 0.5;
1091 Kcy *= 0.5;
1092 Kcz *= 0.5;
1093 Kcw *= 0.5;
1094
1095 #ifdef IS_MPI
1096 MPI_Allreduce(MPI_IN_PLACE, &Phx, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1097 MPI_Allreduce(MPI_IN_PLACE, &Phy, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1098 MPI_Allreduce(MPI_IN_PLACE, &Phz, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1099 MPI_Allreduce(MPI_IN_PLACE, &Pcx, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1100 MPI_Allreduce(MPI_IN_PLACE, &Pcy, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1101 MPI_Allreduce(MPI_IN_PLACE, &Pcz, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1102
1103 MPI_Allreduce(MPI_IN_PLACE, &Khx, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1104 MPI_Allreduce(MPI_IN_PLACE, &Khy, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1105 MPI_Allreduce(MPI_IN_PLACE, &Khz, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1106 MPI_Allreduce(MPI_IN_PLACE, &Khw, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1107
1108 MPI_Allreduce(MPI_IN_PLACE, &Kcx, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1109 MPI_Allreduce(MPI_IN_PLACE, &Kcy, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1110 MPI_Allreduce(MPI_IN_PLACE, &Kcz, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1111 MPI_Allreduce(MPI_IN_PLACE, &Kcw, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1112 #endif
1113
1114 //solve coldBin coeff's first
1115 RealType px = Pcx / Phx;
1116 RealType py = Pcy / Phy;
1117 RealType pz = Pcz / Phz;
1118 RealType c, x, y, z;
1119 bool successfulScale = false;
1120 if ((rnemdFluxType_ == rnemdFullKE) ||
1121 (rnemdFluxType_ == rnemdRotKE)) {
1122 //may need sanity check Khw & Kcw > 0
1123
1124 if (rnemdFluxType_ == rnemdFullKE) {
1125 c = 1.0 - kineticTarget_ / (Kcx + Kcy + Kcz + Kcw);
1126 } else {
1127 c = 1.0 - kineticTarget_ / Kcw;
1128 }
1129
1130 if ((c > 0.81) && (c < 1.21)) {//restrict scaling coefficients
1131 c = sqrt(c);
1132
1133 RealType w = 0.0;
1134 if (rnemdFluxType_ == rnemdFullKE) {
1135 x = 1.0 + px * (1.0 - c);
1136 y = 1.0 + py * (1.0 - c);
1137 z = 1.0 + pz * (1.0 - c);
1138 /* more complicated way
1139 w = 1.0 + (Kcw - Kcw * c * c - (c * c * (Kcx + Kcy + Kcz
1140 + Khx * px * px + Khy * py * py + Khz * pz * pz)
1141 - 2.0 * c * (Khx * px * (1.0 + px) + Khy * py * (1.0 + py)
1142 + Khz * pz * (1.0 + pz)) + Khx * px * (2.0 + px)
1143 + Khy * py * (2.0 + py) + Khz * pz * (2.0 + pz)
1144 - Kcx - Kcy - Kcz)) / Khw; the following is simpler
1145 */
1146 if ((fabs(x - 1.0) < 0.1) && (fabs(y - 1.0) < 0.1) &&
1147 (fabs(z - 1.0) < 0.1)) {
1148 w = 1.0 + (kineticTarget_
1149 + Khx * (1.0 - x * x) + Khy * (1.0 - y * y)
1150 + Khz * (1.0 - z * z)) / Khw;
1151 }//no need to calculate w if x, y or z is out of range
1152 } else {
1153 w = 1.0 + kineticTarget_ / Khw;
1154 }
1155 if ((w > 0.81) && (w < 1.21)) {//restrict scaling coefficients
1156 //if w is in the right range, so should be x, y, z.
1157 vector<StuntDouble*>::iterator sdi;
1158 Vector3d vel;
1159 for (sdi = coldBin.begin(); sdi != coldBin.end(); ++sdi) {
1160 if (rnemdFluxType_ == rnemdFullKE) {
1161 vel = (*sdi)->getVel() * c;
1162 (*sdi)->setVel(vel);
1163 }
1164 if ((*sdi)->isDirectional()) {
1165 Vector3d angMom = (*sdi)->getJ() * c;
1166 (*sdi)->setJ(angMom);
1167 }
1168 }
1169 w = sqrt(w);
1170 for (sdi = hotBin.begin(); sdi != hotBin.end(); ++sdi) {
1171 if (rnemdFluxType_ == rnemdFullKE) {
1172 vel = (*sdi)->getVel();
1173 vel.x() *= x;
1174 vel.y() *= y;
1175 vel.z() *= z;
1176 (*sdi)->setVel(vel);
1177 }
1178 if ((*sdi)->isDirectional()) {
1179 Vector3d angMom = (*sdi)->getJ() * w;
1180 (*sdi)->setJ(angMom);
1181 }
1182 }
1183 successfulScale = true;
1184 kineticExchange_ += kineticTarget_;
1185 }
1186 }
1187 } else {
1188 RealType a000, a110, c0, a001, a111, b01, b11, c1;
1189 switch(rnemdFluxType_) {
1190 case rnemdKE :
1191 /* used hotBin coeff's & only scale x & y dimensions
1192 RealType px = Phx / Pcx;
1193 RealType py = Phy / Pcy;
1194 a110 = Khy;
1195 c0 = - Khx - Khy - kineticTarget_;
1196 a000 = Khx;
1197 a111 = Kcy * py * py;
1198 b11 = -2.0 * Kcy * py * (1.0 + py);
1199 c1 = Kcy * py * (2.0 + py) + Kcx * px * ( 2.0 + px) + kineticTarget_;
1200 b01 = -2.0 * Kcx * px * (1.0 + px);
1201 a001 = Kcx * px * px;
1202 */
1203 //scale all three dimensions, let c_x = c_y
1204 a000 = Kcx + Kcy;
1205 a110 = Kcz;
1206 c0 = kineticTarget_ - Kcx - Kcy - Kcz;
1207 a001 = Khx * px * px + Khy * py * py;
1208 a111 = Khz * pz * pz;
1209 b01 = -2.0 * (Khx * px * (1.0 + px) + Khy * py * (1.0 + py));
1210 b11 = -2.0 * Khz * pz * (1.0 + pz);
1211 c1 = Khx * px * (2.0 + px) + Khy * py * (2.0 + py)
1212 + Khz * pz * (2.0 + pz) - kineticTarget_;
1213 break;
1214 case rnemdPx :
1215 c = 1 - momentumTarget_.x() / Pcx;
1216 a000 = Kcy;
1217 a110 = Kcz;
1218 c0 = Kcx * c * c - Kcx - Kcy - Kcz;
1219 a001 = py * py * Khy;
1220 a111 = pz * pz * Khz;
1221 b01 = -2.0 * Khy * py * (1.0 + py);
1222 b11 = -2.0 * Khz * pz * (1.0 + pz);
1223 c1 = Khy * py * (2.0 + py) + Khz * pz * (2.0 + pz)
1224 + Khx * (fastpow(c * px - px - 1.0, 2) - 1.0);
1225 break;
1226 case rnemdPy :
1227 c = 1 - momentumTarget_.y() / Pcy;
1228 a000 = Kcx;
1229 a110 = Kcz;
1230 c0 = Kcy * c * c - Kcx - Kcy - Kcz;
1231 a001 = px * px * Khx;
1232 a111 = pz * pz * Khz;
1233 b01 = -2.0 * Khx * px * (1.0 + px);
1234 b11 = -2.0 * Khz * pz * (1.0 + pz);
1235 c1 = Khx * px * (2.0 + px) + Khz * pz * (2.0 + pz)
1236 + Khy * (fastpow(c * py - py - 1.0, 2) - 1.0);
1237 break;
1238 case rnemdPz ://we don't really do this, do we?
1239 c = 1 - momentumTarget_.z() / Pcz;
1240 a000 = Kcx;
1241 a110 = Kcy;
1242 c0 = Kcz * c * c - Kcx - Kcy - Kcz;
1243 a001 = px * px * Khx;
1244 a111 = py * py * Khy;
1245 b01 = -2.0 * Khx * px * (1.0 + px);
1246 b11 = -2.0 * Khy * py * (1.0 + py);
1247 c1 = Khx * px * (2.0 + px) + Khy * py * (2.0 + py)
1248 + Khz * (fastpow(c * pz - pz - 1.0, 2) - 1.0);
1249 break;
1250 default :
1251 break;
1252 }
1253
1254 RealType v1 = a000 * a111 - a001 * a110;
1255 RealType v2 = a000 * b01;
1256 RealType v3 = a000 * b11;
1257 RealType v4 = a000 * c1 - a001 * c0;
1258 RealType v8 = a110 * b01;
1259 RealType v10 = - b01 * c0;
1260
1261 RealType u0 = v2 * v10 - v4 * v4;
1262 RealType u1 = -2.0 * v3 * v4;
1263 RealType u2 = -v2 * v8 - v3 * v3 - 2.0 * v1 * v4;
1264 RealType u3 = -2.0 * v1 * v3;
1265 RealType u4 = - v1 * v1;
1266 //rescale coefficients
1267 RealType maxAbs = fabs(u0);
1268 if (maxAbs < fabs(u1)) maxAbs = fabs(u1);
1269 if (maxAbs < fabs(u2)) maxAbs = fabs(u2);
1270 if (maxAbs < fabs(u3)) maxAbs = fabs(u3);
1271 if (maxAbs < fabs(u4)) maxAbs = fabs(u4);
1272 u0 /= maxAbs;
1273 u1 /= maxAbs;
1274 u2 /= maxAbs;
1275 u3 /= maxAbs;
1276 u4 /= maxAbs;
1277 //max_element(start, end) is also available.
1278 Polynomial<RealType> poly; //same as DoublePolynomial poly;
1279 poly.setCoefficient(4, u4);
1280 poly.setCoefficient(3, u3);
1281 poly.setCoefficient(2, u2);
1282 poly.setCoefficient(1, u1);
1283 poly.setCoefficient(0, u0);
1284 vector<RealType> realRoots = poly.FindRealRoots();
1285
1286 vector<RealType>::iterator ri;
1287 RealType r1, r2, alpha0;
1288 vector<pair<RealType,RealType> > rps;
1289 for (ri = realRoots.begin(); ri !=realRoots.end(); ++ri) {
1290 r2 = *ri;
1291 //check if FindRealRoots() give the right answer
1292 if ( fabs(u0 + r2 * (u1 + r2 * (u2 + r2 * (u3 + r2 * u4)))) > 1e-6 ) {
1293 sprintf(painCave.errMsg,
1294 "RNEMD Warning: polynomial solve seems to have an error!");
1295 painCave.isFatal = 0;
1296 simError();
1297 failRootCount_++;
1298 }
1299 //might not be useful w/o rescaling coefficients
1300 alpha0 = -c0 - a110 * r2 * r2;
1301 if (alpha0 >= 0.0) {
1302 r1 = sqrt(alpha0 / a000);
1303 if (fabs(c1 + r1 * (b01 + r1 * a001) + r2 * (b11 + r2 * a111))
1304 < 1e-6)
1305 { rps.push_back(make_pair(r1, r2)); }
1306 if (r1 > 1e-6) { //r1 non-negative
1307 r1 = -r1;
1308 if (fabs(c1 + r1 * (b01 + r1 * a001) + r2 * (b11 + r2 * a111))
1309 < 1e-6)
1310 { rps.push_back(make_pair(r1, r2)); }
1311 }
1312 }
1313 }
1314 // Consider combining together the solving pair part w/ the searching
1315 // best solution part so that we don't need the pairs vector
1316 if (!rps.empty()) {
1317 RealType smallestDiff = HONKING_LARGE_VALUE;
1318 RealType diff;
1319 pair<RealType,RealType> bestPair = make_pair(1.0, 1.0);
1320 vector<pair<RealType,RealType> >::iterator rpi;
1321 for (rpi = rps.begin(); rpi != rps.end(); ++rpi) {
1322 r1 = (*rpi).first;
1323 r2 = (*rpi).second;
1324 switch(rnemdFluxType_) {
1325 case rnemdKE :
1326 diff = fastpow(1.0 - r1, 2) + fastpow(1.0 - r2, 2)
1327 + fastpow(r1 * r1 / r2 / r2 - Kcz/Kcx, 2)
1328 + fastpow(r1 * r1 / r2 / r2 - Kcz/Kcy, 2);
1329 break;
1330 case rnemdPx :
1331 diff = fastpow(1.0 - r1, 2) + fastpow(1.0 - r2, 2)
1332 + fastpow(r1 * r1 / r2 / r2 - Kcz/Kcy, 2);
1333 break;
1334 case rnemdPy :
1335 diff = fastpow(1.0 - r1, 2) + fastpow(1.0 - r2, 2)
1336 + fastpow(r1 * r1 / r2 / r2 - Kcz/Kcx, 2);
1337 break;
1338 case rnemdPz :
1339 diff = fastpow(1.0 - r1, 2) + fastpow(1.0 - r2, 2)
1340 + fastpow(r1 * r1 / r2 / r2 - Kcy/Kcx, 2);
1341 default :
1342 break;
1343 }
1344 if (diff < smallestDiff) {
1345 smallestDiff = diff;
1346 bestPair = *rpi;
1347 }
1348 }
1349 #ifdef IS_MPI
1350 if (worldRank == 0) {
1351 #endif
1352 // sprintf(painCave.errMsg,
1353 // "RNEMD: roots r1= %lf\tr2 = %lf\n",
1354 // bestPair.first, bestPair.second);
1355 // painCave.isFatal = 0;
1356 // painCave.severity = OPENMD_INFO;
1357 // simError();
1358 #ifdef IS_MPI
1359 }
1360 #endif
1361
1362 switch(rnemdFluxType_) {
1363 case rnemdKE :
1364 x = bestPair.first;
1365 y = bestPair.first;
1366 z = bestPair.second;
1367 break;
1368 case rnemdPx :
1369 x = c;
1370 y = bestPair.first;
1371 z = bestPair.second;
1372 break;
1373 case rnemdPy :
1374 x = bestPair.first;
1375 y = c;
1376 z = bestPair.second;
1377 break;
1378 case rnemdPz :
1379 x = bestPair.first;
1380 y = bestPair.second;
1381 z = c;
1382 break;
1383 default :
1384 break;
1385 }
1386 vector<StuntDouble*>::iterator sdi;
1387 Vector3d vel;
1388 for (sdi = coldBin.begin(); sdi != coldBin.end(); ++sdi) {
1389 vel = (*sdi)->getVel();
1390 vel.x() *= x;
1391 vel.y() *= y;
1392 vel.z() *= z;
1393 (*sdi)->setVel(vel);
1394 }
1395 //convert to hotBin coefficient
1396 x = 1.0 + px * (1.0 - x);
1397 y = 1.0 + py * (1.0 - y);
1398 z = 1.0 + pz * (1.0 - z);
1399 for (sdi = hotBin.begin(); sdi != hotBin.end(); ++sdi) {
1400 vel = (*sdi)->getVel();
1401 vel.x() *= x;
1402 vel.y() *= y;
1403 vel.z() *= z;
1404 (*sdi)->setVel(vel);
1405 }
1406 successfulScale = true;
1407 switch(rnemdFluxType_) {
1408 case rnemdKE :
1409 kineticExchange_ += kineticTarget_;
1410 break;
1411 case rnemdPx :
1412 case rnemdPy :
1413 case rnemdPz :
1414 momentumExchange_ += momentumTarget_;
1415 break;
1416 default :
1417 break;
1418 }
1419 }
1420 }
1421 if (successfulScale != true) {
1422 sprintf(painCave.errMsg,
1423 "RNEMD::doNIVS exchange NOT performed - roots that solve\n"
1424 "\tthe constraint equations may not exist or there may be\n"
1425 "\tno selected objects in one or both slabs.\n");
1426 painCave.isFatal = 0;
1427 painCave.severity = OPENMD_INFO;
1428 simError();
1429 failTrialCount_++;
1430 }
1431 }
1432
1433 void RNEMD::doVSS(SelectionManager& smanA, SelectionManager& smanB) {
1434 if (!doRNEMD_) return;
1435 int selei;
1436 int selej;
1437
1438 Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
1439 RealType time = currentSnap_->getTime();
1440 Mat3x3d hmat = currentSnap_->getHmat();
1441
1442 StuntDouble* sd;
1443
1444 vector<StuntDouble*> hotBin, coldBin;
1445
1446 Vector3d Ph(V3Zero);
1447 Vector3d Lh(V3Zero);
1448 RealType Mh = 0.0;
1449 Mat3x3d Ih(0.0);
1450 RealType Kh = 0.0;
1451 Vector3d Pc(V3Zero);
1452 Vector3d Lc(V3Zero);
1453 RealType Mc = 0.0;
1454 Mat3x3d Ic(0.0);
1455 RealType Kc = 0.0;
1456
1457 // Constraints can be on only the linear or angular momentum, but
1458 // not both. Usually, the user will specify which they want, but
1459 // in case they don't, the use of periodic boundaries should make
1460 // the choice for us.
1461 bool doLinearPart = false;
1462 bool doAngularPart = false;
1463
1464 switch (rnemdFluxType_) {
1465 case rnemdPx:
1466 case rnemdPy:
1467 case rnemdPz:
1468 case rnemdPvector:
1469 case rnemdKePx:
1470 case rnemdKePy:
1471 case rnemdKePvector:
1472 doLinearPart = true;
1473 break;
1474 case rnemdLx:
1475 case rnemdLy:
1476 case rnemdLz:
1477 case rnemdLvector:
1478 case rnemdKeLx:
1479 case rnemdKeLy:
1480 case rnemdKeLz:
1481 case rnemdKeLvector:
1482 doAngularPart = true;
1483 break;
1484 case rnemdKE:
1485 case rnemdRotKE:
1486 case rnemdFullKE:
1487 default:
1488 if (usePeriodicBoundaryConditions_)
1489 doLinearPart = true;
1490 else
1491 doAngularPart = true;
1492 break;
1493 }
1494
1495 for (sd = smanA.beginSelected(selei); sd != NULL;
1496 sd = smanA.nextSelected(selei)) {
1497
1498 Vector3d pos = sd->getPos();
1499
1500 // wrap the stuntdouble's position back into the box:
1501
1502 if (usePeriodicBoundaryConditions_)
1503 currentSnap_->wrapVector(pos);
1504
1505 RealType mass = sd->getMass();
1506 Vector3d vel = sd->getVel();
1507 Vector3d rPos = sd->getPos() - coordinateOrigin_;
1508 RealType r2;
1509
1510 hotBin.push_back(sd);
1511 Ph += mass * vel;
1512 Mh += mass;
1513 Kh += mass * vel.lengthSquare();
1514 Lh += mass * cross(rPos, vel);
1515 Ih -= outProduct(rPos, rPos) * mass;
1516 r2 = rPos.lengthSquare();
1517 Ih(0, 0) += mass * r2;
1518 Ih(1, 1) += mass * r2;
1519 Ih(2, 2) += mass * r2;
1520
1521 if (rnemdFluxType_ == rnemdFullKE) {
1522 if (sd->isDirectional()) {
1523 Vector3d angMom = sd->getJ();
1524 Mat3x3d I = sd->getI();
1525 if (sd->isLinear()) {
1526 int i = sd->linearAxis();
1527 int j = (i + 1) % 3;
1528 int k = (i + 2) % 3;
1529 Kh += angMom[j] * angMom[j] / I(j, j) +
1530 angMom[k] * angMom[k] / I(k, k);
1531 } else {
1532 Kh += angMom[0] * angMom[0] / I(0, 0) +
1533 angMom[1] * angMom[1] / I(1, 1) +
1534 angMom[2] * angMom[2] / I(2, 2);
1535 }
1536 }
1537 }
1538 }
1539 for (sd = smanB.beginSelected(selej); sd != NULL;
1540 sd = smanB.nextSelected(selej)) {
1541
1542 Vector3d pos = sd->getPos();
1543
1544 // wrap the stuntdouble's position back into the box:
1545
1546 if (usePeriodicBoundaryConditions_)
1547 currentSnap_->wrapVector(pos);
1548
1549 RealType mass = sd->getMass();
1550 Vector3d vel = sd->getVel();
1551 Vector3d rPos = sd->getPos() - coordinateOrigin_;
1552 RealType r2;
1553
1554 coldBin.push_back(sd);
1555 Pc += mass * vel;
1556 Mc += mass;
1557 Kc += mass * vel.lengthSquare();
1558 Lc += mass * cross(rPos, vel);
1559 Ic -= outProduct(rPos, rPos) * mass;
1560 r2 = rPos.lengthSquare();
1561 Ic(0, 0) += mass * r2;
1562 Ic(1, 1) += mass * r2;
1563 Ic(2, 2) += mass * r2;
1564
1565 if (rnemdFluxType_ == rnemdFullKE) {
1566 if (sd->isDirectional()) {
1567 Vector3d angMom = sd->getJ();
1568 Mat3x3d I = sd->getI();
1569 if (sd->isLinear()) {
1570 int i = sd->linearAxis();
1571 int j = (i + 1) % 3;
1572 int k = (i + 2) % 3;
1573 Kc += angMom[j] * angMom[j] / I(j, j) +
1574 angMom[k] * angMom[k] / I(k, k);
1575 } else {
1576 Kc += angMom[0] * angMom[0] / I(0, 0) +
1577 angMom[1] * angMom[1] / I(1, 1) +
1578 angMom[2] * angMom[2] / I(2, 2);
1579 }
1580 }
1581 }
1582 }
1583
1584 Kh *= 0.5;
1585 Kc *= 0.5;
1586
1587 #ifdef IS_MPI
1588 MPI_Allreduce(MPI_IN_PLACE, &Ph[0], 3, MPI_REALTYPE, MPI_SUM,
1589 MPI_COMM_WORLD);
1590 MPI_Allreduce(MPI_IN_PLACE, &Pc[0], 3, MPI_REALTYPE, MPI_SUM,
1591 MPI_COMM_WORLD);
1592 MPI_Allreduce(MPI_IN_PLACE, &Lh[0], 3, MPI_REALTYPE, MPI_SUM,
1593 MPI_COMM_WORLD);
1594 MPI_Allreduce(MPI_IN_PLACE, &Lc[0], 3, MPI_REALTYPE, MPI_SUM,
1595 MPI_COMM_WORLD);
1596 MPI_Allreduce(MPI_IN_PLACE, &Mh, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1597 MPI_Allreduce(MPI_IN_PLACE, &Kh, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1598 MPI_Allreduce(MPI_IN_PLACE, &Mc, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1599 MPI_Allreduce(MPI_IN_PLACE, &Kc, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1600 MPI_Allreduce(MPI_IN_PLACE, Ih.getArrayPointer(), 9,
1601 MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1602 MPI_Allreduce(MPI_IN_PLACE, Ic.getArrayPointer(), 9,
1603 MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1604 #endif
1605
1606
1607 Vector3d ac, acrec, bc, bcrec;
1608 Vector3d ah, ahrec, bh, bhrec;
1609
1610 bool successfulExchange = false;
1611 if ((Mh > 0.0) && (Mc > 0.0)) {//both slabs are not empty
1612 Vector3d vc = Pc / Mc;
1613 ac = -momentumTarget_ / Mc + vc;
1614 acrec = -momentumTarget_ / Mc;
1615
1616 // We now need the inverse of the inertia tensor to calculate the
1617 // angular velocity of the cold slab;
1618 Mat3x3d Ici = Ic.inverse();
1619 Vector3d omegac = Ici * Lc;
1620 bc = -(Ici * angularMomentumTarget_) + omegac;
1621 bcrec = bc - omegac;
1622
1623 RealType cNumerator = Kc - kineticTarget_;
1624 if (doLinearPart)
1625 cNumerator -= 0.5 * Mc * ac.lengthSquare();
1626
1627 if (doAngularPart)
1628 cNumerator -= 0.5 * ( dot(bc, Ic * bc));
1629
1630 if (cNumerator > 0.0) {
1631
1632 RealType cDenominator = Kc;
1633
1634 if (doLinearPart)
1635 cDenominator -= 0.5 * Mc * vc.lengthSquare();
1636
1637 if (doAngularPart)
1638 cDenominator -= 0.5*(dot(omegac, Ic * omegac));
1639
1640 if (cDenominator > 0.0) {
1641 RealType c = sqrt(cNumerator / cDenominator);
1642 if ((c > 0.9) && (c < 1.1)) {//restrict scaling coefficients
1643
1644 Vector3d vh = Ph / Mh;
1645 ah = momentumTarget_ / Mh + vh;
1646 ahrec = momentumTarget_ / Mh;
1647
1648 // We now need the inverse of the inertia tensor to
1649 // calculate the angular velocity of the hot slab;
1650 Mat3x3d Ihi = Ih.inverse();
1651 Vector3d omegah = Ihi * Lh;
1652 bh = (Ihi * angularMomentumTarget_) + omegah;
1653 bhrec = bh - omegah;
1654
1655 RealType hNumerator = Kh + kineticTarget_;
1656 if (doLinearPart)
1657 hNumerator -= 0.5 * Mh * ah.lengthSquare();
1658
1659 if (doAngularPart)
1660 hNumerator -= 0.5 * ( dot(bh, Ih * bh));
1661
1662 if (hNumerator > 0.0) {
1663
1664 RealType hDenominator = Kh;
1665 if (doLinearPart)
1666 hDenominator -= 0.5 * Mh * vh.lengthSquare();
1667 if (doAngularPart)
1668 hDenominator -= 0.5*(dot(omegah, Ih * omegah));
1669
1670 if (hDenominator > 0.0) {
1671 RealType h = sqrt(hNumerator / hDenominator);
1672 if ((h > 0.9) && (h < 1.1)) {
1673
1674 vector<StuntDouble*>::iterator sdi;
1675 Vector3d vel;
1676 Vector3d rPos;
1677
1678 for (sdi = coldBin.begin(); sdi != coldBin.end(); ++sdi) {
1679 //vel = (*sdi)->getVel();
1680 rPos = (*sdi)->getPos() - coordinateOrigin_;
1681 if (doLinearPart)
1682 vel = ((*sdi)->getVel() - vc) * c + ac;
1683 if (doAngularPart)
1684 vel = ((*sdi)->getVel() - cross(omegac, rPos)) * c + cross(bc, rPos);
1685
1686 (*sdi)->setVel(vel);
1687 if (rnemdFluxType_ == rnemdFullKE) {
1688 if ((*sdi)->isDirectional()) {
1689 Vector3d angMom = (*sdi)->getJ() * c;
1690 (*sdi)->setJ(angMom);
1691 }
1692 }
1693 }
1694 for (sdi = hotBin.begin(); sdi != hotBin.end(); ++sdi) {
1695 //vel = (*sdi)->getVel();
1696 rPos = (*sdi)->getPos() - coordinateOrigin_;
1697 if (doLinearPart)
1698 vel = ((*sdi)->getVel() - vh) * h + ah;
1699 if (doAngularPart)
1700 vel = ((*sdi)->getVel() - cross(omegah, rPos)) * h + cross(bh, rPos);
1701
1702 (*sdi)->setVel(vel);
1703 if (rnemdFluxType_ == rnemdFullKE) {
1704 if ((*sdi)->isDirectional()) {
1705 Vector3d angMom = (*sdi)->getJ() * h;
1706 (*sdi)->setJ(angMom);
1707 }
1708 }
1709 }
1710 successfulExchange = true;
1711 kineticExchange_ += kineticTarget_;
1712 momentumExchange_ += momentumTarget_;
1713 angularMomentumExchange_ += angularMomentumTarget_;
1714 }
1715 }
1716 }
1717 }
1718 }
1719 }
1720 }
1721 if (successfulExchange != true) {
1722 sprintf(painCave.errMsg,
1723 "RNEMD::doVSS exchange NOT performed - roots that solve\n"
1724 "\tthe constraint equations may not exist or there may be\n"
1725 "\tno selected objects in one or both slabs.\n");
1726 painCave.isFatal = 0;
1727 painCave.severity = OPENMD_INFO;
1728 simError();
1729 failTrialCount_++;
1730 }
1731 }
1732
1733 RealType RNEMD::getDividingArea() {
1734
1735 if (hasDividingArea_) return dividingArea_;
1736
1737 RealType areaA, areaB;
1738 Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot();
1739
1740 if (hasSelectionA_) {
1741
1742 if (evaluatorA_.hasSurfaceArea())
1743 areaA = evaluatorA_.getSurfaceArea();
1744 else {
1745
1746 cerr << "selection A did not have surface area, recomputing\n";
1747 int isd;
1748 StuntDouble* sd;
1749 vector<StuntDouble*> aSites;
1750 seleManA_.setSelectionSet(evaluatorA_.evaluate());
1751 for (sd = seleManA_.beginSelected(isd); sd != NULL;
1752 sd = seleManA_.nextSelected(isd)) {
1753 aSites.push_back(sd);
1754 }
1755 #if defined(HAVE_QHULL)
1756 ConvexHull* surfaceMeshA = new ConvexHull();
1757 surfaceMeshA->computeHull(aSites);
1758 areaA = surfaceMeshA->getArea();
1759 delete surfaceMeshA;
1760 #else
1761 sprintf( painCave.errMsg,
1762 "RNEMD::getDividingArea : Hull calculation is not possible\n"
1763 "\twithout libqhull. Please rebuild OpenMD with qhull enabled.");
1764 painCave.severity = OPENMD_ERROR;
1765 painCave.isFatal = 1;
1766 simError();
1767 #endif
1768 }
1769
1770 } else {
1771 if (usePeriodicBoundaryConditions_) {
1772 // in periodic boundaries, the surface area is twice the x-y
1773 // area of the current box:
1774 areaA = 2.0 * snap->getXYarea();
1775 } else {
1776 // in non-periodic simulations, without explicitly setting
1777 // selections, the sphere radius sets the surface area of the
1778 // dividing surface:
1779 areaA = 4.0 * M_PI * pow(sphereARadius_, 2);
1780 }
1781 }
1782
1783 if (hasSelectionB_) {
1784 if (evaluatorB_.hasSurfaceArea())
1785 areaB = evaluatorB_.getSurfaceArea();
1786 else {
1787 cerr << "selection B did not have surface area, recomputing\n";
1788
1789 int isd;
1790 StuntDouble* sd;
1791 vector<StuntDouble*> bSites;
1792 seleManB_.setSelectionSet(evaluatorB_.evaluate());
1793 for (sd = seleManB_.beginSelected(isd); sd != NULL;
1794 sd = seleManB_.nextSelected(isd)) {
1795 bSites.push_back(sd);
1796 }
1797
1798 #if defined(HAVE_QHULL)
1799 ConvexHull* surfaceMeshB = new ConvexHull();
1800 surfaceMeshB->computeHull(bSites);
1801 areaB = surfaceMeshB->getArea();
1802 delete surfaceMeshB;
1803 #else
1804 sprintf( painCave.errMsg,
1805 "RNEMD::getDividingArea : Hull calculation is not possible\n"
1806 "\twithout libqhull. Please rebuild OpenMD with qhull enabled.");
1807 painCave.severity = OPENMD_ERROR;
1808 painCave.isFatal = 1;
1809 simError();
1810 #endif
1811 }
1812
1813 } else {
1814 if (usePeriodicBoundaryConditions_) {
1815 // in periodic boundaries, the surface area is twice the x-y
1816 // area of the current box:
1817 areaB = 2.0 * snap->getXYarea();
1818 } else {
1819 // in non-periodic simulations, without explicitly setting
1820 // selections, but if a sphereBradius has been set, just use that:
1821 areaB = 4.0 * M_PI * pow(sphereBRadius_, 2);
1822 }
1823 }
1824
1825 dividingArea_ = min(areaA, areaB);
1826 hasDividingArea_ = true;
1827 return dividingArea_;
1828 }
1829
1830 void RNEMD::doRNEMD() {
1831 if (!doRNEMD_) return;
1832 trialCount_++;
1833
1834 // object evaluator:
1835 evaluator_.loadScriptString(rnemdObjectSelection_);
1836 seleMan_.setSelectionSet(evaluator_.evaluate());
1837
1838 evaluatorA_.loadScriptString(selectionA_);
1839 evaluatorB_.loadScriptString(selectionB_);
1840
1841 seleManA_.setSelectionSet(evaluatorA_.evaluate());
1842 seleManB_.setSelectionSet(evaluatorB_.evaluate());
1843
1844 commonA_ = seleManA_ & seleMan_;
1845 commonB_ = seleManB_ & seleMan_;
1846
1847 // Target exchange quantities (in each exchange) = dividingArea * dt * flux
1848 // dt = exchange time interval
1849 // flux = target flux
1850 // dividingArea = smallest dividing surface between the two regions
1851
1852 hasDividingArea_ = false;
1853 RealType area = getDividingArea();
1854
1855 kineticTarget_ = kineticFlux_ * exchangeTime_ * area;
1856 momentumTarget_ = momentumFluxVector_ * exchangeTime_ * area;
1857 angularMomentumTarget_ = angularMomentumFluxVector_ * exchangeTime_ * area;
1858
1859 switch(rnemdMethod_) {
1860 case rnemdSwap:
1861 doSwap(commonA_, commonB_);
1862 break;
1863 case rnemdNIVS:
1864 doNIVS(commonA_, commonB_);
1865 break;
1866 case rnemdVSS:
1867 doVSS(commonA_, commonB_);
1868 break;
1869 case rnemdUnkownMethod:
1870 default :
1871 break;
1872 }
1873 }
1874
1875 void RNEMD::collectData() {
1876 if (!doRNEMD_) return;
1877 Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
1878
1879 // collectData can be called more frequently than the doRNEMD, so use the
1880 // computed area from the last exchange time:
1881 RealType area = getDividingArea();
1882 areaAccumulator_->add(area);
1883 Mat3x3d hmat = currentSnap_->getHmat();
1884 Vector3d u = angularMomentumFluxVector_;
1885 u.normalize();
1886
1887 seleMan_.setSelectionSet(evaluator_.evaluate());
1888
1889 int selei(0);
1890 StuntDouble* sd;
1891 int binNo;
1892 RealType mass;
1893 Vector3d vel;
1894 Vector3d rPos;
1895 RealType KE;
1896 Vector3d L;
1897 Mat3x3d I;
1898 RealType r2;
1899
1900 vector<RealType> binMass(nBins_, 0.0);
1901 vector<Vector3d> binP(nBins_, V3Zero);
1902 vector<RealType> binOmega(nBins_, 0.0);
1903 vector<Vector3d> binL(nBins_, V3Zero);
1904 vector<Mat3x3d> binI(nBins_);
1905 vector<RealType> binKE(nBins_, 0.0);
1906 vector<int> binDOF(nBins_, 0);
1907 vector<int> binCount(nBins_, 0);
1908
1909 // alternative approach, track all molecules instead of only those
1910 // selected for scaling/swapping:
1911 /*
1912 SimInfo::MoleculeIterator miter;
1913 vector<StuntDouble*>::iterator iiter;
1914 Molecule* mol;
1915 StuntDouble* sd;
1916 for (mol = info_->beginMolecule(miter); mol != NULL;
1917 mol = info_->nextMolecule(miter))
1918 sd is essentially sd
1919 for (sd = mol->beginIntegrableObject(iiter);
1920 sd != NULL;
1921 sd = mol->nextIntegrableObject(iiter))
1922 */
1923
1924 for (sd = seleMan_.beginSelected(selei); sd != NULL;
1925 sd = seleMan_.nextSelected(selei)) {
1926
1927 Vector3d pos = sd->getPos();
1928
1929 // wrap the stuntdouble's position back into the box:
1930
1931 if (usePeriodicBoundaryConditions_) {
1932 currentSnap_->wrapVector(pos);
1933 // which bin is this stuntdouble in?
1934 // wrapped positions are in the range [-0.5*hmat(2,2), +0.5*hmat(2,2)]
1935 // Shift molecules by half a box to have bins start at 0
1936 // The modulo operator is used to wrap the case when we are
1937 // beyond the end of the bins back to the beginning.
1938 binNo = int(nBins_ * (pos.z() / hmat(2,2) + 0.5)) % nBins_;
1939 } else {
1940 Vector3d rPos = pos - coordinateOrigin_;
1941 binNo = int(rPos.length() / binWidth_);
1942 }
1943
1944 mass = sd->getMass();
1945 vel = sd->getVel();
1946 rPos = sd->getPos() - coordinateOrigin_;
1947 KE = 0.5 * mass * vel.lengthSquare();
1948 L = mass * cross(rPos, vel);
1949 I = outProduct(rPos, rPos) * mass;
1950 r2 = rPos.lengthSquare();
1951 I(0, 0) += mass * r2;
1952 I(1, 1) += mass * r2;
1953 I(2, 2) += mass * r2;
1954
1955 // Project the relative position onto a plane perpendicular to
1956 // the angularMomentumFluxVector:
1957 // Vector3d rProj = rPos - dot(rPos, u) * u;
1958 // Project the velocity onto a plane perpendicular to the
1959 // angularMomentumFluxVector:
1960 // Vector3d vProj = vel - dot(vel, u) * u;
1961 // Compute angular velocity vector (should be nearly parallel to
1962 // angularMomentumFluxVector
1963 // Vector3d aVel = cross(rProj, vProj);
1964
1965 if (binNo >= 0 && binNo < nBins_) {
1966 binCount[binNo]++;
1967 binMass[binNo] += mass;
1968 binP[binNo] += mass*vel;
1969 binKE[binNo] += KE;
1970 binI[binNo] += I;
1971 binL[binNo] += L;
1972 binDOF[binNo] += 3;
1973
1974 if (sd->isDirectional()) {
1975 Vector3d angMom = sd->getJ();
1976 Mat3x3d Ia = sd->getI();
1977 if (sd->isLinear()) {
1978 int i = sd->linearAxis();
1979 int j = (i + 1) % 3;
1980 int k = (i + 2) % 3;
1981 binKE[binNo] += 0.5 * (angMom[j] * angMom[j] / Ia(j, j) +
1982 angMom[k] * angMom[k] / Ia(k, k));
1983 binDOF[binNo] += 2;
1984 } else {
1985 binKE[binNo] += 0.5 * (angMom[0] * angMom[0] / Ia(0, 0) +
1986 angMom[1] * angMom[1] / Ia(1, 1) +
1987 angMom[2] * angMom[2] / Ia(2, 2));
1988 binDOF[binNo] += 3;
1989 }
1990 }
1991 }
1992 }
1993
1994 #ifdef IS_MPI
1995
1996 for (int i = 0; i < nBins_; i++) {
1997
1998 MPI_Allreduce(MPI_IN_PLACE, &binCount[i],
1999 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
2000 MPI_Allreduce(MPI_IN_PLACE, &binMass[i],
2001 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
2002 MPI_Allreduce(MPI_IN_PLACE, &binP[i],
2003 3, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
2004 MPI_Allreduce(MPI_IN_PLACE, &binL[i],
2005 3, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
2006 MPI_Allreduce(MPI_IN_PLACE, &binI[i],
2007 9, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
2008 MPI_Allreduce(MPI_IN_PLACE, &binKE[i],
2009 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
2010 MPI_Allreduce(MPI_IN_PLACE, &binDOF[i],
2011 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
2012 //MPI_Allreduce(MPI_IN_PLACE, &binOmega[i],
2013 // 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
2014 }
2015
2016 #endif
2017
2018 Vector3d omega;
2019 RealType den;
2020 RealType temp;
2021 RealType z;
2022 RealType r;
2023 for (int i = 0; i < nBins_; i++) {
2024 if (usePeriodicBoundaryConditions_) {
2025 z = (((RealType)i + 0.5) / (RealType)nBins_) * hmat(2,2);
2026 den = binMass[i] * nBins_ * PhysicalConstants::densityConvert
2027 / currentSnap_->getVolume() ;
2028 } else {
2029 r = (((RealType)i + 0.5) * binWidth_);
2030 RealType rinner = (RealType)i * binWidth_;
2031 RealType router = (RealType)(i+1) * binWidth_;
2032 den = binMass[i] * 3.0 * PhysicalConstants::densityConvert
2033 / (4.0 * M_PI * (pow(router,3) - pow(rinner,3)));
2034 }
2035 vel = binP[i] / binMass[i];
2036
2037 omega = binI[i].inverse() * binL[i];
2038
2039 // omega = binOmega[i] / binCount[i];
2040
2041 if (binCount[i] > 0) {
2042 // only add values if there are things to add
2043 temp = 2.0 * binKE[i] / (binDOF[i] * PhysicalConstants::kb *
2044 PhysicalConstants::energyConvert);
2045
2046 for (unsigned int j = 0; j < outputMask_.size(); ++j) {
2047 if(outputMask_[j]) {
2048 switch(j) {
2049 case Z:
2050 dynamic_cast<Accumulator *>(data_[j].accumulator[i])->add(z);
2051 break;
2052 case R:
2053 dynamic_cast<Accumulator *>(data_[j].accumulator[i])->add(r);
2054 break;
2055 case TEMPERATURE:
2056 dynamic_cast<Accumulator *>(data_[j].accumulator[i])->add(temp);
2057 break;
2058 case VELOCITY:
2059 dynamic_cast<VectorAccumulator *>(data_[j].accumulator[i])->add(vel);
2060 break;
2061 case ANGULARVELOCITY:
2062 dynamic_cast<VectorAccumulator *>(data_[j].accumulator[i])->add(omega);
2063 break;
2064 case DENSITY:
2065 dynamic_cast<Accumulator *>(data_[j].accumulator[i])->add(den);
2066 break;
2067 }
2068 }
2069 }
2070 }
2071 }
2072 hasData_ = true;
2073 }
2074
2075 void RNEMD::getStarted() {
2076 if (!doRNEMD_) return;
2077 hasDividingArea_ = false;
2078 collectData();
2079 writeOutputFile();
2080 }
2081
2082 void RNEMD::parseOutputFileFormat(const std::string& format) {
2083 if (!doRNEMD_) return;
2084 StringTokenizer tokenizer(format, " ,;|\t\n\r");
2085
2086 while(tokenizer.hasMoreTokens()) {
2087 std::string token(tokenizer.nextToken());
2088 toUpper(token);
2089 OutputMapType::iterator i = outputMap_.find(token);
2090 if (i != outputMap_.end()) {
2091 outputMask_.set(i->second);
2092 } else {
2093 sprintf( painCave.errMsg,
2094 "RNEMD::parseOutputFileFormat: %s is not a recognized\n"
2095 "\toutputFileFormat keyword.\n", token.c_str() );
2096 painCave.isFatal = 0;
2097 painCave.severity = OPENMD_ERROR;
2098 simError();
2099 }
2100 }
2101 }
2102
2103 void RNEMD::writeOutputFile() {
2104 if (!doRNEMD_) return;
2105 if (!hasData_) return;
2106
2107 #ifdef IS_MPI
2108 // If we're the root node, should we print out the results
2109 int worldRank;
2110 MPI_Comm_rank( MPI_COMM_WORLD, &worldRank);
2111
2112 if (worldRank == 0) {
2113 #endif
2114 rnemdFile_.open(rnemdFileName_.c_str(), std::ios::out | std::ios::trunc );
2115
2116 if( !rnemdFile_ ){
2117 sprintf( painCave.errMsg,
2118 "Could not open \"%s\" for RNEMD output.\n",
2119 rnemdFileName_.c_str());
2120 painCave.isFatal = 1;
2121 simError();
2122 }
2123
2124 Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
2125
2126 RealType time = currentSnap_->getTime();
2127 RealType avgArea;
2128 areaAccumulator_->getAverage(avgArea);
2129
2130 RealType Jz(0.0);
2131 Vector3d JzP(V3Zero);
2132 Vector3d JzL(V3Zero);
2133 if (time >= info_->getSimParams()->getDt()) {
2134 Jz = kineticExchange_ / (time * avgArea)
2135 / PhysicalConstants::energyConvert;
2136 JzP = momentumExchange_ / (time * avgArea);
2137 JzL = angularMomentumExchange_ / (time * avgArea);
2138 }
2139
2140 rnemdFile_ << "#######################################################\n";
2141 rnemdFile_ << "# RNEMD {\n";
2142
2143 map<string, RNEMDMethod>::iterator mi;
2144 for(mi = stringToMethod_.begin(); mi != stringToMethod_.end(); ++mi) {
2145 if ( (*mi).second == rnemdMethod_)
2146 rnemdFile_ << "# exchangeMethod = \"" << (*mi).first << "\";\n";
2147 }
2148 map<string, RNEMDFluxType>::iterator fi;
2149 for(fi = stringToFluxType_.begin(); fi != stringToFluxType_.end(); ++fi) {
2150 if ( (*fi).second == rnemdFluxType_)
2151 rnemdFile_ << "# fluxType = \"" << (*fi).first << "\";\n";
2152 }
2153
2154 rnemdFile_ << "# exchangeTime = " << exchangeTime_ << ";\n";
2155
2156 rnemdFile_ << "# objectSelection = \""
2157 << rnemdObjectSelection_ << "\";\n";
2158 rnemdFile_ << "# selectionA = \"" << selectionA_ << "\";\n";
2159 rnemdFile_ << "# selectionB = \"" << selectionB_ << "\";\n";
2160 rnemdFile_ << "# }\n";
2161 rnemdFile_ << "#######################################################\n";
2162 rnemdFile_ << "# RNEMD report:\n";
2163 rnemdFile_ << "# running time = " << time << " fs\n";
2164 rnemdFile_ << "# Target flux:\n";
2165 rnemdFile_ << "# kinetic = "
2166 << kineticFlux_ / PhysicalConstants::energyConvert
2167 << " (kcal/mol/A^2/fs)\n";
2168 rnemdFile_ << "# momentum = " << momentumFluxVector_
2169 << " (amu/A/fs^2)\n";
2170 rnemdFile_ << "# angular momentum = " << angularMomentumFluxVector_
2171 << " (amu/A^2/fs^2)\n";
2172 rnemdFile_ << "# Target one-time exchanges:\n";
2173 rnemdFile_ << "# kinetic = "
2174 << kineticTarget_ / PhysicalConstants::energyConvert
2175 << " (kcal/mol)\n";
2176 rnemdFile_ << "# momentum = " << momentumTarget_
2177 << " (amu*A/fs)\n";
2178 rnemdFile_ << "# angular momentum = " << angularMomentumTarget_
2179 << " (amu*A^2/fs)\n";
2180 rnemdFile_ << "# Actual exchange totals:\n";
2181 rnemdFile_ << "# kinetic = "
2182 << kineticExchange_ / PhysicalConstants::energyConvert
2183 << " (kcal/mol)\n";
2184 rnemdFile_ << "# momentum = " << momentumExchange_
2185 << " (amu*A/fs)\n";
2186 rnemdFile_ << "# angular momentum = " << angularMomentumExchange_
2187 << " (amu*A^2/fs)\n";
2188 rnemdFile_ << "# Actual flux:\n";
2189 rnemdFile_ << "# kinetic = " << Jz
2190 << " (kcal/mol/A^2/fs)\n";
2191 rnemdFile_ << "# momentum = " << JzP
2192 << " (amu/A/fs^2)\n";
2193 rnemdFile_ << "# angular momentum = " << JzL
2194 << " (amu/A^2/fs^2)\n";
2195 rnemdFile_ << "# Exchange statistics:\n";
2196 rnemdFile_ << "# attempted = " << trialCount_ << "\n";
2197 rnemdFile_ << "# failed = " << failTrialCount_ << "\n";
2198 if (rnemdMethod_ == rnemdNIVS) {
2199 rnemdFile_ << "# NIVS root-check errors = "
2200 << failRootCount_ << "\n";
2201 }
2202 rnemdFile_ << "#######################################################\n";
2203
2204
2205
2206 //write title
2207 rnemdFile_ << "#";
2208 for (unsigned int i = 0; i < outputMask_.size(); ++i) {
2209 if (outputMask_[i]) {
2210 rnemdFile_ << "\t" << data_[i].title <<
2211 "(" << data_[i].units << ")";
2212 // add some extra tabs for column alignment
2213 if (data_[i].dataType == "Vector3d") rnemdFile_ << "\t\t";
2214 }
2215 }
2216 rnemdFile_ << std::endl;
2217
2218 rnemdFile_.precision(8);
2219
2220 for (int j = 0; j < nBins_; j++) {
2221
2222 for (unsigned int i = 0; i < outputMask_.size(); ++i) {
2223 if (outputMask_[i]) {
2224 if (data_[i].dataType == "RealType")
2225 writeReal(i,j);
2226 else if (data_[i].dataType == "Vector3d")
2227 writeVector(i,j);
2228 else {
2229 sprintf( painCave.errMsg,
2230 "RNEMD found an unknown data type for: %s ",
2231 data_[i].title.c_str());
2232 painCave.isFatal = 1;
2233 simError();
2234 }
2235 }
2236 }
2237 rnemdFile_ << std::endl;
2238
2239 }
2240
2241 rnemdFile_ << "#######################################################\n";
2242 rnemdFile_ << "# Standard Deviations in those quantities follow:\n";
2243 rnemdFile_ << "#######################################################\n";
2244
2245
2246 for (int j = 0; j < nBins_; j++) {
2247 rnemdFile_ << "#";
2248 for (unsigned int i = 0; i < outputMask_.size(); ++i) {
2249 if (outputMask_[i]) {
2250 if (data_[i].dataType == "RealType")
2251 writeRealStdDev(i,j);
2252 else if (data_[i].dataType == "Vector3d")
2253 writeVectorStdDev(i,j);
2254 else {
2255 sprintf( painCave.errMsg,
2256 "RNEMD found an unknown data type for: %s ",
2257 data_[i].title.c_str());
2258 painCave.isFatal = 1;
2259 simError();
2260 }
2261 }
2262 }
2263 rnemdFile_ << std::endl;
2264
2265 }
2266
2267 rnemdFile_.flush();
2268 rnemdFile_.close();
2269
2270 #ifdef IS_MPI
2271 }
2272 #endif
2273
2274 }
2275
2276 void RNEMD::writeReal(int index, unsigned int bin) {
2277 if (!doRNEMD_) return;
2278 assert(index >=0 && index < ENDINDEX);
2279 assert(int(bin) < nBins_);
2280 RealType s;
2281 int count;
2282
2283 count = data_[index].accumulator[bin]->count();
2284 if (count == 0) return;
2285
2286 dynamic_cast<Accumulator *>(data_[index].accumulator[bin])->getAverage(s);
2287
2288 if (! isinf(s) && ! isnan(s)) {
2289 rnemdFile_ << "\t" << s;
2290 } else{
2291 sprintf( painCave.errMsg,
2292 "RNEMD detected a numerical error writing: %s for bin %u",
2293 data_[index].title.c_str(), bin);
2294 painCave.isFatal = 1;
2295 simError();
2296 }
2297 }
2298
2299 void RNEMD::writeVector(int index, unsigned int bin) {
2300 if (!doRNEMD_) return;
2301 assert(index >=0 && index < ENDINDEX);
2302 assert(int(bin) < nBins_);
2303 Vector3d s;
2304 int count;
2305
2306 count = data_[index].accumulator[bin]->count();
2307
2308 if (count == 0) return;
2309
2310 dynamic_cast<VectorAccumulator*>(data_[index].accumulator[bin])->getAverage(s);
2311 if (isinf(s[0]) || isnan(s[0]) ||
2312 isinf(s[1]) || isnan(s[1]) ||
2313 isinf(s[2]) || isnan(s[2]) ) {
2314 sprintf( painCave.errMsg,
2315 "RNEMD detected a numerical error writing: %s for bin %u",
2316 data_[index].title.c_str(), bin);
2317 painCave.isFatal = 1;
2318 simError();
2319 } else {
2320 rnemdFile_ << "\t" << s[0] << "\t" << s[1] << "\t" << s[2];
2321 }
2322 }
2323
2324 void RNEMD::writeRealStdDev(int index, unsigned int bin) {
2325 if (!doRNEMD_) return;
2326 assert(index >=0 && index < ENDINDEX);
2327 assert(int(bin) < nBins_);
2328 RealType s;
2329 int count;
2330
2331 count = data_[index].accumulator[bin]->count();
2332 if (count == 0) return;
2333
2334 dynamic_cast<Accumulator *>(data_[index].accumulator[bin])->getStdDev(s);
2335
2336 if (! isinf(s) && ! isnan(s)) {
2337 rnemdFile_ << "\t" << s;
2338 } else{
2339 sprintf( painCave.errMsg,
2340 "RNEMD detected a numerical error writing: %s std. dev. for bin %u",
2341 data_[index].title.c_str(), bin);
2342 painCave.isFatal = 1;
2343 simError();
2344 }
2345 }
2346
2347 void RNEMD::writeVectorStdDev(int index, unsigned int bin) {
2348 if (!doRNEMD_) return;
2349 assert(index >=0 && index < ENDINDEX);
2350 assert(int(bin) < nBins_);
2351 Vector3d s;
2352 int count;
2353
2354 count = data_[index].accumulator[bin]->count();
2355 if (count == 0) return;
2356
2357 dynamic_cast<VectorAccumulator*>(data_[index].accumulator[bin])->getStdDev(s);
2358 if (isinf(s[0]) || isnan(s[0]) ||
2359 isinf(s[1]) || isnan(s[1]) ||
2360 isinf(s[2]) || isnan(s[2]) ) {
2361 sprintf( painCave.errMsg,
2362 "RNEMD detected a numerical error writing: %s std. dev. for bin %u",
2363 data_[index].title.c_str(), bin);
2364 painCave.isFatal = 1;
2365 simError();
2366 } else {
2367 rnemdFile_ << "\t" << s[0] << "\t" << s[1] << "\t" << s[2];
2368 }
2369 }
2370 }
2371

Properties

Name Value
svn:keywords Author Id Revision Date