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/integrators/RNEMD.cpp (file contents):
Revision 1723 by gezelter, Thu May 24 20:59:54 2012 UTC vs.
Revision 1728 by jmarr, Wed May 30 16:07:03 2012 UTC

# Line 85 | Line 85 | namespace OpenMD {
85      stringToEnumMap_["ShiftScaleVAM"] = rnemdShiftScaleVAM;
86      stringToEnumMap_["Unknown"] = rnemdUnknown;
87  
88 +    runTime_ = simParams->getRunTime();
89 +    statusTime_ = simParams->getStatusTime();
90 +
91      rnemdObjectSelection_ = simParams->getRNEMD_objectSelection();
92      evaluator_.loadScriptString(rnemdObjectSelection_);
93      seleMan_.setSelectionSet(evaluator_.evaluate());
# Line 155 | Line 158 | namespace OpenMD {
158      if (simParams->haveRNEMD_outputRotTemperature()) {
159        outputRotTemp_ = simParams->getRNEMD_outputRotTemperature();
160      }
161 +    // James put this in.
162 +    outputDen_ = false;
163 +    if (simParams->haveRNEMD_outputDen()) {
164 +      outputDen_ = simParams->getRNEMD_outputDen();
165 +    }
166 +    outputAh_ = false;
167 +    if (simParams->haveRNEMD_outputAh()) {
168 +      outputAh_ = simParams->getRNEMD_outputAh();
169 +    }    
170 +    outputVz_ = false;
171 +    if (simParams->haveRNEMD_outputVz()) {
172 +      outputVz_ = simParams->getRNEMD_outputVz();
173 +    } else if ((rnemdType_ == rnemdPz) || (rnemdType_ == rnemdPzScale)) {
174 +      outputVz_ = true;
175 +    }
176 +    
177  
178   #ifdef IS_MPI
179      if (worldRank == 0) {
# Line 188 | Line 207 | namespace OpenMD {
207          rnemdFileName = "temperatureR.log";
208          rotTempLog_.open(rnemdFileName.c_str());
209        }
210 <
210 >      
211 >      //James put this in
212 >      if (outputDen_) {
213 >        rnemdFileName = "Density.log";
214 >        denLog_.open(rnemdFileName.c_str());
215 >      }
216 >      if (outputAh_) {
217 >        rnemdFileName = "Ah.log";
218 >        AhLog_.open(rnemdFileName.c_str());
219 >      }
220 >      if (outputVz_) {
221 >        rnemdFileName = "velocityZ.log";
222 >        vzzLog_.open(rnemdFileName.c_str());
223 >      }
224 >      logFrameCount_ = 0;
225   #ifdef IS_MPI
226      }
227   #endif
# Line 233 | Line 266 | namespace OpenMD {
266      xyzTempCount_.resize(rnemdLogWidth_, 0);
267      rotTempHist_.resize(rnemdLogWidth_, 0.0);
268      rotTempCount_.resize(rnemdLogWidth_, 0);
269 +    // James put this in
270 +    DenHist_.resize(rnemdLogWidth_, 0.0);
271 +    pzzHist_.resize(rnemdLogWidth_, 0.0);
272  
273      set_RNEMD_exchange_total(0.0);
274      if (simParams->haveRNEMD_targetFlux()) {
# Line 297 | Line 333 | namespace OpenMD {
333        painCave.isFatal = 0;
334        painCave.severity = OPENMD_INFO;
335        simError();
336 <
336 >      
337        if (outputTemp_) tempLog_.close();
338        if (outputVx_)   vxzLog_.close();
339        if (outputVy_)   vyzLog_.close();
# Line 317 | Line 353 | namespace OpenMD {
353          zTempLog_.close();
354        }
355        if (outputRotTemp_) rotTempLog_.close();
356 <
356 >      // James put this in
357 >      if (outputDen_) denLog_.close();
358 >      if (outputAh_)  AhLog_.close();
359 >      if (outputVz_)  vzzLog_.close();
360 >      
361   #ifdef IS_MPI
362      }
363   #endif
# Line 454 | Line 494 | namespace OpenMD {
494          RealType val;
495          int rank;
496        } max_vals, min_vals;
497 <    
497 >      
498        if (my_min_found) {
499          min_vals.val = min_val;
500        } else {
# Line 760 | Line 800 | namespace OpenMD {
800      Kcz *= 0.5;
801      Kcw *= 0.5;
802  
803 <    std::cerr << "Khx= " << Khx << "\tKhy= " << Khy << "\tKhz= " << Khz
804 <              << "\tKhw= " << Khw << "\tKcx= " << Kcx << "\tKcy= " << Kcy
805 <              << "\tKcz= " << Kcz << "\tKcw= " << Kcw << "\n";
806 <    std::cerr << "Phx= " << Phx << "\tPhy= " << Phy << "\tPhz= " << Phz
807 <              << "\tPcx= " << Pcx << "\tPcy= " << Pcy << "\tPcz= " <<Pcz<<"\n";
803 >    // std::cerr << "Khx= " << Khx << "\tKhy= " << Khy << "\tKhz= " << Khz
804 >    //        << "\tKhw= " << Khw << "\tKcx= " << Kcx << "\tKcy= " << Kcy
805 >    //        << "\tKcz= " << Kcz << "\tKcw= " << Kcw << "\n";
806 >    // std::cerr << "Phx= " << Phx << "\tPhy= " << Phy << "\tPhz= " << Phz
807 >    //        << "\tPcx= " << Pcx << "\tPcy= " << Pcy << "\tPcz= " <<Pcz<<"\n";
808  
809   #ifdef IS_MPI
810      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Phx, 1, MPI::REALTYPE, MPI::SUM);
# Line 1105 | Line 1145 | namespace OpenMD {
1145    void RNEMD::doShiftScale() {
1146  
1147      Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
1148 +    RealType time = currentSnap_->getTime();    
1149      Mat3x3d hmat = currentSnap_->getHmat();
1150  
1151      seleMan_.setSelectionSet(evaluator_.evaluate());
# Line 1121 | Line 1162 | namespace OpenMD {
1162      Vector3d Pc(V3Zero);
1163      RealType Mc = 0.0;
1164      RealType Kc = 0.0;
1165 +    
1166  
1167      for (sd = seleMan_.beginSelected(selei); sd != NULL;
1168           sd = seleMan_.nextSelected(selei)) {
# Line 1198 | Line 1240 | namespace OpenMD {
1240      Kh *= 0.5;
1241      Kc *= 0.5;
1242  
1243 <    std::cerr << "Mh= " << Mh << "\tKh= " << Kh << "\tMc= " << Mc
1244 <              << "\tKc= " << Kc << endl;
1245 <    std::cerr << "Ph= " << Ph << "\tPc= " << Pc << endl;
1246 <
1243 >    // std::cerr << "Mh= " << Mh << "\tKh= " << Kh << "\tMc= " << Mc
1244 >    //        << "\tKc= " << Kc << endl;
1245 >    // std::cerr << "Ph= " << Ph << "\tPc= " << Pc << endl;
1246 >    
1247   #ifdef IS_MPI
1248      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Ph[0], 3, MPI::REALTYPE, MPI::SUM);
1249      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Pc[0], 3, MPI::REALTYPE, MPI::SUM);
# Line 1215 | Line 1257 | namespace OpenMD {
1257      if ((Mh > 0.0) && (Mc > 0.0)) {//both slabs are not empty
1258        Vector3d vc = Pc / Mc;
1259        Vector3d ac = njzp_ / Mc + vc;
1260 +      Vector3d acrec = njzp_ / Mc;
1261        RealType cNumerator = Kc - targetJzKE_ - 0.5 * Mc * ac.lengthSquare();
1262        if (cNumerator > 0.0) {
1263          RealType cDenominator = Kc - 0.5 * Mc * vc.lengthSquare();
# Line 1223 | Line 1266 | namespace OpenMD {
1266            if ((c > 0.9) && (c < 1.1)) {//restrict scaling coefficients
1267              Vector3d vh = Ph / Mh;
1268              Vector3d ah = jzp_ / Mh + vh;
1269 +            Vector3d ahrec = jzp_ / Mh;
1270              RealType hNumerator = Kh + targetJzKE_
1271                - 0.5 * Mh * ah.lengthSquare();
1272              if (hNumerator > 0.0) {
# Line 1230 | Line 1274 | namespace OpenMD {
1274                if (hDenominator > 0.0) {
1275                  RealType h = sqrt(hNumerator / hDenominator);
1276                  if ((h > 0.9) && (h < 1.1)) {
1277 <                  std::cerr << "cold slab scaling coefficient: " << c << "\n";
1278 <                  std::cerr << "hot slab scaling coefficient: " << h << "\n";
1277 >                  // std::cerr << "cold slab scaling coefficient: " << c << "\n";
1278 >                  // std::cerr << "hot slab scaling coefficient: " << h <<  "\n";
1279                    vector<StuntDouble*>::iterator sdi;
1280                    Vector3d vel;
1281                    for (sdi = coldBin.begin(); sdi != coldBin.end(); sdi++) {
# Line 1261 | Line 1305 | namespace OpenMD {
1305                    // this is a redundant variable for doShiftScale() so that
1306                    // RNEMD can output one exchange quantity needed in a job.
1307                    // need a better way to do this.
1308 +                  //cerr << "acx =" << ac.x() << "ahx =" << ah.x() << '\n';
1309 +                  //cerr << "acy =" << ac.y() << "ahy =" << ah.y() << '\n';
1310 +                  //cerr << "acz =" << ac.z() << "ahz =" << ah.z() << '\n';
1311 +                  Asum_ += (ahrec.z() - acrec.z());
1312 +                  Jsum_ += (jzp_.z()*((1/Mh)+(1/Mc)));
1313 +                  AhCount_ = ahrec.z();
1314 +                  if (outputAh_) {
1315 +                    AhLog_ << time << "   ";
1316 +                    AhLog_ << AhCount_;
1317 +                    AhLog_ << endl;
1318 +                  }              
1319                  }
1320                }
1321              }
# Line 1269 | Line 1324 | namespace OpenMD {
1324        }
1325      }
1326      if (successfulExchange != true) {
1327 <      sprintf(painCave.errMsg,
1328 <              "RNEMD: exchange NOT performed!\n");
1329 <      painCave.isFatal = 0;
1330 <      painCave.severity = OPENMD_INFO;
1331 <      simError();        
1327 >      //   sprintf(painCave.errMsg,
1328 >      //              "RNEMD: exchange NOT performed!\n");
1329 >      //   painCave.isFatal = 0;
1330 >      //   painCave.severity = OPENMD_INFO;
1331 >      //   simError();        
1332        failTrialCount_++;
1333      }
1334    }
# Line 1316 | Line 1371 | namespace OpenMD {
1371      StuntDouble* sd;
1372      int idx;
1373  
1374 +    logFrameCount_++;
1375 +
1376      // alternative approach, track all molecules instead of only those
1377      // selected for scaling/swapping:
1378      /*
# Line 1324 | Line 1381 | namespace OpenMD {
1381      Molecule* mol;
1382      StuntDouble* integrableObject;
1383      for (mol = info_->beginMolecule(miter); mol != NULL;
1384 <         mol = info_->nextMolecule(miter))
1384 >      mol = info_->nextMolecule(miter))
1385        integrableObject is essentially sd
1386          for (integrableObject = mol->beginIntegrableObject(iiter);
1387               integrableObject != NULL;
# Line 1427 | Line 1484 | namespace OpenMD {
1484            / PhysicalConstants::kb;//may move to getStatus()
1485          rotTempHist_[binNo] += value;
1486        }
1487 <
1487 >      // James put this in.
1488 >      if (outputDen_) {
1489 >        //value = 1.0;
1490 >        DenHist_[binNo] += 1;
1491 >      }
1492 >      if (outputVz_) {
1493 >        value = mass * vel[2];
1494 >        //vyzCount_[binNo]++;
1495 >        pzzHist_[binNo] += value;
1496 >      }    
1497      }
1498    }
1499  
# Line 1493 | Line 1559 | namespace OpenMD {
1559        MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &rotTempCount_[0],
1560                                  rnemdLogWidth_, MPI::INT, MPI::SUM);
1561      }
1562 <
1562 >    // James put this in
1563 >    if (outputDen_) {
1564 >      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &DenHist_[0],
1565 >                                rnemdLogWidth_, MPI::REALTYPE, MPI::SUM);
1566 >    }
1567 >    if (outputAh_) {
1568 >      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &AhCount_,
1569 >                                1, MPI::REALTYPE, MPI::SUM);
1570 >    }
1571 >    if (outputVz_) {
1572 >      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &pzzHist_[0],
1573 >                                rnemdLogWidth_, MPI::REALTYPE, MPI::SUM);
1574 >    }
1575 >    
1576      // If we're the root node, should we print out the results
1577      int worldRank = MPI::COMM_WORLD.Get_rank();
1578      if (worldRank == 0) {
# Line 1550 | Line 1629 | namespace OpenMD {
1629          }
1630          rotTempLog_ << endl;
1631        }
1632 <
1632 >      // James put this in.
1633 >      Mat3x3d hmat = currentSnap_->getHmat();
1634 >      if (outputDen_) {
1635 >        denLog_ << time;
1636 >        for (j = 0; j < rnemdLogWidth_; j++) {
1637 >          
1638 >          RealType binVol = hmat(0,0) * hmat(1,1) * (hmat(2,2) / float(nBins_));
1639 >          denLog_ << "\t" << DenHist_[j] / (float(logFrameCount_) * binVol);
1640 >        }
1641 >        denLog_ << endl;
1642 >      }
1643 >      if (outputVz_) {
1644 >        vzzLog_ << time;
1645 >        for (j = 0; j < rnemdLogWidth_; j++) {
1646 >          vzzLog_ << "\t" << pzzHist_[j] / mHist_[j];
1647 >        }
1648 >        vzzLog_ << endl;
1649 >      }      
1650   #ifdef IS_MPI
1651      }
1652   #endif
# Line 1586 | Line 1682 | namespace OpenMD {
1682          rotTempCount_[j] = 0;
1683          rotTempHist_[j] = 0.0;
1684        }
1685 +    // James put this in
1686 +    if (outputDen_)
1687 +      for (j = 0; j < rnemdLogWidth_; j++) {
1688 +        //pyzCount_[j] = 0;
1689 +        DenHist_[j] = 0.0;
1690 +      }
1691 +    if (outputVz_)
1692 +      for (j = 0; j < rnemdLogWidth_; j++) {
1693 +        //pyzCount_[j] = 0;
1694 +        pzzHist_[j] = 0.0;
1695 +      }    
1696 +     // reset the counter
1697 +    
1698 +    Numcount_++;
1699 +    if (Numcount_ > int(runTime_/statusTime_))
1700 +      cerr << "time =" << time << "  Asum =" << Asum_ << '\n';
1701 +    if (Numcount_ > int(runTime_/statusTime_))
1702 +      cerr << "time =" << time << "  Jsum =" << Jsum_ << '\n';
1703 +    
1704 +    logFrameCount_ = 0;
1705    }
1706   }
1707  

Diff Legend

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