--- branches/development/src/nonbonded/InteractionManager.cpp 2011/05/25 16:20:37 1568 +++ branches/development/src/nonbonded/InteractionManager.cpp 2013/05/15 15:09:35 1874 @@ -35,47 +35,46 @@ * * [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" namespace OpenMD { - InteractionManager* InteractionManager::_instance = NULL; - SimInfo* InteractionManager::info_ = NULL; - bool InteractionManager::initialized_ = false; + InteractionManager::InteractionManager() { - RealType InteractionManager::rCut_ = 0.0; - RealType InteractionManager::rSwitch_ = 0.0; - CutoffMethod InteractionManager::cutoffMethod_ = SHIFTED_FORCE; - SwitchingFunctionType InteractionManager::sft_ = cubic; - RealType InteractionManager::vdwScale_[4] = {1.0, 0.0, 0.0, 0.0}; - RealType InteractionManager::electrostaticScale_[4] = {1.0, 0.0, 0.0, 0.0}; + initialized_ = false; + + lj_ = new LJ(); + gb_ = new GB(); + sticky_ = new Sticky(); + morse_ = new Morse(); + repulsivePower_ = new RepulsivePower(); + eam_ = new EAM(); + sc_ = new SC(); + electrostatic_ = new Electrostatic(); + maw_ = new MAW(); + } - map InteractionManager::typeMap_; - map, set > InteractionManager::interactions_; - - LJ* InteractionManager::lj_ = new LJ(); - GB* InteractionManager::gb_ = new GB(); - Sticky* InteractionManager::sticky_ = new Sticky(); - Morse* InteractionManager::morse_ = new Morse(); - EAM* InteractionManager::eam_ = new EAM(); - SC* InteractionManager::sc_ = new SC(); - Electrostatic* InteractionManager::electrostatic_ = new Electrostatic(); - MAW* InteractionManager::maw_ = new MAW(); - SwitchingFunction* InteractionManager::switcher_ = new SwitchingFunction(); - - InteractionManager* InteractionManager::Instance() { - if (!_instance) { - _instance = new InteractionManager(); - } - return _instance; + 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_); @@ -84,48 +83,28 @@ namespace OpenMD { eam_->setForceField(forceField_); sc_->setForceField(forceField_); morse_->setForceField(forceField_); + electrostatic_->setSimInfo(info_); electrostatic_->setForceField(forceField_); maw_->setForceField(forceField_); - - ForceFieldOptions& fopts = forceField_->getForceFieldOptions(); - - // Force fields can set options on how to scale van der Waals and electrostatic - // interactions for atoms connected via bonds, bends and torsions - // in this case the topological distance between atoms is: - // 0 = topologically unconnected - // 1 = bonded together - // 2 = connected via a bend - // 3 = connected via a torsion - - vdwScale_[0] = 1.0; - vdwScale_[1] = fopts.getvdw12scale(); - vdwScale_[2] = fopts.getvdw13scale(); - vdwScale_[3] = fopts.getvdw14scale(); - - electrostaticScale_[0] = 1.0; - electrostaticScale_[1] = fopts.getelectrostatic12scale(); - electrostaticScale_[2] = fopts.getelectrostatic13scale(); - electrostaticScale_[3] = fopts.getelectrostatic14scale(); + 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(); @@ -140,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()) { @@ -179,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: @@ -214,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: @@ -278,190 +281,51 @@ namespace OpenMD { } - // make sure every pair of atom types in this simulation has a - // non-bonded interaction: + // Make sure every pair of atom types in this simulation has a + // non-bonded interaction. If not, just inform the user. set simTypes = info_->getSimulatedAtomTypes(); set::iterator it, jt; + for (it = simTypes.begin(); it != simTypes.end(); ++it) { atype1 = (*it); - for (jt = simTypes.begin(); jt != simTypes.end(); ++jt) { + for (jt = it; jt != simTypes.end(); ++jt) { atype2 = (*jt); key = make_pair(atype1, atype2); if (interactions_[key].size() == 0) { sprintf( painCave.errMsg, - "InteractionManager unable to find an appropriate non-bonded\n" - "\tinteraction for atom types %s - %s\n", + "InteractionManager could not find a matching non-bonded\n" + "\tinteraction for atom types %s - %s\n" + "\tProceeding without this interaction.\n", atype1->getName().c_str(), atype2->getName().c_str()); painCave.severity = OPENMD_INFO; - painCave.isFatal = 1; + painCave.isFatal = 0; simError(); } } } - setupCutoffs(); - setupSwitching(); - - //int ljsp = cutoffMethod_ == SHIFTED_POTENTIAL ? 1 : 0; - //int ljsf = cutoffMethod_ == SHIFTED_FORCE ? 1 : 0; - //notifyFortranCutoffs(&rCut_, &rSwitch_, &ljsp, &ljsf); - initialized_ = true; } - - /** - * setupCutoffs - * - * Sets the values of cutoffRadius and cutoffMethod - * - * cutoffRadius : realType - * If the cutoffRadius was explicitly set, use that value. - * If the cutoffRadius was not explicitly set: - * Are there electrostatic atoms? Use 12.0 Angstroms. - * No electrostatic atoms? Poll the atom types present in the - * simulation for suggested cutoff values (e.g. 2.5 * sigma). - * Use the maximum suggested value that was found. - * - * cutoffMethod : (one of HARD, SWITCHED, SHIFTED_FORCE, SHIFTED_POTENTIAL) - * If cutoffMethod was explicitly set, use that choice. - * If cutoffMethod was not explicitly set, use SHIFTED_FORCE - */ - void InteractionManager::setupCutoffs() { - - Globals* simParams_ = info_->getSimParams(); - - if (simParams_->haveCutoffRadius()) { - rCut_ = simParams_->getCutoffRadius(); - } else { - if (info_->usesElectrostaticAtoms()) { - sprintf(painCave.errMsg, - "InteractionManager::setupCutoffs: No value was set for the cutoffRadius.\n" - "\tOpenMD will use a default value of 12.0 angstroms" - "\tfor the cutoffRadius.\n"); - painCave.isFatal = 0; - painCave.severity = OPENMD_INFO; - simError(); - rCut_ = 12.0; - } else { - RealType thisCut; - set::iterator i; - set atomTypes; - atomTypes = info_->getSimulatedAtomTypes(); - for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { - thisCut = getSuggestedCutoffRadius((*i)); - rCut_ = max(thisCut, rCut_); - } - sprintf(painCave.errMsg, - "InteractionManager::setupCutoffs: No value was set for the cutoffRadius.\n" - "\tOpenMD will use %lf angstroms.\n", - rCut_); - painCave.isFatal = 0; - painCave.severity = OPENMD_INFO; - simError(); - } - } - map stringToCutoffMethod; - stringToCutoffMethod["HARD"] = HARD; - stringToCutoffMethod["SWITCHED"] = SWITCHED; - stringToCutoffMethod["SHIFTED_POTENTIAL"] = SHIFTED_POTENTIAL; - stringToCutoffMethod["SHIFTED_FORCE"] = SHIFTED_FORCE; - - if (simParams_->haveCutoffMethod()) { - string cutMeth = toUpperCopy(simParams_->getCutoffMethod()); - map::iterator i; - i = stringToCutoffMethod.find(cutMeth); - if (i == stringToCutoffMethod.end()) { - sprintf(painCave.errMsg, - "InteractionManager::setupCutoffs: Could not find chosen cutoffMethod %s\n" - "\tShould be one of: " - "HARD, SWITCHED, SHIFTED_POTENTIAL, or SHIFTED_FORCE\n", - cutMeth.c_str()); - painCave.isFatal = 1; - painCave.severity = OPENMD_ERROR; - simError(); - } else { - cutoffMethod_ = i->second; - } - } else { - sprintf(painCave.errMsg, - "InteractionManager::setupCutoffs: No value was set for the cutoffMethod.\n" - "\tOpenMD will use SHIFTED_FORCE.\n"); - painCave.isFatal = 0; - painCave.severity = OPENMD_INFO; - simError(); - cutoffMethod_ = SHIFTED_FORCE; - } - } - - - /** - * setupSwitching - * - * Sets the values of switchingRadius and - * If the switchingRadius was explicitly set, use that value (but check it) - * If the switchingRadius was not explicitly set: use 0.85 * cutoffRadius_ - */ - void InteractionManager::setupSwitching() { - Globals* simParams_ = info_->getSimParams(); - - if (simParams_->haveSwitchingRadius()) { - rSwitch_ = simParams_->getSwitchingRadius(); - if (rSwitch_ > rCut_) { - sprintf(painCave.errMsg, - "InteractionManager::setupSwitching: switchingRadius (%f) is larger than cutoffRadius(%f)\n", - rSwitch_, rCut_); - painCave.isFatal = 1; - painCave.severity = OPENMD_ERROR; - simError(); - } - } else { - rSwitch_ = 0.85 * rCut_; - sprintf(painCave.errMsg, - "InteractionManager::setupSwitching: No value was set for the switchingRadius.\n" - "\tOpenMD will use a default value of 85 percent of the cutoffRadius.\n" - "\tswitchingRadius = %f. for this simulation\n", rSwitch_); - painCave.isFatal = 0; - painCave.severity = OPENMD_WARNING; - simError(); - } + void InteractionManager::setCutoffRadius(RealType rcut) { - if (simParams_->haveSwitchingFunctionType()) { - string funcType = simParams_->getSwitchingFunctionType(); - toUpper(funcType); - if (funcType == "CUBIC") { - sft_ = cubic; - } else { - if (funcType == "FIFTH_ORDER_POLYNOMIAL") { - sft_ = fifth_order_poly; - } else { - // throw error - sprintf( painCave.errMsg, - "InteractionManager::setupSwitching : Unknown switchingFunctionType. (Input file specified %s .)\n" - "\tswitchingFunctionType must be one of: " - "\"cubic\" or \"fifth_order_polynomial\".", - funcType.c_str() ); - painCave.isFatal = 1; - painCave.severity = OPENMD_ERROR; - simError(); - } - } - } - - switcher_->setSwitchType(sft_); - switcher_->setSwitch(rSwitch_, rCut_); + electrostatic_->setCutoffRadius(rcut); + eam_->setCutoffRadius(rcut); } void InteractionManager::doPrePair(InteractionData idat){ if (!initialized_) initialize(); + // excluded interaction, so just return + if (idat.excluded) return; + set::iterator it; - for (it = interactions_[ *(idat.atypes) ].begin(); - it != interactions_[ *(idat.atypes) ].end(); ++it){ + for (it = interactions_[ idat.atypes ].begin(); + it != interactions_[ idat.atypes ].end(); ++it){ if ((*it)->getFamily() == METALLIC_FAMILY) { dynamic_cast(*it)->calcDensity(idat); } @@ -489,26 +353,17 @@ 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; - } + for (it = interactions_[ idat.atypes ].begin(); + it != interactions_[ idat.atypes ].end(); ++it) { - void InteractionManager::doSkipCorrection(InteractionData idat){ + // electrostatics still has to worry about indirect + // contributions from excluded pairs of atoms: - if (!initialized_) initialize(); - - set::iterator it; - - for (it = interactions_[ *(idat.atypes) ].begin(); - it != interactions_[ *(idat.atypes) ].end(); ++it){ - if ((*it)->getFamily() == ELECTROSTATIC_FAMILY) { - dynamic_cast(*it)->calcSkipCorrection(idat); + if (!idat.excluded || (*it)->getFamily() == ELECTROSTATIC_FAMILY) { + (*it)->calcForce(idat); } } @@ -533,7 +388,7 @@ namespace OpenMD { RealType InteractionManager::getSuggestedCutoffRadius(int *atid) { if (!initialized_) initialize(); - + AtomType* atype = typeMap_[*atid]; pair key = make_pair(atype, atype);