--- branches/development/src/nonbonded/Morse.cpp 2010/09/15 19:32:10 1501 +++ 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,92 +46,67 @@ #include #include "nonbonded/Morse.hpp" #include "utils/simError.h" -#include "types/NonBondedInteractionType.hpp" +#include "types/MorseInteractionType.hpp" using namespace std; namespace OpenMD { - bool Morse::initialized_ = false; - bool Morse::shiftedPot_ = false; - bool Morse::shiftedFrc_ = false; - ForceField* Morse::forceField_ = NULL; - map Morse::MorseMap; - map, MorseInteractionData> Morse::MixingMap; - map Morse::stringToEnumMap_; + Morse::Morse() : name_("Morse"), initialized_(false), forceField_(NULL) {} - Morse* Morse::_instance = NULL; - - Morse* Morse::Instance() { - if (!_instance) { - _instance = new Morse(); - } - return _instance; - } - 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; @@ -138,165 +114,138 @@ namespace OpenMD { void Morse::addExplicitInteraction(AtomType* atype1, AtomType* atype2, RealType De, RealType Re, RealType beta, - MorseInteractionType mit) { + MorseType mt) { - AtomTypeProperties atp = atype1->getATP(); - MorseMap.insert( pair(atp.ident, atype1) ); - - atp = atype2->getATP(); - MorseMap.insert( pair(atp.ident, atype2) ); - 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(AtomType* at1, AtomType* at2, Vector3d d, - RealType r, RealType r2, RealType rcut, RealType sw, - RealType vdwMult, RealType &vpair, RealType &pot, - Vector3d &f1) { + 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(at1, at2); - 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*(r - 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*(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 * (r - 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 * (r - 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 = vdwMult * (myPot - myPotC); - vpair += pot_temp; - - RealType dudr = sw * vdwMult * (myDeriv - myDerivC); + RealType pot_temp = *(idat.vdwMult) * (myPot - myPotC); + *(idat.vpair) += pot_temp; - pot += sw * pot_temp; - f1 = d * dudr / r; + RealType dudr = *(idat.sw) * *(idat.vdwMult) * (myDeriv - myDerivC); - return; - - } - - void Morse::do_morse_pair(int *atid1, int *atid2, RealType *d, - RealType *rij, RealType *r2, RealType *rcut, - RealType *sw, RealType *vdwMult, - RealType *vpair, RealType *pot, RealType *f1) { - - if (!initialized_) initialize(); - AtomType* atype1 = MorseMap[*atid1]; - AtomType* atype2 = MorseMap[*atid2]; + (*(idat.pot))[VANDERWAALS_FAMILY] += *(idat.sw) * pot_temp; + *(idat.f1) = *(idat.d) * dudr / *(idat.rij); - Vector3d disp(d[0], d[1], d[2]); - Vector3d frc(f1[0], f1[1], f1[2]); - - calcForce(atype1, atype2, disp, *rij, *r2, *rcut, *sw, *vdwMult, *vpair, - *pot, frc); - - f1[0] = frc.x(); - f1[1] = frc.y(); - f1[2] = frc.z(); - return; } - - void Morse::setMorseDefaultCutoff(RealType *thisRcut, int *sP, int *sF) { - shiftedPot_ = (bool)(*sP); - shiftedFrc_ = (bool)(*sF); - } -} -extern "C" { - -#define fortranSetMorseCutoff FC_FUNC(setmorsedefaultcutoff, SETMORSEDEFAULTCUTOFF) -#define fortranDoMorsePair FC_FUNC(do_morse_pair, DO_MORSE_PAIR) - - void fortranSetMorseCutoff(RealType *rcut, int *shiftedPot, int *shiftedFrc) { - return OpenMD::Morse::Instance()->setMorseDefaultCutoff(rcut, shiftedPot, - shiftedFrc); - } - void fortranDoMorsePair(int *atid1, int *atid2, RealType *d, RealType *rij, - RealType *r2, RealType *rcut, RealType *sw, - RealType *vdwMult, RealType* vpair, RealType* pot, - RealType *f1){ + RealType Morse::getSuggestedCutoffRadius(pair atypes) { + if (!initialized_) initialize(); - return OpenMD::Morse::Instance()->do_morse_pair(atid1, atid2, d, rij, r2, - rcut, sw, vdwMult, vpair, - pot, f1); + 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; + } } } +