--- branches/development/src/nonbonded/Morse.cpp 2011/05/27 16:45:44 1571 +++ 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,21 +114,25 @@ 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; } } @@ -146,100 +140,104 @@ namespace OpenMD { if (!initialized_) initialize(); - map, MorseInteractionData>::iterator it; - it = MixingMap.find( idat.atypes ); - if (it != MixingMap.end()) { - MorseInteractionData mixer = (*it).second; + 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; + + RealType De = mixer.De; + RealType Re = mixer.Re; + RealType beta = mixer.beta; + MorseType variant = mixer.variant; + + // V(r) = D_e exp(-a(r-re)(exp(-a(r-re))-2) + + 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 (idat.shiftedPot || idat.shiftedForce) { + exptC = -beta*( *(idat.rcut) - Re); + expfncC = exp(exptC); + expfnc2C = expfncC*expfncC; + } + + + switch(variant) { + case mtShifted : { - RealType myPot = 0.0; - RealType myPotC = 0.0; - RealType myDeriv = 0.0; - RealType myDerivC = 0.0; + myPot = De * (expfnc2 - 2.0 * expfnc); + myDeriv = 2.0 * De * beta * (expfnc - expfnc2); - RealType De = mixer.De; - RealType Re = mixer.Re; - RealType beta = mixer.beta; - MorseInteractionType interactionType = mixer.interactionType; - - // V(r) = D_e exp(-a(r-re)(exp(-a(r-re))-2) - - 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); - expfncC = exp(exptC); - expfnc2C = expfncC*expfncC; + if (idat.shiftedPot) { + myPotC = De * (expfnc2C - 2.0 * expfncC); + myDerivC = 0.0; + } else if (idat.shiftedForce) { + myPotC = De * (expfnc2C - 2.0 * expfncC); + myDerivC = 2.0 * De * beta * (expfncC - expfnc2C); + myPotC += myDerivC * ( *(idat.rij) - *(idat.rcut) ); + } else { + myPotC = 0.0; + myDerivC = 0.0; } - - switch(interactionType) { - case shiftedMorse : { + break; + } + case mtRepulsive : { - myPot = De * (expfnc2 - 2.0 * expfnc); - myDeriv = 2.0 * De * beta * (expfnc - expfnc2); - - if (Morse::shiftedPot_) { - myPotC = De * (expfnc2C - 2.0 * expfncC); - myDerivC = 0.0; - } else if (Morse::shiftedFrc_) { - myPotC = De * (expfnc2C - 2.0 * expfncC); - myDerivC = 2.0 * De * beta * (expfnc2C - expfnc2C); - myPotC += myDerivC * ( *(idat.rij) - *(idat.rcut) ); - } else { - myPotC = 0.0; - myDerivC = 0.0; - } - - break; + myPot = De * expfnc2; + myDeriv = -2.0 * De * beta * expfnc2; + + if (idat.shiftedPot) { + myPotC = De * expfnc2C; + myDerivC = 0.0; + } else if (idat.shiftedForce) { + myPotC = De * expfnc2C; + myDerivC = -2.0 * De * beta * expfnc2C; + myPotC += myDerivC * ( *(idat.rij) - *(idat.rcut)); + } else { + myPotC = 0.0; + myDerivC = 0.0; } - case repulsiveMorse : { - - myPot = De * expfnc2; - myDeriv = -2.0 * De * beta * expfnc2; - - if (Morse::shiftedPot_) { - myPotC = De * expfnc2C; - myDerivC = 0.0; - } else if (Morse::shiftedFrc_) { - myPotC = De * expfnc2C; - myDerivC = -2.0 * De * beta * expfnc2C; - myPotC += myDerivC * ( *(idat.rij) - *(idat.rcut)); - } else { - myPotC = 0.0; - myDerivC = 0.0; - } - - break; - } - } - RealType pot_temp = *(idat.vdwMult) * (myPot - myPotC); - *(idat.vpair) += pot_temp; - - RealType dudr = *(idat.sw) * *(idat.vdwMult) * (myDeriv - myDerivC); - - idat.pot[VANDERWAALS_FAMILY] += *(idat.sw) * pot_temp; - *(idat.f1) = *(idat.d) * dudr / *(idat.rij); + break; } - return; + 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); + + + (*(idat.pot))[VANDERWAALS_FAMILY] += *(idat.sw) * pot_temp; + *(idat.f1) = *(idat.d) * dudr / *(idat.rij); + + return; + } + RealType Morse::getSuggestedCutoffRadius(pair atypes) { if (!initialized_) initialize(); - map, MorseInteractionData>::iterator it; - it = MixingMap.find(atypes); - if (it == MixingMap.end()) - return 0.0; - else { - MorseInteractionData mixer = (*it).second; - + + 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