--- trunk/src/rnemd/RNEMD.cpp 2013/06/16 15:15:42 1879 +++ trunk/src/rnemd/RNEMD.cpp 2015/03/05 15:40:58 2068 @@ -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)) @@ -291,7 +291,18 @@ namespace OpenMD { kineticFlux_ = 0.0; } if (hasMomentumFluxVector) { - momentumFluxVector_ = rnemdParams->getMomentumFluxVector(); + std::vector mf = rnemdParams->getMomentumFluxVector(); + if (mf.size() != 3) { + sprintf(painCave.errMsg, + "RNEMD: Incorrect number of parameters specified for momentumFluxVector.\n" + "\tthere should be 3 parameters, but %lu were specified.\n", + mf.size()); + painCave.isFatal = 1; + simError(); + } + momentumFluxVector_.x() = mf[0]; + momentumFluxVector_.y() = mf[1]; + momentumFluxVector_.z() = mf[2]; } else { momentumFluxVector_ = V3Zero; if (hasMomentumFlux) { @@ -316,269 +327,292 @@ namespace OpenMD { break; } } - if (hasAngularMomentumFluxVector) { - angularMomentumFluxVector_ = rnemdParams->getAngularMomentumFluxVector(); - } else { - angularMomentumFluxVector_ = V3Zero; - if (hasAngularMomentumFlux) { - RealType angularMomentumFlux = rnemdParams->getAngularMomentumFlux(); - switch (rnemdFluxType_) { - case rnemdLx: - angularMomentumFluxVector_.x() = angularMomentumFlux; - break; - case rnemdLy: - angularMomentumFluxVector_.y() = angularMomentumFlux; - break; - case rnemdLz: - angularMomentumFluxVector_.z() = angularMomentumFlux; - break; - case rnemdKeLx: - angularMomentumFluxVector_.x() = angularMomentumFlux; - break; - case rnemdKeLy: - angularMomentumFluxVector_.y() = angularMomentumFlux; - break; - case rnemdKeLz: - angularMomentumFluxVector_.z() = angularMomentumFlux; - break; - default: - break; - } - } - } - - if (hasCoordinateOrigin) { - coordinateOrigin_ = rnemdParams->getCoordinateOrigin(); - } else { - coordinateOrigin_ = V3Zero; - } - - // do some sanity checking - - int selectionCount = seleMan_.getSelectionCount(); - - int nIntegrable = info->getNGlobalIntegrableObjects(); - - if (selectionCount > nIntegrable) { - sprintf(painCave.errMsg, - "RNEMD: The current objectSelection,\n" - "\t\t%s\n" - "\thas resulted in %d selected objects. However,\n" - "\tthe total number of integrable objects in the system\n" - "\tis only %d. This is almost certainly not what you want\n" - "\tto do. A likely cause of this is forgetting the _RB_0\n" - "\tselector in the selection script!\n", - rnemdObjectSelection_.c_str(), - selectionCount, nIntegrable); - painCave.isFatal = 0; - painCave.severity = OPENMD_WARNING; - simError(); + } + if (hasAngularMomentumFluxVector) { + std::vector amf = rnemdParams->getAngularMomentumFluxVector(); + if (amf.size() != 3) { + sprintf(painCave.errMsg, + "RNEMD: Incorrect number of parameters specified for angularMomentumFluxVector.\n" + "\tthere should be 3 parameters, but %lu were specified.\n", + amf.size()); + painCave.isFatal = 1; + simError(); } - - areaAccumulator_ = new Accumulator(); - - nBins_ = rnemdParams->getOutputBins(); - binWidth_ = rnemdParams->getOutputBinWidth(); - - data_.resize(RNEMD::ENDINDEX); - OutputData z; - z.units = "Angstroms"; - z.title = "Z"; - z.dataType = "RealType"; - z.accumulator.reserve(nBins_); - for (int i = 0; i < nBins_; i++) - z.accumulator.push_back( new Accumulator() ); - data_[Z] = z; - outputMap_["Z"] = Z; - - OutputData r; - r.units = "Angstroms"; - r.title = "R"; - r.dataType = "RealType"; - r.accumulator.reserve(nBins_); - for (int i = 0; i < nBins_; i++) - r.accumulator.push_back( new Accumulator() ); - data_[R] = r; - outputMap_["R"] = R; - - OutputData temperature; - temperature.units = "K"; - temperature.title = "Temperature"; - temperature.dataType = "RealType"; - temperature.accumulator.reserve(nBins_); - for (int i = 0; i < nBins_; i++) - temperature.accumulator.push_back( new Accumulator() ); - data_[TEMPERATURE] = temperature; - outputMap_["TEMPERATURE"] = TEMPERATURE; - - OutputData velocity; - velocity.units = "angstroms/fs"; - velocity.title = "Velocity"; - velocity.dataType = "Vector3d"; - velocity.accumulator.reserve(nBins_); - for (int i = 0; i < nBins_; i++) - 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"; - density.title = "Density"; - density.dataType = "RealType"; - density.accumulator.reserve(nBins_); - for (int i = 0; i < nBins_; i++) - density.accumulator.push_back( new Accumulator() ); - data_[DENSITY] = density; - outputMap_["DENSITY"] = DENSITY; - - if (hasOutputFields) { - parseOutputFileFormat(rnemdParams->getOutputFields()); - } else { - if (usePeriodicBoundaryConditions_) - outputMask_.set(Z); - else - outputMask_.set(R); + angularMomentumFluxVector_.x() = amf[0]; + angularMomentumFluxVector_.y() = amf[1]; + angularMomentumFluxVector_.z() = amf[2]; + } else { + angularMomentumFluxVector_ = V3Zero; + if (hasAngularMomentumFlux) { + RealType angularMomentumFlux = rnemdParams->getAngularMomentumFlux(); switch (rnemdFluxType_) { - case rnemdKE: - case rnemdRotKE: - case rnemdFullKE: - outputMask_.set(TEMPERATURE); - break; - case rnemdPx: - case rnemdPy: - outputMask_.set(VELOCITY); - break; - case rnemdPz: - case rnemdPvector: - outputMask_.set(VELOCITY); - outputMask_.set(DENSITY); - break; case rnemdLx: + angularMomentumFluxVector_.x() = angularMomentumFlux; + break; case rnemdLy: + angularMomentumFluxVector_.y() = angularMomentumFlux; + break; case rnemdLz: - case rnemdLvector: - outputMask_.set(ANGULARVELOCITY); + angularMomentumFluxVector_.z() = angularMomentumFlux; break; case rnemdKeLx: + angularMomentumFluxVector_.x() = angularMomentumFlux; + break; case rnemdKeLy: + angularMomentumFluxVector_.y() = angularMomentumFlux; + break; case rnemdKeLz: - case rnemdKeLvector: - outputMask_.set(TEMPERATURE); - outputMask_.set(ANGULARVELOCITY); + angularMomentumFluxVector_.z() = angularMomentumFlux; break; - case rnemdKePx: - case rnemdKePy: - outputMask_.set(TEMPERATURE); - outputMask_.set(VELOCITY); - break; - case rnemdKePvector: - outputMask_.set(TEMPERATURE); - outputMask_.set(VELOCITY); - outputMask_.set(DENSITY); - break; default: break; } - } - - if (hasOutputFileName) { - rnemdFileName_ = rnemdParams->getOutputFileName(); - } else { - rnemdFileName_ = getPrefix(info->getFinalConfigFileName()) + ".rnemd"; - } + } + } - exchangeTime_ = rnemdParams->getExchangeTime(); - - Snapshot* currentSnap_ = info->getSnapshotManager()->getCurrentSnapshot(); - // total exchange sums are zeroed out at the beginning: - - kineticExchange_ = 0.0; - momentumExchange_ = V3Zero; - angularMomentumExchange_ = V3Zero; - - std::ostringstream selectionAstream; - std::ostringstream selectionBstream; - - if (hasSelectionA_) { - selectionA_ = rnemdParams->getSelectionA(); - } else { - if (usePeriodicBoundaryConditions_) { - Mat3x3d hmat = currentSnap_->getHmat(); + if (hasCoordinateOrigin) { + std::vector co = rnemdParams->getCoordinateOrigin(); + if (co.size() != 3) { + sprintf(painCave.errMsg, + "RNEMD: Incorrect number of parameters specified for coordinateOrigin.\n" + "\tthere should be 3 parameters, but %lu were specified.\n", + co.size()); + painCave.isFatal = 1; + simError(); + } + coordinateOrigin_.x() = co[0]; + coordinateOrigin_.y() = co[1]; + coordinateOrigin_.z() = co[2]; + } else { + coordinateOrigin_ = V3Zero; + } + + // do some sanity checking + + int selectionCount = seleMan_.getSelectionCount(); + int nIntegrable = info->getNGlobalIntegrableObjects(); + if (selectionCount > nIntegrable) { + sprintf(painCave.errMsg, + "RNEMD: The current objectSelection,\n" + "\t\t%s\n" + "\thas resulted in %d selected objects. However,\n" + "\tthe total number of integrable objects in the system\n" + "\tis only %d. This is almost certainly not what you want\n" + "\tto do. A likely cause of this is forgetting the _RB_0\n" + "\tselector in the selection script!\n", + rnemdObjectSelection_.c_str(), + selectionCount, nIntegrable); + painCave.isFatal = 0; + painCave.severity = OPENMD_WARNING; + simError(); + } + + areaAccumulator_ = new Accumulator(); + + nBins_ = rnemdParams->getOutputBins(); + binWidth_ = rnemdParams->getOutputBinWidth(); + + data_.resize(RNEMD::ENDINDEX); + OutputData z; + z.units = "Angstroms"; + z.title = "Z"; + z.dataType = "RealType"; + z.accumulator.reserve(nBins_); + for (int i = 0; i < nBins_; i++) + z.accumulator.push_back( new Accumulator() ); + data_[Z] = z; + outputMap_["Z"] = Z; + + OutputData r; + r.units = "Angstroms"; + r.title = "R"; + r.dataType = "RealType"; + r.accumulator.reserve(nBins_); + for (int i = 0; i < nBins_; i++) + r.accumulator.push_back( new Accumulator() ); + data_[R] = r; + outputMap_["R"] = R; + + OutputData temperature; + temperature.units = "K"; + temperature.title = "Temperature"; + temperature.dataType = "RealType"; + temperature.accumulator.reserve(nBins_); + for (int i = 0; i < nBins_; i++) + temperature.accumulator.push_back( new Accumulator() ); + data_[TEMPERATURE] = temperature; + outputMap_["TEMPERATURE"] = TEMPERATURE; + + OutputData velocity; + velocity.units = "angstroms/fs"; + velocity.title = "Velocity"; + velocity.dataType = "Vector3d"; + velocity.accumulator.reserve(nBins_); + for (int i = 0; i < nBins_; i++) + 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"; + density.title = "Density"; + density.dataType = "RealType"; + density.accumulator.reserve(nBins_); + for (int i = 0; i < nBins_; i++) + density.accumulator.push_back( new Accumulator() ); + data_[DENSITY] = density; + outputMap_["DENSITY"] = DENSITY; + + if (hasOutputFields) { + parseOutputFileFormat(rnemdParams->getOutputFields()); + } else { + if (usePeriodicBoundaryConditions_) + outputMask_.set(Z); + else + outputMask_.set(R); + switch (rnemdFluxType_) { + case rnemdKE: + case rnemdRotKE: + case rnemdFullKE: + outputMask_.set(TEMPERATURE); + break; + case rnemdPx: + case rnemdPy: + outputMask_.set(VELOCITY); + break; + case rnemdPz: + case rnemdPvector: + outputMask_.set(VELOCITY); + outputMask_.set(DENSITY); + break; + case rnemdLx: + case rnemdLy: + case rnemdLz: + case rnemdLvector: + outputMask_.set(ANGULARVELOCITY); + break; + case rnemdKeLx: + case rnemdKeLy: + case rnemdKeLz: + case rnemdKeLvector: + outputMask_.set(TEMPERATURE); + outputMask_.set(ANGULARVELOCITY); + break; + case rnemdKePx: + case rnemdKePy: + outputMask_.set(TEMPERATURE); + outputMask_.set(VELOCITY); + break; + case rnemdKePvector: + outputMask_.set(TEMPERATURE); + outputMask_.set(VELOCITY); + outputMask_.set(DENSITY); + break; + default: + break; + } + } + + if (hasOutputFileName) { + rnemdFileName_ = rnemdParams->getOutputFileName(); + } else { + rnemdFileName_ = getPrefix(info->getFinalConfigFileName()) + ".rnemd"; + } + + exchangeTime_ = rnemdParams->getExchangeTime(); + + Snapshot* currentSnap_ = info->getSnapshotManager()->getCurrentSnapshot(); + // total exchange sums are zeroed out at the beginning: + + kineticExchange_ = 0.0; + momentumExchange_ = V3Zero; + angularMomentumExchange_ = V3Zero; + + std::ostringstream selectionAstream; + std::ostringstream selectionBstream; + + if (hasSelectionA_) { + selectionA_ = rnemdParams->getSelectionA(); + } else { + if (usePeriodicBoundaryConditions_) { + Mat3x3d hmat = currentSnap_->getHmat(); - if (hasSlabWidth) - slabWidth_ = rnemdParams->getSlabWidth(); - else - slabWidth_ = hmat(2,2) / 10.0; + if (hasSlabWidth) + slabWidth_ = rnemdParams->getSlabWidth(); + else + slabWidth_ = hmat(2,2) / 10.0; - if (hasSlabACenter) - slabACenter_ = rnemdParams->getSlabACenter(); - else - slabACenter_ = 0.0; + if (hasSlabACenter) + slabACenter_ = rnemdParams->getSlabACenter(); + else + slabACenter_ = 0.0; - selectionAstream << "select wrappedz > " - << slabACenter_ - 0.5*slabWidth_ - << " && wrappedz < " - << slabACenter_ + 0.5*slabWidth_; - selectionA_ = selectionAstream.str(); - } else { - if (hasSphereARadius) - sphereARadius_ = rnemdParams->getSphereARadius(); - else { - // use an initial guess to the size of the inner slab to be 1/10 the - // radius of an approximately spherical hull: - Thermo thermo(info); - RealType hVol = thermo.getHullVolume(); - sphereARadius_ = 0.1 * pow((3.0 * hVol / (4.0 * M_PI)), 1.0/3.0); - } - selectionAstream << "select r < " << sphereARadius_; - selectionA_ = selectionAstream.str(); + selectionAstream << "select wrappedz > " + << slabACenter_ - 0.5*slabWidth_ + << " && wrappedz < " + << slabACenter_ + 0.5*slabWidth_; + selectionA_ = selectionAstream.str(); + } else { + if (hasSphereARadius) + sphereARadius_ = rnemdParams->getSphereARadius(); + else { + // use an initial guess to the size of the inner slab to be 1/10 the + // radius of an approximately spherical hull: + Thermo thermo(info); + RealType hVol = thermo.getHullVolume(); + sphereARadius_ = 0.1 * pow((3.0 * hVol / (4.0 * M_PI)), 1.0/3.0); } + selectionAstream << "select r < " << sphereARadius_; + selectionA_ = selectionAstream.str(); } + } - if (hasSelectionB_) { - selectionB_ = rnemdParams->getSelectionB(); - } else { - if (usePeriodicBoundaryConditions_) { - Mat3x3d hmat = currentSnap_->getHmat(); + if (hasSelectionB_) { + selectionB_ = rnemdParams->getSelectionB(); + + } else { + if (usePeriodicBoundaryConditions_) { + Mat3x3d hmat = currentSnap_->getHmat(); - if (hasSlabWidth) - slabWidth_ = rnemdParams->getSlabWidth(); - else - slabWidth_ = hmat(2,2) / 10.0; + if (hasSlabWidth) + slabWidth_ = rnemdParams->getSlabWidth(); + else + slabWidth_ = hmat(2,2) / 10.0; - if (hasSlabBCenter) - slabBCenter_ = rnemdParams->getSlabBCenter(); - else - slabBCenter_ = hmat(2,2) / 2.0; + if (hasSlabBCenter) + slabBCenter_ = rnemdParams->getSlabBCenter(); + else + slabBCenter_ = hmat(2,2) / 2.0; - selectionBstream << "select wrappedz > " - << slabBCenter_ - 0.5*slabWidth_ - << " && wrappedz < " - << slabBCenter_ + 0.5*slabWidth_; + selectionBstream << "select wrappedz > " + << slabBCenter_ - 0.5*slabWidth_ + << " && wrappedz < " + << slabBCenter_ + 0.5*slabWidth_; + selectionB_ = selectionBstream.str(); + } else { + if (hasSphereBRadius_) { + sphereBRadius_ = rnemdParams->getSphereBRadius(); + selectionBstream << "select r > " << sphereBRadius_; selectionB_ = selectionBstream.str(); } else { - if (hasSphereBRadius_) { - sphereBRadius_ = rnemdParams->getSphereBRadius(); - selectionBstream << "select r > " << sphereBRadius_; - selectionB_ = selectionBstream.str(); - } else { - selectionB_ = "select hull"; - hasSelectionB_ = true; - } + selectionB_ = "select hull"; + BisHull_ = true; + hasSelectionB_ = true; } } } - + + // object evaluator: evaluator_.loadScriptString(rnemdObjectSelection_); seleMan_.setSelectionSet(evaluator_.evaluate()); @@ -620,13 +654,13 @@ namespace OpenMD { StuntDouble* sd; - RealType min_val; - bool min_found = false; - StuntDouble* min_sd; + RealType min_val(0.0); + int min_found = 0; + StuntDouble* min_sd = NULL; - RealType max_val; - bool max_found = false; - StuntDouble* max_sd; + RealType max_val(0.0); + int max_found = 0; + StuntDouble* max_sd = NULL; for (sd = seleManA_.beginSelected(selei); sd != NULL; sd = seleManA_.nextSelected(selei)) { @@ -680,7 +714,7 @@ namespace OpenMD { if (!max_found) { max_val = value; max_sd = sd; - max_found = true; + max_found = 1; } else { if (max_val < value) { max_val = value; @@ -742,7 +776,7 @@ namespace OpenMD { if (!min_found) { min_val = value; min_sd = sd; - min_found = true; + min_found = 1; } else { if (min_val > value) { min_val = value; @@ -752,15 +786,18 @@ namespace OpenMD { } #ifdef IS_MPI - int worldRank = MPI::COMM_WORLD.Get_rank(); - - bool my_min_found = min_found; - bool my_max_found = max_found; + int worldRank; + MPI_Comm_rank( MPI_COMM_WORLD, &worldRank); + + int my_min_found = min_found; + int my_max_found = max_found; // Even if we didn't find a minimum, did someone else? - MPI::COMM_WORLD.Allreduce(&my_min_found, &min_found, 1, MPI::BOOL, MPI::LOR); + MPI_Allreduce(&my_min_found, &min_found, 1, MPI_INT, MPI_LOR, + MPI_COMM_WORLD); // Even if we didn't find a maximum, did someone else? - MPI::COMM_WORLD.Allreduce(&my_max_found, &max_found, 1, MPI::BOOL, MPI::LOR); + MPI_Allreduce(&my_max_found, &max_found, 1, MPI_INT, MPI_LOR, + MPI_COMM_WORLD); #endif if (max_found && min_found) { @@ -779,8 +816,8 @@ namespace OpenMD { min_vals.rank = worldRank; // Who had the minimum? - MPI::COMM_WORLD.Allreduce(&min_vals, &min_vals, - 1, MPI::REALTYPE_INT, MPI::MINLOC); + MPI_Allreduce(&min_vals, &min_vals, + 1, MPI_REALTYPE_INT, MPI_MINLOC, MPI_COMM_WORLD); min_val = min_vals.val; if (my_max_found) { @@ -791,8 +828,8 @@ namespace OpenMD { max_vals.rank = worldRank; // Who had the maximum? - MPI::COMM_WORLD.Allreduce(&max_vals, &max_vals, - 1, MPI::REALTYPE_INT, MPI::MAXLOC); + MPI_Allreduce(&max_vals, &max_vals, + 1, MPI_REALTYPE_INT, MPI_MAXLOC, MPI_COMM_WORLD); max_val = max_vals.val; #endif @@ -852,13 +889,13 @@ namespace OpenMD { Vector3d min_vel; Vector3d max_vel = max_sd->getVel(); - MPI::Status status; + MPI_Status status; // point-to-point swap of the velocity vector - MPI::COMM_WORLD.Sendrecv(max_vel.getArrayPointer(), 3, MPI::REALTYPE, - min_vals.rank, 0, - min_vel.getArrayPointer(), 3, MPI::REALTYPE, - min_vals.rank, 0, status); + MPI_Sendrecv(max_vel.getArrayPointer(), 3, MPI_REALTYPE, + min_vals.rank, 0, + min_vel.getArrayPointer(), 3, MPI_REALTYPE, + min_vals.rank, 0, MPI_COMM_WORLD, &status); switch(rnemdFluxType_) { case rnemdKE : @@ -869,11 +906,11 @@ namespace OpenMD { Vector3d max_angMom = max_sd->getJ(); // point-to-point swap of the angular momentum vector - MPI::COMM_WORLD.Sendrecv(max_angMom.getArrayPointer(), 3, - MPI::REALTYPE, min_vals.rank, 1, - min_angMom.getArrayPointer(), 3, - MPI::REALTYPE, min_vals.rank, 1, - status); + MPI_Sendrecv(max_angMom.getArrayPointer(), 3, + MPI_REALTYPE, min_vals.rank, 1, + min_angMom.getArrayPointer(), 3, + MPI_REALTYPE, min_vals.rank, 1, + MPI_COMM_WORLD, &status); max_sd->setJ(min_angMom); } @@ -898,13 +935,13 @@ namespace OpenMD { Vector3d max_vel; Vector3d min_vel = min_sd->getVel(); - MPI::Status status; + MPI_Status status; // point-to-point swap of the velocity vector - MPI::COMM_WORLD.Sendrecv(min_vel.getArrayPointer(), 3, MPI::REALTYPE, - max_vals.rank, 0, - max_vel.getArrayPointer(), 3, MPI::REALTYPE, - max_vals.rank, 0, status); + MPI_Sendrecv(min_vel.getArrayPointer(), 3, MPI_REALTYPE, + max_vals.rank, 0, + max_vel.getArrayPointer(), 3, MPI_REALTYPE, + max_vals.rank, 0, MPI_COMM_WORLD, &status); switch(rnemdFluxType_) { case rnemdKE : @@ -915,11 +952,11 @@ namespace OpenMD { Vector3d max_angMom; // point-to-point swap of the angular momentum vector - MPI::COMM_WORLD.Sendrecv(min_angMom.getArrayPointer(), 3, - MPI::REALTYPE, max_vals.rank, 1, - max_angMom.getArrayPointer(), 3, - MPI::REALTYPE, max_vals.rank, 1, - status); + MPI_Sendrecv(min_angMom.getArrayPointer(), 3, + MPI_REALTYPE, max_vals.rank, 1, + max_angMom.getArrayPointer(), 3, + MPI_REALTYPE, max_vals.rank, 1, + MPI_COMM_WORLD, &status); min_sd->setJ(max_angMom); } @@ -983,7 +1020,6 @@ namespace OpenMD { int selej; Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot(); - RealType time = currentSnap_->getTime(); Mat3x3d hmat = currentSnap_->getHmat(); StuntDouble* sd; @@ -1088,22 +1124,22 @@ namespace OpenMD { Kcw *= 0.5; #ifdef IS_MPI - MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Phx, 1, MPI::REALTYPE, MPI::SUM); - MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Phy, 1, MPI::REALTYPE, MPI::SUM); - MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Phz, 1, MPI::REALTYPE, MPI::SUM); - MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Pcx, 1, MPI::REALTYPE, MPI::SUM); - MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Pcy, 1, MPI::REALTYPE, MPI::SUM); - MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Pcz, 1, MPI::REALTYPE, MPI::SUM); + MPI_Allreduce(MPI_IN_PLACE, &Phx, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, &Phy, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, &Phz, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, &Pcx, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, &Pcy, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, &Pcz, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); - MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Khx, 1, MPI::REALTYPE, MPI::SUM); - MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Khy, 1, MPI::REALTYPE, MPI::SUM); - MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Khz, 1, MPI::REALTYPE, MPI::SUM); - MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Khw, 1, MPI::REALTYPE, MPI::SUM); + MPI_Allreduce(MPI_IN_PLACE, &Khx, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, &Khy, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, &Khz, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, &Khw, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); - MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Kcx, 1, MPI::REALTYPE, MPI::SUM); - MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Kcy, 1, MPI::REALTYPE, MPI::SUM); - MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Kcz, 1, MPI::REALTYPE, MPI::SUM); - MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Kcw, 1, MPI::REALTYPE, MPI::SUM); + MPI_Allreduce(MPI_IN_PLACE, &Kcx, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, &Kcy, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, &Kcz, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, &Kcw, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); #endif //solve coldBin coeff's first @@ -1431,7 +1467,6 @@ namespace OpenMD { int selej; Snapshot* currentSnap_ = info_->getSnapshotManager()->getCurrentSnapshot(); - RealType time = currentSnap_->getTime(); Mat3x3d hmat = currentSnap_->getHmat(); StuntDouble* sd; @@ -1580,18 +1615,22 @@ namespace OpenMD { Kc *= 0.5; #ifdef IS_MPI - MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Ph[0], 3, MPI::REALTYPE, MPI::SUM); - MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Pc[0], 3, MPI::REALTYPE, MPI::SUM); - MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Lh[0], 3, MPI::REALTYPE, MPI::SUM); - MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Lc[0], 3, MPI::REALTYPE, MPI::SUM); - MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Mh, 1, MPI::REALTYPE, MPI::SUM); - MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Kh, 1, MPI::REALTYPE, MPI::SUM); - MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Mc, 1, MPI::REALTYPE, MPI::SUM); - MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &Kc, 1, MPI::REALTYPE, MPI::SUM); - MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, Ih.getArrayPointer(), 9, - MPI::REALTYPE, MPI::SUM); - MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, Ic.getArrayPointer(), 9, - MPI::REALTYPE, MPI::SUM); + MPI_Allreduce(MPI_IN_PLACE, &Ph[0], 3, MPI_REALTYPE, MPI_SUM, + MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, &Pc[0], 3, MPI_REALTYPE, MPI_SUM, + MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, &Lh[0], 3, MPI_REALTYPE, MPI_SUM, + MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, &Lc[0], 3, MPI_REALTYPE, MPI_SUM, + MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, &Mh, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, &Kh, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, &Mc, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, &Kc, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, Ih.getArrayPointer(), 9, + MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, Ic.getArrayPointer(), 9, + MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); #endif @@ -1729,27 +1768,33 @@ 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 (evaluatorA_.hasSurfaceArea()) + areaA = evaluatorA_.getSurfaceArea(); + else { + + 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, + 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(); + "\twithout libqhull. Please rebuild OpenMD with qhull enabled."); + painCave.severity = OPENMD_ERROR; + painCave.isFatal = 1; + simError(); #endif + } } else { if (usePeriodicBoundaryConditions_) { @@ -1765,30 +1810,34 @@ 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 { + 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 +1849,7 @@ namespace OpenMD { areaB = 4.0 * M_PI * pow(sphereBRadius_, 2); } } - + dividingArea_ = min(areaA, areaB); hasDividingArea_ = true; return dividingArea_; @@ -1854,25 +1903,33 @@ 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: 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,68 +1969,81 @@ 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; } } } } - + #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_Allreduce(MPI_IN_PLACE, &binCount[i], + 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, &binMass[i], + 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, binP[i].getArrayPointer(), + 3, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, binL[i].getArrayPointer(), + 3, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, binI[i].getArrayPointer(), + 9, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, &binKE[i], + 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, &binDOF[i], + 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD); + //MPI_Allreduce(MPI_IN_PLACE, &binOmega[i], + // 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); + } + #endif - Vector3d vel; - Vector3d aVel; + Vector3d omega; RealType den; RealType temp; RealType z; @@ -1990,13 +2060,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 +2087,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); @@ -2056,7 +2125,7 @@ namespace OpenMD { painCave.severity = OPENMD_ERROR; simError(); } - } + } } void RNEMD::writeOutputFile() { @@ -2065,7 +2134,9 @@ namespace OpenMD { #ifdef IS_MPI // If we're the root node, should we print out the results - int worldRank = MPI::COMM_WORLD.Get_rank(); + int worldRank; + MPI_Comm_rank( MPI_COMM_WORLD, &worldRank); + if (worldRank == 0) { #endif rnemdFile_.open(rnemdFileName_.c_str(), std::ios::out | std::ios::trunc ); @@ -2196,7 +2267,7 @@ namespace OpenMD { } rnemdFile_ << "#######################################################\n"; - rnemdFile_ << "# Standard Deviations in those quantities follow:\n"; + rnemdFile_ << "# 95% confidence intervals in those quantities follow:\n"; rnemdFile_ << "#######################################################\n"; @@ -2205,9 +2276,9 @@ namespace OpenMD { for (unsigned int i = 0; i < outputMask_.size(); ++i) { if (outputMask_[i]) { if (data_[i].dataType == "RealType") - writeRealStdDev(i,j); + writeRealErrorBars(i,j); else if (data_[i].dataType == "Vector3d") - writeVectorStdDev(i,j); + writeVectorErrorBars(i,j); else { sprintf( painCave.errMsg, "RNEMD found an unknown data type for: %s ", @@ -2278,7 +2349,7 @@ namespace OpenMD { } } - void RNEMD::writeRealStdDev(int index, unsigned int bin) { + void RNEMD::writeRealErrorBars(int index, unsigned int bin) { if (!doRNEMD_) return; assert(index >=0 && index < ENDINDEX); assert(int(bin) < nBins_); @@ -2288,7 +2359,7 @@ namespace OpenMD { count = data_[index].accumulator[bin]->count(); if (count == 0) return; - dynamic_cast(data_[index].accumulator[bin])->getStdDev(s); + dynamic_cast(data_[index].accumulator[bin])->get95percentConfidenceInterval(s); if (! isinf(s) && ! isnan(s)) { rnemdFile_ << "\t" << s; @@ -2301,7 +2372,7 @@ namespace OpenMD { } } - void RNEMD::writeVectorStdDev(int index, unsigned int bin) { + void RNEMD::writeVectorErrorBars(int index, unsigned int bin) { if (!doRNEMD_) return; assert(index >=0 && index < ENDINDEX); assert(int(bin) < nBins_); @@ -2311,7 +2382,7 @@ namespace OpenMD { count = data_[index].accumulator[bin]->count(); if (count == 0) return; - dynamic_cast(data_[index].accumulator[bin])->getStdDev(s); + dynamic_cast(data_[index].accumulator[bin])->get95percentConfidenceInterval(s); if (isinf(s[0]) || isnan(s[0]) || isinf(s[1]) || isnan(s[1]) || isinf(s[2]) || isnan(s[2]) ) {