--- branches/development/src/nonbonded/InteractionManager.cpp 2010/10/02 19:53:32 1502 +++ branches/development/src/nonbonded/InteractionManager.cpp 2011/05/27 16:45:44 1571 @@ -42,12 +42,30 @@ namespace OpenMD { #include "nonbonded/InteractionManager.hpp" namespace OpenMD { - + + InteractionManager* InteractionManager::_instance = NULL; + SimInfo* InteractionManager::info_ = NULL; bool InteractionManager::initialized_ = false; - ForceField* InteractionManager::forceField_ = NULL; - InteractionManager* InteractionManager::_instance = NULL; + + 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}; + 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) { @@ -58,14 +76,8 @@ namespace OpenMD { void InteractionManager::initialize() { - lj_ = new LJ(); - gb_ = new GB(); - sticky_ = new Sticky(); - eam_ = new EAM(); - sc_ = new SC(); - morse_ = new Morse(); - electrostatic_ = new Electrostatic(); - + ForceField* forceField_ = info_->getForceField(); + lj_->setForceField(forceField_); gb_->setForceField(forceField_); sticky_->setForceField(forceField_); @@ -73,7 +85,28 @@ namespace OpenMD { sc_->setForceField(forceField_); morse_->setForceField(forceField_); 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(); + ForceField::AtomTypeContainer* atomTypes = forceField_->getAtomTypes(); ForceField::AtomTypeContainer::MapTypeIterator i1, i2; AtomType* atype1; @@ -143,80 +176,117 @@ 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); + + 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; } - 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); - } - 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(); - } - // 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: + + 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) { + atype2 = (*jt); key = make_pair(atype1, atype2); if (interactions_[key].size() == 0) { @@ -230,231 +300,260 @@ namespace OpenMD { } } } - } + setupCutoffs(); + setupSwitching(); - void InteractionManager::doPrePair(AtomType* atype1, - AtomType* atype2, - RealType rij, - RealType &rho_i_at_j, - RealType &rho_j_at_i) { - + //int ljsp = cutoffMethod_ == SHIFTED_POTENTIAL ? 1 : 0; + //int ljsf = cutoffMethod_ == SHIFTED_FORCE ? 1 : 0; + //notifyFortranCutoffs(&rCut_, &rSwitch_, &ljsp, &ljsf); + + initialized_ = true; } - void InteractionManager::doPreForce(AtomType* atype, - RealType rho, - RealType &frho, - RealType &dfrhodrho) { + /** + * 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; + } } - void InteractionManager::doSkipCorrection(AtomType* atype1, - AtomType* atype2, - Vector3d d, - RealType rij, - RealType &skippedCharge1, - RealType &skippedCharge2, - RealType sw, - RealType electroMult, - RealType &pot, - RealType &vpair, - Vector3d &f1, - Mat3x3d eFrame1, - Mat3x3d eFrame2, - Vector3d &t1, - Vector3d &t2) { + + /** + * 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(); + } + + 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_); } - - void InteractionManager::doSelfCorrection(AtomType* atype, - Mat3x3d eFrame, - RealType skippedCharge, - RealType &pot, - Vector3d &t) { - } - void InteractionManager::do_prepair(int *atid1, int *atid2, RealType *rij, RealType *rho_i_at_j, RealType *rho_j_at_i){ + void InteractionManager::doPrePair(InteractionData idat){ if (!initialized_) initialize(); - AtomType* atype1 = typeMap_[*atid1]; - AtomType* atype2 = typeMap_[*atid2]; + + 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); + } + } - doPrePair(atype1, atype2, *rij, *rho_i_at_j, *rho_j_at_i); - return; } + + void InteractionManager::doPreForce(SelfData sdat){ - void InteractionManager::do_preforce(int *atid, RealType *rho, RealType *frho, RealType *dfrhodrho){ - - if (!initialized_) initialize(); - AtomType* atype = typeMap_[*atid]; - - doPreForce(atype, *rho, *frho, *dfrhodrho); + if (!initialized_) initialize(); + 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(sdat); + } + } + return; } - void InteractionManager::do_pair(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) + for (it = interactions_[ idat.atypes ].begin(); + it != interactions_[ idat.atypes ].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(); - return; } - void InteractionManager::do_skip_correction(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){ + void InteractionManager::doSkipCorrection(InteractionData idat){ - if (!initialized_) initialize(); - - AtomType* atype1 = typeMap_[*atid1]; - AtomType* atype2 = typeMap_[*atid2]; - Vector3d disp(d); - Vector3d frc(f1); - Vector3d trq1(t1); - Vector3d trq2(t2); - Mat3x3d eFi(eFrame1); - Mat3x3d eFj(eFrame2); - - doSkipCorrection(atype1, atype2, disp, *r, *skippedCharge1, *skippedCharge2, *sw, - *electroMult, *pot, *vpair, frc, eFi, eFj, trq1, trq2); - - f1[0] = frc.x(); - f1[1] = frc.y(); - f1[2] = frc.z(); + if (!initialized_) initialize(); - t1[0] = trq1.x(); - t1[1] = trq1.y(); - t1[2] = trq1.z(); - - t2[0] = trq2.x(); - t2[1] = trq2.y(); - t2[2] = trq2.z(); + set::iterator it; - return; + for (it = interactions_[ idat.atypes ].begin(); + it != interactions_[ idat.atypes ].end(); ++it){ + if ((*it)->getFamily() == ELECTROSTATIC_FAMILY) { + dynamic_cast(*it)->calcSkipCorrection(idat); + } + } + + return; } - void InteractionManager::do_self_correction(int *atid, RealType *eFrame, RealType *skippedCharge, RealType *pot, RealType *t){ + void InteractionManager::doSelfCorrection(SelfData sdat){ if (!initialized_) initialize(); - - AtomType* atype = typeMap_[*atid]; - Mat3x3d eFi(eFrame); - Vector3d trq1(t); - - doSelfCorrection(atype, eFi, *skippedCharge, *pot, trq1); + + pair key = make_pair(sdat.atype, sdat.atype); + set::iterator it; - t[0] = trq1.x(); - t[1] = trq1.y(); - t[2] = trq1.z(); - - return; + for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it){ + if ((*it)->getFamily() == ELECTROSTATIC_FAMILY) { + dynamic_cast(*it)->calcSelfCorrection(sdat); + } + } + + return; } -} //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()->do_prepair(atid1, atid2, rij, - rho_i_at_j, - rho_j_at_i); - } - void fortranDoPreforce(int *atid, RealType *rho, RealType *frho, - RealType *dfrhodrho) { + RealType InteractionManager::getSuggestedCutoffRadius(int *atid) { + if (!initialized_) initialize(); - return OpenMD::InteractionManager::Instance()->do_preforce(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){ - - return OpenMD::InteractionManager::Instance()->do_pair(atid1, atid2, d, r, r2, rcut, - sw, vdwMult, electroMult, pot, - vpair, f1, eFrame1, eFrame2, - A1, A2, t1, t2, rho1, rho2, - dfrho1, dfrho2, fshift1, fshift2); - } + AtomType* atype = typeMap_[*atid]; - 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){ + pair key = make_pair(atype, atype); + set::iterator it; + RealType cutoff = 0.0; - return OpenMD::InteractionManager::Instance()->do_skip_correction(atid1, atid2, d, r, - skippedCharge1, - skippedCharge2, - sw, electroMult, pot, - vpair, f1, eFrame1, - eFrame2, t1, t2); + for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it) + cutoff = max(cutoff, (*it)->getSuggestedCutoffRadius(key)); + return cutoff; } - void fortranDoSelfCorrection(int *atid, RealType *eFrame, RealType *skippedCharge, - RealType *pot, RealType *t) { + RealType InteractionManager::getSuggestedCutoffRadius(AtomType* atype) { + if (!initialized_) initialize(); - return OpenMD::InteractionManager::Instance()->do_self_correction(atid, eFrame, - skippedCharge, - pot, t); + 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)); + return cutoff; } - RealType fortranGetCutoff() { - return OpenMD::InteractionManager::Instance()->getCutoff(); - } -} +} //end namespace OpenMD