--- branches/development/src/nonbonded/RepulsivePower.cpp 2011/09/13 14:12:54 1624 +++ trunk/src/nonbonded/RepulsivePower.cpp 2014/11/01 14:12:16 2033 @@ -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,7 +46,7 @@ #include #include "nonbonded/RepulsivePower.hpp" #include "utils/simError.h" -#include "types/NonBondedInteractionType.hpp" +#include "types/RepulsivePowerInteractionType.hpp" using namespace std; @@ -56,48 +57,55 @@ namespace OpenMD { void RepulsivePower::initialize() { + RPtypes.clear(); + RPtids.clear(); + MixingMap.clear(); + RPtids.resize( forceField_->getNAtomType(), -1); + ForceField::NonBondedInteractionTypeContainer* nbiTypes = forceField_->getNonBondedInteractionTypes(); ForceField::NonBondedInteractionTypeContainer::MapTypeIterator j; + ForceField::NonBondedInteractionTypeContainer::KeyType keys; NonBondedInteractionType* nbt; + int rptid1, rptid2; for (nbt = nbiTypes->beginType(j); nbt != NULL; nbt = nbiTypes->nextType(j)) { if (nbt->isRepulsivePower()) { - - pair atypes = nbt->getAtomTypes(); - - GenericData* data = nbt->getPropertyByName("RepulsivePower"); - if (data == NULL) { - sprintf( painCave.errMsg, "RepulsivePower::initialize could not\n" - "\tfind RepulsivePower 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 (RPtids[atid1] == -1) { + rptid1 = RPtypes.size(); + RPtypes.insert(atid1); + RPtids[atid1] = rptid1; } - - RepulsivePowerData* rpData = dynamic_cast(data); - if (rpData == NULL) { + int atid2 = at2->getIdent(); + if (RPtids[atid2] == -1) { + rptid2 = RPtypes.size(); + RPtypes.insert(atid2); + RPtids[atid2] = rptid2; + } + + RepulsivePowerInteractionType* rpit = dynamic_cast(nbt); + if (rpit == NULL) { sprintf( painCave.errMsg, - "RepulsivePower::initialize could not convert GenericData\n" - "\tto RepulsivePowerData for %s - %s interaction.\n", - atypes.first->getName().c_str(), - atypes.second->getName().c_str()); + "RepulsivePower::initialize could not convert NonBondedInteractionType\n" + "\tto RepulsivePowerInteractionType for %s - %s interaction.\n", + at1->getName().c_str(), + at2->getName().c_str()); painCave.severity = OPENMD_ERROR; painCave.isFatal = 1; simError(); } - RepulsivePowerParam rpParam = rpData->getData(); + RealType sigma = rpit->getSigma(); + RealType epsilon = rpit->getEpsilon(); + int nRep = rpit->getNrep(); - RealType sigma = rpParam.sigma; - RealType epsilon = rpParam.epsilon; - int nRep = rpParam.nRep; - - addExplicitInteraction(atypes.first, atypes.second, - sigma, epsilon, nRep); + addExplicitInteraction(at1, at2, sigma, epsilon, nRep); } } initialized_ = true; @@ -115,13 +123,17 @@ namespace OpenMD { mixer.sigmai = 1.0 / mixer.sigma; mixer.nRep = nRep; - pair key1, key2; - key1 = make_pair(atype1, atype2); - key2 = make_pair(atype2, atype1); + int rptid1 = RPtids[atype1->getIdent()]; + int rptid2 = RPtids[atype2->getIdent()]; + int nRP = RPtypes.size(); + + MixingMap.resize(nRP); + MixingMap[rptid1].resize(nRP); - MixingMap[key1] = mixer; - if (key2 != key1) { - MixingMap[key2] = mixer; + MixingMap[rptid1][rptid2] = mixer; + if (rptid2 != rptid1) { + MixingMap[rptid2].resize(nRP); + MixingMap[rptid2][rptid1] = mixer; } } @@ -129,53 +141,49 @@ namespace OpenMD { if (!initialized_) initialize(); - map, RPInteractionData>::iterator it; - it = MixingMap.find( idat.atypes ); - - if (it != MixingMap.end()) { - - RPInteractionData mixer = (*it).second; - RealType sigmai = mixer.sigmai; - RealType epsilon = mixer.epsilon; - int nRep = mixer.nRep; - - RealType ros; - RealType rcos; - RealType myPot = 0.0; - RealType myPotC = 0.0; - RealType myDeriv = 0.0; - RealType myDerivC = 0.0; - - ros = *(idat.rij) * sigmai; - - getNRepulsionFunc(ros, nRep, myPot, myDeriv); - - if (idat.shiftedPot) { - rcos = *(idat.rcut) * sigmai; - getNRepulsionFunc(rcos, nRep, myPotC, myDerivC); - myDerivC = 0.0; - } else if (idat.shiftedForce) { - rcos = *(idat.rcut) * sigmai; - getNRepulsionFunc(rcos, nRep, myPotC, myDerivC); - myPotC = myPotC + myDerivC * (*(idat.rij) - *(idat.rcut)) * sigmai; - } else { - myPotC = 0.0; - myDerivC = 0.0; - } - - RealType pot_temp = *(idat.vdwMult) * epsilon * (myPot - myPotC); - *(idat.vpair) += pot_temp; - - RealType dudr = *(idat.sw) * *(idat.vdwMult) * epsilon * (myDeriv - - myDerivC)*sigmai; - - (*(idat.pot))[VANDERWAALS_FAMILY] += *(idat.sw) * pot_temp; - *(idat.f1) = *(idat.d) * dudr / *(idat.rij); + RPInteractionData &mixer = MixingMap[RPtids[idat.atid1]][RPtids[idat.atid2]]; + RealType sigmai = mixer.sigmai; + RealType epsilon = mixer.epsilon; + int nRep = mixer.nRep; + + RealType ros; + RealType rcos; + RealType myPot = 0.0; + RealType myPotC = 0.0; + RealType myDeriv = 0.0; + RealType myDerivC = 0.0; + + ros = *(idat.rij) * sigmai; + + getNRepulsionFunc(ros, nRep, myPot, myDeriv); + + if (idat.shiftedPot) { + rcos = *(idat.rcut) * sigmai; + getNRepulsionFunc(rcos, nRep, myPotC, myDerivC); + myDerivC = 0.0; + } else if (idat.shiftedForce) { + rcos = *(idat.rcut) * sigmai; + getNRepulsionFunc(rcos, nRep, myPotC, myDerivC); + myPotC = myPotC + myDerivC * (*(idat.rij) - *(idat.rcut)) * sigmai; + } else { + myPotC = 0.0; + myDerivC = 0.0; } + + RealType pot_temp = *(idat.vdwMult) * epsilon * (myPot - myPotC); + *(idat.vpair) += pot_temp; + + RealType dudr = *(idat.sw) * *(idat.vdwMult) * epsilon * (myDeriv - + myDerivC)*sigmai; + + (*(idat.pot))[VANDERWAALS_FAMILY] += *(idat.sw) * pot_temp; + *(idat.f1) = *(idat.d) * dudr / *(idat.rij); + return; } - void RepulsivePower::getNRepulsionFunc(RealType r, int n, RealType &pot, RealType &deriv) { + void RepulsivePower::getNRepulsionFunc(const RealType &r, int &n, + RealType &pot, RealType &deriv) { RealType ri = 1.0 / r; RealType rin = pow(ri, n); @@ -190,15 +198,17 @@ namespace OpenMD { RealType RepulsivePower::getSuggestedCutoffRadius(pair atypes) { if (!initialized_) initialize(); - map, RPInteractionData>::iterator it; - it = MixingMap.find(atypes); - if (it == MixingMap.end()) - return 0.0; - else { - RPInteractionData mixer = (*it).second; + + int atid1 = atypes.first->getIdent(); + int atid2 = atypes.second->getIdent(); + int rptid1 = RPtids[atid1]; + int rptid2 = RPtids[atid2]; + + if (rptid1 == -1 || rptid2 == -1) return 0.0; + else { + RPInteractionData mixer = MixingMap[rptid1][rptid2]; return 2.5 * mixer.sigma; } - } - + } }