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 1803 by gezelter, Wed Oct 3 14:20:07 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 270 | Line 277 | namespace OpenMD {
277      // do some sanity checking
278  
279      int selectionCount = seleMan_.getSelectionCount();
280 +
281      int nIntegrable = info->getNGlobalIntegrableObjects();
282  
283      if (selectionCount > nIntegrable) {
# Line 288 | Line 296 | namespace OpenMD {
296        simError();
297      }
298  
299 +    areaAccumulator_ = new Accumulator();
300 +
301      nBins_ = rnemdParams->getOutputBins();
302  
303      data_.resize(RNEMD::ENDINDEX);
# Line 312 | Line 322 | namespace OpenMD {
322      outputMap_["TEMPERATURE"] =  TEMPERATURE;
323  
324      OutputData velocity;
325 <    velocity.units = "amu/fs";
325 >    velocity.units = "angstroms/fs";
326      velocity.title =  "Velocity";  
327      velocity.dataType = "Vector3d";
328      velocity.accumulator.reserve(nBins_);
# Line 381 | Line 391 | namespace OpenMD {
391      // dt = exchange time interval
392      // flux = target flux
393  
394 <    kineticTarget_ = 2.0*kineticFlux_*exchangeTime_*hmat(0,0)*hmat(1,1);
395 <    momentumTarget_ = 2.0*momentumFluxVector_*exchangeTime_*hmat(0,0)*hmat(1,1);
394 >    RealType area = currentSnap_->getXYarea();
395 >    kineticTarget_ = 2.0 * kineticFlux_ * exchangeTime_ * area;
396 >    momentumTarget_ = 2.0 * momentumFluxVector_ * exchangeTime_ * area;
397  
398      // total exchange sums are zeroed out at the beginning:
399  
# Line 407 | Line 418 | namespace OpenMD {
418    }
419    
420    RNEMD::~RNEMD() {
421 <    
421 >    if (!doRNEMD_) return;
422   #ifdef IS_MPI
423      if (worldRank == 0) {
424   #endif
# Line 429 | Line 440 | namespace OpenMD {
440    }
441  
442    void RNEMD::doSwap() {
443 <
443 >    if (!doRNEMD_) return;
444      Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
445      Mat3x3d hmat = currentSnap_->getHmat();
446  
# Line 488 | Line 499 | namespace OpenMD {
499                  + angMom[2]*angMom[2]/I(2, 2);
500              }
501            } //angular momenta exchange enabled
491 –          //energyConvert temporarily disabled
492 –          //make kineticExchange_ comparable between swap & scale
493 –          //value = value * 0.5 / PhysicalConstants::energyConvert;
502            value *= 0.5;
503            break;
504          case rnemdPx :
# Line 728 | Line 736 | namespace OpenMD {
736          
737          switch(rnemdFluxType_) {
738          case rnemdKE:
731 –          cerr << "KE\n";
739            kineticExchange_ += max_val - min_val;
740            break;
741          case rnemdPx:
# Line 741 | Line 748 | namespace OpenMD {
748            momentumExchange_.z() += max_val - min_val;
749            break;
750          default:
744 –          cerr << "default\n";
751            break;
752          }
753        } else {        
# 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 908 | Line 914 | namespace OpenMD {
914  
915        if ((c > 0.81) && (c < 1.21)) {//restrict scaling coefficients
916          c = sqrt(c);
917 <        //std::cerr << "cold slab scaling coefficient: " << c << endl;
912 <        //now convert to hotBin coefficient
917 >
918          RealType w = 0.0;
919          if (rnemdFluxType_ ==  rnemdFullKE) {
920            x = 1.0 + px * (1.0 - c);
# Line 947 | Line 952 | namespace OpenMD {
952              }
953            }
954            w = sqrt(w);
950 –          // std::cerr << "xh= " << x << "\tyh= " << y << "\tzh= " << z
951 –          //           << "\twh= " << w << endl;
955            for (sdi = hotBin.begin(); sdi != hotBin.end(); sdi++) {
956              if (rnemdFluxType_ == rnemdFullKE) {
957                vel = (*sdi)->getVel();
# Line 1213 | Line 1216 | namespace OpenMD {
1216    }
1217  
1218    void RNEMD::doVSS() {
1219 <
1219 >    if (!doRNEMD_) return;
1220      Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
1221      RealType time = currentSnap_->getTime();    
1222      Mat3x3d hmat = currentSnap_->getHmat();
# Line 1257 | Line 1260 | namespace OpenMD {
1260        
1261          if (inA) {
1262            hotBin.push_back(sd);
1260 –          //std::cerr << "before, velocity = " << vel << endl;
1263            Ph += mass * vel;
1262 –          //std::cerr << "after, velocity = " << vel << endl;
1264            Mh += mass;
1265            Kh += mass * vel.lengthSquare();
1266            if (rnemdFluxType_ == rnemdFullKE) {
# Line 1307 | Line 1308 | namespace OpenMD {
1308      
1309      Kh *= 0.5;
1310      Kc *= 0.5;
1310 –
1311 –    // std::cerr << "Mh= " << Mh << "\tKh= " << Kh << "\tMc= " << Mc
1312 –    //        << "\tKc= " << Kc << endl;
1313 –    // std::cerr << "Ph= " << Ph << "\tPc= " << Pc << endl;
1311      
1312   #ifdef IS_MPI
1313      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Ph[0], 3, MPI::REALTYPE, MPI::SUM);
# Line 1342 | Line 1339 | namespace OpenMD {
1339                if (hDenominator > 0.0) {
1340                  RealType h = sqrt(hNumerator / hDenominator);
1341                  if ((h > 0.9) && (h < 1.1)) {
1342 <                  // std::cerr << "cold slab scaling coefficient: " << c << "\n";
1346 <                  // std::cerr << "hot slab scaling coefficient: " << h <<  "\n";
1342 >
1343                    vector<StuntDouble*>::iterator sdi;
1344                    Vector3d vel;
1345                    for (sdi = coldBin.begin(); sdi != coldBin.end(); sdi++) {
# Line 1391 | Line 1387 | namespace OpenMD {
1387    }
1388  
1389    void RNEMD::doRNEMD() {
1390 <
1390 >    if (!doRNEMD_) return;
1391      trialCount_++;
1392      switch(rnemdMethod_) {
1393      case rnemdSwap:
# Line 1410 | Line 1406 | namespace OpenMD {
1406    }
1407  
1408    void RNEMD::collectData() {
1409 <
1409 >    if (!doRNEMD_) return;
1410      Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
1411      Mat3x3d hmat = currentSnap_->getHmat();
1412  
1413 +    areaAccumulator_->add(currentSnap_->getXYarea());
1414 +
1415      seleMan_.setSelectionSet(evaluator_.evaluate());
1416  
1417 <    int selei;
1417 >    int selei(0);
1418      StuntDouble* sd;
1419      int idx;
1420  
# Line 1442 | Line 1440 | namespace OpenMD {
1440               sd != NULL;
1441               sd = mol->nextIntegrableObject(iiter))
1442      */
1443 +
1444      for (sd = seleMan_.beginSelected(selei); sd != NULL;
1445           sd = seleMan_.nextSelected(selei)) {
1446 <      
1446 >    
1447        idx = sd->getLocalIndex();
1448        
1449        Vector3d pos = sd->getPos();
# Line 1454 | Line 1453 | namespace OpenMD {
1453        if (usePeriodicBoundaryConditions_)
1454          currentSnap_->wrapVector(pos);
1455  
1456 +
1457        // which bin is this stuntdouble in?
1458        // wrapped positions are in the range [-0.5*hmat(2,2), +0.5*hmat(2,2)]
1459        // Shift molecules by half a box to have bins start at 0
1460        // The modulo operator is used to wrap the case when we are
1461        // beyond the end of the bins back to the beginning.
1462        int binNo = int(nBins_ * (pos.z() / hmat(2,2) + 0.5)) % nBins_;
1463 <    
1463 >
1464        RealType mass = sd->getMass();
1465        Vector3d vel = sd->getVel();
1466  
# Line 1491 | Line 1491 | namespace OpenMD {
1491        }
1492      }
1493      
1494 –
1494   #ifdef IS_MPI
1495      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binCount[0],
1496                                nBins_, MPI::INT, MPI::SUM);
# Line 1518 | Line 1517 | namespace OpenMD {
1517        vel.x() = binPx[i] / binMass[i];
1518        vel.y() = binPy[i] / binMass[i];
1519        vel.z() = binPz[i] / binMass[i];
1520 <      den = binCount[i] * nBins_ / (hmat(0,0) * hmat(1,1) * hmat(2,2));
1520 >
1521 >      den = binMass[i] * nBins_ * PhysicalConstants::densityConvert
1522 >        / currentSnap_->getVolume() ;
1523 >
1524        temp = 2.0 * binKE[i] / (binDOF[i] * PhysicalConstants::kb *
1525                                 PhysicalConstants::energyConvert);
1526 <
1526 >  
1527        for (unsigned int j = 0; j < outputMask_.size(); ++j) {
1528          if(outputMask_[j]) {
1529            switch(j) {
# Line 1544 | Line 1546 | namespace OpenMD {
1546    }
1547  
1548    void RNEMD::getStarted() {
1549 +    if (!doRNEMD_) return;
1550      collectData();
1551      writeOutputFile();
1552    }
1553  
1554    void RNEMD::parseOutputFileFormat(const std::string& format) {
1555 +    if (!doRNEMD_) return;
1556      StringTokenizer tokenizer(format, " ,;|\t\n\r");
1557      
1558      while(tokenizer.hasMoreTokens()) {
# Line 1569 | Line 1573 | namespace OpenMD {
1573    }
1574    
1575    void RNEMD::writeOutputFile() {
1576 +    if (!doRNEMD_) return;
1577      
1578   #ifdef IS_MPI
1579      // If we're the root node, should we print out the results
# Line 1588 | Line 1593 | namespace OpenMD {
1593        Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
1594  
1595        RealType time = currentSnap_->getTime();
1596 <      
1597 <      
1596 >      RealType avgArea;
1597 >      areaAccumulator_->getAverage(avgArea);
1598 >      RealType Jz = kineticExchange_ / (2.0 * time * avgArea)
1599 >        / PhysicalConstants::energyConvert;
1600 >      Vector3d JzP = momentumExchange_ / (2.0 * time * avgArea);      
1601 >
1602        rnemdFile_ << "#######################################################\n";
1603        rnemdFile_ << "# RNEMD {\n";
1604  
1605        map<string, RNEMDMethod>::iterator mi;
1606        for(mi = stringToMethod_.begin(); mi != stringToMethod_.end(); ++mi) {
1607          if ( (*mi).second == rnemdMethod_)
1608 <          rnemdFile_ << "#    exchangeMethod  = " << (*mi).first << "\n";
1608 >          rnemdFile_ << "#    exchangeMethod  = \"" << (*mi).first << "\";\n";
1609        }
1610        map<string, RNEMDFluxType>::iterator fi;
1611        for(fi = stringToFluxType_.begin(); fi != stringToFluxType_.end(); ++fi) {
1612          if ( (*fi).second == rnemdFluxType_)
1613 <          rnemdFile_ << "#    fluxType  = " << (*fi).first << "\n";
1613 >          rnemdFile_ << "#    fluxType  = \"" << (*fi).first << "\";\n";
1614        }
1615        
1616 <      rnemdFile_ << "#    exchangeTime = " << exchangeTime_ << " fs\n";
1616 >      rnemdFile_ << "#    exchangeTime = " << exchangeTime_ << ";\n";
1617  
1618        rnemdFile_ << "#    objectSelection = \""
1619 <                 << rnemdObjectSelection_ << "\"\n";
1620 <      rnemdFile_ << "#    slabWidth = " << slabWidth_ << " angstroms\n";
1621 <      rnemdFile_ << "#    slabAcenter = " << slabACenter_ << " angstroms\n";
1622 <      rnemdFile_ << "#    slabBcenter = " << slabBCenter_ << " angstroms\n";
1619 >                 << rnemdObjectSelection_ << "\";\n";
1620 >      rnemdFile_ << "#    slabWidth = " << slabWidth_ << ";\n";
1621 >      rnemdFile_ << "#    slabAcenter = " << slabACenter_ << ";\n";
1622 >      rnemdFile_ << "#    slabBcenter = " << slabBCenter_ << ";\n";
1623        rnemdFile_ << "# }\n";
1624        rnemdFile_ << "#######################################################\n";
1625 <      
1626 <      rnemdFile_ << "# running time = " << time << " fs\n";
1627 <      rnemdFile_ << "# target kinetic flux = " << kineticFlux_ << "\n";
1628 <      rnemdFile_ << "# target momentum flux = " << momentumFluxVector_ << "\n";
1629 <      
1630 <      rnemdFile_ << "# target one-time kinetic exchange = " << kineticTarget_
1631 <                 << "\n";
1632 <      rnemdFile_ << "# target one-time momentum exchange = " << momentumTarget_
1633 <                 << "\n";
1634 <      
1635 <      rnemdFile_ << "# actual kinetic exchange = " << kineticExchange_ << "\n";
1636 <      rnemdFile_ << "# actual momentum exchange = " << momentumExchange_
1637 <                 << "\n";
1638 <      
1639 <      rnemdFile_ << "# attempted exchanges: " << trialCount_ << "\n";
1640 <      rnemdFile_ << "# failed exchanges: " << failTrialCount_ << "\n";
1641 <
1642 <      
1625 >      rnemdFile_ << "# RNEMD report:\n";      
1626 >      rnemdFile_ << "#     running time = " << time << " fs\n";
1627 >      rnemdFile_ << "#     target flux:\n";
1628 >      rnemdFile_ << "#         kinetic = "
1629 >                 << kineticFlux_ / PhysicalConstants::energyConvert
1630 >                 << " (kcal/mol/A^2/fs)\n";
1631 >      rnemdFile_ << "#         momentum = " << momentumFluxVector_
1632 >                 << " (amu/A/fs^2)\n";
1633 >      rnemdFile_ << "#     target one-time exchanges:\n";
1634 >      rnemdFile_ << "#         kinetic = "
1635 >                 << kineticTarget_ / PhysicalConstants::energyConvert
1636 >                 << " (kcal/mol)\n";
1637 >      rnemdFile_ << "#         momentum = " << momentumTarget_
1638 >                 << " (amu*A/fs)\n";
1639 >      rnemdFile_ << "#     actual exchange totals:\n";
1640 >      rnemdFile_ << "#         kinetic = "
1641 >                 << kineticExchange_ / PhysicalConstants::energyConvert
1642 >                 << " (kcal/mol)\n";
1643 >      rnemdFile_ << "#         momentum = " << momentumExchange_
1644 >                 << " (amu*A/fs)\n";      
1645 >      rnemdFile_ << "#     actual flux:\n";
1646 >      rnemdFile_ << "#         kinetic = " << Jz
1647 >                 << " (kcal/mol/A^2/fs)\n";
1648 >      rnemdFile_ << "#         momentum = " << JzP
1649 >                 << " (amu/A/fs^2)\n";
1650 >      rnemdFile_ << "#     exchange statistics:\n";
1651 >      rnemdFile_ << "#         attempted = " << trialCount_ << "\n";
1652 >      rnemdFile_ << "#         failed = " << failTrialCount_ << "\n";    
1653        if (rnemdMethod_ == rnemdNIVS) {
1654 <        rnemdFile_ << "# NIVS root-check warnings: " << failRootCount_ << "\n";
1654 >        rnemdFile_ << "#         NIVS root-check errors = "
1655 >                   << failRootCount_ << "\n";
1656        }
1637 –
1657        rnemdFile_ << "#######################################################\n";
1658        
1659        
# Line 1645 | Line 1664 | namespace OpenMD {
1664          if (outputMask_[i]) {
1665            rnemdFile_ << "\t" << data_[i].title <<
1666              "(" << data_[i].units << ")";
1667 +          // add some extra tabs for column alignment
1668 +          if (data_[i].dataType == "Vector3d") rnemdFile_ << "\t\t";
1669          }
1670        }
1671        rnemdFile_ << std::endl;
# Line 1671 | Line 1692 | namespace OpenMD {
1692          rnemdFile_ << std::endl;
1693          
1694        }        
1695 +
1696 +      rnemdFile_ << "#######################################################\n";
1697 +      rnemdFile_ << "# Standard Deviations in those quantities follow:\n";
1698 +      rnemdFile_ << "#######################################################\n";
1699 +
1700 +
1701 +      for (unsigned int j = 0; j < nBins_; j++) {        
1702 +        rnemdFile_ << "#";
1703 +        for (unsigned int i = 0; i < outputMask_.size(); ++i) {
1704 +          if (outputMask_[i]) {
1705 +            if (data_[i].dataType == "RealType")
1706 +              writeRealStdDev(i,j);
1707 +            else if (data_[i].dataType == "Vector3d")
1708 +              writeVectorStdDev(i,j);
1709 +            else {
1710 +              sprintf( painCave.errMsg,
1711 +                       "RNEMD found an unknown data type for: %s ",
1712 +                       data_[i].title.c_str());
1713 +              painCave.isFatal = 1;
1714 +              simError();
1715 +            }
1716 +          }
1717 +        }
1718 +        rnemdFile_ << std::endl;
1719 +        
1720 +      }        
1721        
1722        rnemdFile_.flush();
1723        rnemdFile_.close();
# Line 1682 | Line 1729 | namespace OpenMD {
1729    }
1730    
1731    void RNEMD::writeReal(int index, unsigned int bin) {
1732 +    if (!doRNEMD_) return;
1733      assert(index >=0 && index < ENDINDEX);
1734 <    assert(bin >=0 && bin < nBins_);
1734 >    assert(bin < nBins_);
1735      RealType s;
1736      
1737      data_[index].accumulator[bin]->getAverage(s);
# Line 1700 | Line 1748 | namespace OpenMD {
1748    }
1749    
1750    void RNEMD::writeVector(int index, unsigned int bin) {
1751 +    if (!doRNEMD_) return;
1752      assert(index >=0 && index < ENDINDEX);
1753 <    assert(bin >=0 && bin < nBins_);
1753 >    assert(bin < nBins_);
1754      Vector3d s;
1755      dynamic_cast<VectorAccumulator*>(data_[index].accumulator[bin])->getAverage(s);
1756      if (isinf(s[0]) || isnan(s[0]) ||
# Line 1716 | Line 1765 | namespace OpenMD {
1765        rnemdFile_ << "\t" << s[0] << "\t" << s[1] << "\t" << s[2];
1766      }
1767    }  
1768 +
1769 +  void RNEMD::writeRealStdDev(int index, unsigned int bin) {
1770 +    if (!doRNEMD_) return;
1771 +    assert(index >=0 && index < ENDINDEX);
1772 +    assert(bin < nBins_);
1773 +    RealType s;
1774 +    
1775 +    data_[index].accumulator[bin]->getStdDev(s);
1776 +    
1777 +    if (! isinf(s) && ! isnan(s)) {
1778 +      rnemdFile_ << "\t" << s;
1779 +    } else{
1780 +      sprintf( painCave.errMsg,
1781 +               "RNEMD detected a numerical error writing: %s std. dev. for bin %d",
1782 +               data_[index].title.c_str(), bin);
1783 +      painCave.isFatal = 1;
1784 +      simError();
1785 +    }    
1786 +  }
1787 +  
1788 +  void RNEMD::writeVectorStdDev(int index, unsigned int bin) {
1789 +    if (!doRNEMD_) return;
1790 +    assert(index >=0 && index < ENDINDEX);
1791 +    assert(bin < nBins_);
1792 +    Vector3d s;
1793 +    dynamic_cast<VectorAccumulator*>(data_[index].accumulator[bin])->getStdDev(s);
1794 +    if (isinf(s[0]) || isnan(s[0]) ||
1795 +        isinf(s[1]) || isnan(s[1]) ||
1796 +        isinf(s[2]) || isnan(s[2]) ) {      
1797 +      sprintf( painCave.errMsg,
1798 +               "RNEMD detected a numerical error writing: %s std. dev. for bin %d",
1799 +               data_[index].title.c_str(), bin);
1800 +      painCave.isFatal = 1;
1801 +      simError();
1802 +    } else {
1803 +      rnemdFile_ << "\t" << s[0] << "\t" << s[1] << "\t" << s[2];
1804 +    }
1805 +  }  
1806   }
1807  

Diff Legend

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