--- branches/development/src/rnemd/RNEMD.cpp 2013/04/02 18:31:51 1855 +++ branches/development/src/rnemd/RNEMD.cpp 2013/04/29 17:53:48 1867 @@ -556,7 +556,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 +577,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_; } @@ -1445,6 +1443,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 +1589,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 +1667,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,9 +1683,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); - 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 +1730,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 @@ -1685,18 +1753,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 @@ -1718,7 +1790,6 @@ namespace OpenMD { if (!doRNEMD_) return; trialCount_++; - cerr << "trialCount = " << trialCount_ << "\n"; // object evaluator: evaluator_.loadScriptString(rnemdObjectSelection_); seleMan_.setSelectionSet(evaluator_.evaluate()); @@ -1764,7 +1835,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 +1973,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