--- branches/development/src/rnemd/RNEMD.cpp 2013/04/09 19:45:54 1861 +++ branches/development/src/rnemd/RNEMD.cpp 2013/05/15 15:09:35 1874 @@ -71,7 +71,8 @@ namespace OpenMD { RNEMD::RNEMD(SimInfo* info) : info_(info), evaluator_(info), seleMan_(info), evaluatorA_(info), seleManA_(info), commonA_(info), evaluatorB_(info), - seleManB_(info), commonB_(info), + seleManB_(info), commonB_(info), + hasData_(false), hasDividingArea_(false), usePeriodicBoundaryConditions_(info->getSimParams()->getUsePeriodicBoundaryConditions()) { trialCount_ = 0; @@ -556,7 +557,7 @@ namespace OpenMD { slabWidth_ = hmat(2,2) / 10.0; if (hasSlabBCenter) - slabBCenter_ = rnemdParams->getSlabACenter(); + slabBCenter_ = rnemdParams->getSlabBCenter(); else slabBCenter_ = hmat(2,2) / 2.0; @@ -577,18 +578,16 @@ namespace OpenMD { } } } + // object evaluator: evaluator_.loadScriptString(rnemdObjectSelection_); seleMan_.setSelectionSet(evaluator_.evaluate()); - evaluatorA_.loadScriptString(selectionA_); evaluatorB_.loadScriptString(selectionB_); - seleManA_.setSelectionSet(evaluatorA_.evaluate()); seleManB_.setSelectionSet(evaluatorB_.evaluate()); - commonA_ = seleManA_ & seleMan_; - commonB_ = seleManB_ & seleMan_; + commonB_ = seleManB_ & seleMan_; } @@ -605,6 +604,10 @@ namespace OpenMD { #ifdef IS_MPI } #endif + + // delete all of the objects we created: + delete areaAccumulator_; + data_.clear(); } void RNEMD::doSwap(SelectionManager& smanA, SelectionManager& smanB) { @@ -1148,7 +1151,7 @@ namespace OpenMD { //if w is in the right range, so should be x, y, z. vector::iterator sdi; Vector3d vel; - for (sdi = coldBin.begin(); sdi != coldBin.end(); sdi++) { + for (sdi = coldBin.begin(); sdi != coldBin.end(); ++sdi) { if (rnemdFluxType_ == rnemdFullKE) { vel = (*sdi)->getVel() * c; (*sdi)->setVel(vel); @@ -1159,7 +1162,7 @@ namespace OpenMD { } } w = sqrt(w); - for (sdi = hotBin.begin(); sdi != hotBin.end(); sdi++) { + for (sdi = hotBin.begin(); sdi != hotBin.end(); ++sdi) { if (rnemdFluxType_ == rnemdFullKE) { vel = (*sdi)->getVel(); vel.x() *= x; @@ -1278,7 +1281,7 @@ namespace OpenMD { vector::iterator ri; RealType r1, r2, alpha0; vector > rps; - for (ri = realRoots.begin(); ri !=realRoots.end(); ri++) { + 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 ) { @@ -1310,7 +1313,7 @@ namespace OpenMD { RealType diff; pair bestPair = make_pair(1.0, 1.0); vector >::iterator rpi; - for (rpi = rps.begin(); rpi != rps.end(); rpi++) { + for (rpi = rps.begin(); rpi != rps.end(); ++rpi) { r1 = (*rpi).first; r2 = (*rpi).second; switch(rnemdFluxType_) { @@ -1377,7 +1380,7 @@ namespace OpenMD { } vector::iterator sdi; Vector3d vel; - for (sdi = coldBin.begin(); sdi != coldBin.end(); sdi++) { + for (sdi = coldBin.begin(); sdi != coldBin.end(); ++sdi) { vel = (*sdi)->getVel(); vel.x() *= x; vel.y() *= y; @@ -1388,7 +1391,7 @@ namespace OpenMD { x = 1.0 + px * (1.0 - x); y = 1.0 + py * (1.0 - y); z = 1.0 + pz * (1.0 - z); - for (sdi = hotBin.begin(); sdi != hotBin.end(); sdi++) { + for (sdi = hotBin.begin(); sdi != hotBin.end(); ++sdi) { vel = (*sdi)->getVel(); vel.x() *= x; vel.y() *= y; @@ -1666,7 +1669,7 @@ namespace OpenMD { Vector3d vel; Vector3d rPos; - for (sdi = coldBin.begin(); sdi != coldBin.end(); sdi++) { + for (sdi = coldBin.begin(); sdi != coldBin.end(); ++sdi) { //vel = (*sdi)->getVel(); rPos = (*sdi)->getPos() - coordinateOrigin_; if (doLinearPart) @@ -1682,7 +1685,7 @@ namespace OpenMD { } } } - for (sdi = hotBin.begin(); sdi != hotBin.end(); sdi++) { + for (sdi = hotBin.begin(); sdi != hotBin.end(); ++sdi) { //vel = (*sdi)->getVel(); rPos = (*sdi)->getPos() - coordinateOrigin_; if (doLinearPart) @@ -1732,14 +1735,16 @@ namespace OpenMD { int isd; StuntDouble* sd; vector aSites; - ConvexHull* surfaceMeshA = new ConvexHull(); seleManA_.setSelectionSet(evaluatorA_.evaluate()); for (sd = seleManA_.beginSelected(isd); sd != NULL; sd = seleManA_.nextSelected(isd)) { aSites.push_back(sd); } + ConvexHull* surfaceMeshA = new ConvexHull(); surfaceMeshA->computeHull(aSites); areaA = surfaceMeshA->getArea(); + delete surfaceMeshA; + } else { if (usePeriodicBoundaryConditions_) { // in periodic boundaries, the surface area is twice the x-y @@ -1753,18 +1758,22 @@ namespace OpenMD { } } + + if (hasSelectionB_) { int isd; StuntDouble* sd; vector bSites; - ConvexHull* surfaceMeshB = new ConvexHull(); seleManB_.setSelectionSet(evaluatorB_.evaluate()); for (sd = seleManB_.beginSelected(isd); sd != NULL; sd = seleManB_.nextSelected(isd)) { bSites.push_back(sd); } + ConvexHull* surfaceMeshB = new ConvexHull(); surfaceMeshB->computeHull(bSites); areaB = surfaceMeshB->getArea(); + delete surfaceMeshB; + } else { if (usePeriodicBoundaryConditions_) { // in periodic boundaries, the surface area is twice the x-y @@ -2004,6 +2013,7 @@ namespace OpenMD { } } } + hasData_ = true; } void RNEMD::getStarted() { @@ -2036,6 +2046,7 @@ namespace OpenMD { void RNEMD::writeOutputFile() { if (!doRNEMD_) return; + if (!hasData_) return; #ifdef IS_MPI // If we're the root node, should we print out the results @@ -2057,10 +2068,16 @@ namespace OpenMD { RealType time = currentSnap_->getTime(); RealType avgArea; areaAccumulator_->getAverage(avgArea); - RealType Jz = kineticExchange_ / (time * avgArea) - / PhysicalConstants::energyConvert; - Vector3d JzP = momentumExchange_ / (time * avgArea); - Vector3d JzL = angularMomentumExchange_ / (time * avgArea); + + RealType Jz(0.0); + Vector3d JzP(V3Zero); + Vector3d JzL(V3Zero); + if (time >= info_->getSimParams()->getDt()) { + Jz = kineticExchange_ / (time * avgArea) + / PhysicalConstants::energyConvert; + JzP = momentumExchange_ / (time * avgArea); + JzL = angularMomentumExchange_ / (time * avgArea); + } rnemdFile_ << "#######################################################\n"; rnemdFile_ << "# RNEMD {\n"; @@ -2214,7 +2231,7 @@ namespace OpenMD { rnemdFile_ << "\t" << s; } else{ sprintf( painCave.errMsg, - "RNEMD detected a numerical error writing: %s for bin %d", + "RNEMD detected a numerical error writing: %s for bin %u", data_[index].title.c_str(), bin); painCave.isFatal = 1; simError(); @@ -2237,7 +2254,7 @@ namespace OpenMD { isinf(s[1]) || isnan(s[1]) || isinf(s[2]) || isnan(s[2]) ) { sprintf( painCave.errMsg, - "RNEMD detected a numerical error writing: %s for bin %d", + "RNEMD detected a numerical error writing: %s for bin %u", data_[index].title.c_str(), bin); painCave.isFatal = 1; simError(); @@ -2262,7 +2279,7 @@ namespace OpenMD { rnemdFile_ << "\t" << s; } else{ sprintf( painCave.errMsg, - "RNEMD detected a numerical error writing: %s std. dev. for bin %d", + "RNEMD detected a numerical error writing: %s std. dev. for bin %u", data_[index].title.c_str(), bin); painCave.isFatal = 1; simError(); @@ -2284,7 +2301,7 @@ namespace OpenMD { isinf(s[1]) || isnan(s[1]) || isinf(s[2]) || isnan(s[2]) ) { sprintf( painCave.errMsg, - "RNEMD detected a numerical error writing: %s std. dev. for bin %d", + "RNEMD detected a numerical error writing: %s std. dev. for bin %u", data_[index].title.c_str(), bin); painCave.isFatal = 1; simError();