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 1812 by gezelter, Fri Nov 16 21:18:42 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 1453 | Line 1452 | namespace OpenMD {
1452        
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)]
# Line 1460 | 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 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];
1521 –      den = binCount[i] * nBins_ / (hmat(0,0) * hmat(1,1) * hmat(2,2));
1522 –      temp = 2.0 * binKE[i] / (binDOF[i] * PhysicalConstants::kb *
1523 –                               PhysicalConstants::energyConvert);
1520  
1521 <      for (unsigned int j = 0; j < outputMask_.size(); ++j) {
1522 <        if(outputMask_[j]) {
1523 <          switch(j) {
1524 <          case Z:
1525 <            (data_[j].accumulator[i])->add(z);
1526 <            break;
1527 <          case TEMPERATURE:
1528 <            data_[j].accumulator[i]->add(temp);
1529 <            break;
1530 <          case VELOCITY:
1531 <            dynamic_cast<VectorAccumulator *>(data_[j].accumulator[i])->add(vel);
1532 <            break;
1533 <          case DENSITY:
1534 <            data_[j].accumulator[i]->add(den);
1535 <            break;
1521 >      den = binMass[i] * nBins_ * PhysicalConstants::densityConvert
1522 >        / currentSnap_->getVolume() ;
1523 >
1524 >      if (binCount[i] > 0) {
1525 >        // only add values if there are things to add
1526 >        temp = 2.0 * binKE[i] / (binDOF[i] * PhysicalConstants::kb *
1527 >                                 PhysicalConstants::energyConvert);
1528 >        
1529 >        for (unsigned int j = 0; j < outputMask_.size(); ++j) {
1530 >          if(outputMask_[j]) {
1531 >            switch(j) {
1532 >            case Z:
1533 >              dynamic_cast<Accumulator *>(data_[j].accumulator[i])->add(z);
1534 >              break;
1535 >            case TEMPERATURE:
1536 >              dynamic_cast<Accumulator *>(data_[j].accumulator[i])->add(temp);
1537 >              break;
1538 >            case VELOCITY:
1539 >              dynamic_cast<VectorAccumulator *>(data_[j].accumulator[i])->add(vel);
1540 >              break;
1541 >            case DENSITY:
1542 >              dynamic_cast<Accumulator *>(data_[j].accumulator[i])->add(den);
1543 >              break;
1544 >            }
1545            }
1546          }
1547        }
# Line 1544 | Line 1549 | namespace OpenMD {
1549    }
1550  
1551    void RNEMD::getStarted() {
1552 +    if (!doRNEMD_) return;
1553      collectData();
1554      writeOutputFile();
1555    }
1556  
1557    void RNEMD::parseOutputFileFormat(const std::string& format) {
1558 +    if (!doRNEMD_) return;
1559      StringTokenizer tokenizer(format, " ,;|\t\n\r");
1560      
1561      while(tokenizer.hasMoreTokens()) {
# Line 1569 | Line 1576 | namespace OpenMD {
1576    }
1577    
1578    void RNEMD::writeOutputFile() {
1579 +    if (!doRNEMD_) return;
1580      
1581   #ifdef IS_MPI
1582      // If we're the root node, should we print out the results
# Line 1588 | Line 1596 | namespace OpenMD {
1596        Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
1597  
1598        RealType time = currentSnap_->getTime();
1599 <      
1600 <      
1599 >      RealType avgArea;
1600 >      areaAccumulator_->getAverage(avgArea);
1601 >      RealType Jz = kineticExchange_ / (2.0 * time * avgArea)
1602 >        / PhysicalConstants::energyConvert;
1603 >      Vector3d JzP = momentumExchange_ / (2.0 * time * avgArea);      
1604 >
1605        rnemdFile_ << "#######################################################\n";
1606        rnemdFile_ << "# RNEMD {\n";
1607  
1608        map<string, RNEMDMethod>::iterator mi;
1609        for(mi = stringToMethod_.begin(); mi != stringToMethod_.end(); ++mi) {
1610          if ( (*mi).second == rnemdMethod_)
1611 <          rnemdFile_ << "#    exchangeMethod  = " << (*mi).first << "\n";
1611 >          rnemdFile_ << "#    exchangeMethod  = \"" << (*mi).first << "\";\n";
1612        }
1613        map<string, RNEMDFluxType>::iterator fi;
1614        for(fi = stringToFluxType_.begin(); fi != stringToFluxType_.end(); ++fi) {
1615          if ( (*fi).second == rnemdFluxType_)
1616 <          rnemdFile_ << "#    fluxType  = " << (*fi).first << "\n";
1616 >          rnemdFile_ << "#    fluxType  = \"" << (*fi).first << "\";\n";
1617        }
1618        
1619 <      rnemdFile_ << "#    exchangeTime = " << exchangeTime_ << " fs\n";
1619 >      rnemdFile_ << "#    exchangeTime = " << exchangeTime_ << ";\n";
1620  
1621        rnemdFile_ << "#    objectSelection = \""
1622 <                 << rnemdObjectSelection_ << "\"\n";
1623 <      rnemdFile_ << "#    slabWidth = " << slabWidth_ << " angstroms\n";
1624 <      rnemdFile_ << "#    slabAcenter = " << slabACenter_ << " angstroms\n";
1625 <      rnemdFile_ << "#    slabBcenter = " << slabBCenter_ << " angstroms\n";
1622 >                 << rnemdObjectSelection_ << "\";\n";
1623 >      rnemdFile_ << "#    slabWidth = " << slabWidth_ << ";\n";
1624 >      rnemdFile_ << "#    slabAcenter = " << slabACenter_ << ";\n";
1625 >      rnemdFile_ << "#    slabBcenter = " << slabBCenter_ << ";\n";
1626        rnemdFile_ << "# }\n";
1627        rnemdFile_ << "#######################################################\n";
1628 <      
1629 <      rnemdFile_ << "# running time = " << time << " fs\n";
1630 <      rnemdFile_ << "# target kinetic flux = " << kineticFlux_ << "\n";
1631 <      rnemdFile_ << "# target momentum flux = " << momentumFluxVector_ << "\n";
1632 <      
1633 <      rnemdFile_ << "# target one-time kinetic exchange = " << kineticTarget_
1634 <                 << "\n";
1635 <      rnemdFile_ << "# target one-time momentum exchange = " << momentumTarget_
1636 <                 << "\n";
1637 <      
1638 <      rnemdFile_ << "# actual kinetic exchange = " << kineticExchange_ << "\n";
1639 <      rnemdFile_ << "# actual momentum exchange = " << momentumExchange_
1640 <                 << "\n";
1641 <      
1642 <      rnemdFile_ << "# attempted exchanges: " << trialCount_ << "\n";
1643 <      rnemdFile_ << "# failed exchanges: " << failTrialCount_ << "\n";
1644 <
1645 <      
1628 >      rnemdFile_ << "# RNEMD report:\n";      
1629 >      rnemdFile_ << "#     running time = " << time << " fs\n";
1630 >      rnemdFile_ << "#     target flux:\n";
1631 >      rnemdFile_ << "#         kinetic = "
1632 >                 << kineticFlux_ / PhysicalConstants::energyConvert
1633 >                 << " (kcal/mol/A^2/fs)\n";
1634 >      rnemdFile_ << "#         momentum = " << momentumFluxVector_
1635 >                 << " (amu/A/fs^2)\n";
1636 >      rnemdFile_ << "#     target one-time exchanges:\n";
1637 >      rnemdFile_ << "#         kinetic = "
1638 >                 << kineticTarget_ / PhysicalConstants::energyConvert
1639 >                 << " (kcal/mol)\n";
1640 >      rnemdFile_ << "#         momentum = " << momentumTarget_
1641 >                 << " (amu*A/fs)\n";
1642 >      rnemdFile_ << "#     actual exchange totals:\n";
1643 >      rnemdFile_ << "#         kinetic = "
1644 >                 << kineticExchange_ / PhysicalConstants::energyConvert
1645 >                 << " (kcal/mol)\n";
1646 >      rnemdFile_ << "#         momentum = " << momentumExchange_
1647 >                 << " (amu*A/fs)\n";      
1648 >      rnemdFile_ << "#     actual flux:\n";
1649 >      rnemdFile_ << "#         kinetic = " << Jz
1650 >                 << " (kcal/mol/A^2/fs)\n";
1651 >      rnemdFile_ << "#         momentum = " << JzP
1652 >                 << " (amu/A/fs^2)\n";
1653 >      rnemdFile_ << "#     exchange statistics:\n";
1654 >      rnemdFile_ << "#         attempted = " << trialCount_ << "\n";
1655 >      rnemdFile_ << "#         failed = " << failTrialCount_ << "\n";    
1656        if (rnemdMethod_ == rnemdNIVS) {
1657 <        rnemdFile_ << "# NIVS root-check warnings: " << failRootCount_ << "\n";
1657 >        rnemdFile_ << "#         NIVS root-check errors = "
1658 >                   << failRootCount_ << "\n";
1659        }
1637 –
1660        rnemdFile_ << "#######################################################\n";
1661        
1662        
# Line 1645 | Line 1667 | namespace OpenMD {
1667          if (outputMask_[i]) {
1668            rnemdFile_ << "\t" << data_[i].title <<
1669              "(" << data_[i].units << ")";
1670 +          // add some extra tabs for column alignment
1671 +          if (data_[i].dataType == "Vector3d") rnemdFile_ << "\t\t";
1672          }
1673        }
1674        rnemdFile_ << std::endl;
# Line 1671 | Line 1695 | namespace OpenMD {
1695          rnemdFile_ << std::endl;
1696          
1697        }        
1698 +
1699 +      rnemdFile_ << "#######################################################\n";
1700 +      rnemdFile_ << "# Standard Deviations in those quantities follow:\n";
1701 +      rnemdFile_ << "#######################################################\n";
1702 +
1703 +
1704 +      for (unsigned int j = 0; j < nBins_; j++) {        
1705 +        rnemdFile_ << "#";
1706 +        for (unsigned int i = 0; i < outputMask_.size(); ++i) {
1707 +          if (outputMask_[i]) {
1708 +            if (data_[i].dataType == "RealType")
1709 +              writeRealStdDev(i,j);
1710 +            else if (data_[i].dataType == "Vector3d")
1711 +              writeVectorStdDev(i,j);
1712 +            else {
1713 +              sprintf( painCave.errMsg,
1714 +                       "RNEMD found an unknown data type for: %s ",
1715 +                       data_[i].title.c_str());
1716 +              painCave.isFatal = 1;
1717 +              simError();
1718 +            }
1719 +          }
1720 +        }
1721 +        rnemdFile_ << std::endl;
1722 +        
1723 +      }        
1724        
1725        rnemdFile_.flush();
1726        rnemdFile_.close();
# Line 1682 | Line 1732 | namespace OpenMD {
1732    }
1733    
1734    void RNEMD::writeReal(int index, unsigned int bin) {
1735 +    if (!doRNEMD_) return;
1736      assert(index >=0 && index < ENDINDEX);
1737 <    assert(bin >=0 && bin < nBins_);
1737 >    assert(bin < nBins_);
1738      RealType s;
1739      
1740      data_[index].accumulator[bin]->getAverage(s);
# Line 1700 | Line 1751 | namespace OpenMD {
1751    }
1752    
1753    void RNEMD::writeVector(int index, unsigned int bin) {
1754 +    if (!doRNEMD_) return;
1755      assert(index >=0 && index < ENDINDEX);
1756 <    assert(bin >=0 && bin < nBins_);
1756 >    assert(bin < nBins_);
1757      Vector3d s;
1758      dynamic_cast<VectorAccumulator*>(data_[index].accumulator[bin])->getAverage(s);
1759      if (isinf(s[0]) || isnan(s[0]) ||
# Line 1716 | Line 1768 | namespace OpenMD {
1768        rnemdFile_ << "\t" << s[0] << "\t" << s[1] << "\t" << s[2];
1769      }
1770    }  
1771 +
1772 +  void RNEMD::writeRealStdDev(int index, unsigned int bin) {
1773 +    if (!doRNEMD_) return;
1774 +    assert(index >=0 && index < ENDINDEX);
1775 +    assert(bin < nBins_);
1776 +    RealType s;
1777 +    
1778 +    data_[index].accumulator[bin]->getStdDev(s);
1779 +    
1780 +    if (! isinf(s) && ! isnan(s)) {
1781 +      rnemdFile_ << "\t" << s;
1782 +    } else{
1783 +      sprintf( painCave.errMsg,
1784 +               "RNEMD detected a numerical error writing: %s std. dev. for bin %d",
1785 +               data_[index].title.c_str(), bin);
1786 +      painCave.isFatal = 1;
1787 +      simError();
1788 +    }    
1789 +  }
1790 +  
1791 +  void RNEMD::writeVectorStdDev(int index, unsigned int bin) {
1792 +    if (!doRNEMD_) return;
1793 +    assert(index >=0 && index < ENDINDEX);
1794 +    assert(bin < nBins_);
1795 +    Vector3d s;
1796 +    dynamic_cast<VectorAccumulator*>(data_[index].accumulator[bin])->getStdDev(s);
1797 +    if (isinf(s[0]) || isnan(s[0]) ||
1798 +        isinf(s[1]) || isnan(s[1]) ||
1799 +        isinf(s[2]) || isnan(s[2]) ) {      
1800 +      sprintf( painCave.errMsg,
1801 +               "RNEMD detected a numerical error writing: %s std. dev. for bin %d",
1802 +               data_[index].title.c_str(), bin);
1803 +      painCave.isFatal = 1;
1804 +      simError();
1805 +    } else {
1806 +      rnemdFile_ << "\t" << s[0] << "\t" << s[1] << "\t" << s[2];
1807 +    }
1808 +  }  
1809   }
1810  

Diff Legend

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