--- trunk/src/integrators/RNEMD.cpp 2009/04/23 18:22:30 1338 +++ trunk/src/integrators/RNEMD.cpp 2009/05/08 19:47:05 1341 @@ -53,12 +53,6 @@ #include "math/ParallelRandNumGen.hpp" #endif -/* Remove me after testing*/ -/* -#include -#include -*/ -/*End remove me*/ namespace oopse { @@ -74,10 +68,30 @@ namespace oopse { stringToEnumMap_["Unknown"] = rnemdUnknown; rnemdObjectSelection_ = simParams->getRNEMD_objectSelection(); - - std::cerr << "calling evaluator with " << rnemdObjectSelection_ << "\n"; evaluator_.loadScriptString(rnemdObjectSelection_); - std::cerr << "done\n"; + seleMan_.setSelectionSet(evaluator_.evaluate()); + + + // do some sanity checking + + int selectionCount = seleMan_.getSelectionCount(); + int nIntegrable = info->getNGlobalIntegrableObjects(); + + if (selectionCount > nIntegrable) { + sprintf(painCave.errMsg, + "RNEMD warning: 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" + "\tis only %d. This is almost certainly not what you want\n" + "\tto do. A likely cause of this is forgetting the _RB_0\n" + "\tselector in the selection script!\n", + rnemdObjectSelection_.c_str(), + selectionCount, nIntegrable); + painCave.isFatal = 0; + simError(); + + } const std::string st = simParams->getRNEMD_swapType(); @@ -89,7 +103,6 @@ namespace oopse { set_RNEMD_nBins(simParams->getRNEMD_nBins()); exchangeSum_ = 0.0; counter_ = 0; //added by shenyu - //profile_.open("profile", std::ios::out); #ifndef IS_MPI if (simParams->haveSeed()) { @@ -110,27 +123,16 @@ namespace oopse { RNEMD::~RNEMD() { delete randNumGen_; - //profile_.close(); } void RNEMD::doSwap() { - //std::cerr << "in RNEMD!\n"; - //std::cerr << "nBins = " << nBins_ << "\n"; int midBin = nBins_ / 2; - //std::cerr << "midBin = " << midBin << "\n"; - //std::cerr << "swapTime = " << swapTime_ << "\n"; - //std::cerr << "swapType = " << rnemdType_ << "\n"; - //std::cerr << "selection = " << rnemdObjectSelection_ << "\n"; Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot(); Mat3x3d hmat = currentSnap_->getHmat(); - //std::cerr << "hmat = " << hmat << "\n"; - seleMan_.setSelectionSet(evaluator_.evaluate()); - //std::cerr << "selectionCount = " << seleMan_.getSelectionCount() << "\n\n"; - int selei; StuntDouble* sd; int idx; @@ -158,9 +160,8 @@ namespace oopse { // which bin is this stuntdouble in? // wrapped positions are in the range [-0.5*hmat(2,2), +0.5*hmat(2,2)] - int binNo = midBin + int(nBins_ * (pos.z()) / hmat(2,2)); + int binNo = int(nBins_ * (pos.z() / hmat(2,2) + 0.5)) % nBins_; - //std::cerr << "pos.z() = " << pos.z() << " bin = " << binNo << "\n"; // if we're in bin 0 or the middleBin if (binNo == 0 || binNo == midBin) { @@ -174,7 +175,6 @@ namespace oopse { value = mass * (vel[0]*vel[0] + vel[1]*vel[1] + vel[2]*vel[2]); - if (sd->isDirectional()) { Vector3d angMom = sd->getJ(); Mat3x3d I = sd->getI(); @@ -232,26 +232,27 @@ namespace oopse { } } } - //std::cerr << "smallest value = " << min_val << "\n"; - //std::cerr << "largest value = " << max_val << "\n"; - + // missing: swap information in parallel if (max_found && min_found) { if (min_val< max_val) { + Vector3d min_vel = min_sd->getVel(); Vector3d max_vel = max_sd->getVel(); RealType temp_vel; + switch(rnemdType_) { case rnemdKinetic : min_sd->setVel(max_vel); - max_sd->setVel(min_vel); + max_sd->setVel(min_vel); + 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); - } + } break; case rnemdPx : temp_vel = min_vel.x(); @@ -285,42 +286,25 @@ namespace oopse { } else { std::cerr << "exchange NOT performed.\none of the two slabs empty.\n"; } - std::cerr << "exchangeSum = " << exchangeSum_ << "\n"; } void RNEMD::getStatus() { - //std::cerr << "in RNEMD!\n"; - //std::cerr << "nBins = " << nBins_ << "\n"; - int midBin = nBins_ / 2; - //std::cerr << "midBin = " << midBin << "\n"; - //std::cerr << "swapTime = " << swapTime_ << "\n"; - //std::cerr << "exchangeSum = " << exchangeSum_ << "\n"; - //std::cerr << "swapType = " << rnemdType_ << "\n"; - //std::cerr << "selection = " << rnemdObjectSelection_ << "\n"; Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot(); Mat3x3d hmat = currentSnap_->getHmat(); + Stats& stat = currentSnap_->statData; - //std::cerr << "hmat = " << hmat << "\n"; + stat[Stats::RNEMD_SWAP_TOTAL] = exchangeSum_; seleMan_.setSelectionSet(evaluator_.evaluate()); - //std::cerr << "selectionCount = " << seleMan_.getSelectionCount() << "\n\n"; - int selei; StuntDouble* sd; int idx; - /* - RealType min_val; - bool min_found = false; - StuntDouble* min_sd; - RealType max_val; - bool max_found = false; - StuntDouble* max_sd; - */ - std::vector valueHist; - std::vector valueCount; + std::vector valueHist; // keeps track of what's being averaged + std::vector valueCount; // keeps track of the number of degrees of + // freedom being averaged valueHist.resize(nBins_); valueCount.resize(nBins_); //do they initialize themselves to zero automatically? @@ -339,7 +323,7 @@ namespace oopse { // which bin is this stuntdouble in? // wrapped positions are in the range [-0.5*hmat(2,2), +0.5*hmat(2,2)] - int binNo = midBin + int(nBins_ * (pos.z()) / hmat(2,2)); + int binNo = int(nBins_ * (pos.z() / hmat(2,2) + 0.5)) % nBins_; //std::cerr << "pos.z() = " << pos.z() << " bin = " << binNo << "\n"; @@ -354,8 +338,12 @@ namespace oopse { value = mass * (vel[0]*vel[0] + vel[1]*vel[1] + vel[2]*vel[2]); + valueCount[binNo] += 3; + //std::cerr <<"starting value = " << value << "\n"; if (sd->isDirectional()) { + //std::cerr << "angMom calculated.\n"; Vector3d angMom = sd->getJ(); + //std::cerr << "current angMom: " << angMom << "\n"; Mat3x3d I = sd->getI(); if (sd->isLinear()) { @@ -364,24 +352,35 @@ namespace oopse { int k = (i + 2) % 3; value += angMom[j] * angMom[j] / I(j, j) + angMom[k] * angMom[k] / I(k, k); - } else { + + valueCount[binNo] +=2; + + } else { + //std::cerr << "non-linear molecule.\n"; value += angMom[0]*angMom[0]/I(0, 0) + angMom[1]*angMom[1]/I(1, 1) + angMom[2]*angMom[2]/I(2, 2); + valueCount[binNo] +=3; + } } - //std::cerr <<"this value = " << value << "\n"; - value = value * 0.5 / OOPSEConstant::energyConvert; - //std::cerr <<"this value = " << value << "\n"; + //std::cerr <<"total value = " << value << "\n"; + //value *= 0.5 / OOPSEConstant::energyConvert; // get it in kcal / mol + //value *= 2.0 / OOPSEConstant::kb; // convert to temperature + value = value / OOPSEConstant::energyConvert / OOPSEConstant::kb; + //std::cerr <<"value = " << value << "\n"; break; case rnemdPx : value = mass * vel[0]; + valueCount[binNo]++; break; case rnemdPy : value = mass * vel[1]; + valueCount[binNo]++; break; case rnemdPz : value = mass * vel[2]; + valueCount[binNo]++; break; case rnemdUnknown : default : @@ -389,7 +388,6 @@ namespace oopse { } //std::cerr << "bin = " << binNo << " value = " << value ; valueHist[binNo] += value; - valueCount[binNo]++; //std::cerr << " hist = " << valueHist[binNo] << " count = " << valueCount[binNo] << "\n"; }