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 1776 by gezelter, Thu Aug 9 15:52:59 2012 UTC vs.
Revision 1812 by gezelter, Fri Nov 16 21:18:42 2012 UTC

# Line 80 | 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 101 | 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 200 | 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 227 | Line 229 | namespace OpenMD {
229      }
230      if (!hasCorrectFlux) {
231        sprintf(painCave.errMsg,
232 <              "RNEMD: The current method,\n"
231 <              "\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 238 | 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 273 | 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 317 | 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 494 | Line 499 | namespace OpenMD {
499                  + angMom[2]*angMom[2]/I(2, 2);
500              }
501            } //angular momenta exchange enabled
497 –          //energyConvert temporarily disabled
498 –          //make kineticExchange_ comparable between swap & scale
499 –          //value = value * 0.5 / PhysicalConstants::energyConvert;
502            value *= 0.5;
503            break;
504          case rnemdPx :
# Line 734 | Line 736 | namespace OpenMD {
736          
737          switch(rnemdFluxType_) {
738          case rnemdKE:
737 –          cerr << "KE\n";
739            kineticExchange_ += max_val - min_val;
740            break;
741          case rnemdPx:
# Line 747 | Line 748 | namespace OpenMD {
748            momentumExchange_.z() += max_val - min_val;
749            break;
750          default:
750 –          cerr << "default\n";
751            break;
752          }
753        } else {        
# Line 914 | 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;
918 <        //now convert to hotBin coefficient
917 >
918          RealType w = 0.0;
919          if (rnemdFluxType_ ==  rnemdFullKE) {
920            x = 1.0 + px * (1.0 - c);
# Line 953 | Line 952 | namespace OpenMD {
952              }
953            }
954            w = sqrt(w);
956 –          // std::cerr << "xh= " << x << "\tyh= " << y << "\tzh= " << z
957 –          //           << "\twh= " << w << endl;
955            for (sdi = hotBin.begin(); sdi != hotBin.end(); sdi++) {
956              if (rnemdFluxType_ == rnemdFullKE) {
957                vel = (*sdi)->getVel();
# Line 1263 | Line 1260 | namespace OpenMD {
1260        
1261          if (inA) {
1262            hotBin.push_back(sd);
1266 –          //std::cerr << "before, velocity = " << vel << endl;
1263            Ph += mass * vel;
1268 –          //std::cerr << "after, velocity = " << vel << endl;
1264            Mh += mass;
1265            Kh += mass * vel.lengthSquare();
1266            if (rnemdFluxType_ == rnemdFullKE) {
# Line 1313 | Line 1308 | namespace OpenMD {
1308      
1309      Kh *= 0.5;
1310      Kc *= 0.5;
1316 –
1317 –    // std::cerr << "Mh= " << Mh << "\tKh= " << Kh << "\tMc= " << Mc
1318 –    //        << "\tKc= " << Kc << endl;
1319 –    // 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 1348 | 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";
1352 <                  // 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 1424 | 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 1450 | 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 1469 | 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 1500 | Line 1491 | namespace OpenMD {
1491        }
1492      }
1493      
1503 –
1494   #ifdef IS_MPI
1495      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binCount[0],
1496                                nBins_, MPI::INT, MPI::SUM);
# Line 1527 | 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];
1530 –      den = binCount[i] * nBins_ / currentSnap_->getVolume();
1531 –      temp = 2.0 * binKE[i] / (binDOF[i] * PhysicalConstants::kb *
1532 –                               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 1602 | Line 1598 | namespace OpenMD {
1598        RealType time = currentSnap_->getTime();
1599        RealType avgArea;
1600        areaAccumulator_->getAverage(avgArea);
1601 <      RealType Jz = kineticExchange_ / (2.0 * time * avgArea);
1601 >      RealType Jz = kineticExchange_ / (2.0 * time * avgArea)
1602 >        / PhysicalConstants::energyConvert;
1603        Vector3d JzP = momentumExchange_ / (2.0 * time * avgArea);      
1604  
1605        rnemdFile_ << "#######################################################\n";
# Line 1631 | Line 1628 | namespace OpenMD {
1628        rnemdFile_ << "# RNEMD report:\n";      
1629        rnemdFile_ << "#     running time = " << time << " fs\n";
1630        rnemdFile_ << "#     target flux:\n";
1631 <      rnemdFile_ << "#         kinetic = " << kineticFlux_ << "\n";
1632 <      rnemdFile_ << "#         momentum = " << momentumFluxVector_ << "\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 = " << kineticTarget_ << "\n";
1638 <      rnemdFile_ << "#         momentum = " << momentumTarget_ << "\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 = " << kineticExchange_ << "\n";
1644 <      rnemdFile_ << "#         momentum = " << momentumExchange_  << "\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 << "\n";
1650 <      rnemdFile_ << "#         momentum = " << JzP  << "\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";    
# Line 1659 | 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;

Diff Legend

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