--- branches/development/src/nonbonded/InteractionManager.cpp 2010/10/03 22:18:59 1505 +++ branches/development/src/nonbonded/InteractionManager.cpp 2013/02/20 15:39:39 1850 @@ -35,44 +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 { - - bool InteractionManager::initialized_ = false; - ForceField* InteractionManager::forceField_ = NULL; - InteractionManager* InteractionManager::_instance = NULL; - map InteractionManager::typeMap_; - map, set > InteractionManager::interactions_; - - InteractionManager* InteractionManager::Instance() { - if (!_instance) { - _instance = new InteractionManager(); - } - return _instance; - } - void InteractionManager::initialize() { - + InteractionManager::InteractionManager() { + + initialized_ = false; + lj_ = new LJ(); gb_ = new GB(); sticky_ = new Sticky(); + morse_ = new Morse(); + repulsivePower_ = new RepulsivePower(); eam_ = new EAM(); sc_ = new SC(); - morse_ = new Morse(); electrostatic_ = new Electrostatic(); + maw_ = new MAW(); + } + void InteractionManager::initialize() { + + if (initialized_) return; + + ForceField* forceField_ = info_->getForceField(); + lj_->setForceField(forceField_); gb_->setForceField(forceField_); sticky_->setForceField(forceField_); eam_->setForceField(forceField_); sc_->setForceField(forceField_); morse_->setForceField(forceField_); + electrostatic_->setSimInfo(info_); electrostatic_->setForceField(forceField_); + maw_->setForceField(forceField_); + repulsivePower_->setForceField(forceField_); ForceField::AtomTypeContainer* atomTypes = forceField_->getAtomTypes(); ForceField::AtomTypeContainer::MapTypeIterator i1, i2; @@ -85,14 +87,13 @@ namespace OpenMD { 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(); @@ -110,7 +111,7 @@ namespace OpenMD { bool vdwExplicit = false; bool metExplicit = false; - bool hbExplicit = false; + // bool hbExplicit = false; key = make_pair(atype1, atype2); @@ -143,272 +144,240 @@ namespace OpenMD { // look for an explicitly-set non-bonded interaction type using the // two atom types. NonBondedInteractionType* nbiType = forceField_->getNonBondedInteractionType(atype1->getName(), atype2->getName()); + + if (nbiType != NULL) { - 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) { - InteractionFamily ifam = (*it)->getFamily(); - if (ifam == VANDERWAALS_FAMILY) interactions_[key].erase(*it); + 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) { + InteractionFamily ifam = (*it)->getFamily(); + if (ifam == VANDERWAALS_FAMILY) interactions_[key].erase(*it); + } + interactions_[key].insert(lj_); + vdwExplicit = true; } - interactions_[key].insert(lj_); - vdwExplicit = true; - } - - if (nbiType->isMorse()) { - 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(); + + if (nbiType->isMorse()) { + 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 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) { + InteractionFamily ifam = (*it)->getFamily(); + if (ifam == VANDERWAALS_FAMILY) interactions_[key].erase(*it); + } + interactions_[key].insert(morse_); + vdwExplicit = true; } - // 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) { - InteractionFamily ifam = (*it)->getFamily(); - if (ifam == VANDERWAALS_FAMILY) interactions_[key].erase(*it); - } - interactions_[key].insert(morse_); - vdwExplicit = true; - } - if (nbiType->isEAM()) { - // 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) { - InteractionFamily ifam = (*it)->getFamily(); - if (ifam == METALLIC_FAMILY) interactions_[key].erase(*it); + 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; } - interactions_[key].insert(eam_); - metExplicit = true; - } - - if (nbiType->isSC()) { - if (metExplicit) { - sprintf( painCave.errMsg, - "InteractionManager::initialize found more than one explicit\n" - "\tmetallic interaction for atom types %s - %s\n", - atype1->getName().c_str(), atype2->getName().c_str()); - painCave.severity = OPENMD_ERROR; - painCave.isFatal = 1; - simError(); + + + if (nbiType->isEAM()) { + // 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) { + InteractionFamily ifam = (*it)->getFamily(); + if (ifam == METALLIC_FAMILY) interactions_[key].erase(*it); + } + interactions_[key].insert(eam_); + metExplicit = true; } - // 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) { - InteractionFamily ifam = (*it)->getFamily(); - if (ifam == METALLIC_FAMILY) interactions_[key].erase(*it); + + if (nbiType->isSC()) { + if (metExplicit) { + sprintf( painCave.errMsg, + "InteractionManager::initialize found more than one " + "explicit\n" + "\tmetallic 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 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) { + InteractionFamily ifam = (*it)->getFamily(); + if (ifam == METALLIC_FAMILY) interactions_[key].erase(*it); + } + interactions_[key].insert(sc_); + metExplicit = true; } - interactions_[key].insert(sc_); - metExplicit = true; + + if (nbiType->isMAW()) { + 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 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) { + InteractionFamily ifam = (*it)->getFamily(); + if (ifam == VANDERWAALS_FAMILY) interactions_[key].erase(*it); + } + interactions_[key].insert(maw_); + vdwExplicit = true; + } } } } - // make sure every pair of atom types has a non-bonded interaction: - for (atype1 = atomTypes->beginType(i1); atype1 != NULL; - atype1 = atomTypes->nextType(i1)) { - for (atype2 = atomTypes->beginType(i2); atype2 != NULL; - atype2 = atomTypes->nextType(i2)) { + + // 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 = 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(); } } } - } - - void InteractionManager::doPrePair(int *atid1, int *atid2, RealType *rij, RealType *rho_i_at_j, RealType *rho_j_at_i){ + + initialized_ = true; + } + + void InteractionManager::setCutoffRadius(RealType rcut) { + electrostatic_->setCutoffRadius(rcut); + eam_->setCutoffRadius(rcut); + } + + void InteractionManager::doPrePair(InteractionData idat){ + if (!initialized_) initialize(); - - DensityData ddat; + + // excluded interaction, so just return + if (idat.excluded) return; - ddat.atype1 = typeMap_[*atid1]; - ddat.atype2 = typeMap_[*atid2]; - ddat.rij = *rij; - ddat.rho_i_at_j = *rho_i_at_j; - ddat.rho_j_at_i = *rho_j_at_i; - - pair key = make_pair(ddat.atype1, ddat.atype2); set::iterator it; - for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it){ + for (it = interactions_[ idat.atypes ].begin(); + it != interactions_[ idat.atypes ].end(); ++it){ if ((*it)->getFamily() == METALLIC_FAMILY) { - dynamic_cast(*it)->calcDensity(ddat); + dynamic_cast(*it)->calcDensity(idat); } } return; } - void InteractionManager::doPreForce(int *atid, RealType *rho, RealType *frho, RealType *dfrhodrho){ + void InteractionManager::doPreForce(SelfData sdat){ if (!initialized_) initialize(); - - FunctionalData fdat; - - fdat.atype = typeMap_[*atid]; - fdat.rho = *rho; - fdat.frho = *frho; - fdat.dfrhodrho = *dfrhodrho; - - pair key = make_pair(fdat.atype, fdat.atype); + + pair key = make_pair(sdat.atype, sdat.atype); set::iterator it; for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it){ if ((*it)->getFamily() == METALLIC_FAMILY) { - dynamic_cast(*it)->calcFunctional(fdat); + dynamic_cast(*it)->calcFunctional(sdat); } } return; } - void InteractionManager::doPair(int *atid1, int *atid2, RealType *d, RealType *r, RealType *r2, RealType *rcut, RealType *sw, RealType *vdwMult,RealType *electroMult, RealType *pot, RealType *vpair, RealType *f1, RealType *eFrame1, RealType *eFrame2, RealType *A1, RealType *A2, RealType *t1, RealType *t2, RealType *rho1, RealType *rho2, RealType *dfrho1, RealType *dfrho2, RealType *fshift1, RealType *fshift2){ + void InteractionManager::doPair(InteractionData idat){ if (!initialized_) initialize(); - - InteractionData idat; - - idat.atype1 = typeMap_[*atid1]; - idat.atype2 = typeMap_[*atid2]; - idat.d = Vector3d(d); - idat.rij = *r; - idat.r2 = *r2; - idat.rcut = *rcut; - idat.sw = *sw; - idat.vdwMult = *vdwMult; - idat.electroMult = *electroMult; - idat.pot = *pot; - idat.vpair = *vpair; - idat.f1 = Vector3d(f1); - idat.eFrame1 = Mat3x3d(eFrame1); - idat.eFrame2 = Mat3x3d(eFrame2); - idat.A1 = RotMat3x3d(A1); - idat.A2 = RotMat3x3d(A2); - idat.t1 = Vector3d(t1); - idat.t2 = Vector3d(t2); - idat.rho1 = *rho1; - idat.rho2 = *rho2; - idat.dfrho1 = *dfrho1; - idat.dfrho2 = *dfrho2; - idat.fshift1 = *fshift1; - idat.fshift2 = *fshift2; - pair key = make_pair(idat.atype1, idat.atype2); set::iterator it; - for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it) - (*it)->calcForce(idat); - - f1[0] = idat.f1.x(); - f1[1] = idat.f1.y(); - f1[2] = idat.f1.z(); - - t1[0] = idat.t1.x(); - t1[1] = idat.t1.y(); - t1[2] = idat.t1.z(); - - t2[0] = idat.t2.x(); - t2[1] = idat.t2.y(); - t2[2] = idat.t2.z(); + for (it = interactions_[ idat.atypes ].begin(); + it != interactions_[ idat.atypes ].end(); ++it) { - return; - } + // electrostatics still has to worry about indirect + // contributions from excluded pairs of atoms: - void InteractionManager::doSkipCorrection(int *atid1, int *atid2, RealType *d, RealType *r, RealType *skippedCharge1, RealType *skippedCharge2, RealType *sw, RealType *electroMult, RealType *pot, RealType *vpair, RealType *f1, RealType *eFrame1, RealType *eFrame2, RealType *t1, RealType *t2){ - - if (!initialized_) initialize(); - - SkipCorrectionData skdat; - - skdat.atype1 = typeMap_[*atid1]; - skdat.atype2 = typeMap_[*atid2]; - skdat.d = Vector3d(d); - skdat.rij = *r; - skdat.skippedCharge1 = *skippedCharge1; - skdat.skippedCharge2 = *skippedCharge2; - skdat.sw = *sw; - skdat.electroMult = *electroMult; - skdat.pot = *pot; - skdat.vpair = *vpair; - skdat.f1 = Vector3d(f1); - skdat.eFrame1 = Mat3x3d(eFrame1); - skdat.eFrame2 = Mat3x3d(eFrame2); - skdat.t1 = Vector3d(t1); - skdat.t2 = Vector3d(t2); - - pair key = make_pair(skdat.atype1, skdat.atype2); - set::iterator it; - - for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it){ - if ((*it)->getFamily() == ELECTROSTATIC_FAMILY) { - dynamic_cast(*it)->calcSkipCorrection(skdat); + if (!idat.excluded || (*it)->getFamily() == ELECTROSTATIC_FAMILY) { + (*it)->calcForce(idat); } } - f1[0] = skdat.f1.x(); - f1[1] = skdat.f1.y(); - f1[2] = skdat.f1.z(); - - t1[0] = skdat.t1.x(); - t1[1] = skdat.t1.y(); - t1[2] = skdat.t1.z(); - - t2[0] = skdat.t2.x(); - t2[1] = skdat.t2.y(); - t2[2] = skdat.t2.z(); - return; } - void InteractionManager::doSelfCorrection(int *atid, RealType *eFrame, RealType *skippedCharge, RealType *pot, RealType *t){ + void InteractionManager::doSelfCorrection(SelfData sdat){ if (!initialized_) initialize(); - SelfCorrectionData scdat; - - scdat.atype = typeMap_[*atid]; - scdat.eFrame = Mat3x3d(eFrame); - scdat.skippedCharge = *skippedCharge; - scdat.pot = *pot; - scdat.t = Vector3d(t); - - pair key = make_pair(scdat.atype, scdat.atype); + 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(scdat); + dynamic_cast(*it)->calcSelfCorrection(sdat); } } - - t[0] = scdat.t.x(); - t[1] = scdat.t.y(); - t[2] = scdat.t.z(); - + return; } - RealType InteractionManager::getSuggestedCutoffRadius(int *atid) { if (!initialized_) initialize(); - + AtomType* atype = typeMap_[*atid]; pair key = make_pair(atype, atype); @@ -416,81 +385,19 @@ namespace OpenMD { RealType cutoff = 0.0; for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it) - cutoff = max(cutoff, (*it)->getSuggestedCutoffRadius(atype, atype)); + cutoff = max(cutoff, (*it)->getSuggestedCutoffRadius(key)); return cutoff; } -} //end namespace OpenMD - -extern "C" { - -#define fortranDoPrePair FC_FUNC(do_prepair, DO_PREPAIR) -#define fortranDoPreForce FC_FUNC(do_preforce, DO_PREFORCE) -#define fortranDoPair FC_FUNC(do_pair, DO_PAIR) -#define fortranDoSkipCorrection FC_FUNC(do_skip_correction, DO_SKIP_CORRECTION) -#define fortranDoSelfCorrection FC_FUNC(do_self_correction, DO_SELF_CORRECTION) -#define fortranGetCutoff FC_FUNC(get_cutoff, GET_CUTOFF) - - void fortranDoPrePair(int *atid1, int *atid2, RealType *rij, - RealType *rho_i_at_j, RealType *rho_j_at_i) { - - return OpenMD::InteractionManager::Instance()->doPrePair(atid1, atid2, rij, - rho_i_at_j, - rho_j_at_i); - } - void fortranDoPreforce(int *atid, RealType *rho, RealType *frho, - RealType *dfrhodrho) { + RealType InteractionManager::getSuggestedCutoffRadius(AtomType* atype) { + if (!initialized_) initialize(); - return OpenMD::InteractionManager::Instance()->doPreForce(atid, rho, frho, - dfrhodrho); - } - - void fortranDoPair(int *atid1, int *atid2, RealType *d, RealType *r, - RealType *r2, RealType *rcut, RealType *sw, - RealType *vdwMult, RealType *electroMult, RealType *pot, - RealType *vpair, RealType *f1, RealType *eFrame1, - RealType *eFrame2, RealType *A1, RealType *A2, - RealType *t1, RealType *t2, RealType *rho1, RealType *rho2, - RealType *dfrho1, RealType *dfrho2, RealType *fshift1, - RealType *fshift2){ + pair key = make_pair(atype, atype); + set::iterator it; + RealType cutoff = 0.0; - return OpenMD::InteractionManager::Instance()->doPair(atid1, atid2, d, r, - r2, rcut, sw, - vdwMult, electroMult, - pot, vpair, f1, - eFrame1, eFrame2, - A1, A2, t1, t2, rho1, - rho2, dfrho1, dfrho2, - fshift1, fshift2); + for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it) + cutoff = max(cutoff, (*it)->getSuggestedCutoffRadius(key)); + return cutoff; } - - void fortranDoSkipCorrection(int *atid1, int *atid2, RealType *d, - RealType *r, RealType *skippedCharge1, - RealType *skippedCharge2, RealType *sw, - RealType *electroMult, RealType *pot, - RealType *vpair, RealType *f1, - RealType *eFrame1, RealType *eFrame2, - RealType *t1, RealType *t2){ - - return OpenMD::InteractionManager::Instance()->doSkipCorrection(atid1, - atid2, d, - r, - skippedCharge1, - skippedCharge2, - sw, electroMult, pot, - vpair, f1, eFrame1, - eFrame2, t1, t2); - } - - void fortranDoSelfCorrection(int *atid, RealType *eFrame, RealType *skippedCharge, - RealType *pot, RealType *t) { - - return OpenMD::InteractionManager::Instance()->doSelfCorrection(atid, - eFrame, - skippedCharge, - pot, t); - } - RealType fortranGetCutoff(int *atid) { - return OpenMD::InteractionManager::Instance()->getSuggestedCutoffRadius(atid); - } -} +} //end namespace OpenMD