ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/rnemd/RNEMD.cpp
Revision: 1854
Committed: Thu Mar 28 20:54:06 2013 UTC (12 years, 1 month ago) by gezelter
File size: 72707 byte(s)
Log Message: 
First pass at angular momentum aware version of VSS-RNEMD.

File Contents

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

Properties

Name Value
svn:keywords Author Id Revision Date