--- branches/development/src/nonbonded/InteractionManager.cpp 2011/06/21 06:34:35 1586 +++ branches/development/src/nonbonded/InteractionManager.cpp 2013/06/06 15:43:35 1877 @@ -35,8 +35,9 @@ * * [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). */ #include "nonbonded/InteractionManager.hpp" @@ -51,14 +52,29 @@ namespace OpenMD { gb_ = new GB(); sticky_ = new Sticky(); morse_ = new Morse(); + repulsivePower_ = new RepulsivePower(); eam_ = new EAM(); sc_ = new SC(); electrostatic_ = new Electrostatic(); maw_ = new MAW(); } + InteractionManager::~InteractionManager() { + delete lj_; + delete gb_; + delete sticky_; + delete morse_; + delete repulsivePower_; + delete eam_; + delete sc_; + delete electrostatic_; + delete maw_; + } + void InteractionManager::initialize() { - + + if (initialized_) return; + ForceField* forceField_ = info_->getForceField(); lj_->setForceField(forceField_); @@ -70,26 +86,25 @@ namespace OpenMD { electrostatic_->setSimInfo(info_); electrostatic_->setForceField(forceField_); maw_->setForceField(forceField_); + repulsivePower_->setForceField(forceField_); ForceField::AtomTypeContainer* atomTypes = forceField_->getAtomTypes(); ForceField::AtomTypeContainer::MapTypeIterator i1, i2; AtomType* atype1; AtomType* atype2; pair key; - pair::iterator, bool> ret; for (atype1 = atomTypes->beginType(i1); atype1 != NULL; atype1 = atomTypes->nextType(i1)) { // add it to the map: - AtomTypeProperties atp = atype1->getATP(); pair::iterator,bool> ret; - ret = typeMap_.insert( pair(atp.ident, atype1) ); + ret = typeMap_.insert( pair(atype1->getIdent(), atype1) ); if (ret.second == false) { sprintf( painCave.errMsg, "InteractionManager already had a previous entry with ident %d\n", - atp.ident); + atype1->getIdent()); painCave.severity = OPENMD_INFO; painCave.isFatal = 0; simError(); @@ -104,37 +119,43 @@ namespace OpenMD { for( it2 = typeMap_.begin(); it2 != typeMap_.end(); ++it2) { atype2 = (*it2).second; - - bool vdwExplicit = false; - bool metExplicit = false; - bool hbExplicit = false; - + key = make_pair(atype1, atype2); + + iHash_[key] = 0; if (atype1->isLennardJones() && atype2->isLennardJones()) { interactions_[key].insert(lj_); + iHash_[key] |= LJ_PAIR; } if (atype1->isElectrostatic() && atype2->isElectrostatic() ) { interactions_[key].insert(electrostatic_); + iHash_[key] |= ELECTROSTATIC_PAIR; } if (atype1->isSticky() && atype2->isSticky() ) { interactions_[key].insert(sticky_); + iHash_[key] |= STICKY_PAIR; } if (atype1->isStickyPower() && atype2->isStickyPower() ) { interactions_[key].insert(sticky_); + iHash_[key] |= STICKY_PAIR; } if (atype1->isEAM() && atype2->isEAM() ) { interactions_[key].insert(eam_); + iHash_[key] |= EAM_PAIR; } if (atype1->isSC() && atype2->isSC() ) { interactions_[key].insert(sc_); + iHash_[key] |= SC_PAIR; } if (atype1->isGayBerne() && atype2->isGayBerne() ) { interactions_[key].insert(gb_); + iHash_[key] |= GB_PAIR; } if ((atype1->isGayBerne() && atype2->isLennardJones()) || (atype1->isLennardJones() && atype2->isGayBerne())) { interactions_[key].insert(gb_); + iHash_[key] |= GB_PAIR; } // look for an explicitly-set non-bonded interaction type using the @@ -143,6 +164,10 @@ namespace OpenMD { if (nbiType != NULL) { + bool vdwExplicit = false; + bool metExplicit = false; + // bool hbExplicit = false; + if (nbiType->isLennardJones()) { // We found an explicit Lennard-Jones interaction. // override all other vdw entries for this pair of atom types: @@ -150,9 +175,13 @@ namespace OpenMD { for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it) { InteractionFamily ifam = (*it)->getFamily(); - if (ifam == VANDERWAALS_FAMILY) interactions_[key].erase(*it); + if (ifam == VANDERWAALS_FAMILY) { + interactions_[key].erase(*it); + // work on iHash here; + } } interactions_[key].insert(lj_); + iHash_[key] |= LJ_PAIR; vdwExplicit = true; } @@ -173,12 +202,44 @@ namespace OpenMD { for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it) { InteractionFamily ifam = (*it)->getFamily(); - if (ifam == VANDERWAALS_FAMILY) interactions_[key].erase(*it); + if (ifam == VANDERWAALS_FAMILY) { + interactions_[key].erase(*it); + // work on iHash here; + } } interactions_[key].insert(morse_); + iHash_[key] |= MORSE_PAIR; vdwExplicit = true; } + + if (nbiType->isRepulsivePower()) { + if (vdwExplicit) { + sprintf( painCave.errMsg, + "InteractionManager::initialize found more than one " + "explicit \n" + "\tvan der Waals interaction for atom types %s - %s\n", + atype1->getName().c_str(), atype2->getName().c_str()); + painCave.severity = OPENMD_ERROR; + painCave.isFatal = 1; + simError(); + } + // We found an explicit RepulsivePower interaction. + // override all other vdw entries for this pair of atom types: + set::iterator it; + for (it = interactions_[key].begin(); + it != interactions_[key].end(); ++it) { + InteractionFamily ifam = (*it)->getFamily(); + if (ifam == VANDERWAALS_FAMILY) { + interactions_[key].erase(*it); + // work on iHash here; + } + } + interactions_[key].insert(repulsivePower_); + iHash_[key] |= REPULSIVEPOWER_PAIR; + vdwExplicit = true; + } + if (nbiType->isEAM()) { // We found an explicit EAM interaction. // override all other metallic entries for this pair of atom types: @@ -186,9 +247,13 @@ namespace OpenMD { for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it) { InteractionFamily ifam = (*it)->getFamily(); - if (ifam == METALLIC_FAMILY) interactions_[key].erase(*it); + if (ifam == METALLIC_FAMILY) { + interactions_[key].erase(*it); + // work on iHash here; + } } interactions_[key].insert(eam_); + iHash_[key] |= EAM_PAIR; metExplicit = true; } @@ -209,9 +274,13 @@ namespace OpenMD { for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it) { InteractionFamily ifam = (*it)->getFamily(); - if (ifam == METALLIC_FAMILY) interactions_[key].erase(*it); + if (ifam == METALLIC_FAMILY) { + interactions_[key].erase(*it); + // work on iHash here; + } } interactions_[key].insert(sc_); + iHash_[key] |= SC_PAIR; metExplicit = true; } @@ -232,9 +301,13 @@ namespace OpenMD { for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it) { InteractionFamily ifam = (*it)->getFamily(); - if (ifam == VANDERWAALS_FAMILY) interactions_[key].erase(*it); + if (ifam == VANDERWAALS_FAMILY) { + interactions_[key].erase(*it); + // work on iHash here; + } } interactions_[key].insert(maw_); + iHash_[key] |= MAW_PAIR; vdwExplicit = true; } } @@ -271,26 +344,31 @@ namespace OpenMD { } void InteractionManager::setCutoffRadius(RealType rcut) { + electrostatic_->setCutoffRadius(rcut); eam_->setCutoffRadius(rcut); } - void InteractionManager::setSwitchingRadius(RealType rswitch) { - electrostatic_->setSwitchingRadius(rswitch); - } - void InteractionManager::doPrePair(InteractionData idat){ if (!initialized_) initialize(); - set::iterator it; + // excluded interaction, so just return + if (idat.excluded) return; - for (it = interactions_[ idat.atypes ].begin(); - it != interactions_[ idat.atypes ].end(); ++it){ - if ((*it)->getFamily() == METALLIC_FAMILY) { - dynamic_cast(*it)->calcDensity(idat); - } - } + int& iHash = iHash_[idat.atypes]; + + if ((iHash & EAM_PAIR) != 0) eam_->calcDensity(idat); + if ((iHash & SC_PAIR) != 0) sc_->calcDensity(idat); + + // set::iterator it; + + // for (it = interactions_[ idat.atypes ].begin(); + // it != interactions_[ idat.atypes ].end(); ++it){ + // if ((*it)->getFamily() == METALLIC_FAMILY) { + // dynamic_cast(*it)->calcDensity(idat); + // } + // } return; } @@ -299,14 +377,20 @@ namespace OpenMD { if (!initialized_) initialize(); - pair key = make_pair(sdat.atype, sdat.atype); - set::iterator it; + // pair key = make_pair(sdat.atype, sdat.atype); + + int& iHash = iHash_[ make_pair(sdat.atype, sdat.atype) ]; + + if ((iHash & EAM_PAIR) != 0) eam_->calcFunctional(sdat); + if ((iHash & SC_PAIR) != 0) sc_->calcFunctional(sdat); + + // set::iterator it; - for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it){ - if ((*it)->getFamily() == METALLIC_FAMILY) { - dynamic_cast(*it)->calcFunctional(sdat); - } - } + // for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it){ + // if ((*it)->getFamily() == METALLIC_FAMILY) { + // dynamic_cast(*it)->calcFunctional(sdat); + // } + // } return; } @@ -314,28 +398,37 @@ namespace OpenMD { void InteractionManager::doPair(InteractionData idat){ if (!initialized_) initialize(); - - set::iterator it; - for (it = interactions_[ idat.atypes ].begin(); - it != interactions_[ idat.atypes ].end(); ++it) - (*it)->calcForce(idat); - - return; - } + int& iHash = iHash_[idat.atypes]; - void InteractionManager::doSkipCorrection(InteractionData idat){ + if ((iHash & ELECTROSTATIC_PAIR) != 0) electrostatic_->calcForce(idat); + + // electrostatics still has to worry about indirect + // contributions from excluded pairs of atoms, but nothing else does: - if (!initialized_) initialize(); + if (idat.excluded) return; - set::iterator it; + if ((iHash & LJ_PAIR) != 0) lj_->calcForce(idat); + if ((iHash & GB_PAIR) != 0) gb_->calcForce(idat); + if ((iHash & STICKY_PAIR) != 0) sticky_->calcForce(idat); + if ((iHash & MORSE_PAIR) != 0) morse_->calcForce(idat); + if ((iHash & REPULSIVEPOWER_PAIR) != 0) repulsivePower_->calcForce(idat); + if ((iHash & EAM_PAIR) != 0) eam_->calcForce(idat); + if ((iHash & SC_PAIR) != 0) sc_->calcForce(idat); + if ((iHash & MAW_PAIR) != 0) maw_->calcForce(idat); - for (it = interactions_[ idat.atypes ].begin(); - it != interactions_[ idat.atypes ].end(); ++it){ - if ((*it)->getFamily() == ELECTROSTATIC_FAMILY) { - dynamic_cast(*it)->calcSkipCorrection(idat); - } - } + // set::iterator it; + + // for (it = interactions_[ idat.atypes ].begin(); + // it != interactions_[ idat.atypes ].end(); ++it) { + + // // electrostatics still has to worry about indirect + // // contributions from excluded pairs of atoms: + + // if (!idat.excluded || (*it)->getFamily() == ELECTROSTATIC_FAMILY) { + // (*it)->calcForce(idat); + // } + // } return; } @@ -344,21 +437,26 @@ namespace OpenMD { if (!initialized_) initialize(); - pair key = make_pair(sdat.atype, sdat.atype); - set::iterator it; + int& iHash = iHash_[ make_pair(sdat.atype, sdat.atype) ]; - for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it){ - if ((*it)->getFamily() == ELECTROSTATIC_FAMILY) { - dynamic_cast(*it)->calcSelfCorrection(sdat); - } - } + if ((iHash & ELECTROSTATIC_PAIR) != 0) electrostatic_->calcSelfCorrection(sdat); + + + // pair key = make_pair(sdat.atype, sdat.atype); + // set::iterator it; + + // for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it){ + // if ((*it)->getFamily() == ELECTROSTATIC_FAMILY) { + // dynamic_cast(*it)->calcSelfCorrection(sdat); + // } + // } return; } RealType InteractionManager::getSuggestedCutoffRadius(int *atid) { if (!initialized_) initialize(); - + AtomType* atype = typeMap_[*atid]; pair key = make_pair(atype, atype);