--- branches/development/src/rnemd/RNEMD.cpp 2013/04/29 17:53:48 1867 +++ branches/development/src/rnemd/RNEMD.cpp 2013/05/17 17:10:11 1876 @@ -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; @@ -603,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) { @@ -1146,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); @@ -1157,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; @@ -1276,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 ) { @@ -1308,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_) { @@ -1375,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; @@ -1386,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; @@ -1592,10 +1597,7 @@ namespace OpenMD { Vector3d ac, acrec, bc, bcrec; Vector3d ah, ahrec, bh, bhrec; - RealType cNumerator, cDenominator; - RealType hNumerator, hDenominator; - bool successfulExchange = false; if ((Mh > 0.0) && (Mc > 0.0)) {//both slabs are not empty Vector3d vc = Pc / Mc; @@ -1609,7 +1611,7 @@ namespace OpenMD { bc = -(Ici * angularMomentumTarget_) + omegac; bcrec = bc - omegac; - cNumerator = Kc - kineticTarget_; + RealType cNumerator = Kc - kineticTarget_; if (doLinearPart) cNumerator -= 0.5 * Mc * ac.lengthSquare(); @@ -1618,7 +1620,7 @@ namespace OpenMD { if (cNumerator > 0.0) { - cDenominator = Kc; + RealType cDenominator = Kc; if (doLinearPart) cDenominator -= 0.5 * Mc * vc.lengthSquare(); @@ -1641,7 +1643,7 @@ namespace OpenMD { bh = (Ihi * angularMomentumTarget_) + omegah; bhrec = bh - omegah; - hNumerator = Kh + kineticTarget_; + RealType hNumerator = Kh + kineticTarget_; if (doLinearPart) hNumerator -= 0.5 * Mh * ah.lengthSquare(); @@ -1650,7 +1652,7 @@ namespace OpenMD { if (hNumerator > 0.0) { - hDenominator = Kh; + RealType hDenominator = Kh; if (doLinearPart) hDenominator -= 0.5 * Mh * vh.lengthSquare(); if (doAngularPart) @@ -1664,7 +1666,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) @@ -1680,7 +1682,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) @@ -1735,10 +1737,19 @@ namespace OpenMD { sd = seleManA_.nextSelected(isd)) { aSites.push_back(sd); } +#if defined(HAVE_QHULL) ConvexHull* surfaceMeshA = new ConvexHull(); surfaceMeshA->computeHull(aSites); areaA = surfaceMeshA->getArea(); delete surfaceMeshA; +#else + sprintf( painCave.errMsg, + "RNEMD::getDividingArea : Hull calculation is not possible\n" + "\twithout libqhull. Please rebuild OpenMD with qhull enabled."); + painCave.severity = OPENMD_ERROR; + painCave.isFatal = 1; + simError(); +#endif } else { if (usePeriodicBoundaryConditions_) { @@ -1753,8 +1764,6 @@ namespace OpenMD { } } - - if (hasSelectionB_) { int isd; StuntDouble* sd; @@ -1764,11 +1773,22 @@ namespace OpenMD { sd = seleManB_.nextSelected(isd)) { bSites.push_back(sd); } + +#if defined(HAVE_QHULL) ConvexHull* surfaceMeshB = new ConvexHull(); surfaceMeshB->computeHull(bSites); areaB = surfaceMeshB->getArea(); delete surfaceMeshB; +#else + sprintf( painCave.errMsg, + "RNEMD::getDividingArea : Hull calculation is not possible\n" + "\twithout libqhull. Please rebuild OpenMD with qhull enabled."); + painCave.severity = OPENMD_ERROR; + painCave.isFatal = 1; + simError(); +#endif + } else { if (usePeriodicBoundaryConditions_) { // in periodic boundaries, the surface area is twice the x-y @@ -2008,6 +2028,7 @@ namespace OpenMD { } } } + hasData_ = true; } void RNEMD::getStarted() { @@ -2040,6 +2061,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 @@ -2061,11 +2083,17 @@ 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"; @@ -2218,7 +2246,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(); @@ -2241,7 +2269,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(); @@ -2266,7 +2294,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(); @@ -2288,7 +2316,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();