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.
branches/development/src/rnemd/RNEMD.cpp (file contents), Revision 1764 by gezelter, Tue Jul 3 18:32:27 2012 UTC

# Line 40 | Line 40
40   */
41  
42   #include <cmath>
43 < #include "integrators/RNEMD.hpp"
43 > #include "rnemd/RNEMD.hpp"
44   #include "math/Vector3.hpp"
45   #include "math/Vector.hpp"
46   #include "math/SquareMatrix3.hpp"
# Line 70 | Line 70 | namespace OpenMD {
70  
71      int seedValue;
72      Globals * simParams = info->getSimParams();
73 +    RNEMDParameters* rnemdParams = simParams->getRNEMDParameters();
74  
75      stringToEnumMap_["KineticSwap"] = rnemdKineticSwap;
76      stringToEnumMap_["KineticScale"] = rnemdKineticScale;
# Line 85 | Line 86 | namespace OpenMD {
86      stringToEnumMap_["ShiftScaleVAM"] = rnemdShiftScaleVAM;
87      stringToEnumMap_["Unknown"] = rnemdUnknown;
88  
89 <    rnemdObjectSelection_ = simParams->getRNEMD_objectSelection();
89 >    runTime_ = simParams->getRunTime();
90 >    statusTime_ = simParams->getStatusTime();
91 >
92 >    rnemdObjectSelection_ = rnemdParams->getObjectSelection();
93      evaluator_.loadScriptString(rnemdObjectSelection_);
94      seleMan_.setSelectionSet(evaluator_.evaluate());
95  
# Line 110 | Line 114 | namespace OpenMD {
114        simError();
115      }
116      
117 <    const string st = simParams->getRNEMD_exchangeType();
117 >    const string st = rnemdParams->getExchangeType();
118  
119      map<string, RNEMDTypeEnum>::iterator i;
120      i = stringToEnumMap_.find(st);
# Line 127 | Line 131 | namespace OpenMD {
131      }
132      
133      outputTemp_ = false;
134 <    if (simParams->haveRNEMD_outputTemperature()) {
135 <      outputTemp_ = simParams->getRNEMD_outputTemperature();
134 >    if (rnemdParams->haveOutputTemperature()) {
135 >      outputTemp_ = rnemdParams->getOutputTemperature();
136      } else if ((rnemdType_ == rnemdKineticSwap) ||
137                 (rnemdType_ == rnemdKineticScale) ||
138                 (rnemdType_ == rnemdKineticScaleVAM) ||
# Line 136 | Line 140 | namespace OpenMD {
140        outputTemp_ = true;
141      }
142      outputVx_ = false;
143 <    if (simParams->haveRNEMD_outputVx()) {
144 <      outputVx_ = simParams->getRNEMD_outputVx();
143 >    if (rnemdParams->haveOutputVx()) {
144 >      outputVx_ = rnemdParams->getOutputVx();
145      } else if ((rnemdType_ == rnemdPx) || (rnemdType_ == rnemdPxScale)) {
146        outputVx_ = true;
147      }
148      outputVy_ = false;
149 <    if (simParams->haveRNEMD_outputVy()) {
150 <      outputVy_ = simParams->getRNEMD_outputVy();
149 >    if (rnemdParams->haveOutputVy()) {
150 >      outputVy_ = rnemdParams->getOutputVy();
151      } else if ((rnemdType_ == rnemdPy) || (rnemdType_ == rnemdPyScale)) {
152        outputVy_ = true;
153      }
154      output3DTemp_ = false;
155 <    if (simParams->haveRNEMD_outputXyzTemperature()) {
156 <      output3DTemp_ = simParams->getRNEMD_outputXyzTemperature();
155 >    if (rnemdParams->haveOutputXyzTemperature()) {
156 >      output3DTemp_ = rnemdParams->getOutputXyzTemperature();
157      }
158      outputRotTemp_ = false;
159 <    if (simParams->haveRNEMD_outputRotTemperature()) {
160 <      outputRotTemp_ = simParams->getRNEMD_outputRotTemperature();
159 >    if (rnemdParams->haveOutputRotTemperature()) {
160 >      outputRotTemp_ = rnemdParams->getOutputRotTemperature();
161      }
162 +    // James put this in.
163 +    outputDen_ = false;
164 +    if (rnemdParams->haveOutputDen()) {
165 +      outputDen_ = rnemdParams->getOutputDen();
166 +    }
167 +    outputAh_ = false;
168 +    if (rnemdParams->haveOutputAh()) {
169 +      outputAh_ = rnemdParams->getOutputAh();
170 +    }    
171 +    outputVz_ = false;
172 +    if (rnemdParams->haveOutputVz()) {
173 +      outputVz_ = rnemdParams->getOutputVz();
174 +    } else if ((rnemdType_ == rnemdPz) || (rnemdType_ == rnemdPzScale)) {
175 +      outputVz_ = true;
176 +    }
177 +    
178  
179   #ifdef IS_MPI
180      if (worldRank == 0) {
# Line 188 | Line 208 | namespace OpenMD {
208          rnemdFileName = "temperatureR.log";
209          rotTempLog_.open(rnemdFileName.c_str());
210        }
211 <
211 >      
212 >      //James put this in
213 >      if (outputDen_) {
214 >        rnemdFileName = "Density.log";
215 >        denLog_.open(rnemdFileName.c_str());
216 >      }
217 >      if (outputAh_) {
218 >        rnemdFileName = "Ah.log";
219 >        AhLog_.open(rnemdFileName.c_str());
220 >      }
221 >      if (outputVz_) {
222 >        rnemdFileName = "velocityZ.log";
223 >        vzzLog_.open(rnemdFileName.c_str());
224 >      }
225 >      logFrameCount_ = 0;
226   #ifdef IS_MPI
227      }
228   #endif
229  
230 <    set_RNEMD_exchange_time(simParams->getRNEMD_exchangeTime());
231 <    set_RNEMD_nBins(simParams->getRNEMD_nBins());
230 >    set_RNEMD_exchange_time(rnemdParams->getExchangeTime());
231 >    set_RNEMD_nBins(rnemdParams->getNbins());
232      midBin_ = nBins_ / 2;
233 <    if (simParams->haveRNEMD_binShift()) {
234 <      if (simParams->getRNEMD_binShift()) {
233 >    if (rnemdParams->haveBinShift()) {
234 >      if (rnemdParams->getBinShift()) {
235          zShift_ = 0.5 / (RealType)(nBins_);
236        } else {
237          zShift_ = 0.0;
# Line 208 | Line 242 | namespace OpenMD {
242      //cerr << "I shift slabs by " << zShift_ << " Lz\n";
243      //shift slabs by half slab width, maybe useful in heterogeneous systems
244      //set to 0.0 if not using it; N/A in status output yet
245 <    if (simParams->haveRNEMD_logWidth()) {
246 <      set_RNEMD_logWidth(simParams->getRNEMD_logWidth());
245 >    if (rnemdParams->haveLogWidth()) {
246 >      set_RNEMD_logWidth(rnemdParams->getLogWidth());
247        /*arbitary rnemdLogWidth_, no checking;
248        if (rnemdLogWidth_ != nBins_ && rnemdLogWidth_ != midBin_ + 1) {
249          cerr << "WARNING! RNEMD_logWidth has abnormal value!\n";
# Line 233 | Line 267 | namespace OpenMD {
267      xyzTempCount_.resize(rnemdLogWidth_, 0);
268      rotTempHist_.resize(rnemdLogWidth_, 0.0);
269      rotTempCount_.resize(rnemdLogWidth_, 0);
270 +    // James put this in
271 +    DenHist_.resize(rnemdLogWidth_, 0.0);
272 +    pzzHist_.resize(rnemdLogWidth_, 0.0);
273  
274      set_RNEMD_exchange_total(0.0);
275 <    if (simParams->haveRNEMD_targetFlux()) {
276 <      set_RNEMD_target_flux(simParams->getRNEMD_targetFlux());
275 >    if (rnemdParams->haveTargetFlux()) {
276 >      set_RNEMD_target_flux(rnemdParams->getTargetFlux());
277      } else {
278        set_RNEMD_target_flux(0.0);
279      }
280 <    if (simParams->haveRNEMD_targetJzKE()) {
281 <      set_RNEMD_target_JzKE(simParams->getRNEMD_targetJzKE());
280 >    if (rnemdParams->haveTargetJzKE()) {
281 >      set_RNEMD_target_JzKE(rnemdParams->getTargetJzKE());
282      } else {
283        set_RNEMD_target_JzKE(0.0);
284      }
285 <    if (simParams->haveRNEMD_targetJzpx()) {
286 <      set_RNEMD_target_jzpx(simParams->getRNEMD_targetJzpx());
285 >    if (rnemdParams->haveTargetJzpx()) {
286 >      set_RNEMD_target_jzpx(rnemdParams->getTargetJzpx());
287      } else {
288        set_RNEMD_target_jzpx(0.0);
289      }
290      jzp_.x() = targetJzpx_;
291      njzp_.x() = -targetJzpx_;
292 <    if (simParams->haveRNEMD_targetJzpy()) {
293 <      set_RNEMD_target_jzpy(simParams->getRNEMD_targetJzpy());
292 >    if (rnemdParams->haveTargetJzpy()) {
293 >      set_RNEMD_target_jzpy(rnemdParams->getTargetJzpy());
294      } else {
295        set_RNEMD_target_jzpy(0.0);
296      }
297      jzp_.y() = targetJzpy_;
298      njzp_.y() = -targetJzpy_;
299 <    if (simParams->haveRNEMD_targetJzpz()) {
300 <      set_RNEMD_target_jzpz(simParams->getRNEMD_targetJzpz());
299 >    if (rnemdParams->haveTargetJzpz()) {
300 >      set_RNEMD_target_jzpz(rnemdParams->getTargetJzpz());
301      } else {
302        set_RNEMD_target_jzpz(0.0);
303      }
# Line 297 | Line 334 | namespace OpenMD {
334        painCave.isFatal = 0;
335        painCave.severity = OPENMD_INFO;
336        simError();
337 <
337 >      
338        if (outputTemp_) tempLog_.close();
339        if (outputVx_)   vxzLog_.close();
340        if (outputVy_)   vyzLog_.close();
# Line 317 | Line 354 | namespace OpenMD {
354          zTempLog_.close();
355        }
356        if (outputRotTemp_) rotTempLog_.close();
357 <
357 >      // James put this in
358 >      if (outputDen_) denLog_.close();
359 >      if (outputAh_)  AhLog_.close();
360 >      if (outputVz_)  vzzLog_.close();
361 >      
362   #ifdef IS_MPI
363      }
364   #endif
# Line 454 | Line 495 | namespace OpenMD {
495          RealType val;
496          int rank;
497        } max_vals, min_vals;
498 <    
498 >      
499        if (my_min_found) {
500          min_vals.val = min_val;
501        } else {
# Line 760 | Line 801 | namespace OpenMD {
801      Kcz *= 0.5;
802      Kcw *= 0.5;
803  
804 <    std::cerr << "Khx= " << Khx << "\tKhy= " << Khy << "\tKhz= " << Khz
805 <              << "\tKhw= " << Khw << "\tKcx= " << Kcx << "\tKcy= " << Kcy
806 <              << "\tKcz= " << Kcz << "\tKcw= " << Kcw << "\n";
807 <    std::cerr << "Phx= " << Phx << "\tPhy= " << Phy << "\tPhz= " << Phz
808 <              << "\tPcx= " << Pcx << "\tPcy= " << Pcy << "\tPcz= " <<Pcz<<"\n";
804 >    // std::cerr << "Khx= " << Khx << "\tKhy= " << Khy << "\tKhz= " << Khz
805 >    //        << "\tKhw= " << Khw << "\tKcx= " << Kcx << "\tKcy= " << Kcy
806 >    //        << "\tKcz= " << Kcz << "\tKcw= " << Kcw << "\n";
807 >    // std::cerr << "Phx= " << Phx << "\tPhy= " << Phy << "\tPhz= " << Phz
808 >    //        << "\tPcx= " << Pcx << "\tPcy= " << Pcy << "\tPcz= " <<Pcz<<"\n";
809  
810   #ifdef IS_MPI
811      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Phx, 1, MPI::REALTYPE, MPI::SUM);
# Line 1105 | Line 1146 | namespace OpenMD {
1146    void RNEMD::doShiftScale() {
1147  
1148      Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
1149 +    RealType time = currentSnap_->getTime();    
1150      Mat3x3d hmat = currentSnap_->getHmat();
1151  
1152      seleMan_.setSelectionSet(evaluator_.evaluate());
# Line 1121 | Line 1163 | namespace OpenMD {
1163      Vector3d Pc(V3Zero);
1164      RealType Mc = 0.0;
1165      RealType Kc = 0.0;
1166 +    
1167  
1168      for (sd = seleMan_.beginSelected(selei); sd != NULL;
1169           sd = seleMan_.nextSelected(selei)) {
# Line 1198 | Line 1241 | namespace OpenMD {
1241      Kh *= 0.5;
1242      Kc *= 0.5;
1243  
1244 <    std::cerr << "Mh= " << Mh << "\tKh= " << Kh << "\tMc= " << Mc
1245 <              << "\tKc= " << Kc << endl;
1246 <    std::cerr << "Ph= " << Ph << "\tPc= " << Pc << endl;
1247 <
1244 >    // std::cerr << "Mh= " << Mh << "\tKh= " << Kh << "\tMc= " << Mc
1245 >    //        << "\tKc= " << Kc << endl;
1246 >    // std::cerr << "Ph= " << Ph << "\tPc= " << Pc << endl;
1247 >    
1248   #ifdef IS_MPI
1249      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Ph[0], 3, MPI::REALTYPE, MPI::SUM);
1250      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Pc[0], 3, MPI::REALTYPE, MPI::SUM);
# Line 1215 | Line 1258 | namespace OpenMD {
1258      if ((Mh > 0.0) && (Mc > 0.0)) {//both slabs are not empty
1259        Vector3d vc = Pc / Mc;
1260        Vector3d ac = njzp_ / Mc + vc;
1261 +      Vector3d acrec = njzp_ / Mc;
1262        RealType cNumerator = Kc - targetJzKE_ - 0.5 * Mc * ac.lengthSquare();
1263        if (cNumerator > 0.0) {
1264          RealType cDenominator = Kc - 0.5 * Mc * vc.lengthSquare();
# Line 1223 | Line 1267 | namespace OpenMD {
1267            if ((c > 0.9) && (c < 1.1)) {//restrict scaling coefficients
1268              Vector3d vh = Ph / Mh;
1269              Vector3d ah = jzp_ / Mh + vh;
1270 +            Vector3d ahrec = jzp_ / Mh;
1271              RealType hNumerator = Kh + targetJzKE_
1272                - 0.5 * Mh * ah.lengthSquare();
1273              if (hNumerator > 0.0) {
# Line 1230 | Line 1275 | namespace OpenMD {
1275                if (hDenominator > 0.0) {
1276                  RealType h = sqrt(hNumerator / hDenominator);
1277                  if ((h > 0.9) && (h < 1.1)) {
1278 <                  std::cerr << "cold slab scaling coefficient: " << c << "\n";
1279 <                  std::cerr << "hot slab scaling coefficient: " << h << "\n";
1278 >                  // std::cerr << "cold slab scaling coefficient: " << c << "\n";
1279 >                  // std::cerr << "hot slab scaling coefficient: " << h <<  "\n";
1280                    vector<StuntDouble*>::iterator sdi;
1281                    Vector3d vel;
1282                    for (sdi = coldBin.begin(); sdi != coldBin.end(); sdi++) {
# Line 1261 | Line 1306 | namespace OpenMD {
1306                    // this is a redundant variable for doShiftScale() so that
1307                    // RNEMD can output one exchange quantity needed in a job.
1308                    // need a better way to do this.
1309 +                  //cerr << "acx =" << ac.x() << "ahx =" << ah.x() << '\n';
1310 +                  //cerr << "acy =" << ac.y() << "ahy =" << ah.y() << '\n';
1311 +                  //cerr << "acz =" << ac.z() << "ahz =" << ah.z() << '\n';
1312 +                  Asum_ += (ahrec.z() - acrec.z());
1313 +                  Jsum_ += (jzp_.z()*((1/Mh)+(1/Mc)));
1314 +                  AhCount_ = ahrec.z();
1315 +                  if (outputAh_) {
1316 +                    AhLog_ << time << "   ";
1317 +                    AhLog_ << AhCount_;
1318 +                    AhLog_ << endl;
1319 +                  }              
1320                  }
1321                }
1322              }
# Line 1269 | Line 1325 | namespace OpenMD {
1325        }
1326      }
1327      if (successfulExchange != true) {
1328 <      sprintf(painCave.errMsg,
1329 <              "RNEMD: exchange NOT performed!\n");
1330 <      painCave.isFatal = 0;
1331 <      painCave.severity = OPENMD_INFO;
1332 <      simError();        
1328 >      //   sprintf(painCave.errMsg,
1329 >      //              "RNEMD: exchange NOT performed!\n");
1330 >      //   painCave.isFatal = 0;
1331 >      //   painCave.severity = OPENMD_INFO;
1332 >      //   simError();        
1333        failTrialCount_++;
1334      }
1335    }
# Line 1316 | Line 1372 | namespace OpenMD {
1372      StuntDouble* sd;
1373      int idx;
1374  
1375 +    logFrameCount_++;
1376 +
1377      // alternative approach, track all molecules instead of only those
1378      // selected for scaling/swapping:
1379      /*
# Line 1324 | Line 1382 | namespace OpenMD {
1382      Molecule* mol;
1383      StuntDouble* integrableObject;
1384      for (mol = info_->beginMolecule(miter); mol != NULL;
1385 <         mol = info_->nextMolecule(miter))
1385 >      mol = info_->nextMolecule(miter))
1386        integrableObject is essentially sd
1387          for (integrableObject = mol->beginIntegrableObject(iiter);
1388               integrableObject != NULL;
# Line 1427 | Line 1485 | namespace OpenMD {
1485            / PhysicalConstants::kb;//may move to getStatus()
1486          rotTempHist_[binNo] += value;
1487        }
1488 <
1488 >      // James put this in.
1489 >      if (outputDen_) {
1490 >        //value = 1.0;
1491 >        DenHist_[binNo] += 1;
1492 >      }
1493 >      if (outputVz_) {
1494 >        value = mass * vel[2];
1495 >        //vyzCount_[binNo]++;
1496 >        pzzHist_[binNo] += value;
1497 >      }    
1498      }
1499    }
1500  
# Line 1445 | Line 1512 | namespace OpenMD {
1512    void RNEMD::getStatus() {
1513  
1514      Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot();
1448 –    Stats& stat = currentSnap_->statData;
1515      RealType time = currentSnap_->getTime();
1450 –
1451 –    stat[Stats::RNEMD_EXCHANGE_TOTAL] = exchangeSum_;
1516      //or to be more meaningful, define another item as exchangeSum_ / time
1517      int j;
1518  
# Line 1493 | Line 1557 | namespace OpenMD {
1557        MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &rotTempCount_[0],
1558                                  rnemdLogWidth_, MPI::INT, MPI::SUM);
1559      }
1560 <
1560 >    // James put this in
1561 >    if (outputDen_) {
1562 >      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &DenHist_[0],
1563 >                                rnemdLogWidth_, MPI::REALTYPE, MPI::SUM);
1564 >    }
1565 >    if (outputAh_) {
1566 >      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &AhCount_,
1567 >                                1, MPI::REALTYPE, MPI::SUM);
1568 >    }
1569 >    if (outputVz_) {
1570 >      MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &pzzHist_[0],
1571 >                                rnemdLogWidth_, MPI::REALTYPE, MPI::SUM);
1572 >    }
1573 >    
1574      // If we're the root node, should we print out the results
1575      int worldRank = MPI::COMM_WORLD.Get_rank();
1576      if (worldRank == 0) {
# Line 1550 | Line 1627 | namespace OpenMD {
1627          }
1628          rotTempLog_ << endl;
1629        }
1630 <
1630 >      // James put this in.
1631 >      Mat3x3d hmat = currentSnap_->getHmat();
1632 >      if (outputDen_) {
1633 >        denLog_ << time;
1634 >        for (j = 0; j < rnemdLogWidth_; j++) {
1635 >          
1636 >          RealType binVol = hmat(0,0) * hmat(1,1) * (hmat(2,2) / float(nBins_));
1637 >          denLog_ << "\t" << DenHist_[j] / (float(logFrameCount_) * binVol);
1638 >        }
1639 >        denLog_ << endl;
1640 >      }
1641 >      if (outputVz_) {
1642 >        vzzLog_ << time;
1643 >        for (j = 0; j < rnemdLogWidth_; j++) {
1644 >          vzzLog_ << "\t" << pzzHist_[j] / mHist_[j];
1645 >        }
1646 >        vzzLog_ << endl;
1647 >      }      
1648   #ifdef IS_MPI
1649      }
1650   #endif
# Line 1586 | Line 1680 | namespace OpenMD {
1680          rotTempCount_[j] = 0;
1681          rotTempHist_[j] = 0.0;
1682        }
1683 +    // James put this in
1684 +    if (outputDen_)
1685 +      for (j = 0; j < rnemdLogWidth_; j++) {
1686 +        //pyzCount_[j] = 0;
1687 +        DenHist_[j] = 0.0;
1688 +      }
1689 +    if (outputVz_)
1690 +      for (j = 0; j < rnemdLogWidth_; j++) {
1691 +        //pyzCount_[j] = 0;
1692 +        pzzHist_[j] = 0.0;
1693 +      }    
1694 +     // reset the counter
1695 +    
1696 +    Numcount_++;
1697 +    if (Numcount_ > int(runTime_/statusTime_))
1698 +      cerr << "time =" << time << "  Asum =" << Asum_ << '\n';
1699 +    if (Numcount_ > int(runTime_/statusTime_))
1700 +      cerr << "time =" << time << "  Jsum =" << Jsum_ << '\n';
1701 +    
1702 +    logFrameCount_ = 0;
1703    }
1704   }
1705  

Diff Legend

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