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 1938 by gezelter, Thu Oct 31 15:32:17 2013 UTC vs.
Revision 2071 by gezelter, Sat Mar 7 21:41:51 2015 UTC

# Line 68 | Line 68 | namespace OpenMD {
68   using namespace std;
69   namespace OpenMD {
70    
71 <  RNEMD::RNEMD(SimInfo* info) : info_(info), evaluator_(info), seleMan_(info),
72 <                                evaluatorA_(info), seleManA_(info),
73 <                                commonA_(info), evaluatorB_(info),
74 <                                seleManB_(info), commonB_(info),
75 <                                hasData_(false), hasDividingArea_(false),
76 <                                usePeriodicBoundaryConditions_(info->getSimParams()->getUsePeriodicBoundaryConditions()) {
71 >  RNEMD::RNEMD(SimInfo* info) : info_(info),
72 >                                evaluator_(info_), seleMan_(info_),
73 >                                evaluatorA_(info_), seleManA_(info_),
74 >                                evaluatorB_(info_), seleManB_(info_),
75 >                                commonA_(info_), commonB_(info_),
76 >                                usePeriodicBoundaryConditions_(info_->getSimParams()->getUsePeriodicBoundaryConditions()),
77 >                                hasDividingArea_(false),
78 >                                hasData_(false) {
79  
80      trialCount_ = 0;
81      failTrialCount_ = 0;
# Line 291 | Line 293 | namespace OpenMD {
293        kineticFlux_ = 0.0;
294      }
295      if (hasMomentumFluxVector) {
296 <      momentumFluxVector_ = rnemdParams->getMomentumFluxVector();
296 >      std::vector<RealType> mf = rnemdParams->getMomentumFluxVector();
297 >      if (mf.size() != 3) {
298 >          sprintf(painCave.errMsg,
299 >                  "RNEMD: Incorrect number of parameters specified for momentumFluxVector.\n"
300 >                  "\tthere should be 3 parameters, but %lu were specified.\n",
301 >                  mf.size());
302 >          painCave.isFatal = 1;
303 >          simError();      
304 >      }
305 >      momentumFluxVector_.x() = mf[0];
306 >      momentumFluxVector_.y() = mf[1];
307 >      momentumFluxVector_.z() = mf[2];
308      } else {
309        momentumFluxVector_ = V3Zero;
310        if (hasMomentumFlux) {
# Line 316 | Line 329 | namespace OpenMD {
329            break;
330          }
331        }
332 <      if (hasAngularMomentumFluxVector) {
333 <        angularMomentumFluxVector_ = rnemdParams->getAngularMomentumFluxVector();
334 <      } else {
335 <        angularMomentumFluxVector_ = V3Zero;
336 <        if (hasAngularMomentumFlux) {
337 <          RealType angularMomentumFlux = rnemdParams->getAngularMomentumFlux();
338 <          switch (rnemdFluxType_) {
339 <          case rnemdLx:
340 <            angularMomentumFluxVector_.x() = angularMomentumFlux;
341 <            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 <        }        
332 >    }
333 >    if (hasAngularMomentumFluxVector) {
334 >      std::vector<RealType> amf = rnemdParams->getAngularMomentumFluxVector();
335 >      if (amf.size() != 3) {
336 >        sprintf(painCave.errMsg,
337 >                "RNEMD: Incorrect number of parameters specified for angularMomentumFluxVector.\n"
338 >                "\tthere should be 3 parameters, but %lu were specified.\n",
339 >                amf.size());
340 >        painCave.isFatal = 1;
341 >        simError();      
342        }
343 <
344 <      if (hasCoordinateOrigin) {
345 <        coordinateOrigin_ = rnemdParams->getCoordinateOrigin();
346 <      } else {
347 <        coordinateOrigin_ = V3Zero;
348 <      }
349 <
356 <      // do some sanity checking
357 <
358 <      int selectionCount = seleMan_.getSelectionCount();
359 <
360 <      int nIntegrable = info->getNGlobalIntegrableObjects();
361 <
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);
343 >      angularMomentumFluxVector_.x() = amf[0];
344 >      angularMomentumFluxVector_.y() = amf[1];
345 >      angularMomentumFluxVector_.z() = amf[2];
346 >    } else {
347 >      angularMomentumFluxVector_ = V3Zero;
348 >      if (hasAngularMomentumFlux) {
349 >        RealType angularMomentumFlux = rnemdParams->getAngularMomentumFlux();
350          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;
351          case rnemdLx:
352 +          angularMomentumFluxVector_.x() = angularMomentumFlux;
353 +          break;
354          case rnemdLy:
355 +          angularMomentumFluxVector_.y() = angularMomentumFlux;
356 +          break;
357          case rnemdLz:
358 <        case rnemdLvector:
470 <          outputMask_.set(ANGULARVELOCITY);
358 >          angularMomentumFluxVector_.z() = angularMomentumFlux;
359            break;
360          case rnemdKeLx:
361 +          angularMomentumFluxVector_.x() = angularMomentumFlux;
362 +          break;
363          case rnemdKeLy:
364 +          angularMomentumFluxVector_.y() = angularMomentumFlux;
365 +          break;
366          case rnemdKeLz:
367 <        case rnemdKeLvector:
476 <          outputMask_.set(TEMPERATURE);
477 <          outputMask_.set(ANGULARVELOCITY);
367 >          angularMomentumFluxVector_.z() = angularMomentumFlux;
368            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;
369          default:
370            break;
371          }
372 +      }        
373 +    }
374 +
375 +    if (hasCoordinateOrigin) {
376 +      std::vector<RealType> co = rnemdParams->getCoordinateOrigin();
377 +      if (co.size() != 3) {
378 +        sprintf(painCave.errMsg,
379 +                "RNEMD: Incorrect number of parameters specified for coordinateOrigin.\n"
380 +                "\tthere should be 3 parameters, but %lu were specified.\n",
381 +                co.size());
382 +        painCave.isFatal = 1;
383 +        simError();      
384        }
385 <      
386 <      if (hasOutputFileName) {
387 <        rnemdFileName_ = rnemdParams->getOutputFileName();
388 <      } else {
389 <        rnemdFileName_ = getPrefix(info->getFinalConfigFileName()) + ".rnemd";
390 <      }          
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;
385 >      coordinateOrigin_.x() = co[0];
386 >      coordinateOrigin_.y() = co[1];
387 >      coordinateOrigin_.z() = co[2];
388 >    } else {
389 >      coordinateOrigin_ = V3Zero;
390 >    }
391      
392 <      if (hasSelectionA_) {
393 <        selectionA_ = rnemdParams->getSelectionA();
394 <      } else {
395 <        if (usePeriodicBoundaryConditions_) {    
396 <          Mat3x3d hmat = currentSnap_->getHmat();
392 >    // do some sanity checking
393 >    
394 >    int selectionCount = seleMan_.getSelectionCount();    
395 >    int nIntegrable = info->getNGlobalIntegrableObjects();
396 >    if (selectionCount > nIntegrable) {
397 >      sprintf(painCave.errMsg,
398 >              "RNEMD: The current objectSelection,\n"
399 >              "\t\t%s\n"
400 >              "\thas resulted in %d selected objects.  However,\n"
401 >              "\tthe total number of integrable objects in the system\n"
402 >              "\tis only %d.  This is almost certainly not what you want\n"
403 >              "\tto do.  A likely cause of this is forgetting the _RB_0\n"
404 >              "\tselector in the selection script!\n",
405 >              rnemdObjectSelection_.c_str(),
406 >              selectionCount, nIntegrable);
407 >      painCave.isFatal = 0;
408 >      painCave.severity = OPENMD_WARNING;
409 >      simError();
410 >    }
411 >    
412 >    areaAccumulator_ = new Accumulator();
413 >    
414 >    nBins_ = rnemdParams->getOutputBins();
415 >    binWidth_ = rnemdParams->getOutputBinWidth();
416 >    
417 >    data_.resize(RNEMD::ENDINDEX);
418 >    OutputData z;
419 >    z.units =  "Angstroms";
420 >    z.title =  "Z";
421 >    z.dataType = "RealType";
422 >    z.accumulator.reserve(nBins_);
423 >    for (int i = 0; i < nBins_; i++)
424 >      z.accumulator.push_back( new Accumulator() );
425 >    data_[Z] = z;
426 >    outputMap_["Z"] =  Z;
427 >    
428 >    OutputData r;
429 >    r.units =  "Angstroms";
430 >    r.title =  "R";
431 >    r.dataType = "RealType";
432 >    r.accumulator.reserve(nBins_);
433 >    for (int i = 0; i < nBins_; i++)
434 >      r.accumulator.push_back( new Accumulator() );
435 >    data_[R] = r;
436 >    outputMap_["R"] =  R;
437 >    
438 >    OutputData temperature;
439 >    temperature.units =  "K";
440 >    temperature.title =  "Temperature";
441 >    temperature.dataType = "RealType";
442 >    temperature.accumulator.reserve(nBins_);
443 >    for (int i = 0; i < nBins_; i++)
444 >      temperature.accumulator.push_back( new Accumulator() );
445 >    data_[TEMPERATURE] = temperature;
446 >    outputMap_["TEMPERATURE"] =  TEMPERATURE;
447 >    
448 >    OutputData velocity;
449 >    velocity.units = "angstroms/fs";
450 >    velocity.title =  "Velocity";  
451 >    velocity.dataType = "Vector3d";
452 >    velocity.accumulator.reserve(nBins_);
453 >    for (int i = 0; i < nBins_; i++)
454 >      velocity.accumulator.push_back( new VectorAccumulator() );
455 >    data_[VELOCITY] = velocity;
456 >    outputMap_["VELOCITY"] = VELOCITY;
457 >    
458 >    OutputData angularVelocity;
459 >    angularVelocity.units = "angstroms^2/fs";
460 >    angularVelocity.title =  "AngularVelocity";  
461 >    angularVelocity.dataType = "Vector3d";
462 >    angularVelocity.accumulator.reserve(nBins_);
463 >    for (int i = 0; i < nBins_; i++)
464 >      angularVelocity.accumulator.push_back( new VectorAccumulator() );
465 >    data_[ANGULARVELOCITY] = angularVelocity;
466 >    outputMap_["ANGULARVELOCITY"] = ANGULARVELOCITY;
467 >    
468 >    OutputData density;
469 >    density.units =  "g cm^-3";
470 >    density.title =  "Density";
471 >    density.dataType = "RealType";
472 >    density.accumulator.reserve(nBins_);
473 >    for (int i = 0; i < nBins_; i++)
474 >      density.accumulator.push_back( new Accumulator() );
475 >    data_[DENSITY] = density;
476 >    outputMap_["DENSITY"] =  DENSITY;
477 >    
478 >    if (hasOutputFields) {
479 >      parseOutputFileFormat(rnemdParams->getOutputFields());
480 >    } else {
481 >      if (usePeriodicBoundaryConditions_)
482 >        outputMask_.set(Z);
483 >      else
484 >        outputMask_.set(R);
485 >      switch (rnemdFluxType_) {
486 >      case rnemdKE:
487 >      case rnemdRotKE:
488 >      case rnemdFullKE:
489 >        outputMask_.set(TEMPERATURE);
490 >        break;
491 >      case rnemdPx:
492 >      case rnemdPy:
493 >        outputMask_.set(VELOCITY);
494 >        break;
495 >      case rnemdPz:        
496 >      case rnemdPvector:
497 >        outputMask_.set(VELOCITY);
498 >        outputMask_.set(DENSITY);
499 >        break;
500 >      case rnemdLx:
501 >      case rnemdLy:
502 >      case rnemdLz:
503 >      case rnemdLvector:
504 >        outputMask_.set(ANGULARVELOCITY);
505 >        break;
506 >      case rnemdKeLx:
507 >      case rnemdKeLy:
508 >      case rnemdKeLz:
509 >      case rnemdKeLvector:
510 >        outputMask_.set(TEMPERATURE);
511 >        outputMask_.set(ANGULARVELOCITY);
512 >        break;
513 >      case rnemdKePx:
514 >      case rnemdKePy:
515 >        outputMask_.set(TEMPERATURE);
516 >        outputMask_.set(VELOCITY);
517 >        break;
518 >      case rnemdKePvector:
519 >        outputMask_.set(TEMPERATURE);
520 >        outputMask_.set(VELOCITY);
521 >        outputMask_.set(DENSITY);        
522 >        break;
523 >      default:
524 >        break;
525 >      }
526 >    }
527 >    
528 >    if (hasOutputFileName) {
529 >      rnemdFileName_ = rnemdParams->getOutputFileName();
530 >    } else {
531 >      rnemdFileName_ = getPrefix(info->getFinalConfigFileName()) + ".rnemd";
532 >    }          
533 >    
534 >    exchangeTime_ = rnemdParams->getExchangeTime();
535 >    
536 >    Snapshot* currentSnap_ = info->getSnapshotManager()->getCurrentSnapshot();
537 >    // total exchange sums are zeroed out at the beginning:
538 >    
539 >    kineticExchange_ = 0.0;
540 >    momentumExchange_ = V3Zero;
541 >    angularMomentumExchange_ = V3Zero;
542 >    
543 >    std::ostringstream selectionAstream;
544 >    std::ostringstream selectionBstream;
545 >    
546 >    if (hasSelectionA_) {
547 >      selectionA_ = rnemdParams->getSelectionA();
548 >    } else {
549 >      if (usePeriodicBoundaryConditions_) {    
550 >        Mat3x3d hmat = currentSnap_->getHmat();
551          
552 <          if (hasSlabWidth)
553 <            slabWidth_ = rnemdParams->getSlabWidth();
554 <          else
555 <            slabWidth_ = hmat(2,2) / 10.0;
552 >        if (hasSlabWidth)
553 >          slabWidth_ = rnemdParams->getSlabWidth();
554 >        else
555 >          slabWidth_ = hmat(2,2) / 10.0;
556          
557 <          if (hasSlabACenter)
558 <            slabACenter_ = rnemdParams->getSlabACenter();
559 <          else
560 <            slabACenter_ = 0.0;
557 >        if (hasSlabACenter)
558 >          slabACenter_ = rnemdParams->getSlabACenter();
559 >        else
560 >          slabACenter_ = 0.0;
561          
562 <          selectionAstream << "select wrappedz > "
563 <                           << slabACenter_ - 0.5*slabWidth_
564 <                           <<  " && wrappedz < "
565 <                           << slabACenter_ + 0.5*slabWidth_;
566 <          selectionA_ = selectionAstream.str();
567 <        } else {
568 <          if (hasSphereARadius)
569 <            sphereARadius_ = rnemdParams->getSphereARadius();
570 <          else {
571 <            // use an initial guess to the size of the inner slab to be 1/10 the
572 <            // radius of an approximately spherical hull:
573 <            Thermo thermo(info);
574 <            RealType hVol = thermo.getHullVolume();
575 <            sphereARadius_ = 0.1 * pow((3.0 * hVol / (4.0 * M_PI)), 1.0/3.0);
542 <          }
543 <          selectionAstream << "select r < " << sphereARadius_;
544 <          selectionA_ = selectionAstream.str();
562 >        selectionAstream << "select wrappedz > "
563 >                         << slabACenter_ - 0.5*slabWidth_
564 >                         <<  " && wrappedz < "
565 >                         << slabACenter_ + 0.5*slabWidth_;
566 >        selectionA_ = selectionAstream.str();
567 >      } else {
568 >        if (hasSphereARadius)
569 >          sphereARadius_ = rnemdParams->getSphereARadius();
570 >        else {
571 >          // use an initial guess to the size of the inner slab to be 1/10 the
572 >          // radius of an approximately spherical hull:
573 >          Thermo thermo(info);
574 >          RealType hVol = thermo.getHullVolume();
575 >          sphereARadius_ = 0.1 * pow((3.0 * hVol / (4.0 * M_PI)), 1.0/3.0);
576          }
577 +        selectionAstream << "select r < " << sphereARadius_;
578 +        selectionA_ = selectionAstream.str();
579        }
580 +    }
581      
582 <      if (hasSelectionB_) {
583 <        selectionB_ = rnemdParams->getSelectionB();
584 <
585 <      } else {
586 <        if (usePeriodicBoundaryConditions_) {    
587 <          Mat3x3d hmat = currentSnap_->getHmat();
582 >    if (hasSelectionB_) {
583 >      selectionB_ = rnemdParams->getSelectionB();
584 >      
585 >    } else {
586 >      if (usePeriodicBoundaryConditions_) {    
587 >        Mat3x3d hmat = currentSnap_->getHmat();
588          
589 <          if (hasSlabWidth)
590 <            slabWidth_ = rnemdParams->getSlabWidth();
591 <          else
592 <            slabWidth_ = hmat(2,2) / 10.0;
589 >        if (hasSlabWidth)
590 >          slabWidth_ = rnemdParams->getSlabWidth();
591 >        else
592 >          slabWidth_ = hmat(2,2) / 10.0;
593          
594 <          if (hasSlabBCenter)
595 <            slabBCenter_ = rnemdParams->getSlabBCenter();
596 <          else
597 <            slabBCenter_ = hmat(2,2) / 2.0;
594 >        if (hasSlabBCenter)
595 >          slabBCenter_ = rnemdParams->getSlabBCenter();
596 >        else
597 >          slabBCenter_ = hmat(2,2) / 2.0;
598          
599 <          selectionBstream << "select wrappedz > "
600 <                           << slabBCenter_ - 0.5*slabWidth_
601 <                           <<  " && wrappedz < "
602 <                           << slabBCenter_ + 0.5*slabWidth_;
599 >        selectionBstream << "select wrappedz > "
600 >                         << slabBCenter_ - 0.5*slabWidth_
601 >                         <<  " && wrappedz < "
602 >                         << slabBCenter_ + 0.5*slabWidth_;
603 >        selectionB_ = selectionBstream.str();
604 >      } else {
605 >        if (hasSphereBRadius_) {
606 >          sphereBRadius_ = rnemdParams->getSphereBRadius();
607 >          selectionBstream << "select r > " << sphereBRadius_;
608            selectionB_ = selectionBstream.str();
609          } else {
610 <          if (hasSphereBRadius_) {
611 <            sphereBRadius_ = rnemdParams->getSphereBRadius();
612 <            selectionBstream << "select r > " << sphereBRadius_;
574 <            selectionB_ = selectionBstream.str();
575 <          } else {
576 <            selectionB_ = "select hull";
577 <            BisHull_ = true;
578 <            hasSelectionB_ = true;
579 <          }
610 >          selectionB_ = "select hull";
611 >          BisHull_ = true;
612 >          hasSelectionB_ = true;
613          }
614        }
615      }
616 <
616 >  
617 >  
618      // object evaluator:
619      evaluator_.loadScriptString(rnemdObjectSelection_);
620      seleMan_.setSelectionSet(evaluator_.evaluate());
# Line 622 | Line 656 | namespace OpenMD {
656  
657      StuntDouble* sd;
658  
659 <    RealType min_val;
660 <    bool min_found = false;  
661 <    StuntDouble* min_sd;
659 >    RealType min_val(0.0);
660 >    int min_found = 0;  
661 >    StuntDouble* min_sd = NULL;
662  
663 <    RealType max_val;
664 <    bool max_found = false;
665 <    StuntDouble* max_sd;
663 >    RealType max_val(0.0);
664 >    int max_found = 0;
665 >    StuntDouble* max_sd = NULL;
666  
667      for (sd = seleManA_.beginSelected(selei); sd != NULL;
668           sd = seleManA_.nextSelected(selei)) {
# Line 642 | Line 676 | namespace OpenMD {
676        
677        RealType mass = sd->getMass();
678        Vector3d vel = sd->getVel();
679 <      RealType value;
679 >      RealType value(0.0);
680        
681        switch(rnemdFluxType_) {
682        case rnemdKE :
# Line 682 | Line 716 | namespace OpenMD {
716        if (!max_found) {
717          max_val = value;
718          max_sd = sd;
719 <        max_found = true;
719 >        max_found = 1;
720        } else {
721          if (max_val < value) {
722            max_val = value;
# Line 703 | Line 737 | namespace OpenMD {
737        
738        RealType mass = sd->getMass();
739        Vector3d vel = sd->getVel();
740 <      RealType value;
740 >      RealType value(0.0);
741        
742        switch(rnemdFluxType_) {
743        case rnemdKE :
# Line 744 | Line 778 | namespace OpenMD {
778        if (!min_found) {
779          min_val = value;
780          min_sd = sd;
781 <        min_found = true;
781 >        min_found = 1;
782        } else {
783          if (min_val > value) {
784            min_val = value;
# Line 754 | Line 788 | namespace OpenMD {
788      }
789      
790   #ifdef IS_MPI    
791 <    int worldRank = MPI::COMM_WORLD.Get_rank();
792 <    
793 <    bool my_min_found = min_found;
794 <    bool my_max_found = max_found;
791 >    int worldRank;
792 >    MPI_Comm_rank( MPI_COMM_WORLD, &worldRank);
793 >        
794 >    int my_min_found = min_found;
795 >    int my_max_found = max_found;
796  
797      // Even if we didn't find a minimum, did someone else?
798 <    MPI::COMM_WORLD.Allreduce(&my_min_found, &min_found, 1, MPI::BOOL, MPI::LOR);
798 >    MPI_Allreduce(&my_min_found, &min_found, 1, MPI_INT, MPI_LOR,
799 >                  MPI_COMM_WORLD);
800      // Even if we didn't find a maximum, did someone else?
801 <    MPI::COMM_WORLD.Allreduce(&my_max_found, &max_found, 1, MPI::BOOL, MPI::LOR);
801 >    MPI_Allreduce(&my_max_found, &max_found, 1, MPI_INT, MPI_LOR,
802 >                  MPI_COMM_WORLD);
803   #endif
804  
805      if (max_found && min_found) {
# Line 781 | Line 818 | namespace OpenMD {
818        min_vals.rank = worldRank;    
819        
820        // Who had the minimum?
821 <      MPI::COMM_WORLD.Allreduce(&min_vals, &min_vals,
822 <                                1, MPI::REALTYPE_INT, MPI::MINLOC);
821 >      MPI_Allreduce(&min_vals, &min_vals,
822 >                    1, MPI_REALTYPE_INT, MPI_MINLOC, MPI_COMM_WORLD);
823        min_val = min_vals.val;
824        
825        if (my_max_found) {
# Line 793 | Line 830 | namespace OpenMD {
830        max_vals.rank = worldRank;    
831        
832        // Who had the maximum?
833 <      MPI::COMM_WORLD.Allreduce(&max_vals, &max_vals,
834 <                                1, MPI::REALTYPE_INT, MPI::MAXLOC);
833 >      MPI_Allreduce(&max_vals, &max_vals,
834 >                    1, MPI_REALTYPE_INT, MPI_MAXLOC, MPI_COMM_WORLD);
835        max_val = max_vals.val;
836   #endif
837        
# Line 854 | Line 891 | namespace OpenMD {
891            
892            Vector3d min_vel;
893            Vector3d max_vel = max_sd->getVel();
894 <          MPI::Status status;
894 >          MPI_Status status;
895  
896            // point-to-point swap of the velocity vector
897 <          MPI::COMM_WORLD.Sendrecv(max_vel.getArrayPointer(), 3, MPI::REALTYPE,
898 <                                   min_vals.rank, 0,
899 <                                   min_vel.getArrayPointer(), 3, MPI::REALTYPE,
900 <                                   min_vals.rank, 0, status);
897 >          MPI_Sendrecv(max_vel.getArrayPointer(), 3, MPI_REALTYPE,
898 >                       min_vals.rank, 0,
899 >                       min_vel.getArrayPointer(), 3, MPI_REALTYPE,
900 >                       min_vals.rank, 0, MPI_COMM_WORLD, &status);
901            
902            switch(rnemdFluxType_) {
903            case rnemdKE :
# Line 871 | Line 908 | namespace OpenMD {
908                Vector3d max_angMom = max_sd->getJ();
909                
910                // point-to-point swap of the angular momentum vector
911 <              MPI::COMM_WORLD.Sendrecv(max_angMom.getArrayPointer(), 3,
912 <                                       MPI::REALTYPE, min_vals.rank, 1,
913 <                                       min_angMom.getArrayPointer(), 3,
914 <                                       MPI::REALTYPE, min_vals.rank, 1,
915 <                                       status);
911 >              MPI_Sendrecv(max_angMom.getArrayPointer(), 3,
912 >                           MPI_REALTYPE, min_vals.rank, 1,
913 >                           min_angMom.getArrayPointer(), 3,
914 >                           MPI_REALTYPE, min_vals.rank, 1,
915 >                           MPI_COMM_WORLD, &status);
916                
917                max_sd->setJ(min_angMom);
918              }
# Line 900 | Line 937 | namespace OpenMD {
937            
938            Vector3d max_vel;
939            Vector3d min_vel = min_sd->getVel();
940 <          MPI::Status status;
940 >          MPI_Status status;
941            
942            // point-to-point swap of the velocity vector
943 <          MPI::COMM_WORLD.Sendrecv(min_vel.getArrayPointer(), 3, MPI::REALTYPE,
944 <                                   max_vals.rank, 0,
945 <                                   max_vel.getArrayPointer(), 3, MPI::REALTYPE,
946 <                                   max_vals.rank, 0, status);
943 >          MPI_Sendrecv(min_vel.getArrayPointer(), 3, MPI_REALTYPE,
944 >                       max_vals.rank, 0,
945 >                       max_vel.getArrayPointer(), 3, MPI_REALTYPE,
946 >                       max_vals.rank, 0, MPI_COMM_WORLD, &status);
947            
948            switch(rnemdFluxType_) {
949            case rnemdKE :
# Line 917 | Line 954 | namespace OpenMD {
954                Vector3d max_angMom;
955                
956                // point-to-point swap of the angular momentum vector
957 <              MPI::COMM_WORLD.Sendrecv(min_angMom.getArrayPointer(), 3,
958 <                                       MPI::REALTYPE, max_vals.rank, 1,
959 <                                       max_angMom.getArrayPointer(), 3,
960 <                                       MPI::REALTYPE, max_vals.rank, 1,
961 <                                       status);
957 >              MPI_Sendrecv(min_angMom.getArrayPointer(), 3,
958 >                           MPI_REALTYPE, max_vals.rank, 1,
959 >                           max_angMom.getArrayPointer(), 3,
960 >                           MPI_REALTYPE, max_vals.rank, 1,
961 >                           MPI_COMM_WORLD, &status);
962                
963                min_sd->setJ(max_angMom);
964              }
# Line 985 | Line 1022 | namespace OpenMD {
1022      int selej;
1023  
1024      Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
988 –    RealType time = currentSnap_->getTime();    
1025      Mat3x3d hmat = currentSnap_->getHmat();
1026  
1027      StuntDouble* sd;
# Line 1090 | Line 1126 | namespace OpenMD {
1126      Kcw *= 0.5;
1127  
1128   #ifdef IS_MPI
1129 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Phx, 1, MPI::REALTYPE, MPI::SUM);
1130 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Phy, 1, MPI::REALTYPE, MPI::SUM);
1131 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Phz, 1, MPI::REALTYPE, MPI::SUM);
1132 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Pcx, 1, MPI::REALTYPE, MPI::SUM);
1133 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Pcy, 1, MPI::REALTYPE, MPI::SUM);
1134 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Pcz, 1, MPI::REALTYPE, MPI::SUM);
1129 >    MPI_Allreduce(MPI_IN_PLACE, &Phx, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1130 >    MPI_Allreduce(MPI_IN_PLACE, &Phy, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1131 >    MPI_Allreduce(MPI_IN_PLACE, &Phz, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1132 >    MPI_Allreduce(MPI_IN_PLACE, &Pcx, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1133 >    MPI_Allreduce(MPI_IN_PLACE, &Pcy, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1134 >    MPI_Allreduce(MPI_IN_PLACE, &Pcz, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1135  
1136 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Khx, 1, MPI::REALTYPE, MPI::SUM);
1137 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Khy, 1, MPI::REALTYPE, MPI::SUM);
1138 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Khz, 1, MPI::REALTYPE, MPI::SUM);
1139 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Khw, 1, MPI::REALTYPE, MPI::SUM);
1136 >    MPI_Allreduce(MPI_IN_PLACE, &Khx, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1137 >    MPI_Allreduce(MPI_IN_PLACE, &Khy, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1138 >    MPI_Allreduce(MPI_IN_PLACE, &Khz, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1139 >    MPI_Allreduce(MPI_IN_PLACE, &Khw, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1140  
1141 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Kcx, 1, MPI::REALTYPE, MPI::SUM);
1142 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Kcy, 1, MPI::REALTYPE, MPI::SUM);
1143 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Kcz, 1, MPI::REALTYPE, MPI::SUM);
1144 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Kcw, 1, MPI::REALTYPE, MPI::SUM);
1141 >    MPI_Allreduce(MPI_IN_PLACE, &Kcx, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1142 >    MPI_Allreduce(MPI_IN_PLACE, &Kcy, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1143 >    MPI_Allreduce(MPI_IN_PLACE, &Kcz, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1144 >    MPI_Allreduce(MPI_IN_PLACE, &Kcw, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1145   #endif
1146  
1147      //solve coldBin coeff's first
1148      RealType px = Pcx / Phx;
1149      RealType py = Pcy / Phy;
1150      RealType pz = Pcz / Phz;
1151 <    RealType c, x, y, z;
1151 >    RealType c(0.0), x(0.0), y(0.0), z(0.0);
1152      bool successfulScale = false;
1153      if ((rnemdFluxType_ == rnemdFullKE) ||
1154          (rnemdFluxType_ == rnemdRotKE)) {
# Line 1182 | Line 1218 | namespace OpenMD {
1218          }
1219        }
1220      } else {
1221 <      RealType a000, a110, c0, a001, a111, b01, b11, c1;
1221 >      RealType a000(0.0), a110(0.0), c0(0.0);
1222 >      RealType a001(0.0), a111(0.0), b01(0.0), b11(0.0), c1(0.0);
1223        switch(rnemdFluxType_) {
1224        case rnemdKE :
1225          /* used hotBin coeff's & only scale x & y dimensions
# Line 1285 | Line 1322 | namespace OpenMD {
1322        vector<pair<RealType,RealType> > rps;
1323        for (ri = realRoots.begin(); ri !=realRoots.end(); ++ri) {
1324          r2 = *ri;
1325 <        //check if FindRealRoots() give the right answer
1325 >        // Check to see if FindRealRoots() gave the right answer:
1326          if ( fabs(u0 + r2 * (u1 + r2 * (u2 + r2 * (u3 + r2 * u4)))) > 1e-6 ) {
1327            sprintf(painCave.errMsg,
1328                    "RNEMD Warning: polynomial solve seems to have an error!");
# Line 1293 | Line 1330 | namespace OpenMD {
1330            simError();
1331            failRootCount_++;
1332          }
1333 <        //might not be useful w/o rescaling coefficients
1333 >        // Might not be useful w/o rescaling coefficients
1334          alpha0 = -c0 - a110 * r2 * r2;
1335          if (alpha0 >= 0.0) {
1336            r1 = sqrt(alpha0 / a000);
# Line 1308 | Line 1345 | namespace OpenMD {
1345            }
1346          }
1347        }
1348 <      // Consider combining together the solving pair part w/ the searching
1349 <      // best solution part so that we don't need the pairs vector
1348 >      // Consider combining together the part for solving for the pair
1349 >      // w/ the searching for the best solution part so that we don't
1350 >      // need the pairs vector:
1351        if (!rps.empty()) {
1352          RealType smallestDiff = HONKING_LARGE_VALUE;
1353 <        RealType diff;
1354 <        pair<RealType,RealType> bestPair = make_pair(1.0, 1.0);
1355 <        vector<pair<RealType,RealType> >::iterator rpi;
1353 >        RealType diff(0.0);
1354 >        std::pair<RealType,RealType> bestPair = std::make_pair(1.0, 1.0);
1355 >        std::vector<std::pair<RealType,RealType> >::iterator rpi;
1356          for (rpi = rps.begin(); rpi != rps.end(); ++rpi) {
1357            r1 = (*rpi).first;
1358            r2 = (*rpi).second;
# Line 1433 | Line 1471 | namespace OpenMD {
1471      int selej;
1472  
1473      Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
1436 –    RealType time = currentSnap_->getTime();    
1474      Mat3x3d hmat = currentSnap_->getHmat();
1475  
1476      StuntDouble* sd;
# Line 1582 | Line 1619 | namespace OpenMD {
1619      Kc *= 0.5;
1620      
1621   #ifdef IS_MPI
1622 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Ph[0], 3, MPI::REALTYPE, MPI::SUM);
1623 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Pc[0], 3, MPI::REALTYPE, MPI::SUM);
1624 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Lh[0], 3, MPI::REALTYPE, MPI::SUM);
1625 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Lc[0], 3, MPI::REALTYPE, MPI::SUM);
1626 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Mh, 1, MPI::REALTYPE, MPI::SUM);
1627 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Kh, 1, MPI::REALTYPE, MPI::SUM);
1628 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Mc, 1, MPI::REALTYPE, MPI::SUM);
1629 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Kc, 1, MPI::REALTYPE, MPI::SUM);
1630 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, Ih.getArrayPointer(), 9,
1631 <                              MPI::REALTYPE, MPI::SUM);
1632 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, Ic.getArrayPointer(), 9,
1633 <                              MPI::REALTYPE, MPI::SUM);
1622 >    MPI_Allreduce(MPI_IN_PLACE, &Ph[0], 3, MPI_REALTYPE, MPI_SUM,
1623 >                  MPI_COMM_WORLD);
1624 >    MPI_Allreduce(MPI_IN_PLACE, &Pc[0], 3, MPI_REALTYPE, MPI_SUM,
1625 >                  MPI_COMM_WORLD);
1626 >    MPI_Allreduce(MPI_IN_PLACE, &Lh[0], 3, MPI_REALTYPE, MPI_SUM,
1627 >                  MPI_COMM_WORLD);
1628 >    MPI_Allreduce(MPI_IN_PLACE, &Lc[0], 3, MPI_REALTYPE, MPI_SUM,
1629 >                  MPI_COMM_WORLD);
1630 >    MPI_Allreduce(MPI_IN_PLACE, &Mh, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1631 >    MPI_Allreduce(MPI_IN_PLACE, &Kh, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1632 >    MPI_Allreduce(MPI_IN_PLACE, &Mc, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1633 >    MPI_Allreduce(MPI_IN_PLACE, &Kc, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1634 >    MPI_Allreduce(MPI_IN_PLACE, Ih.getArrayPointer(), 9,
1635 >                  MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1636 >    MPI_Allreduce(MPI_IN_PLACE, Ic.getArrayPointer(), 9,
1637 >                  MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
1638   #endif
1639      
1640  
# Line 1736 | Line 1777 | namespace OpenMD {
1777          areaA = evaluatorA_.getSurfaceArea();
1778        else {
1779          
1739 –        cerr << "selection A did not have surface area, recomputing\n";
1780          int isd;
1781          StuntDouble* sd;
1782          vector<StuntDouble*> aSites;
# Line 1774 | Line 1814 | namespace OpenMD {
1814      }
1815  
1816      if (hasSelectionB_) {
1817 <      if (evaluatorB_.hasSurfaceArea())
1817 >      if (evaluatorB_.hasSurfaceArea()) {
1818          areaB = evaluatorB_.getSurfaceArea();
1819 <      else {
1780 <        cerr << "selection B did not have surface area, recomputing\n";
1819 >      } else {
1820  
1821          int isd;
1822          StuntDouble* sd;
# Line 1868 | Line 1907 | namespace OpenMD {
1907    void RNEMD::collectData() {
1908      if (!doRNEMD_) return;
1909      Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
1910 <    
1910 >
1911      // collectData can be called more frequently than the doRNEMD, so use the
1912      // computed area from the last exchange time:
1913      RealType area = getDividingArea();
1914      areaAccumulator_->add(area);
1915      Mat3x3d hmat = currentSnap_->getHmat();
1916 +    Vector3d u = angularMomentumFluxVector_;
1917 +    u.normalize();
1918 +
1919      seleMan_.setSelectionSet(evaluator_.evaluate());
1920  
1921      int selei(0);
1922      StuntDouble* sd;
1923      int binNo;
1924 +    RealType mass;
1925 +    Vector3d vel;
1926 +    Vector3d rPos;
1927 +    RealType KE;
1928 +    Vector3d L;
1929 +    Mat3x3d I;
1930 +    RealType r2;
1931  
1932      vector<RealType> binMass(nBins_, 0.0);
1933 <    vector<RealType> binPx(nBins_, 0.0);
1934 <    vector<RealType> binPy(nBins_, 0.0);
1935 <    vector<RealType> binPz(nBins_, 0.0);
1936 <    vector<RealType> binOmegax(nBins_, 0.0);
1888 <    vector<RealType> binOmegay(nBins_, 0.0);
1889 <    vector<RealType> binOmegaz(nBins_, 0.0);
1933 >    vector<Vector3d> binP(nBins_, V3Zero);
1934 >    vector<RealType> binOmega(nBins_, 0.0);
1935 >    vector<Vector3d> binL(nBins_, V3Zero);
1936 >    vector<Mat3x3d>  binI(nBins_);
1937      vector<RealType> binKE(nBins_, 0.0);
1938      vector<int> binDOF(nBins_, 0);
1939      vector<int> binCount(nBins_, 0);
# Line 1926 | Line 1973 | namespace OpenMD {
1973          binNo = int(rPos.length() / binWidth_);
1974        }
1975  
1976 <      RealType mass = sd->getMass();
1977 <      Vector3d vel = sd->getVel();
1978 <      Vector3d rPos = sd->getPos() - coordinateOrigin_;
1979 <      Vector3d aVel = cross(rPos, vel);
1980 <      
1976 >      mass = sd->getMass();
1977 >      vel = sd->getVel();
1978 >      rPos = sd->getPos() - coordinateOrigin_;
1979 >      KE = 0.5 * mass * vel.lengthSquare();
1980 >      L = mass * cross(rPos, vel);
1981 >      I = outProduct(rPos, rPos) * mass;
1982 >      r2 = rPos.lengthSquare();
1983 >      I(0, 0) += mass * r2;
1984 >      I(1, 1) += mass * r2;
1985 >      I(2, 2) += mass * r2;
1986 >
1987 >      // Project the relative position onto a plane perpendicular to
1988 >      // the angularMomentumFluxVector:
1989 >      // Vector3d rProj = rPos - dot(rPos, u) * u;
1990 >      // Project the velocity onto a plane perpendicular to the
1991 >      // angularMomentumFluxVector:
1992 >      // Vector3d vProj = vel  - dot(vel, u) * u;
1993 >      // Compute angular velocity vector (should be nearly parallel to
1994 >      // angularMomentumFluxVector
1995 >      // Vector3d aVel = cross(rProj, vProj);
1996 >
1997        if (binNo >= 0 && binNo < nBins_)  {
1998          binCount[binNo]++;
1999          binMass[binNo] += mass;
2000 <        binPx[binNo] += mass*vel.x();
2001 <        binPy[binNo] += mass*vel.y();
2002 <        binPz[binNo] += mass*vel.z();
2003 <        binOmegax[binNo] += aVel.x();
1941 <        binOmegay[binNo] += aVel.y();
1942 <        binOmegaz[binNo] += aVel.z();
1943 <        binKE[binNo] += 0.5 * (mass * vel.lengthSquare());
2000 >        binP[binNo] += mass*vel;
2001 >        binKE[binNo] += KE;
2002 >        binI[binNo] += I;
2003 >        binL[binNo] += L;
2004          binDOF[binNo] += 3;
2005          
2006          if (sd->isDirectional()) {
2007            Vector3d angMom = sd->getJ();
2008 <          Mat3x3d I = sd->getI();
2008 >          Mat3x3d Ia = sd->getI();
2009            if (sd->isLinear()) {
2010              int i = sd->linearAxis();
2011              int j = (i + 1) % 3;
2012              int k = (i + 2) % 3;
2013 <            binKE[binNo] += 0.5 * (angMom[j] * angMom[j] / I(j, j) +
2014 <                                   angMom[k] * angMom[k] / I(k, k));
2013 >            binKE[binNo] += 0.5 * (angMom[j] * angMom[j] / Ia(j, j) +
2014 >                                   angMom[k] * angMom[k] / Ia(k, k));
2015              binDOF[binNo] += 2;
2016            } else {
2017 <            binKE[binNo] += 0.5 * (angMom[0] * angMom[0] / I(0, 0) +
2018 <                                   angMom[1] * angMom[1] / I(1, 1) +
2019 <                                   angMom[2] * angMom[2] / I(2, 2));
2017 >            binKE[binNo] += 0.5 * (angMom[0] * angMom[0] / Ia(0, 0) +
2018 >                                   angMom[1] * angMom[1] / Ia(1, 1) +
2019 >                                   angMom[2] * angMom[2] / Ia(2, 2));
2020              binDOF[binNo] += 3;
2021            }
2022          }
2023        }
2024      }
2025 <    
2025 >
2026   #ifdef IS_MPI
2027 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binCount[0],
2028 <                              nBins_, MPI::INT, MPI::SUM);
2029 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binMass[0],
2030 <                              nBins_, MPI::REALTYPE, MPI::SUM);
2031 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binPx[0],
2032 <                              nBins_, MPI::REALTYPE, MPI::SUM);
2033 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binPy[0],
2034 <                              nBins_, MPI::REALTYPE, MPI::SUM);
2035 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binPz[0],
2036 <                              nBins_, MPI::REALTYPE, MPI::SUM);
2037 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binOmegax[0],
2038 <                              nBins_, MPI::REALTYPE, MPI::SUM);
2039 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binOmegay[0],
2040 <                              nBins_, MPI::REALTYPE, MPI::SUM);
2041 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binOmegaz[0],
2042 <                              nBins_, MPI::REALTYPE, MPI::SUM);
2043 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binKE[0],
2044 <                              nBins_, MPI::REALTYPE, MPI::SUM);
2045 <    MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binDOF[0],
2046 <                              nBins_, MPI::INT, MPI::SUM);
2027 >
2028 >    for (int i = 0; i < nBins_; i++) {
2029 >      
2030 >      MPI_Allreduce(MPI_IN_PLACE, &binCount[i],
2031 >                    1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
2032 >      MPI_Allreduce(MPI_IN_PLACE, &binMass[i],
2033 >                    1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
2034 >      MPI_Allreduce(MPI_IN_PLACE, binP[i].getArrayPointer(),
2035 >                    3, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
2036 >      MPI_Allreduce(MPI_IN_PLACE, binL[i].getArrayPointer(),
2037 >                    3, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
2038 >      MPI_Allreduce(MPI_IN_PLACE, binI[i].getArrayPointer(),
2039 >                    9, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
2040 >      MPI_Allreduce(MPI_IN_PLACE, &binKE[i],
2041 >                    1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
2042 >      MPI_Allreduce(MPI_IN_PLACE, &binDOF[i],
2043 >                    1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
2044 >      //MPI_Allreduce(MPI_IN_PLACE, &binOmega[i],
2045 >      //                          1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
2046 >    }
2047 >    
2048   #endif
2049  
2050 <    Vector3d vel;
1990 <    Vector3d aVel;
2050 >    Vector3d omega;
2051      RealType den;
2052      RealType temp;
2053      RealType z;
# Line 2004 | Line 2064 | namespace OpenMD {
2064          den = binMass[i] * 3.0 * PhysicalConstants::densityConvert
2065            / (4.0 * M_PI * (pow(router,3) - pow(rinner,3)));
2066        }
2067 <      vel.x() = binPx[i] / binMass[i];
2068 <      vel.y() = binPy[i] / binMass[i];
2069 <      vel.z() = binPz[i] / binMass[i];
2070 <      aVel.x() = binOmegax[i] / binCount[i];
2071 <      aVel.y() = binOmegay[i] / binCount[i];
2012 <      aVel.z() = binOmegaz[i] / binCount[i];
2067 >      vel = binP[i] / binMass[i];
2068 >
2069 >      omega = binI[i].inverse() * binL[i];
2070 >
2071 >      // omega = binOmega[i] / binCount[i];
2072  
2073        if (binCount[i] > 0) {
2074          // only add values if there are things to add
# Line 2032 | Line 2091 | namespace OpenMD {
2091                dynamic_cast<VectorAccumulator *>(data_[j].accumulator[i])->add(vel);
2092                break;
2093              case ANGULARVELOCITY:  
2094 <              dynamic_cast<VectorAccumulator *>(data_[j].accumulator[i])->add(aVel);
2094 >              dynamic_cast<VectorAccumulator *>(data_[j].accumulator[i])->add(omega);
2095                break;
2096              case DENSITY:
2097                dynamic_cast<Accumulator *>(data_[j].accumulator[i])->add(den);
# Line 2070 | Line 2129 | namespace OpenMD {
2129          painCave.severity = OPENMD_ERROR;
2130          simError();            
2131        }
2132 <    }  
2132 >    }
2133    }
2134    
2135    void RNEMD::writeOutputFile() {
# Line 2079 | Line 2138 | namespace OpenMD {
2138      
2139   #ifdef IS_MPI
2140      // If we're the root node, should we print out the results
2141 <    int worldRank = MPI::COMM_WORLD.Get_rank();
2141 >    int worldRank;
2142 >    MPI_Comm_rank( MPI_COMM_WORLD, &worldRank);
2143 >
2144      if (worldRank == 0) {
2145   #endif
2146        rnemdFile_.open(rnemdFileName_.c_str(), std::ios::out | std::ios::trunc );
# Line 2210 | Line 2271 | namespace OpenMD {
2271        }        
2272  
2273        rnemdFile_ << "#######################################################\n";
2274 <      rnemdFile_ << "# Standard Deviations in those quantities follow:\n";
2274 >      rnemdFile_ << "# 95% confidence intervals in those quantities follow:\n";
2275        rnemdFile_ << "#######################################################\n";
2276  
2277  
# Line 2219 | Line 2280 | namespace OpenMD {
2280          for (unsigned int i = 0; i < outputMask_.size(); ++i) {
2281            if (outputMask_[i]) {
2282              if (data_[i].dataType == "RealType")
2283 <              writeRealStdDev(i,j);
2283 >              writeRealErrorBars(i,j);
2284              else if (data_[i].dataType == "Vector3d")
2285 <              writeVectorStdDev(i,j);
2285 >              writeVectorErrorBars(i,j);
2286              else {
2287                sprintf( painCave.errMsg,
2288                         "RNEMD found an unknown data type for: %s ",
# Line 2292 | Line 2353 | namespace OpenMD {
2353      }
2354    }  
2355  
2356 <  void RNEMD::writeRealStdDev(int index, unsigned int bin) {
2356 >  void RNEMD::writeRealErrorBars(int index, unsigned int bin) {
2357      if (!doRNEMD_) return;
2358      assert(index >=0 && index < ENDINDEX);
2359      assert(int(bin) < nBins_);
# Line 2302 | Line 2363 | namespace OpenMD {
2363      count = data_[index].accumulator[bin]->count();
2364      if (count == 0) return;
2365      
2366 <    dynamic_cast<Accumulator *>(data_[index].accumulator[bin])->getStdDev(s);
2366 >    dynamic_cast<Accumulator *>(data_[index].accumulator[bin])->get95percentConfidenceInterval(s);
2367      
2368      if (! isinf(s) && ! isnan(s)) {
2369        rnemdFile_ << "\t" << s;
# Line 2315 | Line 2376 | namespace OpenMD {
2376      }    
2377    }
2378    
2379 <  void RNEMD::writeVectorStdDev(int index, unsigned int bin) {
2379 >  void RNEMD::writeVectorErrorBars(int index, unsigned int bin) {
2380      if (!doRNEMD_) return;
2381      assert(index >=0 && index < ENDINDEX);
2382      assert(int(bin) < nBins_);
# Line 2325 | Line 2386 | namespace OpenMD {
2386      count = data_[index].accumulator[bin]->count();
2387      if (count == 0) return;
2388  
2389 <    dynamic_cast<VectorAccumulator*>(data_[index].accumulator[bin])->getStdDev(s);
2389 >    dynamic_cast<VectorAccumulator*>(data_[index].accumulator[bin])->get95percentConfidenceInterval(s);
2390      if (isinf(s[0]) || isnan(s[0]) ||
2391          isinf(s[1]) || isnan(s[1]) ||
2392          isinf(s[2]) || isnan(s[2]) ) {      

Diff Legend

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