ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/rnemd/RNEMD.cpp
Revision: 1855
Committed: Tue Apr 2 18:31:51 2013 UTC (12 years, 1 month ago) by gezelter
File size: 73219 byte(s)
Log Message: 
Fixed a bunch of bugs in CubicSpline debug sections, ForceMatrix Decomp 
computing costhetas for non-existent rotation matrices, Hidden accumulator counts, and SMIPD non-coupling

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

Properties

Name Value
svn:keywords Author Id Revision Date