ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/rnemd/RNEMD.cpp
(Generate patch)

Comparing branches/development/src/rnemd/RNEMD.cpp (file contents):
Revision 1803 by gezelter, Wed Oct 3 14:20:07 2012 UTC vs.
Revision 1856 by gezelter, Tue Apr 2 21:30:34 2013 UTC

# Line 35 | Line 35
35   *                                                                      
36   * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).            
37   * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
38 < * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
38 > * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).          
39   * [4]  Vardeman & Gezelter, in progress (2009).                        
40   */
41  
42   #include <cmath>
43 + #include <sstream>
44 + #include <string>
45 +
46   #include "rnemd/RNEMD.hpp"
47   #include "math/Vector3.hpp"
48   #include "math/Vector.hpp"
# Line 49 | Line 52
52   #include "primitives/StuntDouble.hpp"
53   #include "utils/PhysicalConstants.hpp"
54   #include "utils/Tuple.hpp"
55 + #include "brains/Thermo.hpp"
56 + #include "math/ConvexHull.hpp"
57   #ifdef IS_MPI
58   #include <mpi.h>
59   #endif
60  
61 + #ifdef _MSC_VER
62 + #define isnan(x) _isnan((x))
63 + #define isinf(x) (!_finite(x) && !_isnan(x))
64 + #endif
65 +
66   #define HONKING_LARGE_VALUE 1.0e10
67  
68   using namespace std;
69   namespace OpenMD {
70    
71    RNEMD::RNEMD(SimInfo* info) : info_(info), evaluator_(info), seleMan_(info),
72 +                                evaluatorA_(info), seleManA_(info),
73 +                                commonA_(info), evaluatorB_(info),
74 +                                seleManB_(info), commonB_(info),
75                                  usePeriodicBoundaryConditions_(info->getSimParams()->getUsePeriodicBoundaryConditions()) {
76  
77      trialCount_ = 0;
78      failTrialCount_ = 0;
79      failRootCount_ = 0;
80  
81 <    int seedValue;
69 <    Globals * simParams = info->getSimParams();
81 >    Globals* simParams = info->getSimParams();
82      RNEMDParameters* rnemdParams = simParams->getRNEMDParameters();
83  
84      doRNEMD_ = rnemdParams->getUseRNEMD();
# Line 81 | Line 93 | namespace OpenMD {
93      stringToFluxType_["Py"]  = rnemdPy;
94      stringToFluxType_["Pz"]  = rnemdPz;
95      stringToFluxType_["Pvector"]  = rnemdPvector;
96 +    stringToFluxType_["Lx"]  = rnemdLx;
97 +    stringToFluxType_["Ly"]  = rnemdLy;
98 +    stringToFluxType_["Lz"]  = rnemdLz;
99 +    stringToFluxType_["Lvector"]  = rnemdLvector;
100      stringToFluxType_["KE+Px"]  = rnemdKePx;
101      stringToFluxType_["KE+Py"]  = rnemdKePy;
102      stringToFluxType_["KE+Pvector"]  = rnemdKePvector;
103 +    stringToFluxType_["KE+Lx"]  = rnemdKeLx;
104 +    stringToFluxType_["KE+Ly"]  = rnemdKeLy;
105 +    stringToFluxType_["KE+Lz"]  = rnemdKeLz;
106 +    stringToFluxType_["KE+Lvector"]  = rnemdKeLvector;
107  
108      runTime_ = simParams->getRunTime();
109      statusTime_ = simParams->getStatusTime();
110  
91    rnemdObjectSelection_ = rnemdParams->getObjectSelection();
92    evaluator_.loadScriptString(rnemdObjectSelection_);
93    seleMan_.setSelectionSet(evaluator_.evaluate());
94
111      const string methStr = rnemdParams->getMethod();
112      bool hasFluxType = rnemdParams->haveFluxType();
113  
114 +    rnemdObjectSelection_ = rnemdParams->getObjectSelection();
115 +
116      string fluxStr;
117      if (hasFluxType) {
118        fluxStr = rnemdParams->getFluxType();
# Line 102 | Line 120 | namespace OpenMD {
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 <              "\tKE, Px, Py, Pz, Pvector, KE+Px, KE+Py, KE+Pvector\n"
123 >              "\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                "\tmust be set to use RNEMD\n");
126        painCave.isFatal = 1;
127        painCave.severity = OPENMD_ERROR;
# Line 112 | Line 131 | namespace OpenMD {
131      bool hasKineticFlux = rnemdParams->haveKineticFlux();
132      bool hasMomentumFlux = rnemdParams->haveMomentumFlux();
133      bool hasMomentumFluxVector = rnemdParams->haveMomentumFluxVector();
134 +    bool hasAngularMomentumFlux = rnemdParams->haveAngularMomentumFlux();
135 +    bool hasAngularMomentumFluxVector = rnemdParams->haveAngularMomentumFluxVector();
136 +    hasSelectionA_ = rnemdParams->haveSelectionA();
137 +    hasSelectionB_ = rnemdParams->haveSelectionB();
138      bool hasSlabWidth = rnemdParams->haveSlabWidth();
139      bool hasSlabACenter = rnemdParams->haveSlabACenter();
140      bool hasSlabBCenter = rnemdParams->haveSlabBCenter();
141 +    bool hasSphereARadius = rnemdParams->haveSphereARadius();
142 +    hasSphereBRadius_ = rnemdParams->haveSphereBRadius();
143 +    bool hasCoordinateOrigin = rnemdParams->haveCoordinateOrigin();
144      bool hasOutputFileName = rnemdParams->haveOutputFileName();
145      bool hasOutputFields = rnemdParams->haveOutputFields();
146      
# Line 199 | Line 225 | namespace OpenMD {
225        case rnemdPz:
226          hasCorrectFlux = hasMomentumFlux;
227          break;
228 +      case rnemdLx:
229 +      case rnemdLy:
230 +      case rnemdLz:
231 +        hasCorrectFlux = hasAngularMomentumFlux;
232 +        break;
233        case rnemdPvector:
234          hasCorrectFlux = hasMomentumFluxVector;
235          break;
236 +      case rnemdLvector:
237 +        hasCorrectFlux = hasAngularMomentumFluxVector;
238 +        break;
239        case rnemdKePx:
240        case rnemdKePy:
241          hasCorrectFlux = hasMomentumFlux && hasKineticFlux;
242          break;
243 +      case rnemdKeLx:
244 +      case rnemdKeLy:
245 +      case rnemdKeLz:
246 +        hasCorrectFlux = hasAngularMomentumFlux && hasKineticFlux;
247 +        break;
248        case rnemdKePvector:
249          hasCorrectFlux = hasMomentumFluxVector && hasKineticFlux;
250          break;
251 +      case rnemdKeLvector:
252 +        hasCorrectFlux = hasAngularMomentumFluxVector && hasKineticFlux;
253 +        break;
254        default:
255          methodFluxMismatch = true;
256          break;
# Line 231 | Line 273 | namespace OpenMD {
273        sprintf(painCave.errMsg,
274                "RNEMD: The current method, %s, and flux type, %s,\n"
275                "\tdid not have the correct flux value specified. Options\n"
276 <              "\tinclude: kineticFlux, momentumFlux, and momentumFluxVector\n",
276 >              "\tinclude: kineticFlux, momentumFlux, angularMomentumFlux,\n"
277 >              "\tmomentumFluxVector, and angularMomentumFluxVector.\n",
278                methStr.c_str(), fluxStr.c_str());
279        painCave.isFatal = 1;
280        painCave.severity = OPENMD_ERROR;
# Line 271 | Line 314 | namespace OpenMD {
314          default:
315            break;
316          }
317 <      }    
318 <    }
317 >      }
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  
349 <    // do some sanity checking
350 <
351 <    int selectionCount = seleMan_.getSelectionCount();
349 >      if (hasCoordinateOrigin) {
350 >        coordinateOrigin_ = rnemdParams->getCoordinateOrigin();
351 >      } else {
352 >        coordinateOrigin_ = V3Zero;
353 >      }
354  
355 <    int nIntegrable = info->getNGlobalIntegrableObjects();
355 >      // do some sanity checking
356  
357 <    if (selectionCount > nIntegrable) {
284 <      sprintf(painCave.errMsg,
285 <              "RNEMD: The current objectSelection,\n"
286 <              "\t\t%s\n"
287 <              "\thas resulted in %d selected objects.  However,\n"
288 <              "\tthe total number of integrable objects in the system\n"
289 <              "\tis only %d.  This is almost certainly not what you want\n"
290 <              "\tto do.  A likely cause of this is forgetting the _RB_0\n"
291 <              "\tselector in the selection script!\n",
292 <              rnemdObjectSelection_.c_str(),
293 <              selectionCount, nIntegrable);
294 <      painCave.isFatal = 0;
295 <      painCave.severity = OPENMD_WARNING;
296 <      simError();
297 <    }
357 >      int selectionCount = seleMan_.getSelectionCount();
358  
359 <    areaAccumulator_ = new Accumulator();
359 >      int nIntegrable = info->getNGlobalIntegrableObjects();
360  
361 <    nBins_ = rnemdParams->getOutputBins();
361 >      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  
377 <    data_.resize(RNEMD::ENDINDEX);
304 <    OutputData z;
305 <    z.units =  "Angstroms";
306 <    z.title =  "Z";
307 <    z.dataType = "RealType";
308 <    z.accumulator.reserve(nBins_);
309 <    for (unsigned int i = 0; i < nBins_; i++)
310 <      z.accumulator.push_back( new Accumulator() );
311 <    data_[Z] = z;
312 <    outputMap_["Z"] =  Z;
377 >      areaAccumulator_ = new Accumulator();
378  
379 <    OutputData temperature;
380 <    temperature.units =  "K";
316 <    temperature.title =  "Temperature";
317 <    temperature.dataType = "RealType";
318 <    temperature.accumulator.reserve(nBins_);
319 <    for (unsigned int i = 0; i < nBins_; i++)
320 <      temperature.accumulator.push_back( new Accumulator() );
321 <    data_[TEMPERATURE] = temperature;
322 <    outputMap_["TEMPERATURE"] =  TEMPERATURE;
379 >      nBins_ = rnemdParams->getOutputBins();
380 >      binWidth_ = rnemdParams->getOutputBinWidth();
381  
382 <    OutputData velocity;
383 <    velocity.units = "angstroms/fs";
384 <    velocity.title =  "Velocity";  
385 <    velocity.dataType = "Vector3d";
386 <    velocity.accumulator.reserve(nBins_);
387 <    for (unsigned int i = 0; i < nBins_; i++)
388 <      velocity.accumulator.push_back( new VectorAccumulator() );
389 <    data_[VELOCITY] = velocity;
390 <    outputMap_["VELOCITY"] = VELOCITY;
382 >      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 >
393 >      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  
403 <    OutputData density;
404 <    density.units =  "g cm^-3";
405 <    density.title =  "Density";
406 <    density.dataType = "RealType";
407 <    density.accumulator.reserve(nBins_);
408 <    for (unsigned int i = 0; i < nBins_; i++)
409 <      density.accumulator.push_back( new Accumulator() );
410 <    data_[DENSITY] = density;
411 <    outputMap_["DENSITY"] =  DENSITY;
403 >      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  
413 <    if (hasOutputFields) {
414 <      parseOutputFileFormat(rnemdParams->getOutputFields());
415 <    } else {
416 <      outputMask_.set(Z);
417 <      switch (rnemdFluxType_) {
418 <      case rnemdKE:
419 <      case rnemdRotKE:
420 <      case rnemdFullKE:
421 <        outputMask_.set(TEMPERATURE);
353 <        break;
354 <      case rnemdPx:
355 <      case rnemdPy:
356 <        outputMask_.set(VELOCITY);
357 <        break;
358 <      case rnemdPz:        
359 <      case rnemdPvector:
360 <        outputMask_.set(VELOCITY);
361 <        outputMask_.set(DENSITY);
362 <        break;
363 <      case rnemdKePx:
364 <      case rnemdKePy:
365 <        outputMask_.set(TEMPERATURE);
366 <        outputMask_.set(VELOCITY);
367 <        break;
368 <      case rnemdKePvector:
369 <        outputMask_.set(TEMPERATURE);
370 <        outputMask_.set(VELOCITY);
371 <        outputMask_.set(DENSITY);        
372 <        break;
373 <      default:
374 <        break;
375 <      }
376 <    }
377 <      
378 <    if (hasOutputFileName) {
379 <      rnemdFileName_ = rnemdParams->getOutputFileName();
380 <    } else {
381 <      rnemdFileName_ = getPrefix(info->getFinalConfigFileName()) + ".rnemd";
382 <    }          
413 >      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 <    exchangeTime_ = rnemdParams->getExchangeTime();
423 >      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 <    Snapshot* currentSnap_ = info->getSnapshotManager()->getCurrentSnapshot();
434 <    Mat3x3d hmat = currentSnap_->getHmat();
435 <  
436 <    // Target exchange quantities (in each exchange) =  2 Lx Ly dt flux
437 <    // Lx, Ly = box dimensions in x & y
438 <    // dt = exchange time interval
439 <    // flux = target flux
433 >      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 <    RealType area = currentSnap_->getXYarea();
444 <    kineticTarget_ = 2.0 * kineticFlux_ * exchangeTime_ * area;
445 <    momentumTarget_ = 2.0 * momentumFluxVector_ * exchangeTime_ * area;
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 >      }
492 >      
493 >      if (hasOutputFileName) {
494 >        rnemdFileName_ = rnemdParams->getOutputFileName();
495 >      } else {
496 >        rnemdFileName_ = getPrefix(info->getFinalConfigFileName()) + ".rnemd";
497 >      }          
498  
499 <    // total exchange sums are zeroed out at the beginning:
499 >      exchangeTime_ = rnemdParams->getExchangeTime();
500  
501 <    kineticExchange_ = 0.0;
502 <    momentumExchange_ = V3Zero;
501 >      Snapshot* currentSnap_ = info->getSnapshotManager()->getCurrentSnapshot();
502 >      // total exchange sums are zeroed out at the beginning:
503  
504 <    if (hasSlabWidth)
505 <      slabWidth_ = rnemdParams->getSlabWidth();
506 <    else
507 <      slabWidth_ = hmat(2,2) / 10.0;
508 <  
509 <    if (hasSlabACenter)
409 <      slabACenter_ = rnemdParams->getSlabACenter();
410 <    else
411 <      slabACenter_ = 0.0;
504 >      kineticExchange_ = 0.0;
505 >      momentumExchange_ = V3Zero;
506 >      angularMomentumExchange_ = V3Zero;
507 >
508 >      std::ostringstream selectionAstream;
509 >      std::ostringstream selectionBstream;
510      
511 <    if (hasSlabBCenter)
512 <      slabBCenter_ = rnemdParams->getSlabBCenter();
513 <    else
514 <      slabBCenter_ = hmat(2,2) / 2.0;
511 >      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      
547 +      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    }
593 <  
593 >  
594 >    
595    RNEMD::~RNEMD() {
596      if (!doRNEMD_) return;
597   #ifdef IS_MPI
# Line 432 | Line 607 | namespace OpenMD {
607   #endif
608    }
609    
610 <  bool RNEMD::inSlabA(Vector3d pos) {
436 <    return (abs(pos.z() - slabACenter_) < 0.5*slabWidth_);
437 <  }
438 <  bool RNEMD::inSlabB(Vector3d pos) {
439 <    return (abs(pos.z() - slabBCenter_) < 0.5*slabWidth_);
440 <  }
441 <
442 <  void RNEMD::doSwap() {
610 >  void RNEMD::doSwap(SelectionManager& smanA, SelectionManager& smanB) {
611      if (!doRNEMD_) return;
612 +    int selei;
613 +    int selej;
614 +
615      Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
616      Mat3x3d hmat = currentSnap_->getHmat();
617  
447    seleMan_.setSelectionSet(evaluator_.evaluate());
448
449    int selei;
618      StuntDouble* sd;
451    int idx;
619  
620      RealType min_val;
621      bool min_found = false;  
# Line 458 | Line 625 | namespace OpenMD {
625      bool max_found = false;
626      StuntDouble* max_sd;
627  
628 <    for (sd = seleMan_.beginSelected(selei); sd != NULL;
629 <         sd = seleMan_.nextSelected(selei)) {
628 >    for (sd = seleManA_.beginSelected(selei); sd != NULL;
629 >         sd = seleManA_.nextSelected(selei)) {
630  
464      idx = sd->getLocalIndex();
465
631        Vector3d pos = sd->getPos();
632 <
632 >      
633        // wrap the stuntdouble's position back into the box:
634 <
634 >      
635        if (usePeriodicBoundaryConditions_)
636          currentSnap_->wrapVector(pos);
637 <      bool inA = inSlabA(pos);
638 <      bool inB = inSlabB(pos);
639 <
640 <      if (inA || inB) {
637 >      
638 >      RealType mass = sd->getMass();
639 >      Vector3d vel = sd->getVel();
640 >      RealType value;
641 >      
642 >      switch(rnemdFluxType_) {
643 >      case rnemdKE :
644          
645 <        RealType mass = sd->getMass();
646 <        Vector3d vel = sd->getVel();
647 <        RealType value;
648 <        
649 <        switch(rnemdFluxType_) {
650 <        case rnemdKE :
651 <          
652 <          value = mass * vel.lengthSquare();
653 <          
654 <          if (sd->isDirectional()) {
655 <            Vector3d angMom = sd->getJ();
656 <            Mat3x3d I = sd->getI();
657 <            
658 <            if (sd->isLinear()) {
659 <              int i = sd->linearAxis();
660 <              int j = (i + 1) % 3;
661 <              int k = (i + 2) % 3;
662 <              value += angMom[j] * angMom[j] / I(j, j) +
663 <                angMom[k] * angMom[k] / I(k, k);
664 <            } else {                        
665 <              value += angMom[0]*angMom[0]/I(0, 0)
666 <                + angMom[1]*angMom[1]/I(1, 1)
667 <                + angMom[2]*angMom[2]/I(2, 2);
668 <            }
669 <          } //angular momenta exchange enabled
670 <          value *= 0.5;
671 <          break;
672 <        case rnemdPx :
673 <          value = mass * vel[0];
674 <          break;
675 <        case rnemdPy :
676 <          value = mass * vel[1];
677 <          break;
678 <        case rnemdPz :
679 <          value = mass * vel[2];
680 <          break;
681 <        default :
682 <          break;
645 >        value = mass * vel.lengthSquare();
646 >        
647 >        if (sd->isDirectional()) {
648 >          Vector3d angMom = sd->getJ();
649 >          Mat3x3d I = sd->getI();
650 >          
651 >          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          }
686 +      }  
687 +    }
688          
689 <        if (inA == 0) {
690 <          if (!min_found) {
691 <            min_val = value;
692 <            min_sd = sd;
693 <            min_found = true;
694 <          } else {
695 <            if (min_val > value) {
696 <              min_val = value;
697 <              min_sd = sd;
698 <            }
699 <          }
700 <        } else {
701 <          if (!max_found) {
702 <            max_val = value;
703 <            max_sd = sd;
704 <            max_found = true;
705 <          } else {
706 <            if (max_val < value) {
707 <              max_val = value;
708 <              max_sd = sd;
709 <            }
710 <          }      
711 <        }
689 >    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        }
738 +      
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      }
750      
751 < #ifdef IS_MPI
752 <    int nProc, worldRank;
751 > #ifdef IS_MPI    
752 >    int worldRank = MPI::COMM_WORLD.Get_rank();
753      
546    nProc = MPI::COMM_WORLD.Get_size();
547    worldRank = MPI::COMM_WORLD.Get_rank();
548
754      bool my_min_found = min_found;
755      bool my_max_found = max_found;
756  
# Line 769 | Line 974 | namespace OpenMD {
974      }    
975    }
976    
977 <  void RNEMD::doNIVS() {
977 >  void RNEMD::doNIVS(SelectionManager& smanA, SelectionManager& smanB) {
978      if (!doRNEMD_) return;
979 +    int selei;
980 +    int selej;
981 +
982      Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
983 +    RealType time = currentSnap_->getTime();    
984      Mat3x3d hmat = currentSnap_->getHmat();
985  
777    seleMan_.setSelectionSet(evaluator_.evaluate());
778
779    int selei;
986      StuntDouble* sd;
781    int idx;
987  
988      vector<StuntDouble*> hotBin, coldBin;
989  
# Line 797 | Line 1002 | namespace OpenMD {
1002      RealType Kcz = 0.0;
1003      RealType Kcw = 0.0;
1004  
1005 <    for (sd = seleMan_.beginSelected(selei); sd != NULL;
1006 <         sd = seleMan_.nextSelected(selei)) {
1005 >    for (sd = smanA.beginSelected(selei); sd != NULL;
1006 >         sd = smanA.nextSelected(selei)) {
1007  
803      idx = sd->getLocalIndex();
804
1008        Vector3d pos = sd->getPos();
1009 <
1009 >      
1010        // wrap the stuntdouble's position back into the box:
1011 <
1011 >      
1012        if (usePeriodicBoundaryConditions_)
1013          currentSnap_->wrapVector(pos);
1014 <
1015 <      // which bin is this stuntdouble in?
1016 <      bool inA = inSlabA(pos);
1017 <      bool inB = inSlabB(pos);
1014 >      
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  
1054 <      if (inA || inB) {
1055 <              
1056 <        RealType mass = sd->getMass();
1057 <        Vector3d vel = sd->getVel();
1058 <      
1059 <        if (inA) {
1060 <          hotBin.push_back(sd);
1061 <          Phx += mass * vel.x();
1062 <          Phy += mass * vel.y();
1063 <          Phz += mass * vel.z();
1064 <          Khx += mass * vel.x() * vel.x();
1065 <          Khy += mass * vel.y() * vel.y();
1066 <          Khz += mass * vel.z() * vel.z();
1067 <          if (sd->isDirectional()) {
1068 <            Vector3d angMom = sd->getJ();
1069 <            Mat3x3d I = sd->getI();
1070 <            if (sd->isLinear()) {
1071 <              int i = sd->linearAxis();
1072 <              int j = (i + 1) % 3;
1073 <              int k = (i + 2) % 3;
1074 <              Khw += angMom[j] * angMom[j] / I(j, j) +
837 <                angMom[k] * angMom[k] / I(k, k);
838 <            } else {
839 <              Khw += angMom[0]*angMom[0]/I(0, 0)
840 <                + angMom[1]*angMom[1]/I(1, 1)
841 <                + angMom[2]*angMom[2]/I(2, 2);
842 <            }
843 <          }
844 <        } else {
845 <          coldBin.push_back(sd);
846 <          Pcx += mass * vel.x();
847 <          Pcy += mass * vel.y();
848 <          Pcz += mass * vel.z();
849 <          Kcx += mass * vel.x() * vel.x();
850 <          Kcy += mass * vel.y() * vel.y();
851 <          Kcz += mass * vel.z() * vel.z();
852 <          if (sd->isDirectional()) {
853 <            Vector3d angMom = sd->getJ();
854 <            Mat3x3d I = sd->getI();
855 <            if (sd->isLinear()) {
856 <              int i = sd->linearAxis();
857 <              int j = (i + 1) % 3;
858 <              int k = (i + 2) % 3;
859 <              Kcw += angMom[j] * angMom[j] / I(j, j) +
860 <                angMom[k] * angMom[k] / I(k, k);
861 <            } else {
862 <              Kcw += angMom[0]*angMom[0]/I(0, 0)
863 <                + angMom[1]*angMom[1]/I(1, 1)
864 <                + angMom[2]*angMom[2]/I(2, 2);
865 <            }
866 <          }
867 <        }
1054 >      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        }
1076      }
1077      
# Line 1214 | Line 1421 | namespace OpenMD {
1421        failTrialCount_++;
1422      }
1423    }
1424 <
1425 <  void RNEMD::doVSS() {
1424 >  
1425 >  void RNEMD::doVSS(SelectionManager& smanA, SelectionManager& smanB) {
1426      if (!doRNEMD_) return;
1427 +    int selei;
1428 +    int selej;
1429 +
1430      Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
1431      RealType time = currentSnap_->getTime();    
1432      Mat3x3d hmat = currentSnap_->getHmat();
1433  
1224    seleMan_.setSelectionSet(evaluator_.evaluate());
1225
1226    int selei;
1434      StuntDouble* sd;
1228    int idx;
1435  
1436      vector<StuntDouble*> hotBin, coldBin;
1437  
1438      Vector3d Ph(V3Zero);
1439 +    Vector3d Lh(V3Zero);
1440      RealType Mh = 0.0;
1441 +    Mat3x3d Ih(0.0);
1442      RealType Kh = 0.0;
1443      Vector3d Pc(V3Zero);
1444 +    Vector3d Lc(V3Zero);
1445      RealType Mc = 0.0;
1446 +    Mat3x3d Ic(0.0);
1447      RealType Kc = 0.0;
1448      
1449 +    for (sd = smanA.beginSelected(selei); sd != NULL;
1450 +         sd = smanA.nextSelected(selei)) {
1451  
1240    for (sd = seleMan_.beginSelected(selei); sd != NULL;
1241         sd = seleMan_.nextSelected(selei)) {
1242
1243      idx = sd->getLocalIndex();
1244
1452        Vector3d pos = sd->getPos();
1453  
1454        // wrap the stuntdouble's position back into the box:
1455 +      
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  
1496 +      Vector3d pos = sd->getPos();
1497 +      
1498 +      // wrap the stuntdouble's position back into the box:
1499 +      
1500        if (usePeriodicBoundaryConditions_)
1501          currentSnap_->wrapVector(pos);
1502 +      
1503 +      RealType mass = sd->getMass();
1504 +      Vector3d vel = sd->getVel();
1505 +      Vector3d rPos = sd->getPos() - coordinateOrigin_;
1506 +      RealType r2;
1507  
1508 <      // which bin is this stuntdouble in?
1509 <      bool inA = inSlabA(pos);
1510 <      bool inB = inSlabB(pos);
1508 >      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        
1519 <      if (inA || inB) {
1520 <        
1521 <        RealType mass = sd->getMass();
1522 <        Vector3d vel = sd->getVel();
1523 <      
1524 <        if (inA) {
1525 <          hotBin.push_back(sd);
1526 <          Ph += mass * vel;
1527 <          Mh += mass;
1528 <          Kh += mass * vel.lengthSquare();
1529 <          if (rnemdFluxType_ == rnemdFullKE) {
1530 <            if (sd->isDirectional()) {
1531 <              Vector3d angMom = sd->getJ();
1532 <              Mat3x3d I = sd->getI();
1533 <              if (sd->isLinear()) {
1534 <                int i = sd->linearAxis();
1272 <                int j = (i + 1) % 3;
1273 <                int k = (i + 2) % 3;
1274 <                Kh += angMom[j] * angMom[j] / I(j, j) +
1275 <                  angMom[k] * angMom[k] / I(k, k);
1276 <              } else {
1277 <                Kh += angMom[0] * angMom[0] / I(0, 0) +
1278 <                  angMom[1] * angMom[1] / I(1, 1) +
1279 <                  angMom[2] * angMom[2] / I(2, 2);
1280 <              }
1281 <            }
1282 <          }
1283 <        } else { //midBin_
1284 <          coldBin.push_back(sd);
1285 <          Pc += mass * vel;
1286 <          Mc += mass;
1287 <          Kc += mass * vel.lengthSquare();
1288 <          if (rnemdFluxType_ == rnemdFullKE) {
1289 <            if (sd->isDirectional()) {
1290 <              Vector3d angMom = sd->getJ();
1291 <              Mat3x3d I = sd->getI();
1292 <              if (sd->isLinear()) {
1293 <                int i = sd->linearAxis();
1294 <                int j = (i + 1) % 3;
1295 <                int k = (i + 2) % 3;
1296 <                Kc += angMom[j] * angMom[j] / I(j, j) +
1297 <                  angMom[k] * angMom[k] / I(k, k);
1298 <              } else {
1299 <                Kc += angMom[0] * angMom[0] / I(0, 0) +
1300 <                  angMom[1] * angMom[1] / I(1, 1) +
1301 <                  angMom[2] * angMom[2] / I(2, 2);
1302 <              }
1303 <            }
1304 <          }
1305 <        }
1519 >      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        }
1536      }
1537      
# Line 1312 | Line 1541 | namespace OpenMD {
1541   #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 +    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      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 +    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   #endif
1555 <
1555 >    
1556      bool successfulExchange = false;
1557      if ((Mh > 0.0) && (Mc > 0.0)) {//both slabs are not empty
1558        Vector3d vc = Pc / Mc;
1559        Vector3d ac = -momentumTarget_ / Mc + vc;
1560        Vector3d acrec = -momentumTarget_ / Mc;
1561 <      RealType cNumerator = Kc - kineticTarget_ - 0.5 * Mc * ac.lengthSquare();
1561 >      
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        if (cNumerator > 0.0) {
1572 <        RealType cDenominator = Kc - 0.5 * Mc * vc.lengthSquare();
1572 >        
1573 >        RealType cDenominator = Kc - 0.5 * Mc * vc.lengthSquare()
1574 >          - 0.5*(dot(omegac, Ic * omegac));
1575 >        
1576          if (cDenominator > 0.0) {
1577            RealType c = sqrt(cNumerator / cDenominator);
1578            if ((c > 0.9) && (c < 1.1)) {//restrict scaling coefficients
1579 +            
1580              Vector3d vh = Ph / Mh;
1581              Vector3d ah = momentumTarget_ / Mh + vh;
1582              Vector3d ahrec = momentumTarget_ / Mh;
1583 <            RealType hNumerator = Kh + kineticTarget_
1584 <              - 0.5 * Mh * ah.lengthSquare();
1585 <            if (hNumerator > 0.0) {
1586 <              RealType hDenominator = Kh - 0.5 * Mh * vh.lengthSquare();
1583 >            
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                if (hDenominator > 0.0) {
1599                  RealType h = sqrt(hNumerator / hDenominator);
1600                  if ((h > 0.9) && (h < 1.1)) {
1601 <
1601 >                  
1602                    vector<StuntDouble*>::iterator sdi;
1603                    Vector3d vel;
1604 +                  Vector3d rPos;
1605 +                  
1606                    for (sdi = coldBin.begin(); sdi != coldBin.end(); sdi++) {
1607                      //vel = (*sdi)->getVel();
1608 <                    vel = ((*sdi)->getVel() - vc) * c + ac;
1608 >                    rPos = (*sdi)->getPos() - coordinateOrigin_;
1609 >                    vel = ((*sdi)->getVel() - vc - cross(omegac, rPos)) * c
1610 >                      + ac + cross(bc, rPos);
1611                      (*sdi)->setVel(vel);
1612                      if (rnemdFluxType_ == rnemdFullKE) {
1613                        if ((*sdi)->isDirectional()) {
# Line 1355 | Line 1618 | namespace OpenMD {
1618                    }
1619                    for (sdi = hotBin.begin(); sdi != hotBin.end(); sdi++) {
1620                      //vel = (*sdi)->getVel();
1621 <                    vel = ((*sdi)->getVel() - vh) * h + ah;
1621 >                    rPos = (*sdi)->getPos() - coordinateOrigin_;
1622 >                    vel = ((*sdi)->getVel() - vh - cross(omegah, rPos)) * h
1623 >                      + ah + cross(bh, rPos);    
1624                      (*sdi)->setVel(vel);
1625                      if (rnemdFluxType_ == rnemdFullKE) {
1626                        if ((*sdi)->isDirectional()) {
# Line 1367 | Line 1632 | namespace OpenMD {
1632                    successfulExchange = true;
1633                    kineticExchange_ += kineticTarget_;
1634                    momentumExchange_ += momentumTarget_;
1635 +                  angularMomentumExchange_ += angularMomentumTarget_;
1636                  }
1637                }
1638              }
# Line 1386 | Line 1652 | namespace OpenMD {
1652      }
1653    }
1654  
1655 +  RealType RNEMD::getDividingArea() {
1656 +
1657 +    if (hasDividingArea_) return dividingArea_;
1658 +
1659 +    RealType areaA, areaB;
1660 +    Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot();
1661 +
1662 +    if (hasSelectionA_) {
1663 +      int isd;
1664 +      StuntDouble* sd;
1665 +      vector<StuntDouble*> aSites;
1666 +      ConvexHull* surfaceMeshA = new ConvexHull();
1667 +      seleManA_.setSelectionSet(evaluatorA_.evaluate());
1668 +      for (sd = seleManA_.beginSelected(isd); sd != NULL;
1669 +           sd = seleManA_.nextSelected(isd)) {
1670 +        aSites.push_back(sd);
1671 +      }
1672 +      surfaceMeshA->computeHull(aSites);
1673 +      areaA = surfaceMeshA->getArea();
1674 +    } else {
1675 +      if (usePeriodicBoundaryConditions_) {
1676 +        // in periodic boundaries, the surface area is twice the x-y
1677 +        // area of the current box:
1678 +        areaA = 2.0 * snap->getXYarea();
1679 +      } else {
1680 +        // in non-periodic simulations, without explicitly setting
1681 +        // selections, the sphere radius sets the surface area of the
1682 +        // dividing surface:
1683 +        areaA = 4.0 * M_PI * pow(sphereARadius_, 2);
1684 +      }
1685 +    }
1686 +
1687 +    if (hasSelectionB_) {
1688 +      int isd;
1689 +      StuntDouble* sd;
1690 +      vector<StuntDouble*> bSites;
1691 +      ConvexHull* surfaceMeshB = new ConvexHull();
1692 +      seleManB_.setSelectionSet(evaluatorB_.evaluate());
1693 +      for (sd = seleManB_.beginSelected(isd); sd != NULL;
1694 +           sd = seleManB_.nextSelected(isd)) {
1695 +        bSites.push_back(sd);
1696 +      }
1697 +      surfaceMeshB->computeHull(bSites);
1698 +      areaB = surfaceMeshB->getArea();
1699 +    } else {
1700 +      if (usePeriodicBoundaryConditions_) {
1701 +        // in periodic boundaries, the surface area is twice the x-y
1702 +        // area of the current box:
1703 +        areaB = 2.0 * snap->getXYarea();
1704 +      } else {
1705 +        // in non-periodic simulations, without explicitly setting
1706 +        // selections, but if a sphereBradius has been set, just use that:
1707 +        areaB = 4.0 * M_PI * pow(sphereBRadius_, 2);
1708 +      }
1709 +    }
1710 +    
1711 +    dividingArea_ = min(areaA, areaB);
1712 +    hasDividingArea_ = true;
1713 +    return dividingArea_;
1714 +  }
1715 +  
1716    void RNEMD::doRNEMD() {
1717      if (!doRNEMD_) return;
1718      trialCount_++;
1719 +
1720 +    // object evaluator:
1721 +    evaluator_.loadScriptString(rnemdObjectSelection_);
1722 +    seleMan_.setSelectionSet(evaluator_.evaluate());
1723 +
1724 +    evaluatorA_.loadScriptString(selectionA_);
1725 +    evaluatorB_.loadScriptString(selectionB_);
1726 +
1727 +    seleManA_.setSelectionSet(evaluatorA_.evaluate());
1728 +    seleManB_.setSelectionSet(evaluatorB_.evaluate());
1729 +
1730 +    commonA_ = seleManA_ & seleMan_;
1731 +    commonB_ = seleManB_ & seleMan_;
1732 +
1733 +    // Target exchange quantities (in each exchange) = dividingArea * dt * flux
1734 +    // dt = exchange time interval
1735 +    // flux = target flux
1736 +    // dividingArea = smallest dividing surface between the two regions
1737 +
1738 +    hasDividingArea_ = false;
1739 +    RealType area = getDividingArea();
1740 +
1741 +    kineticTarget_ = kineticFlux_ * exchangeTime_ * area;
1742 +    momentumTarget_ = momentumFluxVector_ * exchangeTime_ * area;
1743 +    angularMomentumTarget_ = angularMomentumFluxVector_ * exchangeTime_ * area;
1744 +
1745      switch(rnemdMethod_) {
1746      case rnemdSwap:
1747 <      doSwap();
1747 >      doSwap(commonA_, commonB_);
1748        break;
1749      case rnemdNIVS:
1750 <      doNIVS();
1750 >      doNIVS(commonA_, commonB_);
1751        break;
1752      case rnemdVSS:
1753 <      doVSS();
1753 >      doVSS(commonA_, commonB_);
1754        break;
1755      case rnemdUnkownMethod:
1756      default :
# Line 1408 | Line 1761 | namespace OpenMD {
1761    void RNEMD::collectData() {
1762      if (!doRNEMD_) return;
1763      Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
1764 +    
1765 +    // collectData can be called more frequently than the doRNEMD, so use the
1766 +    // computed area from the last exchange time:
1767 +    RealType area = getDividingArea();
1768 +    areaAccumulator_->add(area);
1769      Mat3x3d hmat = currentSnap_->getHmat();
1412
1413    areaAccumulator_->add(currentSnap_->getXYarea());
1414
1770      seleMan_.setSelectionSet(evaluator_.evaluate());
1771  
1772      int selei(0);
1773      StuntDouble* sd;
1774 <    int idx;
1774 >    int binNo;
1775  
1776      vector<RealType> binMass(nBins_, 0.0);
1777      vector<RealType> binPx(nBins_, 0.0);
1778      vector<RealType> binPy(nBins_, 0.0);
1779      vector<RealType> binPz(nBins_, 0.0);
1780 +    vector<RealType> binOmegax(nBins_, 0.0);
1781 +    vector<RealType> binOmegay(nBins_, 0.0);
1782 +    vector<RealType> binOmegaz(nBins_, 0.0);
1783      vector<RealType> binKE(nBins_, 0.0);
1784      vector<int> binDOF(nBins_, 0);
1785      vector<int> binCount(nBins_, 0);
# Line 1429 | Line 1787 | namespace OpenMD {
1787      // alternative approach, track all molecules instead of only those
1788      // selected for scaling/swapping:
1789      /*
1790 <    SimInfo::MoleculeIterator miter;
1791 <    vector<StuntDouble*>::iterator iiter;
1792 <    Molecule* mol;
1793 <    StuntDouble* sd;
1794 <    for (mol = info_->beginMolecule(miter); mol != NULL;
1790 >      SimInfo::MoleculeIterator miter;
1791 >      vector<StuntDouble*>::iterator iiter;
1792 >      Molecule* mol;
1793 >      StuntDouble* sd;
1794 >      for (mol = info_->beginMolecule(miter); mol != NULL;
1795        mol = info_->nextMolecule(miter))
1796        sd is essentially sd
1797 <        for (sd = mol->beginIntegrableObject(iiter);
1798 <             sd != NULL;
1799 <             sd = mol->nextIntegrableObject(iiter))
1797 >      for (sd = mol->beginIntegrableObject(iiter);
1798 >      sd != NULL;
1799 >      sd = mol->nextIntegrableObject(iiter))
1800      */
1801  
1802      for (sd = seleMan_.beginSelected(selei); sd != NULL;
1803 <         sd = seleMan_.nextSelected(selei)) {
1803 >         sd = seleMan_.nextSelected(selei)) {    
1804      
1447      idx = sd->getLocalIndex();
1448      
1805        Vector3d pos = sd->getPos();
1806  
1807        // wrap the stuntdouble's position back into the box:
1808        
1809 <      if (usePeriodicBoundaryConditions_)
1809 >      if (usePeriodicBoundaryConditions_) {
1810          currentSnap_->wrapVector(pos);
1811 +        // which bin is this stuntdouble in?
1812 +        // wrapped positions are in the range [-0.5*hmat(2,2), +0.5*hmat(2,2)]
1813 +        // Shift molecules by half a box to have bins start at 0
1814 +        // The modulo operator is used to wrap the case when we are
1815 +        // beyond the end of the bins back to the beginning.
1816 +        binNo = int(nBins_ * (pos.z() / hmat(2,2) + 0.5)) % nBins_;
1817 +      } else {
1818 +        Vector3d rPos = pos - coordinateOrigin_;
1819 +        binNo = int(rPos.length() / binWidth_);
1820 +      }
1821  
1456
1457      // which bin is this stuntdouble in?
1458      // wrapped positions are in the range [-0.5*hmat(2,2), +0.5*hmat(2,2)]
1459      // Shift molecules by half a box to have bins start at 0
1460      // The modulo operator is used to wrap the case when we are
1461      // beyond the end of the bins back to the beginning.
1462      int binNo = int(nBins_ * (pos.z() / hmat(2,2) + 0.5)) % nBins_;
1463
1822        RealType mass = sd->getMass();
1823        Vector3d vel = sd->getVel();
1824 <
1825 <      binCount[binNo]++;
1826 <      binMass[binNo] += mass;
1827 <      binPx[binNo] += mass*vel.x();
1828 <      binPy[binNo] += mass*vel.y();
1829 <      binPz[binNo] += mass*vel.z();
1830 <      binKE[binNo] += 0.5 * (mass * vel.lengthSquare());
1831 <      binDOF[binNo] += 3;
1832 <
1833 <      if (sd->isDirectional()) {
1834 <        Vector3d angMom = sd->getJ();
1835 <        Mat3x3d I = sd->getI();
1836 <        if (sd->isLinear()) {
1837 <          int i = sd->linearAxis();
1838 <          int j = (i + 1) % 3;
1839 <          int k = (i + 2) % 3;
1840 <          binKE[binNo] += 0.5 * (angMom[j] * angMom[j] / I(j, j) +
1841 <                                 angMom[k] * angMom[k] / I(k, k));
1842 <          binDOF[binNo] += 2;
1843 <        } else {
1844 <          binKE[binNo] += 0.5 * (angMom[0] * angMom[0] / I(0, 0) +
1845 <                                 angMom[1] * angMom[1] / I(1, 1) +
1846 <                                 angMom[2] * angMom[2] / I(2, 2));
1847 <          binDOF[binNo] += 3;
1824 >      Vector3d rPos = sd->getPos() - coordinateOrigin_;
1825 >      Vector3d aVel = cross(rPos, vel);
1826 >      
1827 >      if (binNo >= 0 && binNo < nBins_)  {
1828 >        binCount[binNo]++;
1829 >        binMass[binNo] += mass;
1830 >        binPx[binNo] += mass*vel.x();
1831 >        binPy[binNo] += mass*vel.y();
1832 >        binPz[binNo] += mass*vel.z();
1833 >        binOmegax[binNo] += aVel.x();
1834 >        binOmegay[binNo] += aVel.y();
1835 >        binOmegaz[binNo] += aVel.z();
1836 >        binKE[binNo] += 0.5 * (mass * vel.lengthSquare());
1837 >        binDOF[binNo] += 3;
1838 >        
1839 >        if (sd->isDirectional()) {
1840 >          Vector3d angMom = sd->getJ();
1841 >          Mat3x3d I = sd->getI();
1842 >          if (sd->isLinear()) {
1843 >            int i = sd->linearAxis();
1844 >            int j = (i + 1) % 3;
1845 >            int k = (i + 2) % 3;
1846 >            binKE[binNo] += 0.5 * (angMom[j] * angMom[j] / I(j, j) +
1847 >                                   angMom[k] * angMom[k] / I(k, k));
1848 >            binDOF[binNo] += 2;
1849 >          } else {
1850 >            binKE[binNo] += 0.5 * (angMom[0] * angMom[0] / I(0, 0) +
1851 >                                   angMom[1] * angMom[1] / I(1, 1) +
1852 >                                   angMom[2] * angMom[2] / I(2, 2));
1853 >            binDOF[binNo] += 3;
1854 >          }
1855          }
1856        }
1857      }
# Line 1502 | Line 1867 | namespace OpenMD {
1867                                nBins_, MPI::REALTYPE, MPI::SUM);
1868      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binPz[0],
1869                                nBins_, MPI::REALTYPE, MPI::SUM);
1870 +    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binOmegax[0],
1871 +                              nBins_, MPI::REALTYPE, MPI::SUM);
1872 +    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binOmegay[0],
1873 +                              nBins_, MPI::REALTYPE, MPI::SUM);
1874 +    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binOmegaz[0],
1875 +                              nBins_, MPI::REALTYPE, MPI::SUM);
1876      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binKE[0],
1877                                nBins_, MPI::REALTYPE, MPI::SUM);
1878      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binDOF[0],
# Line 1509 | Line 1880 | namespace OpenMD {
1880   #endif
1881  
1882      Vector3d vel;
1883 +    Vector3d aVel;
1884      RealType den;
1885      RealType temp;
1886      RealType z;
1887 +    RealType r;
1888      for (int i = 0; i < nBins_; i++) {
1889 <      z = (((RealType)i + 0.5) / (RealType)nBins_) * hmat(2,2);
1889 >      if (usePeriodicBoundaryConditions_) {
1890 >        z = (((RealType)i + 0.5) / (RealType)nBins_) * hmat(2,2);
1891 >        den = binMass[i] * nBins_ * PhysicalConstants::densityConvert
1892 >          / currentSnap_->getVolume() ;
1893 >      } else {
1894 >        r = (((RealType)i + 0.5) * binWidth_);
1895 >        RealType rinner = (RealType)i * binWidth_;
1896 >        RealType router = (RealType)(i+1) * binWidth_;
1897 >        den = binMass[i] * 3.0 * PhysicalConstants::densityConvert
1898 >          / (4.0 * M_PI * (pow(router,3) - pow(rinner,3)));
1899 >      }
1900        vel.x() = binPx[i] / binMass[i];
1901        vel.y() = binPy[i] / binMass[i];
1902        vel.z() = binPz[i] / binMass[i];
1903 +      aVel.x() = binOmegax[i];
1904 +      aVel.y() = binOmegay[i];
1905 +      aVel.z() = binOmegaz[i];
1906  
1907 <      den = binMass[i] * nBins_ * PhysicalConstants::densityConvert
1908 <        / currentSnap_->getVolume() ;
1909 <
1910 <      temp = 2.0 * binKE[i] / (binDOF[i] * PhysicalConstants::kb *
1911 <                               PhysicalConstants::energyConvert);
1912 <  
1913 <      for (unsigned int j = 0; j < outputMask_.size(); ++j) {
1914 <        if(outputMask_[j]) {
1915 <          switch(j) {
1916 <          case Z:
1917 <            (data_[j].accumulator[i])->add(z);
1918 <            break;
1919 <          case TEMPERATURE:
1920 <            data_[j].accumulator[i]->add(temp);
1921 <            break;
1922 <          case VELOCITY:
1923 <            dynamic_cast<VectorAccumulator *>(data_[j].accumulator[i])->add(vel);
1924 <            break;
1925 <          case DENSITY:
1926 <            data_[j].accumulator[i]->add(den);
1927 <            break;
1907 >      if (binCount[i] > 0) {
1908 >        // only add values if there are things to add
1909 >        temp = 2.0 * binKE[i] / (binDOF[i] * PhysicalConstants::kb *
1910 >                                 PhysicalConstants::energyConvert);
1911 >        
1912 >        for (unsigned int j = 0; j < outputMask_.size(); ++j) {
1913 >          if(outputMask_[j]) {
1914 >            switch(j) {
1915 >            case Z:
1916 >              dynamic_cast<Accumulator *>(data_[j].accumulator[i])->add(z);
1917 >              break;
1918 >            case R:
1919 >              dynamic_cast<Accumulator *>(data_[j].accumulator[i])->add(r);
1920 >              break;
1921 >            case TEMPERATURE:
1922 >              dynamic_cast<Accumulator *>(data_[j].accumulator[i])->add(temp);
1923 >              break;
1924 >            case VELOCITY:
1925 >              dynamic_cast<VectorAccumulator *>(data_[j].accumulator[i])->add(vel);
1926 >              break;
1927 >            case ANGULARVELOCITY:  
1928 >              dynamic_cast<VectorAccumulator *>(data_[j].accumulator[i])->add(aVel);
1929 >              break;
1930 >            case DENSITY:
1931 >              dynamic_cast<Accumulator *>(data_[j].accumulator[i])->add(den);
1932 >              break;
1933 >            }
1934            }
1935          }
1936        }
# Line 1547 | Line 1939 | namespace OpenMD {
1939  
1940    void RNEMD::getStarted() {
1941      if (!doRNEMD_) return;
1942 +    hasDividingArea_ = false;
1943      collectData();
1944      writeOutputFile();
1945    }
# Line 1595 | Line 1988 | namespace OpenMD {
1988        RealType time = currentSnap_->getTime();
1989        RealType avgArea;
1990        areaAccumulator_->getAverage(avgArea);
1991 <      RealType Jz = kineticExchange_ / (2.0 * time * avgArea)
1991 >      RealType Jz = kineticExchange_ / (time * avgArea)
1992          / PhysicalConstants::energyConvert;
1993 <      Vector3d JzP = momentumExchange_ / (2.0 * time * avgArea);      
1993 >      Vector3d JzP = momentumExchange_ / (time * avgArea);      
1994 >      Vector3d JzL = angularMomentumExchange_ / (time * avgArea);      
1995  
1996        rnemdFile_ << "#######################################################\n";
1997        rnemdFile_ << "# RNEMD {\n";
# Line 1617 | Line 2011 | namespace OpenMD {
2011  
2012        rnemdFile_ << "#    objectSelection = \""
2013                   << rnemdObjectSelection_ << "\";\n";
2014 <      rnemdFile_ << "#    slabWidth = " << slabWidth_ << ";\n";
2015 <      rnemdFile_ << "#    slabAcenter = " << slabACenter_ << ";\n";
1622 <      rnemdFile_ << "#    slabBcenter = " << slabBCenter_ << ";\n";
2014 >      rnemdFile_ << "#    selectionA = \"" << selectionA_ << "\";\n";
2015 >      rnemdFile_ << "#    selectionB = \"" << selectionB_ << "\";\n";
2016        rnemdFile_ << "# }\n";
2017        rnemdFile_ << "#######################################################\n";
2018        rnemdFile_ << "# RNEMD report:\n";      
2019 <      rnemdFile_ << "#     running time = " << time << " fs\n";
2020 <      rnemdFile_ << "#     target flux:\n";
2021 <      rnemdFile_ << "#         kinetic = "
2019 >      rnemdFile_ << "#      running time = " << time << " fs\n";
2020 >      rnemdFile_ << "# Target flux:\n";
2021 >      rnemdFile_ << "#           kinetic = "
2022                   << kineticFlux_ / PhysicalConstants::energyConvert
2023                   << " (kcal/mol/A^2/fs)\n";
2024 <      rnemdFile_ << "#         momentum = " << momentumFluxVector_
2024 >      rnemdFile_ << "#          momentum = " << momentumFluxVector_
2025                   << " (amu/A/fs^2)\n";
2026 <      rnemdFile_ << "#     target one-time exchanges:\n";
2027 <      rnemdFile_ << "#         kinetic = "
2026 >      rnemdFile_ << "#  angular momentum = " << angularMomentumFluxVector_
2027 >                 << " (amu/A^2/fs^2)\n";
2028 >      rnemdFile_ << "# Target one-time exchanges:\n";
2029 >      rnemdFile_ << "#          kinetic = "
2030                   << kineticTarget_ / PhysicalConstants::energyConvert
2031                   << " (kcal/mol)\n";
2032 <      rnemdFile_ << "#         momentum = " << momentumTarget_
2032 >      rnemdFile_ << "#          momentum = " << momentumTarget_
2033                   << " (amu*A/fs)\n";
2034 <      rnemdFile_ << "#     actual exchange totals:\n";
2035 <      rnemdFile_ << "#         kinetic = "
2034 >      rnemdFile_ << "#  angular momentum = " << angularMomentumTarget_
2035 >                 << " (amu*A^2/fs)\n";
2036 >      rnemdFile_ << "# Actual exchange totals:\n";
2037 >      rnemdFile_ << "#          kinetic = "
2038                   << kineticExchange_ / PhysicalConstants::energyConvert
2039                   << " (kcal/mol)\n";
2040 <      rnemdFile_ << "#         momentum = " << momentumExchange_
2040 >      rnemdFile_ << "#          momentum = " << momentumExchange_
2041                   << " (amu*A/fs)\n";      
2042 <      rnemdFile_ << "#     actual flux:\n";
2043 <      rnemdFile_ << "#         kinetic = " << Jz
2042 >      rnemdFile_ << "#  angular momentum = " << angularMomentumExchange_
2043 >                 << " (amu*A^2/fs)\n";      
2044 >      rnemdFile_ << "# Actual flux:\n";
2045 >      rnemdFile_ << "#          kinetic = " << Jz
2046                   << " (kcal/mol/A^2/fs)\n";
2047 <      rnemdFile_ << "#         momentum = " << JzP
2047 >      rnemdFile_ << "#          momentum = " << JzP
2048                   << " (amu/A/fs^2)\n";
2049 <      rnemdFile_ << "#     exchange statistics:\n";
2050 <      rnemdFile_ << "#         attempted = " << trialCount_ << "\n";
2051 <      rnemdFile_ << "#         failed = " << failTrialCount_ << "\n";    
2049 >      rnemdFile_ << "#  angular momentum = " << JzL
2050 >                 << " (amu/A^2/fs^2)\n";
2051 >      rnemdFile_ << "# Exchange statistics:\n";
2052 >      rnemdFile_ << "#               attempted = " << trialCount_ << "\n";
2053 >      rnemdFile_ << "#                  failed = " << failTrialCount_ << "\n";
2054        if (rnemdMethod_ == rnemdNIVS) {
2055 <        rnemdFile_ << "#         NIVS root-check errors = "
2055 >        rnemdFile_ << "#  NIVS root-check errors = "
2056                     << failRootCount_ << "\n";
2057        }
2058        rnemdFile_ << "#######################################################\n";
# Line 1672 | Line 2073 | namespace OpenMD {
2073        
2074        rnemdFile_.precision(8);
2075        
2076 <      for (unsigned int j = 0; j < nBins_; j++) {        
2076 >      for (int j = 0; j < nBins_; j++) {        
2077          
2078          for (unsigned int i = 0; i < outputMask_.size(); ++i) {
2079            if (outputMask_[i]) {
2080              if (data_[i].dataType == "RealType")
2081                writeReal(i,j);
2082 <            else if (data_[i].dataType == "Vector3d")
2082 >            else if (data_[i].dataType == "Vector3d")
2083                writeVector(i,j);
2084              else {
2085                sprintf( painCave.errMsg,
# Line 1698 | Line 2099 | namespace OpenMD {
2099        rnemdFile_ << "#######################################################\n";
2100  
2101  
2102 <      for (unsigned int j = 0; j < nBins_; j++) {        
2102 >      for (int j = 0; j < nBins_; j++) {        
2103          rnemdFile_ << "#";
2104          for (unsigned int i = 0; i < outputMask_.size(); ++i) {
2105            if (outputMask_[i]) {
# Line 1731 | Line 2132 | namespace OpenMD {
2132    void RNEMD::writeReal(int index, unsigned int bin) {
2133      if (!doRNEMD_) return;
2134      assert(index >=0 && index < ENDINDEX);
2135 <    assert(bin < nBins_);
2135 >    assert(int(bin) < nBins_);
2136      RealType s;
2137 +    int count;
2138      
2139 <    data_[index].accumulator[bin]->getAverage(s);
2139 >    count = data_[index].accumulator[bin]->count();
2140 >    if (count == 0) return;
2141      
2142 +    dynamic_cast<Accumulator *>(data_[index].accumulator[bin])->getAverage(s);
2143 +    
2144      if (! isinf(s) && ! isnan(s)) {
2145        rnemdFile_ << "\t" << s;
2146      } else{
# Line 1750 | Line 2155 | namespace OpenMD {
2155    void RNEMD::writeVector(int index, unsigned int bin) {
2156      if (!doRNEMD_) return;
2157      assert(index >=0 && index < ENDINDEX);
2158 <    assert(bin < nBins_);
2158 >    assert(int(bin) < nBins_);
2159      Vector3d s;
2160 +    int count;
2161 +    
2162 +    count = data_[index].accumulator[bin]->count();
2163 +
2164 +    if (count == 0) return;
2165 +
2166      dynamic_cast<VectorAccumulator*>(data_[index].accumulator[bin])->getAverage(s);
2167      if (isinf(s[0]) || isnan(s[0]) ||
2168          isinf(s[1]) || isnan(s[1]) ||
# Line 1769 | Line 2180 | namespace OpenMD {
2180    void RNEMD::writeRealStdDev(int index, unsigned int bin) {
2181      if (!doRNEMD_) return;
2182      assert(index >=0 && index < ENDINDEX);
2183 <    assert(bin < nBins_);
2183 >    assert(int(bin) < nBins_);
2184      RealType s;
2185 +    int count;
2186      
2187 <    data_[index].accumulator[bin]->getStdDev(s);
2187 >    count = data_[index].accumulator[bin]->count();
2188 >    if (count == 0) return;
2189      
2190 +    dynamic_cast<Accumulator *>(data_[index].accumulator[bin])->getStdDev(s);
2191 +    
2192      if (! isinf(s) && ! isnan(s)) {
2193        rnemdFile_ << "\t" << s;
2194      } else{
# Line 1788 | Line 2203 | namespace OpenMD {
2203    void RNEMD::writeVectorStdDev(int index, unsigned int bin) {
2204      if (!doRNEMD_) return;
2205      assert(index >=0 && index < ENDINDEX);
2206 <    assert(bin < nBins_);
2206 >    assert(int(bin) < nBins_);
2207      Vector3d s;
2208 +    int count;
2209 +    
2210 +    count = data_[index].accumulator[bin]->count();
2211 +    if (count == 0) return;
2212 +
2213      dynamic_cast<VectorAccumulator*>(data_[index].accumulator[bin])->getStdDev(s);
2214      if (isinf(s[0]) || isnan(s[0]) ||
2215          isinf(s[1]) || isnan(s[1]) ||

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines