--- branches/development/src/nonbonded/Morse.cpp 2010/10/02 19:53:32 1502 +++ trunk/src/nonbonded/Morse.cpp 2014/09/02 18:31:44 2017 @@ -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 @@ -45,78 +46,67 @@ #include #include "nonbonded/Morse.hpp" #include "utils/simError.h" -#include "types/NonBondedInteractionType.hpp" +#include "types/MorseInteractionType.hpp" using namespace std; namespace OpenMD { - Morse::Morse() : name_("Morse"), initialized_(false), forceField_(NULL), - shiftedPot_(false), shiftedFrc_(false) {} + Morse::Morse() : name_("Morse"), initialized_(false), forceField_(NULL) {} void Morse::initialize() { - stringToEnumMap_["shiftedMorse"] = shiftedMorse; - stringToEnumMap_["repulsiveMorse"] = repulsiveMorse; + Mtypes.clear(); + Mtids.clear(); + MixingMap.clear(); + Mtids.resize( forceField_->getNAtomType(), -1); ForceField::NonBondedInteractionTypeContainer* nbiTypes = forceField_->getNonBondedInteractionTypes(); ForceField::NonBondedInteractionTypeContainer::MapTypeIterator j; + ForceField::NonBondedInteractionTypeContainer::KeyType keys; NonBondedInteractionType* nbt; + int mtid1, mtid2; for (nbt = nbiTypes->beginType(j); nbt != NULL; nbt = nbiTypes->nextType(j)) { if (nbt->isMorse()) { - - pair atypes = nbt->getAtomTypes(); - - GenericData* data = nbt->getPropertyByName("Morse"); - if (data == NULL) { - sprintf( painCave.errMsg, "Morse::initialize could not find\n" - "\tMorse parameters for %s - %s interaction.\n", - atypes.first->getName().c_str(), - atypes.second->getName().c_str()); - painCave.severity = OPENMD_ERROR; - painCave.isFatal = 1; - simError(); + keys = nbiTypes->getKeys(j); + AtomType* at1 = forceField_->getAtomType(keys[0]); + AtomType* at2 = forceField_->getAtomType(keys[1]); + + int atid1 = at1->getIdent(); + if (Mtids[atid1] == -1) { + mtid1 = Mtypes.size(); + Mtypes.insert(atid1); + Mtids[atid1] = mtid1; } - - MorseData* morseData = dynamic_cast(data); - if (morseData == NULL) { - sprintf( painCave.errMsg, - "Morse::initialize could not convert GenericData to\n" - "\tMorseData for %s - %s interaction.\n", - atypes.first->getName().c_str(), - atypes.second->getName().c_str()); - painCave.severity = OPENMD_ERROR; - painCave.isFatal = 1; - simError(); + int atid2 = at2->getIdent(); + if (Mtids[atid2] == -1) { + mtid2 = Mtypes.size(); + Mtypes.insert(atid2); + Mtids[atid2] = mtid2; } - - MorseParam morseParam = morseData->getData(); + + MorseInteractionType* mit = dynamic_cast(nbt); - RealType De = morseParam.De; - RealType Re = morseParam.Re; - RealType beta = morseParam.beta; - string interactionType = morseParam.interactionType; - - toUpper(interactionType); - map::iterator i; - i = stringToEnumMap_.find(interactionType); - if (i != stringToEnumMap_.end()) { - addExplicitInteraction(atypes.first, atypes.second, - De, Re, beta, i->second ); - } else { + if (mit == NULL) { sprintf( painCave.errMsg, - "Morse::initialize found unknown Morse interaction type\n" - "\t(%s) for %s - %s interaction.\n", - morseParam.interactionType.c_str(), - atypes.first->getName().c_str(), - atypes.second->getName().c_str()); + "Morse::initialize could not convert NonBondedInteractionType\n" + "\tto MorseInteractionType for %s - %s interaction.\n", + at1->getName().c_str(), + at2->getName().c_str()); painCave.severity = OPENMD_ERROR; painCave.isFatal = 1; simError(); } + + RealType De = mit->getD(); + RealType Re = mit->getR(); + RealType beta = mit->getBeta(); + + MorseType variant = mit->getInteractionType(); + addExplicitInteraction(at1, at2, De, Re, beta, variant ); } } initialized_ = true; @@ -124,110 +114,138 @@ namespace OpenMD { void Morse::addExplicitInteraction(AtomType* atype1, AtomType* atype2, RealType De, RealType Re, RealType beta, - MorseInteractionType mit) { + MorseType mt) { MorseInteractionData mixer; mixer.De = De; mixer.Re = Re; mixer.beta = beta; - mixer.interactionType = mit; + mixer.variant = mt; - pair key1, key2; - key1 = make_pair(atype1, atype2); - key2 = make_pair(atype2, atype1); + int mtid1 = Mtids[atype1->getIdent()]; + int mtid2 = Mtids[atype2->getIdent()]; + int nM = Mtypes.size(); + + MixingMap.resize(nM); + MixingMap[mtid1].resize(nM); - MixingMap[key1] = mixer; - if (key2 != key1) { - MixingMap[key2] = mixer; + MixingMap[mtid1][mtid2] = mixer; + if (mtid2 != mtid1) { + MixingMap[mtid2].resize(nM); + MixingMap[mtid2][mtid1] = mixer; } } - void Morse::calcForce(InteractionData idat) { + void Morse::calcForce(InteractionData &idat) { if (!initialized_) initialize(); + MorseInteractionData &mixer = MixingMap[Mtids[idat.atid1]][Mtids[idat.atid2]]; + RealType myPot = 0.0; RealType myPotC = 0.0; RealType myDeriv = 0.0; RealType myDerivC = 0.0; - - pair key = make_pair(idat.atype1, idat.atype2); - MorseInteractionData mixer = MixingMap[key]; - + RealType De = mixer.De; RealType Re = mixer.Re; RealType beta = mixer.beta; - MorseInteractionType interactionType = mixer.interactionType; - + MorseType variant = mixer.variant; + // V(r) = D_e exp(-a(r-re)(exp(-a(r-re))-2) - RealType expt = -beta*(idat.rij - Re); + RealType expt = -beta*( *(idat.rij) - Re); RealType expfnc = exp(expt); RealType expfnc2 = expfnc*expfnc; - + RealType exptC = 0.0; RealType expfncC = 0.0; RealType expfnc2C = 0.0; - - if (Morse::shiftedPot_ || Morse::shiftedFrc_) { - exptC = -beta*(idat.rcut - Re); + + if (idat.shiftedPot || idat.shiftedForce) { + exptC = -beta*( *(idat.rcut) - Re); expfncC = exp(exptC); expfnc2C = expfncC*expfncC; } - - - switch(interactionType) { - case shiftedMorse : { - + + + switch(variant) { + case mtShifted : { + myPot = De * (expfnc2 - 2.0 * expfnc); myDeriv = 2.0 * De * beta * (expfnc - expfnc2); - - if (Morse::shiftedPot_) { + + if (idat.shiftedPot) { myPotC = De * (expfnc2C - 2.0 * expfncC); myDerivC = 0.0; - } else if (Morse::shiftedFrc_) { + } else if (idat.shiftedForce) { myPotC = De * (expfnc2C - 2.0 * expfncC); - myDerivC = 2.0 * De * beta * (expfnc2C - expfnc2C); - myPotC += myDerivC * (idat.rij - idat.rcut); + myDerivC = 2.0 * De * beta * (expfncC - expfnc2C); + myPotC += myDerivC * ( *(idat.rij) - *(idat.rcut) ); } else { myPotC = 0.0; myDerivC = 0.0; } - + break; } - case repulsiveMorse : { - + case mtRepulsive : { + myPot = De * expfnc2; myDeriv = -2.0 * De * beta * expfnc2; - - if (Morse::shiftedPot_) { + + if (idat.shiftedPot) { myPotC = De * expfnc2C; myDerivC = 0.0; - } else if (Morse::shiftedFrc_) { + } else if (idat.shiftedForce) { myPotC = De * expfnc2C; myDerivC = -2.0 * De * beta * expfnc2C; - myPotC += myDerivC * (idat.rij - idat.rcut); + myPotC += myDerivC * ( *(idat.rij) - *(idat.rcut)); } else { myPotC = 0.0; myDerivC = 0.0; } - + break; } + case mtUnknown: { + // don't know what to do so don't do anything + break; } + } + - RealType pot_temp = idat.vdwMult * (myPot - myPotC); - idat.vpair += pot_temp; - - RealType dudr = idat.sw * idat.vdwMult * (myDeriv - myDerivC); + RealType pot_temp = *(idat.vdwMult) * (myPot - myPotC); + *(idat.vpair) += pot_temp; - idat.pot += idat.sw * pot_temp; - idat.f1 = idat.d * dudr / idat.rij; + RealType dudr = *(idat.sw) * *(idat.vdwMult) * (myDeriv - myDerivC); - return; - + + (*(idat.pot))[VANDERWAALS_FAMILY] += *(idat.sw) * pot_temp; + *(idat.f1) = *(idat.d) * dudr / *(idat.rij); + + return; } + + RealType Morse::getSuggestedCutoffRadius(pair atypes) { + if (!initialized_) initialize(); + int atid1 = atypes.first->getIdent(); + int atid2 = atypes.second->getIdent(); + int mtid1 = Mtids[atid1]; + int mtid2 = Mtids[atid2]; + + if ( mtid1 == -1 || mtid2 == -1) return 0.0; + else { + MorseInteractionData mixer = MixingMap[mtid1][mtid2]; + RealType Re = mixer.Re; + RealType beta = mixer.beta; + // This value of the r corresponds to an energy about 1.48% of + // the energy at the bottom of the Morse well. For comparison, the + // Lennard-Jones function is about 1.63% of it's minimum value at + // a distance of 2.5 sigma. + return (4.9 + beta * Re) / beta; + } + } }