--- branches/development/src/nonbonded/LJ.cpp 2011/04/27 18:38:15 1549 +++ branches/development/src/nonbonded/LJ.cpp 2013/02/20 15:39:39 1850 @@ -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,56 +46,19 @@ #include #include "nonbonded/LJ.hpp" #include "utils/simError.h" +#include "types/LennardJonesAdapter.hpp" +#include "types/LennardJonesInteractionType.hpp" namespace OpenMD { - LJ::LJ() : name_("LJ"), initialized_(false), shiftedPot_(false), - shiftedFrc_(false), forceField_(NULL) {} + LJ::LJ() : name_("LJ"), initialized_(false), forceField_(NULL) {} - LJParam LJ::getLJParam(AtomType* atomType) { - - // Do sanity checking on the AtomType we were passed before - // building any data structures: - if (!atomType->isLennardJones()) { - sprintf( painCave.errMsg, - "LJ::getLJParam was passed an atomType (%s) that does not\n" - "\tappear to be a Lennard-Jones atom.\n", - atomType->getName().c_str()); - painCave.severity = OPENMD_ERROR; - painCave.isFatal = 1; - simError(); - } - - GenericData* data = atomType->getPropertyByName("LennardJones"); - if (data == NULL) { - sprintf( painCave.errMsg, "LJ::getLJParam could not find Lennard-Jones\n" - "\tparameters for atomType %s.\n", atomType->getName().c_str()); - painCave.severity = OPENMD_ERROR; - painCave.isFatal = 1; - simError(); - } - - LJParamGenericData* ljData = dynamic_cast(data); - if (ljData == NULL) { - sprintf( painCave.errMsg, - "LJ::getLJParam could not convert GenericData to LJParam for\n" - "\tatom type %s\n", atomType->getName().c_str()); - painCave.severity = OPENMD_ERROR; - painCave.isFatal = 1; - simError(); - } - - return ljData->getData(); - } + RealType LJ::getSigma(AtomType* atomType1, AtomType* atomType2) { - RealType LJ::getSigma(AtomType* atomType) { - LJParam ljParam = getLJParam(atomType); - return ljParam.sigma; - } - - RealType LJ::getSigma(AtomType* atomType1, AtomType* atomType2) { - RealType sigma1 = getSigma(atomType1); - RealType sigma2 = getSigma(atomType2); + LennardJonesAdapter lja1 = LennardJonesAdapter(atomType1); + LennardJonesAdapter lja2 = LennardJonesAdapter(atomType2); + RealType sigma1 = lja1.getSigma(); + RealType sigma2 = lja2.getSigma(); ForceFieldOptions& fopts = forceField_->getForceFieldOptions(); string DistanceMix = fopts.getDistanceMixingRule(); @@ -106,14 +70,12 @@ namespace OpenMD { return 0.5 * (sigma1 + sigma2); } - RealType LJ::getEpsilon(AtomType* atomType) { - LJParam ljParam = getLJParam(atomType); - return ljParam.epsilon; - } - - RealType LJ::getEpsilon(AtomType* atomType1, AtomType* atomType2) { - RealType epsilon1 = getEpsilon(atomType1); - RealType epsilon2 = getEpsilon(atomType2); + RealType LJ::getEpsilon(AtomType* atomType1, AtomType* atomType2) { + LennardJonesAdapter lja1 = LennardJonesAdapter(atomType1); + LennardJonesAdapter lja2 = LennardJonesAdapter(atomType2); + + RealType epsilon1 = lja1.getEpsilon(); + RealType epsilon2 = lja2.getEpsilon(); return sqrt(epsilon1 * epsilon2); } @@ -124,51 +86,41 @@ namespace OpenMD { for (at = atomTypes->beginType(i); at != NULL; at = atomTypes->nextType(i)) { - - if (at->isLennardJones()) - addType(at); + LennardJonesAdapter lja = LennardJonesAdapter(at); + if (lja.isLennardJones()){ + addType(at); + } } - ForceField::NonBondedInteractionTypeContainer* nbiTypes = forceField_->getNonBondedInteractionTypes(); ForceField::NonBondedInteractionTypeContainer::MapTypeIterator j; NonBondedInteractionType* nbt; + ForceField::NonBondedInteractionTypeContainer::KeyType keys; + for (nbt = nbiTypes->beginType(j); nbt != NULL; nbt = nbiTypes->nextType(j)) { - + if (nbt->isLennardJones()) { - - pair atypes = nbt->getAtomTypes(); - - GenericData* data = nbt->getPropertyByName("LennardJones"); - if (data == NULL) { - sprintf( painCave.errMsg, "LJ::rebuildMixingMap could not find\n" - "\tLennard-Jones parameters for %s - %s interaction.\n", - atypes.first->getName().c_str(), - atypes.second->getName().c_str()); - painCave.severity = OPENMD_ERROR; - painCave.isFatal = 1; - simError(); - } - - LJParamGenericData* ljData = dynamic_cast(data); - if (ljData == NULL) { + keys = nbiTypes->getKeys(j); + AtomType* at1 = forceField_->getAtomType(keys[0]); + AtomType* at2 = forceField_->getAtomType(keys[1]); + + LennardJonesInteractionType* ljit = dynamic_cast(nbt); + + if (ljit == NULL) { sprintf( painCave.errMsg, - "LJ::rebuildMixingMap could not convert GenericData to\n" - "\tLJParam for %s - %s interaction.\n", - atypes.first->getName().c_str(), - atypes.second->getName().c_str()); + "LJ::initialize could not convert NonBondedInteractionType\n" + "\tto LennardJonesInteractionType for %s - %s interaction.\n", + at1->getName().c_str(), + at2->getName().c_str()); painCave.severity = OPENMD_ERROR; painCave.isFatal = 1; simError(); } - - LJParam ljParam = ljData->getData(); - RealType sigma = ljParam.sigma; - RealType epsilon = ljParam.epsilon; - - addExplicitInteraction(atypes.first, atypes.second, sigma, epsilon); + RealType sigma = ljit->getSigma(); + RealType epsilon = ljit->getEpsilon(); + addExplicitInteraction(at1, at2, sigma, epsilon); } } initialized_ = true; @@ -177,18 +129,19 @@ namespace OpenMD { void LJ::addType(AtomType* atomType){ - RealType sigma1 = getSigma(atomType); - RealType epsilon1 = getEpsilon(atomType); + LennardJonesAdapter lja1 = LennardJonesAdapter(atomType); + + RealType sigma1 = lja1.getSigma(); + RealType epsilon1 = lja1.getEpsilon(); // add it to the map: - AtomTypeProperties atp = atomType->getATP(); pair::iterator,bool> ret; - ret = LJMap.insert( pair(atp.ident, atomType) ); + ret = LJMap.insert( pair(atomType->getIdent(), atomType) ); if (ret.second == false) { sprintf( painCave.errMsg, "LJ already had a previous entry with ident %d\n", - atp.ident); + atomType->getIdent()); painCave.severity = OPENMD_INFO; painCave.isFatal = 0; simError(); @@ -220,10 +173,6 @@ namespace OpenMD { void LJ::addExplicitInteraction(AtomType* atype1, AtomType* atype2, RealType sigma, RealType epsilon){ - // in case these weren't already in the map - addType(atype1); - addType(atype2); - LJInteractionData mixer; mixer.sigma = sigma; mixer.epsilon = epsilon; @@ -241,11 +190,9 @@ namespace OpenMD { } void LJ::calcForce(InteractionData &idat) { - if (!initialized_) initialize(); - map, LJInteractionData>::iterator it; - it = MixingMap.find(idat.atypes); + it = MixingMap.find( idat.atypes ); if (it != MixingMap.end()) { @@ -261,31 +208,31 @@ namespace OpenMD { RealType myDeriv = 0.0; RealType myDerivC = 0.0; - ros = idat.rij * sigmai; + ros = *(idat.rij) * sigmai; getLJfunc(ros, myPot, myDeriv); - if (shiftedPot_) { - rcos = idat.rcut * sigmai; + if (idat.shiftedPot) { + rcos = *(idat.rcut) * sigmai; getLJfunc(rcos, myPotC, myDerivC); myDerivC = 0.0; - } else if (LJ::shiftedFrc_) { - rcos = idat.rcut * sigmai; + } else if (idat.shiftedForce) { + rcos = *(idat.rcut) * sigmai; getLJfunc(rcos, myPotC, myDerivC); - myPotC = myPotC + myDerivC * (idat.rij - idat.rcut) * sigmai; + 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 pot_temp = *(idat.vdwMult) * epsilon * (myPot - myPotC); + *(idat.vpair) += pot_temp; - RealType dudr = idat.sw * idat.vdwMult * epsilon * (myDeriv - - myDerivC)*sigmai; - - idat.pot[0] += idat.sw * pot_temp; - idat.f1 = idat.d * dudr / idat.rij; + 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; }