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 1775 by gezelter, Wed Aug 8 18:45:52 2012 UTC

# Line 287 | Line 287 | namespace OpenMD {
287        painCave.severity = OPENMD_WARNING;
288        simError();
289      }
290 +
291 +    areaAccumulator_ = new Accumulator();
292  
293      nBins_ = rnemdParams->getOutputBins();
294  
# Line 381 | Line 383 | namespace OpenMD {
383      // dt = exchange time interval
384      // flux = target flux
385  
386 <    kineticTarget_ = 2.0*kineticFlux_*exchangeTime_*hmat(0,0)*hmat(1,1);
387 <    momentumTarget_ = 2.0*momentumFluxVector_*exchangeTime_*hmat(0,0)*hmat(1,1);
386 >    RealType area = currentSnap_->getXYarea();
387 >    kineticTarget_ = 2.0 * kineticFlux_ * exchangeTime_ * area;
388 >    momentumTarget_ = 2.0 * momentumFluxVector_ * exchangeTime_ * area;
389  
390      // total exchange sums are zeroed out at the beginning:
391  
# Line 1414 | Line 1417 | namespace OpenMD {
1417      Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
1418      Mat3x3d hmat = currentSnap_->getHmat();
1419  
1420 +    areaAccumulator_->add(currentSnap_->getXYarea());
1421 +
1422      seleMan_.setSelectionSet(evaluator_.evaluate());
1423  
1424      int selei;
# Line 1454 | Line 1459 | namespace OpenMD {
1459        if (usePeriodicBoundaryConditions_)
1460          currentSnap_->wrapVector(pos);
1461  
1462 +
1463        // which bin is this stuntdouble in?
1464        // wrapped positions are in the range [-0.5*hmat(2,2), +0.5*hmat(2,2)]
1465        // Shift molecules by half a box to have bins start at 0
# Line 1518 | Line 1524 | namespace OpenMD {
1524        vel.x() = binPx[i] / binMass[i];
1525        vel.y() = binPy[i] / binMass[i];
1526        vel.z() = binPz[i] / binMass[i];
1527 <      den = binCount[i] * nBins_ / (hmat(0,0) * hmat(1,1) * hmat(2,2));
1527 >      den = binCount[i] * nBins_ / currentSnap_->getVolume();
1528        temp = 2.0 * binKE[i] / (binDOF[i] * PhysicalConstants::kb *
1529                                 PhysicalConstants::energyConvert);
1530  
# Line 1588 | Line 1594 | namespace OpenMD {
1594        Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
1595  
1596        RealType time = currentSnap_->getTime();
1597 <      
1598 <      
1597 >      RealType avgArea;
1598 >      areaAccumulator_->getAverage(avgArea);
1599 >      RealType Jz = kineticExchange_ / (2.0 * time * avgArea);
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 <
1633 <      
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";
1630 >      rnemdFile_ << "#     target one-time exchanges:\n";
1631 >      rnemdFile_ << "#         kinetic = " << kineticTarget_ << "\n";
1632 >      rnemdFile_ << "#         momentum = " << momentumTarget_ << "\n";
1633 >      rnemdFile_ << "#     actual exchange totals:\n";
1634 >      rnemdFile_ << "#         kinetic = " << kineticExchange_ << "\n";
1635 >      rnemdFile_ << "#         momentum = " << momentumExchange_  << "\n";
1636 >      rnemdFile_ << "#     actual flux:\n";
1637 >      rnemdFile_ << "#         kinetic = " << Jz << "\n";
1638 >      rnemdFile_ << "#         momentum = " << JzP  << "\n";
1639 >      rnemdFile_ << "#     exchange statistics:\n";
1640 >      rnemdFile_ << "#         attempted = " << trialCount_ << "\n";
1641 >      rnemdFile_ << "#         failed = " << failTrialCount_ << "\n";    
1642        if (rnemdMethod_ == rnemdNIVS) {
1643 <        rnemdFile_ << "# NIVS root-check warnings: " << failRootCount_ << "\n";
1643 >        rnemdFile_ << "#         NIVS root-check errors = "
1644 >                   << failRootCount_ << "\n";
1645        }
1637 –
1646        rnemdFile_ << "#######################################################\n";
1647        
1648        
# Line 1671 | Line 1679 | namespace OpenMD {
1679          rnemdFile_ << std::endl;
1680          
1681        }        
1682 +
1683 +      rnemdFile_ << "#######################################################\n";
1684 +      rnemdFile_ << "# Standard Deviations in those quantities follow:\n";
1685 +      rnemdFile_ << "#######################################################\n";
1686 +
1687 +
1688 +      for (unsigned int j = 0; j < nBins_; j++) {        
1689 +        rnemdFile_ << "#";
1690 +        for (unsigned int i = 0; i < outputMask_.size(); ++i) {
1691 +          if (outputMask_[i]) {
1692 +            if (data_[i].dataType == "RealType")
1693 +              writeRealStdDev(i,j);
1694 +            else if (data_[i].dataType == "Vector3d")
1695 +              writeVectorStdDev(i,j);
1696 +            else {
1697 +              sprintf( painCave.errMsg,
1698 +                       "RNEMD found an unknown data type for: %s ",
1699 +                       data_[i].title.c_str());
1700 +              painCave.isFatal = 1;
1701 +              simError();
1702 +            }
1703 +          }
1704 +        }
1705 +        rnemdFile_ << std::endl;
1706 +        
1707 +      }        
1708        
1709        rnemdFile_.flush();
1710        rnemdFile_.close();
# Line 1683 | Line 1717 | namespace OpenMD {
1717    
1718    void RNEMD::writeReal(int index, unsigned int bin) {
1719      assert(index >=0 && index < ENDINDEX);
1720 <    assert(bin >=0 && bin < nBins_);
1720 >    assert(bin < nBins_);
1721      RealType s;
1722      
1723      data_[index].accumulator[bin]->getAverage(s);
# Line 1701 | Line 1735 | namespace OpenMD {
1735    
1736    void RNEMD::writeVector(int index, unsigned int bin) {
1737      assert(index >=0 && index < ENDINDEX);
1738 <    assert(bin >=0 && bin < nBins_);
1738 >    assert(bin < nBins_);
1739      Vector3d s;
1740      dynamic_cast<VectorAccumulator*>(data_[index].accumulator[bin])->getAverage(s);
1741      if (isinf(s[0]) || isnan(s[0]) ||
# Line 1716 | Line 1750 | namespace OpenMD {
1750        rnemdFile_ << "\t" << s[0] << "\t" << s[1] << "\t" << s[2];
1751      }
1752    }  
1753 +
1754 +  void RNEMD::writeRealStdDev(int index, unsigned int bin) {
1755 +    assert(index >=0 && index < ENDINDEX);
1756 +    assert(bin < nBins_);
1757 +    RealType s;
1758 +    
1759 +    data_[index].accumulator[bin]->getStdDev(s);
1760 +    
1761 +    if (! isinf(s) && ! isnan(s)) {
1762 +      rnemdFile_ << "\t" << s;
1763 +    } else{
1764 +      sprintf( painCave.errMsg,
1765 +               "RNEMD detected a numerical error writing: %s std. dev. for bin %d",
1766 +               data_[index].title.c_str(), bin);
1767 +      painCave.isFatal = 1;
1768 +      simError();
1769 +    }    
1770 +  }
1771 +  
1772 +  void RNEMD::writeVectorStdDev(int index, unsigned int bin) {
1773 +    assert(index >=0 && index < ENDINDEX);
1774 +    assert(bin < nBins_);
1775 +    Vector3d s;
1776 +    dynamic_cast<VectorAccumulator*>(data_[index].accumulator[bin])->getStdDev(s);
1777 +    if (isinf(s[0]) || isnan(s[0]) ||
1778 +        isinf(s[1]) || isnan(s[1]) ||
1779 +        isinf(s[2]) || isnan(s[2]) ) {      
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 +    } else {
1786 +      rnemdFile_ << "\t" << s[0] << "\t" << s[1] << "\t" << s[2];
1787 +    }
1788 +  }  
1789   }
1790  

Diff Legend

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