--- branches/development/src/rnemd/RNEMD.cpp 2012/08/07 18:26:40 1773 +++ branches/development/src/rnemd/RNEMD.cpp 2013/01/09 22:02:30 1830 @@ -53,6 +53,11 @@ #include #endif +#ifdef _MSC_VER +#define isnan(x) _isnan((x)) +#define isinf(x) (!_finite(x) && !_isnan(x)) +#endif + #define HONKING_LARGE_VALUE 1.0e10 using namespace std; @@ -65,10 +70,12 @@ namespace OpenMD { failTrialCount_ = 0; failRootCount_ = 0; - int seedValue; Globals * simParams = info->getSimParams(); RNEMDParameters* rnemdParams = simParams->getRNEMDParameters(); + doRNEMD_ = rnemdParams->getUseRNEMD(); + if (!doRNEMD_) return; + stringToMethod_["Swap"] = rnemdSwap; stringToMethod_["NIVS"] = rnemdNIVS; stringToMethod_["VSS"] = rnemdVSS; @@ -77,6 +84,7 @@ namespace OpenMD { stringToFluxType_["Px"] = rnemdPx; stringToFluxType_["Py"] = rnemdPy; stringToFluxType_["Pz"] = rnemdPz; + stringToFluxType_["Pvector"] = rnemdPvector; stringToFluxType_["KE+Px"] = rnemdKePx; stringToFluxType_["KE+Py"] = rnemdKePy; stringToFluxType_["KE+Pvector"] = rnemdKePvector; @@ -98,8 +106,8 @@ namespace OpenMD { 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"); + "\tKE, Px, Py, Pz, Pvector, KE+Px, KE+Py, KE+Pvector\n" + "\tmust be set to use RNEMD\n"); painCave.isFatal = 1; painCave.severity = OPENMD_ERROR; simError(); @@ -197,6 +205,7 @@ namespace OpenMD { break; case rnemdPvector: hasCorrectFlux = hasMomentumFluxVector; + break; case rnemdKePx: case rnemdKePy: hasCorrectFlux = hasMomentumFlux && hasKineticFlux; @@ -224,8 +233,7 @@ namespace OpenMD { } if (!hasCorrectFlux) { sprintf(painCave.errMsg, - "RNEMD: The current method,\n" - "\t%s, and flux type %s\n" + "RNEMD: The current method, %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()); @@ -235,7 +243,10 @@ namespace OpenMD { } if (hasKineticFlux) { - kineticFlux_ = rnemdParams->getKineticFlux(); + // convert the kcal / mol / Angstroms^2 / fs values in the md file + // into amu / fs^3: + kineticFlux_ = rnemdParams->getKineticFlux() + * PhysicalConstants::energyConvert; } else { kineticFlux_ = 0.0; } @@ -270,6 +281,7 @@ namespace OpenMD { // do some sanity checking int selectionCount = seleMan_.getSelectionCount(); + int nIntegrable = info->getNGlobalIntegrableObjects(); if (selectionCount > nIntegrable) { @@ -288,6 +300,8 @@ namespace OpenMD { simError(); } + areaAccumulator_ = new Accumulator(); + nBins_ = rnemdParams->getOutputBins(); data_.resize(RNEMD::ENDINDEX); @@ -296,7 +310,7 @@ namespace OpenMD { z.title = "Z"; z.dataType = "RealType"; z.accumulator.reserve(nBins_); - for (unsigned int i = 0; i < nBins_; i++) + for (int i = 0; i < nBins_; i++) z.accumulator.push_back( new Accumulator() ); data_[Z] = z; outputMap_["Z"] = Z; @@ -306,17 +320,17 @@ namespace OpenMD { temperature.title = "Temperature"; temperature.dataType = "RealType"; temperature.accumulator.reserve(nBins_); - for (unsigned int i = 0; i < nBins_; i++) + for (int i = 0; i < nBins_; i++) temperature.accumulator.push_back( new Accumulator() ); data_[TEMPERATURE] = temperature; outputMap_["TEMPERATURE"] = TEMPERATURE; OutputData velocity; - velocity.units = "amu/fs"; + velocity.units = "angstroms/fs"; velocity.title = "Velocity"; velocity.dataType = "Vector3d"; velocity.accumulator.reserve(nBins_); - for (unsigned int i = 0; i < nBins_; i++) + for (int i = 0; i < nBins_; i++) velocity.accumulator.push_back( new VectorAccumulator() ); data_[VELOCITY] = velocity; outputMap_["VELOCITY"] = VELOCITY; @@ -326,7 +340,7 @@ namespace OpenMD { density.title = "Density"; density.dataType = "RealType"; density.accumulator.reserve(nBins_); - for (unsigned int i = 0; i < nBins_; i++) + for (int i = 0; i < nBins_; i++) density.accumulator.push_back( new Accumulator() ); data_[DENSITY] = density; outputMap_["DENSITY"] = DENSITY; @@ -381,8 +395,9 @@ namespace OpenMD { // dt = exchange time interval // flux = target flux - kineticTarget_ = 2.0*kineticFlux_*exchangeTime_*hmat(0,0)*hmat(1,1); - momentumTarget_ = 2.0*momentumFluxVector_*exchangeTime_*hmat(0,0)*hmat(1,1); + RealType area = currentSnap_->getXYarea(); + kineticTarget_ = 2.0 * kineticFlux_ * exchangeTime_ * area; + momentumTarget_ = 2.0 * momentumFluxVector_ * exchangeTime_ * area; // total exchange sums are zeroed out at the beginning: @@ -407,7 +422,7 @@ namespace OpenMD { } RNEMD::~RNEMD() { - + if (!doRNEMD_) return; #ifdef IS_MPI if (worldRank == 0) { #endif @@ -429,7 +444,7 @@ namespace OpenMD { } void RNEMD::doSwap() { - + if (!doRNEMD_) return; Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot(); Mat3x3d hmat = currentSnap_->getHmat(); @@ -437,7 +452,6 @@ namespace OpenMD { int selei; StuntDouble* sd; - int idx; RealType min_val; bool min_found = false; @@ -450,8 +464,6 @@ namespace OpenMD { for (sd = seleMan_.beginSelected(selei); sd != NULL; sd = seleMan_.nextSelected(selei)) { - idx = sd->getLocalIndex(); - Vector3d pos = sd->getPos(); // wrap the stuntdouble's position back into the box: @@ -488,9 +500,6 @@ namespace OpenMD { + angMom[2]*angMom[2]/I(2, 2); } } //angular momenta exchange enabled - //energyConvert temporarily disabled - //make kineticExchange_ comparable between swap & scale - //value = value * 0.5 / PhysicalConstants::energyConvert; value *= 0.5; break; case rnemdPx : @@ -532,12 +541,9 @@ namespace OpenMD { } } -#ifdef IS_MPI - int nProc, worldRank; +#ifdef IS_MPI + int worldRank = MPI::COMM_WORLD.Get_rank(); - nProc = MPI::COMM_WORLD.Get_size(); - worldRank = MPI::COMM_WORLD.Get_rank(); - bool my_min_found = min_found; bool my_max_found = max_found; @@ -728,7 +734,6 @@ namespace OpenMD { switch(rnemdFluxType_) { case rnemdKE: - cerr << "KE\n"; kineticExchange_ += max_val - min_val; break; case rnemdPx: @@ -741,7 +746,6 @@ namespace OpenMD { momentumExchange_.z() += max_val - min_val; break; default: - cerr << "default\n"; break; } } else { @@ -764,7 +768,7 @@ namespace OpenMD { } void RNEMD::doNIVS() { - + if (!doRNEMD_) return; Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot(); Mat3x3d hmat = currentSnap_->getHmat(); @@ -772,7 +776,6 @@ namespace OpenMD { int selei; StuntDouble* sd; - int idx; vector hotBin, coldBin; @@ -794,8 +797,6 @@ namespace OpenMD { for (sd = seleMan_.beginSelected(selei); sd != NULL; sd = seleMan_.nextSelected(selei)) { - idx = sd->getLocalIndex(); - Vector3d pos = sd->getPos(); // wrap the stuntdouble's position back into the box: @@ -908,8 +909,7 @@ namespace OpenMD { if ((c > 0.81) && (c < 1.21)) {//restrict scaling coefficients c = sqrt(c); - //std::cerr << "cold slab scaling coefficient: " << c << endl; - //now convert to hotBin coefficient + RealType w = 0.0; if (rnemdFluxType_ == rnemdFullKE) { x = 1.0 + px * (1.0 - c); @@ -947,8 +947,6 @@ namespace OpenMD { } } w = sqrt(w); - // std::cerr << "xh= " << x << "\tyh= " << y << "\tzh= " << z - // << "\twh= " << w << endl; for (sdi = hotBin.begin(); sdi != hotBin.end(); sdi++) { if (rnemdFluxType_ == rnemdFullKE) { vel = (*sdi)->getVel(); @@ -1213,7 +1211,7 @@ namespace OpenMD { } void RNEMD::doVSS() { - + if (!doRNEMD_) return; Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot(); RealType time = currentSnap_->getTime(); Mat3x3d hmat = currentSnap_->getHmat(); @@ -1222,7 +1220,6 @@ namespace OpenMD { int selei; StuntDouble* sd; - int idx; vector hotBin, coldBin; @@ -1237,8 +1234,6 @@ namespace OpenMD { for (sd = seleMan_.beginSelected(selei); sd != NULL; sd = seleMan_.nextSelected(selei)) { - idx = sd->getLocalIndex(); - Vector3d pos = sd->getPos(); // wrap the stuntdouble's position back into the box: @@ -1257,9 +1252,7 @@ namespace OpenMD { 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 (rnemdFluxType_ == rnemdFullKE) { @@ -1307,10 +1300,6 @@ 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; #ifdef IS_MPI MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Ph[0], 3, MPI::REALTYPE, MPI::SUM); @@ -1342,8 +1331,7 @@ 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"; + vector::iterator sdi; Vector3d vel; for (sdi = coldBin.begin(); sdi != coldBin.end(); sdi++) { @@ -1391,7 +1379,7 @@ namespace OpenMD { } void RNEMD::doRNEMD() { - + if (!doRNEMD_) return; trialCount_++; switch(rnemdMethod_) { case rnemdSwap: @@ -1410,15 +1398,16 @@ namespace OpenMD { } void RNEMD::collectData() { - + if (!doRNEMD_) return; Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot(); Mat3x3d hmat = currentSnap_->getHmat(); + areaAccumulator_->add(currentSnap_->getXYarea()); + seleMan_.setSelectionSet(evaluator_.evaluate()); - int selei; + int selei(0); StuntDouble* sd; - int idx; vector binMass(nBins_, 0.0); vector binPx(nBins_, 0.0); @@ -1442,11 +1431,10 @@ namespace OpenMD { sd != NULL; sd = mol->nextIntegrableObject(iiter)) */ + for (sd = seleMan_.beginSelected(selei); sd != NULL; - sd = seleMan_.nextSelected(selei)) { - - idx = sd->getLocalIndex(); - + sd = seleMan_.nextSelected(selei)) { + Vector3d pos = sd->getPos(); // wrap the stuntdouble's position back into the box: @@ -1454,13 +1442,14 @@ 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)] // 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(); Vector3d vel = sd->getVel(); @@ -1491,7 +1480,6 @@ namespace OpenMD { } } - #ifdef IS_MPI MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binCount[0], nBins_, MPI::INT, MPI::SUM); @@ -1518,25 +1506,31 @@ namespace OpenMD { 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; + den = binMass[i] * nBins_ * PhysicalConstants::densityConvert + / currentSnap_->getVolume() ; + + if (binCount[i] > 0) { + // only add values if there are things to add + 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: + dynamic_cast(data_[j].accumulator[i])->add(z); + break; + case TEMPERATURE: + dynamic_cast(data_[j].accumulator[i])->add(temp); + break; + case VELOCITY: + dynamic_cast(data_[j].accumulator[i])->add(vel); + break; + case DENSITY: + dynamic_cast(data_[j].accumulator[i])->add(den); + break; + } } } } @@ -1544,11 +1538,13 @@ namespace OpenMD { } void RNEMD::getStarted() { + if (!doRNEMD_) return; collectData(); writeOutputFile(); } void RNEMD::parseOutputFileFormat(const std::string& format) { + if (!doRNEMD_) return; StringTokenizer tokenizer(format, " ,;|\t\n\r"); while(tokenizer.hasMoreTokens()) { @@ -1569,6 +1565,7 @@ namespace OpenMD { } void RNEMD::writeOutputFile() { + if (!doRNEMD_) return; #ifdef IS_MPI // If we're the root node, should we print out the results @@ -1588,53 +1585,67 @@ namespace OpenMD { Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot(); RealType time = currentSnap_->getTime(); - - + RealType avgArea; + areaAccumulator_->getAverage(avgArea); + RealType Jz = kineticExchange_ / (2.0 * time * avgArea) + / PhysicalConstants::energyConvert; + Vector3d JzP = momentumExchange_ / (2.0 * time * avgArea); + 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"; + rnemdFile_ << "# exchangeMethod = \"" << (*mi).first << "\";\n"; } map::iterator fi; for(fi = stringToFluxType_.begin(); fi != stringToFluxType_.end(); ++fi) { if ( (*fi).second == rnemdFluxType_) - rnemdFile_ << "# fluxType = " << (*fi).first << "\n"; + rnemdFile_ << "# fluxType = \"" << (*fi).first << "\";\n"; } - rnemdFile_ << "# exchangeTime = " << exchangeTime_ << " fs\n"; + rnemdFile_ << "# exchangeTime = " << exchangeTime_ << ";\n"; rnemdFile_ << "# objectSelection = \"" - << rnemdObjectSelection_ << "\"\n"; - rnemdFile_ << "# slabWidth = " << slabWidth_ << " angstroms\n"; - rnemdFile_ << "# slabAcenter = " << slabACenter_ << " angstroms\n"; - rnemdFile_ << "# slabBcenter = " << slabBCenter_ << " angstroms\n"; + << rnemdObjectSelection_ << "\";\n"; + rnemdFile_ << "# slabWidth = " << slabWidth_ << ";\n"; + rnemdFile_ << "# slabAcenter = " << slabACenter_ << ";\n"; + rnemdFile_ << "# slabBcenter = " << slabBCenter_ << ";\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"; - - + rnemdFile_ << "# RNEMD report:\n"; + rnemdFile_ << "# running time = " << time << " fs\n"; + rnemdFile_ << "# target flux:\n"; + rnemdFile_ << "# kinetic = " + << kineticFlux_ / PhysicalConstants::energyConvert + << " (kcal/mol/A^2/fs)\n"; + rnemdFile_ << "# momentum = " << momentumFluxVector_ + << " (amu/A/fs^2)\n"; + rnemdFile_ << "# target one-time exchanges:\n"; + rnemdFile_ << "# kinetic = " + << kineticTarget_ / PhysicalConstants::energyConvert + << " (kcal/mol)\n"; + rnemdFile_ << "# momentum = " << momentumTarget_ + << " (amu*A/fs)\n"; + rnemdFile_ << "# actual exchange totals:\n"; + rnemdFile_ << "# kinetic = " + << kineticExchange_ / PhysicalConstants::energyConvert + << " (kcal/mol)\n"; + rnemdFile_ << "# momentum = " << momentumExchange_ + << " (amu*A/fs)\n"; + rnemdFile_ << "# actual flux:\n"; + rnemdFile_ << "# kinetic = " << Jz + << " (kcal/mol/A^2/fs)\n"; + rnemdFile_ << "# momentum = " << JzP + << " (amu/A/fs^2)\n"; + rnemdFile_ << "# exchange statistics:\n"; + rnemdFile_ << "# attempted = " << trialCount_ << "\n"; + rnemdFile_ << "# failed = " << failTrialCount_ << "\n"; if (rnemdMethod_ == rnemdNIVS) { - rnemdFile_ << "# NIVS root-check warnings: " << failRootCount_ << "\n"; + rnemdFile_ << "# NIVS root-check errors = " + << failRootCount_ << "\n"; } - rnemdFile_ << "#######################################################\n"; @@ -1645,13 +1656,15 @@ namespace OpenMD { if (outputMask_[i]) { rnemdFile_ << "\t" << data_[i].title << "(" << data_[i].units << ")"; + // add some extra tabs for column alignment + if (data_[i].dataType == "Vector3d") rnemdFile_ << "\t\t"; } } rnemdFile_ << std::endl; rnemdFile_.precision(8); - for (unsigned int j = 0; j < nBins_; j++) { + for (int j = 0; j < nBins_; j++) { for (unsigned int i = 0; i < outputMask_.size(); ++i) { if (outputMask_[i]) { @@ -1671,6 +1684,32 @@ namespace OpenMD { rnemdFile_ << std::endl; } + + rnemdFile_ << "#######################################################\n"; + rnemdFile_ << "# Standard Deviations in those quantities follow:\n"; + rnemdFile_ << "#######################################################\n"; + + + for (int j = 0; j < nBins_; j++) { + rnemdFile_ << "#"; + for (unsigned int i = 0; i < outputMask_.size(); ++i) { + if (outputMask_[i]) { + if (data_[i].dataType == "RealType") + writeRealStdDev(i,j); + else if (data_[i].dataType == "Vector3d") + writeVectorStdDev(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(); @@ -1682,11 +1721,12 @@ namespace OpenMD { } void RNEMD::writeReal(int index, unsigned int bin) { + if (!doRNEMD_) return; assert(index >=0 && index < ENDINDEX); - assert(bin >=0 && bin < nBins_); + assert(int(bin) < nBins_); RealType s; - data_[index].accumulator[bin]->getAverage(s); + dynamic_cast(data_[index].accumulator[bin])->getAverage(s); if (! isinf(s) && ! isnan(s)) { rnemdFile_ << "\t" << s; @@ -1700,8 +1740,9 @@ namespace OpenMD { } void RNEMD::writeVector(int index, unsigned int bin) { + if (!doRNEMD_) return; assert(index >=0 && index < ENDINDEX); - assert(bin >=0 && bin < nBins_); + assert(int(bin) < nBins_); Vector3d s; dynamic_cast(data_[index].accumulator[bin])->getAverage(s); if (isinf(s[0]) || isnan(s[0]) || @@ -1716,5 +1757,43 @@ namespace OpenMD { rnemdFile_ << "\t" << s[0] << "\t" << s[1] << "\t" << s[2]; } } + + void RNEMD::writeRealStdDev(int index, unsigned int bin) { + if (!doRNEMD_) return; + assert(index >=0 && index < ENDINDEX); + assert(int(bin) < nBins_); + RealType s; + + dynamic_cast(data_[index].accumulator[bin])->getStdDev(s); + + if (! isinf(s) && ! isnan(s)) { + rnemdFile_ << "\t" << s; + } else{ + sprintf( painCave.errMsg, + "RNEMD detected a numerical error writing: %s std. dev. for bin %d", + data_[index].title.c_str(), bin); + painCave.isFatal = 1; + simError(); + } + } + + void RNEMD::writeVectorStdDev(int index, unsigned int bin) { + if (!doRNEMD_) return; + assert(index >=0 && index < ENDINDEX); + assert(int(bin) < nBins_); + Vector3d s; + dynamic_cast(data_[index].accumulator[bin])->getStdDev(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 std. dev. for bin %d", + data_[index].title.c_str(), bin); + painCave.isFatal = 1; + simError(); + } else { + rnemdFile_ << "\t" << s[0] << "\t" << s[1] << "\t" << s[2]; + } + } }