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 1980 by gezelter, Sat Apr 5 21:17:12 2014 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 315 | Line 328 | namespace OpenMD {
328          default:
329            break;
330          }
318 –      }
319 –      if (hasAngularMomentumFluxVector) {
320 –        angularMomentumFluxVector_ = rnemdParams->getAngularMomentumFluxVector();
321 –      } else {
322 –        angularMomentumFluxVector_ = V3Zero;
323 –        if (hasAngularMomentumFlux) {
324 –          RealType angularMomentumFlux = rnemdParams->getAngularMomentumFlux();
325 –          switch (rnemdFluxType_) {
326 –          case rnemdLx:
327 –            angularMomentumFluxVector_.x() = angularMomentumFlux;
328 –            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 –        }        
331        }
332 <
333 <      if (hasCoordinateOrigin) {
334 <        coordinateOrigin_ = rnemdParams->getCoordinateOrigin();
335 <      } else {
336 <        coordinateOrigin_ = V3Zero;
337 <      }
338 <
339 <      // do some sanity checking
340 <
341 <      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();
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 <      areaAccumulator_ = new Accumulator();
345 <
346 <      nBins_ = rnemdParams->getOutputBins();
347 <      binWidth_ = rnemdParams->getOutputBinWidth();
348 <
349 <      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 <      
494 <      if (hasOutputFileName) {
495 <        rnemdFileName_ = rnemdParams->getOutputFileName();
496 <      } else {
497 <        rnemdFileName_ = getPrefix(info->getFinalConfigFileName()) + ".rnemd";
498 <      }          
372 >      }        
373 >    }
374  
375 <      exchangeTime_ = rnemdParams->getExchangeTime();
376 <
377 <      Snapshot* currentSnap_ = info->getSnapshotManager()->getCurrentSnapshot();
378 <      // total exchange sums are zeroed out at the beginning:
379 <
380 <      kineticExchange_ = 0.0;
381 <      momentumExchange_ = V3Zero;
382 <      angularMomentumExchange_ = V3Zero;
383 <
384 <      std::ostringstream selectionAstream;
385 <      std::ostringstream selectionBstream;
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 >      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();
397 <        
398 <          if (hasSlabWidth)
399 <            slabWidth_ = rnemdParams->getSlabWidth();
400 <          else
401 <            slabWidth_ = hmat(2,2) / 10.0;
402 <        
403 <          if (hasSlabACenter)
404 <            slabACenter_ = rnemdParams->getSlabACenter();
405 <          else
406 <            slabACenter_ = 0.0;
407 <        
408 <          selectionAstream << "select wrappedz > "
409 <                           << slabACenter_ - 0.5*slabWidth_
410 <                           <<  " && wrappedz < "
411 <                           << slabACenter_ + 0.5*slabWidth_;
412 <          selectionA_ = selectionAstream.str();
413 <        } else {
414 <          if (hasSphereARadius)
415 <            sphereARadius_ = rnemdParams->getSphereARadius();
416 <          else {
417 <            // use an initial guess to the size of the inner slab to be 1/10 the
418 <            // radius of an approximately spherical hull:
419 <            Thermo thermo(info);
420 <            RealType hVol = thermo.getHullVolume();
421 <            sphereARadius_ = 0.1 * pow((3.0 * hVol / (4.0 * M_PI)), 1.0/3.0);
422 <          }
423 <          selectionAstream << "select r < " << sphereARadius_;
424 <          selectionA_ = selectionAstream.str();
425 <        }
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 (hasSelectionB_) {
529 <        selectionB_ = rnemdParams->getSelectionB();
530 <
531 <      } else {
532 <        if (usePeriodicBoundaryConditions_) {    
533 <          Mat3x3d hmat = currentSnap_->getHmat();
534 <        
535 <          if (hasSlabWidth)
536 <            slabWidth_ = rnemdParams->getSlabWidth();
537 <          else
538 <            slabWidth_ = hmat(2,2) / 10.0;
539 <        
540 <          if (hasSlabBCenter)
541 <            slabBCenter_ = rnemdParams->getSlabBCenter();
542 <          else
543 <            slabBCenter_ = hmat(2,2) / 2.0;
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 <          selectionBstream << "select wrappedz > "
553 <                           << slabBCenter_ - 0.5*slabWidth_
554 <                           <<  " && wrappedz < "
555 <                           << slabBCenter_ + 0.5*slabWidth_;
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;
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);
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();
588 >        
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;
598 >        
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;
659 >    RealType min_val(0.0);
660      int min_found = 0;  
661 <    StuntDouble* min_sd;
661 >    StuntDouble* min_sd = NULL;
662  
663 <    RealType max_val;
663 >    RealType max_val(0.0);
664      int max_found = 0;
665 <    StuntDouble* max_sd;
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 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 988 | Line 1022 | namespace OpenMD {
1022      int selej;
1023  
1024      Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
991 –    RealType time = currentSnap_->getTime();    
1025      Mat3x3d hmat = currentSnap_->getHmat();
1026  
1027      StuntDouble* sd;
# Line 1115 | Line 1148 | namespace OpenMD {
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 1185 | 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 1288 | 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 1296 | 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 1311 | 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 1436 | Line 1471 | namespace OpenMD {
1471      int selej;
1472  
1473      Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
1439 –    RealType time = currentSnap_->getTime();    
1474      Mat3x3d hmat = currentSnap_->getHmat();
1475  
1476      StuntDouble* sd;
# Line 1743 | Line 1777 | namespace OpenMD {
1777          areaA = evaluatorA_.getSurfaceArea();
1778        else {
1779          
1746 –        cerr << "selection A did not have surface area, recomputing\n";
1780          int isd;
1781          StuntDouble* sd;
1782          vector<StuntDouble*> aSites;
# Line 1781 | Line 1814 | namespace OpenMD {
1814      }
1815  
1816      if (hasSelectionB_) {
1817 <      if (evaluatorB_.hasSurfaceArea())
1817 >      if (evaluatorB_.hasSurfaceArea()) {
1818          areaB = evaluatorB_.getSurfaceArea();
1819 <      else {
1787 <        cerr << "selection B did not have surface area, recomputing\n";
1819 >      } else {
1820  
1821          int isd;
1822          StuntDouble* sd;
# Line 1875 | 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();
# Line 1990 | Line 2022 | namespace OpenMD {
2022          }
2023        }
2024      }
2025 <    
2025 >
2026   #ifdef IS_MPI
2027  
2028      for (int i = 0; i < nBins_; i++) {
# Line 2097 | Line 2129 | namespace OpenMD {
2129          painCave.severity = OPENMD_ERROR;
2130          simError();            
2131        }
2132 <    }  
2132 >    }
2133    }
2134    
2135    void RNEMD::writeOutputFile() {

Diff Legend

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