--- branches/development/src/rnemd/RNEMD.cpp 2013/03/28 20:54:06 1854 +++ 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; @@ -419,6 +420,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"; @@ -546,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; @@ -567,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_; } @@ -595,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) { @@ -1435,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)) { @@ -1543,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); @@ -1596,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()) { @@ -1609,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()) { @@ -1653,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 @@ -1674,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 @@ -1751,11 +1839,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 +1902,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 +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 @@ -1914,7 +2002,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: @@ -1925,6 +2013,7 @@ namespace OpenMD { } } } + hasData_ = true; } void RNEMD::getStarted() { @@ -1957,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 @@ -1978,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"; @@ -2069,7 +2165,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 +2222,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 +2245,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 +2270,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 +2293,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);