--- trunk/src/rnemd/RNEMD.cpp 2013/06/16 15:15:42 1879 +++ trunk/src/rnemd/RNEMD.cpp 2013/11/01 19:31:41 1940 @@ -38,6 +38,9 @@ * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008). * [4] Vardeman & Gezelter, in progress (2009). */ +#ifdef IS_MPI +#include +#endif #include #include @@ -54,9 +57,6 @@ #include "utils/Tuple.hpp" #include "brains/Thermo.hpp" #include "math/ConvexHull.hpp" -#ifdef IS_MPI -#include -#endif #ifdef _MSC_VER #define isnan(x) _isnan((x)) @@ -424,10 +424,10 @@ namespace OpenMD { OutputData angularVelocity; angularVelocity.units = "angstroms^2/fs"; angularVelocity.title = "AngularVelocity"; - angularVelocity.dataType = "Vector3d"; + angularVelocity.dataType = "RealType"; angularVelocity.accumulator.reserve(nBins_); for (int i = 0; i < nBins_; i++) - angularVelocity.accumulator.push_back( new VectorAccumulator() ); + angularVelocity.accumulator.push_back( new Accumulator() ); data_[ANGULARVELOCITY] = angularVelocity; outputMap_["ANGULARVELOCITY"] = ANGULARVELOCITY; @@ -547,6 +547,7 @@ namespace OpenMD { if (hasSelectionB_) { selectionB_ = rnemdParams->getSelectionB(); + } else { if (usePeriodicBoundaryConditions_) { Mat3x3d hmat = currentSnap_->getHmat(); @@ -573,6 +574,7 @@ namespace OpenMD { selectionB_ = selectionBstream.str(); } else { selectionB_ = "select hull"; + BisHull_ = true; hasSelectionB_ = true; } } @@ -1729,27 +1731,34 @@ namespace OpenMD { Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot(); if (hasSelectionA_) { - int isd; - StuntDouble* sd; - vector aSites; - 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(); + + if (evaluatorA_.hasSurfaceArea()) + areaA = evaluatorA_.getSurfaceArea(); + else { + + cerr << "selection A did not have surface area, recomputing\n"; + int isd; + StuntDouble* sd; + vector aSites; + 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_) { @@ -1765,30 +1774,35 @@ namespace OpenMD { } if (hasSelectionB_) { - int isd; - StuntDouble* sd; - vector bSites; - seleManB_.setSelectionSet(evaluatorB_.evaluate()); - for (sd = seleManB_.beginSelected(isd); sd != NULL; - sd = seleManB_.nextSelected(isd)) { - bSites.push_back(sd); - } + if (evaluatorB_.hasSurfaceArea()) + areaB = evaluatorB_.getSurfaceArea(); + else { + cerr << "selection B did not have surface area, recomputing\n"; -#if defined(HAVE_QHULL) - ConvexHull* surfaceMeshB = new ConvexHull(); - surfaceMeshB->computeHull(bSites); - areaB = surfaceMeshB->getArea(); - delete surfaceMeshB; + int isd; + StuntDouble* sd; + vector bSites; + 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(); + 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 @@ -1800,7 +1814,7 @@ namespace OpenMD { areaB = 4.0 * M_PI * pow(sphereBRadius_, 2); } } - + dividingArea_ = min(areaA, areaB); hasDividingArea_ = true; return dividingArea_; @@ -1860,6 +1874,9 @@ namespace OpenMD { RealType area = getDividingArea(); areaAccumulator_->add(area); Mat3x3d hmat = currentSnap_->getHmat(); + Vector3d u = angularMomentumFluxVector_; + u.normalize(); + seleMan_.setSelectionSet(evaluator_.evaluate()); int selei(0); @@ -1870,9 +1887,7 @@ namespace OpenMD { vector binPx(nBins_, 0.0); vector binPy(nBins_, 0.0); vector binPz(nBins_, 0.0); - vector binOmegax(nBins_, 0.0); - vector binOmegay(nBins_, 0.0); - vector binOmegaz(nBins_, 0.0); + vector binOmega(nBins_, 0.0); vector binKE(nBins_, 0.0); vector binDOF(nBins_, 0); vector binCount(nBins_, 0); @@ -1912,10 +1927,19 @@ namespace OpenMD { binNo = int(rPos.length() / binWidth_); } + RealType mass = sd->getMass(); Vector3d vel = sd->getVel(); Vector3d rPos = sd->getPos() - coordinateOrigin_; - Vector3d aVel = cross(rPos, vel); + // Project the relative position onto a plane perpendicular to + // the angularMomentumFluxVector: + Vector3d rProj = rPos - dot(rPos, u) * u; + // Project the velocity onto a plane perpendicular to the + // angularMomentumFluxVector: + Vector3d vProj = vel - dot(vel, u) * u; + // Compute angular velocity vector (should be nearly parallel to + // angularMomentumFluxVector + Vector3d aVel = cross(rProj, vProj); if (binNo >= 0 && binNo < nBins_) { binCount[binNo]++; @@ -1923,9 +1947,7 @@ namespace OpenMD { binPx[binNo] += mass*vel.x(); binPy[binNo] += mass*vel.y(); binPz[binNo] += mass*vel.z(); - binOmegax[binNo] += aVel.x(); - binOmegay[binNo] += aVel.y(); - binOmegaz[binNo] += aVel.z(); + binOmega[binNo] += dot(aVel, u); binKE[binNo] += 0.5 * (mass * vel.lengthSquare()); binDOF[binNo] += 3; @@ -1960,11 +1982,7 @@ namespace OpenMD { nBins_, MPI::REALTYPE, MPI::SUM); MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binPz[0], nBins_, MPI::REALTYPE, MPI::SUM); - MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binOmegax[0], - nBins_, MPI::REALTYPE, MPI::SUM); - MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binOmegay[0], - nBins_, MPI::REALTYPE, MPI::SUM); - MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binOmegaz[0], + MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binOmega[0], nBins_, MPI::REALTYPE, MPI::SUM); MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binKE[0], nBins_, MPI::REALTYPE, MPI::SUM); @@ -1973,7 +1991,7 @@ namespace OpenMD { #endif Vector3d vel; - Vector3d aVel; + RealType omega; RealType den; RealType temp; RealType z; @@ -1993,9 +2011,7 @@ namespace OpenMD { vel.x() = binPx[i] / binMass[i]; vel.y() = binPy[i] / binMass[i]; vel.z() = binPz[i] / binMass[i]; - aVel.x() = binOmegax[i] / binCount[i]; - aVel.y() = binOmegay[i] / binCount[i]; - aVel.z() = binOmegaz[i] / binCount[i]; + omega = binOmega[i] / binCount[i]; if (binCount[i] > 0) { // only add values if there are things to add @@ -2018,7 +2034,7 @@ namespace OpenMD { dynamic_cast(data_[j].accumulator[i])->add(vel); break; case ANGULARVELOCITY: - dynamic_cast(data_[j].accumulator[i])->add(aVel); + dynamic_cast(data_[j].accumulator[i])->add(omega); break; case DENSITY: dynamic_cast(data_[j].accumulator[i])->add(den);