--- trunk/src/restraints/RestraintForceManager.cpp 2009/10/07 20:49:50 1364 +++ trunk/src/restraints/RestraintForceManager.cpp 2014/10/16 19:13:51 2024 @@ -6,19 +6,10 @@ * redistribute this software in source and binary code form, provided * that the following conditions are met: * - * 1. Acknowledgement of the program authors must be made in any - * publication of scientific results based in part on use of the - * program. An acceptable form of acknowledgement is citation of - * the article in which the program was described (Matthew - * A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher - * J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented - * Parallel Simulation Engine for Molecular Dynamics," - * J. Comput. Chem. 26, pp. 252-271 (2005)) - * - * 2. Redistributions of source code must retain the above copyright + * 1. Redistributions of source code must retain the above copyright * notice, this list of conditions and the following disclaimer. * - * 3. Redistributions in binary form must reproduce the above copyright + * 2. Redistributions in binary form must reproduce the above copyright * notice, this list of conditions and the following disclaimer in the * documentation and/or other materials provided with the * distribution. @@ -37,25 +28,37 @@ * arising out of the use of or inability to use software, even if the * University of Notre Dame has been advised of the possibility of * such damages. + * + * SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your + * research, please cite the appropriate papers when you publish your + * work. Good starting points are: + * + * [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005). + * [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006). + * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008). + * [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). + * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). */ +#ifdef IS_MPI +#include +#endif + +#include "config.h" #include + #include "restraints/RestraintForceManager.hpp" #include "restraints/MolecularRestraint.hpp" #include "restraints/ObjectRestraint.hpp" #include "io/RestReader.hpp" #include "utils/simError.h" -#include "utils/OOPSEConstant.hpp" +#include "utils/PhysicalConstants.hpp" #include "utils/StringUtils.hpp" #include "selection/SelectionEvaluator.hpp" #include "selection/SelectionManager.hpp" -#ifdef IS_MPI -#include -#endif +namespace OpenMD { -namespace oopse { - RestraintForceManager::RestraintForceManager(SimInfo* info): ForceManager(info) { // order of affairs: @@ -89,6 +92,8 @@ namespace oopse { int nRestraintStamps = simParam->getNRestraintStamps(); std::vector stamp = simParam->getRestraintStamps(); + std::vector stuntDoubleIndex; + for (int i = 0; i < nRestraintStamps; i++){ std::string myType = toUpperCopy(stamp[i]->getType()); @@ -96,7 +101,6 @@ namespace oopse { if (myType.compare("MOLECULAR")==0){ int molIndex; - std::vector ref; Vector3d refCom; if (!stamp[i]->haveMolIndex()) { @@ -109,17 +113,59 @@ namespace oopse { molIndex = stamp[i]->getMolIndex(); } - Molecule* mol = info_->getMoleculeByGlobalIndex(molIndex); - - if (mol == NULL) { + if (molIndex < 0) { sprintf(painCave.errMsg, - "Restraint Error: A molecular restraint was specified, but\n" - "\tno molecule was found with global index %d.\n", - molIndex); + "Restraint Error: A molecular restraint was specified\n" + "\twith a molIndex that was less than 0\n"); painCave.isFatal = 1; simError(); + } + if (molIndex >= info_->getNGlobalMolecules()) { + sprintf(painCave.errMsg, + "Restraint Error: A molecular restraint was specified with\n" + "\ta molIndex that was greater than the total number of molecules\n"); + painCave.isFatal = 1; + simError(); + } + + Molecule* mol = info_->getMoleculeByGlobalIndex(molIndex); + + if (mol == NULL) { +#ifdef IS_MPI + // getMoleculeByGlobalIndex returns a NULL in parallel if + // this proc doesn't have the molecule. Do a quick check to + // make sure another processor is supposed to have it. + + int myrank; + MPI_Comm_rank( MPI_COMM_WORLD, &myrank); + + if (info_->getMolToProc(molIndex) == myrank) { + + // If we were supposed to have it but got a null, then freak out. +#endif + + sprintf(painCave.errMsg, + "Restraint Error: A molecular restraint was specified, but\n" + "\tno molecule was found with global index %d.\n", + molIndex); + painCave.isFatal = 1; + simError(); + +#ifdef IS_MPI + } +#endif } + +#ifdef IS_MPI + // only handle this molecular restraint if this processor owns the + // molecule + int myrank; + MPI_Comm_rank( MPI_COMM_WORLD, &myrank); + if (info_->getMolToProc(molIndex) == myrank) { + +#endif + MolecularRestraint* rest = new MolecularRestraint(); std::string restPre("mol_"); @@ -155,7 +201,9 @@ namespace oopse { restraints_.push_back(rest); mol->addProperty(new RestraintData("Restraint", rest)); restrainedMols_.push_back(mol); - +#ifdef IS_MPI + } +#endif } else if (myType.compare("OBJECT") == 0) { std::string objectSelection; @@ -177,12 +225,18 @@ namespace oopse { evaluator.loadScriptString(objectSelection); seleMan.setSelectionSet(evaluator.evaluate()); int selectionCount = seleMan.getSelectionCount(); - + +#ifdef IS_MPI + MPI_Allreduce(MPI_IN_PLACE, &selectionCount, 1, MPI_INT, MPI_SUM, + MPI_COMM_WORLD); +#endif + sprintf(painCave.errMsg, "Restraint Info: The specified restraint objectSelection,\n" "\t\t%s\n" "\twill result in %d integrable objects being\n" "\trestrained.\n", objectSelection.c_str(), selectionCount); + painCave.severity = OPENMD_INFO; painCave.isFatal = 0; simError(); @@ -191,7 +245,8 @@ namespace oopse { for (sd = seleMan.beginSelected(selei); sd != NULL; sd = seleMan.nextSelected(selei)) { - + stuntDoubleIndex.push_back(sd->getGlobalIntegrableObjectIndex()); + ObjectRestraint* rest = new ObjectRestraint(); if (stamp[i]->haveDisplacementSpringConstant()) { @@ -231,16 +286,14 @@ namespace oopse { // are times when it won't use restraints at all, so only open the // restraint file if we are actually using restraints: - if (simParam->getUseRestraints()) { + if (simParam->getUseRestraints()) { std::string refFile = simParam->getRestraint_file(); - RestReader* rr = new RestReader(info, refFile); - + RestReader* rr = new RestReader(info, refFile, stuntDoubleIndex); rr->readReferenceStructure(); } restOutput_ = getPrefix(info_->getFinalConfigFileName()) + ".rest"; restOut = new RestWriter(info_, restOutput_.c_str(), restraints_); - if(!restOut){ sprintf(painCave.errMsg, "Restraint error: Failed to create RestWriter\n"); painCave.isFatal = 1; @@ -263,22 +316,20 @@ namespace oopse { currRestTime_ = currSnapshot_->getTime(); } - void RestraintForceManager::calcForces(bool needPotential, bool needStress){ + void RestraintForceManager::calcForces(){ - ForceManager::calcForces(needPotential, needStress); - RealType restPot_local, restPot; + ForceManager::calcForces(); + RealType restPot(0.0); - restPot_local = doRestraints(1.0); + restPot = doRestraints(1.0); #ifdef IS_MPI - MPI::COMM_WORLD.Allreduce(&restPot_local, &restPot, 1, - MPI::REALTYPE, MPI::SUM); -#else - restPot = restPot_local; + MPI_Allreduce(MPI_IN_PLACE, &restPot, 1, MPI_REALTYPE, MPI_SUM, + MPI_COMM_WORLD); #endif + currSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot(); - currSnapshot_->statData[Stats::LONG_RANGE_POTENTIAL] += restPot; - currSnapshot_->statData[Stats::VHARM] = restPot; + currSnapshot_->setRestraintPotential(restPot); //write out forces and current positions of restrained molecules if (currSnapshot_->getTime() >= currRestTime_){ @@ -293,7 +344,6 @@ namespace oopse { Molecule::IntegrableObjectIterator ioi; MolecularRestraint* mRest; StuntDouble* sd; - RealType pTot; std::vector::const_iterator ro; ObjectRestraint* oRest; @@ -318,20 +368,20 @@ namespace oopse { if (mRest == NULL) { sprintf( painCave.errMsg, "Can not cast RestraintData to MolecularRestraint\n"); - painCave.severity = OOPSE_ERROR; + painCave.severity = OPENMD_ERROR; painCave.isFatal = 1; simError(); } } else { sprintf( painCave.errMsg, "Can not cast GenericData to RestraintData\n"); - painCave.severity = OOPSE_ERROR; + painCave.severity = OPENMD_ERROR; painCave.isFatal = 1; simError(); } } else { sprintf( painCave.errMsg, "Can not find Restraint for RestrainedObject\n"); - painCave.severity = OOPSE_ERROR; + painCave.severity = OPENMD_ERROR; painCave.isFatal = 1; simError(); } @@ -343,9 +393,10 @@ namespace oopse { std::vector struc; std::vector forces; + for(sd = (*rm)->beginIntegrableObject(ioi); sd != NULL; - sd = (*rm)->nextIntegrableObject(ioi)) { + sd = (*rm)->nextIntegrableObject(ioi)) { struc.push_back(sd->getPos()); } @@ -361,12 +412,11 @@ namespace oopse { index++; } - unscaledPotential_ += mRest->getUnscaledPotential(); + unscaledPotential_ += mRest->getUnscaledPotential(); - restInfo = mRest->getRestraintInfo(); - // only collect data on restraints that we're going to print: if (mRest->getPrintRestraint()) + restInfo = mRest->getRestraintInfo(); restInfo_.push_back(restInfo); } @@ -383,27 +433,26 @@ namespace oopse { if (oRest == NULL) { sprintf( painCave.errMsg, "Can not cast RestraintData to ObjectRestraint\n"); - painCave.severity = OOPSE_ERROR; + painCave.severity = OPENMD_ERROR; painCave.isFatal = 1; simError(); } } else { sprintf( painCave.errMsg, "Can not cast GenericData to RestraintData\n"); - painCave.severity = OOPSE_ERROR; + painCave.severity = OPENMD_ERROR; painCave.isFatal = 1; simError(); } } else { sprintf( painCave.errMsg, "Can not find Restraint for RestrainedObject\n"); - painCave.severity = OOPSE_ERROR; + painCave.severity = OPENMD_ERROR; painCave.isFatal = 1; simError(); } // phew. At this point, we should have the pointer to the // correct Object restraint in the variable oRest. - oRest->setScaleFactor(scalingFactor); Vector3d pos = (*ro)->getPos(); @@ -427,10 +476,9 @@ namespace oopse { unscaledPotential_ += oRest->getUnscaledPotential(); - restInfo = oRest->getRestraintInfo(); - // only collect data on restraints that we're going to print: if (oRest->getPrintRestraint()) + restInfo = oRest->getRestraintInfo(); restInfo_.push_back(restInfo); }