--- branches/development/src/rnemd/RNEMD.cpp 2013/04/02 18:31:51 1855 +++ 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; @@ -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; @@ -1445,6 +1448,44 @@ namespace OpenMD { RealType Mc = 0.0; Mat3x3d Ic(0.0); RealType Kc = 0.0; + + // Constraints can be on only the linear or angular momentum, but + // not both. Usually, the user will specify which they want, but + // in case they don't, the use of periodic boundaries should make + // the choice for us. + bool doLinearPart = false; + bool doAngularPart = false; + + switch (rnemdFluxType_) { + case rnemdPx: + case rnemdPy: + case rnemdPz: + case rnemdPvector: + case rnemdKePx: + case rnemdKePy: + case rnemdKePvector: + doLinearPart = true; + break; + case rnemdLx: + case rnemdLy: + case rnemdLz: + case rnemdLvector: + case rnemdKeLx: + case rnemdKeLy: + case rnemdKeLz: + case rnemdKeLvector: + doAngularPart = true; + break; + case rnemdKE: + case rnemdRotKE: + case rnemdFullKE: + default: + if (usePeriodicBoundaryConditions_) + doLinearPart = true; + else + doAngularPart = true; + break; + } for (sd = smanA.beginSelected(selei); sd != NULL; sd = smanA.nextSelected(selei)) { @@ -1553,47 +1594,69 @@ namespace OpenMD { MPI::REALTYPE, MPI::SUM); #endif + + Vector3d ac, acrec, bc, bcrec; + Vector3d ah, ahrec, bh, bhrec; + bool successfulExchange = false; if ((Mh > 0.0) && (Mc > 0.0)) {//both slabs are not empty Vector3d vc = Pc / Mc; - Vector3d ac = -momentumTarget_ / Mc + vc; - Vector3d acrec = -momentumTarget_ / Mc; + ac = -momentumTarget_ / Mc + vc; + acrec = -momentumTarget_ / Mc; // We now need the inverse of the inertia tensor to calculate the // angular velocity of the cold slab; Mat3x3d Ici = Ic.inverse(); Vector3d omegac = Ici * Lc; - Vector3d bc = -(Ici * angularMomentumTarget_) + omegac; - Vector3d bcrec = bc - omegac; + bc = -(Ici * angularMomentumTarget_) + omegac; + bcrec = bc - omegac; - RealType cNumerator = Kc - kineticTarget_ - - 0.5 * Mc * ac.lengthSquare() - 0.5 * ( dot(bc, Ic * bc)); + RealType cNumerator = Kc - kineticTarget_; + if (doLinearPart) + cNumerator -= 0.5 * Mc * ac.lengthSquare(); + + if (doAngularPart) + cNumerator -= 0.5 * ( dot(bc, Ic * bc)); + if (cNumerator > 0.0) { - RealType cDenominator = Kc - 0.5 * Mc * vc.lengthSquare() - - 0.5*(dot(omegac, Ic * omegac)); + RealType cDenominator = Kc; + + if (doLinearPart) + cDenominator -= 0.5 * Mc * vc.lengthSquare(); + + if (doAngularPart) + cDenominator -= 0.5*(dot(omegac, Ic * omegac)); if (cDenominator > 0.0) { RealType c = sqrt(cNumerator / cDenominator); if ((c > 0.9) && (c < 1.1)) {//restrict scaling coefficients Vector3d vh = Ph / Mh; - Vector3d ah = momentumTarget_ / Mh + vh; - Vector3d ahrec = momentumTarget_ / Mh; + ah = momentumTarget_ / Mh + vh; + ahrec = momentumTarget_ / Mh; // We now need the inverse of the inertia tensor to // calculate the angular velocity of the hot slab; Mat3x3d Ihi = Ih.inverse(); Vector3d omegah = Ihi * Lh; - Vector3d bh = (Ihi * angularMomentumTarget_) + omegah; - Vector3d bhrec = bh - omegah; + bh = (Ihi * angularMomentumTarget_) + omegah; + bhrec = bh - omegah; - RealType hNumerator = Kh + kineticTarget_ - - 0.5 * Mh * ah.lengthSquare() - 0.5 * ( dot(bh, Ih * bh));; + RealType hNumerator = Kh + kineticTarget_; + if (doLinearPart) + hNumerator -= 0.5 * Mh * ah.lengthSquare(); + + if (doAngularPart) + hNumerator -= 0.5 * ( dot(bh, Ih * bh)); + if (hNumerator > 0.0) { - RealType hDenominator = Kh - 0.5 * Mh * vh.lengthSquare() - - 0.5*(dot(omegah, Ih * omegah)); + RealType hDenominator = Kh; + if (doLinearPart) + hDenominator -= 0.5 * Mh * vh.lengthSquare(); + if (doAngularPart) + hDenominator -= 0.5*(dot(omegah, Ih * omegah)); if (hDenominator > 0.0) { RealType h = sqrt(hNumerator / hDenominator); @@ -1603,11 +1666,14 @@ 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_; - vel = ((*sdi)->getVel() - vc - cross(omegac, rPos)) * c - + ac + cross(bc, rPos); + if (doLinearPart) + vel = ((*sdi)->getVel() - vc) * c + ac; + if (doAngularPart) + vel = ((*sdi)->getVel() - cross(omegac, rPos)) * c + cross(bc, rPos); + (*sdi)->setVel(vel); if (rnemdFluxType_ == rnemdFullKE) { if ((*sdi)->isDirectional()) { @@ -1616,12 +1682,14 @@ 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_; - vel = ((*sdi)->getVel() - vh - cross(omegah, rPos)) * h - + ah + cross(bh, rPos); - cerr << "setting vel to " << vel << "\n"; + if (doLinearPart) + vel = ((*sdi)->getVel() - vh) * h + ah; + if (doAngularPart) + vel = ((*sdi)->getVel() - cross(omegah, rPos)) * h + cross(bh, rPos); + (*sdi)->setVel(vel); if (rnemdFluxType_ == rnemdFullKE) { if ((*sdi)->isDirectional()) { @@ -1664,14 +1732,25 @@ 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); } +#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_) { // in periodic boundaries, the surface area is twice the x-y @@ -1689,14 +1768,27 @@ namespace OpenMD { 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); } + +#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 @@ -1718,7 +1810,6 @@ namespace OpenMD { if (!doRNEMD_) return; trialCount_++; - cerr << "trialCount = " << trialCount_ << "\n"; // object evaluator: evaluator_.loadScriptString(rnemdObjectSelection_); seleMan_.setSelectionSet(evaluator_.evaluate()); @@ -1764,7 +1855,6 @@ namespace OpenMD { if (!doRNEMD_) return; Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot(); - cerr << "collecting data\n"; // collectData can be called more frequently than the doRNEMD, so use the // computed area from the last exchange time: RealType area = getDividingArea(); @@ -1903,9 +1993,9 @@ namespace OpenMD { vel.x() = binPx[i] / binMass[i]; vel.y() = binPy[i] / binMass[i]; vel.z() = binPz[i] / binMass[i]; - aVel.x() = binOmegax[i]; - aVel.y() = binOmegay[i]; - aVel.z() = binOmegaz[i]; + aVel.x() = binOmegax[i] / binCount[i]; + aVel.y() = binOmegay[i] / binCount[i]; + aVel.z() = binOmegaz[i] / binCount[i]; if (binCount[i] > 0) { // only add values if there are things to add @@ -1938,6 +2028,7 @@ namespace OpenMD { } } } + hasData_ = true; } void RNEMD::getStarted() { @@ -1970,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 @@ -1991,10 +2083,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"; @@ -2148,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(); @@ -2171,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(); @@ -2196,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(); @@ -2218,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();