--- trunk/src/integrators/RNEMD.cpp 2009/03/19 21:03:36 1330 +++ trunk/src/integrators/RNEMD.cpp 2009/04/02 19:55:37 1334 @@ -40,9 +40,12 @@ */ #include "integrators/RNEMD.hpp" +#include "math/Vector3.hpp" #include "math/SquareMatrix3.hpp" #include "primitives/Molecule.hpp" #include "primitives/StuntDouble.hpp" +#include "utils/OOPSEConstant.hpp" +#include "utils/Tuple.hpp" #ifndef IS_MPI #include "math/SeqRandNumGen.hpp" @@ -59,7 +62,7 @@ namespace oopse { namespace oopse { - RNEMD::RNEMD(SimInfo* info) : info_(info) { + RNEMD::RNEMD(SimInfo* info) : info_(info), evaluator_(info), seleMan_(info), usePeriodicBoundaryConditions_(info->getSimParams()->getUsePeriodicBoundaryConditions()) { int seedValue; Globals * simParams = info->getSimParams(); @@ -70,13 +73,18 @@ namespace oopse { stringToEnumMap_["Pz"] = rnemdPz; stringToEnumMap_["Unknown"] = rnemdUnknown; + rnemdObjectSelection_ = simParams->getRNEMD_objectSelection(); + + std::cerr << "calling evaluator with " << rnemdObjectSelection_ << "\n"; + evaluator_.loadScriptString(rnemdObjectSelection_); + std::cerr << "done"; + const std::string st = simParams->getRNEMD_swapType(); std::map::iterator i; i = stringToEnumMap_.find(st); rnemdType_ = (i == stringToEnumMap_.end()) ? RNEMD::rnemdUnknown : i->second; - set_RNEMD_swapTime(simParams->getRNEMD_swapTime()); set_RNEMD_nBins(simParams->getRNEMD_nBins()); exchangeSum_ = 0.0; @@ -105,8 +113,116 @@ namespace oopse { 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 << "exchangeSum = " << exchangeSum_ << "\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; + + std::vector > endSlice; + std::vector > midSlice; + + 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: + + 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)] + + int binNo = midBin + int(nBins_ * (pos.z()) / hmat(2,2)); + + std::cerr << "pos.z() = " << pos.z() << " bin = " << binNo << "\n"; + + // if we're in bin 0 or the middleBin + if (binNo == 0 || binNo == midBin) { + + RealType mass = sd->getMass(); + Vector3d vel = sd->getVel(); + RealType value; + + switch(rnemdType_) { + case rnemdKinetic : + + 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(); + + 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); + } else { + value += angMom[0]*angMom[0]/I(0, 0) + + angMom[1]*angMom[1]/I(1, 1) + + angMom[2]*angMom[2]/I(2, 2); + } + } + value = value * 0.5 / OOPSEConstant::energyConvert; + break; + case rnemdPx : + value = mass * vel[0]; + break; + case rnemdPy : + value = mass * vel[1]; + break; + case rnemdPz : + value = mass * vel[2]; + break; + case rnemdUnknown : + default : + break; + } + + if (binNo == 0) + endSlice.push_back( make_tuple3(value, idx, sd)); + + if (binNo == midBin) + midSlice.push_back( make_tuple3(value, idx, sd)); + } + } + + // find smallest value in endSlice: + std::sort(endSlice.begin(), endSlice.end()); + RealType min_val = endSlice[0].first; + int min_idx = endSlice[0].second; + StuntDouble* min_sd = endSlice[0].third; + + std::cerr << "smallest value = " << min_val << " idx = " << min_idx << "\n"; + + // find largest value in midSlice: + std::sort(midSlice.rbegin(), midSlice.rend()); + RealType max_val = midSlice[0].first; + int max_idx = midSlice[0].second; + StuntDouble* max_sd = midSlice[0].third; + + std::cerr << "largest value = " << max_val << " idx = " << max_idx << "\n"; + + } +}