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.
Revision 1776 by gezelter, Thu Aug 9 15:52:59 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 287 | Line 290 | namespace OpenMD {
290        painCave.severity = OPENMD_WARNING;
291        simError();
292      }
293 +
294 +    areaAccumulator_ = new Accumulator();
295  
296      nBins_ = rnemdParams->getOutputBins();
297  
# Line 381 | Line 386 | namespace OpenMD {
386      // dt = exchange time interval
387      // flux = target flux
388  
389 <    kineticTarget_ = 2.0*kineticFlux_*exchangeTime_*hmat(0,0)*hmat(1,1);
390 <    momentumTarget_ = 2.0*momentumFluxVector_*exchangeTime_*hmat(0,0)*hmat(1,1);
389 >    RealType area = currentSnap_->getXYarea();
390 >    kineticTarget_ = 2.0 * kineticFlux_ * exchangeTime_ * area;
391 >    momentumTarget_ = 2.0 * momentumFluxVector_ * exchangeTime_ * area;
392  
393      // total exchange sums are zeroed out at the beginning:
394  
# Line 407 | Line 413 | namespace OpenMD {
413    }
414    
415    RNEMD::~RNEMD() {
416 <    
416 >    if (!doRNEMD_) return;
417   #ifdef IS_MPI
418      if (worldRank == 0) {
419   #endif
# Line 429 | Line 435 | namespace OpenMD {
435    }
436  
437    void RNEMD::doSwap() {
438 <
438 >    if (!doRNEMD_) return;
439      Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
440      Mat3x3d hmat = currentSnap_->getHmat();
441  
# Line 764 | Line 770 | namespace OpenMD {
770    }
771    
772    void RNEMD::doNIVS() {
773 <
773 >    if (!doRNEMD_) return;
774      Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
775      Mat3x3d hmat = currentSnap_->getHmat();
776  
# Line 1213 | Line 1219 | namespace OpenMD {
1219    }
1220  
1221    void RNEMD::doVSS() {
1222 <
1222 >    if (!doRNEMD_) return;
1223      Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
1224      RealType time = currentSnap_->getTime();    
1225      Mat3x3d hmat = currentSnap_->getHmat();
# Line 1391 | Line 1397 | namespace OpenMD {
1397    }
1398  
1399    void RNEMD::doRNEMD() {
1400 <
1400 >    if (!doRNEMD_) return;
1401      trialCount_++;
1402      switch(rnemdMethod_) {
1403      case rnemdSwap:
# Line 1410 | Line 1416 | namespace OpenMD {
1416    }
1417  
1418    void RNEMD::collectData() {
1419 <
1419 >    if (!doRNEMD_) return;
1420      Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
1421      Mat3x3d hmat = currentSnap_->getHmat();
1422  
1423 +    areaAccumulator_->add(currentSnap_->getXYarea());
1424 +
1425      seleMan_.setSelectionSet(evaluator_.evaluate());
1426  
1427      int selei;
# Line 1454 | Line 1462 | namespace OpenMD {
1462        if (usePeriodicBoundaryConditions_)
1463          currentSnap_->wrapVector(pos);
1464  
1465 +
1466        // which bin is this stuntdouble in?
1467        // wrapped positions are in the range [-0.5*hmat(2,2), +0.5*hmat(2,2)]
1468        // Shift molecules by half a box to have bins start at 0
# Line 1518 | Line 1527 | namespace OpenMD {
1527        vel.x() = binPx[i] / binMass[i];
1528        vel.y() = binPy[i] / binMass[i];
1529        vel.z() = binPz[i] / binMass[i];
1530 <      den = binCount[i] * nBins_ / (hmat(0,0) * hmat(1,1) * hmat(2,2));
1530 >      den = binCount[i] * nBins_ / currentSnap_->getVolume();
1531        temp = 2.0 * binKE[i] / (binDOF[i] * PhysicalConstants::kb *
1532                                 PhysicalConstants::energyConvert);
1533  
# Line 1544 | Line 1553 | namespace OpenMD {
1553    }
1554  
1555    void RNEMD::getStarted() {
1556 +    if (!doRNEMD_) return;
1557      collectData();
1558      writeOutputFile();
1559    }
1560  
1561    void RNEMD::parseOutputFileFormat(const std::string& format) {
1562 +    if (!doRNEMD_) return;
1563      StringTokenizer tokenizer(format, " ,;|\t\n\r");
1564      
1565      while(tokenizer.hasMoreTokens()) {
# Line 1569 | Line 1580 | namespace OpenMD {
1580    }
1581    
1582    void RNEMD::writeOutputFile() {
1583 +    if (!doRNEMD_) return;
1584      
1585   #ifdef IS_MPI
1586      // If we're the root node, should we print out the results
# Line 1588 | Line 1600 | namespace OpenMD {
1600        Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
1601  
1602        RealType time = currentSnap_->getTime();
1603 <      
1604 <      
1603 >      RealType avgArea;
1604 >      areaAccumulator_->getAverage(avgArea);
1605 >      RealType Jz = kineticExchange_ / (2.0 * time * avgArea);
1606 >      Vector3d JzP = momentumExchange_ / (2.0 * time * avgArea);      
1607 >
1608        rnemdFile_ << "#######################################################\n";
1609        rnemdFile_ << "# RNEMD {\n";
1610  
1611        map<string, RNEMDMethod>::iterator mi;
1612        for(mi = stringToMethod_.begin(); mi != stringToMethod_.end(); ++mi) {
1613          if ( (*mi).second == rnemdMethod_)
1614 <          rnemdFile_ << "#    exchangeMethod  = " << (*mi).first << "\n";
1614 >          rnemdFile_ << "#    exchangeMethod  = \"" << (*mi).first << "\";\n";
1615        }
1616        map<string, RNEMDFluxType>::iterator fi;
1617        for(fi = stringToFluxType_.begin(); fi != stringToFluxType_.end(); ++fi) {
1618          if ( (*fi).second == rnemdFluxType_)
1619 <          rnemdFile_ << "#    fluxType  = " << (*fi).first << "\n";
1619 >          rnemdFile_ << "#    fluxType  = \"" << (*fi).first << "\";\n";
1620        }
1621        
1622 <      rnemdFile_ << "#    exchangeTime = " << exchangeTime_ << " fs\n";
1622 >      rnemdFile_ << "#    exchangeTime = " << exchangeTime_ << ";\n";
1623  
1624        rnemdFile_ << "#    objectSelection = \""
1625 <                 << rnemdObjectSelection_ << "\"\n";
1626 <      rnemdFile_ << "#    slabWidth = " << slabWidth_ << " angstroms\n";
1627 <      rnemdFile_ << "#    slabAcenter = " << slabACenter_ << " angstroms\n";
1628 <      rnemdFile_ << "#    slabBcenter = " << slabBCenter_ << " angstroms\n";
1625 >                 << rnemdObjectSelection_ << "\";\n";
1626 >      rnemdFile_ << "#    slabWidth = " << slabWidth_ << ";\n";
1627 >      rnemdFile_ << "#    slabAcenter = " << slabACenter_ << ";\n";
1628 >      rnemdFile_ << "#    slabBcenter = " << slabBCenter_ << ";\n";
1629        rnemdFile_ << "# }\n";
1630        rnemdFile_ << "#######################################################\n";
1631 <      
1632 <      rnemdFile_ << "# running time = " << time << " fs\n";
1633 <      rnemdFile_ << "# target kinetic flux = " << kineticFlux_ << "\n";
1634 <      rnemdFile_ << "# target momentum flux = " << momentumFluxVector_ << "\n";
1635 <      
1636 <      rnemdFile_ << "# target one-time kinetic exchange = " << kineticTarget_
1637 <                 << "\n";
1638 <      rnemdFile_ << "# target one-time momentum exchange = " << momentumTarget_
1639 <                 << "\n";
1640 <      
1641 <      rnemdFile_ << "# actual kinetic exchange = " << kineticExchange_ << "\n";
1642 <      rnemdFile_ << "# actual momentum exchange = " << momentumExchange_
1643 <                 << "\n";
1644 <      
1645 <      rnemdFile_ << "# attempted exchanges: " << trialCount_ << "\n";
1646 <      rnemdFile_ << "# failed exchanges: " << failTrialCount_ << "\n";
1647 <
1633 <      
1631 >      rnemdFile_ << "# RNEMD report:\n";      
1632 >      rnemdFile_ << "#     running time = " << time << " fs\n";
1633 >      rnemdFile_ << "#     target flux:\n";
1634 >      rnemdFile_ << "#         kinetic = " << kineticFlux_ << "\n";
1635 >      rnemdFile_ << "#         momentum = " << momentumFluxVector_ << "\n";
1636 >      rnemdFile_ << "#     target one-time exchanges:\n";
1637 >      rnemdFile_ << "#         kinetic = " << kineticTarget_ << "\n";
1638 >      rnemdFile_ << "#         momentum = " << momentumTarget_ << "\n";
1639 >      rnemdFile_ << "#     actual exchange totals:\n";
1640 >      rnemdFile_ << "#         kinetic = " << kineticExchange_ << "\n";
1641 >      rnemdFile_ << "#         momentum = " << momentumExchange_  << "\n";
1642 >      rnemdFile_ << "#     actual flux:\n";
1643 >      rnemdFile_ << "#         kinetic = " << Jz << "\n";
1644 >      rnemdFile_ << "#         momentum = " << JzP  << "\n";
1645 >      rnemdFile_ << "#     exchange statistics:\n";
1646 >      rnemdFile_ << "#         attempted = " << trialCount_ << "\n";
1647 >      rnemdFile_ << "#         failed = " << failTrialCount_ << "\n";    
1648        if (rnemdMethod_ == rnemdNIVS) {
1649 <        rnemdFile_ << "# NIVS root-check warnings: " << failRootCount_ << "\n";
1649 >        rnemdFile_ << "#         NIVS root-check errors = "
1650 >                   << failRootCount_ << "\n";
1651        }
1637 –
1652        rnemdFile_ << "#######################################################\n";
1653        
1654        
# Line 1671 | Line 1685 | namespace OpenMD {
1685          rnemdFile_ << std::endl;
1686          
1687        }        
1688 +
1689 +      rnemdFile_ << "#######################################################\n";
1690 +      rnemdFile_ << "# Standard Deviations in those quantities follow:\n";
1691 +      rnemdFile_ << "#######################################################\n";
1692 +
1693 +
1694 +      for (unsigned int j = 0; j < nBins_; j++) {        
1695 +        rnemdFile_ << "#";
1696 +        for (unsigned int i = 0; i < outputMask_.size(); ++i) {
1697 +          if (outputMask_[i]) {
1698 +            if (data_[i].dataType == "RealType")
1699 +              writeRealStdDev(i,j);
1700 +            else if (data_[i].dataType == "Vector3d")
1701 +              writeVectorStdDev(i,j);
1702 +            else {
1703 +              sprintf( painCave.errMsg,
1704 +                       "RNEMD found an unknown data type for: %s ",
1705 +                       data_[i].title.c_str());
1706 +              painCave.isFatal = 1;
1707 +              simError();
1708 +            }
1709 +          }
1710 +        }
1711 +        rnemdFile_ << std::endl;
1712 +        
1713 +      }        
1714        
1715        rnemdFile_.flush();
1716        rnemdFile_.close();
# Line 1682 | Line 1722 | namespace OpenMD {
1722    }
1723    
1724    void RNEMD::writeReal(int index, unsigned int bin) {
1725 +    if (!doRNEMD_) return;
1726      assert(index >=0 && index < ENDINDEX);
1727 <    assert(bin >=0 && bin < nBins_);
1727 >    assert(bin < nBins_);
1728      RealType s;
1729      
1730      data_[index].accumulator[bin]->getAverage(s);
# Line 1700 | Line 1741 | namespace OpenMD {
1741    }
1742    
1743    void RNEMD::writeVector(int index, unsigned int bin) {
1744 +    if (!doRNEMD_) return;
1745      assert(index >=0 && index < ENDINDEX);
1746 <    assert(bin >=0 && bin < nBins_);
1746 >    assert(bin < nBins_);
1747      Vector3d s;
1748      dynamic_cast<VectorAccumulator*>(data_[index].accumulator[bin])->getAverage(s);
1749      if (isinf(s[0]) || isnan(s[0]) ||
# Line 1716 | Line 1758 | namespace OpenMD {
1758        rnemdFile_ << "\t" << s[0] << "\t" << s[1] << "\t" << s[2];
1759      }
1760    }  
1761 +
1762 +  void RNEMD::writeRealStdDev(int index, unsigned int bin) {
1763 +    if (!doRNEMD_) return;
1764 +    assert(index >=0 && index < ENDINDEX);
1765 +    assert(bin < nBins_);
1766 +    RealType s;
1767 +    
1768 +    data_[index].accumulator[bin]->getStdDev(s);
1769 +    
1770 +    if (! isinf(s) && ! isnan(s)) {
1771 +      rnemdFile_ << "\t" << s;
1772 +    } else{
1773 +      sprintf( painCave.errMsg,
1774 +               "RNEMD detected a numerical error writing: %s std. dev. for bin %d",
1775 +               data_[index].title.c_str(), bin);
1776 +      painCave.isFatal = 1;
1777 +      simError();
1778 +    }    
1779 +  }
1780 +  
1781 +  void RNEMD::writeVectorStdDev(int index, unsigned int bin) {
1782 +    if (!doRNEMD_) return;
1783 +    assert(index >=0 && index < ENDINDEX);
1784 +    assert(bin < nBins_);
1785 +    Vector3d s;
1786 +    dynamic_cast<VectorAccumulator*>(data_[index].accumulator[bin])->getStdDev(s);
1787 +    if (isinf(s[0]) || isnan(s[0]) ||
1788 +        isinf(s[1]) || isnan(s[1]) ||
1789 +        isinf(s[2]) || isnan(s[2]) ) {      
1790 +      sprintf( painCave.errMsg,
1791 +               "RNEMD detected a numerical error writing: %s std. dev. for bin %d",
1792 +               data_[index].title.c_str(), bin);
1793 +      painCave.isFatal = 1;
1794 +      simError();
1795 +    } else {
1796 +      rnemdFile_ << "\t" << s[0] << "\t" << s[1] << "\t" << s[2];
1797 +    }
1798 +  }  
1799   }
1800  

Diff Legend

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