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 1773 by gezelter, Tue Aug 7 18:26:40 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 288 | Line 295 | namespace OpenMD {
295        simError();
296      }
297  
298 +    areaAccumulator_ = new Accumulator();
299 +
300      nBins_ = rnemdParams->getOutputBins();
301  
302      data_.resize(RNEMD::ENDINDEX);
# Line 312 | 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 381 | Line 390 | namespace OpenMD {
390      // dt = exchange time interval
391      // flux = target flux
392  
393 <    kineticTarget_ = 2.0*kineticFlux_*exchangeTime_*hmat(0,0)*hmat(1,1);
394 <    momentumTarget_ = 2.0*momentumFluxVector_*exchangeTime_*hmat(0,0)*hmat(1,1);
393 >    RealType area = currentSnap_->getXYarea();
394 >    kineticTarget_ = 2.0 * kineticFlux_ * exchangeTime_ * area;
395 >    momentumTarget_ = 2.0 * momentumFluxVector_ * exchangeTime_ * area;
396  
397      // total exchange sums are zeroed out at the beginning:
398  
# Line 407 | 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 429 | 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 488 | Line 498 | namespace OpenMD {
498                  + angMom[2]*angMom[2]/I(2, 2);
499              }
500            } //angular momenta exchange enabled
491 –          //energyConvert temporarily disabled
492 –          //make kineticExchange_ comparable between swap & scale
493 –          //value = value * 0.5 / PhysicalConstants::energyConvert;
501            value *= 0.5;
502            break;
503          case rnemdPx :
# Line 728 | Line 735 | namespace OpenMD {
735          
736          switch(rnemdFluxType_) {
737          case rnemdKE:
731 –          cerr << "KE\n";
738            kineticExchange_ += max_val - min_val;
739            break;
740          case rnemdPx:
# Line 741 | Line 747 | namespace OpenMD {
747            momentumExchange_.z() += max_val - min_val;
748            break;
749          default:
744 –          cerr << "default\n";
750            break;
751          }
752        } else {        
# Line 764 | 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 1213 | 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 1391 | Line 1396 | namespace OpenMD {
1396    }
1397  
1398    void RNEMD::doRNEMD() {
1399 <
1399 >    if (!doRNEMD_) return;
1400      trialCount_++;
1401      switch(rnemdMethod_) {
1402      case rnemdSwap:
# Line 1410 | 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  
1422 +    areaAccumulator_->add(currentSnap_->getXYarea());
1423 +
1424      seleMan_.setSelectionSet(evaluator_.evaluate());
1425  
1426      int selei;
# Line 1454 | Line 1461 | namespace OpenMD {
1461        if (usePeriodicBoundaryConditions_)
1462          currentSnap_->wrapVector(pos);
1463  
1464 +
1465        // which bin is this stuntdouble in?
1466        // wrapped positions are in the range [-0.5*hmat(2,2), +0.5*hmat(2,2)]
1467        // Shift molecules by half a box to have bins start at 0
# Line 1518 | 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_ / (hmat(0,0) * hmat(1,1) * hmat(2,2));
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 1544 | 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 1569 | 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 1588 | Line 1602 | namespace OpenMD {
1602        Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
1603  
1604        RealType time = currentSnap_->getTime();
1605 <      
1606 <      
1605 >      RealType avgArea;
1606 >      areaAccumulator_->getAverage(avgArea);
1607 >      RealType Jz = kineticExchange_ / (2.0 * time * avgArea)
1608 >        / PhysicalConstants::energyConvert;
1609 >      Vector3d JzP = momentumExchange_ / (2.0 * time * avgArea);      
1610 >
1611        rnemdFile_ << "#######################################################\n";
1612        rnemdFile_ << "# RNEMD {\n";
1613  
1614        map<string, RNEMDMethod>::iterator mi;
1615        for(mi = stringToMethod_.begin(); mi != stringToMethod_.end(); ++mi) {
1616          if ( (*mi).second == rnemdMethod_)
1617 <          rnemdFile_ << "#    exchangeMethod  = " << (*mi).first << "\n";
1617 >          rnemdFile_ << "#    exchangeMethod  = \"" << (*mi).first << "\";\n";
1618        }
1619        map<string, RNEMDFluxType>::iterator fi;
1620        for(fi = stringToFluxType_.begin(); fi != stringToFluxType_.end(); ++fi) {
1621          if ( (*fi).second == rnemdFluxType_)
1622 <          rnemdFile_ << "#    fluxType  = " << (*fi).first << "\n";
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";
1629 <      rnemdFile_ << "#    slabWidth = " << slabWidth_ << " angstroms\n";
1630 <      rnemdFile_ << "#    slabAcenter = " << slabACenter_ << " angstroms\n";
1631 <      rnemdFile_ << "#    slabBcenter = " << slabBCenter_ << " angstroms\n";
1628 >                 << rnemdObjectSelection_ << "\";\n";
1629 >      rnemdFile_ << "#    slabWidth = " << slabWidth_ << ";\n";
1630 >      rnemdFile_ << "#    slabAcenter = " << slabACenter_ << ";\n";
1631 >      rnemdFile_ << "#    slabBcenter = " << slabBCenter_ << ";\n";
1632        rnemdFile_ << "# }\n";
1633        rnemdFile_ << "#######################################################\n";
1634 <      
1635 <      rnemdFile_ << "# running time = " << time << " fs\n";
1636 <      rnemdFile_ << "# target kinetic flux = " << kineticFlux_ << "\n";
1637 <      rnemdFile_ << "# target momentum flux = " << momentumFluxVector_ << "\n";
1638 <      
1639 <      rnemdFile_ << "# target one-time kinetic exchange = " << kineticTarget_
1640 <                 << "\n";
1641 <      rnemdFile_ << "# target one-time momentum exchange = " << momentumTarget_
1642 <                 << "\n";
1643 <      
1644 <      rnemdFile_ << "# actual kinetic exchange = " << kineticExchange_ << "\n";
1645 <      rnemdFile_ << "# actual momentum exchange = " << momentumExchange_
1646 <                 << "\n";
1647 <      
1648 <      rnemdFile_ << "# attempted exchanges: " << trialCount_ << "\n";
1649 <      rnemdFile_ << "# failed exchanges: " << failTrialCount_ << "\n";
1650 <
1651 <      
1634 >      rnemdFile_ << "# RNEMD report:\n";      
1635 >      rnemdFile_ << "#     running time = " << time << " fs\n";
1636 >      rnemdFile_ << "#     target flux:\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 = "
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 = "
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
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";    
1662        if (rnemdMethod_ == rnemdNIVS) {
1663 <        rnemdFile_ << "# NIVS root-check warnings: " << failRootCount_ << "\n";
1663 >        rnemdFile_ << "#         NIVS root-check errors = "
1664 >                   << failRootCount_ << "\n";
1665        }
1637 –
1666        rnemdFile_ << "#######################################################\n";
1667        
1668        
# Line 1645 | 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 1671 | Line 1701 | namespace OpenMD {
1701          rnemdFile_ << std::endl;
1702          
1703        }        
1704 +
1705 +      rnemdFile_ << "#######################################################\n";
1706 +      rnemdFile_ << "# Standard Deviations in those quantities follow:\n";
1707 +      rnemdFile_ << "#######################################################\n";
1708 +
1709 +
1710 +      for (unsigned int j = 0; j < nBins_; j++) {        
1711 +        rnemdFile_ << "#";
1712 +        for (unsigned int i = 0; i < outputMask_.size(); ++i) {
1713 +          if (outputMask_[i]) {
1714 +            if (data_[i].dataType == "RealType")
1715 +              writeRealStdDev(i,j);
1716 +            else if (data_[i].dataType == "Vector3d")
1717 +              writeVectorStdDev(i,j);
1718 +            else {
1719 +              sprintf( painCave.errMsg,
1720 +                       "RNEMD found an unknown data type for: %s ",
1721 +                       data_[i].title.c_str());
1722 +              painCave.isFatal = 1;
1723 +              simError();
1724 +            }
1725 +          }
1726 +        }
1727 +        rnemdFile_ << std::endl;
1728 +        
1729 +      }        
1730        
1731        rnemdFile_.flush();
1732        rnemdFile_.close();
# Line 1682 | 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 >=0 && bin < nBins_);
1743 >    assert(bin < nBins_);
1744      RealType s;
1745      
1746      data_[index].accumulator[bin]->getAverage(s);
# Line 1700 | 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 >=0 && bin < nBins_);
1762 >    assert(bin < nBins_);
1763      Vector3d s;
1764      dynamic_cast<VectorAccumulator*>(data_[index].accumulator[bin])->getAverage(s);
1765      if (isinf(s[0]) || isnan(s[0]) ||
# Line 1716 | Line 1774 | namespace OpenMD {
1774        rnemdFile_ << "\t" << s[0] << "\t" << s[1] << "\t" << s[2];
1775      }
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;
1783 +    
1784 +    data_[index].accumulator[bin]->getStdDev(s);
1785 +    
1786 +    if (! isinf(s) && ! isnan(s)) {
1787 +      rnemdFile_ << "\t" << s;
1788 +    } else{
1789 +      sprintf( painCave.errMsg,
1790 +               "RNEMD detected a numerical error writing: %s std. dev. for bin %d",
1791 +               data_[index].title.c_str(), bin);
1792 +      painCave.isFatal = 1;
1793 +      simError();
1794 +    }    
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;
1802 +    dynamic_cast<VectorAccumulator*>(data_[index].accumulator[bin])->getStdDev(s);
1803 +    if (isinf(s[0]) || isnan(s[0]) ||
1804 +        isinf(s[1]) || isnan(s[1]) ||
1805 +        isinf(s[2]) || isnan(s[2]) ) {      
1806 +      sprintf( painCave.errMsg,
1807 +               "RNEMD detected a numerical error writing: %s std. dev. for bin %d",
1808 +               data_[index].title.c_str(), bin);
1809 +      painCave.isFatal = 1;
1810 +      simError();
1811 +    } else {
1812 +      rnemdFile_ << "\t" << s[0] << "\t" << s[1] << "\t" << s[2];
1813 +    }
1814 +  }  
1815   }
1816  

Diff Legend

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