--- branches/development/src/integrators/RNEMD.cpp 2012/05/24 20:59:54 1723 +++ branches/development/src/rnemd/RNEMD.cpp 2012/07/09 14:15:52 1769 @@ -40,7 +40,7 @@ */ #include -#include "integrators/RNEMD.hpp" +#include "rnemd/RNEMD.hpp" #include "math/Vector3.hpp" #include "math/Vector.hpp" #include "math/SquareMatrix3.hpp" @@ -70,6 +70,7 @@ namespace OpenMD { int seedValue; Globals * simParams = info->getSimParams(); + RNEMDParameters* rnemdParams = simParams->getRNEMDParameters(); stringToEnumMap_["KineticSwap"] = rnemdKineticSwap; stringToEnumMap_["KineticScale"] = rnemdKineticScale; @@ -85,7 +86,10 @@ namespace OpenMD { stringToEnumMap_["ShiftScaleVAM"] = rnemdShiftScaleVAM; stringToEnumMap_["Unknown"] = rnemdUnknown; - rnemdObjectSelection_ = simParams->getRNEMD_objectSelection(); + runTime_ = simParams->getRunTime(); + statusTime_ = simParams->getStatusTime(); + + rnemdObjectSelection_ = rnemdParams->getObjectSelection(); evaluator_.loadScriptString(rnemdObjectSelection_); seleMan_.setSelectionSet(evaluator_.evaluate()); @@ -110,7 +114,7 @@ namespace OpenMD { simError(); } - const string st = simParams->getRNEMD_exchangeType(); + const string st = rnemdParams->getExchangeType(); map::iterator i; i = stringToEnumMap_.find(st); @@ -127,8 +131,8 @@ namespace OpenMD { } outputTemp_ = false; - if (simParams->haveRNEMD_outputTemperature()) { - outputTemp_ = simParams->getRNEMD_outputTemperature(); + if (rnemdParams->haveOutputTemperature()) { + outputTemp_ = rnemdParams->getOutputTemperature(); } else if ((rnemdType_ == rnemdKineticSwap) || (rnemdType_ == rnemdKineticScale) || (rnemdType_ == rnemdKineticScaleVAM) || @@ -136,25 +140,41 @@ namespace OpenMD { outputTemp_ = true; } outputVx_ = false; - if (simParams->haveRNEMD_outputVx()) { - outputVx_ = simParams->getRNEMD_outputVx(); + if (rnemdParams->haveOutputVx()) { + outputVx_ = rnemdParams->getOutputVx(); } else if ((rnemdType_ == rnemdPx) || (rnemdType_ == rnemdPxScale)) { outputVx_ = true; } outputVy_ = false; - if (simParams->haveRNEMD_outputVy()) { - outputVy_ = simParams->getRNEMD_outputVy(); + if (rnemdParams->haveOutputVy()) { + outputVy_ = rnemdParams->getOutputVy(); } else if ((rnemdType_ == rnemdPy) || (rnemdType_ == rnemdPyScale)) { outputVy_ = true; } output3DTemp_ = false; - if (simParams->haveRNEMD_outputXyzTemperature()) { - output3DTemp_ = simParams->getRNEMD_outputXyzTemperature(); + if (rnemdParams->haveOutputXyzTemperature()) { + output3DTemp_ = rnemdParams->getOutputXyzTemperature(); } outputRotTemp_ = false; - if (simParams->haveRNEMD_outputRotTemperature()) { - outputRotTemp_ = simParams->getRNEMD_outputRotTemperature(); + if (rnemdParams->haveOutputRotTemperature()) { + outputRotTemp_ = rnemdParams->getOutputRotTemperature(); } + // James put this in. + outputDen_ = false; + if (rnemdParams->haveOutputDen()) { + outputDen_ = rnemdParams->getOutputDen(); + } + outputAh_ = false; + if (rnemdParams->haveOutputAh()) { + outputAh_ = rnemdParams->getOutputAh(); + } + outputVz_ = false; + if (rnemdParams->haveOutputVz()) { + outputVz_ = rnemdParams->getOutputVz(); + } else if ((rnemdType_ == rnemdPz) || (rnemdType_ == rnemdPzScale)) { + outputVz_ = true; + } + #ifdef IS_MPI if (worldRank == 0) { @@ -188,16 +208,30 @@ namespace OpenMD { rnemdFileName = "temperatureR.log"; rotTempLog_.open(rnemdFileName.c_str()); } - + + //James put this in + if (outputDen_) { + rnemdFileName = "Density.log"; + denLog_.open(rnemdFileName.c_str()); + } + if (outputAh_) { + rnemdFileName = "Ah.log"; + AhLog_.open(rnemdFileName.c_str()); + } + if (outputVz_) { + rnemdFileName = "velocityZ.log"; + vzzLog_.open(rnemdFileName.c_str()); + } + logFrameCount_ = 0; #ifdef IS_MPI } #endif - set_RNEMD_exchange_time(simParams->getRNEMD_exchangeTime()); - set_RNEMD_nBins(simParams->getRNEMD_nBins()); + set_RNEMD_exchange_time(rnemdParams->getExchangeTime()); + set_RNEMD_nBins(rnemdParams->getNbins()); midBin_ = nBins_ / 2; - if (simParams->haveRNEMD_binShift()) { - if (simParams->getRNEMD_binShift()) { + if (rnemdParams->haveBinShift()) { + if (rnemdParams->getBinShift()) { zShift_ = 0.5 / (RealType)(nBins_); } else { zShift_ = 0.0; @@ -208,8 +242,8 @@ namespace OpenMD { //cerr << "I shift slabs by " << zShift_ << " Lz\n"; //shift slabs by half slab width, maybe useful in heterogeneous systems //set to 0.0 if not using it; N/A in status output yet - if (simParams->haveRNEMD_logWidth()) { - set_RNEMD_logWidth(simParams->getRNEMD_logWidth()); + if (rnemdParams->haveLogWidth()) { + set_RNEMD_logWidth(rnemdParams->getLogWidth()); /*arbitary rnemdLogWidth_, no checking; if (rnemdLogWidth_ != nBins_ && rnemdLogWidth_ != midBin_ + 1) { cerr << "WARNING! RNEMD_logWidth has abnormal value!\n"; @@ -233,34 +267,37 @@ namespace OpenMD { xyzTempCount_.resize(rnemdLogWidth_, 0); rotTempHist_.resize(rnemdLogWidth_, 0.0); rotTempCount_.resize(rnemdLogWidth_, 0); + // James put this in + DenHist_.resize(rnemdLogWidth_, 0.0); + pzzHist_.resize(rnemdLogWidth_, 0.0); set_RNEMD_exchange_total(0.0); - if (simParams->haveRNEMD_targetFlux()) { - set_RNEMD_target_flux(simParams->getRNEMD_targetFlux()); + if (rnemdParams->haveTargetFlux()) { + set_RNEMD_target_flux(rnemdParams->getTargetFlux()); } else { set_RNEMD_target_flux(0.0); } - if (simParams->haveRNEMD_targetJzKE()) { - set_RNEMD_target_JzKE(simParams->getRNEMD_targetJzKE()); + if (rnemdParams->haveTargetJzKE()) { + set_RNEMD_target_JzKE(rnemdParams->getTargetJzKE()); } else { set_RNEMD_target_JzKE(0.0); } - if (simParams->haveRNEMD_targetJzpx()) { - set_RNEMD_target_jzpx(simParams->getRNEMD_targetJzpx()); + if (rnemdParams->haveTargetJzpx()) { + set_RNEMD_target_jzpx(rnemdParams->getTargetJzpx()); } else { set_RNEMD_target_jzpx(0.0); } jzp_.x() = targetJzpx_; njzp_.x() = -targetJzpx_; - if (simParams->haveRNEMD_targetJzpy()) { - set_RNEMD_target_jzpy(simParams->getRNEMD_targetJzpy()); + if (rnemdParams->haveTargetJzpy()) { + set_RNEMD_target_jzpy(rnemdParams->getTargetJzpy()); } else { set_RNEMD_target_jzpy(0.0); } jzp_.y() = targetJzpy_; njzp_.y() = -targetJzpy_; - if (simParams->haveRNEMD_targetJzpz()) { - set_RNEMD_target_jzpz(simParams->getRNEMD_targetJzpz()); + if (rnemdParams->haveTargetJzpz()) { + set_RNEMD_target_jzpz(rnemdParams->getTargetJzpz()); } else { set_RNEMD_target_jzpz(0.0); } @@ -297,7 +334,7 @@ namespace OpenMD { painCave.isFatal = 0; painCave.severity = OPENMD_INFO; simError(); - + if (outputTemp_) tempLog_.close(); if (outputVx_) vxzLog_.close(); if (outputVy_) vyzLog_.close(); @@ -317,7 +354,11 @@ namespace OpenMD { zTempLog_.close(); } if (outputRotTemp_) rotTempLog_.close(); - + // James put this in + if (outputDen_) denLog_.close(); + if (outputAh_) AhLog_.close(); + if (outputVz_) vzzLog_.close(); + #ifdef IS_MPI } #endif @@ -454,7 +495,7 @@ namespace OpenMD { RealType val; int rank; } max_vals, min_vals; - + if (my_min_found) { min_vals.val = min_val; } else { @@ -760,11 +801,11 @@ namespace OpenMD { Kcz *= 0.5; Kcw *= 0.5; - std::cerr << "Khx= " << Khx << "\tKhy= " << Khy << "\tKhz= " << Khz - << "\tKhw= " << Khw << "\tKcx= " << Kcx << "\tKcy= " << Kcy - << "\tKcz= " << Kcz << "\tKcw= " << Kcw << "\n"; - std::cerr << "Phx= " << Phx << "\tPhy= " << Phy << "\tPhz= " << Phz - << "\tPcx= " << Pcx << "\tPcy= " << Pcy << "\tPcz= " <getSnapshotManager()->getCurrentSnapshot(); + RealType time = currentSnap_->getTime(); Mat3x3d hmat = currentSnap_->getHmat(); seleMan_.setSelectionSet(evaluator_.evaluate()); @@ -1121,6 +1163,7 @@ namespace OpenMD { Vector3d Pc(V3Zero); RealType Mc = 0.0; RealType Kc = 0.0; + for (sd = seleMan_.beginSelected(selei); sd != NULL; sd = seleMan_.nextSelected(selei)) { @@ -1198,10 +1241,10 @@ namespace OpenMD { Kh *= 0.5; Kc *= 0.5; - std::cerr << "Mh= " << Mh << "\tKh= " << Kh << "\tMc= " << Mc - << "\tKc= " << Kc << endl; - std::cerr << "Ph= " << Ph << "\tPc= " << Pc << endl; - + // std::cerr << "Mh= " << Mh << "\tKh= " << Kh << "\tMc= " << Mc + // << "\tKc= " << Kc << endl; + // std::cerr << "Ph= " << Ph << "\tPc= " << Pc << endl; + #ifdef IS_MPI MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Ph[0], 3, MPI::REALTYPE, MPI::SUM); MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Pc[0], 3, MPI::REALTYPE, MPI::SUM); @@ -1215,6 +1258,7 @@ namespace OpenMD { if ((Mh > 0.0) && (Mc > 0.0)) {//both slabs are not empty Vector3d vc = Pc / Mc; Vector3d ac = njzp_ / Mc + vc; + Vector3d acrec = njzp_ / Mc; RealType cNumerator = Kc - targetJzKE_ - 0.5 * Mc * ac.lengthSquare(); if (cNumerator > 0.0) { RealType cDenominator = Kc - 0.5 * Mc * vc.lengthSquare(); @@ -1223,6 +1267,7 @@ namespace OpenMD { if ((c > 0.9) && (c < 1.1)) {//restrict scaling coefficients Vector3d vh = Ph / Mh; Vector3d ah = jzp_ / Mh + vh; + Vector3d ahrec = jzp_ / Mh; RealType hNumerator = Kh + targetJzKE_ - 0.5 * Mh * ah.lengthSquare(); if (hNumerator > 0.0) { @@ -1230,8 +1275,8 @@ namespace OpenMD { if (hDenominator > 0.0) { RealType h = sqrt(hNumerator / hDenominator); if ((h > 0.9) && (h < 1.1)) { - std::cerr << "cold slab scaling coefficient: " << c << "\n"; - std::cerr << "hot slab scaling coefficient: " << h << "\n"; + // std::cerr << "cold slab scaling coefficient: " << c << "\n"; + // std::cerr << "hot slab scaling coefficient: " << h << "\n"; vector::iterator sdi; Vector3d vel; for (sdi = coldBin.begin(); sdi != coldBin.end(); sdi++) { @@ -1261,6 +1306,17 @@ namespace OpenMD { // this is a redundant variable for doShiftScale() so that // RNEMD can output one exchange quantity needed in a job. // need a better way to do this. + //cerr << "acx =" << ac.x() << "ahx =" << ah.x() << '\n'; + //cerr << "acy =" << ac.y() << "ahy =" << ah.y() << '\n'; + //cerr << "acz =" << ac.z() << "ahz =" << ah.z() << '\n'; + Asum_ += (ahrec.z() - acrec.z()); + Jsum_ += (jzp_.z()*((1/Mh)+(1/Mc))); + AhCount_ = ahrec.z(); + if (outputAh_) { + AhLog_ << time << " "; + AhLog_ << AhCount_; + AhLog_ << endl; + } } } } @@ -1269,11 +1325,11 @@ namespace OpenMD { } } if (successfulExchange != true) { - sprintf(painCave.errMsg, - "RNEMD: exchange NOT performed!\n"); - painCave.isFatal = 0; - painCave.severity = OPENMD_INFO; - simError(); + // sprintf(painCave.errMsg, + // "RNEMD: exchange NOT performed!\n"); + // painCave.isFatal = 0; + // painCave.severity = OPENMD_INFO; + // simError(); failTrialCount_++; } } @@ -1316,19 +1372,21 @@ namespace OpenMD { StuntDouble* sd; int idx; + logFrameCount_++; + // alternative approach, track all molecules instead of only those // selected for scaling/swapping: /* SimInfo::MoleculeIterator miter; vector::iterator iiter; Molecule* mol; - StuntDouble* integrableObject; + StuntDouble* sd; for (mol = info_->beginMolecule(miter); mol != NULL; - mol = info_->nextMolecule(miter)) - integrableObject is essentially sd - for (integrableObject = mol->beginIntegrableObject(iiter); - integrableObject != NULL; - integrableObject = mol->nextIntegrableObject(iiter)) + mol = info_->nextMolecule(miter)) + sd is essentially sd + for (sd = mol->beginIntegrableObject(iiter); + sd != NULL; + sd = mol->nextIntegrableObject(iiter)) */ for (sd = seleMan_.beginSelected(selei); sd != NULL; sd = seleMan_.nextSelected(selei)) { @@ -1427,7 +1485,16 @@ namespace OpenMD { / PhysicalConstants::kb;//may move to getStatus() rotTempHist_[binNo] += value; } - + // James put this in. + if (outputDen_) { + //value = 1.0; + DenHist_[binNo] += 1; + } + if (outputVz_) { + value = mass * vel[2]; + //vyzCount_[binNo]++; + pzzHist_[binNo] += value; + } } } @@ -1445,10 +1512,7 @@ namespace OpenMD { void RNEMD::getStatus() { Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot(); - Stats& stat = currentSnap_->statData; RealType time = currentSnap_->getTime(); - - stat[Stats::RNEMD_EXCHANGE_TOTAL] = exchangeSum_; //or to be more meaningful, define another item as exchangeSum_ / time int j; @@ -1493,7 +1557,20 @@ namespace OpenMD { MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &rotTempCount_[0], rnemdLogWidth_, MPI::INT, MPI::SUM); } - + // James put this in + if (outputDen_) { + MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &DenHist_[0], + rnemdLogWidth_, MPI::REALTYPE, MPI::SUM); + } + if (outputAh_) { + MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &AhCount_, + 1, MPI::REALTYPE, MPI::SUM); + } + if (outputVz_) { + MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &pzzHist_[0], + rnemdLogWidth_, MPI::REALTYPE, MPI::SUM); + } + // If we're the root node, should we print out the results int worldRank = MPI::COMM_WORLD.Get_rank(); if (worldRank == 0) { @@ -1550,7 +1627,24 @@ namespace OpenMD { } rotTempLog_ << endl; } - + // James put this in. + Mat3x3d hmat = currentSnap_->getHmat(); + if (outputDen_) { + denLog_ << time; + for (j = 0; j < rnemdLogWidth_; j++) { + + RealType binVol = hmat(0,0) * hmat(1,1) * (hmat(2,2) / float(nBins_)); + denLog_ << "\t" << DenHist_[j] / (float(logFrameCount_) * binVol); + } + denLog_ << endl; + } + if (outputVz_) { + vzzLog_ << time; + for (j = 0; j < rnemdLogWidth_; j++) { + vzzLog_ << "\t" << pzzHist_[j] / mHist_[j]; + } + vzzLog_ << endl; + } #ifdef IS_MPI } #endif @@ -1586,6 +1680,26 @@ namespace OpenMD { rotTempCount_[j] = 0; rotTempHist_[j] = 0.0; } + // James put this in + if (outputDen_) + for (j = 0; j < rnemdLogWidth_; j++) { + //pyzCount_[j] = 0; + DenHist_[j] = 0.0; + } + if (outputVz_) + for (j = 0; j < rnemdLogWidth_; j++) { + //pyzCount_[j] = 0; + pzzHist_[j] = 0.0; + } + // reset the counter + + Numcount_++; + if (Numcount_ > int(runTime_/statusTime_)) + cerr << "time =" << time << " Asum =" << Asum_ << '\n'; + if (Numcount_ > int(runTime_/statusTime_)) + cerr << "time =" << time << " Jsum =" << Jsum_ << '\n'; + + logFrameCount_ = 0; } }