--- branches/development/src/UseTheForce/ForceField.cpp 2010/12/31 18:31:56 1535 +++ trunk/src/brains/ForceField.cpp 2013/06/17 18:28:30 1880 @@ -35,25 +35,52 @@ * * [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). + * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008). + * [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). + * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). */ /** * @file ForceField.cpp * @author tlin * @date 11/04/2004 - * @time 22:51am * @version 1.0 */ #include -#include "UseTheForce/ForceField.hpp" +#include "brains/ForceField.hpp" #include "utils/simError.h" -#include "utils/Tuple.hpp" + +#include "io/OptionSectionParser.hpp" +#include "io/BaseAtomTypesSectionParser.hpp" +#include "io/DirectionalAtomTypesSectionParser.hpp" +#include "io/AtomTypesSectionParser.hpp" +#include "io/BendTypesSectionParser.hpp" +#include "io/BondTypesSectionParser.hpp" +#include "io/ChargeAtomTypesSectionParser.hpp" +#include "io/EAMAtomTypesSectionParser.hpp" +#include "io/FluctuatingChargeAtomTypesSectionParser.hpp" +#include "io/GayBerneAtomTypesSectionParser.hpp" +#include "io/InversionTypesSectionParser.hpp" +#include "io/LennardJonesAtomTypesSectionParser.hpp" +#include "io/MultipoleAtomTypesSectionParser.hpp" +#include "io/NonBondedInteractionsSectionParser.hpp" +#include "io/PolarizableAtomTypesSectionParser.hpp" +#include "io/SCAtomTypesSectionParser.hpp" +#include "io/ShapeAtomTypesSectionParser.hpp" +#include "io/StickyAtomTypesSectionParser.hpp" +#include "io/StickyPowerAtomTypesSectionParser.hpp" +#include "io/TorsionTypesSectionParser.hpp" + +#include "types/LennardJonesAdapter.hpp" +#include "types/EAMAdapter.hpp" +#include "types/SuttonChenAdapter.hpp" +#include "types/GayBerneAdapter.hpp" +#include "types/StickyAdapter.hpp" + namespace OpenMD { - ForceField::ForceField() { + ForceField::ForceField(std::string ffName) { char* tempPath; tempPath = getenv("FORCE_PARAM_PATH"); @@ -63,7 +90,84 @@ namespace OpenMD { STR_DEFINE(ffPath_, FRC_PATH ); } else { ffPath_ = tempPath; + } + + setForceFieldFileName(ffName + ".frc"); + + /** + * The order of adding section parsers is important. + * + * OptionSectionParser must come first to set options for other + * parsers + * + * DirectionalAtomTypesSectionParser should be added before + * AtomTypesSectionParser, and these two section parsers will + * actually create "real" AtomTypes (AtomTypesSectionParser will + * create AtomType and DirectionalAtomTypesSectionParser will + * create DirectionalAtomType, which is a subclass of AtomType and + * should come first). + * + * Other AtomTypes Section Parsers will not create the "real" + * AtomType, they only add and set some attributes of the AtomType + * (via the Adapters). Thus ordering of these is not important. + * AtomTypesSectionParser should be added before other atom type + * + * The order of BondTypesSectionParser, BendTypesSectionParser and + * TorsionTypesSectionParser, etc. are not important. + */ + + spMan_.push_back(new OptionSectionParser(forceFieldOptions_)); + spMan_.push_back(new BaseAtomTypesSectionParser()); + spMan_.push_back(new DirectionalAtomTypesSectionParser(forceFieldOptions_)); + spMan_.push_back(new AtomTypesSectionParser()); + + spMan_.push_back(new LennardJonesAtomTypesSectionParser(forceFieldOptions_)); + spMan_.push_back(new ChargeAtomTypesSectionParser(forceFieldOptions_)); + spMan_.push_back(new MultipoleAtomTypesSectionParser(forceFieldOptions_)); + spMan_.push_back(new FluctuatingChargeAtomTypesSectionParser(forceFieldOptions_)); + spMan_.push_back(new PolarizableAtomTypesSectionParser(forceFieldOptions_)); + spMan_.push_back(new GayBerneAtomTypesSectionParser(forceFieldOptions_)); + spMan_.push_back(new EAMAtomTypesSectionParser(forceFieldOptions_)); + spMan_.push_back(new SCAtomTypesSectionParser(forceFieldOptions_)); + spMan_.push_back(new ShapeAtomTypesSectionParser(forceFieldOptions_)); + spMan_.push_back(new StickyAtomTypesSectionParser(forceFieldOptions_)); + spMan_.push_back(new StickyPowerAtomTypesSectionParser(forceFieldOptions_)); + + spMan_.push_back(new BondTypesSectionParser(forceFieldOptions_)); + spMan_.push_back(new BendTypesSectionParser(forceFieldOptions_)); + spMan_.push_back(new TorsionTypesSectionParser(forceFieldOptions_)); + spMan_.push_back(new InversionTypesSectionParser(forceFieldOptions_)); + + spMan_.push_back(new NonBondedInteractionsSectionParser(forceFieldOptions_)); + } + + void ForceField::parse(const std::string& filename) { + ifstrstream* ffStream; + + ffStream = openForceFieldFile(filename); + + spMan_.parse(*ffStream, *this); + + ForceField::AtomTypeContainer::MapTypeIterator i; + AtomType* at; + + for (at = atomTypeCont_.beginType(i); at != NULL; + at = atomTypeCont_.nextType(i)) { + + // useBase sets the responsibilities, and these have to be done + // after the atomTypes and Base types have all been scanned: + + std::vector ayb = at->allYourBase(); + if (ayb.size() > 1) { + for (int j = ayb.size()-1; j > 0; j--) { + + ayb[j-1]->useBase(ayb[j]); + + } + } } + + delete ffStream; } /** @@ -123,9 +227,9 @@ namespace OpenMD { std::vector > > foundBonds; - for (i = at1Chain.begin(); i != at1Chain.end(); i++) { + for (i = at1Chain.begin(); i != at1Chain.end(); ++i) { jj = 0; - for (j = at2Chain.begin(); j != at2Chain.end(); j++) { + for (j = at2Chain.begin(); j != at2Chain.end(); ++j) { bondTypeScore = ii + jj; @@ -143,11 +247,10 @@ namespace OpenMD { } - if (foundBonds.size() > 0) { + if (!foundBonds.empty()) { // 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); @@ -205,11 +308,11 @@ namespace OpenMD { std::vector > > foundBends; - for (j = at2Chain.begin(); j != at2Chain.end(); j++) { + for (j = at2Chain.begin(); j != at2Chain.end(); ++j) { ii = 0; - for (i = at1Chain.begin(); i != at1Chain.end(); i++) { + for (i = at1Chain.begin(); i != at1Chain.end(); ++i) { kk = 0; - for (k = at3Chain.begin(); k != at3Chain.end(); k++) { + for (k = at3Chain.begin(); k != at3Chain.end(); ++k) { IKscore = ii + kk; @@ -229,10 +332,8 @@ namespace OpenMD { jj++; } - if (foundBends.size() > 0) { + if (!foundBends.empty()) { 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); @@ -301,13 +402,13 @@ namespace OpenMD { std::vector > > foundTorsions; - for (j = at2Chain.begin(); j != at2Chain.end(); j++) { + for (j = at2Chain.begin(); j != at2Chain.end(); ++j) { kk = 0; - for (k = at3Chain.begin(); k != at3Chain.end(); k++) { + for (k = at3Chain.begin(); k != at3Chain.end(); ++k) { ii = 0; - for (i = at1Chain.begin(); i != at1Chain.end(); i++) { + for (i = at1Chain.begin(); i != at1Chain.end(); ++i) { ll = 0; - for (l = at4Chain.begin(); l != at4Chain.end(); l++) { + for (l = at4Chain.begin(); l != at4Chain.end(); ++l) { ILscore = ii + ll; JKscore = jj + kk; @@ -331,10 +432,8 @@ namespace OpenMD { jj++; } - if (foundTorsions.size() > 0) { + if (!foundTorsions.empty()) { 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); @@ -402,13 +501,13 @@ namespace OpenMD { std::vector > > foundInversions; - for (j = at2Chain.begin(); j != at2Chain.end(); j++) { + for (j = at2Chain.begin(); j != at2Chain.end(); ++j) { kk = 0; - for (k = at3Chain.begin(); k != at3Chain.end(); k++) { + for (k = at3Chain.begin(); k != at3Chain.end(); ++k) { ii = 0; - for (i = at1Chain.begin(); i != at1Chain.end(); i++) { + for (i = at1Chain.begin(); i != at1Chain.end(); ++i) { ll = 0; - for (l = at4Chain.begin(); l != at4Chain.end(); l++) { + for (l = at4Chain.begin(); l != at4Chain.end(); ++l) { Iscore = ii; JKLscore = jj + kk + ll; @@ -432,10 +531,8 @@ namespace OpenMD { jj++; } - if (foundInversions.size() > 0) { + if (!foundInversions.empty()) { 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); @@ -448,18 +545,71 @@ 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.empty()) { + // sort the foundNBI by the score: + std::sort(foundNBI.begin(), foundNBI.end()); + 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, @@ -579,35 +729,29 @@ namespace OpenMD { } RealType ForceField::getRcutFromAtomType(AtomType* at) { - /**@todo */ - GenericData* data; - RealType rcut = 0.0; + 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 = OPENMD_ERROR; - painCave.isFatal = 1; - simError(); - } - } else { - sprintf( painCave.errMsg, "Can not find Parameters for LennardJones\n"); - painCave.severity = OPENMD_ERROR; - painCave.isFatal = 1; - simError(); - } + LennardJonesAdapter lja = LennardJonesAdapter(at); + if (lja.isLennardJones()) { + rcut = 2.5 * lja.getSigma(); } + EAMAdapter ea = EAMAdapter(at); + if (ea.isEAM()) { + rcut = max(rcut, ea.getRcut()); + } + SuttonChenAdapter sca = SuttonChenAdapter(at); + if (sca.isSuttonChen()) { + rcut = max(rcut, 2.0 * sca.getAlpha()); + } + GayBerneAdapter gba = GayBerneAdapter(at); + if (gba.isGayBerne()) { + rcut = max(rcut, 2.5 * sqrt(2.0) * max(gba.getD(), gba.getL())); + } + StickyAdapter sa = StickyAdapter(at); + if (sa.isSticky()) { + rcut = max(rcut, max(sa.getRu(), sa.getRup())); + } + return rcut; }