--- branches/development/src/integrators/RNEMD.cpp 2012/05/24 20:59:54 1723 +++ branches/development/src/rnemd/RNEMD.cpp 2012/08/07 18:26:40 1773 @@ -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" @@ -49,11 +49,7 @@ #include "primitives/StuntDouble.hpp" #include "utils/PhysicalConstants.hpp" #include "utils/Tuple.hpp" - -#ifndef IS_MPI -#include "math/SeqRandNumGen.hpp" -#else -#include "math/ParallelRandNumGen.hpp" +#ifdef IS_MPI #include #endif @@ -65,30 +61,212 @@ namespace OpenMD { RNEMD::RNEMD(SimInfo* info) : info_(info), evaluator_(info), seleMan_(info), usePeriodicBoundaryConditions_(info->getSimParams()->getUsePeriodicBoundaryConditions()) { + trialCount_ = 0; failTrialCount_ = 0; failRootCount_ = 0; int seedValue; Globals * simParams = info->getSimParams(); + RNEMDParameters* rnemdParams = simParams->getRNEMDParameters(); - stringToEnumMap_["KineticSwap"] = rnemdKineticSwap; - stringToEnumMap_["KineticScale"] = rnemdKineticScale; - stringToEnumMap_["KineticScaleVAM"] = rnemdKineticScaleVAM; - stringToEnumMap_["KineticScaleAM"] = rnemdKineticScaleAM; - stringToEnumMap_["PxScale"] = rnemdPxScale; - stringToEnumMap_["PyScale"] = rnemdPyScale; - stringToEnumMap_["PzScale"] = rnemdPzScale; - stringToEnumMap_["Px"] = rnemdPx; - stringToEnumMap_["Py"] = rnemdPy; - stringToEnumMap_["Pz"] = rnemdPz; - stringToEnumMap_["ShiftScaleV"] = rnemdShiftScaleV; - stringToEnumMap_["ShiftScaleVAM"] = rnemdShiftScaleVAM; - stringToEnumMap_["Unknown"] = rnemdUnknown; + stringToMethod_["Swap"] = rnemdSwap; + stringToMethod_["NIVS"] = rnemdNIVS; + stringToMethod_["VSS"] = rnemdVSS; - rnemdObjectSelection_ = simParams->getRNEMD_objectSelection(); + stringToFluxType_["KE"] = rnemdKE; + stringToFluxType_["Px"] = rnemdPx; + stringToFluxType_["Py"] = rnemdPy; + stringToFluxType_["Pz"] = rnemdPz; + stringToFluxType_["KE+Px"] = rnemdKePx; + stringToFluxType_["KE+Py"] = rnemdKePy; + stringToFluxType_["KE+Pvector"] = rnemdKePvector; + + runTime_ = simParams->getRunTime(); + statusTime_ = simParams->getStatusTime(); + + rnemdObjectSelection_ = rnemdParams->getObjectSelection(); evaluator_.loadScriptString(rnemdObjectSelection_); seleMan_.setSelectionSet(evaluator_.evaluate()); + + const string methStr = rnemdParams->getMethod(); + bool hasFluxType = rnemdParams->haveFluxType(); + + string fluxStr; + if (hasFluxType) { + fluxStr = rnemdParams->getFluxType(); + } else { + sprintf(painCave.errMsg, + "RNEMD: No fluxType was set in the md file. This parameter,\n" + "\twhich must be one of the following values:\n" + "\tKE, Px, Py, Pz, KE+Px, KE+Py, KE+Pvector, must be set to\n" + "\tuse RNEMD\n"); + painCave.isFatal = 1; + painCave.severity = OPENMD_ERROR; + simError(); + } + + bool hasKineticFlux = rnemdParams->haveKineticFlux(); + bool hasMomentumFlux = rnemdParams->haveMomentumFlux(); + bool hasMomentumFluxVector = rnemdParams->haveMomentumFluxVector(); + bool hasSlabWidth = rnemdParams->haveSlabWidth(); + bool hasSlabACenter = rnemdParams->haveSlabACenter(); + bool hasSlabBCenter = rnemdParams->haveSlabBCenter(); + bool hasOutputFileName = rnemdParams->haveOutputFileName(); + bool hasOutputFields = rnemdParams->haveOutputFields(); + + map::iterator i; + i = stringToMethod_.find(methStr); + if (i != stringToMethod_.end()) + rnemdMethod_ = i->second; + else { + sprintf(painCave.errMsg, + "RNEMD: The current method,\n" + "\t\t%s is not one of the recognized\n" + "\texchange methods: Swap, NIVS, or VSS\n", + methStr.c_str()); + painCave.isFatal = 1; + painCave.severity = OPENMD_ERROR; + simError(); + } + + map::iterator j; + j = stringToFluxType_.find(fluxStr); + if (j != stringToFluxType_.end()) + rnemdFluxType_ = j->second; + else { + sprintf(painCave.errMsg, + "RNEMD: The current fluxType,\n" + "\t\t%s\n" + "\tis not one of the recognized flux types.\n", + fluxStr.c_str()); + painCave.isFatal = 1; + painCave.severity = OPENMD_ERROR; + simError(); + } + + bool methodFluxMismatch = false; + bool hasCorrectFlux = false; + switch(rnemdMethod_) { + case rnemdSwap: + switch (rnemdFluxType_) { + case rnemdKE: + hasCorrectFlux = hasKineticFlux; + break; + case rnemdPx: + case rnemdPy: + case rnemdPz: + hasCorrectFlux = hasMomentumFlux; + break; + default : + methodFluxMismatch = true; + break; + } + break; + case rnemdNIVS: + switch (rnemdFluxType_) { + case rnemdKE: + case rnemdRotKE: + case rnemdFullKE: + hasCorrectFlux = hasKineticFlux; + break; + case rnemdPx: + case rnemdPy: + case rnemdPz: + hasCorrectFlux = hasMomentumFlux; + break; + case rnemdKePx: + case rnemdKePy: + hasCorrectFlux = hasMomentumFlux && hasKineticFlux; + break; + default: + methodFluxMismatch = true; + break; + } + break; + case rnemdVSS: + switch (rnemdFluxType_) { + case rnemdKE: + case rnemdRotKE: + case rnemdFullKE: + hasCorrectFlux = hasKineticFlux; + break; + case rnemdPx: + case rnemdPy: + case rnemdPz: + hasCorrectFlux = hasMomentumFlux; + break; + case rnemdPvector: + hasCorrectFlux = hasMomentumFluxVector; + case rnemdKePx: + case rnemdKePy: + hasCorrectFlux = hasMomentumFlux && hasKineticFlux; + break; + case rnemdKePvector: + hasCorrectFlux = hasMomentumFluxVector && hasKineticFlux; + break; + default: + methodFluxMismatch = true; + break; + } + default: + break; + } + + if (methodFluxMismatch) { + sprintf(painCave.errMsg, + "RNEMD: The current method,\n" + "\t\t%s\n" + "\tcannot be used with the current flux type, %s\n", + methStr.c_str(), fluxStr.c_str()); + painCave.isFatal = 1; + painCave.severity = OPENMD_ERROR; + simError(); + } + if (!hasCorrectFlux) { + sprintf(painCave.errMsg, + "RNEMD: The current method,\n" + "\t%s, and flux type %s\n" + "\tdid not have the correct flux value specified. Options\n" + "\tinclude: kineticFlux, momentumFlux, and momentumFluxVector\n", + methStr.c_str(), fluxStr.c_str()); + painCave.isFatal = 1; + painCave.severity = OPENMD_ERROR; + simError(); + } + if (hasKineticFlux) { + kineticFlux_ = rnemdParams->getKineticFlux(); + } else { + kineticFlux_ = 0.0; + } + if (hasMomentumFluxVector) { + momentumFluxVector_ = rnemdParams->getMomentumFluxVector(); + } else { + momentumFluxVector_ = V3Zero; + if (hasMomentumFlux) { + RealType momentumFlux = rnemdParams->getMomentumFlux(); + switch (rnemdFluxType_) { + case rnemdPx: + momentumFluxVector_.x() = momentumFlux; + break; + case rnemdPy: + momentumFluxVector_.y() = momentumFlux; + break; + case rnemdPz: + momentumFluxVector_.z() = momentumFlux; + break; + case rnemdKePx: + momentumFluxVector_.x() = momentumFlux; + break; + case rnemdKePy: + momentumFluxVector_.y() = momentumFlux; + break; + default: + break; + } + } + } + // do some sanity checking int selectionCount = seleMan_.getSelectionCount(); @@ -96,7 +274,7 @@ namespace OpenMD { if (selectionCount > nIntegrable) { sprintf(painCave.errMsg, - "RNEMD: The current RNEMD_objectSelection,\n" + "RNEMD: The current objectSelection,\n" "\t\t%s\n" "\thas resulted in %d selected objects. However,\n" "\tthe total number of integrable objects in the system\n" @@ -109,219 +287,146 @@ namespace OpenMD { painCave.severity = OPENMD_WARNING; simError(); } - - const string st = simParams->getRNEMD_exchangeType(); - map::iterator i; - i = stringToEnumMap_.find(st); - rnemdType_ = (i == stringToEnumMap_.end()) ? RNEMD::rnemdUnknown : i->second; - if (rnemdType_ == rnemdUnknown) { - sprintf(painCave.errMsg, - "RNEMD: The current RNEMD_exchangeType,\n" - "\t\t%s\n" - "\tis not one of the recognized exchange types.\n", - st.c_str()); - painCave.isFatal = 1; - painCave.severity = OPENMD_ERROR; - simError(); - } - - outputTemp_ = false; - if (simParams->haveRNEMD_outputTemperature()) { - outputTemp_ = simParams->getRNEMD_outputTemperature(); - } else if ((rnemdType_ == rnemdKineticSwap) || - (rnemdType_ == rnemdKineticScale) || - (rnemdType_ == rnemdKineticScaleVAM) || - (rnemdType_ == rnemdKineticScaleAM)) { - outputTemp_ = true; - } - outputVx_ = false; - if (simParams->haveRNEMD_outputVx()) { - outputVx_ = simParams->getRNEMD_outputVx(); - } else if ((rnemdType_ == rnemdPx) || (rnemdType_ == rnemdPxScale)) { - outputVx_ = true; - } - outputVy_ = false; - if (simParams->haveRNEMD_outputVy()) { - outputVy_ = simParams->getRNEMD_outputVy(); - } else if ((rnemdType_ == rnemdPy) || (rnemdType_ == rnemdPyScale)) { - outputVy_ = true; - } - output3DTemp_ = false; - if (simParams->haveRNEMD_outputXyzTemperature()) { - output3DTemp_ = simParams->getRNEMD_outputXyzTemperature(); - } - outputRotTemp_ = false; - if (simParams->haveRNEMD_outputRotTemperature()) { - outputRotTemp_ = simParams->getRNEMD_outputRotTemperature(); - } + nBins_ = rnemdParams->getOutputBins(); -#ifdef IS_MPI - if (worldRank == 0) { -#endif + data_.resize(RNEMD::ENDINDEX); + OutputData z; + z.units = "Angstroms"; + z.title = "Z"; + z.dataType = "RealType"; + z.accumulator.reserve(nBins_); + for (unsigned int i = 0; i < nBins_; i++) + z.accumulator.push_back( new Accumulator() ); + data_[Z] = z; + outputMap_["Z"] = Z; - //may have rnemdWriter separately - string rnemdFileName; + OutputData temperature; + temperature.units = "K"; + temperature.title = "Temperature"; + temperature.dataType = "RealType"; + temperature.accumulator.reserve(nBins_); + for (unsigned int i = 0; i < nBins_; i++) + temperature.accumulator.push_back( new Accumulator() ); + data_[TEMPERATURE] = temperature; + outputMap_["TEMPERATURE"] = TEMPERATURE; - if (outputTemp_) { - rnemdFileName = "temperature.log"; - tempLog_.open(rnemdFileName.c_str()); - } - if (outputVx_) { - rnemdFileName = "velocityX.log"; - vxzLog_.open(rnemdFileName.c_str()); - } - if (outputVy_) { - rnemdFileName = "velocityY.log"; - vyzLog_.open(rnemdFileName.c_str()); - } + OutputData velocity; + velocity.units = "amu/fs"; + velocity.title = "Velocity"; + velocity.dataType = "Vector3d"; + velocity.accumulator.reserve(nBins_); + for (unsigned int i = 0; i < nBins_; i++) + velocity.accumulator.push_back( new VectorAccumulator() ); + data_[VELOCITY] = velocity; + outputMap_["VELOCITY"] = VELOCITY; - if (output3DTemp_) { - rnemdFileName = "temperatureX.log"; - xTempLog_.open(rnemdFileName.c_str()); - rnemdFileName = "temperatureY.log"; - yTempLog_.open(rnemdFileName.c_str()); - rnemdFileName = "temperatureZ.log"; - zTempLog_.open(rnemdFileName.c_str()); - } - if (outputRotTemp_) { - rnemdFileName = "temperatureR.log"; - rotTempLog_.open(rnemdFileName.c_str()); - } + OutputData density; + density.units = "g cm^-3"; + density.title = "Density"; + density.dataType = "RealType"; + density.accumulator.reserve(nBins_); + for (unsigned int i = 0; i < nBins_; i++) + density.accumulator.push_back( new Accumulator() ); + data_[DENSITY] = density; + outputMap_["DENSITY"] = DENSITY; -#ifdef IS_MPI - } -#endif - - set_RNEMD_exchange_time(simParams->getRNEMD_exchangeTime()); - set_RNEMD_nBins(simParams->getRNEMD_nBins()); - midBin_ = nBins_ / 2; - if (simParams->haveRNEMD_binShift()) { - if (simParams->getRNEMD_binShift()) { - zShift_ = 0.5 / (RealType)(nBins_); - } else { - zShift_ = 0.0; - } + if (hasOutputFields) { + parseOutputFileFormat(rnemdParams->getOutputFields()); } else { - zShift_ = 0.0; + outputMask_.set(Z); + switch (rnemdFluxType_) { + case rnemdKE: + case rnemdRotKE: + case rnemdFullKE: + outputMask_.set(TEMPERATURE); + break; + case rnemdPx: + case rnemdPy: + outputMask_.set(VELOCITY); + break; + case rnemdPz: + case rnemdPvector: + outputMask_.set(VELOCITY); + outputMask_.set(DENSITY); + break; + case rnemdKePx: + case rnemdKePy: + outputMask_.set(TEMPERATURE); + outputMask_.set(VELOCITY); + break; + case rnemdKePvector: + outputMask_.set(TEMPERATURE); + outputMask_.set(VELOCITY); + outputMask_.set(DENSITY); + break; + default: + break; + } } - //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()); - /*arbitary rnemdLogWidth_, no checking; - if (rnemdLogWidth_ != nBins_ && rnemdLogWidth_ != midBin_ + 1) { - cerr << "WARNING! RNEMD_logWidth has abnormal value!\n"; - cerr << "Automaically set back to default.\n"; - rnemdLogWidth_ = nBins_; - }*/ + + if (hasOutputFileName) { + rnemdFileName_ = rnemdParams->getOutputFileName(); } else { - set_RNEMD_logWidth(nBins_); - } - tempHist_.resize(rnemdLogWidth_, 0.0); - tempCount_.resize(rnemdLogWidth_, 0); - pxzHist_.resize(rnemdLogWidth_, 0.0); - //vxzCount_.resize(rnemdLogWidth_, 0); - pyzHist_.resize(rnemdLogWidth_, 0.0); - //vyzCount_.resize(rnemdLogWidth_, 0); + rnemdFileName_ = getPrefix(info->getFinalConfigFileName()) + ".rnemd"; + } - mHist_.resize(rnemdLogWidth_, 0.0); - xTempHist_.resize(rnemdLogWidth_, 0.0); - yTempHist_.resize(rnemdLogWidth_, 0.0); - zTempHist_.resize(rnemdLogWidth_, 0.0); - xyzTempCount_.resize(rnemdLogWidth_, 0); - rotTempHist_.resize(rnemdLogWidth_, 0.0); - rotTempCount_.resize(rnemdLogWidth_, 0); + exchangeTime_ = rnemdParams->getExchangeTime(); - set_RNEMD_exchange_total(0.0); - if (simParams->haveRNEMD_targetFlux()) { - set_RNEMD_target_flux(simParams->getRNEMD_targetFlux()); - } else { - set_RNEMD_target_flux(0.0); - } - if (simParams->haveRNEMD_targetJzKE()) { - set_RNEMD_target_JzKE(simParams->getRNEMD_targetJzKE()); - } else { - set_RNEMD_target_JzKE(0.0); - } - if (simParams->haveRNEMD_targetJzpx()) { - set_RNEMD_target_jzpx(simParams->getRNEMD_targetJzpx()); - } else { - set_RNEMD_target_jzpx(0.0); - } - jzp_.x() = targetJzpx_; - njzp_.x() = -targetJzpx_; - if (simParams->haveRNEMD_targetJzpy()) { - set_RNEMD_target_jzpy(simParams->getRNEMD_targetJzpy()); - } else { - set_RNEMD_target_jzpy(0.0); - } - jzp_.y() = targetJzpy_; - njzp_.y() = -targetJzpy_; - if (simParams->haveRNEMD_targetJzpz()) { - set_RNEMD_target_jzpz(simParams->getRNEMD_targetJzpz()); - } else { - set_RNEMD_target_jzpz(0.0); - } - jzp_.z() = targetJzpz_; - njzp_.z() = -targetJzpz_; + Snapshot* currentSnap_ = info->getSnapshotManager()->getCurrentSnapshot(); + Mat3x3d hmat = currentSnap_->getHmat(); + + // Target exchange quantities (in each exchange) = 2 Lx Ly dt flux + // Lx, Ly = box dimensions in x & y + // dt = exchange time interval + // flux = target flux -#ifndef IS_MPI - if (simParams->haveSeed()) { - seedValue = simParams->getSeed(); - randNumGen_ = new SeqRandNumGen(seedValue); - }else { - randNumGen_ = new SeqRandNumGen(); - } -#else - if (simParams->haveSeed()) { - seedValue = simParams->getSeed(); - randNumGen_ = new ParallelRandNumGen(seedValue); - }else { - randNumGen_ = new ParallelRandNumGen(); - } -#endif - } + kineticTarget_ = 2.0*kineticFlux_*exchangeTime_*hmat(0,0)*hmat(1,1); + momentumTarget_ = 2.0*momentumFluxVector_*exchangeTime_*hmat(0,0)*hmat(1,1); + + // total exchange sums are zeroed out at the beginning: + + kineticExchange_ = 0.0; + momentumExchange_ = V3Zero; + + if (hasSlabWidth) + slabWidth_ = rnemdParams->getSlabWidth(); + else + slabWidth_ = hmat(2,2) / 10.0; + if (hasSlabACenter) + slabACenter_ = rnemdParams->getSlabACenter(); + else + slabACenter_ = 0.0; + + if (hasSlabBCenter) + slabBCenter_ = rnemdParams->getSlabBCenter(); + else + slabBCenter_ = hmat(2,2) / 2.0; + + } + RNEMD::~RNEMD() { - delete randNumGen_; #ifdef IS_MPI if (worldRank == 0) { #endif - - sprintf(painCave.errMsg, - "RNEMD: total failed trials: %d\n", - failTrialCount_); - painCave.isFatal = 0; - painCave.severity = OPENMD_INFO; - simError(); - if (outputTemp_) tempLog_.close(); - if (outputVx_) vxzLog_.close(); - if (outputVy_) vyzLog_.close(); + writeOutputFile(); - if (rnemdType_ == rnemdKineticScale || rnemdType_ == rnemdPxScale || - rnemdType_ == rnemdPyScale) { - sprintf(painCave.errMsg, - "RNEMD: total root-checking warnings: %d\n", - failRootCount_); - painCave.isFatal = 0; - painCave.severity = OPENMD_INFO; - simError(); - } - if (output3DTemp_) { - xTempLog_.close(); - yTempLog_.close(); - zTempLog_.close(); - } - if (outputRotTemp_) rotTempLog_.close(); - + rnemdFile_.close(); + #ifdef IS_MPI } #endif } + + bool RNEMD::inSlabA(Vector3d pos) { + return (abs(pos.z() - slabACenter_) < 0.5*slabWidth_); + } + bool RNEMD::inSlabB(Vector3d pos) { + return (abs(pos.z() - slabBCenter_) < 0.5*slabWidth_); + } void RNEMD::doSwap() { @@ -353,22 +458,17 @@ namespace OpenMD { if (usePeriodicBoundaryConditions_) currentSnap_->wrapVector(pos); + bool inA = inSlabA(pos); + bool inB = inSlabB(pos); - // which bin is this stuntdouble in? - // wrapped positions are in the range [-0.5*hmat(2,2), +0.5*hmat(2,2)] - - int binNo = int(nBins_ * (pos.z() / hmat(2,2) + zShift_ + 0.5)) % nBins_; - - - // if we're in bin 0 or the middleBin - if (binNo == 0 || binNo == midBin_) { + if (inA || inB) { RealType mass = sd->getMass(); Vector3d vel = sd->getVel(); RealType value; - - switch(rnemdType_) { - case rnemdKineticSwap : + + switch(rnemdFluxType_) { + case rnemdKE : value = mass * vel.lengthSquare(); @@ -389,7 +489,7 @@ namespace OpenMD { } } //angular momenta exchange enabled //energyConvert temporarily disabled - //make exchangeSum_ comparable between swap & scale + //make kineticExchange_ comparable between swap & scale //value = value * 0.5 / PhysicalConstants::energyConvert; value *= 0.5; break; @@ -406,7 +506,7 @@ namespace OpenMD { break; } - if (binNo == 0) { + if (inA == 0) { if (!min_found) { min_val = value; min_sd = sd; @@ -417,7 +517,7 @@ namespace OpenMD { min_sd = sd; } } - } else { //midBin_ + } else { if (!max_found) { max_val = value; max_sd = sd; @@ -431,10 +531,10 @@ namespace OpenMD { } } } - + #ifdef IS_MPI int nProc, worldRank; - + nProc = MPI::COMM_WORLD.Get_size(); worldRank = MPI::COMM_WORLD.Get_rank(); @@ -454,7 +554,7 @@ namespace OpenMD { RealType val; int rank; } max_vals, min_vals; - + if (my_min_found) { min_vals.val = min_val; } else { @@ -492,8 +592,8 @@ namespace OpenMD { Vector3d max_vel = max_sd->getVel(); RealType temp_vel; - switch(rnemdType_) { - case rnemdKineticSwap : + switch(rnemdFluxType_) { + case rnemdKE : min_sd->setVel(max_vel); max_sd->setVel(min_vel); if (min_sd->isDirectional() && max_sd->isDirectional()) { @@ -544,8 +644,8 @@ namespace OpenMD { min_vel.getArrayPointer(), 3, MPI::REALTYPE, min_vals.rank, 0, status); - switch(rnemdType_) { - case rnemdKineticSwap : + switch(rnemdFluxType_) { + case rnemdKE : max_sd->setVel(min_vel); //angular momenta exchange enabled if (max_sd->isDirectional()) { @@ -590,8 +690,8 @@ namespace OpenMD { max_vel.getArrayPointer(), 3, MPI::REALTYPE, max_vals.rank, 0, status); - switch(rnemdType_) { - case rnemdKineticSwap : + switch(rnemdFluxType_) { + case rnemdKE : min_sd->setVel(max_vel); //angular momenta exchange enabled if (min_sd->isDirectional()) { @@ -625,10 +725,28 @@ namespace OpenMD { } } #endif - exchangeSum_ += max_val - min_val; + + switch(rnemdFluxType_) { + case rnemdKE: + cerr << "KE\n"; + kineticExchange_ += max_val - min_val; + break; + case rnemdPx: + momentumExchange_.x() += max_val - min_val; + break; + case rnemdPy: + momentumExchange_.y() += max_val - min_val; + break; + case rnemdPz: + momentumExchange_.z() += max_val - min_val; + break; + default: + cerr << "default\n"; + break; + } } else { sprintf(painCave.errMsg, - "RNEMD: exchange NOT performed because min_val > max_val\n"); + "RNEMD::doSwap exchange NOT performed because min_val > max_val\n"); painCave.isFatal = 0; painCave.severity = OPENMD_INFO; simError(); @@ -636,17 +754,16 @@ namespace OpenMD { } } else { sprintf(painCave.errMsg, - "RNEMD: exchange NOT performed because selected object\n" - "\tnot present in at least one of the two slabs.\n"); + "RNEMD::doSwap exchange NOT performed because selected object\n" + "\twas not present in at least one of the two slabs.\n"); painCave.isFatal = 0; painCave.severity = OPENMD_INFO; simError(); failTrialCount_++; - } - + } } - void RNEMD::doScale() { + void RNEMD::doNIVS() { Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot(); Mat3x3d hmat = currentSnap_->getHmat(); @@ -687,17 +804,15 @@ namespace OpenMD { currentSnap_->wrapVector(pos); // which bin is this stuntdouble in? - // wrapped positions are in the range [-0.5*hmat(2,2), +0.5*hmat(2,2)] + bool inA = inSlabA(pos); + bool inB = inSlabB(pos); - int binNo = int(nBins_ * (pos.z() / hmat(2,2) + zShift_ + 0.5)) % nBins_; - - // if we're in bin 0 or the middleBin - if (binNo == 0 || binNo == midBin_) { - + if (inA || inB) { + RealType mass = sd->getMass(); Vector3d vel = sd->getVel(); - if (binNo == 0) { + if (inA) { hotBin.push_back(sd); Phx += mass * vel.x(); Phy += mass * vel.y(); @@ -705,7 +820,6 @@ namespace OpenMD { Khx += mass * vel.x() * vel.x(); Khy += mass * vel.y() * vel.y(); Khz += mass * vel.z() * vel.z(); - //if (rnemdType_ == rnemdKineticScaleVAM) { if (sd->isDirectional()) { Vector3d angMom = sd->getJ(); Mat3x3d I = sd->getI(); @@ -721,8 +835,7 @@ namespace OpenMD { + angMom[2]*angMom[2]/I(2, 2); } } - //} - } else { //midBin_ + } else { coldBin.push_back(sd); Pcx += mass * vel.x(); Pcy += mass * vel.y(); @@ -730,7 +843,6 @@ namespace OpenMD { Kcx += mass * vel.x() * vel.x(); Kcy += mass * vel.y() * vel.y(); Kcz += mass * vel.z() * vel.z(); - //if (rnemdType_ == rnemdKineticScaleVAM) { if (sd->isDirectional()) { Vector3d angMom = sd->getJ(); Mat3x3d I = sd->getI(); @@ -746,7 +858,6 @@ namespace OpenMD { + angMom[2]*angMom[2]/I(2, 2); } } - //} } } } @@ -760,12 +871,6 @@ 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= " < 0 - if (rnemdType_ == rnemdKineticScaleVAM) { - c = 1.0 - targetFlux_ / (Kcx + Kcy + Kcz + Kcw); + if (rnemdFluxType_ == rnemdFullKE) { + c = 1.0 - kineticTarget_ / (Kcx + Kcy + Kcz + Kcw); } else { - c = 1.0 - targetFlux_ / Kcw; + c = 1.0 - kineticTarget_ / Kcw; } if ((c > 0.81) && (c < 1.21)) {//restrict scaling coefficients c = sqrt(c); - std::cerr << "cold slab scaling coefficient: " << c << endl; + //std::cerr << "cold slab scaling coefficient: " << c << endl; //now convert to hotBin coefficient RealType w = 0.0; - if (rnemdType_ == rnemdKineticScaleVAM) { + if (rnemdFluxType_ == rnemdFullKE) { x = 1.0 + px * (1.0 - c); y = 1.0 + py * (1.0 - c); z = 1.0 + pz * (1.0 - c); @@ -820,37 +925,32 @@ namespace OpenMD { */ if ((fabs(x - 1.0) < 0.1) && (fabs(y - 1.0) < 0.1) && (fabs(z - 1.0) < 0.1)) { - w = 1.0 + (targetFlux_ + Khx * (1.0 - x * x) + Khy * (1.0 - y * y) + w = 1.0 + (kineticTarget_ + + Khx * (1.0 - x * x) + Khy * (1.0 - y * y) + Khz * (1.0 - z * z)) / Khw; }//no need to calculate w if x, y or z is out of range } else { - w = 1.0 + targetFlux_ / Khw; + w = 1.0 + kineticTarget_ / Khw; } if ((w > 0.81) && (w < 1.21)) {//restrict scaling coefficients //if w is in the right range, so should be x, y, z. vector::iterator sdi; Vector3d vel; for (sdi = coldBin.begin(); sdi != coldBin.end(); sdi++) { - if (rnemdType_ == rnemdKineticScaleVAM) { + if (rnemdFluxType_ == rnemdFullKE) { vel = (*sdi)->getVel() * c; - //vel.x() *= c; - //vel.y() *= c; - //vel.z() *= c; (*sdi)->setVel(vel); } if ((*sdi)->isDirectional()) { Vector3d angMom = (*sdi)->getJ() * c; - //angMom[0] *= c; - //angMom[1] *= c; - //angMom[2] *= c; (*sdi)->setJ(angMom); } } w = sqrt(w); - std::cerr << "xh= " << x << "\tyh= " << y << "\tzh= " << z - << "\twh= " << w << endl; + // std::cerr << "xh= " << x << "\tyh= " << y << "\tzh= " << z + // << "\twh= " << w << endl; for (sdi = hotBin.begin(); sdi != hotBin.end(); sdi++) { - if (rnemdType_ == rnemdKineticScaleVAM) { + if (rnemdFluxType_ == rnemdFullKE) { vel = (*sdi)->getVel(); vel.x() *= x; vel.y() *= y; @@ -859,45 +959,42 @@ namespace OpenMD { } if ((*sdi)->isDirectional()) { Vector3d angMom = (*sdi)->getJ() * w; - //angMom[0] *= w; - //angMom[1] *= w; - //angMom[2] *= w; (*sdi)->setJ(angMom); } } successfulScale = true; - exchangeSum_ += targetFlux_; + kineticExchange_ += kineticTarget_; } } } else { RealType a000, a110, c0, a001, a111, b01, b11, c1; - switch(rnemdType_) { - case rnemdKineticScale : + switch(rnemdFluxType_) { + case rnemdKE : /* used hotBin coeff's & only scale x & y dimensions RealType px = Phx / Pcx; RealType py = Phy / Pcy; a110 = Khy; - c0 = - Khx - Khy - targetFlux_; + c0 = - Khx - Khy - kineticTarget_; a000 = Khx; a111 = Kcy * py * py; b11 = -2.0 * Kcy * py * (1.0 + py); - c1 = Kcy * py * (2.0 + py) + Kcx * px * ( 2.0 + px) + targetFlux_; + c1 = Kcy * py * (2.0 + py) + Kcx * px * ( 2.0 + px) + kineticTarget_; b01 = -2.0 * Kcx * px * (1.0 + px); a001 = Kcx * px * px; */ //scale all three dimensions, let c_x = c_y a000 = Kcx + Kcy; a110 = Kcz; - c0 = targetFlux_ - Kcx - Kcy - Kcz; + c0 = kineticTarget_ - Kcx - Kcy - Kcz; a001 = Khx * px * px + Khy * py * py; a111 = Khz * pz * pz; b01 = -2.0 * (Khx * px * (1.0 + px) + Khy * py * (1.0 + py)); b11 = -2.0 * Khz * pz * (1.0 + pz); c1 = Khx * px * (2.0 + px) + Khy * py * (2.0 + py) - + Khz * pz * (2.0 + pz) - targetFlux_; + + Khz * pz * (2.0 + pz) - kineticTarget_; break; - case rnemdPxScale : - c = 1 - targetFlux_ / Pcx; + case rnemdPx : + c = 1 - momentumTarget_.x() / Pcx; a000 = Kcy; a110 = Kcz; c0 = Kcx * c * c - Kcx - Kcy - Kcz; @@ -908,8 +1005,8 @@ namespace OpenMD { c1 = Khy * py * (2.0 + py) + Khz * pz * (2.0 + pz) + Khx * (fastpow(c * px - px - 1.0, 2) - 1.0); break; - case rnemdPyScale : - c = 1 - targetFlux_ / Pcy; + case rnemdPy : + c = 1 - momentumTarget_.y() / Pcy; a000 = Kcx; a110 = Kcz; c0 = Kcy * c * c - Kcx - Kcy - Kcz; @@ -920,8 +1017,8 @@ namespace OpenMD { c1 = Khx * px * (2.0 + px) + Khz * pz * (2.0 + pz) + Khy * (fastpow(c * py - py - 1.0, 2) - 1.0); break; - case rnemdPzScale ://we don't really do this, do we? - c = 1 - targetFlux_ / Pcz; + case rnemdPz ://we don't really do this, do we? + c = 1 - momentumTarget_.z() / Pcz; a000 = Kcx; a110 = Kcy; c0 = Kcz * c * c - Kcx - Kcy - Kcz; @@ -1006,21 +1103,21 @@ namespace OpenMD { for (rpi = rps.begin(); rpi != rps.end(); rpi++) { r1 = (*rpi).first; r2 = (*rpi).second; - switch(rnemdType_) { - case rnemdKineticScale : + switch(rnemdFluxType_) { + case rnemdKE : diff = fastpow(1.0 - r1, 2) + fastpow(1.0 - r2, 2) + fastpow(r1 * r1 / r2 / r2 - Kcz/Kcx, 2) + fastpow(r1 * r1 / r2 / r2 - Kcz/Kcy, 2); break; - case rnemdPxScale : + case rnemdPx : diff = fastpow(1.0 - r1, 2) + fastpow(1.0 - r2, 2) + fastpow(r1 * r1 / r2 / r2 - Kcz/Kcy, 2); break; - case rnemdPyScale : + case rnemdPy : diff = fastpow(1.0 - r1, 2) + fastpow(1.0 - r2, 2) + fastpow(r1 * r1 / r2 / r2 - Kcz/Kcx, 2); break; - case rnemdPzScale : + case rnemdPz : diff = fastpow(1.0 - r1, 2) + fastpow(1.0 - r2, 2) + fastpow(r1 * r1 / r2 / r2 - Kcy/Kcx, 2); default : @@ -1034,33 +1131,33 @@ namespace OpenMD { #ifdef IS_MPI if (worldRank == 0) { #endif - sprintf(painCave.errMsg, - "RNEMD: roots r1= %lf\tr2 = %lf\n", - bestPair.first, bestPair.second); - painCave.isFatal = 0; - painCave.severity = OPENMD_INFO; - simError(); + // sprintf(painCave.errMsg, + // "RNEMD: roots r1= %lf\tr2 = %lf\n", + // bestPair.first, bestPair.second); + // painCave.isFatal = 0; + // painCave.severity = OPENMD_INFO; + // simError(); #ifdef IS_MPI } #endif - switch(rnemdType_) { - case rnemdKineticScale : + switch(rnemdFluxType_) { + case rnemdKE : x = bestPair.first; y = bestPair.first; z = bestPair.second; break; - case rnemdPxScale : + case rnemdPx : x = c; y = bestPair.first; z = bestPair.second; break; - case rnemdPyScale : + case rnemdPy : x = bestPair.first; y = c; z = bestPair.second; break; - case rnemdPzScale : + case rnemdPz : x = bestPair.first; y = bestPair.second; z = c; @@ -1089,12 +1186,25 @@ namespace OpenMD { (*sdi)->setVel(vel); } successfulScale = true; - exchangeSum_ += targetFlux_; + switch(rnemdFluxType_) { + case rnemdKE : + kineticExchange_ += kineticTarget_; + break; + case rnemdPx : + case rnemdPy : + case rnemdPz : + momentumExchange_ += momentumTarget_; + break; + default : + break; + } } } if (successfulScale != true) { sprintf(painCave.errMsg, - "RNEMD: exchange NOT performed!\n"); + "RNEMD::doNIVS exchange NOT performed - roots that solve\n" + "\tthe constraint equations may not exist or there may be\n" + "\tno selected objects in one or both slabs.\n"); painCave.isFatal = 0; painCave.severity = OPENMD_INFO; simError(); @@ -1102,9 +1212,10 @@ namespace OpenMD { } } - void RNEMD::doShiftScale() { + void RNEMD::doVSS() { Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot(); + RealType time = currentSnap_->getTime(); Mat3x3d hmat = currentSnap_->getHmat(); seleMan_.setSelectionSet(evaluator_.evaluate()); @@ -1121,6 +1232,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)) { @@ -1135,24 +1247,22 @@ namespace OpenMD { currentSnap_->wrapVector(pos); // which bin is this stuntdouble in? - // wrapped positions are in the range [-0.5*hmat(2,2), +0.5*hmat(2,2)] - - int binNo = int(nBins_ * (pos.z() / hmat(2,2) + zShift_ + 0.5)) % nBins_; - - // if we're in bin 0 or the middleBin - if (binNo == 0 || binNo == midBin_) { + bool inA = inSlabA(pos); + bool inB = inSlabB(pos); + + if (inA || inB) { RealType mass = sd->getMass(); Vector3d vel = sd->getVel(); - if (binNo == 0) { + if (inA) { hotBin.push_back(sd); //std::cerr << "before, velocity = " << vel << endl; Ph += mass * vel; //std::cerr << "after, velocity = " << vel << endl; Mh += mass; Kh += mass * vel.lengthSquare(); - if (rnemdType_ == rnemdShiftScaleVAM) { + if (rnemdFluxType_ == rnemdFullKE) { if (sd->isDirectional()) { Vector3d angMom = sd->getJ(); Mat3x3d I = sd->getI(); @@ -1174,7 +1284,7 @@ namespace OpenMD { Pc += mass * vel; Mc += mass; Kc += mass * vel.lengthSquare(); - if (rnemdType_ == rnemdShiftScaleVAM) { + if (rnemdFluxType_ == rnemdFullKE) { if (sd->isDirectional()) { Vector3d angMom = sd->getJ(); Mat3x3d I = sd->getI(); @@ -1198,10 +1308,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); @@ -1214,31 +1324,33 @@ namespace OpenMD { bool successfulExchange = false; if ((Mh > 0.0) && (Mc > 0.0)) {//both slabs are not empty Vector3d vc = Pc / Mc; - Vector3d ac = njzp_ / Mc + vc; - RealType cNumerator = Kc - targetJzKE_ - 0.5 * Mc * ac.lengthSquare(); + Vector3d ac = -momentumTarget_ / Mc + vc; + Vector3d acrec = -momentumTarget_ / Mc; + RealType cNumerator = Kc - kineticTarget_ - 0.5 * Mc * ac.lengthSquare(); if (cNumerator > 0.0) { RealType cDenominator = Kc - 0.5 * Mc * vc.lengthSquare(); if (cDenominator > 0.0) { RealType c = sqrt(cNumerator / cDenominator); if ((c > 0.9) && (c < 1.1)) {//restrict scaling coefficients Vector3d vh = Ph / Mh; - Vector3d ah = jzp_ / Mh + vh; - RealType hNumerator = Kh + targetJzKE_ + Vector3d ah = momentumTarget_ / Mh + vh; + Vector3d ahrec = momentumTarget_ / Mh; + RealType hNumerator = Kh + kineticTarget_ - 0.5 * Mh * ah.lengthSquare(); if (hNumerator > 0.0) { RealType hDenominator = Kh - 0.5 * Mh * vh.lengthSquare(); 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++) { //vel = (*sdi)->getVel(); vel = ((*sdi)->getVel() - vc) * c + ac; (*sdi)->setVel(vel); - if (rnemdType_ == rnemdShiftScaleVAM) { + if (rnemdFluxType_ == rnemdFullKE) { if ((*sdi)->isDirectional()) { Vector3d angMom = (*sdi)->getJ() * c; (*sdi)->setJ(angMom); @@ -1249,7 +1361,7 @@ namespace OpenMD { //vel = (*sdi)->getVel(); vel = ((*sdi)->getVel() - vh) * h + ah; (*sdi)->setVel(vel); - if (rnemdType_ == rnemdShiftScaleVAM) { + if (rnemdFluxType_ == rnemdFullKE) { if ((*sdi)->isDirectional()) { Vector3d angMom = (*sdi)->getJ() * h; (*sdi)->setJ(angMom); @@ -1257,10 +1369,8 @@ namespace OpenMD { } } successfulExchange = true; - exchangeSum_ += targetFlux_; - // 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. + kineticExchange_ += kineticTarget_; + momentumExchange_ += momentumTarget_; } } } @@ -1270,7 +1380,9 @@ namespace OpenMD { } if (successfulExchange != true) { sprintf(painCave.errMsg, - "RNEMD: exchange NOT performed!\n"); + "RNEMD::doVSS exchange NOT performed - roots that solve\n" + "\tthe constraint equations may not exist or there may be\n" + "\tno selected objects in one or both slabs.\n"); painCave.isFatal = 0; painCave.severity = OPENMD_INFO; simError(); @@ -1280,26 +1392,18 @@ namespace OpenMD { void RNEMD::doRNEMD() { - switch(rnemdType_) { - case rnemdKineticScale : - case rnemdKineticScaleVAM : - case rnemdKineticScaleAM : - case rnemdPxScale : - case rnemdPyScale : - case rnemdPzScale : - doScale(); - break; - case rnemdKineticSwap : - case rnemdPx : - case rnemdPy : - case rnemdPz : + trialCount_++; + switch(rnemdMethod_) { + case rnemdSwap: doSwap(); break; - case rnemdShiftScaleV : - case rnemdShiftScaleVAM : - doShiftScale(); + case rnemdNIVS: + doNIVS(); break; - case rnemdUnknown : + case rnemdVSS: + doVSS(); + break; + case rnemdUnkownMethod: default : break; } @@ -1316,19 +1420,27 @@ namespace OpenMD { StuntDouble* sd; int idx; + vector binMass(nBins_, 0.0); + vector binPx(nBins_, 0.0); + vector binPy(nBins_, 0.0); + vector binPz(nBins_, 0.0); + vector binKE(nBins_, 0.0); + vector binDOF(nBins_, 0); + vector binCount(nBins_, 0); + // 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)) { @@ -1341,251 +1453,268 @@ namespace OpenMD { if (usePeriodicBoundaryConditions_) currentSnap_->wrapVector(pos); - + // which bin is this stuntdouble in? // wrapped positions are in the range [-0.5*hmat(2,2), +0.5*hmat(2,2)] - - int binNo = int(rnemdLogWidth_ * (pos.z() / hmat(2,2) + 0.5)) % - rnemdLogWidth_; - // no symmetrization allowed due to arbitary rnemdLogWidth_ - /* - if (rnemdLogWidth_ == midBin_ + 1) - if (binNo > midBin_) - binNo = nBins_ - binNo; - */ + // Shift molecules by half a box to have bins start at 0 + // The modulo operator is used to wrap the case when we are + // beyond the end of the bins back to the beginning. + int binNo = int(nBins_ * (pos.z() / hmat(2,2) + 0.5)) % nBins_; + RealType mass = sd->getMass(); - mHist_[binNo] += mass; Vector3d vel = sd->getVel(); - RealType value; - //RealType xVal, yVal, zVal; - if (outputTemp_) { - value = mass * vel.lengthSquare(); - tempCount_[binNo] += 3; - if (sd->isDirectional()) { - Vector3d angMom = sd->getJ(); - Mat3x3d I = sd->getI(); - if (sd->isLinear()) { - int i = sd->linearAxis(); - int j = (i + 1) % 3; - int k = (i + 2) % 3; - value += angMom[j] * angMom[j] / I(j, j) + - angMom[k] * angMom[k] / I(k, k); - tempCount_[binNo] +=2; - } else { - value += angMom[0] * angMom[0] / I(0, 0) + - angMom[1]*angMom[1]/I(1, 1) + - angMom[2]*angMom[2]/I(2, 2); - tempCount_[binNo] +=3; - } - } - value = value / PhysicalConstants::energyConvert - / PhysicalConstants::kb;//may move to getStatus() - tempHist_[binNo] += value; - } - if (outputVx_) { - value = mass * vel[0]; - //vxzCount_[binNo]++; - pxzHist_[binNo] += value; - } - if (outputVy_) { - value = mass * vel[1]; - //vyzCount_[binNo]++; - pyzHist_[binNo] += value; - } + binCount[binNo]++; + binMass[binNo] += mass; + binPx[binNo] += mass*vel.x(); + binPy[binNo] += mass*vel.y(); + binPz[binNo] += mass*vel.z(); + binKE[binNo] += 0.5 * (mass * vel.lengthSquare()); + binDOF[binNo] += 3; - if (output3DTemp_) { - value = mass * vel.x() * vel.x(); - xTempHist_[binNo] += value; - value = mass * vel.y() * vel.y() / PhysicalConstants::energyConvert - / PhysicalConstants::kb; - yTempHist_[binNo] += value; - value = mass * vel.z() * vel.z() / PhysicalConstants::energyConvert - / PhysicalConstants::kb; - zTempHist_[binNo] += value; - xyzTempCount_[binNo]++; + if (sd->isDirectional()) { + Vector3d angMom = sd->getJ(); + Mat3x3d I = sd->getI(); + if (sd->isLinear()) { + int i = sd->linearAxis(); + int j = (i + 1) % 3; + int k = (i + 2) % 3; + binKE[binNo] += 0.5 * (angMom[j] * angMom[j] / I(j, j) + + angMom[k] * angMom[k] / I(k, k)); + binDOF[binNo] += 2; + } else { + binKE[binNo] += 0.5 * (angMom[0] * angMom[0] / I(0, 0) + + angMom[1] * angMom[1] / I(1, 1) + + angMom[2] * angMom[2] / I(2, 2)); + binDOF[binNo] += 3; + } } - if (outputRotTemp_) { - if (sd->isDirectional()) { - Vector3d angMom = sd->getJ(); - Mat3x3d I = sd->getI(); - if (sd->isLinear()) { - int i = sd->linearAxis(); - int j = (i + 1) % 3; - int k = (i + 2) % 3; - value = angMom[j] * angMom[j] / I(j, j) + - angMom[k] * angMom[k] / I(k, k); - rotTempCount_[binNo] +=2; - } else { - value = angMom[0] * angMom[0] / I(0, 0) + - angMom[1] * angMom[1] / I(1, 1) + - angMom[2] * angMom[2] / I(2, 2); - rotTempCount_[binNo] +=3; - } - } - value = value / PhysicalConstants::energyConvert - / PhysicalConstants::kb;//may move to getStatus() - rotTempHist_[binNo] += value; - } + } + +#ifdef IS_MPI + MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binCount[0], + nBins_, MPI::INT, MPI::SUM); + MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binMass[0], + nBins_, MPI::REALTYPE, MPI::SUM); + MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binPx[0], + nBins_, MPI::REALTYPE, MPI::SUM); + MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binPy[0], + nBins_, MPI::REALTYPE, MPI::SUM); + MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binPz[0], + nBins_, MPI::REALTYPE, MPI::SUM); + MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binKE[0], + nBins_, MPI::REALTYPE, MPI::SUM); + MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binDOF[0], + nBins_, MPI::INT, MPI::SUM); +#endif + + Vector3d vel; + RealType den; + RealType temp; + RealType z; + for (int i = 0; i < nBins_; i++) { + z = (((RealType)i + 0.5) / (RealType)nBins_) * hmat(2,2); + vel.x() = binPx[i] / binMass[i]; + vel.y() = binPy[i] / binMass[i]; + vel.z() = binPz[i] / binMass[i]; + den = binCount[i] * nBins_ / (hmat(0,0) * hmat(1,1) * hmat(2,2)); + temp = 2.0 * binKE[i] / (binDOF[i] * PhysicalConstants::kb * + PhysicalConstants::energyConvert); + + for (unsigned int j = 0; j < outputMask_.size(); ++j) { + if(outputMask_[j]) { + switch(j) { + case Z: + (data_[j].accumulator[i])->add(z); + break; + case TEMPERATURE: + data_[j].accumulator[i]->add(temp); + break; + case VELOCITY: + dynamic_cast(data_[j].accumulator[i])->add(vel); + break; + case DENSITY: + data_[j].accumulator[i]->add(den); + break; + } + } + } } } void RNEMD::getStarted() { collectData(); - /*now can output profile in step 0, but might not be useful; - Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot(); - Stats& stat = currentSnap_->statData; - stat[Stats::RNEMD_EXCHANGE_TOTAL] = exchangeSum_; - */ - //may output a header for the log file here - getStatus(); + writeOutputFile(); } - 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; - + void RNEMD::parseOutputFileFormat(const std::string& format) { + StringTokenizer tokenizer(format, " ,;|\t\n\r"); + + while(tokenizer.hasMoreTokens()) { + std::string token(tokenizer.nextToken()); + toUpper(token); + OutputMapType::iterator i = outputMap_.find(token); + if (i != outputMap_.end()) { + outputMask_.set(i->second); + } else { + sprintf( painCave.errMsg, + "RNEMD::parseOutputFileFormat: %s is not a recognized\n" + "\toutputFileFormat keyword.\n", token.c_str() ); + painCave.isFatal = 0; + painCave.severity = OPENMD_ERROR; + simError(); + } + } + } + + void RNEMD::writeOutputFile() { + #ifdef IS_MPI - - // all processors have the same number of bins, and STL vectors pack their - // arrays, so in theory, this should be safe: - - MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &mHist_[0], - rnemdLogWidth_, MPI::REALTYPE, MPI::SUM); - if (outputTemp_) { - MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &tempHist_[0], - rnemdLogWidth_, MPI::REALTYPE, MPI::SUM); - MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &tempCount_[0], - rnemdLogWidth_, MPI::INT, MPI::SUM); - } - if (outputVx_) { - MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &pxzHist_[0], - rnemdLogWidth_, MPI::REALTYPE, MPI::SUM); - //MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &vxzCount_[0], - // rnemdLogWidth_, MPI::INT, MPI::SUM); - } - if (outputVy_) { - MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &pyzHist_[0], - rnemdLogWidth_, MPI::REALTYPE, MPI::SUM); - //MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &vyzCount_[0], - // rnemdLogWidth_, MPI::INT, MPI::SUM); - } - if (output3DTemp_) { - MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &xTempHist_[0], - rnemdLogWidth_, MPI::REALTYPE, MPI::SUM); - MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &yTempHist_[0], - rnemdLogWidth_, MPI::REALTYPE, MPI::SUM); - MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &zTempHist_[0], - rnemdLogWidth_, MPI::REALTYPE, MPI::SUM); - MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &xyzTempCount_[0], - rnemdLogWidth_, MPI::INT, MPI::SUM); - } - if (outputRotTemp_) { - MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &rotTempHist_[0], - rnemdLogWidth_, MPI::REALTYPE, MPI::SUM); - MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &rotTempCount_[0], - rnemdLogWidth_, MPI::INT, MPI::SUM); - } - // If we're the root node, should we print out the results int worldRank = MPI::COMM_WORLD.Get_rank(); if (worldRank == 0) { #endif + rnemdFile_.open(rnemdFileName_.c_str(), std::ios::out | std::ios::trunc ); + + if( !rnemdFile_ ){ + sprintf( painCave.errMsg, + "Could not open \"%s\" for RNEMD output.\n", + rnemdFileName_.c_str()); + painCave.isFatal = 1; + simError(); + } - if (outputTemp_) { - tempLog_ << time; - for (j = 0; j < rnemdLogWidth_; j++) { - tempLog_ << "\t" << tempHist_[j] / (RealType)tempCount_[j]; - } - tempLog_ << endl; + Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot(); + + RealType time = currentSnap_->getTime(); + + + rnemdFile_ << "#######################################################\n"; + rnemdFile_ << "# RNEMD {\n"; + + map::iterator mi; + for(mi = stringToMethod_.begin(); mi != stringToMethod_.end(); ++mi) { + if ( (*mi).second == rnemdMethod_) + rnemdFile_ << "# exchangeMethod = " << (*mi).first << "\n"; } - if (outputVx_) { - vxzLog_ << time; - for (j = 0; j < rnemdLogWidth_; j++) { - vxzLog_ << "\t" << pxzHist_[j] / mHist_[j]; - } - vxzLog_ << endl; + map::iterator fi; + for(fi = stringToFluxType_.begin(); fi != stringToFluxType_.end(); ++fi) { + if ( (*fi).second == rnemdFluxType_) + rnemdFile_ << "# fluxType = " << (*fi).first << "\n"; } - if (outputVy_) { - vyzLog_ << time; - for (j = 0; j < rnemdLogWidth_; j++) { - vyzLog_ << "\t" << pyzHist_[j] / mHist_[j]; - } - vyzLog_ << endl; + + rnemdFile_ << "# exchangeTime = " << exchangeTime_ << " fs\n"; + + rnemdFile_ << "# objectSelection = \"" + << rnemdObjectSelection_ << "\"\n"; + rnemdFile_ << "# slabWidth = " << slabWidth_ << " angstroms\n"; + rnemdFile_ << "# slabAcenter = " << slabACenter_ << " angstroms\n"; + rnemdFile_ << "# slabBcenter = " << slabBCenter_ << " angstroms\n"; + rnemdFile_ << "# }\n"; + rnemdFile_ << "#######################################################\n"; + + rnemdFile_ << "# running time = " << time << " fs\n"; + rnemdFile_ << "# target kinetic flux = " << kineticFlux_ << "\n"; + rnemdFile_ << "# target momentum flux = " << momentumFluxVector_ << "\n"; + + rnemdFile_ << "# target one-time kinetic exchange = " << kineticTarget_ + << "\n"; + rnemdFile_ << "# target one-time momentum exchange = " << momentumTarget_ + << "\n"; + + rnemdFile_ << "# actual kinetic exchange = " << kineticExchange_ << "\n"; + rnemdFile_ << "# actual momentum exchange = " << momentumExchange_ + << "\n"; + + rnemdFile_ << "# attempted exchanges: " << trialCount_ << "\n"; + rnemdFile_ << "# failed exchanges: " << failTrialCount_ << "\n"; + + + if (rnemdMethod_ == rnemdNIVS) { + rnemdFile_ << "# NIVS root-check warnings: " << failRootCount_ << "\n"; } - if (output3DTemp_) { - RealType temp; - xTempLog_ << time; - for (j = 0; j < rnemdLogWidth_; j++) { - if (outputVx_) - xTempHist_[j] -= pxzHist_[j] * pxzHist_[j] / mHist_[j]; - temp = xTempHist_[j] / (RealType)xyzTempCount_[j] - / PhysicalConstants::energyConvert / PhysicalConstants::kb; - xTempLog_ << "\t" << temp; + rnemdFile_ << "#######################################################\n"; + + + + //write title + rnemdFile_ << "#"; + for (unsigned int i = 0; i < outputMask_.size(); ++i) { + if (outputMask_[i]) { + rnemdFile_ << "\t" << data_[i].title << + "(" << data_[i].units << ")"; } - xTempLog_ << endl; - yTempLog_ << time; - for (j = 0; j < rnemdLogWidth_; j++) { - yTempLog_ << "\t" << yTempHist_[j] / (RealType)xyzTempCount_[j]; - } - yTempLog_ << endl; - zTempLog_ << time; - for (j = 0; j < rnemdLogWidth_; j++) { - zTempLog_ << "\t" << zTempHist_[j] / (RealType)xyzTempCount_[j]; - } - zTempLog_ << endl; } - if (outputRotTemp_) { - rotTempLog_ << time; - for (j = 0; j < rnemdLogWidth_; j++) { - rotTempLog_ << "\t" << rotTempHist_[j] / (RealType)rotTempCount_[j]; - } - rotTempLog_ << endl; - } - + rnemdFile_ << std::endl; + + rnemdFile_.precision(8); + + for (unsigned int j = 0; j < nBins_; j++) { + + for (unsigned int i = 0; i < outputMask_.size(); ++i) { + if (outputMask_[i]) { + if (data_[i].dataType == "RealType") + writeReal(i,j); + else if (data_[i].dataType == "Vector3d") + writeVector(i,j); + else { + sprintf( painCave.errMsg, + "RNEMD found an unknown data type for: %s ", + data_[i].title.c_str()); + painCave.isFatal = 1; + simError(); + } + } + } + rnemdFile_ << std::endl; + + } + + rnemdFile_.flush(); + rnemdFile_.close(); + #ifdef IS_MPI } #endif - - for (j = 0; j < rnemdLogWidth_; j++) { - mHist_[j] = 0.0; - } - if (outputTemp_) - for (j = 0; j < rnemdLogWidth_; j++) { - tempCount_[j] = 0; - tempHist_[j] = 0.0; - } - if (outputVx_) - for (j = 0; j < rnemdLogWidth_; j++) { - //pxzCount_[j] = 0; - pxzHist_[j] = 0.0; - } - if (outputVy_) - for (j = 0; j < rnemdLogWidth_; j++) { - //pyzCount_[j] = 0; - pyzHist_[j] = 0.0; - } - - if (output3DTemp_) - for (j = 0; j < rnemdLogWidth_; j++) { - xTempHist_[j] = 0.0; - yTempHist_[j] = 0.0; - zTempHist_[j] = 0.0; - xyzTempCount_[j] = 0; - } - if (outputRotTemp_) - for (j = 0; j < rnemdLogWidth_; j++) { - rotTempCount_[j] = 0; - rotTempHist_[j] = 0.0; - } + } + + void RNEMD::writeReal(int index, unsigned int bin) { + assert(index >=0 && index < ENDINDEX); + assert(bin >=0 && bin < nBins_); + RealType s; + + data_[index].accumulator[bin]->getAverage(s); + + if (! isinf(s) && ! isnan(s)) { + rnemdFile_ << "\t" << s; + } else{ + sprintf( painCave.errMsg, + "RNEMD detected a numerical error writing: %s for bin %d", + data_[index].title.c_str(), bin); + painCave.isFatal = 1; + simError(); + } + } + + void RNEMD::writeVector(int index, unsigned int bin) { + assert(index >=0 && index < ENDINDEX); + assert(bin >=0 && bin < nBins_); + Vector3d s; + dynamic_cast(data_[index].accumulator[bin])->getAverage(s); + if (isinf(s[0]) || isnan(s[0]) || + isinf(s[1]) || isnan(s[1]) || + isinf(s[2]) || isnan(s[2]) ) { + sprintf( painCave.errMsg, + "RNEMD detected a numerical error writing: %s for bin %d", + data_[index].title.c_str(), bin); + painCave.isFatal = 1; + simError(); + } else { + rnemdFile_ << "\t" << s[0] << "\t" << s[1] << "\t" << s[2]; + } + } }