--- branches/development/src/rnemd/RNEMD.cpp 2012/08/09 15:52:59 1776 +++ trunk/src/rnemd/RNEMD.cpp 2012/08/31 21:16:10 1793 @@ -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,7 +70,6 @@ namespace OpenMD { failTrialCount_ = 0; failRootCount_ = 0; - int seedValue; Globals * simParams = info->getSimParams(); RNEMDParameters* rnemdParams = simParams->getRNEMDParameters(); @@ -80,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; @@ -101,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(); @@ -200,6 +205,7 @@ namespace OpenMD { break; case rnemdPvector: hasCorrectFlux = hasMomentumFluxVector; + break; case rnemdKePx: case rnemdKePy: hasCorrectFlux = hasMomentumFlux && hasKineticFlux; @@ -227,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()); @@ -238,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; } @@ -301,7 +309,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; @@ -311,17 +319,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; @@ -331,7 +339,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; @@ -443,7 +451,6 @@ namespace OpenMD { int selei; StuntDouble* sd; - int idx; RealType min_val; bool min_found = false; @@ -456,8 +463,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: @@ -494,9 +499,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 : @@ -538,12 +540,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; @@ -734,7 +733,6 @@ namespace OpenMD { switch(rnemdFluxType_) { case rnemdKE: - cerr << "KE\n"; kineticExchange_ += max_val - min_val; break; case rnemdPx: @@ -747,7 +745,6 @@ namespace OpenMD { momentumExchange_.z() += max_val - min_val; break; default: - cerr << "default\n"; break; } } else { @@ -778,7 +775,6 @@ namespace OpenMD { int selei; StuntDouble* sd; - int idx; vector hotBin, coldBin; @@ -800,8 +796,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: @@ -1228,7 +1222,6 @@ namespace OpenMD { int selei; StuntDouble* sd; - int idx; vector hotBin, coldBin; @@ -1243,8 +1236,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: @@ -1426,7 +1417,6 @@ namespace OpenMD { int selei; StuntDouble* sd; - int idx; vector binMass(nBins_, 0.0); vector binPx(nBins_, 0.0); @@ -1451,10 +1441,8 @@ namespace OpenMD { sd = mol->nextIntegrableObject(iiter)) */ for (sd = seleMan_.beginSelected(selei); sd != NULL; - sd = seleMan_.nextSelected(selei)) { + sd = seleMan_.nextSelected(selei)) { - idx = sd->getLocalIndex(); - Vector3d pos = sd->getPos(); // wrap the stuntdouble's position back into the box: @@ -1527,7 +1515,10 @@ namespace OpenMD { vel.x() = binPx[i] / binMass[i]; vel.y() = binPy[i] / binMass[i]; vel.z() = binPz[i] / binMass[i]; - den = binCount[i] * nBins_ / currentSnap_->getVolume(); + + den = binMass[i] * nBins_ * PhysicalConstants::densityConvert + / currentSnap_->getVolume() ; + temp = 2.0 * binKE[i] / (binDOF[i] * PhysicalConstants::kb * PhysicalConstants::energyConvert); @@ -1535,16 +1526,16 @@ namespace OpenMD { if(outputMask_[j]) { switch(j) { case Z: - (data_[j].accumulator[i])->add(z); + dynamic_cast(data_[j].accumulator[i])->add(z); break; case TEMPERATURE: - data_[j].accumulator[i]->add(temp); + dynamic_cast(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); + dynamic_cast(data_[j].accumulator[i])->add(den); break; } } @@ -1602,7 +1593,8 @@ namespace OpenMD { RealType time = currentSnap_->getTime(); RealType avgArea; areaAccumulator_->getAverage(avgArea); - RealType Jz = kineticExchange_ / (2.0 * time * avgArea); + RealType Jz = kineticExchange_ / (2.0 * time * avgArea) + / PhysicalConstants::energyConvert; Vector3d JzP = momentumExchange_ / (2.0 * time * avgArea); rnemdFile_ << "#######################################################\n"; @@ -1631,17 +1623,28 @@ namespace OpenMD { rnemdFile_ << "# RNEMD report:\n"; rnemdFile_ << "# running time = " << time << " fs\n"; rnemdFile_ << "# target flux:\n"; - rnemdFile_ << "# kinetic = " << kineticFlux_ << "\n"; - rnemdFile_ << "# momentum = " << momentumFluxVector_ << "\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_ << "\n"; - rnemdFile_ << "# momentum = " << momentumTarget_ << "\n"; - rnemdFile_ << "# actual exchange totals:\n"; - rnemdFile_ << "# kinetic = " << kineticExchange_ << "\n"; - rnemdFile_ << "# momentum = " << momentumExchange_ << "\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 << "\n"; - rnemdFile_ << "# momentum = " << JzP << "\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"; @@ -1659,13 +1662,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]) { @@ -1691,7 +1696,7 @@ namespace OpenMD { rnemdFile_ << "#######################################################\n"; - for (unsigned int j = 0; j < nBins_; j++) { + for (int j = 0; j < nBins_; j++) { rnemdFile_ << "#"; for (unsigned int i = 0; i < outputMask_.size(); ++i) { if (outputMask_[i]) { @@ -1727,7 +1732,7 @@ namespace OpenMD { assert(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; @@ -1765,7 +1770,7 @@ namespace OpenMD { assert(bin < nBins_); RealType s; - data_[index].accumulator[bin]->getStdDev(s); + dynamic_cast(data_[index].accumulator[bin])->getStdDev(s); if (! isinf(s) && ! isnan(s)) { rnemdFile_ << "\t" << s;