--- trunk/src/restraints/RestraintForceManager.cpp 2009/09/07 16:31:51 1360 +++ trunk/src/restraints/RestraintForceManager.cpp 2009/11/25 20:02:06 1390 @@ -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) { @@ -78,8 +78,8 @@ namespace oopse { restTime_ = simParam->getStatusTime(); } else { sprintf(painCave.errMsg, - "Restraint warning: If you use restraints without setting\n", - "\tstatusTime, no restraint data will be written to the rest\n", + "Restraint warning: If you use restraints without setting\n" + "\tstatusTime, no restraint data will be written to the rest\n" "\tfile.\n"); painCave.isFatal = 0; simError(); @@ -109,15 +109,44 @@ 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 } MolecularRestraint* rest = new MolecularRestraint(); @@ -148,6 +177,9 @@ namespace oopse { if (stamp[i]->haveRestrainedSwingXAngle()) { rest->setRestrainedSwingXAngle(stamp[i]->getRestrainedSwingXAngle() * M_PI/180.0); } + if (stamp[i]->havePrint()) { + rest->setPrintRestraint(stamp[i]->getPrint()); + } restraints_.push_back(rest); mol->addProperty(new RestraintData("Restraint", rest)); @@ -170,7 +202,7 @@ namespace oopse { SelectionEvaluator evaluator(info); SelectionManager seleMan(info); - + evaluator.loadScriptString(objectSelection); seleMan.setSelectionSet(evaluator.evaluate()); int selectionCount = seleMan.getSelectionCount(); @@ -211,7 +243,11 @@ namespace oopse { } if (stamp[i]->haveRestrainedSwingYAngle()) { rest->setRestrainedSwingYAngle(stamp[i]->getRestrainedSwingYAngle()); - } + } + if (stamp[i]->havePrint()) { + rest->setPrintRestraint(stamp[i]->getPrint()); + } + restraints_.push_back(rest); sd->addProperty(new RestraintData("Restraint", rest)); restrainedObjs_.push_back(sd); @@ -311,20 +347,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(); } @@ -357,7 +393,10 @@ namespace oopse { unscaledPotential_ += mRest->getUnscaledPotential(); restInfo = mRest->getRestraintInfo(); - restInfo_.push_back(restInfo); + + // only collect data on restraints that we're going to print: + if (mRest->getPrintRestraint()) + restInfo_.push_back(restInfo); } for(ro=restrainedObjs_.begin(); ro != restrainedObjs_.end(); ++ro){ @@ -373,20 +412,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(); } @@ -418,7 +457,10 @@ namespace oopse { unscaledPotential_ += oRest->getUnscaledPotential(); restInfo = oRest->getRestraintInfo(); - restInfo_.push_back(restInfo); + + // only collect data on restraints that we're going to print: + if (oRest->getPrintRestraint()) + restInfo_.push_back(restInfo); } return unscaledPotential_ * scalingFactor;