--- trunk/src/restraints/RestraintForceManager.cpp 2009/10/07 20:49:50 1364 +++ trunk/src/restraints/RestraintForceManager.cpp 2010/07/09 19:29:05 1464 @@ -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,6 +28,15 @@ * 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, 24107 (2008). + * [4] Vardeman & Gezelter, in progress (2009). */ #include @@ -45,7 +45,7 @@ #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" @@ -54,7 +54,7 @@ #endif -namespace oopse { +namespace OpenMD { RestraintForceManager::RestraintForceManager(SimInfo* info): ForceManager(info) { @@ -89,6 +89,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()); @@ -109,17 +111,56 @@ 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_WORLD.Get_rank(); + 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_WORLD.Get_rank(); + if (info_->getMolToProc(molIndex) == myrank) { + +#endif + MolecularRestraint* rest = new MolecularRestraint(); std::string restPre("mol_"); @@ -155,7 +196,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; @@ -191,7 +234,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 +275,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,9 +305,9 @@ namespace oopse { currRestTime_ = currSnapshot_->getTime(); } - void RestraintForceManager::calcForces(bool needPotential, bool needStress){ + void RestraintForceManager::calcForces(){ - ForceManager::calcForces(needPotential, needStress); + ForceManager::calcForces(); RealType restPot_local, restPot; restPot_local = doRestraints(1.0); @@ -318,20 +360,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(); } @@ -345,7 +387,7 @@ namespace oopse { std::vector forces; for(sd = (*rm)->beginIntegrableObject(ioi); sd != NULL; - sd = (*rm)->nextIntegrableObject(ioi)) { + sd = (*rm)->nextIntegrableObject(ioi)) { struc.push_back(sd->getPos()); } @@ -383,20 +425,20 @@ 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(); }