--- branches/development/src/rnemd/RNEMD.cpp 2013/03/28 20:54:06 1854 +++ branches/development/src/rnemd/RNEMD.cpp 2013/04/09 19:45:54 1861 @@ -419,6 +419,16 @@ namespace OpenMD { velocity.accumulator.push_back( new VectorAccumulator() ); data_[VELOCITY] = velocity; outputMap_["VELOCITY"] = VELOCITY; + + OutputData angularVelocity; + angularVelocity.units = "angstroms^2/fs"; + angularVelocity.title = "AngularVelocity"; + angularVelocity.dataType = "Vector3d"; + angularVelocity.accumulator.reserve(nBins_); + for (int i = 0; i < nBins_; i++) + angularVelocity.accumulator.push_back( new VectorAccumulator() ); + data_[ANGULARVELOCITY] = angularVelocity; + outputMap_["ANGULARVELOCITY"] = ANGULARVELOCITY; OutputData density; density.units = "g cm^-3"; @@ -1435,6 +1445,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)) { @@ -1543,47 +1591,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); @@ -1596,8 +1669,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()) { @@ -1609,8 +1685,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()) { @@ -1751,11 +1830,11 @@ namespace OpenMD { void RNEMD::collectData() { if (!doRNEMD_) return; Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot(); - + // collectData can be called more frequently than the doRNEMD, so use the // computed area from the last exchange time: - - areaAccumulator_->add(getDividingArea()); + RealType area = getDividingArea(); + areaAccumulator_->add(area); Mat3x3d hmat = currentSnap_->getHmat(); seleMan_.setSelectionSet(evaluator_.evaluate()); @@ -1814,7 +1893,7 @@ namespace OpenMD { Vector3d rPos = sd->getPos() - coordinateOrigin_; Vector3d aVel = cross(rPos, vel); - if (binNo < nBins_) { + if (binNo >= 0 && binNo < nBins_) { binCount[binNo]++; binMass[binNo] += mass; binPx[binNo] += mass*vel.x(); @@ -1890,9 +1969,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 @@ -1914,7 +1993,7 @@ namespace OpenMD { case VELOCITY: dynamic_cast(data_[j].accumulator[i])->add(vel); break; - case ANGULARVELOCITY: + case ANGULARVELOCITY: dynamic_cast(data_[j].accumulator[i])->add(aVel); break; case DENSITY: @@ -2069,7 +2148,7 @@ namespace OpenMD { if (outputMask_[i]) { if (data_[i].dataType == "RealType") writeReal(i,j); - else if (data_[i].dataType == "Vector3d") + else if (data_[i].dataType == "Vector3d") writeVector(i,j); else { sprintf( painCave.errMsg, @@ -2126,7 +2205,7 @@ namespace OpenMD { RealType s; int count; - count = dynamic_cast(data_[index].accumulator[bin])->count(); + count = data_[index].accumulator[bin]->count(); if (count == 0) return; dynamic_cast(data_[index].accumulator[bin])->getAverage(s); @@ -2149,7 +2228,8 @@ namespace OpenMD { Vector3d s; int count; - count = dynamic_cast(data_[index].accumulator[bin])->count(); + count = data_[index].accumulator[bin]->count(); + if (count == 0) return; dynamic_cast(data_[index].accumulator[bin])->getAverage(s); @@ -2173,7 +2253,7 @@ namespace OpenMD { RealType s; int count; - count = dynamic_cast(data_[index].accumulator[bin])->count(); + count = data_[index].accumulator[bin]->count(); if (count == 0) return; dynamic_cast(data_[index].accumulator[bin])->getStdDev(s); @@ -2196,7 +2276,7 @@ namespace OpenMD { Vector3d s; int count; - count = dynamic_cast(data_[index].accumulator[bin])->count(); + count = data_[index].accumulator[bin]->count(); if (count == 0) return; dynamic_cast(data_[index].accumulator[bin])->getStdDev(s);