--- trunk/src/integrators/RNEMD.cpp 2009/11/25 20:02:06 1390 +++ branches/development/src/integrators/RNEMD.cpp 2011/11/22 20:38:56 1665 @@ -36,7 +36,8 @@ * [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). - * [4] Vardeman & Gezelter, in progress (2009). + * [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). + * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). */ #include @@ -52,14 +53,17 @@ #ifndef IS_MPI #include "math/SeqRandNumGen.hpp" #else +#include #include "math/ParallelRandNumGen.hpp" #endif #define HONKING_LARGE_VALUE 1.0e10 +using namespace std; namespace OpenMD { - RNEMD::RNEMD(SimInfo* info) : info_(info), evaluator_(info), seleMan_(info), usePeriodicBoundaryConditions_(info->getSimParams()->getUsePeriodicBoundaryConditions()) { + RNEMD::RNEMD(SimInfo* info) : info_(info), evaluator_(info), seleMan_(info), + usePeriodicBoundaryConditions_(info->getSimParams()->getUsePeriodicBoundaryConditions()) { failTrialCount_ = 0; failRootCount_ = 0; @@ -88,7 +92,7 @@ namespace OpenMD { if (selectionCount > nIntegrable) { sprintf(painCave.errMsg, - "RNEMD warning: The current RNEMD_objectSelection,\n" + "RNEMD: The current RNEMD_objectSelection,\n" "\t\t%s\n" "\thas resulted in %d selected objects. However,\n" "\tthe total number of integrable objects in the system\n" @@ -98,27 +102,36 @@ namespace OpenMD { rnemdObjectSelection_.c_str(), selectionCount, nIntegrable); painCave.isFatal = 0; + painCave.severity = OPENMD_WARNING; simError(); - } - const std::string st = simParams->getRNEMD_exchangeType(); + const string st = simParams->getRNEMD_exchangeType(); - std::map::iterator i; + map::iterator i; i = stringToEnumMap_.find(st); rnemdType_ = (i == stringToEnumMap_.end()) ? RNEMD::rnemdUnknown : i->second; if (rnemdType_ == rnemdUnknown) { - std::cerr << "WARNING! RNEMD Type Unknown!\n"; + 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(); } + + output3DTemp_ = false; + if (simParams->haveRNEMD_outputDimensionalTemperature()) { + output3DTemp_ = simParams->getRNEMD_outputDimensionalTemperature(); + } #ifdef IS_MPI if (worldRank == 0) { #endif - std::string rnemdFileName; - std::string xTempFileName; - std::string yTempFileName; - std::string zTempFileName; + string rnemdFileName; switch(rnemdType_) { case rnemdKineticSwap : case rnemdKineticScale : @@ -129,12 +142,6 @@ namespace OpenMD { case rnemdPy : case rnemdPyScale : rnemdFileName = "momemtum.log"; - xTempFileName = "temperatureX.log"; - yTempFileName = "temperatureY.log"; - zTempFileName = "temperatureZ.log"; - xTempLog_.open(xTempFileName.c_str()); - yTempLog_.open(yTempFileName.c_str()); - zTempLog_.open(zTempFileName.c_str()); break; case rnemdPz : case rnemdPzScale : @@ -145,6 +152,18 @@ namespace OpenMD { } rnemdLog_.open(rnemdFileName.c_str()); + string xTempFileName; + string yTempFileName; + string zTempFileName; + if (output3DTemp_) { + xTempFileName = "temperatureX.log"; + yTempFileName = "temperatureY.log"; + zTempFileName = "temperatureZ.log"; + xTempLog_.open(xTempFileName.c_str()); + yTempLog_.open(yTempFileName.c_str()); + zTempLog_.open(zTempFileName.c_str()); + } + #ifdef IS_MPI } #endif @@ -152,21 +171,35 @@ namespace OpenMD { 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; + } + } else { + zShift_ = 0.0; + } + //cerr << "we have zShift_ = " << zShift_ << "\n"; + //shift slabs by half slab width, might be useful in heterogeneous systems + //set to 0.0 if not using it; can NOT be used in status output yet if (simParams->haveRNEMD_logWidth()) { - rnemdLogWidth_ = simParams->getRNEMD_logWidth(); - if (rnemdLogWidth_ != nBins_ || rnemdLogWidth_ != midBin_ + 1) { - std::cerr << "WARNING! RNEMD_logWidth has abnormal value!\n"; - std::cerr << "Automaically set back to default.\n"; + 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_; - } + }*/ } else { - rnemdLogWidth_ = nBins_; + set_RNEMD_logWidth(nBins_); } valueHist_.resize(rnemdLogWidth_, 0.0); valueCount_.resize(rnemdLogWidth_, 0); xTempHist_.resize(rnemdLogWidth_, 0.0); yTempHist_.resize(rnemdLogWidth_, 0.0); zTempHist_.resize(rnemdLogWidth_, 0.0); + xyzTempCount_.resize(rnemdLogWidth_, 0); set_RNEMD_exchange_total(0.0); if (simParams->haveRNEMD_targetFlux()) { @@ -194,15 +227,28 @@ namespace OpenMD { RNEMD::~RNEMD() { delete randNumGen_; - - std::cerr << "total fail trials: " << failTrialCount_ << "\n"; + #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(); + rnemdLog_.close(); - if (rnemdType_ == rnemdKineticScale || rnemdType_ == rnemdPxScale || rnemdType_ == rnemdPyScale) - std::cerr<< "total root-checking warnings: " << failRootCount_ << "\n"; - if (rnemdType_ == rnemdPx || rnemdType_ == rnemdPxScale || rnemdType_ == rnemdPy || rnemdType_ == rnemdPyScale) { + 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(); @@ -246,7 +292,7 @@ namespace OpenMD { // 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) + 0.5)) % nBins_; + int binNo = int(nBins_ * (pos.z() / hmat(2,2) + zShift_ + 0.5)) % nBins_; // if we're in bin 0 or the middleBin @@ -261,22 +307,24 @@ namespace OpenMD { value = mass * (vel[0]*vel[0] + vel[1]*vel[1] + vel[2]*vel[2]); - if (sd->isDirectional()) { + /* + 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); + 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); } else { - value += angMom[0]*angMom[0]/I(0, 0) - + angMom[1]*angMom[1]/I(1, 1) - + angMom[2]*angMom[2]/I(2, 2); + value += angMom[0]*angMom[0]/I(0, 0) + + angMom[1]*angMom[1]/I(1, 1) + + angMom[2]*angMom[2]/I(2, 2); } - } + } no exchange of angular momenta + */ //make exchangeSum_ comparable between swap & scale //temporarily without using energyConvert //value = value * 0.5 / PhysicalConstants::energyConvert; @@ -331,13 +379,9 @@ namespace OpenMD { bool my_max_found = max_found; // Even if we didn't find a minimum, did someone else? - MPI::COMM_WORLD.Allreduce(&my_min_found, &min_found, - 1, MPI::BOOL, MPI::LAND); - + MPI::COMM_WORLD.Allreduce(&my_min_found, &min_found, 1, MPI::BOOL, MPI::LOR); // Even if we didn't find a maximum, did someone else? - MPI::COMM_WORLD.Allreduce(&my_max_found, &max_found, - 1, MPI::BOOL, MPI::LAND); - + MPI::COMM_WORLD.Allreduce(&my_max_found, &max_found, 1, MPI::BOOL, MPI::LOR); struct { RealType val; int rank; @@ -373,14 +417,14 @@ namespace OpenMD { #endif if (max_found && min_found) { - if (min_val< max_val) { + if (min_val < max_val) { #ifdef IS_MPI if (max_vals.rank == worldRank && min_vals.rank == worldRank) { // I have both maximum and minimum, so proceed like a single // processor version: #endif - // objects to be swapped: velocity & angular velocity + // objects to be swapped: velocity ONLY Vector3d min_vel = min_sd->getVel(); Vector3d max_vel = max_sd->getVel(); RealType temp_vel; @@ -389,12 +433,14 @@ namespace OpenMD { case rnemdKineticSwap : min_sd->setVel(max_vel); max_sd->setVel(min_vel); - if (min_sd->isDirectional() && max_sd->isDirectional()) { + /* + if (min_sd->isDirectional() && max_sd->isDirectional()) { Vector3d min_angMom = min_sd->getJ(); Vector3d max_angMom = max_sd->getJ(); min_sd->setJ(max_angMom); max_sd->setJ(min_angMom); - } + } no angular momentum exchange + */ break; case rnemdPx : temp_vel = min_vel.x(); @@ -438,20 +484,22 @@ namespace OpenMD { switch(rnemdType_) { case rnemdKineticSwap : max_sd->setVel(min_vel); - + //no angular momentum exchange for now + /* if (max_sd->isDirectional()) { Vector3d min_angMom; Vector3d max_angMom = max_sd->getJ(); - + // point-to-point swap of the angular momentum vector MPI::COMM_WORLD.Sendrecv(max_angMom.getArrayPointer(), 3, MPI::REALTYPE, min_vals.rank, 1, min_angMom.getArrayPointer(), 3, MPI::REALTYPE, min_vals.rank, 1, status); - + max_sd->setJ(min_angMom); - } + } + */ break; case rnemdPx : max_vel.x() = min_vel.x(); @@ -484,20 +532,22 @@ namespace OpenMD { switch(rnemdType_) { case rnemdKineticSwap : min_sd->setVel(max_vel); - + // no angular momentum exchange for now + /* if (min_sd->isDirectional()) { Vector3d min_angMom = min_sd->getJ(); Vector3d max_angMom; - + // point-to-point swap of the angular momentum vector MPI::COMM_WORLD.Sendrecv(min_angMom.getArrayPointer(), 3, MPI::REALTYPE, max_vals.rank, 1, max_angMom.getArrayPointer(), 3, MPI::REALTYPE, max_vals.rank, 1, status); - + min_sd->setJ(max_angMom); } + */ break; case rnemdPx : min_vel.x() = max_vel.x(); @@ -517,13 +567,21 @@ namespace OpenMD { } #endif exchangeSum_ += max_val - min_val; - } else { - std::cerr << "exchange NOT performed!\nmin_val > max_val.\n"; + } else { + sprintf(painCave.errMsg, + "RNEMD: exchange NOT performed because min_val > max_val\n"); + painCave.isFatal = 0; + painCave.severity = OPENMD_INFO; + simError(); failTrialCount_++; } } else { - std::cerr << "exchange NOT performed!\n"; - std::cerr << "at least one of the two slabs empty.\n"; + sprintf(painCave.errMsg, + "RNEMD: exchange NOT performed because at least one\n" + "\tof the two slabs is empty\n"); + painCave.isFatal = 0; + painCave.severity = OPENMD_INFO; + simError(); failTrialCount_++; } @@ -540,7 +598,7 @@ namespace OpenMD { StuntDouble* sd; int idx; - std::vector hotBin, coldBin; + vector hotBin, coldBin; RealType Phx = 0.0; RealType Phy = 0.0; @@ -570,7 +628,7 @@ namespace OpenMD { // 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) + 0.5)) % nBins_; + 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_) { @@ -629,20 +687,20 @@ namespace OpenMD { RealType a000, a110, c0, a001, a111, b01, b11, c1, c; switch(rnemdType_) { case rnemdKineticScale : - /*used hotBin coeff's & only scale x & y dimensions + // used hotBin coeff's & only scale x & y dimensions + /* RealType px = Phx / Pcx; RealType py = Phy / Pcy; a110 = Khy; c0 = - Khx - Khy - targetFlux_; a000 = Khx; - a111 = Kcy * py * py + a111 = Kcy * py * py; b11 = -2.0 * Kcy * py * (1.0 + py); c1 = Kcy * py * (2.0 + py) + Kcx * px * ( 2.0 + px) + targetFlux_; b01 = -2.0 * Kcx * px * (1.0 + px); a001 = Kcx * px * px; - */ - - //scale all three dimensions, let x = y + */ + //scale all three dimensions, let c_x = c_y a000 = Kcx + Kcy; a110 = Kcz; c0 = targetFlux_ - Kcx - Kcy - Kcz; @@ -651,7 +709,7 @@ namespace OpenMD { 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) - targetFlux_; break; case rnemdPxScale : c = 1 - targetFlux_ / Pcx; @@ -663,7 +721,7 @@ namespace OpenMD { b01 = -2.0 * Khy * py * (1.0 + py); b11 = -2.0 * Khz * pz * (1.0 + pz); c1 = Khy * py * (2.0 + py) + Khz * pz * (2.0 + pz) - + Khx * (fastpow(c * px - px - 1.0, 2) - 1.0); + + Khx * (fastpow(c * px - px - 1.0, 2) - 1.0); break; case rnemdPyScale : c = 1 - targetFlux_ / Pcy; @@ -675,9 +733,20 @@ namespace OpenMD { b01 = -2.0 * Khx * px * (1.0 + px); b11 = -2.0 * Khz * pz * (1.0 + pz); c1 = Khx * px * (2.0 + px) + Khz * pz * (2.0 + pz) - + Khy * (fastpow(c * py - py - 1.0, 2) - 1.0); + + 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; + a000 = Kcx; + a110 = Kcy; + c0 = Kcz * c * c - Kcx - Kcy - Kcz; + a001 = px * px * Khx; + a111 = py * py * Khy; + b01 = -2.0 * Khx * px * (1.0 + px); + b11 = -2.0 * Khy * py * (1.0 + py); + c1 = Khx * px * (2.0 + px) + Khy * py * (2.0 + py) + + Khz * (fastpow(c * pz - pz - 1.0, 2) - 1.0); + break; default : break; } @@ -712,16 +781,19 @@ namespace OpenMD { poly.setCoefficient(2, u2); poly.setCoefficient(1, u1); poly.setCoefficient(0, u0); - std::vector realRoots = poly.FindRealRoots(); + vector realRoots = poly.FindRealRoots(); - std::vector::iterator ri; + vector::iterator ri; RealType r1, r2, alpha0; - std::vector > rps; + vector > rps; for (ri = realRoots.begin(); ri !=realRoots.end(); ri++) { r2 = *ri; //check if FindRealRoots() give the right answer if ( fabs(u0 + r2 * (u1 + r2 * (u2 + r2 * (u3 + r2 * u4)))) > 1e-6 ) { - std::cerr << "WARNING! eq solvers might have mistakes!\n"; + sprintf(painCave.errMsg, + "RNEMD Warning: polynomial solve seems to have an error!"); + painCave.isFatal = 0; + simError(); failRootCount_++; } //might not be useful w/o rescaling coefficients @@ -729,21 +801,21 @@ namespace OpenMD { if (alpha0 >= 0.0) { r1 = sqrt(alpha0 / a000); if (fabs(c1 + r1 * (b01 + r1 * a001) + r2 * (b11 + r2 * a111)) < 1e-6) - { rps.push_back(std::make_pair(r1, r2)); } + { rps.push_back(make_pair(r1, r2)); } if (r1 > 1e-6) { //r1 non-negative r1 = -r1; if (fabs(c1 + r1 * (b01 + r1 * a001) + r2 * (b11 + r2 * a111)) <1e-6) - { rps.push_back(std::make_pair(r1, r2)); } + { rps.push_back(make_pair(r1, r2)); } } } } - // Consider combininig together the solving pair part w/ the searching + // Consider combining together the solving pair part w/ the searching // best solution part so that we don't need the pairs vector if (!rps.empty()) { RealType smallestDiff = HONKING_LARGE_VALUE; RealType diff; - std::pair bestPair = std::make_pair(1.0, 1.0); - std::vector >::iterator rpi; + pair bestPair = make_pair(1.0, 1.0); + vector >::iterator rpi; for (rpi = rps.begin(); rpi != rps.end(); rpi++) { r1 = (*rpi).first; r2 = (*rpi).second; @@ -762,6 +834,8 @@ namespace OpenMD { + fastpow(r1 * r1 / r2 / r2 - Kcz/Kcx, 2); break; case rnemdPzScale : + diff = fastpow(1.0 - r1, 2) + fastpow(1.0 - r2, 2) + + fastpow(r1 * r1 / r2 / r2 - Kcy/Kcx, 2); default : break; } @@ -773,34 +847,42 @@ namespace OpenMD { #ifdef IS_MPI if (worldRank == 0) { #endif - std::cerr << "we choose r1 = " << bestPair.first - << " and r2 = " << bestPair.second << "\n"; + 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 - + RealType x, y, z; - switch(rnemdType_) { - case rnemdKineticScale : - x = bestPair.first; - y = bestPair.first; - z = bestPair.second; - break; - case rnemdPxScale : - x = c; - y = bestPair.first; - z = bestPair.second; - break; - case rnemdPyScale : - x = bestPair.first; - y = c; - z = bestPair.second; - break; - case rnemdPzScale : - default : - break; - } - std::vector::iterator sdi; + switch(rnemdType_) { + case rnemdKineticScale : + x = bestPair.first; + y = bestPair.first; + z = bestPair.second; + break; + case rnemdPxScale : + x = c; + y = bestPair.first; + z = bestPair.second; + break; + case rnemdPyScale : + x = bestPair.first; + y = c; + z = bestPair.second; + break; + case rnemdPzScale : + x = bestPair.first; + y = bestPair.second; + z = c; + break; + default : + break; + } + vector::iterator sdi; Vector3d vel; for (sdi = coldBin.begin(); sdi != coldBin.end(); sdi++) { vel = (*sdi)->getVel(); @@ -823,7 +905,11 @@ namespace OpenMD { exchangeSum_ += targetFlux_; //we may want to check whether the exchange has been successful } else { - std::cerr << "exchange NOT performed!\n"; + sprintf(painCave.errMsg, + "RNEMD: exchange NOT performed!\n"); + painCave.isFatal = 0; + painCave.severity = OPENMD_INFO; + simError(); failTrialCount_++; } @@ -861,6 +947,20 @@ namespace OpenMD { StuntDouble* sd; int idx; + // alternative approach, track all molecules instead of only those + // selected for scaling/swapping: + /* + SimInfo::MoleculeIterator miter; + vector::iterator iiter; + Molecule* mol; + StuntDouble* integrableObject; + 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)) + */ for (sd = seleMan_.beginSelected(selei); sd != NULL; sd = seleMan_.nextSelected(selei)) { @@ -876,12 +976,14 @@ namespace OpenMD { // 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) + 0.5)) % nBins_; - + int binNo = int(rnemdLogWidth_ * (pos.z() / hmat(2,2) + 0.5)) % + rnemdLogWidth_; + // no symmetrization allowed due to arbitary rnemdLogWidth_ value + /* if (rnemdLogWidth_ == midBin_ + 1) if (binNo > midBin_) binNo = nBins_ - binNo; - + */ RealType mass = sd->getMass(); Vector3d vel = sd->getVel(); RealType value; @@ -891,8 +993,7 @@ namespace OpenMD { case rnemdKineticSwap : case rnemdKineticScale : - value = mass * (vel[0]*vel[0] + vel[1]*vel[1] + - vel[2]*vel[2]); + value = mass * (vel[0]*vel[0] + vel[1]*vel[1] + vel[2]*vel[2]); valueCount_[binNo] += 3; if (sd->isDirectional()) { @@ -905,9 +1006,9 @@ namespace OpenMD { int k = (i + 2) % 3; value += angMom[j] * angMom[j] / I(j, j) + angMom[k] * angMom[k] / I(k, k); - + valueCount_[binNo] +=2; - + } else { value += angMom[0]*angMom[0]/I(0, 0) + angMom[1]*angMom[1]/I(1, 1) @@ -916,21 +1017,12 @@ namespace OpenMD { } } value = value / PhysicalConstants::energyConvert / PhysicalConstants::kb; - + break; case rnemdPx : case rnemdPxScale : value = mass * vel[0]; valueCount_[binNo]++; - xVal = mass * vel.x() * vel.x() / PhysicalConstants::energyConvert - / PhysicalConstants::kb; - yVal = mass * vel.y() * vel.y() / PhysicalConstants::energyConvert - / PhysicalConstants::kb; - zVal = mass * vel.z() * vel.z() / PhysicalConstants::energyConvert - / PhysicalConstants::kb; - xTempHist_[binNo] += xVal; - yTempHist_[binNo] += yVal; - zTempHist_[binNo] += zVal; break; case rnemdPy : case rnemdPyScale : @@ -939,22 +1031,40 @@ namespace OpenMD { break; case rnemdPz : case rnemdPzScale : - value = mass * vel[2]; + value = pos.z(); //temporarily for homogeneous systems ONLY valueCount_[binNo]++; break; case rnemdUnknown : default : + value = 1.0; + valueCount_[binNo]++; break; } valueHist_[binNo] += value; - } + if (output3DTemp_) { + xVal = mass * vel.x() * vel.x() / PhysicalConstants::energyConvert + / PhysicalConstants::kb; + yVal = mass * vel.y() * vel.y() / PhysicalConstants::energyConvert + / PhysicalConstants::kb; + zVal = mass * vel.z() * vel.z() / PhysicalConstants::energyConvert + / PhysicalConstants::kb; + xTempHist_[binNo] += xVal; + yTempHist_[binNo] += yVal; + zTempHist_[binNo] += zVal; + xyzTempCount_[binNo]++; + } + } } void RNEMD::getStarted() { - Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot(); - Stats& stat = currentSnap_->statData; - stat[Stats::RNEMD_EXCHANGE_TOTAL] = exchangeSum_; + collectData(); + /* now should be able to 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_; + */ + getStatus(); } void RNEMD::getStatus() { @@ -965,8 +1075,8 @@ namespace OpenMD { stat[Stats::RNEMD_EXCHANGE_TOTAL] = exchangeSum_; //or to be more meaningful, define another item as exchangeSum_ / time + int j; - #ifdef IS_MPI // all processors have the same number of bins, and STL vectors pack their @@ -976,44 +1086,55 @@ namespace OpenMD { rnemdLogWidth_, MPI::REALTYPE, MPI::SUM); MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &valueCount_[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 we're the root node, should we print out the results int worldRank = MPI::COMM_WORLD.Get_rank(); if (worldRank == 0) { #endif - int j; rnemdLog_ << time; for (j = 0; j < rnemdLogWidth_; j++) { rnemdLog_ << "\t" << valueHist_[j] / (RealType)valueCount_[j]; - valueHist_[j] = 0.0; } rnemdLog_ << "\n"; - if (rnemdType_ == rnemdPx || rnemdType_ == rnemdPxScale ) { + if (output3DTemp_) { xTempLog_ << time; for (j = 0; j < rnemdLogWidth_; j++) { - xTempLog_ << "\t" << xTempHist_[j] / (RealType)valueCount_[j]; - xTempHist_[j] = 0.0; + xTempLog_ << "\t" << xTempHist_[j] / (RealType)xyzTempCount_[j]; } xTempLog_ << "\n"; yTempLog_ << time; for (j = 0; j < rnemdLogWidth_; j++) { - yTempLog_ << "\t" << yTempHist_[j] / (RealType)valueCount_[j]; - yTempHist_[j] = 0.0; + yTempLog_ << "\t" << yTempHist_[j] / (RealType)xyzTempCount_[j]; } yTempLog_ << "\n"; zTempLog_ << time; for (j = 0; j < rnemdLogWidth_; j++) { - zTempLog_ << "\t" << zTempHist_[j] / (RealType)valueCount_[j]; - zTempHist_[j] = 0.0; + zTempLog_ << "\t" << zTempHist_[j] / (RealType)xyzTempCount_[j]; } zTempLog_ << "\n"; } - for (j = 0; j < rnemdLogWidth_; j++) valueCount_[j] = 0; #ifdef IS_MPI - } + } #endif - - + for (j = 0; j < rnemdLogWidth_; j++) { + valueCount_[j] = 0; + valueHist_[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; + } } - }