--- branches/development/src/UseTheForce/ForceField.cpp 2010/12/31 18:31:56 1535 +++ branches/development/src/UseTheForce/ForceField.cpp 2011/11/22 20:38:56 1665 @@ -36,7 +36,8 @@ * [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). + * [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). + * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). */ /** @@ -448,18 +449,73 @@ namespace OpenMD { } 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_); - } + 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,