--- trunk/src/UseTheForce/ForceField.cpp 2008/10/13 21:35:22 1303 +++ branches/development/src/UseTheForce/ForceField.cpp 2012/05/18 21:44:02 1710 @@ -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,16 @@ * 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] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). + * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). */ /** @@ -51,35 +52,46 @@ #include "UseTheForce/ForceField.hpp" #include "utils/simError.h" #include "utils/Tuple.hpp" -#include "UseTheForce/DarkSide/atype_interface.h" -#include "UseTheForce/DarkSide/fForceOptions_interface.h" -#include "UseTheForce/DarkSide/switcheroo_interface.h" -namespace oopse { +#include "types/LennardJonesAdapter.hpp" +namespace OpenMD { + ForceField::ForceField() { + char* tempPath; tempPath = getenv("FORCE_PARAM_PATH"); - + if (tempPath == NULL) { //convert a macro from compiler to a string in c++ STR_DEFINE(ffPath_, FRC_PATH ); } else { ffPath_ = tempPath; } - } - - - ForceField::~ForceField() { - deleteAtypes(); - deleteSwitch(); } + /** + * getAtomType by string + * + * finds the requested atom type in this force field using the string + * name of the atom type. + */ AtomType* ForceField::getAtomType(const std::string &at) { std::vector keys; keys.push_back(at); return atomTypeCont_.find(keys); } + /** + * getAtomType by ident + * + * finds the requested atom type in this force field using the + * integer ident instead of the string name of the atom type. + */ + AtomType* ForceField::getAtomType(int ident) { + std::string at = atypeIdentToName.find(ident)->second; + return getAtomType(at); + } + BondType* ForceField::getBondType(const std::string &at1, const std::string &at2) { std::vector keys; @@ -439,6 +451,7 @@ namespace oopse { } NonBondedInteractionType* ForceField::getNonBondedInteractionType(const std::string &at1, const std::string &at2) { + std::vector keys; keys.push_back(at1); keys.push_back(at2); @@ -448,9 +461,63 @@ namespace oopse { if (nbiType) { return nbiType; } else { - //if no exact match found, try wild card match - return nonBondedInteractionTypeCont_.find(keys, wildCardAtomTypeName_); - } + AtomType* atype1; + AtomType* atype2; + std::vector at1key; + at1key.push_back(at1); + atype1 = atomTypeCont_.find(at1key); + + std::vector at2key; + at2key.push_back(at2); + atype2 = atomTypeCont_.find(at2key); + + // query atom types for their chains of responsibility + std::vector at1Chain = atype1->allYourBase(); + std::vector at2Chain = atype2->allYourBase(); + + std::vector::iterator i; + std::vector::iterator j; + + int ii = 0; + int jj = 0; + int nbiTypeScore; + + std::vector > > foundNBI; + + for (i = at1Chain.begin(); i != at1Chain.end(); i++) { + jj = 0; + for (j = at2Chain.begin(); j != at2Chain.end(); j++) { + + nbiTypeScore = ii + jj; + + std::vector myKeys; + myKeys.push_back((*i)->getName()); + myKeys.push_back((*j)->getName()); + + NonBondedInteractionType* nbiType = nonBondedInteractionTypeCont_.find(myKeys); + if (nbiType) { + foundNBI.push_back(std::make_pair(nbiTypeScore, myKeys)); + } + jj++; + } + ii++; + } + + + if (foundNBI.size() > 0) { + // sort the foundNBI by the score: + std::sort(foundNBI.begin(), foundNBI.end()); + + int bestScore = foundNBI[0].first; + std::vector theKeys = foundNBI[0].second; + + NonBondedInteractionType* bestType = nonBondedInteractionTypeCont_.find(theKeys); + return bestType; + } else { + //if no exact match found, try wild card match + return nonBondedInteractionTypeCont_.find(keys, wildCardAtomTypeName_); + } + } } BondType* ForceField::getExactBondType(const std::string &at1, @@ -506,12 +573,14 @@ namespace oopse { bool ForceField::addAtomType(const std::string &at, AtomType* atomType) { std::vector keys; keys.push_back(at); + atypeIdentToName[atomType->getIdent()] = at; return atomTypeCont_.add(keys, atomType); } bool ForceField::replaceAtomType(const std::string &at, AtomType* atomType) { std::vector keys; keys.push_back(at); + atypeIdentToName[atomType->getIdent()] = at; return atomTypeCont_.replace(keys, atomType); } @@ -568,34 +637,11 @@ namespace oopse { } RealType ForceField::getRcutFromAtomType(AtomType* at) { - /**@todo */ - GenericData* data; RealType rcut = 0.0; - if (at->isLennardJones()) { - data = at->getPropertyByName("LennardJones"); - if (data != NULL) { - LJParamGenericData* ljData = dynamic_cast(data); - - if (ljData != NULL) { - LJParam ljParam = ljData->getData(); - - //by default use 2.5*sigma as cutoff radius - rcut = 2.5 * ljParam.sigma; - - } else { - sprintf( painCave.errMsg, - "Can not cast GenericData to LJParam\n"); - painCave.severity = OOPSE_ERROR; - painCave.isFatal = 1; - simError(); - } - } else { - sprintf( painCave.errMsg, "Can not find Parameters for LennardJones\n"); - painCave.severity = OOPSE_ERROR; - painCave.isFatal = 1; - simError(); - } + LennardJonesAdapter lja = LennardJonesAdapter(at); + if (lja.isLennardJones()) { + rcut = 2.5 * lja.getSigma(); } return rcut; } @@ -622,7 +668,7 @@ namespace oopse { "\tHave you tried setting the FORCE_PARAM_PATH environment " "variable?\n", forceFieldFilename.c_str() ); - painCave.severity = OOPSE_ERROR; + painCave.severity = OPENMD_ERROR; painCave.isFatal = 1; simError(); } @@ -630,9 +676,4 @@ namespace oopse { return ffStream; } - void ForceField::setFortranForceOptions(){ - ForceOptions theseFortranOptions; - forceFieldOptions_.makeFortranOptions(theseFortranOptions); - setfForceOptions(&theseFortranOptions); - } -} //end namespace oopse +} //end namespace OpenMD