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

Comparing branches/development/src/rnemd/RNEMD.cpp (file contents):
Revision 1774 by gezelter, Wed Aug 8 16:03:02 2012 UTC vs.
Revision 1777 by gezelter, Thu Aug 9 18:35:09 2012 UTC

# Line 69 | Line 69 | namespace OpenMD {
69      Globals * simParams = info->getSimParams();
70      RNEMDParameters* rnemdParams = simParams->getRNEMDParameters();
71  
72 +    doRNEMD_ = rnemdParams->getUseRNEMD();
73 +    if (!doRNEMD_) return;
74 +
75      stringToMethod_["Swap"]  = rnemdSwap;
76      stringToMethod_["NIVS"]  = rnemdNIVS;
77      stringToMethod_["VSS"]   = rnemdVSS;
# Line 77 | Line 80 | namespace OpenMD {
80      stringToFluxType_["Px"]  = rnemdPx;
81      stringToFluxType_["Py"]  = rnemdPy;
82      stringToFluxType_["Pz"]  = rnemdPz;
83 +    stringToFluxType_["Pvector"]  = rnemdPvector;
84      stringToFluxType_["KE+Px"]  = rnemdKePx;
85      stringToFluxType_["KE+Py"]  = rnemdKePy;
86      stringToFluxType_["KE+Pvector"]  = rnemdKePvector;
# Line 98 | Line 102 | namespace OpenMD {
102        sprintf(painCave.errMsg,
103                "RNEMD: No fluxType was set in the md file.  This parameter,\n"
104                "\twhich must be one of the following values:\n"
105 <              "\tKE, Px, Py, Pz, KE+Px, KE+Py, KE+Pvector, must be set to\n"
106 <              "\tuse RNEMD\n");
105 >              "\tKE, Px, Py, Pz, Pvector, KE+Px, KE+Py, KE+Pvector\n"
106 >              "\tmust be set to use RNEMD\n");
107        painCave.isFatal = 1;
108        painCave.severity = OPENMD_ERROR;
109        simError();
# Line 197 | Line 201 | namespace OpenMD {
201          break;
202        case rnemdPvector:
203          hasCorrectFlux = hasMomentumFluxVector;
204 +        break;
205        case rnemdKePx:
206        case rnemdKePy:
207          hasCorrectFlux = hasMomentumFlux && hasKineticFlux;
# Line 224 | Line 229 | namespace OpenMD {
229      }
230      if (!hasCorrectFlux) {
231        sprintf(painCave.errMsg,
232 <              "RNEMD: The current method,\n"
228 <              "\t%s, and flux type %s\n"
232 >              "RNEMD: The current method, %s, and flux type, %s,\n"
233                "\tdid not have the correct flux value specified. Options\n"
234                "\tinclude: kineticFlux, momentumFlux, and momentumFluxVector\n",
235                methStr.c_str(), fluxStr.c_str());
# Line 235 | Line 239 | namespace OpenMD {
239      }
240  
241      if (hasKineticFlux) {
242 <      kineticFlux_ = rnemdParams->getKineticFlux();
242 >      // convert the kcal / mol / Angstroms^2 / fs values in the md file
243 >      // into  amu / fs^3:
244 >      kineticFlux_ = rnemdParams->getKineticFlux()
245 >        * PhysicalConstants::energyConvert;
246      } else {
247        kineticFlux_ = 0.0;
248      }
# Line 314 | Line 321 | namespace OpenMD {
321      outputMap_["TEMPERATURE"] =  TEMPERATURE;
322  
323      OutputData velocity;
324 <    velocity.units = "amu/fs";
324 >    velocity.units = "angstroms/fs";
325      velocity.title =  "Velocity";  
326      velocity.dataType = "Vector3d";
327      velocity.accumulator.reserve(nBins_);
# Line 410 | Line 417 | namespace OpenMD {
417    }
418    
419    RNEMD::~RNEMD() {
420 <    
420 >    if (!doRNEMD_) return;
421   #ifdef IS_MPI
422      if (worldRank == 0) {
423   #endif
# Line 432 | Line 439 | namespace OpenMD {
439    }
440  
441    void RNEMD::doSwap() {
442 <
442 >    if (!doRNEMD_) return;
443      Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
444      Mat3x3d hmat = currentSnap_->getHmat();
445  
# Line 491 | Line 498 | namespace OpenMD {
498                  + angMom[2]*angMom[2]/I(2, 2);
499              }
500            } //angular momenta exchange enabled
494 –          //energyConvert temporarily disabled
495 –          //make kineticExchange_ comparable between swap & scale
496 –          //value = value * 0.5 / PhysicalConstants::energyConvert;
501            value *= 0.5;
502            break;
503          case rnemdPx :
# Line 731 | Line 735 | namespace OpenMD {
735          
736          switch(rnemdFluxType_) {
737          case rnemdKE:
734 –          cerr << "KE\n";
738            kineticExchange_ += max_val - min_val;
739            break;
740          case rnemdPx:
# Line 744 | Line 747 | namespace OpenMD {
747            momentumExchange_.z() += max_val - min_val;
748            break;
749          default:
747 –          cerr << "default\n";
750            break;
751          }
752        } else {        
# Line 767 | Line 769 | namespace OpenMD {
769    }
770    
771    void RNEMD::doNIVS() {
772 <
772 >    if (!doRNEMD_) return;
773      Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
774      Mat3x3d hmat = currentSnap_->getHmat();
775  
# Line 1216 | Line 1218 | namespace OpenMD {
1218    }
1219  
1220    void RNEMD::doVSS() {
1221 <
1221 >    if (!doRNEMD_) return;
1222      Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
1223      RealType time = currentSnap_->getTime();    
1224      Mat3x3d hmat = currentSnap_->getHmat();
# Line 1394 | Line 1396 | namespace OpenMD {
1396    }
1397  
1398    void RNEMD::doRNEMD() {
1399 <
1399 >    if (!doRNEMD_) return;
1400      trialCount_++;
1401      switch(rnemdMethod_) {
1402      case rnemdSwap:
# Line 1413 | Line 1415 | namespace OpenMD {
1415    }
1416  
1417    void RNEMD::collectData() {
1418 <
1418 >    if (!doRNEMD_) return;
1419      Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
1420      Mat3x3d hmat = currentSnap_->getHmat();
1421  
# Line 1524 | Line 1526 | namespace OpenMD {
1526        vel.x() = binPx[i] / binMass[i];
1527        vel.y() = binPy[i] / binMass[i];
1528        vel.z() = binPz[i] / binMass[i];
1529 <      den = binCount[i] * nBins_ / currentSnap_->getVolume();
1529 >
1530 >      den = binMass[i] * nBins_ * PhysicalConstants::densityConvert
1531 >        / currentSnap_->getVolume() ;
1532 >
1533        temp = 2.0 * binKE[i] / (binDOF[i] * PhysicalConstants::kb *
1534                                 PhysicalConstants::energyConvert);
1535  
# Line 1550 | Line 1555 | namespace OpenMD {
1555    }
1556  
1557    void RNEMD::getStarted() {
1558 +    if (!doRNEMD_) return;
1559      collectData();
1560      writeOutputFile();
1561    }
1562  
1563    void RNEMD::parseOutputFileFormat(const std::string& format) {
1564 +    if (!doRNEMD_) return;
1565      StringTokenizer tokenizer(format, " ,;|\t\n\r");
1566      
1567      while(tokenizer.hasMoreTokens()) {
# Line 1575 | Line 1582 | namespace OpenMD {
1582    }
1583    
1584    void RNEMD::writeOutputFile() {
1585 +    if (!doRNEMD_) return;
1586      
1587   #ifdef IS_MPI
1588      // If we're the root node, should we print out the results
# Line 1596 | Line 1604 | namespace OpenMD {
1604        RealType time = currentSnap_->getTime();
1605        RealType avgArea;
1606        areaAccumulator_->getAverage(avgArea);
1607 <      RealType Jz = kineticExchange_ / (2.0 * time * avgArea);
1607 >      RealType Jz = kineticExchange_ / (2.0 * time * avgArea)
1608 >        / PhysicalConstants::energyConvert;
1609        Vector3d JzP = momentumExchange_ / (2.0 * time * avgArea);      
1610  
1611        rnemdFile_ << "#######################################################\n";
# Line 1613 | Line 1622 | namespace OpenMD {
1622            rnemdFile_ << "#    fluxType  = \"" << (*fi).first << "\";\n";
1623        }
1624        
1625 <      rnemdFile_ << "#    exchangeTime = " << exchangeTime_ << "; fs\n";
1625 >      rnemdFile_ << "#    exchangeTime = " << exchangeTime_ << ";\n";
1626  
1627        rnemdFile_ << "#    objectSelection = \""
1628                   << rnemdObjectSelection_ << "\";\n";
# Line 1625 | Line 1634 | namespace OpenMD {
1634        rnemdFile_ << "# RNEMD report:\n";      
1635        rnemdFile_ << "#     running time = " << time << " fs\n";
1636        rnemdFile_ << "#     target flux:\n";
1637 <      rnemdFile_ << "#         kinetic = " << kineticFlux_ << "\n";
1638 <      rnemdFile_ << "#         momentum = " << momentumFluxVector_ << "\n";
1637 >      rnemdFile_ << "#         kinetic = "
1638 >                 << kineticFlux_ / PhysicalConstants::energyConvert
1639 >                 << " (kcal/mol/A^2/fs)\n";
1640 >      rnemdFile_ << "#         momentum = " << momentumFluxVector_
1641 >                 << " (amu/A/fs^2)\n";
1642        rnemdFile_ << "#     target one-time exchanges:\n";
1643 <      rnemdFile_ << "#         kinetic = " << kineticTarget_ << "\n";
1644 <      rnemdFile_ << "#         momentum = " << momentumTarget_ << "\n";
1643 >      rnemdFile_ << "#         kinetic = "
1644 >                 << kineticTarget_ / PhysicalConstants::energyConvert
1645 >                 << " (kcal/mol)\n";
1646 >      rnemdFile_ << "#         momentum = " << momentumTarget_
1647 >                 << " (amu*A/fs)\n";
1648        rnemdFile_ << "#     actual exchange totals:\n";
1649 <      rnemdFile_ << "#         kinetic = " << kineticExchange_ << "\n";
1650 <      rnemdFile_ << "#         momentum = " << momentumExchange_  << "\n";
1649 >      rnemdFile_ << "#         kinetic = "
1650 >                 << kineticExchange_ / PhysicalConstants::energyConvert
1651 >                 << " (kcal/mol)\n";
1652 >      rnemdFile_ << "#         momentum = " << momentumExchange_
1653 >                 << " (amu*A/fs)\n";      
1654        rnemdFile_ << "#     actual flux:\n";
1655 <      rnemdFile_ << "#         kinetic = " << Jz << "\n";
1656 <      rnemdFile_ << "#         momentum = " << JzP  << "\n";
1655 >      rnemdFile_ << "#         kinetic = " << Jz
1656 >                 << " (kcal/mol/A^2/fs)\n";
1657 >      rnemdFile_ << "#         momentum = " << JzP
1658 >                 << " (amu/A/fs^2)\n";
1659        rnemdFile_ << "#     exchange statistics:\n";
1660        rnemdFile_ << "#         attempted = " << trialCount_ << "\n";
1661        rnemdFile_ << "#         failed = " << failTrialCount_ << "\n";    
# Line 1653 | Line 1673 | namespace OpenMD {
1673          if (outputMask_[i]) {
1674            rnemdFile_ << "\t" << data_[i].title <<
1675              "(" << data_[i].units << ")";
1676 +          // add some extra tabs for column alignment
1677 +          if (data_[i].dataType == "Vector3d") rnemdFile_ << "\t\t";
1678          }
1679        }
1680        rnemdFile_ << std::endl;
# Line 1716 | Line 1738 | namespace OpenMD {
1738    }
1739    
1740    void RNEMD::writeReal(int index, unsigned int bin) {
1741 +    if (!doRNEMD_) return;
1742      assert(index >=0 && index < ENDINDEX);
1743      assert(bin < nBins_);
1744      RealType s;
# Line 1734 | Line 1757 | namespace OpenMD {
1757    }
1758    
1759    void RNEMD::writeVector(int index, unsigned int bin) {
1760 +    if (!doRNEMD_) return;
1761      assert(index >=0 && index < ENDINDEX);
1762      assert(bin < nBins_);
1763      Vector3d s;
# Line 1752 | Line 1776 | namespace OpenMD {
1776    }  
1777  
1778    void RNEMD::writeRealStdDev(int index, unsigned int bin) {
1779 +    if (!doRNEMD_) return;
1780      assert(index >=0 && index < ENDINDEX);
1781      assert(bin < nBins_);
1782      RealType s;
# Line 1770 | Line 1795 | namespace OpenMD {
1795    }
1796    
1797    void RNEMD::writeVectorStdDev(int index, unsigned int bin) {
1798 +    if (!doRNEMD_) return;
1799      assert(index >=0 && index < ENDINDEX);
1800      assert(bin < nBins_);
1801      Vector3d s;

Diff Legend

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