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

Comparing trunk/src/rnemd/RNEMD.cpp (file contents):
Revision 1903 by gezelter, Tue Jul 16 18:58:08 2013 UTC vs.
Revision 2068 by gezelter, Thu Mar 5 15:40:58 2015 UTC

# Line 38 | Line 38
38   * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).          
39   * [4]  Vardeman & Gezelter, in progress (2009).                        
40   */
41 + #ifdef IS_MPI
42 + #include <mpi.h>
43 + #endif
44  
45   #include <cmath>
46   #include <sstream>
# Line 54 | Line 57
57   #include "utils/Tuple.hpp"
58   #include "brains/Thermo.hpp"
59   #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))
# Line 291 | Line 291 | namespace OpenMD {
291        kineticFlux_ = 0.0;
292      }
293      if (hasMomentumFluxVector) {
294 <      momentumFluxVector_ = rnemdParams->getMomentumFluxVector();
294 >      std::vector<RealType> mf = rnemdParams->getMomentumFluxVector();
295 >      if (mf.size() != 3) {
296 >          sprintf(painCave.errMsg,
297 >                  "RNEMD: Incorrect number of parameters specified for momentumFluxVector.\n"
298 >                  "\tthere should be 3 parameters, but %lu were specified.\n",
299 >                  mf.size());
300 >          painCave.isFatal = 1;
301 >          simError();      
302 >      }
303 >      momentumFluxVector_.x() = mf[0];
304 >      momentumFluxVector_.y() = mf[1];
305 >      momentumFluxVector_.z() = mf[2];
306      } else {
307        momentumFluxVector_ = V3Zero;
308        if (hasMomentumFlux) {
# Line 316 | Line 327 | namespace OpenMD {
327            break;
328          }
329        }
330 <      if (hasAngularMomentumFluxVector) {
331 <        angularMomentumFluxVector_ = rnemdParams->getAngularMomentumFluxVector();
332 <      } else {
333 <        angularMomentumFluxVector_ = V3Zero;
334 <        if (hasAngularMomentumFlux) {
335 <          RealType angularMomentumFlux = rnemdParams->getAngularMomentumFlux();
336 <          switch (rnemdFluxType_) {
337 <          case rnemdLx:
338 <            angularMomentumFluxVector_.x() = angularMomentumFlux;
339 <            break;
329 <          case rnemdLy:
330 <            angularMomentumFluxVector_.y() = angularMomentumFlux;
331 <            break;
332 <          case rnemdLz:
333 <            angularMomentumFluxVector_.z() = angularMomentumFlux;
334 <            break;
335 <          case rnemdKeLx:
336 <            angularMomentumFluxVector_.x() = angularMomentumFlux;
337 <            break;
338 <          case rnemdKeLy:
339 <            angularMomentumFluxVector_.y() = angularMomentumFlux;
340 <            break;
341 <          case rnemdKeLz:
342 <            angularMomentumFluxVector_.z() = angularMomentumFlux;
343 <            break;
344 <          default:
345 <            break;
346 <          }
347 <        }        
348 <      }
349 <
350 <      if (hasCoordinateOrigin) {
351 <        coordinateOrigin_ = rnemdParams->getCoordinateOrigin();
352 <      } else {
353 <        coordinateOrigin_ = V3Zero;
330 >    }
331 >    if (hasAngularMomentumFluxVector) {
332 >      std::vector<RealType> amf = rnemdParams->getAngularMomentumFluxVector();
333 >      if (amf.size() != 3) {
334 >        sprintf(painCave.errMsg,
335 >                "RNEMD: Incorrect number of parameters specified for angularMomentumFluxVector.\n"
336 >                "\tthere should be 3 parameters, but %lu were specified.\n",
337 >                amf.size());
338 >        painCave.isFatal = 1;
339 >        simError();      
340        }
341 <
342 <      // do some sanity checking
343 <
344 <      int selectionCount = seleMan_.getSelectionCount();
345 <
346 <      int nIntegrable = info->getNGlobalIntegrableObjects();
347 <
362 <      if (selectionCount > nIntegrable) {
363 <        sprintf(painCave.errMsg,
364 <                "RNEMD: The current objectSelection,\n"
365 <                "\t\t%s\n"
366 <                "\thas resulted in %d selected objects.  However,\n"
367 <                "\tthe total number of integrable objects in the system\n"
368 <                "\tis only %d.  This is almost certainly not what you want\n"
369 <                "\tto do.  A likely cause of this is forgetting the _RB_0\n"
370 <                "\tselector in the selection script!\n",
371 <                rnemdObjectSelection_.c_str(),
372 <                selectionCount, nIntegrable);
373 <        painCave.isFatal = 0;
374 <        painCave.severity = OPENMD_WARNING;
375 <        simError();
376 <      }
377 <
378 <      areaAccumulator_ = new Accumulator();
379 <
380 <      nBins_ = rnemdParams->getOutputBins();
381 <      binWidth_ = rnemdParams->getOutputBinWidth();
382 <
383 <      data_.resize(RNEMD::ENDINDEX);
384 <      OutputData z;
385 <      z.units =  "Angstroms";
386 <      z.title =  "Z";
387 <      z.dataType = "RealType";
388 <      z.accumulator.reserve(nBins_);
389 <      for (int i = 0; i < nBins_; i++)
390 <        z.accumulator.push_back( new Accumulator() );
391 <      data_[Z] = z;
392 <      outputMap_["Z"] =  Z;
393 <
394 <      OutputData r;
395 <      r.units =  "Angstroms";
396 <      r.title =  "R";
397 <      r.dataType = "RealType";
398 <      r.accumulator.reserve(nBins_);
399 <      for (int i = 0; i < nBins_; i++)
400 <        r.accumulator.push_back( new Accumulator() );
401 <      data_[R] = r;
402 <      outputMap_["R"] =  R;
403 <
404 <      OutputData temperature;
405 <      temperature.units =  "K";
406 <      temperature.title =  "Temperature";
407 <      temperature.dataType = "RealType";
408 <      temperature.accumulator.reserve(nBins_);
409 <      for (int i = 0; i < nBins_; i++)
410 <        temperature.accumulator.push_back( new Accumulator() );
411 <      data_[TEMPERATURE] = temperature;
412 <      outputMap_["TEMPERATURE"] =  TEMPERATURE;
413 <
414 <      OutputData velocity;
415 <      velocity.units = "angstroms/fs";
416 <      velocity.title =  "Velocity";  
417 <      velocity.dataType = "Vector3d";
418 <      velocity.accumulator.reserve(nBins_);
419 <      for (int i = 0; i < nBins_; i++)
420 <        velocity.accumulator.push_back( new VectorAccumulator() );
421 <      data_[VELOCITY] = velocity;
422 <      outputMap_["VELOCITY"] = VELOCITY;
423 <
424 <      OutputData angularVelocity;
425 <      angularVelocity.units = "angstroms^2/fs";
426 <      angularVelocity.title =  "AngularVelocity";  
427 <      angularVelocity.dataType = "Vector3d";
428 <      angularVelocity.accumulator.reserve(nBins_);
429 <      for (int i = 0; i < nBins_; i++)
430 <        angularVelocity.accumulator.push_back( new VectorAccumulator() );
431 <      data_[ANGULARVELOCITY] = angularVelocity;
432 <      outputMap_["ANGULARVELOCITY"] = ANGULARVELOCITY;
433 <
434 <      OutputData density;
435 <      density.units =  "g cm^-3";
436 <      density.title =  "Density";
437 <      density.dataType = "RealType";
438 <      density.accumulator.reserve(nBins_);
439 <      for (int i = 0; i < nBins_; i++)
440 <        density.accumulator.push_back( new Accumulator() );
441 <      data_[DENSITY] = density;
442 <      outputMap_["DENSITY"] =  DENSITY;
443 <
444 <      if (hasOutputFields) {
445 <        parseOutputFileFormat(rnemdParams->getOutputFields());
446 <      } else {
447 <        if (usePeriodicBoundaryConditions_)
448 <          outputMask_.set(Z);
449 <        else
450 <          outputMask_.set(R);
341 >      angularMomentumFluxVector_.x() = amf[0];
342 >      angularMomentumFluxVector_.y() = amf[1];
343 >      angularMomentumFluxVector_.z() = amf[2];
344 >    } else {
345 >      angularMomentumFluxVector_ = V3Zero;
346 >      if (hasAngularMomentumFlux) {
347 >        RealType angularMomentumFlux = rnemdParams->getAngularMomentumFlux();
348          switch (rnemdFluxType_) {
452 –        case rnemdKE:
453 –        case rnemdRotKE:
454 –        case rnemdFullKE:
455 –          outputMask_.set(TEMPERATURE);
456 –          break;
457 –        case rnemdPx:
458 –        case rnemdPy:
459 –          outputMask_.set(VELOCITY);
460 –          break;
461 –        case rnemdPz:        
462 –        case rnemdPvector:
463 –          outputMask_.set(VELOCITY);
464 –          outputMask_.set(DENSITY);
465 –          break;
349          case rnemdLx:
350 +          angularMomentumFluxVector_.x() = angularMomentumFlux;
351 +          break;
352          case rnemdLy:
353 +          angularMomentumFluxVector_.y() = angularMomentumFlux;
354 +          break;
355          case rnemdLz:
356 <        case rnemdLvector:
470 <          outputMask_.set(ANGULARVELOCITY);
356 >          angularMomentumFluxVector_.z() = angularMomentumFlux;
357            break;
358          case rnemdKeLx:
359 +          angularMomentumFluxVector_.x() = angularMomentumFlux;
360 +          break;
361          case rnemdKeLy:
362 +          angularMomentumFluxVector_.y() = angularMomentumFlux;
363 +          break;
364          case rnemdKeLz:
365 <        case rnemdKeLvector:
476 <          outputMask_.set(TEMPERATURE);
477 <          outputMask_.set(ANGULARVELOCITY);
365 >          angularMomentumFluxVector_.z() = angularMomentumFlux;
366            break;
479 –        case rnemdKePx:
480 –        case rnemdKePy:
481 –          outputMask_.set(TEMPERATURE);
482 –          outputMask_.set(VELOCITY);
483 –          break;
484 –        case rnemdKePvector:
485 –          outputMask_.set(TEMPERATURE);
486 –          outputMask_.set(VELOCITY);
487 –          outputMask_.set(DENSITY);        
488 –          break;
367          default:
368            break;
369          }
370 +      }        
371 +    }
372 +
373 +    if (hasCoordinateOrigin) {
374 +      std::vector<RealType> co = rnemdParams->getCoordinateOrigin();
375 +      if (co.size() != 3) {
376 +        sprintf(painCave.errMsg,
377 +                "RNEMD: Incorrect number of parameters specified for coordinateOrigin.\n"
378 +                "\tthere should be 3 parameters, but %lu were specified.\n",
379 +                co.size());
380 +        painCave.isFatal = 1;
381 +        simError();      
382        }
383 <      
384 <      if (hasOutputFileName) {
385 <        rnemdFileName_ = rnemdParams->getOutputFileName();
386 <      } else {
387 <        rnemdFileName_ = getPrefix(info->getFinalConfigFileName()) + ".rnemd";
388 <      }          
499 <
500 <      exchangeTime_ = rnemdParams->getExchangeTime();
501 <
502 <      Snapshot* currentSnap_ = info->getSnapshotManager()->getCurrentSnapshot();
503 <      // total exchange sums are zeroed out at the beginning:
504 <
505 <      kineticExchange_ = 0.0;
506 <      momentumExchange_ = V3Zero;
507 <      angularMomentumExchange_ = V3Zero;
508 <
509 <      std::ostringstream selectionAstream;
510 <      std::ostringstream selectionBstream;
383 >      coordinateOrigin_.x() = co[0];
384 >      coordinateOrigin_.y() = co[1];
385 >      coordinateOrigin_.z() = co[2];
386 >    } else {
387 >      coordinateOrigin_ = V3Zero;
388 >    }
389      
390 <      if (hasSelectionA_) {
391 <        selectionA_ = rnemdParams->getSelectionA();
392 <      } else {
393 <        if (usePeriodicBoundaryConditions_) {    
394 <          Mat3x3d hmat = currentSnap_->getHmat();
395 <        
396 <          if (hasSlabWidth)
397 <            slabWidth_ = rnemdParams->getSlabWidth();
398 <          else
399 <            slabWidth_ = hmat(2,2) / 10.0;
400 <        
401 <          if (hasSlabACenter)
402 <            slabACenter_ = rnemdParams->getSlabACenter();
403 <          else
404 <            slabACenter_ = 0.0;
405 <        
406 <          selectionAstream << "select wrappedz > "
407 <                           << slabACenter_ - 0.5*slabWidth_
408 <                           <<  " && wrappedz < "
409 <                           << slabACenter_ + 0.5*slabWidth_;
410 <          selectionA_ = selectionAstream.str();
411 <        } else {
412 <          if (hasSphereARadius)
413 <            sphereARadius_ = rnemdParams->getSphereARadius();
414 <          else {
415 <            // use an initial guess to the size of the inner slab to be 1/10 the
416 <            // radius of an approximately spherical hull:
417 <            Thermo thermo(info);
418 <            RealType hVol = thermo.getHullVolume();
419 <            sphereARadius_ = 0.1 * pow((3.0 * hVol / (4.0 * M_PI)), 1.0/3.0);
420 <          }
421 <          selectionAstream << "select r < " << sphereARadius_;
422 <          selectionA_ = selectionAstream.str();
423 <        }
390 >    // do some sanity checking
391 >    
392 >    int selectionCount = seleMan_.getSelectionCount();    
393 >    int nIntegrable = info->getNGlobalIntegrableObjects();
394 >    if (selectionCount > nIntegrable) {
395 >      sprintf(painCave.errMsg,
396 >              "RNEMD: The current objectSelection,\n"
397 >              "\t\t%s\n"
398 >              "\thas resulted in %d selected objects.  However,\n"
399 >              "\tthe total number of integrable objects in the system\n"
400 >              "\tis only %d.  This is almost certainly not what you want\n"
401 >              "\tto do.  A likely cause of this is forgetting the _RB_0\n"
402 >              "\tselector in the selection script!\n",
403 >              rnemdObjectSelection_.c_str(),
404 >              selectionCount, nIntegrable);
405 >      painCave.isFatal = 0;
406 >      painCave.severity = OPENMD_WARNING;
407 >      simError();
408 >    }
409 >    
410 >    areaAccumulator_ = new Accumulator();
411 >    
412 >    nBins_ = rnemdParams->getOutputBins();
413 >    binWidth_ = rnemdParams->getOutputBinWidth();
414 >    
415 >    data_.resize(RNEMD::ENDINDEX);
416 >    OutputData z;
417 >    z.units =  "Angstroms";
418 >    z.title =  "Z";
419 >    z.dataType = "RealType";
420 >    z.accumulator.reserve(nBins_);
421 >    for (int i = 0; i < nBins_; i++)
422 >      z.accumulator.push_back( new Accumulator() );
423 >    data_[Z] = z;
424 >    outputMap_["Z"] =  Z;
425 >    
426 >    OutputData r;
427 >    r.units =  "Angstroms";
428 >    r.title =  "R";
429 >    r.dataType = "RealType";
430 >    r.accumulator.reserve(nBins_);
431 >    for (int i = 0; i < nBins_; i++)
432 >      r.accumulator.push_back( new Accumulator() );
433 >    data_[R] = r;
434 >    outputMap_["R"] =  R;
435 >    
436 >    OutputData temperature;
437 >    temperature.units =  "K";
438 >    temperature.title =  "Temperature";
439 >    temperature.dataType = "RealType";
440 >    temperature.accumulator.reserve(nBins_);
441 >    for (int i = 0; i < nBins_; i++)
442 >      temperature.accumulator.push_back( new Accumulator() );
443 >    data_[TEMPERATURE] = temperature;
444 >    outputMap_["TEMPERATURE"] =  TEMPERATURE;
445 >    
446 >    OutputData velocity;
447 >    velocity.units = "angstroms/fs";
448 >    velocity.title =  "Velocity";  
449 >    velocity.dataType = "Vector3d";
450 >    velocity.accumulator.reserve(nBins_);
451 >    for (int i = 0; i < nBins_; i++)
452 >      velocity.accumulator.push_back( new VectorAccumulator() );
453 >    data_[VELOCITY] = velocity;
454 >    outputMap_["VELOCITY"] = VELOCITY;
455 >    
456 >    OutputData angularVelocity;
457 >    angularVelocity.units = "angstroms^2/fs";
458 >    angularVelocity.title =  "AngularVelocity";  
459 >    angularVelocity.dataType = "Vector3d";
460 >    angularVelocity.accumulator.reserve(nBins_);
461 >    for (int i = 0; i < nBins_; i++)
462 >      angularVelocity.accumulator.push_back( new VectorAccumulator() );
463 >    data_[ANGULARVELOCITY] = angularVelocity;
464 >    outputMap_["ANGULARVELOCITY"] = ANGULARVELOCITY;
465 >    
466 >    OutputData density;
467 >    density.units =  "g cm^-3";
468 >    density.title =  "Density";
469 >    density.dataType = "RealType";
470 >    density.accumulator.reserve(nBins_);
471 >    for (int i = 0; i < nBins_; i++)
472 >      density.accumulator.push_back( new Accumulator() );
473 >    data_[DENSITY] = density;
474 >    outputMap_["DENSITY"] =  DENSITY;
475 >    
476 >    if (hasOutputFields) {
477 >      parseOutputFileFormat(rnemdParams->getOutputFields());
478 >    } else {
479 >      if (usePeriodicBoundaryConditions_)
480 >        outputMask_.set(Z);
481 >      else
482 >        outputMask_.set(R);
483 >      switch (rnemdFluxType_) {
484 >      case rnemdKE:
485 >      case rnemdRotKE:
486 >      case rnemdFullKE:
487 >        outputMask_.set(TEMPERATURE);
488 >        break;
489 >      case rnemdPx:
490 >      case rnemdPy:
491 >        outputMask_.set(VELOCITY);
492 >        break;
493 >      case rnemdPz:        
494 >      case rnemdPvector:
495 >        outputMask_.set(VELOCITY);
496 >        outputMask_.set(DENSITY);
497 >        break;
498 >      case rnemdLx:
499 >      case rnemdLy:
500 >      case rnemdLz:
501 >      case rnemdLvector:
502 >        outputMask_.set(ANGULARVELOCITY);
503 >        break;
504 >      case rnemdKeLx:
505 >      case rnemdKeLy:
506 >      case rnemdKeLz:
507 >      case rnemdKeLvector:
508 >        outputMask_.set(TEMPERATURE);
509 >        outputMask_.set(ANGULARVELOCITY);
510 >        break;
511 >      case rnemdKePx:
512 >      case rnemdKePy:
513 >        outputMask_.set(TEMPERATURE);
514 >        outputMask_.set(VELOCITY);
515 >        break;
516 >      case rnemdKePvector:
517 >        outputMask_.set(TEMPERATURE);
518 >        outputMask_.set(VELOCITY);
519 >        outputMask_.set(DENSITY);        
520 >        break;
521 >      default:
522 >        break;
523        }
524 +    }
525      
526 <      if (hasSelectionB_) {
527 <        selectionB_ = rnemdParams->getSelectionB();
528 <
526 >    if (hasOutputFileName) {
527 >      rnemdFileName_ = rnemdParams->getOutputFileName();
528 >    } else {
529 >      rnemdFileName_ = getPrefix(info->getFinalConfigFileName()) + ".rnemd";
530 >    }          
531 >    
532 >    exchangeTime_ = rnemdParams->getExchangeTime();
533 >    
534 >    Snapshot* currentSnap_ = info->getSnapshotManager()->getCurrentSnapshot();
535 >    // total exchange sums are zeroed out at the beginning:
536 >    
537 >    kineticExchange_ = 0.0;
538 >    momentumExchange_ = V3Zero;
539 >    angularMomentumExchange_ = V3Zero;
540 >    
541 >    std::ostringstream selectionAstream;
542 >    std::ostringstream selectionBstream;
543 >    
544 >    if (hasSelectionA_) {
545 >      selectionA_ = rnemdParams->getSelectionA();
546 >    } else {
547 >      if (usePeriodicBoundaryConditions_) {    
548 >        Mat3x3d hmat = currentSnap_->getHmat();
549 >        
550 >        if (hasSlabWidth)
551 >          slabWidth_ = rnemdParams->getSlabWidth();
552 >        else
553 >          slabWidth_ = hmat(2,2) / 10.0;
554 >        
555 >        if (hasSlabACenter)
556 >          slabACenter_ = rnemdParams->getSlabACenter();
557 >        else
558 >          slabACenter_ = 0.0;
559 >        
560 >        selectionAstream << "select wrappedz > "
561 >                         << slabACenter_ - 0.5*slabWidth_
562 >                         <<  " && wrappedz < "
563 >                         << slabACenter_ + 0.5*slabWidth_;
564 >        selectionA_ = selectionAstream.str();
565        } else {
566 <        if (usePeriodicBoundaryConditions_) {    
567 <          Mat3x3d hmat = currentSnap_->getHmat();
566 >        if (hasSphereARadius)
567 >          sphereARadius_ = rnemdParams->getSphereARadius();
568 >        else {
569 >          // use an initial guess to the size of the inner slab to be 1/10 the
570 >          // radius of an approximately spherical hull:
571 >          Thermo thermo(info);
572 >          RealType hVol = thermo.getHullVolume();
573 >          sphereARadius_ = 0.1 * pow((3.0 * hVol / (4.0 * M_PI)), 1.0/3.0);
574 >        }
575 >        selectionAstream << "select r < " << sphereARadius_;
576 >        selectionA_ = selectionAstream.str();
577 >      }
578 >    }
579 >    
580 >    if (hasSelectionB_) {
581 >      selectionB_ = rnemdParams->getSelectionB();
582 >      
583 >    } else {
584 >      if (usePeriodicBoundaryConditions_) {    
585 >        Mat3x3d hmat = currentSnap_->getHmat();
586          
587 <          if (hasSlabWidth)
588 <            slabWidth_ = rnemdParams->getSlabWidth();
589 <          else
590 <            slabWidth_ = hmat(2,2) / 10.0;
587 >        if (hasSlabWidth)
588 >          slabWidth_ = rnemdParams->getSlabWidth();
589 >        else
590 >          slabWidth_ = hmat(2,2) / 10.0;
591          
592 <          if (hasSlabBCenter)
593 <            slabBCenter_ = rnemdParams->getSlabBCenter();
594 <          else
595 <            slabBCenter_ = hmat(2,2) / 2.0;
592 >        if (hasSlabBCenter)
593 >          slabBCenter_ = rnemdParams->getSlabBCenter();
594 >        else
595 >          slabBCenter_ = hmat(2,2) / 2.0;
596          
597 <          selectionBstream << "select wrappedz > "
598 <                           << slabBCenter_ - 0.5*slabWidth_
599 <                           <<  " && wrappedz < "
600 <                           << slabBCenter_ + 0.5*slabWidth_;
597 >        selectionBstream << "select wrappedz > "
598 >                         << slabBCenter_ - 0.5*slabWidth_
599 >                         <<  " && wrappedz < "
600 >                         << slabBCenter_ + 0.5*slabWidth_;
601 >        selectionB_ = selectionBstream.str();
602 >      } else {
603 >        if (hasSphereBRadius_) {
604 >          sphereBRadius_ = rnemdParams->getSphereBRadius();
605 >          selectionBstream << "select r > " << sphereBRadius_;
606            selectionB_ = selectionBstream.str();
607          } else {
608 <          if (hasSphereBRadius_) {
609 <            sphereBRadius_ = rnemdParams->getSphereBRadius();
610 <            selectionBstream << "select r > " << sphereBRadius_;
574 <            selectionB_ = selectionBstream.str();
575 <          } else {
576 <            selectionB_ = "select hull";
577 <            BisHull_ = true;
578 <            hasSelectionB_ = true;
579 <          }
608 >          selectionB_ = "select hull";
609 >          BisHull_ = true;
610 >          hasSelectionB_ = true;
611          }
612        }
613      }
614 <
614 >  
615 >  
616      // object evaluator:
617      evaluator_.loadScriptString(rnemdObjectSelection_);
618      seleMan_.setSelectionSet(evaluator_.evaluate());
# Line 622 | Line 654 | namespace OpenMD {
654  
655      StuntDouble* sd;
656  
657 <    RealType min_val;
658 <    bool min_found = false;  
659 <    StuntDouble* min_sd;
657 >    RealType min_val(0.0);
658 >    int min_found = 0;  
659 >    StuntDouble* min_sd = NULL;
660  
661 <    RealType max_val;
662 <    bool max_found = false;
663 <    StuntDouble* max_sd;
661 >    RealType max_val(0.0);
662 >    int max_found = 0;
663 >    StuntDouble* max_sd = NULL;
664  
665      for (sd = seleManA_.beginSelected(selei); sd != NULL;
666           sd = seleManA_.nextSelected(selei)) {
# Line 682 | Line 714 | namespace OpenMD {
714        if (!max_found) {
715          max_val = value;
716          max_sd = sd;
717 <        max_found = true;
717 >        max_found = 1;
718        } else {
719          if (max_val < value) {
720            max_val = value;
# Line 744 | Line 776 | namespace OpenMD {
776        if (!min_found) {
777          min_val = value;
778          min_sd = sd;
779 <        min_found = true;
779 >        min_found = 1;
780        } else {
781          if (min_val > value) {
782            min_val = value;
# Line 754 | Line 786 | namespace OpenMD {
786      }
787      
788   #ifdef IS_MPI    
789 <    int worldRank = MPI::COMM_WORLD.Get_rank();
790 <    
791 <    bool my_min_found = min_found;
792 <    bool my_max_found = max_found;
789 >    int worldRank;
790 >    MPI_Comm_rank( MPI_COMM_WORLD, &worldRank);
791 >        
792 >    int my_min_found = min_found;
793 >    int my_max_found = max_found;
794  
795      // Even if we didn't find a minimum, did someone else?
796 <    MPI::COMM_WORLD.Allreduce(&my_min_found, &min_found, 1, MPI::BOOL, MPI::LOR);
796 >    MPI_Allreduce(&my_min_found, &min_found, 1, MPI_INT, MPI_LOR,
797 >                  MPI_COMM_WORLD);
798      // Even if we didn't find a maximum, did someone else?
799 <    MPI::COMM_WORLD.Allreduce(&my_max_found, &max_found, 1, MPI::BOOL, MPI::LOR);
799 >    MPI_Allreduce(&my_max_found, &max_found, 1, MPI_INT, MPI_LOR,
800 >                  MPI_COMM_WORLD);
801   #endif
802  
803      if (max_found && min_found) {
# Line 781 | Line 816 | namespace OpenMD {
816        min_vals.rank = worldRank;    
817        
818        // Who had the minimum?
819 <      MPI::COMM_WORLD.Allreduce(&min_vals, &min_vals,
820 <                                1, MPI::REALTYPE_INT, MPI::MINLOC);
819 >      MPI_Allreduce(&min_vals, &min_vals,
820 >                    1, MPI_REALTYPE_INT, MPI_MINLOC, MPI_COMM_WORLD);
821        min_val = min_vals.val;
822        
823        if (my_max_found) {
# Line 793 | Line 828 | namespace OpenMD {
828        max_vals.rank = worldRank;    
829        
830        // Who had the maximum?
831 <      MPI::COMM_WORLD.Allreduce(&max_vals, &max_vals,
832 <                                1, MPI::REALTYPE_INT, MPI::MAXLOC);
831 >      MPI_Allreduce(&max_vals, &max_vals,
832 >                    1, MPI_REALTYPE_INT, MPI_MAXLOC, MPI_COMM_WORLD);
833        max_val = max_vals.val;
834   #endif
835        
# Line 854 | Line 889 | namespace OpenMD {
889            
890            Vector3d min_vel;
891            Vector3d max_vel = max_sd->getVel();
892 <          MPI::Status status;
892 >          MPI_Status status;
893  
894            // point-to-point swap of the velocity vector
895 <          MPI::COMM_WORLD.Sendrecv(max_vel.getArrayPointer(), 3, MPI::REALTYPE,
896 <                                   min_vals.rank, 0,
897 <                                   min_vel.getArrayPointer(), 3, MPI::REALTYPE,
898 <                                   min_vals.rank, 0, status);
895 >          MPI_Sendrecv(max_vel.getArrayPointer(), 3, MPI_REALTYPE,
896 >                       min_vals.rank, 0,
897 >                       min_vel.getArrayPointer(), 3, MPI_REALTYPE,
898 >                       min_vals.rank, 0, MPI_COMM_WORLD, &status);
899            
900            switch(rnemdFluxType_) {
901            case rnemdKE :
# Line 871 | Line 906 | namespace OpenMD {
906                Vector3d max_angMom = max_sd->getJ();
907                
908                // point-to-point swap of the angular momentum vector
909 <              MPI::COMM_WORLD.Sendrecv(max_angMom.getArrayPointer(), 3,
910 <                                       MPI::REALTYPE, min_vals.rank, 1,
911 <                                       min_angMom.getArrayPointer(), 3,
912 <                                       MPI::REALTYPE, min_vals.rank, 1,
913 <                                       status);
909 >              MPI_Sendrecv(max_angMom.getArrayPointer(), 3,
910 >                           MPI_REALTYPE, min_vals.rank, 1,
911 >                           min_angMom.getArrayPointer(), 3,
912 >                           MPI_REALTYPE, min_vals.rank, 1,
913 >                           MPI_COMM_WORLD, &status);
914                
915                max_sd->setJ(min_angMom);
916              }
# Line 900 | Line 935 | namespace OpenMD {
935            
936            Vector3d max_vel;
937            Vector3d min_vel = min_sd->getVel();
938 <          MPI::Status status;
938 >          MPI_Status status;
939            
940            // point-to-point swap of the velocity vector
941 <          MPI::COMM_WORLD.Sendrecv(min_vel.getArrayPointer(), 3, MPI::REALTYPE,
942 <                                   max_vals.rank, 0,
943 <                                   max_vel.getArrayPointer(), 3, MPI::REALTYPE,
944 <                                   max_vals.rank, 0, status);
941 >          MPI_Sendrecv(min_vel.getArrayPointer(), 3, MPI_REALTYPE,
942 >                       max_vals.rank, 0,
943 >                       max_vel.getArrayPointer(), 3, MPI_REALTYPE,
944 >                       max_vals.rank, 0, MPI_COMM_WORLD, &status);
945            
946            switch(rnemdFluxType_) {
947            case rnemdKE :
# Line 917 | Line 952 | namespace OpenMD {
952                Vector3d max_angMom;
953                
954                // point-to-point swap of the angular momentum vector
955 <              MPI::COMM_WORLD.Sendrecv(min_angMom.getArrayPointer(), 3,
956 <                                       MPI::REALTYPE, max_vals.rank, 1,
957 <                                       max_angMom.getArrayPointer(), 3,
958 <                                       MPI::REALTYPE, max_vals.rank, 1,
959 <                                       status);
955 >              MPI_Sendrecv(min_angMom.getArrayPointer(), 3,
956 >                           MPI_REALTYPE, max_vals.rank, 1,
957 >                           max_angMom.getArrayPointer(), 3,
958 >                           MPI_REALTYPE, max_vals.rank, 1,
959 >                           MPI_COMM_WORLD, &status);
960                
961                min_sd->setJ(max_angMom);
962              }
# Line 985 | Line 1020 | namespace OpenMD {
1020      int selej;
1021  
1022      Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
988 –    RealType time = currentSnap_->getTime();    
1023      Mat3x3d hmat = currentSnap_->getHmat();
1024  
1025      StuntDouble* sd;
# Line 1090 | Line 1124 | namespace OpenMD {
1124      Kcw *= 0.5;
1125  
1126   #ifdef IS_MPI
1127 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Phx, 1, MPI::REALTYPE, MPI::SUM);
1128 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Phy, 1, MPI::REALTYPE, MPI::SUM);
1129 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Phz, 1, MPI::REALTYPE, MPI::SUM);
1130 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Pcx, 1, MPI::REALTYPE, MPI::SUM);
1131 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Pcy, 1, MPI::REALTYPE, MPI::SUM);
1132 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Pcz, 1, MPI::REALTYPE, MPI::SUM);
1127 >    MPI_Allreduce(MPI_IN_PLACE, &Phx, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1128 >    MPI_Allreduce(MPI_IN_PLACE, &Phy, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1129 >    MPI_Allreduce(MPI_IN_PLACE, &Phz, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1130 >    MPI_Allreduce(MPI_IN_PLACE, &Pcx, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1131 >    MPI_Allreduce(MPI_IN_PLACE, &Pcy, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1132 >    MPI_Allreduce(MPI_IN_PLACE, &Pcz, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1133  
1134 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Khx, 1, MPI::REALTYPE, MPI::SUM);
1135 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Khy, 1, MPI::REALTYPE, MPI::SUM);
1136 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Khz, 1, MPI::REALTYPE, MPI::SUM);
1137 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Khw, 1, MPI::REALTYPE, MPI::SUM);
1134 >    MPI_Allreduce(MPI_IN_PLACE, &Khx, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1135 >    MPI_Allreduce(MPI_IN_PLACE, &Khy, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1136 >    MPI_Allreduce(MPI_IN_PLACE, &Khz, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1137 >    MPI_Allreduce(MPI_IN_PLACE, &Khw, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1138  
1139 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Kcx, 1, MPI::REALTYPE, MPI::SUM);
1140 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Kcy, 1, MPI::REALTYPE, MPI::SUM);
1141 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Kcz, 1, MPI::REALTYPE, MPI::SUM);
1142 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Kcw, 1, MPI::REALTYPE, MPI::SUM);
1139 >    MPI_Allreduce(MPI_IN_PLACE, &Kcx, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1140 >    MPI_Allreduce(MPI_IN_PLACE, &Kcy, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1141 >    MPI_Allreduce(MPI_IN_PLACE, &Kcz, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1142 >    MPI_Allreduce(MPI_IN_PLACE, &Kcw, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1143   #endif
1144  
1145      //solve coldBin coeff's first
# Line 1433 | Line 1467 | namespace OpenMD {
1467      int selej;
1468  
1469      Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
1436 –    RealType time = currentSnap_->getTime();    
1470      Mat3x3d hmat = currentSnap_->getHmat();
1471  
1472      StuntDouble* sd;
# Line 1582 | Line 1615 | namespace OpenMD {
1615      Kc *= 0.5;
1616      
1617   #ifdef IS_MPI
1618 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Ph[0], 3, MPI::REALTYPE, MPI::SUM);
1619 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Pc[0], 3, MPI::REALTYPE, MPI::SUM);
1620 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Lh[0], 3, MPI::REALTYPE, MPI::SUM);
1621 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Lc[0], 3, MPI::REALTYPE, MPI::SUM);
1622 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Mh, 1, MPI::REALTYPE, MPI::SUM);
1623 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Kh, 1, MPI::REALTYPE, MPI::SUM);
1624 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Mc, 1, MPI::REALTYPE, MPI::SUM);
1625 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Kc, 1, MPI::REALTYPE, MPI::SUM);
1626 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, Ih.getArrayPointer(), 9,
1627 <                              MPI::REALTYPE, MPI::SUM);
1628 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, Ic.getArrayPointer(), 9,
1629 <                              MPI::REALTYPE, MPI::SUM);
1618 >    MPI_Allreduce(MPI_IN_PLACE, &Ph[0], 3, MPI_REALTYPE, MPI_SUM,
1619 >                  MPI_COMM_WORLD);
1620 >    MPI_Allreduce(MPI_IN_PLACE, &Pc[0], 3, MPI_REALTYPE, MPI_SUM,
1621 >                  MPI_COMM_WORLD);
1622 >    MPI_Allreduce(MPI_IN_PLACE, &Lh[0], 3, MPI_REALTYPE, MPI_SUM,
1623 >                  MPI_COMM_WORLD);
1624 >    MPI_Allreduce(MPI_IN_PLACE, &Lc[0], 3, MPI_REALTYPE, MPI_SUM,
1625 >                  MPI_COMM_WORLD);
1626 >    MPI_Allreduce(MPI_IN_PLACE, &Mh, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1627 >    MPI_Allreduce(MPI_IN_PLACE, &Kh, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1628 >    MPI_Allreduce(MPI_IN_PLACE, &Mc, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1629 >    MPI_Allreduce(MPI_IN_PLACE, &Kc, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1630 >    MPI_Allreduce(MPI_IN_PLACE, Ih.getArrayPointer(), 9,
1631 >                  MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1632 >    MPI_Allreduce(MPI_IN_PLACE, Ic.getArrayPointer(), 9,
1633 >                  MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1634   #endif
1635      
1636  
# Line 1736 | Line 1773 | namespace OpenMD {
1773          areaA = evaluatorA_.getSurfaceArea();
1774        else {
1775          
1739 –        cerr << "selection A did not have surface area, recomputing\n";
1776          int isd;
1777          StuntDouble* sd;
1778          vector<StuntDouble*> aSites;
# Line 1774 | Line 1810 | namespace OpenMD {
1810      }
1811  
1812      if (hasSelectionB_) {
1813 <      if (evaluatorB_.hasSurfaceArea())
1813 >      if (evaluatorB_.hasSurfaceArea()) {
1814          areaB = evaluatorB_.getSurfaceArea();
1815 <      else {
1780 <        cerr << "selection B did not have surface area, recomputing\n";
1815 >      } else {
1816  
1817          int isd;
1818          StuntDouble* sd;
# Line 1868 | Line 1903 | namespace OpenMD {
1903    void RNEMD::collectData() {
1904      if (!doRNEMD_) return;
1905      Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
1906 <    
1906 >
1907      // collectData can be called more frequently than the doRNEMD, so use the
1908      // computed area from the last exchange time:
1909      RealType area = getDividingArea();
1910      areaAccumulator_->add(area);
1911      Mat3x3d hmat = currentSnap_->getHmat();
1912 +    Vector3d u = angularMomentumFluxVector_;
1913 +    u.normalize();
1914 +
1915      seleMan_.setSelectionSet(evaluator_.evaluate());
1916  
1917      int selei(0);
1918      StuntDouble* sd;
1919      int binNo;
1920 +    RealType mass;
1921 +    Vector3d vel;
1922 +    Vector3d rPos;
1923 +    RealType KE;
1924 +    Vector3d L;
1925 +    Mat3x3d I;
1926 +    RealType r2;
1927  
1928      vector<RealType> binMass(nBins_, 0.0);
1929 <    vector<RealType> binPx(nBins_, 0.0);
1930 <    vector<RealType> binPy(nBins_, 0.0);
1931 <    vector<RealType> binPz(nBins_, 0.0);
1932 <    vector<RealType> binOmegax(nBins_, 0.0);
1888 <    vector<RealType> binOmegay(nBins_, 0.0);
1889 <    vector<RealType> binOmegaz(nBins_, 0.0);
1929 >    vector<Vector3d> binP(nBins_, V3Zero);
1930 >    vector<RealType> binOmega(nBins_, 0.0);
1931 >    vector<Vector3d> binL(nBins_, V3Zero);
1932 >    vector<Mat3x3d>  binI(nBins_);
1933      vector<RealType> binKE(nBins_, 0.0);
1934      vector<int> binDOF(nBins_, 0);
1935      vector<int> binCount(nBins_, 0);
# Line 1926 | Line 1969 | namespace OpenMD {
1969          binNo = int(rPos.length() / binWidth_);
1970        }
1971  
1972 <      RealType mass = sd->getMass();
1973 <      Vector3d vel = sd->getVel();
1974 <      Vector3d rPos = sd->getPos() - coordinateOrigin_;
1975 <      Vector3d aVel = cross(rPos, vel);
1976 <      
1972 >      mass = sd->getMass();
1973 >      vel = sd->getVel();
1974 >      rPos = sd->getPos() - coordinateOrigin_;
1975 >      KE = 0.5 * mass * vel.lengthSquare();
1976 >      L = mass * cross(rPos, vel);
1977 >      I = outProduct(rPos, rPos) * mass;
1978 >      r2 = rPos.lengthSquare();
1979 >      I(0, 0) += mass * r2;
1980 >      I(1, 1) += mass * r2;
1981 >      I(2, 2) += mass * r2;
1982 >
1983 >      // Project the relative position onto a plane perpendicular to
1984 >      // the angularMomentumFluxVector:
1985 >      // Vector3d rProj = rPos - dot(rPos, u) * u;
1986 >      // Project the velocity onto a plane perpendicular to the
1987 >      // angularMomentumFluxVector:
1988 >      // Vector3d vProj = vel  - dot(vel, u) * u;
1989 >      // Compute angular velocity vector (should be nearly parallel to
1990 >      // angularMomentumFluxVector
1991 >      // Vector3d aVel = cross(rProj, vProj);
1992 >
1993        if (binNo >= 0 && binNo < nBins_)  {
1994          binCount[binNo]++;
1995          binMass[binNo] += mass;
1996 <        binPx[binNo] += mass*vel.x();
1997 <        binPy[binNo] += mass*vel.y();
1998 <        binPz[binNo] += mass*vel.z();
1999 <        binOmegax[binNo] += aVel.x();
1941 <        binOmegay[binNo] += aVel.y();
1942 <        binOmegaz[binNo] += aVel.z();
1943 <        binKE[binNo] += 0.5 * (mass * vel.lengthSquare());
1996 >        binP[binNo] += mass*vel;
1997 >        binKE[binNo] += KE;
1998 >        binI[binNo] += I;
1999 >        binL[binNo] += L;
2000          binDOF[binNo] += 3;
2001          
2002          if (sd->isDirectional()) {
2003            Vector3d angMom = sd->getJ();
2004 <          Mat3x3d I = sd->getI();
2004 >          Mat3x3d Ia = sd->getI();
2005            if (sd->isLinear()) {
2006              int i = sd->linearAxis();
2007              int j = (i + 1) % 3;
2008              int k = (i + 2) % 3;
2009 <            binKE[binNo] += 0.5 * (angMom[j] * angMom[j] / I(j, j) +
2010 <                                   angMom[k] * angMom[k] / I(k, k));
2009 >            binKE[binNo] += 0.5 * (angMom[j] * angMom[j] / Ia(j, j) +
2010 >                                   angMom[k] * angMom[k] / Ia(k, k));
2011              binDOF[binNo] += 2;
2012            } else {
2013 <            binKE[binNo] += 0.5 * (angMom[0] * angMom[0] / I(0, 0) +
2014 <                                   angMom[1] * angMom[1] / I(1, 1) +
2015 <                                   angMom[2] * angMom[2] / I(2, 2));
2013 >            binKE[binNo] += 0.5 * (angMom[0] * angMom[0] / Ia(0, 0) +
2014 >                                   angMom[1] * angMom[1] / Ia(1, 1) +
2015 >                                   angMom[2] * angMom[2] / Ia(2, 2));
2016              binDOF[binNo] += 3;
2017            }
2018          }
2019        }
2020      }
2021 <    
2021 >
2022   #ifdef IS_MPI
2023 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binCount[0],
2024 <                              nBins_, MPI::INT, MPI::SUM);
2025 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binMass[0],
2026 <                              nBins_, MPI::REALTYPE, MPI::SUM);
2027 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binPx[0],
2028 <                              nBins_, MPI::REALTYPE, MPI::SUM);
2029 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binPy[0],
2030 <                              nBins_, MPI::REALTYPE, MPI::SUM);
2031 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binPz[0],
2032 <                              nBins_, MPI::REALTYPE, MPI::SUM);
2033 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binOmegax[0],
2034 <                              nBins_, MPI::REALTYPE, MPI::SUM);
2035 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binOmegay[0],
2036 <                              nBins_, MPI::REALTYPE, MPI::SUM);
2037 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binOmegaz[0],
2038 <                              nBins_, MPI::REALTYPE, MPI::SUM);
2039 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binKE[0],
2040 <                              nBins_, MPI::REALTYPE, MPI::SUM);
2041 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binDOF[0],
2042 <                              nBins_, MPI::INT, MPI::SUM);
2023 >
2024 >    for (int i = 0; i < nBins_; i++) {
2025 >      
2026 >      MPI_Allreduce(MPI_IN_PLACE, &binCount[i],
2027 >                    1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
2028 >      MPI_Allreduce(MPI_IN_PLACE, &binMass[i],
2029 >                    1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
2030 >      MPI_Allreduce(MPI_IN_PLACE, binP[i].getArrayPointer(),
2031 >                    3, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
2032 >      MPI_Allreduce(MPI_IN_PLACE, binL[i].getArrayPointer(),
2033 >                    3, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
2034 >      MPI_Allreduce(MPI_IN_PLACE, binI[i].getArrayPointer(),
2035 >                    9, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
2036 >      MPI_Allreduce(MPI_IN_PLACE, &binKE[i],
2037 >                    1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
2038 >      MPI_Allreduce(MPI_IN_PLACE, &binDOF[i],
2039 >                    1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
2040 >      //MPI_Allreduce(MPI_IN_PLACE, &binOmega[i],
2041 >      //                          1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
2042 >    }
2043 >    
2044   #endif
2045  
2046 <    Vector3d vel;
1990 <    Vector3d aVel;
2046 >    Vector3d omega;
2047      RealType den;
2048      RealType temp;
2049      RealType z;
# Line 2004 | Line 2060 | namespace OpenMD {
2060          den = binMass[i] * 3.0 * PhysicalConstants::densityConvert
2061            / (4.0 * M_PI * (pow(router,3) - pow(rinner,3)));
2062        }
2063 <      vel.x() = binPx[i] / binMass[i];
2064 <      vel.y() = binPy[i] / binMass[i];
2065 <      vel.z() = binPz[i] / binMass[i];
2066 <      aVel.x() = binOmegax[i] / binCount[i];
2067 <      aVel.y() = binOmegay[i] / binCount[i];
2012 <      aVel.z() = binOmegaz[i] / binCount[i];
2063 >      vel = binP[i] / binMass[i];
2064 >
2065 >      omega = binI[i].inverse() * binL[i];
2066 >
2067 >      // omega = binOmega[i] / binCount[i];
2068  
2069        if (binCount[i] > 0) {
2070          // only add values if there are things to add
# Line 2032 | Line 2087 | namespace OpenMD {
2087                dynamic_cast<VectorAccumulator *>(data_[j].accumulator[i])->add(vel);
2088                break;
2089              case ANGULARVELOCITY:  
2090 <              dynamic_cast<VectorAccumulator *>(data_[j].accumulator[i])->add(aVel);
2090 >              dynamic_cast<VectorAccumulator *>(data_[j].accumulator[i])->add(omega);
2091                break;
2092              case DENSITY:
2093                dynamic_cast<Accumulator *>(data_[j].accumulator[i])->add(den);
# Line 2070 | Line 2125 | namespace OpenMD {
2125          painCave.severity = OPENMD_ERROR;
2126          simError();            
2127        }
2128 <    }  
2128 >    }
2129    }
2130    
2131    void RNEMD::writeOutputFile() {
# Line 2079 | Line 2134 | namespace OpenMD {
2134      
2135   #ifdef IS_MPI
2136      // If we're the root node, should we print out the results
2137 <    int worldRank = MPI::COMM_WORLD.Get_rank();
2137 >    int worldRank;
2138 >    MPI_Comm_rank( MPI_COMM_WORLD, &worldRank);
2139 >
2140      if (worldRank == 0) {
2141   #endif
2142        rnemdFile_.open(rnemdFileName_.c_str(), std::ios::out | std::ios::trunc );
# Line 2210 | Line 2267 | namespace OpenMD {
2267        }        
2268  
2269        rnemdFile_ << "#######################################################\n";
2270 <      rnemdFile_ << "# Standard Deviations in those quantities follow:\n";
2270 >      rnemdFile_ << "# 95% confidence intervals in those quantities follow:\n";
2271        rnemdFile_ << "#######################################################\n";
2272  
2273  
# Line 2219 | Line 2276 | namespace OpenMD {
2276          for (unsigned int i = 0; i < outputMask_.size(); ++i) {
2277            if (outputMask_[i]) {
2278              if (data_[i].dataType == "RealType")
2279 <              writeRealStdDev(i,j);
2279 >              writeRealErrorBars(i,j);
2280              else if (data_[i].dataType == "Vector3d")
2281 <              writeVectorStdDev(i,j);
2281 >              writeVectorErrorBars(i,j);
2282              else {
2283                sprintf( painCave.errMsg,
2284                         "RNEMD found an unknown data type for: %s ",
# Line 2292 | Line 2349 | namespace OpenMD {
2349      }
2350    }  
2351  
2352 <  void RNEMD::writeRealStdDev(int index, unsigned int bin) {
2352 >  void RNEMD::writeRealErrorBars(int index, unsigned int bin) {
2353      if (!doRNEMD_) return;
2354      assert(index >=0 && index < ENDINDEX);
2355      assert(int(bin) < nBins_);
# Line 2302 | Line 2359 | namespace OpenMD {
2359      count = data_[index].accumulator[bin]->count();
2360      if (count == 0) return;
2361      
2362 <    dynamic_cast<Accumulator *>(data_[index].accumulator[bin])->getStdDev(s);
2362 >    dynamic_cast<Accumulator *>(data_[index].accumulator[bin])->get95percentConfidenceInterval(s);
2363      
2364      if (! isinf(s) && ! isnan(s)) {
2365        rnemdFile_ << "\t" << s;
# Line 2315 | Line 2372 | namespace OpenMD {
2372      }    
2373    }
2374    
2375 <  void RNEMD::writeVectorStdDev(int index, unsigned int bin) {
2375 >  void RNEMD::writeVectorErrorBars(int index, unsigned int bin) {
2376      if (!doRNEMD_) return;
2377      assert(index >=0 && index < ENDINDEX);
2378      assert(int(bin) < nBins_);
# Line 2325 | Line 2382 | namespace OpenMD {
2382      count = data_[index].accumulator[bin]->count();
2383      if (count == 0) return;
2384  
2385 <    dynamic_cast<VectorAccumulator*>(data_[index].accumulator[bin])->getStdDev(s);
2385 >    dynamic_cast<VectorAccumulator*>(data_[index].accumulator[bin])->get95percentConfidenceInterval(s);
2386      if (isinf(s[0]) || isnan(s[0]) ||
2387          isinf(s[1]) || isnan(s[1]) ||
2388          isinf(s[2]) || isnan(s[2]) ) {      

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines