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 1775 by gezelter, Wed Aug 8 18:45:52 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 314 | 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 410 | 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 432 | 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 491 | Line 499 | namespace OpenMD {
499                  + angMom[2]*angMom[2]/I(2, 2);
500              }
501            } //angular momenta exchange enabled
494 –          //energyConvert temporarily disabled
495 –          //make kineticExchange_ comparable between swap & scale
496 –          //value = value * 0.5 / PhysicalConstants::energyConvert;
502            value *= 0.5;
503            break;
504          case rnemdPx :
# Line 731 | Line 736 | namespace OpenMD {
736          
737          switch(rnemdFluxType_) {
738          case rnemdKE:
734 –          cerr << "KE\n";
739            kineticExchange_ += max_val - min_val;
740            break;
741          case rnemdPx:
# Line 744 | Line 748 | namespace OpenMD {
748            momentumExchange_.z() += max_val - min_val;
749            break;
750          default:
747 –          cerr << "default\n";
751            break;
752          }
753        } else {        
# Line 767 | 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 911 | 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;
915 <        //now convert to hotBin coefficient
917 >
918          RealType w = 0.0;
919          if (rnemdFluxType_ ==  rnemdFullKE) {
920            x = 1.0 + px * (1.0 - c);
# Line 950 | Line 952 | namespace OpenMD {
952              }
953            }
954            w = sqrt(w);
953 –          // std::cerr << "xh= " << x << "\tyh= " << y << "\tzh= " << z
954 –          //           << "\twh= " << w << endl;
955            for (sdi = hotBin.begin(); sdi != hotBin.end(); sdi++) {
956              if (rnemdFluxType_ == rnemdFullKE) {
957                vel = (*sdi)->getVel();
# Line 1216 | 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 1260 | Line 1260 | namespace OpenMD {
1260        
1261          if (inA) {
1262            hotBin.push_back(sd);
1263 –          //std::cerr << "before, velocity = " << vel << endl;
1263            Ph += mass * vel;
1265 –          //std::cerr << "after, velocity = " << vel << endl;
1264            Mh += mass;
1265            Kh += mass * vel.lengthSquare();
1266            if (rnemdFluxType_ == rnemdFullKE) {
# Line 1310 | Line 1308 | namespace OpenMD {
1308      
1309      Kh *= 0.5;
1310      Kc *= 0.5;
1313 –
1314 –    // std::cerr << "Mh= " << Mh << "\tKh= " << Kh << "\tMc= " << Mc
1315 –    //        << "\tKc= " << Kc << endl;
1316 –    // 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 1345 | 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";
1349 <                  // 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 1394 | Line 1387 | namespace OpenMD {
1387    }
1388  
1389    void RNEMD::doRNEMD() {
1390 <
1390 >    if (!doRNEMD_) return;
1391      trialCount_++;
1392      switch(rnemdMethod_) {
1393      case rnemdSwap:
# Line 1413 | 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  
# Line 1421 | Line 1414 | namespace OpenMD {
1414  
1415      seleMan_.setSelectionSet(evaluator_.evaluate());
1416  
1417 <    int selei;
1417 >    int selei(0);
1418      StuntDouble* sd;
1419      int idx;
1420  
# Line 1447 | 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 1466 | Line 1460 | namespace OpenMD {
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 1497 | Line 1491 | namespace OpenMD {
1491        }
1492      }
1493      
1500 –
1494   #ifdef IS_MPI
1495      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binCount[0],
1496                                nBins_, MPI::INT, MPI::SUM);
# Line 1524 | 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_ / currentSnap_->getVolume();
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 1550 | 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 1575 | 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 1596 | Line 1595 | namespace OpenMD {
1595        RealType time = currentSnap_->getTime();
1596        RealType avgArea;
1597        areaAccumulator_->getAverage(avgArea);
1598 <      RealType Jz = kineticExchange_ / (2.0 * time * avgArea);
1598 >      RealType Jz = kineticExchange_ / (2.0 * time * avgArea)
1599 >        / PhysicalConstants::energyConvert;
1600        Vector3d JzP = momentumExchange_ / (2.0 * time * avgArea);      
1601  
1602        rnemdFile_ << "#######################################################\n";
# Line 1625 | Line 1625 | namespace OpenMD {
1625        rnemdFile_ << "# RNEMD report:\n";      
1626        rnemdFile_ << "#     running time = " << time << " fs\n";
1627        rnemdFile_ << "#     target flux:\n";
1628 <      rnemdFile_ << "#         kinetic = " << kineticFlux_ << "\n";
1629 <      rnemdFile_ << "#         momentum = " << momentumFluxVector_ << "\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 = " << kineticTarget_ << "\n";
1635 <      rnemdFile_ << "#         momentum = " << momentumTarget_ << "\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 = " << kineticExchange_ << "\n";
1641 <      rnemdFile_ << "#         momentum = " << momentumExchange_  << "\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 << "\n";
1647 <      rnemdFile_ << "#         momentum = " << JzP  << "\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";    
# Line 1653 | 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 1716 | 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 < nBins_);
1735      RealType s;
# Line 1734 | 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 < nBins_);
1754      Vector3d s;
# Line 1752 | Line 1767 | namespace OpenMD {
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;
# Line 1770 | Line 1786 | namespace OpenMD {
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;

Diff Legend

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