--- branches/development/src/rnemd/RNEMD.cpp 2012/08/09 18:35:09 1777 +++ branches/development/src/rnemd/RNEMD.cpp 2013/02/20 15:39:39 1850 @@ -35,7 +35,7 @@ * * [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005). * [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006). - * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008). + * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008). * [4] Vardeman & Gezelter, in progress (2009). */ @@ -51,6 +51,11 @@ #include "utils/Tuple.hpp" #ifdef IS_MPI #include +#endif + +#ifdef _MSC_VER +#define isnan(x) _isnan((x)) +#define isinf(x) (!_finite(x) && !_isnan(x)) #endif #define HONKING_LARGE_VALUE 1.0e10 @@ -65,7 +70,6 @@ namespace OpenMD { failTrialCount_ = 0; failRootCount_ = 0; - int seedValue; Globals * simParams = info->getSimParams(); RNEMDParameters* rnemdParams = simParams->getRNEMDParameters(); @@ -277,6 +281,7 @@ namespace OpenMD { // do some sanity checking int selectionCount = seleMan_.getSelectionCount(); + int nIntegrable = info->getNGlobalIntegrableObjects(); if (selectionCount > nIntegrable) { @@ -305,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; @@ -315,7 +320,7 @@ 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; @@ -325,7 +330,7 @@ namespace OpenMD { 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; @@ -335,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; @@ -447,7 +452,6 @@ namespace OpenMD { int selei; StuntDouble* sd; - int idx; RealType min_val; bool min_found = false; @@ -460,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: @@ -539,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; @@ -777,7 +776,6 @@ namespace OpenMD { int selei; StuntDouble* sd; - int idx; vector hotBin, coldBin; @@ -799,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: @@ -913,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); @@ -952,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(); @@ -1227,7 +1220,6 @@ namespace OpenMD { int selei; StuntDouble* sd; - int idx; vector hotBin, coldBin; @@ -1242,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: @@ -1262,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) { @@ -1312,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); @@ -1347,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++) { @@ -1423,9 +1406,8 @@ namespace OpenMD { seleMan_.setSelectionSet(evaluator_.evaluate()); - int selei; + int selei(0); StuntDouble* sd; - int idx; vector binMass(nBins_, 0.0); vector binPx(nBins_, 0.0); @@ -1449,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: @@ -1468,7 +1449,7 @@ namespace OpenMD { // 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(); @@ -1499,7 +1480,6 @@ namespace OpenMD { } } - #ifdef IS_MPI MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binCount[0], nBins_, MPI::INT, MPI::SUM); @@ -1530,24 +1510,27 @@ namespace OpenMD { den = binMass[i] * nBins_ * PhysicalConstants::densityConvert / currentSnap_->getVolume() ; - 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; + 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; + } } } } @@ -1681,7 +1664,7 @@ namespace OpenMD { 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]) { @@ -1707,7 +1690,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]) { @@ -1740,10 +1723,10 @@ namespace OpenMD { void RNEMD::writeReal(int index, unsigned int bin) { if (!doRNEMD_) return; assert(index >=0 && index < ENDINDEX); - assert(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; @@ -1759,7 +1742,7 @@ namespace OpenMD { void RNEMD::writeVector(int index, unsigned int bin) { if (!doRNEMD_) return; assert(index >=0 && index < ENDINDEX); - assert(bin < nBins_); + assert(int(bin) < nBins_); Vector3d s; dynamic_cast(data_[index].accumulator[bin])->getAverage(s); if (isinf(s[0]) || isnan(s[0]) || @@ -1778,10 +1761,10 @@ namespace OpenMD { void RNEMD::writeRealStdDev(int index, unsigned int bin) { if (!doRNEMD_) return; assert(index >=0 && index < ENDINDEX); - assert(bin < nBins_); + assert(int(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; @@ -1797,7 +1780,7 @@ namespace OpenMD { void RNEMD::writeVectorStdDev(int index, unsigned int bin) { if (!doRNEMD_) return; assert(index >=0 && index < ENDINDEX); - assert(bin < nBins_); + assert(int(bin) < nBins_); Vector3d s; dynamic_cast(data_[index].accumulator[bin])->getStdDev(s); if (isinf(s[0]) || isnan(s[0]) ||