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 1777 by gezelter, Thu Aug 9 18:35:09 2012 UTC vs.
Revision 1812 by gezelter, Fri Nov 16 21:18:42 2012 UTC

# Line 277 | 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 913 | 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;
917 <        //now convert to hotBin coefficient
917 >
918          RealType w = 0.0;
919          if (rnemdFluxType_ ==  rnemdFullKE) {
920            x = 1.0 + px * (1.0 - c);
# Line 952 | Line 952 | namespace OpenMD {
952              }
953            }
954            w = sqrt(w);
955 –          // std::cerr << "xh= " << x << "\tyh= " << y << "\tzh= " << z
956 –          //           << "\twh= " << w << endl;
955            for (sdi = hotBin.begin(); sdi != hotBin.end(); sdi++) {
956              if (rnemdFluxType_ == rnemdFullKE) {
957                vel = (*sdi)->getVel();
# Line 1262 | Line 1260 | namespace OpenMD {
1260        
1261          if (inA) {
1262            hotBin.push_back(sd);
1265 –          //std::cerr << "before, velocity = " << vel << endl;
1263            Ph += mass * vel;
1267 –          //std::cerr << "after, velocity = " << vel << endl;
1264            Mh += mass;
1265            Kh += mass * vel.lengthSquare();
1266            if (rnemdFluxType_ == rnemdFullKE) {
# Line 1312 | Line 1308 | namespace OpenMD {
1308      
1309      Kh *= 0.5;
1310      Kc *= 0.5;
1315 –
1316 –    // std::cerr << "Mh= " << Mh << "\tKh= " << Kh << "\tMc= " << Mc
1317 –    //        << "\tKc= " << Kc << endl;
1318 –    // 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 1347 | 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";
1351 <                  // 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 1423 | 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 1449 | 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 1468 | 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 1499 | Line 1491 | namespace OpenMD {
1491        }
1492      }
1493      
1502 –
1494   #ifdef IS_MPI
1495      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binCount[0],
1496                                nBins_, MPI::INT, MPI::SUM);
# Line 1530 | Line 1521 | namespace OpenMD {
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 <
1527 <      for (unsigned int j = 0; j < outputMask_.size(); ++j) {
1528 <        if(outputMask_[j]) {
1529 <          switch(j) {
1530 <          case Z:
1531 <            (data_[j].accumulator[i])->add(z);
1532 <            break;
1533 <          case TEMPERATURE:
1534 <            data_[j].accumulator[i]->add(temp);
1535 <            break;
1536 <          case VELOCITY:
1537 <            dynamic_cast<VectorAccumulator *>(data_[j].accumulator[i])->add(vel);
1538 <            break;
1539 <          case DENSITY:
1540 <            data_[j].accumulator[i]->add(den);
1541 <            break;
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        }

Diff Legend

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