--- trunk/src/rnemd/RNEMD.cpp 2013/06/16 15:15:42 1879 +++ trunk/src/rnemd/RNEMD.cpp 2013/11/12 02:18:35 1946 @@ -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)) @@ -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"; + 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; + 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,19 +1874,27 @@ namespace OpenMD { RealType area = getDividingArea(); areaAccumulator_->add(area); Mat3x3d hmat = currentSnap_->getHmat(); + Vector3d u = angularMomentumFluxVector_; + u.normalize(); + seleMan_.setSelectionSet(evaluator_.evaluate()); int selei(0); StuntDouble* sd; int binNo; + RealType mass; + Vector3d vel; + Vector3d rPos; + RealType KE; + Vector3d L; + Mat3x3d I; + RealType r2; vector binMass(nBins_, 0.0); - 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 binP(nBins_, V3Zero); + vector binOmega(nBins_, 0.0); + vector binL(nBins_, V3Zero); + vector binI(nBins_); vector binKE(nBins_, 0.0); vector binDOF(nBins_, 0); vector binCount(nBins_, 0); @@ -1912,37 +1934,50 @@ 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); - + mass = sd->getMass(); + vel = sd->getVel(); + rPos = sd->getPos() - coordinateOrigin_; + KE = 0.5 * mass * vel.lengthSquare(); + L = mass * cross(rPos, vel); + I = outProduct(rPos, rPos) * mass; + r2 = rPos.lengthSquare(); + I(0, 0) += mass * r2; + I(1, 1) += mass * r2; + I(2, 2) += mass * r2; + + // 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]++; binMass[binNo] += mass; - 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(); - binKE[binNo] += 0.5 * (mass * vel.lengthSquare()); + binP[binNo] += mass*vel; + binKE[binNo] += KE; + binI[binNo] += I; + binL[binNo] += L; binDOF[binNo] += 3; if (sd->isDirectional()) { Vector3d angMom = sd->getJ(); - Mat3x3d I = sd->getI(); + Mat3x3d Ia = sd->getI(); if (sd->isLinear()) { int i = sd->linearAxis(); int j = (i + 1) % 3; int k = (i + 2) % 3; - binKE[binNo] += 0.5 * (angMom[j] * angMom[j] / I(j, j) + - angMom[k] * angMom[k] / I(k, k)); + binKE[binNo] += 0.5 * (angMom[j] * angMom[j] / Ia(j, j) + + angMom[k] * angMom[k] / Ia(k, k)); binDOF[binNo] += 2; } else { - binKE[binNo] += 0.5 * (angMom[0] * angMom[0] / I(0, 0) + - angMom[1] * angMom[1] / I(1, 1) + - angMom[2] * angMom[2] / I(2, 2)); + binKE[binNo] += 0.5 * (angMom[0] * angMom[0] / Ia(0, 0) + + angMom[1] * angMom[1] / Ia(1, 1) + + angMom[2] * angMom[2] / Ia(2, 2)); binDOF[binNo] += 3; } } @@ -1950,30 +1985,30 @@ namespace OpenMD { } #ifdef IS_MPI - MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binCount[0], - nBins_, MPI::INT, MPI::SUM); - MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binMass[0], - nBins_, MPI::REALTYPE, MPI::SUM); - MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binPx[0], - nBins_, MPI::REALTYPE, MPI::SUM); - MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binPy[0], - 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], - nBins_, MPI::REALTYPE, MPI::SUM); - MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binKE[0], - nBins_, MPI::REALTYPE, MPI::SUM); - MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binDOF[0], - nBins_, MPI::INT, MPI::SUM); + + for (int i = 0; i < nBins_; i++) { + + MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binCount[i], + 1, MPI::INT, MPI::SUM); + MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binMass[i], + 1, MPI::REALTYPE, MPI::SUM); + MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binP[i], + 3, MPI::REALTYPE, MPI::SUM); + MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binL[i], + 3, MPI::REALTYPE, MPI::SUM); + MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binI[i], + 9, MPI::REALTYPE, MPI::SUM); + MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binKE[i], + 1, MPI::REALTYPE, MPI::SUM); + MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binDOF[i], + 1, MPI::INT, MPI::SUM); + //MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &binOmega[i], + // 1, MPI::REALTYPE, MPI::SUM); + } + #endif - Vector3d vel; - Vector3d aVel; + Vector3d omega; RealType den; RealType temp; RealType z; @@ -1990,13 +2025,12 @@ namespace OpenMD { den = binMass[i] * 3.0 * PhysicalConstants::densityConvert / (4.0 * M_PI * (pow(router,3) - pow(rinner,3))); } - 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]; + vel = binP[i] / binMass[i]; + omega = binI[i].inverse() * binL[i]; + + // omega = binOmega[i] / binCount[i]; + if (binCount[i] > 0) { // only add values if there are things to add temp = 2.0 * binKE[i] / (binDOF[i] * PhysicalConstants::kb * @@ -2018,7 +2052,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);