--- branches/development/src/nonbonded/InteractionManager.cpp 2011/07/08 20:25:32 1587 +++ branches/development/src/nonbonded/InteractionManager.cpp 2013/05/15 15:09:35 1874 @@ -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,11 +119,7 @@ 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); if (atype1->isLennardJones() && atype2->isLennardJones()) { @@ -143,6 +154,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: @@ -178,7 +193,31 @@ namespace OpenMD { interactions_[key].insert(morse_); 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); + } + interactions_[key].insert(repulsivePower_); + vdwExplicit = true; + } + if (nbiType->isEAM()) { // We found an explicit EAM interaction. // override all other metallic entries for this pair of atom types: @@ -276,10 +315,6 @@ namespace OpenMD { eam_->setCutoffRadius(rcut); } - void InteractionManager::setSwitchingRadius(RealType rswitch) { - electrostatic_->setSwitchingRadius(rswitch); - } - void InteractionManager::doPrePair(InteractionData idat){ if (!initialized_) initialize(); @@ -353,7 +388,7 @@ namespace OpenMD { RealType InteractionManager::getSuggestedCutoffRadius(int *atid) { if (!initialized_) initialize(); - + AtomType* atype = typeMap_[*atid]; pair key = make_pair(atype, atype);