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

Comparing:
branches/development/src/rnemd/RNEMD.cpp (file contents), Revision 1773 by gezelter, Tue Aug 7 18:26:40 2012 UTC vs.
trunk/src/rnemd/RNEMD.cpp (file contents), Revision 1789 by gezelter, Wed Aug 29 20:52:19 2012 UTC

# Line 53 | Line 53
53   #include <mpi.h>
54   #endif
55  
56 + #ifdef _MSC_VER
57 + #define isnan(x) _isnan((x))
58 + #define isinf(x) (!_finite(x) && !_isnan(x))
59 + #endif
60 +
61   #define HONKING_LARGE_VALUE 1.0e10
62  
63   using namespace std;
# Line 65 | Line 70 | namespace OpenMD {
70      failTrialCount_ = 0;
71      failRootCount_ = 0;
72  
68 –    int seedValue;
73      Globals * simParams = info->getSimParams();
74      RNEMDParameters* rnemdParams = simParams->getRNEMDParameters();
75  
76 +    doRNEMD_ = rnemdParams->getUseRNEMD();
77 +    if (!doRNEMD_) return;
78 +
79      stringToMethod_["Swap"]  = rnemdSwap;
80      stringToMethod_["NIVS"]  = rnemdNIVS;
81      stringToMethod_["VSS"]   = rnemdVSS;
# Line 77 | Line 84 | namespace OpenMD {
84      stringToFluxType_["Px"]  = rnemdPx;
85      stringToFluxType_["Py"]  = rnemdPy;
86      stringToFluxType_["Pz"]  = rnemdPz;
87 +    stringToFluxType_["Pvector"]  = rnemdPvector;
88      stringToFluxType_["KE+Px"]  = rnemdKePx;
89      stringToFluxType_["KE+Py"]  = rnemdKePy;
90      stringToFluxType_["KE+Pvector"]  = rnemdKePvector;
# Line 98 | Line 106 | namespace OpenMD {
106        sprintf(painCave.errMsg,
107                "RNEMD: No fluxType was set in the md file.  This parameter,\n"
108                "\twhich must be one of the following values:\n"
109 <              "\tKE, Px, Py, Pz, KE+Px, KE+Py, KE+Pvector, must be set to\n"
110 <              "\tuse RNEMD\n");
109 >              "\tKE, Px, Py, Pz, Pvector, KE+Px, KE+Py, KE+Pvector\n"
110 >              "\tmust be set to use RNEMD\n");
111        painCave.isFatal = 1;
112        painCave.severity = OPENMD_ERROR;
113        simError();
# Line 197 | Line 205 | namespace OpenMD {
205          break;
206        case rnemdPvector:
207          hasCorrectFlux = hasMomentumFluxVector;
208 +        break;
209        case rnemdKePx:
210        case rnemdKePy:
211          hasCorrectFlux = hasMomentumFlux && hasKineticFlux;
# Line 224 | Line 233 | namespace OpenMD {
233      }
234      if (!hasCorrectFlux) {
235        sprintf(painCave.errMsg,
236 <              "RNEMD: The current method,\n"
228 <              "\t%s, and flux type %s\n"
236 >              "RNEMD: The current method, %s, and flux type, %s,\n"
237                "\tdid not have the correct flux value specified. Options\n"
238                "\tinclude: kineticFlux, momentumFlux, and momentumFluxVector\n",
239                methStr.c_str(), fluxStr.c_str());
# Line 235 | Line 243 | namespace OpenMD {
243      }
244  
245      if (hasKineticFlux) {
246 <      kineticFlux_ = rnemdParams->getKineticFlux();
246 >      // convert the kcal / mol / Angstroms^2 / fs values in the md file
247 >      // into  amu / fs^3:
248 >      kineticFlux_ = rnemdParams->getKineticFlux()
249 >        * PhysicalConstants::energyConvert;
250      } else {
251        kineticFlux_ = 0.0;
252      }
# Line 288 | Line 299 | namespace OpenMD {
299        simError();
300      }
301  
302 +    areaAccumulator_ = new Accumulator();
303 +
304      nBins_ = rnemdParams->getOutputBins();
305  
306      data_.resize(RNEMD::ENDINDEX);
# Line 296 | Line 309 | namespace OpenMD {
309      z.title =  "Z";
310      z.dataType = "RealType";
311      z.accumulator.reserve(nBins_);
312 <    for (unsigned int i = 0; i < nBins_; i++)
312 >    for (int i = 0; i < nBins_; i++)
313        z.accumulator.push_back( new Accumulator() );
314      data_[Z] = z;
315      outputMap_["Z"] =  Z;
# Line 306 | Line 319 | namespace OpenMD {
319      temperature.title =  "Temperature";
320      temperature.dataType = "RealType";
321      temperature.accumulator.reserve(nBins_);
322 <    for (unsigned int i = 0; i < nBins_; i++)
322 >    for (int i = 0; i < nBins_; i++)
323        temperature.accumulator.push_back( new Accumulator() );
324      data_[TEMPERATURE] = temperature;
325      outputMap_["TEMPERATURE"] =  TEMPERATURE;
326  
327      OutputData velocity;
328 <    velocity.units = "amu/fs";
328 >    velocity.units = "angstroms/fs";
329      velocity.title =  "Velocity";  
330      velocity.dataType = "Vector3d";
331      velocity.accumulator.reserve(nBins_);
332 <    for (unsigned int i = 0; i < nBins_; i++)
332 >    for (int i = 0; i < nBins_; i++)
333        velocity.accumulator.push_back( new VectorAccumulator() );
334      data_[VELOCITY] = velocity;
335      outputMap_["VELOCITY"] = VELOCITY;
# Line 326 | Line 339 | namespace OpenMD {
339      density.title =  "Density";
340      density.dataType = "RealType";
341      density.accumulator.reserve(nBins_);
342 <    for (unsigned int i = 0; i < nBins_; i++)
342 >    for (int i = 0; i < nBins_; i++)
343        density.accumulator.push_back( new Accumulator() );
344      data_[DENSITY] = density;
345      outputMap_["DENSITY"] =  DENSITY;
# Line 381 | Line 394 | namespace OpenMD {
394      // dt = exchange time interval
395      // flux = target flux
396  
397 <    kineticTarget_ = 2.0*kineticFlux_*exchangeTime_*hmat(0,0)*hmat(1,1);
398 <    momentumTarget_ = 2.0*momentumFluxVector_*exchangeTime_*hmat(0,0)*hmat(1,1);
397 >    RealType area = currentSnap_->getXYarea();
398 >    kineticTarget_ = 2.0 * kineticFlux_ * exchangeTime_ * area;
399 >    momentumTarget_ = 2.0 * momentumFluxVector_ * exchangeTime_ * area;
400  
401      // total exchange sums are zeroed out at the beginning:
402  
# Line 407 | Line 421 | namespace OpenMD {
421    }
422    
423    RNEMD::~RNEMD() {
424 <    
424 >    if (!doRNEMD_) return;
425   #ifdef IS_MPI
426      if (worldRank == 0) {
427   #endif
# Line 429 | Line 443 | namespace OpenMD {
443    }
444  
445    void RNEMD::doSwap() {
446 <
446 >    if (!doRNEMD_) return;
447      Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
448      Mat3x3d hmat = currentSnap_->getHmat();
449  
# Line 488 | Line 502 | namespace OpenMD {
502                  + angMom[2]*angMom[2]/I(2, 2);
503              }
504            } //angular momenta exchange enabled
491 –          //energyConvert temporarily disabled
492 –          //make kineticExchange_ comparable between swap & scale
493 –          //value = value * 0.5 / PhysicalConstants::energyConvert;
505            value *= 0.5;
506            break;
507          case rnemdPx :
# Line 728 | Line 739 | namespace OpenMD {
739          
740          switch(rnemdFluxType_) {
741          case rnemdKE:
731 –          cerr << "KE\n";
742            kineticExchange_ += max_val - min_val;
743            break;
744          case rnemdPx:
# Line 741 | Line 751 | namespace OpenMD {
751            momentumExchange_.z() += max_val - min_val;
752            break;
753          default:
744 –          cerr << "default\n";
754            break;
755          }
756        } else {        
# Line 764 | Line 773 | namespace OpenMD {
773    }
774    
775    void RNEMD::doNIVS() {
776 <
776 >    if (!doRNEMD_) return;
777      Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
778      Mat3x3d hmat = currentSnap_->getHmat();
779  
# Line 1213 | Line 1222 | namespace OpenMD {
1222    }
1223  
1224    void RNEMD::doVSS() {
1225 <
1225 >    if (!doRNEMD_) return;
1226      Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
1227      RealType time = currentSnap_->getTime();    
1228      Mat3x3d hmat = currentSnap_->getHmat();
# Line 1391 | Line 1400 | namespace OpenMD {
1400    }
1401  
1402    void RNEMD::doRNEMD() {
1403 <
1403 >    if (!doRNEMD_) return;
1404      trialCount_++;
1405      switch(rnemdMethod_) {
1406      case rnemdSwap:
# Line 1410 | Line 1419 | namespace OpenMD {
1419    }
1420  
1421    void RNEMD::collectData() {
1422 <
1422 >    if (!doRNEMD_) return;
1423      Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
1424      Mat3x3d hmat = currentSnap_->getHmat();
1425  
1426 +    areaAccumulator_->add(currentSnap_->getXYarea());
1427 +
1428      seleMan_.setSelectionSet(evaluator_.evaluate());
1429  
1430      int selei;
# Line 1454 | Line 1465 | namespace OpenMD {
1465        if (usePeriodicBoundaryConditions_)
1466          currentSnap_->wrapVector(pos);
1467  
1468 +
1469        // which bin is this stuntdouble in?
1470        // wrapped positions are in the range [-0.5*hmat(2,2), +0.5*hmat(2,2)]
1471        // Shift molecules by half a box to have bins start at 0
# Line 1518 | Line 1530 | namespace OpenMD {
1530        vel.x() = binPx[i] / binMass[i];
1531        vel.y() = binPy[i] / binMass[i];
1532        vel.z() = binPz[i] / binMass[i];
1533 <      den = binCount[i] * nBins_ / (hmat(0,0) * hmat(1,1) * hmat(2,2));
1533 >
1534 >      den = binMass[i] * nBins_ * PhysicalConstants::densityConvert
1535 >        / currentSnap_->getVolume() ;
1536 >
1537        temp = 2.0 * binKE[i] / (binDOF[i] * PhysicalConstants::kb *
1538                                 PhysicalConstants::energyConvert);
1539  
# Line 1544 | Line 1559 | namespace OpenMD {
1559    }
1560  
1561    void RNEMD::getStarted() {
1562 +    if (!doRNEMD_) return;
1563      collectData();
1564      writeOutputFile();
1565    }
1566  
1567    void RNEMD::parseOutputFileFormat(const std::string& format) {
1568 +    if (!doRNEMD_) return;
1569      StringTokenizer tokenizer(format, " ,;|\t\n\r");
1570      
1571      while(tokenizer.hasMoreTokens()) {
# Line 1569 | Line 1586 | namespace OpenMD {
1586    }
1587    
1588    void RNEMD::writeOutputFile() {
1589 +    if (!doRNEMD_) return;
1590      
1591   #ifdef IS_MPI
1592      // If we're the root node, should we print out the results
# Line 1588 | Line 1606 | namespace OpenMD {
1606        Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
1607  
1608        RealType time = currentSnap_->getTime();
1609 <      
1610 <      
1609 >      RealType avgArea;
1610 >      areaAccumulator_->getAverage(avgArea);
1611 >      RealType Jz = kineticExchange_ / (2.0 * time * avgArea)
1612 >        / PhysicalConstants::energyConvert;
1613 >      Vector3d JzP = momentumExchange_ / (2.0 * time * avgArea);      
1614 >
1615        rnemdFile_ << "#######################################################\n";
1616        rnemdFile_ << "# RNEMD {\n";
1617  
1618        map<string, RNEMDMethod>::iterator mi;
1619        for(mi = stringToMethod_.begin(); mi != stringToMethod_.end(); ++mi) {
1620          if ( (*mi).second == rnemdMethod_)
1621 <          rnemdFile_ << "#    exchangeMethod  = " << (*mi).first << "\n";
1621 >          rnemdFile_ << "#    exchangeMethod  = \"" << (*mi).first << "\";\n";
1622        }
1623        map<string, RNEMDFluxType>::iterator fi;
1624        for(fi = stringToFluxType_.begin(); fi != stringToFluxType_.end(); ++fi) {
1625          if ( (*fi).second == rnemdFluxType_)
1626 <          rnemdFile_ << "#    fluxType  = " << (*fi).first << "\n";
1626 >          rnemdFile_ << "#    fluxType  = \"" << (*fi).first << "\";\n";
1627        }
1628        
1629 <      rnemdFile_ << "#    exchangeTime = " << exchangeTime_ << " fs\n";
1629 >      rnemdFile_ << "#    exchangeTime = " << exchangeTime_ << ";\n";
1630  
1631        rnemdFile_ << "#    objectSelection = \""
1632 <                 << rnemdObjectSelection_ << "\"\n";
1633 <      rnemdFile_ << "#    slabWidth = " << slabWidth_ << " angstroms\n";
1634 <      rnemdFile_ << "#    slabAcenter = " << slabACenter_ << " angstroms\n";
1635 <      rnemdFile_ << "#    slabBcenter = " << slabBCenter_ << " angstroms\n";
1632 >                 << rnemdObjectSelection_ << "\";\n";
1633 >      rnemdFile_ << "#    slabWidth = " << slabWidth_ << ";\n";
1634 >      rnemdFile_ << "#    slabAcenter = " << slabACenter_ << ";\n";
1635 >      rnemdFile_ << "#    slabBcenter = " << slabBCenter_ << ";\n";
1636        rnemdFile_ << "# }\n";
1637        rnemdFile_ << "#######################################################\n";
1638 <      
1639 <      rnemdFile_ << "# running time = " << time << " fs\n";
1640 <      rnemdFile_ << "# target kinetic flux = " << kineticFlux_ << "\n";
1641 <      rnemdFile_ << "# target momentum flux = " << momentumFluxVector_ << "\n";
1642 <      
1643 <      rnemdFile_ << "# target one-time kinetic exchange = " << kineticTarget_
1644 <                 << "\n";
1645 <      rnemdFile_ << "# target one-time momentum exchange = " << momentumTarget_
1646 <                 << "\n";
1647 <      
1648 <      rnemdFile_ << "# actual kinetic exchange = " << kineticExchange_ << "\n";
1649 <      rnemdFile_ << "# actual momentum exchange = " << momentumExchange_
1650 <                 << "\n";
1651 <      
1652 <      rnemdFile_ << "# attempted exchanges: " << trialCount_ << "\n";
1653 <      rnemdFile_ << "# failed exchanges: " << failTrialCount_ << "\n";
1654 <
1655 <      
1638 >      rnemdFile_ << "# RNEMD report:\n";      
1639 >      rnemdFile_ << "#     running time = " << time << " fs\n";
1640 >      rnemdFile_ << "#     target flux:\n";
1641 >      rnemdFile_ << "#         kinetic = "
1642 >                 << kineticFlux_ / PhysicalConstants::energyConvert
1643 >                 << " (kcal/mol/A^2/fs)\n";
1644 >      rnemdFile_ << "#         momentum = " << momentumFluxVector_
1645 >                 << " (amu/A/fs^2)\n";
1646 >      rnemdFile_ << "#     target one-time exchanges:\n";
1647 >      rnemdFile_ << "#         kinetic = "
1648 >                 << kineticTarget_ / PhysicalConstants::energyConvert
1649 >                 << " (kcal/mol)\n";
1650 >      rnemdFile_ << "#         momentum = " << momentumTarget_
1651 >                 << " (amu*A/fs)\n";
1652 >      rnemdFile_ << "#     actual exchange totals:\n";
1653 >      rnemdFile_ << "#         kinetic = "
1654 >                 << kineticExchange_ / PhysicalConstants::energyConvert
1655 >                 << " (kcal/mol)\n";
1656 >      rnemdFile_ << "#         momentum = " << momentumExchange_
1657 >                 << " (amu*A/fs)\n";      
1658 >      rnemdFile_ << "#     actual flux:\n";
1659 >      rnemdFile_ << "#         kinetic = " << Jz
1660 >                 << " (kcal/mol/A^2/fs)\n";
1661 >      rnemdFile_ << "#         momentum = " << JzP
1662 >                 << " (amu/A/fs^2)\n";
1663 >      rnemdFile_ << "#     exchange statistics:\n";
1664 >      rnemdFile_ << "#         attempted = " << trialCount_ << "\n";
1665 >      rnemdFile_ << "#         failed = " << failTrialCount_ << "\n";    
1666        if (rnemdMethod_ == rnemdNIVS) {
1667 <        rnemdFile_ << "# NIVS root-check warnings: " << failRootCount_ << "\n";
1667 >        rnemdFile_ << "#         NIVS root-check errors = "
1668 >                   << failRootCount_ << "\n";
1669        }
1637 –
1670        rnemdFile_ << "#######################################################\n";
1671        
1672        
# Line 1645 | Line 1677 | namespace OpenMD {
1677          if (outputMask_[i]) {
1678            rnemdFile_ << "\t" << data_[i].title <<
1679              "(" << data_[i].units << ")";
1680 +          // add some extra tabs for column alignment
1681 +          if (data_[i].dataType == "Vector3d") rnemdFile_ << "\t\t";
1682          }
1683        }
1684        rnemdFile_ << std::endl;
1685        
1686        rnemdFile_.precision(8);
1687        
1688 <      for (unsigned int j = 0; j < nBins_; j++) {        
1688 >      for (int j = 0; j < nBins_; j++) {        
1689          
1690          for (unsigned int i = 0; i < outputMask_.size(); ++i) {
1691            if (outputMask_[i]) {
# Line 1671 | Line 1705 | namespace OpenMD {
1705          rnemdFile_ << std::endl;
1706          
1707        }        
1708 +
1709 +      rnemdFile_ << "#######################################################\n";
1710 +      rnemdFile_ << "# Standard Deviations in those quantities follow:\n";
1711 +      rnemdFile_ << "#######################################################\n";
1712 +
1713 +
1714 +      for (int j = 0; j < nBins_; j++) {        
1715 +        rnemdFile_ << "#";
1716 +        for (unsigned int i = 0; i < outputMask_.size(); ++i) {
1717 +          if (outputMask_[i]) {
1718 +            if (data_[i].dataType == "RealType")
1719 +              writeRealStdDev(i,j);
1720 +            else if (data_[i].dataType == "Vector3d")
1721 +              writeVectorStdDev(i,j);
1722 +            else {
1723 +              sprintf( painCave.errMsg,
1724 +                       "RNEMD found an unknown data type for: %s ",
1725 +                       data_[i].title.c_str());
1726 +              painCave.isFatal = 1;
1727 +              simError();
1728 +            }
1729 +          }
1730 +        }
1731 +        rnemdFile_ << std::endl;
1732 +        
1733 +      }        
1734        
1735        rnemdFile_.flush();
1736        rnemdFile_.close();
# Line 1682 | Line 1742 | namespace OpenMD {
1742    }
1743    
1744    void RNEMD::writeReal(int index, unsigned int bin) {
1745 +    if (!doRNEMD_) return;
1746      assert(index >=0 && index < ENDINDEX);
1747 <    assert(bin >=0 && bin < nBins_);
1747 >    assert(bin < nBins_);
1748      RealType s;
1749      
1750      data_[index].accumulator[bin]->getAverage(s);
# Line 1700 | Line 1761 | namespace OpenMD {
1761    }
1762    
1763    void RNEMD::writeVector(int index, unsigned int bin) {
1764 +    if (!doRNEMD_) return;
1765      assert(index >=0 && index < ENDINDEX);
1766 <    assert(bin >=0 && bin < nBins_);
1766 >    assert(bin < nBins_);
1767      Vector3d s;
1768      dynamic_cast<VectorAccumulator*>(data_[index].accumulator[bin])->getAverage(s);
1769      if (isinf(s[0]) || isnan(s[0]) ||
# Line 1716 | Line 1778 | namespace OpenMD {
1778        rnemdFile_ << "\t" << s[0] << "\t" << s[1] << "\t" << s[2];
1779      }
1780    }  
1781 +
1782 +  void RNEMD::writeRealStdDev(int index, unsigned int bin) {
1783 +    if (!doRNEMD_) return;
1784 +    assert(index >=0 && index < ENDINDEX);
1785 +    assert(bin < nBins_);
1786 +    RealType s;
1787 +    
1788 +    data_[index].accumulator[bin]->getStdDev(s);
1789 +    
1790 +    if (! isinf(s) && ! isnan(s)) {
1791 +      rnemdFile_ << "\t" << s;
1792 +    } else{
1793 +      sprintf( painCave.errMsg,
1794 +               "RNEMD detected a numerical error writing: %s std. dev. for bin %d",
1795 +               data_[index].title.c_str(), bin);
1796 +      painCave.isFatal = 1;
1797 +      simError();
1798 +    }    
1799 +  }
1800 +  
1801 +  void RNEMD::writeVectorStdDev(int index, unsigned int bin) {
1802 +    if (!doRNEMD_) return;
1803 +    assert(index >=0 && index < ENDINDEX);
1804 +    assert(bin < nBins_);
1805 +    Vector3d s;
1806 +    dynamic_cast<VectorAccumulator*>(data_[index].accumulator[bin])->getStdDev(s);
1807 +    if (isinf(s[0]) || isnan(s[0]) ||
1808 +        isinf(s[1]) || isnan(s[1]) ||
1809 +        isinf(s[2]) || isnan(s[2]) ) {      
1810 +      sprintf( painCave.errMsg,
1811 +               "RNEMD detected a numerical error writing: %s std. dev. for bin %d",
1812 +               data_[index].title.c_str(), bin);
1813 +      painCave.isFatal = 1;
1814 +      simError();
1815 +    } else {
1816 +      rnemdFile_ << "\t" << s[0] << "\t" << s[1] << "\t" << s[2];
1817 +    }
1818 +  }  
1819   }
1820  

Diff Legend

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