--- branches/development/src/nonbonded/InteractionManager.cpp 2011/11/22 20:38:56 1665 +++ trunk/src/nonbonded/InteractionManager.cpp 2013/08/07 15:24:16 1925 @@ -35,7 +35,7 @@ * * [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). + * [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). */ @@ -59,8 +59,22 @@ namespace OpenMD { 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_); @@ -75,24 +89,48 @@ namespace OpenMD { repulsivePower_->setForceField(forceField_); ForceField::AtomTypeContainer* atomTypes = forceField_->getAtomTypes(); + int nTypes = atomTypes->size(); + iHash_.resize(nTypes); + interactions_.resize(nTypes); ForceField::AtomTypeContainer::MapTypeIterator i1, i2; AtomType* atype1; AtomType* atype2; - pair key; - pair::iterator, bool> ret; + int atid1, atid2; + + // We only need to worry about the types that are actually in the simulation: + + set atypes = info_->getSimulatedAtomTypes(); + + lj_->setSimulatedAtomTypes(atypes); + gb_->setSimulatedAtomTypes(atypes); + sticky_->setSimulatedAtomTypes(atypes); + eam_->setSimulatedAtomTypes(atypes); + sc_->setSimulatedAtomTypes(atypes); + morse_->setSimulatedAtomTypes(atypes); + electrostatic_->setSimInfo(info_); + electrostatic_->setSimulatedAtomTypes(atypes); + maw_->setSimulatedAtomTypes(atypes); + repulsivePower_->setSimulatedAtomTypes(atypes); + + set::iterator at; + + for (at = atypes.begin(); at != atypes.end(); ++at) { - for (atype1 = atomTypes->beginType(i1); atype1 != NULL; - atype1 = atomTypes->nextType(i1)) { + //for (atype1 = atomTypes->beginType(i1); atype1 != NULL; + // atype1 = atomTypes->nextType(i1)) { - // add it to the map: - AtomTypeProperties atp = atype1->getATP(); - + atype1 = *at; + atid1 = atype1->getIdent(); + iHash_[atid1].resize(nTypes); + interactions_[atid1].resize(nTypes); + + // add it to the map: pair::iterator,bool> ret; - ret = typeMap_.insert( pair(atp.ident, atype1) ); + ret = typeMap_.insert( pair(atid1, 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,40 +142,46 @@ namespace OpenMD { map::iterator it1, it2; for (it1 = typeMap_.begin(); it1 != typeMap_.end(); ++it1) { atype1 = (*it1).second; + atid1 = atype1->getIdent(); for( it2 = typeMap_.begin(); it2 != typeMap_.end(); ++it2) { atype2 = (*it2).second; + atid2 = atype2->getIdent(); + + iHash_[atid1][atid2] = 0; - bool vdwExplicit = false; - bool metExplicit = false; - bool hbExplicit = false; - - key = make_pair(atype1, atype2); - if (atype1->isLennardJones() && atype2->isLennardJones()) { - interactions_[key].insert(lj_); + interactions_[atid1][atid2].insert(lj_); + iHash_[atid1][atid2] |= LJ_PAIR; } if (atype1->isElectrostatic() && atype2->isElectrostatic() ) { - interactions_[key].insert(electrostatic_); + interactions_[atid1][atid2].insert(electrostatic_); + iHash_[atid1][atid2] |= ELECTROSTATIC_PAIR; } if (atype1->isSticky() && atype2->isSticky() ) { - interactions_[key].insert(sticky_); + interactions_[atid1][atid2].insert(sticky_); + iHash_[atid1][atid2] |= STICKY_PAIR; } if (atype1->isStickyPower() && atype2->isStickyPower() ) { - interactions_[key].insert(sticky_); + interactions_[atid1][atid2].insert(sticky_); + iHash_[atid1][atid2] |= STICKY_PAIR; } if (atype1->isEAM() && atype2->isEAM() ) { - interactions_[key].insert(eam_); + interactions_[atid1][atid2].insert(eam_); + iHash_[atid1][atid2] |= EAM_PAIR; } if (atype1->isSC() && atype2->isSC() ) { - interactions_[key].insert(sc_); + interactions_[atid1][atid2].insert(sc_); + iHash_[atid1][atid2] |= SC_PAIR; } if (atype1->isGayBerne() && atype2->isGayBerne() ) { - interactions_[key].insert(gb_); + interactions_[atid1][atid2].insert(gb_); + iHash_[atid1][atid2] |= GB_PAIR; } if ((atype1->isGayBerne() && atype2->isLennardJones()) || (atype1->isLennardJones() && atype2->isGayBerne())) { - interactions_[key].insert(gb_); + interactions_[atid1][atid2].insert(gb_); + iHash_[atid1][atid2] |= GB_PAIR; } // look for an explicitly-set non-bonded interaction type using the @@ -146,16 +190,24 @@ 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: set::iterator it; - for (it = interactions_[key].begin(); - it != interactions_[key].end(); ++it) { + for (it = interactions_[atid1][atid2].begin(); + it != interactions_[atid1][atid2].end(); ++it) { InteractionFamily ifam = (*it)->getFamily(); - if (ifam == VANDERWAALS_FAMILY) interactions_[key].erase(*it); + if (ifam == VANDERWAALS_FAMILY) { + interactions_[atid1][atid2].erase(*it); + iHash_[atid1][atid2] ^= (*it)->getHash(); + } } - interactions_[key].insert(lj_); + interactions_[atid1][atid2].insert(lj_); + iHash_[atid1][atid2] |= LJ_PAIR; vdwExplicit = true; } @@ -173,12 +225,16 @@ namespace OpenMD { // We found an explicit Morse 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) { + for (it = interactions_[atid1][atid2].begin(); + it != interactions_[atid1][atid2].end(); ++it) { InteractionFamily ifam = (*it)->getFamily(); - if (ifam == VANDERWAALS_FAMILY) interactions_[key].erase(*it); + if (ifam == VANDERWAALS_FAMILY) { + interactions_[atid1][atid2].erase(*it); + iHash_[atid1][atid2] ^= (*it)->getHash(); + } } - interactions_[key].insert(morse_); + interactions_[atid1][atid2].insert(morse_); + iHash_[atid1][atid2] |= MORSE_PAIR; vdwExplicit = true; } @@ -196,12 +252,16 @@ namespace OpenMD { // 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) { + for (it = interactions_[atid1][atid2].begin(); + it != interactions_[atid1][atid2].end(); ++it) { InteractionFamily ifam = (*it)->getFamily(); - if (ifam == VANDERWAALS_FAMILY) interactions_[key].erase(*it); + if (ifam == VANDERWAALS_FAMILY) { + interactions_[atid1][atid2].erase(*it); + iHash_[atid1][atid2] ^= (*it)->getHash(); + } } - interactions_[key].insert(repulsivePower_); + interactions_[atid1][atid2].insert(repulsivePower_); + iHash_[atid1][atid2] |= REPULSIVEPOWER_PAIR; vdwExplicit = true; } @@ -210,12 +270,16 @@ namespace OpenMD { // We found an explicit EAM interaction. // override all other metallic entries for this pair of atom types: set::iterator it; - for (it = interactions_[key].begin(); - it != interactions_[key].end(); ++it) { + for (it = interactions_[atid1][atid2].begin(); + it != interactions_[atid1][atid2].end(); ++it) { InteractionFamily ifam = (*it)->getFamily(); - if (ifam == METALLIC_FAMILY) interactions_[key].erase(*it); + if (ifam == METALLIC_FAMILY) { + interactions_[atid1][atid2].erase(*it); + iHash_[atid1][atid2] ^= (*it)->getHash(); + } } - interactions_[key].insert(eam_); + interactions_[atid1][atid2].insert(eam_); + iHash_[atid1][atid2] |= EAM_PAIR; metExplicit = true; } @@ -233,12 +297,16 @@ namespace OpenMD { // We found an explicit Sutton-Chen interaction. // override all other metallic entries for this pair of atom types: set::iterator it; - for (it = interactions_[key].begin(); - it != interactions_[key].end(); ++it) { + for (it = interactions_[atid1][atid2].begin(); + it != interactions_[atid1][atid2].end(); ++it) { InteractionFamily ifam = (*it)->getFamily(); - if (ifam == METALLIC_FAMILY) interactions_[key].erase(*it); + if (ifam == METALLIC_FAMILY) { + interactions_[atid1][atid2].erase(*it); + iHash_[atid1][atid2] ^= (*it)->getHash(); + } } - interactions_[key].insert(sc_); + interactions_[atid1][atid2].insert(sc_); + iHash_[atid1][atid2] |= SC_PAIR; metExplicit = true; } @@ -256,12 +324,16 @@ namespace OpenMD { // We found an explicit MAW 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) { + for (it = interactions_[atid1][atid2].begin(); + it != interactions_[atid1][atid2].end(); ++it) { InteractionFamily ifam = (*it)->getFamily(); - if (ifam == VANDERWAALS_FAMILY) interactions_[key].erase(*it); + if (ifam == VANDERWAALS_FAMILY) { + interactions_[atid1][atid2].erase(*it); + iHash_[atid1][atid2] ^= (*it)->getHash(); + } } - interactions_[key].insert(maw_); + interactions_[atid1][atid2].insert(maw_); + iHash_[atid1][atid2] |= MAW_PAIR; vdwExplicit = true; } } @@ -277,11 +349,12 @@ namespace OpenMD { for (it = simTypes.begin(); it != simTypes.end(); ++it) { atype1 = (*it); + atid1 = atype1->getIdent(); for (jt = it; jt != simTypes.end(); ++jt) { atype2 = (*jt); - key = make_pair(atype1, atype2); + atid1 = atype1->getIdent(); - if (interactions_[key].size() == 0) { + if (interactions_[atid1][atid2].size() == 0) { sprintf( painCave.errMsg, "InteractionManager could not find a matching non-bonded\n" "\tinteraction for atom types %s - %s\n" @@ -303,104 +376,138 @@ namespace OpenMD { eam_->setCutoffRadius(rcut); } - void InteractionManager::setSwitchingRadius(RealType rswitch) { - electrostatic_->setSwitchingRadius(rswitch); - } - - void InteractionManager::doPrePair(InteractionData idat){ + void InteractionManager::doPrePair(InteractionData &idat){ if (!initialized_) initialize(); // excluded interaction, so just return if (idat.excluded) return; - set::iterator it; + int& iHash = iHash_[idat.atid1][idat.atid2]; - for (it = interactions_[ idat.atypes ].begin(); - it != interactions_[ idat.atypes ].end(); ++it){ - if ((*it)->getFamily() == METALLIC_FAMILY) { - dynamic_cast(*it)->calcDensity(idat); - } - } + 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; } - void InteractionManager::doPreForce(SelfData sdat){ + void InteractionManager::doPreForce(SelfData &sdat){ if (!initialized_) initialize(); - pair key = make_pair(sdat.atype, sdat.atype); - set::iterator it; + int& iHash = iHash_[sdat.atid][sdat.atid]; + + 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_[atid1][atid2].begin(); it != interactions_[atid1][atid2].end(); ++it){ + // if ((*it)->getFamily() == METALLIC_FAMILY) { + // dynamic_cast(*it)->calcFunctional(sdat); + // } + // } return; } - void InteractionManager::doPair(InteractionData idat){ + void InteractionManager::doPair(InteractionData &idat){ if (!initialized_) initialize(); - set::iterator it; + int& iHash = iHash_[idat.atid1][idat.atid2]; - for (it = interactions_[ idat.atypes ].begin(); - it != interactions_[ idat.atypes ].end(); ++it) { + 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: - // electrostatics still has to worry about indirect - // contributions from excluded pairs of atoms: + if (idat.excluded) return; - if (!idat.excluded || (*it)->getFamily() == ELECTROSTATIC_FAMILY) { - (*it)->calcForce(idat); - } - } + 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); + + // 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; } - void InteractionManager::doSelfCorrection(SelfData sdat){ + void InteractionManager::doSelfCorrection(SelfData &sdat){ if (!initialized_) initialize(); - pair key = make_pair(sdat.atype, sdat.atype); - set::iterator it; + int& iHash = iHash_[sdat.atid][sdat.atid]; - 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); + + // set::iterator it; + + // for (it = interactions_[atid1][atid2].begin(); it != interactions_[atid1][atid2].end(); ++it){ + // if ((*it)->getFamily() == ELECTROSTATIC_FAMILY) { + // dynamic_cast(*it)->calcSelfCorrection(sdat); + // } + // } return; } + void InteractionManager::doReciprocalSpaceSum(RealType &pot){ + if (!initialized_) initialize(); + electrostatic_->ReciprocalSpaceSum(pot); + } + RealType InteractionManager::getSuggestedCutoffRadius(int *atid) { if (!initialized_) initialize(); - + AtomType* atype = typeMap_[*atid]; - pair key = make_pair(atype, atype); set::iterator it; RealType cutoff = 0.0; - for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it) - cutoff = max(cutoff, (*it)->getSuggestedCutoffRadius(key)); + for (it = interactions_[*atid][*atid].begin(); + it != interactions_[*atid][*atid].end(); + ++it) + cutoff = max(cutoff, (*it)->getSuggestedCutoffRadius(make_pair(atype, atype))); return cutoff; } RealType InteractionManager::getSuggestedCutoffRadius(AtomType* atype) { if (!initialized_) initialize(); - pair key = make_pair(atype, atype); + int atid = atype->getIdent(); + set::iterator it; RealType cutoff = 0.0; - for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it) - cutoff = max(cutoff, (*it)->getSuggestedCutoffRadius(key)); + for (it = interactions_[atid][atid].begin(); + it != interactions_[atid][atid].end(); ++it) + cutoff = max(cutoff, (*it)->getSuggestedCutoffRadius(make_pair(atype, atype))); return cutoff; } } //end namespace OpenMD