--- branches/development/src/rnemd/RNEMD.cpp 2013/04/02 21:30:34 1856 +++ branches/development/src/rnemd/RNEMD.cpp 2013/04/30 15:56:54 1868 @@ -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), 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) { @@ -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,72 @@ namespace OpenMD { MPI::REALTYPE, MPI::SUM); #endif + + 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; - 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)); + 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)); + 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));; + 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)); + 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); @@ -1606,8 +1672,11 @@ namespace OpenMD { 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()) { @@ -1619,8 +1688,11 @@ namespace OpenMD { 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); + 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()) { @@ -1663,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 @@ -1684,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 @@ -1900,9 +1978,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 @@ -1935,6 +2013,7 @@ namespace OpenMD { } } } + hasData_ = true; } void RNEMD::getStarted() { @@ -1967,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 @@ -1988,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";