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 1774 by gezelter, Wed Aug 8 16:03:02 2012 UTC vs.
trunk/src/rnemd/RNEMD.cpp (file contents), Revision 1793 by gezelter, Fri Aug 31 21:16:10 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 298 | 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 308 | 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 328 | 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 410 | 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 432 | 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 440 | Line 451 | namespace OpenMD {
451  
452      int selei;
453      StuntDouble* sd;
443 –    int idx;
454  
455      RealType min_val;
456      bool min_found = false;  
# Line 453 | Line 463 | namespace OpenMD {
463      for (sd = seleMan_.beginSelected(selei); sd != NULL;
464           sd = seleMan_.nextSelected(selei)) {
465  
456 –      idx = sd->getLocalIndex();
457 –
466        Vector3d pos = sd->getPos();
467  
468        // wrap the stuntdouble's position back into the box:
# 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 535 | Line 540 | namespace OpenMD {
540        }
541      }
542      
543 < #ifdef IS_MPI
544 <    int nProc, worldRank;
543 > #ifdef IS_MPI    
544 >    int worldRank = MPI::COMM_WORLD.Get_rank();
545      
541 –    nProc = MPI::COMM_WORLD.Get_size();
542 –    worldRank = MPI::COMM_WORLD.Get_rank();
543 –
546      bool my_min_found = min_found;
547      bool my_max_found = max_found;
548  
# Line 731 | Line 733 | namespace OpenMD {
733          
734          switch(rnemdFluxType_) {
735          case rnemdKE:
734 –          cerr << "KE\n";
736            kineticExchange_ += max_val - min_val;
737            break;
738          case rnemdPx:
# Line 744 | Line 745 | namespace OpenMD {
745            momentumExchange_.z() += max_val - min_val;
746            break;
747          default:
747 –          cerr << "default\n";
748            break;
749          }
750        } else {        
# Line 767 | Line 767 | namespace OpenMD {
767    }
768    
769    void RNEMD::doNIVS() {
770 <
770 >    if (!doRNEMD_) return;
771      Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
772      Mat3x3d hmat = currentSnap_->getHmat();
773  
# Line 775 | Line 775 | namespace OpenMD {
775  
776      int selei;
777      StuntDouble* sd;
778 –    int idx;
778  
779      vector<StuntDouble*> hotBin, coldBin;
780  
# Line 796 | Line 795 | namespace OpenMD {
795  
796      for (sd = seleMan_.beginSelected(selei); sd != NULL;
797           sd = seleMan_.nextSelected(selei)) {
799 –
800 –      idx = sd->getLocalIndex();
798  
799        Vector3d pos = sd->getPos();
800  
# Line 1216 | Line 1213 | namespace OpenMD {
1213    }
1214  
1215    void RNEMD::doVSS() {
1216 <
1216 >    if (!doRNEMD_) return;
1217      Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
1218      RealType time = currentSnap_->getTime();    
1219      Mat3x3d hmat = currentSnap_->getHmat();
# Line 1225 | Line 1222 | namespace OpenMD {
1222  
1223      int selei;
1224      StuntDouble* sd;
1228 –    int idx;
1225  
1226      vector<StuntDouble*> hotBin, coldBin;
1227  
# Line 1240 | Line 1236 | namespace OpenMD {
1236      for (sd = seleMan_.beginSelected(selei); sd != NULL;
1237           sd = seleMan_.nextSelected(selei)) {
1238  
1243 –      idx = sd->getLocalIndex();
1244 –
1239        Vector3d pos = sd->getPos();
1240  
1241        // wrap the stuntdouble's position back into the box:
# Line 1394 | Line 1388 | namespace OpenMD {
1388    }
1389  
1390    void RNEMD::doRNEMD() {
1391 <
1391 >    if (!doRNEMD_) return;
1392      trialCount_++;
1393      switch(rnemdMethod_) {
1394      case rnemdSwap:
# Line 1413 | Line 1407 | namespace OpenMD {
1407    }
1408  
1409    void RNEMD::collectData() {
1410 <
1410 >    if (!doRNEMD_) return;
1411      Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
1412      Mat3x3d hmat = currentSnap_->getHmat();
1413  
# Line 1423 | Line 1417 | namespace OpenMD {
1417  
1418      int selei;
1419      StuntDouble* sd;
1426 –    int idx;
1420  
1421      vector<RealType> binMass(nBins_, 0.0);
1422      vector<RealType> binPx(nBins_, 0.0);
# Line 1448 | Line 1441 | namespace OpenMD {
1441               sd = mol->nextIntegrableObject(iiter))
1442      */
1443      for (sd = seleMan_.beginSelected(selei); sd != NULL;
1444 <         sd = seleMan_.nextSelected(selei)) {
1452 <      
1453 <      idx = sd->getLocalIndex();
1444 >         sd = seleMan_.nextSelected(selei)) {    
1445        
1446        Vector3d pos = sd->getPos();
1447  
# Line 1524 | Line 1515 | namespace OpenMD {
1515        vel.x() = binPx[i] / binMass[i];
1516        vel.y() = binPy[i] / binMass[i];
1517        vel.z() = binPz[i] / binMass[i];
1518 <      den = binCount[i] * nBins_ / currentSnap_->getVolume();
1518 >
1519 >      den = binMass[i] * nBins_ * PhysicalConstants::densityConvert
1520 >        / currentSnap_->getVolume() ;
1521 >
1522        temp = 2.0 * binKE[i] / (binDOF[i] * PhysicalConstants::kb *
1523                                 PhysicalConstants::energyConvert);
1524  
# Line 1532 | Line 1526 | namespace OpenMD {
1526          if(outputMask_[j]) {
1527            switch(j) {
1528            case Z:
1529 <            (data_[j].accumulator[i])->add(z);
1529 >            dynamic_cast<Accumulator *>(data_[j].accumulator[i])->add(z);
1530              break;
1531            case TEMPERATURE:
1532 <            data_[j].accumulator[i]->add(temp);
1532 >            dynamic_cast<Accumulator *>(data_[j].accumulator[i])->add(temp);
1533              break;
1534            case VELOCITY:
1535              dynamic_cast<VectorAccumulator *>(data_[j].accumulator[i])->add(vel);
1536              break;
1537            case DENSITY:
1538 <            data_[j].accumulator[i]->add(den);
1538 >            dynamic_cast<Accumulator *>(data_[j].accumulator[i])->add(den);
1539              break;
1540            }
1541          }
# Line 1550 | Line 1544 | namespace OpenMD {
1544    }
1545  
1546    void RNEMD::getStarted() {
1547 +    if (!doRNEMD_) return;
1548      collectData();
1549      writeOutputFile();
1550    }
1551  
1552    void RNEMD::parseOutputFileFormat(const std::string& format) {
1553 +    if (!doRNEMD_) return;
1554      StringTokenizer tokenizer(format, " ,;|\t\n\r");
1555      
1556      while(tokenizer.hasMoreTokens()) {
# Line 1575 | Line 1571 | namespace OpenMD {
1571    }
1572    
1573    void RNEMD::writeOutputFile() {
1574 +    if (!doRNEMD_) return;
1575      
1576   #ifdef IS_MPI
1577      // If we're the root node, should we print out the results
# Line 1596 | Line 1593 | namespace OpenMD {
1593        RealType time = currentSnap_->getTime();
1594        RealType avgArea;
1595        areaAccumulator_->getAverage(avgArea);
1596 <      RealType Jz = kineticExchange_ / (2.0 * time * avgArea);
1596 >      RealType Jz = kineticExchange_ / (2.0 * time * avgArea)
1597 >        / PhysicalConstants::energyConvert;
1598        Vector3d JzP = momentumExchange_ / (2.0 * time * avgArea);      
1599  
1600        rnemdFile_ << "#######################################################\n";
# Line 1613 | Line 1611 | namespace OpenMD {
1611            rnemdFile_ << "#    fluxType  = \"" << (*fi).first << "\";\n";
1612        }
1613        
1614 <      rnemdFile_ << "#    exchangeTime = " << exchangeTime_ << "; fs\n";
1614 >      rnemdFile_ << "#    exchangeTime = " << exchangeTime_ << ";\n";
1615  
1616        rnemdFile_ << "#    objectSelection = \""
1617                   << rnemdObjectSelection_ << "\";\n";
# Line 1625 | Line 1623 | namespace OpenMD {
1623        rnemdFile_ << "# RNEMD report:\n";      
1624        rnemdFile_ << "#     running time = " << time << " fs\n";
1625        rnemdFile_ << "#     target flux:\n";
1626 <      rnemdFile_ << "#         kinetic = " << kineticFlux_ << "\n";
1627 <      rnemdFile_ << "#         momentum = " << momentumFluxVector_ << "\n";
1626 >      rnemdFile_ << "#         kinetic = "
1627 >                 << kineticFlux_ / PhysicalConstants::energyConvert
1628 >                 << " (kcal/mol/A^2/fs)\n";
1629 >      rnemdFile_ << "#         momentum = " << momentumFluxVector_
1630 >                 << " (amu/A/fs^2)\n";
1631        rnemdFile_ << "#     target one-time exchanges:\n";
1632 <      rnemdFile_ << "#         kinetic = " << kineticTarget_ << "\n";
1633 <      rnemdFile_ << "#         momentum = " << momentumTarget_ << "\n";
1632 >      rnemdFile_ << "#         kinetic = "
1633 >                 << kineticTarget_ / PhysicalConstants::energyConvert
1634 >                 << " (kcal/mol)\n";
1635 >      rnemdFile_ << "#         momentum = " << momentumTarget_
1636 >                 << " (amu*A/fs)\n";
1637        rnemdFile_ << "#     actual exchange totals:\n";
1638 <      rnemdFile_ << "#         kinetic = " << kineticExchange_ << "\n";
1639 <      rnemdFile_ << "#         momentum = " << momentumExchange_  << "\n";
1638 >      rnemdFile_ << "#         kinetic = "
1639 >                 << kineticExchange_ / PhysicalConstants::energyConvert
1640 >                 << " (kcal/mol)\n";
1641 >      rnemdFile_ << "#         momentum = " << momentumExchange_
1642 >                 << " (amu*A/fs)\n";      
1643        rnemdFile_ << "#     actual flux:\n";
1644 <      rnemdFile_ << "#         kinetic = " << Jz << "\n";
1645 <      rnemdFile_ << "#         momentum = " << JzP  << "\n";
1644 >      rnemdFile_ << "#         kinetic = " << Jz
1645 >                 << " (kcal/mol/A^2/fs)\n";
1646 >      rnemdFile_ << "#         momentum = " << JzP
1647 >                 << " (amu/A/fs^2)\n";
1648        rnemdFile_ << "#     exchange statistics:\n";
1649        rnemdFile_ << "#         attempted = " << trialCount_ << "\n";
1650        rnemdFile_ << "#         failed = " << failTrialCount_ << "\n";    
# Line 1653 | Line 1662 | namespace OpenMD {
1662          if (outputMask_[i]) {
1663            rnemdFile_ << "\t" << data_[i].title <<
1664              "(" << data_[i].units << ")";
1665 +          // add some extra tabs for column alignment
1666 +          if (data_[i].dataType == "Vector3d") rnemdFile_ << "\t\t";
1667          }
1668        }
1669        rnemdFile_ << std::endl;
1670        
1671        rnemdFile_.precision(8);
1672        
1673 <      for (unsigned int j = 0; j < nBins_; j++) {        
1673 >      for (int j = 0; j < nBins_; j++) {        
1674          
1675          for (unsigned int i = 0; i < outputMask_.size(); ++i) {
1676            if (outputMask_[i]) {
# Line 1685 | Line 1696 | namespace OpenMD {
1696        rnemdFile_ << "#######################################################\n";
1697  
1698  
1699 <      for (unsigned int j = 0; j < nBins_; j++) {        
1699 >      for (int j = 0; j < nBins_; j++) {        
1700          rnemdFile_ << "#";
1701          for (unsigned int i = 0; i < outputMask_.size(); ++i) {
1702            if (outputMask_[i]) {
# Line 1716 | Line 1727 | namespace OpenMD {
1727    }
1728    
1729    void RNEMD::writeReal(int index, unsigned int bin) {
1730 +    if (!doRNEMD_) return;
1731      assert(index >=0 && index < ENDINDEX);
1732      assert(bin < nBins_);
1733      RealType s;
1734      
1735 <    data_[index].accumulator[bin]->getAverage(s);
1735 >    dynamic_cast<Accumulator *>(data_[index].accumulator[bin])->getAverage(s);
1736      
1737      if (! isinf(s) && ! isnan(s)) {
1738        rnemdFile_ << "\t" << s;
# Line 1734 | Line 1746 | namespace OpenMD {
1746    }
1747    
1748    void RNEMD::writeVector(int index, unsigned int bin) {
1749 +    if (!doRNEMD_) return;
1750      assert(index >=0 && index < ENDINDEX);
1751      assert(bin < nBins_);
1752      Vector3d s;
# Line 1752 | Line 1765 | namespace OpenMD {
1765    }  
1766  
1767    void RNEMD::writeRealStdDev(int index, unsigned int bin) {
1768 +    if (!doRNEMD_) return;
1769      assert(index >=0 && index < ENDINDEX);
1770      assert(bin < nBins_);
1771      RealType s;
1772      
1773 <    data_[index].accumulator[bin]->getStdDev(s);
1773 >    dynamic_cast<Accumulator *>(data_[index].accumulator[bin])->getStdDev(s);
1774      
1775      if (! isinf(s) && ! isnan(s)) {
1776        rnemdFile_ << "\t" << s;
# Line 1770 | Line 1784 | namespace OpenMD {
1784    }
1785    
1786    void RNEMD::writeVectorStdDev(int index, unsigned int bin) {
1787 +    if (!doRNEMD_) return;
1788      assert(index >=0 && index < ENDINDEX);
1789      assert(bin < nBins_);
1790      Vector3d s;

Diff Legend

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