--- trunk/src/UseTheForce/ForceField.cpp 2007/07/06 18:15:03 1151 +++ trunk/src/UseTheForce/ForceField.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). */ /** @@ -47,12 +47,14 @@ * @version 1.0 */ +#include #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 { +namespace OpenMD { ForceField::ForceField() { char* tempPath; @@ -78,7 +80,8 @@ namespace oopse { return atomTypeCont_.find(keys); } - BondType* ForceField::getBondType(const std::string &at1, const std::string &at2) { + BondType* ForceField::getBondType(const std::string &at1, + const std::string &at2) { std::vector keys; keys.push_back(at1); keys.push_back(at2); @@ -88,13 +91,68 @@ namespace oopse { if (bondType) { return bondType; } else { - //if no exact match found, try wild card match - return bondTypeCont_.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(); - BendType* ForceField::getBendType(const std::string &at1, const std::string &at2, + std::vector::iterator i; + std::vector::iterator j; + + int ii = 0; + int jj = 0; + int bondTypeScore; + + std::vector > > foundBonds; + + for (i = at1Chain.begin(); i != at1Chain.end(); i++) { + jj = 0; + for (j = at2Chain.begin(); j != at2Chain.end(); j++) { + + bondTypeScore = ii + jj; + + std::vector myKeys; + myKeys.push_back((*i)->getName()); + myKeys.push_back((*j)->getName()); + + BondType* bondType = bondTypeCont_.find(myKeys); + if (bondType) { + foundBonds.push_back(std::make_pair(bondTypeScore, myKeys)); + } + jj++; + } + ii++; + } + + + if (foundBonds.size() > 0) { + // sort the foundBonds by the score: + std::sort(foundBonds.begin(), foundBonds.end()); + + int bestScore = foundBonds[0].first; + std::vector theKeys = foundBonds[0].second; + + BondType* bestType = bondTypeCont_.find(theKeys); + + return bestType; + } else { + //if no exact match found, try wild card match + return bondTypeCont_.find(keys, wildCardAtomTypeName_); + } + } + } + + BendType* ForceField::getBendType(const std::string &at1, + const std::string &at2, const std::string &at3) { std::vector keys; keys.push_back(at1); @@ -106,55 +164,305 @@ namespace oopse { if (bendType) { return bendType; } else { - //if no exact match found, try wild card match - return bendTypeCont_.find(keys, wildCardAtomTypeName_); + + AtomType* atype1; + AtomType* atype2; + AtomType* atype3; + std::vector at1key; + at1key.push_back(at1); + atype1 = atomTypeCont_.find(at1key); + + std::vector at2key; + at2key.push_back(at2); + atype2 = atomTypeCont_.find(at2key); + + std::vector at3key; + at3key.push_back(at3); + atype3 = atomTypeCont_.find(at3key); + + // query atom types for their chains of responsibility + std::vector at1Chain = atype1->allYourBase(); + std::vector at2Chain = atype2->allYourBase(); + std::vector at3Chain = atype3->allYourBase(); + + std::vector::iterator i; + std::vector::iterator j; + std::vector::iterator k; + + int ii = 0; + int jj = 0; + int kk = 0; + int IKscore; + + std::vector > > foundBends; + + for (j = at2Chain.begin(); j != at2Chain.end(); j++) { + ii = 0; + for (i = at1Chain.begin(); i != at1Chain.end(); i++) { + kk = 0; + for (k = at3Chain.begin(); k != at3Chain.end(); k++) { + + IKscore = ii + kk; + + std::vector myKeys; + myKeys.push_back((*i)->getName()); + myKeys.push_back((*j)->getName()); + myKeys.push_back((*k)->getName()); + + BendType* bendType = bendTypeCont_.find(myKeys); + if (bendType) { + foundBends.push_back( make_tuple3(jj, IKscore, myKeys) ); + } + kk++; + } + ii++; + } + jj++; + } + + if (foundBends.size() > 0) { + std::sort(foundBends.begin(), foundBends.end()); + int jscore = foundBends[0].first; + int ikscore = foundBends[0].second; + std::vector theKeys = foundBends[0].third; + + BendType* bestType = bendTypeCont_.find(theKeys); + return bestType; + } else { + //if no exact match found, try wild card match + return bendTypeCont_.find(keys, wildCardAtomTypeName_); + } } } - TorsionType* ForceField::getTorsionType(const std::string &at1, const std::string &at2, - const std::string &at3, const std::string &at4) { + TorsionType* ForceField::getTorsionType(const std::string &at1, + const std::string &at2, + const std::string &at3, + const std::string &at4) { std::vector keys; keys.push_back(at1); keys.push_back(at2); keys.push_back(at3); keys.push_back(at4); + + //try exact match first TorsionType* torsionType = torsionTypeCont_.find(keys); if (torsionType) { return torsionType; } else { - //if no exact match found, try wild card match - return torsionTypeCont_.find(keys, wildCardAtomTypeName_); - } - - return torsionTypeCont_.find(keys, wildCardAtomTypeName_); + AtomType* atype1; + AtomType* atype2; + AtomType* atype3; + AtomType* atype4; + std::vector at1key; + at1key.push_back(at1); + atype1 = atomTypeCont_.find(at1key); + + std::vector at2key; + at2key.push_back(at2); + atype2 = atomTypeCont_.find(at2key); + + std::vector at3key; + at3key.push_back(at3); + atype3 = atomTypeCont_.find(at3key); + + std::vector at4key; + at4key.push_back(at4); + atype4 = atomTypeCont_.find(at4key); + + // query atom types for their chains of responsibility + std::vector at1Chain = atype1->allYourBase(); + std::vector at2Chain = atype2->allYourBase(); + std::vector at3Chain = atype3->allYourBase(); + std::vector at4Chain = atype4->allYourBase(); + + std::vector::iterator i; + std::vector::iterator j; + std::vector::iterator k; + std::vector::iterator l; + + int ii = 0; + int jj = 0; + int kk = 0; + int ll = 0; + int ILscore; + int JKscore; + + std::vector > > foundTorsions; + + for (j = at2Chain.begin(); j != at2Chain.end(); j++) { + kk = 0; + for (k = at3Chain.begin(); k != at3Chain.end(); k++) { + ii = 0; + for (i = at1Chain.begin(); i != at1Chain.end(); i++) { + ll = 0; + for (l = at4Chain.begin(); l != at4Chain.end(); l++) { + + ILscore = ii + ll; + JKscore = jj + kk; + + std::vector myKeys; + myKeys.push_back((*i)->getName()); + myKeys.push_back((*j)->getName()); + myKeys.push_back((*k)->getName()); + myKeys.push_back((*l)->getName()); + + TorsionType* torsionType = torsionTypeCont_.find(myKeys); + if (torsionType) { + foundTorsions.push_back( make_tuple3(JKscore, ILscore, myKeys) ); + } + ll++; + } + ii++; + } + kk++; + } + jj++; + } + + if (foundTorsions.size() > 0) { + std::sort(foundTorsions.begin(), foundTorsions.end()); + int jkscore = foundTorsions[0].first; + int ilscore = foundTorsions[0].second; + std::vector theKeys = foundTorsions[0].third; + + TorsionType* bestType = torsionTypeCont_.find(theKeys); + return bestType; + } else { + //if no exact match found, try wild card match + return torsionTypeCont_.find(keys, wildCardAtomTypeName_); + } + } } - NonBondedInteractionType* ForceField::getNonBondedInteractionType(const std::string &at1, const std::string &at2) { + InversionType* ForceField::getInversionType(const std::string &at1, + const std::string &at2, + const std::string &at3, + const std::string &at4) { std::vector keys; keys.push_back(at1); keys.push_back(at2); + keys.push_back(at3); + keys.push_back(at4); //try exact match first + InversionType* inversionType = inversionTypeCont_.permutedFindSkippingFirstElement(keys); + if (inversionType) { + return inversionType; + } else { + + AtomType* atype1; + AtomType* atype2; + AtomType* atype3; + AtomType* atype4; + std::vector at1key; + at1key.push_back(at1); + atype1 = atomTypeCont_.find(at1key); + + std::vector at2key; + at2key.push_back(at2); + atype2 = atomTypeCont_.find(at2key); + + std::vector at3key; + at3key.push_back(at3); + atype3 = atomTypeCont_.find(at3key); + + std::vector at4key; + at4key.push_back(at4); + atype4 = atomTypeCont_.find(at4key); + + // query atom types for their chains of responsibility + std::vector at1Chain = atype1->allYourBase(); + std::vector at2Chain = atype2->allYourBase(); + std::vector at3Chain = atype3->allYourBase(); + std::vector at4Chain = atype4->allYourBase(); + + std::vector::iterator i; + std::vector::iterator j; + std::vector::iterator k; + std::vector::iterator l; + + int ii = 0; + int jj = 0; + int kk = 0; + int ll = 0; + int Iscore; + int JKLscore; + + std::vector > > foundInversions; + + for (j = at2Chain.begin(); j != at2Chain.end(); j++) { + kk = 0; + for (k = at3Chain.begin(); k != at3Chain.end(); k++) { + ii = 0; + for (i = at1Chain.begin(); i != at1Chain.end(); i++) { + ll = 0; + for (l = at4Chain.begin(); l != at4Chain.end(); l++) { + + Iscore = ii; + JKLscore = jj + kk + ll; + + std::vector myKeys; + myKeys.push_back((*i)->getName()); + myKeys.push_back((*j)->getName()); + myKeys.push_back((*k)->getName()); + myKeys.push_back((*l)->getName()); + + InversionType* inversionType = inversionTypeCont_.permutedFindSkippingFirstElement(myKeys); + if (inversionType) { + foundInversions.push_back( make_tuple3(Iscore, JKLscore, myKeys) ); + } + ll++; + } + ii++; + } + kk++; + } + jj++; + } + + if (foundInversions.size() > 0) { + std::sort(foundInversions.begin(), foundInversions.end()); + int iscore = foundInversions[0].first; + int jklscore = foundInversions[0].second; + std::vector theKeys = foundInversions[0].third; + + InversionType* bestType = inversionTypeCont_.permutedFindSkippingFirstElement(theKeys); + return bestType; + } else { + //if no exact match found, try wild card match + return inversionTypeCont_.find(keys, wildCardAtomTypeName_); + } + } + } + + NonBondedInteractionType* ForceField::getNonBondedInteractionType(const std::string &at1, const std::string &at2) { + std::vector keys; + keys.push_back(at1); + keys.push_back(at2); + + //try exact match first NonBondedInteractionType* nbiType = nonBondedInteractionTypeCont_.find(keys); if (nbiType) { return nbiType; } else { //if no exact match found, try wild card match return nonBondedInteractionTypeCont_.find(keys, wildCardAtomTypeName_); - } - + } } - - BondType* ForceField::getExactBondType(const std::string &at1, const std::string &at2){ + + BondType* ForceField::getExactBondType(const std::string &at1, + const std::string &at2){ std::vector keys; keys.push_back(at1); keys.push_back(at2); return bondTypeCont_.find(keys); } - - BendType* ForceField::getExactBendType(const std::string &at1, const std::string &at2, + + BendType* ForceField::getExactBendType(const std::string &at1, + const std::string &at2, const std::string &at3){ std::vector keys; keys.push_back(at1); @@ -162,9 +470,11 @@ namespace oopse { keys.push_back(at3); return bendTypeCont_.find(keys); } - - TorsionType* ForceField::getExactTorsionType(const std::string &at1, const std::string &at2, - const std::string &at3, const std::string &at4){ + + TorsionType* ForceField::getExactTorsionType(const std::string &at1, + const std::string &at2, + const std::string &at3, + const std::string &at4){ std::vector keys; keys.push_back(at1); keys.push_back(at2); @@ -172,29 +482,47 @@ namespace oopse { keys.push_back(at4); return torsionTypeCont_.find(keys); } - + + InversionType* ForceField::getExactInversionType(const std::string &at1, + const std::string &at2, + const std::string &at3, + const std::string &at4){ + std::vector keys; + keys.push_back(at1); + keys.push_back(at2); + keys.push_back(at3); + keys.push_back(at4); + return inversionTypeCont_.find(keys); + } + NonBondedInteractionType* ForceField::getExactNonBondedInteractionType(const std::string &at1, const std::string &at2){ std::vector keys; keys.push_back(at1); keys.push_back(at2); return nonBondedInteractionTypeCont_.find(keys); } + - bool ForceField::addAtomType(const std::string &at, AtomType* atomType) { std::vector keys; keys.push_back(at); return atomTypeCont_.add(keys, atomType); } - bool ForceField::addBondType(const std::string &at1, const std::string &at2, BondType* bondType) { + bool ForceField::replaceAtomType(const std::string &at, AtomType* atomType) { std::vector keys; + keys.push_back(at); + return atomTypeCont_.replace(keys, atomType); + } + + bool ForceField::addBondType(const std::string &at1, const std::string &at2, + BondType* bondType) { + std::vector keys; keys.push_back(at1); keys.push_back(at2); - return bondTypeCont_.add(keys, bondType); - + return bondTypeCont_.add(keys, bondType); } - + bool ForceField::addBendType(const std::string &at1, const std::string &at2, const std::string &at3, BendType* bendType) { std::vector keys; @@ -203,9 +531,12 @@ namespace oopse { keys.push_back(at3); return bendTypeCont_.add(keys, bendType); } - - bool ForceField::addTorsionType(const std::string &at1, const std::string &at2, - const std::string &at3, const std::string &at4, TorsionType* torsionType) { + + bool ForceField::addTorsionType(const std::string &at1, + const std::string &at2, + const std::string &at3, + const std::string &at4, + TorsionType* torsionType) { std::vector keys; keys.push_back(at1); keys.push_back(at2); @@ -214,49 +545,62 @@ namespace oopse { return torsionTypeCont_.add(keys, torsionType); } - bool ForceField::addNonBondedInteractionType(const std::string &at1, const std::string &at2, NonBondedInteractionType* nbiType) { + bool ForceField::addInversionType(const std::string &at1, + const std::string &at2, + const std::string &at3, + const std::string &at4, + InversionType* inversionType) { std::vector keys; keys.push_back(at1); keys.push_back(at2); + keys.push_back(at3); + keys.push_back(at4); + return inversionTypeCont_.add(keys, inversionType); + } + + bool ForceField::addNonBondedInteractionType(const std::string &at1, + const std::string &at2, + NonBondedInteractionType* nbiType) { + std::vector keys; + keys.push_back(at1); + keys.push_back(at2); return nonBondedInteractionTypeCont_.add(keys, nbiType); - } - + 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.severity = OPENMD_ERROR; painCave.isFatal = 1; simError(); } } else { sprintf( painCave.errMsg, "Can not find Parameters for LennardJones\n"); - painCave.severity = OOPSE_ERROR; + painCave.severity = OPENMD_ERROR; painCave.isFatal = 1; simError(); } } - return rcut; } + - ifstrstream* ForceField::openForceFieldFile(const std::string& filename) { std::string forceFieldFilename(filename); ifstrstream* ffStream = new ifstrstream(); @@ -278,14 +622,12 @@ 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(); } } - return ffStream; - } void ForceField::setFortranForceOptions(){ @@ -293,4 +635,4 @@ namespace oopse { forceFieldOptions_.makeFortranOptions(theseFortranOptions); setfForceOptions(&theseFortranOptions); } -} //end namespace oopse +} //end namespace OpenMD