ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/rnemd/RNEMD.cpp
Revision: 2026
Committed: Wed Oct 22 12:23:59 2014 UTC (10 years, 6 months ago) by gezelter
File size: 78086 byte(s)
Log Message: 
Starting to add support for UniformGradient. 
Changed Vector3d input type to a more general std::vector<RealType> input.  This change alters RNEMD and UniformField inputs.

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

Properties

Name Value
svn:keywords Author Id Revision Date