--- branches/development/src/nonbonded/LJ.cpp 2010/07/21 18:31:05 1476 +++ trunk/src/nonbonded/LJ.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,90 +46,22 @@ #include #include "nonbonded/LJ.hpp" #include "utils/simError.h" +#include "types/LennardJonesAdapter.hpp" +#include "types/LennardJonesInteractionType.hpp" - namespace OpenMD { - bool LJ::initialized_ = false; - bool LJ::shiftedPot_ = false; - bool LJ::shiftedFrc_ = false; - ForceField* LJ::forceField_ = NULL; - std::map LJ::LJMap; - std::map, LJInteractionData> LJ::MixingMap; - - LJ* LJ::_instance = NULL; + LJ::LJ() : name_("LJ"), initialized_(false), forceField_(NULL) {} - LJ* LJ::Instance() { - if (!_instance) { - _instance = new LJ(); - } - return _instance; - } + RealType LJ::getSigma(AtomType* atomType1, AtomType* atomType2) { - LJParam LJ::getLJParam(AtomType* atomType) { + LennardJonesAdapter lja1 = LennardJonesAdapter(atomType1); + LennardJonesAdapter lja2 = LennardJonesAdapter(atomType2); + RealType sigma1 = lja1.getSigma(); + RealType sigma2 = lja2.getSigma(); - // 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* atomType) { - LJParam ljParam = getLJParam(atomType); - return ljParam.sigma; - } - - RealType LJ::getSigma(int atid) { - if (!initialized_) initialize(); - std::map :: const_iterator it; - it = LJMap.find(atid); - if (it == LJMap.end()) { - sprintf( painCave.errMsg, - "LJ::getSigma could not find atid %d in LJMap\n", - (atid)); - painCave.severity = OPENMD_ERROR; - painCave.isFatal = 1; - simError(); - } - AtomType* atype = it->second; - - return getSigma(atype); - } - - RealType LJ::getSigma(AtomType* atomType1, AtomType* atomType2) { - RealType sigma1 = getSigma(atomType1); - RealType sigma2 = getSigma(atomType2); - ForceFieldOptions& fopts = forceField_->getForceFieldOptions(); - std::string DistanceMix = fopts.getDistanceMixingRule(); + string DistanceMix = fopts.getDistanceMixingRule(); toUpper(DistanceMix); if (DistanceMix == "GEOMETRIC") @@ -137,228 +70,203 @@ 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) { + LennardJonesAdapter lja1 = LennardJonesAdapter(atomType1); + LennardJonesAdapter lja2 = LennardJonesAdapter(atomType2); + + RealType epsilon1 = lja1.getEpsilon(); + RealType epsilon2 = lja2.getEpsilon(); + return sqrt(epsilon1 * epsilon2); } - RealType LJ::getEpsilon(int atid) { - if (!initialized_) initialize(); - std::map :: const_iterator it; - it = LJMap.find(atid); - if (it == LJMap.end()) { - sprintf( painCave.errMsg, - "LJ::getEpsilon could not find atid %d in LJMap\n", - (atid)); - painCave.severity = OPENMD_ERROR; - painCave.isFatal = 1; - simError(); - } - AtomType* atype = it->second; + void LJ::initialize() { - return getEpsilon(atype); - } + LJtypes.clear(); + LJtids.clear(); + MixingMap.clear(); + nLJ_ = 0; - RealType LJ::getEpsilon(AtomType* atomType1, AtomType* atomType2) { - RealType epsilon1 = getEpsilon(atomType1); - RealType epsilon2 = getEpsilon(atomType2); - return sqrt(epsilon1 * epsilon2); - } + LJtids.resize( forceField_->getNAtomType(), -1); - void LJ::initialize() { - ForceField::AtomTypeContainer* atomTypes = forceField_->getAtomTypes(); - ForceField::AtomTypeContainer::MapTypeIterator i; - AtomType* at; + set::iterator at; + for (at = simTypes_.begin(); at != simTypes_.end(); ++at) { + if ((*at)->isLennardJones()) nLJ_++; + } - for (at = atomTypes->beginType(i); at != NULL; - at = atomTypes->nextType(i)) { - - if (at->isLennardJones()) - addType(at); + MixingMap.resize(nLJ_); + + for (at = simTypes_.begin(); at != simTypes_.end(); ++at) { + if ((*at)->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()) { - - std::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; } - 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(); - - std::pair::iterator,bool> ret; - ret = LJMap.insert( std::pair(atp.ident, atomType) ); + int atid = atomType->getIdent(); + int ljtid = LJtypes.size(); + + pair::iterator,bool> ret; + ret = LJtypes.insert( atid ); if (ret.second == false) { sprintf( painCave.errMsg, "LJ already had a previous entry with ident %d\n", - atp.ident); + atid) ; painCave.severity = OPENMD_INFO; painCave.isFatal = 0; simError(); } - + + LJtids[atid] = ljtid; + MixingMap[ljtid].resize( nLJ_ ); + // Now, iterate over all known types and add to the mixing map: - std::map::iterator it; - for( it = LJMap.begin(); it != LJMap.end(); ++it) { - - AtomType* atype2 = (*it).second; + std::set::iterator it; + for( it = LJtypes.begin(); it != LJtypes.end(); ++it) { + int ljtid2 = LJtids[ (*it) ]; + AtomType* atype2 = forceField_->getAtomType( (*it) ); + LJInteractionData mixer; mixer.sigma = getSigma(atomType, atype2); mixer.epsilon = getEpsilon(atomType, atype2); mixer.sigmai = 1.0 / mixer.sigma; mixer.explicitlySet = false; + MixingMap[ljtid2].resize( nLJ_ ); - std::pair key1, key2; - key1 = std::make_pair(atomType, atype2); - key2 = std::make_pair(atype2, atomType); - - MixingMap[key1] = mixer; - if (key2 != key1) { - MixingMap[key2] = mixer; + MixingMap[ljtid][ljtid2] = mixer; + if (ljtid2 != ljtid) { + MixingMap[ljtid2][ljtid] = mixer; } } } 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; mixer.sigmai = 1.0 / mixer.sigma; mixer.explicitlySet = true; - std::pair key1, key2; - key1 = std::make_pair(atype1, atype2); - key2 = std::make_pair(atype2, atype1); + int atid1 = atype1->getIdent(); + int atid2 = atype2->getIdent(); + + int ljtid1, ljtid2; + + pair::iterator,bool> ret; + ret = LJtypes.insert( atid1 ); + if (ret.second == false) { + // already had this type in the LJMap, just get the ljtid: + ljtid1 = LJtids[ atid1 ]; + } else { + // didn't already have it, so make a new one and assign it: + ljtid1 = nLJ_; + LJtids[atid1] = nLJ_; + nLJ_++; + } + + ret = LJtypes.insert( atid2 ); + if (ret.second == false) { + // already had this type in the LJMap, just get the ljtid: + ljtid2 = LJtids[ atid2 ]; + } else { + // didn't already have it, so make a new one and assign it: + ljtid2 = nLJ_; + LJtids[atid2] = nLJ_; + nLJ_++; + } - MixingMap[key1] = mixer; - if (key2 != key1) { - MixingMap[key2] = mixer; + MixingMap.resize(nLJ_); + MixingMap[ljtid1].resize(nLJ_); + + MixingMap[ljtid1][ljtid2] = mixer; + if (ljtid2 != ljtid1) { + MixingMap[ljtid2].resize(nLJ_); + MixingMap[ljtid2][ljtid1] = mixer; } } - void LJ::calcForce(AtomType* at1, AtomType* at2, Vector3d d, - RealType rij, RealType r2, RealType rcut, RealType sw, - RealType vdwMult, RealType &vpair, RealType &pot, - Vector3d &f1) { - + void LJ::calcForce(InteractionData &idat) { if (!initialized_) initialize(); + LJInteractionData &mixer = MixingMap[LJtids[idat.atid1]][LJtids[idat.atid2]]; + + RealType sigmai = mixer.sigmai; + RealType epsilon = mixer.epsilon; + RealType ros; RealType rcos; RealType myPot = 0.0; RealType myPotC = 0.0; RealType myDeriv = 0.0; RealType myDerivC = 0.0; - - std::pair key = std::make_pair(at1, at2); - LJInteractionData mixer = MixingMap[key]; - - RealType sigmai = mixer.sigmai; - RealType epsilon = mixer.epsilon; - - ros = rij * sigmai; - + ros = *(idat.rij) * sigmai; + getLJfunc(ros, myPot, myDeriv); - - if (shiftedPot_) { - rcos = rcut * sigmai; + + if (idat.shiftedPot) { + rcos = *(idat.rcut) * sigmai; getLJfunc(rcos, myPotC, myDerivC); myDerivC = 0.0; - } else if (LJ::shiftedFrc_) { - rcos = rcut * sigmai; + } else if (idat.shiftedForce) { + rcos = *(idat.rcut) * sigmai; getLJfunc(rcos, myPotC, myDerivC); - myPotC = myPotC + myDerivC * (rij - rcut) * sigmai; + myPotC = myPotC + myDerivC * (*(idat.rij) - *(idat.rcut)) * sigmai; } else { myPotC = 0.0; - myDerivC = 0.0; + myDerivC = 0.0; } - - RealType pot_temp = vdwMult * epsilon * (myPot - myPotC); - vpair += pot_temp; - - RealType dudr = sw * vdwMult * epsilon * (myDeriv - myDerivC)*sigmai; - pot += sw * pot_temp; - f1 = d * dudr / rij; - - return; - - } - - void LJ::do_lj_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(); + RealType pot_temp = *(idat.vdwMult) * epsilon * (myPot - myPotC); + *(idat.vpair) += pot_temp; - AtomType* atype1 = LJMap[*atid1]; - AtomType* atype2 = LJMap[*atid2]; + 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); - 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; + return; } void LJ::getLJfunc(RealType r, RealType &pot, RealType &deriv) { @@ -376,36 +284,19 @@ namespace OpenMD { return; } - - void LJ::setLJDefaultCutoff(RealType *thisRcut, int *sP, int *sF) { - shiftedPot_ = (bool)(*sP); - shiftedFrc_ = (bool)(*sF); - } -} - -extern "C" { - -#define fortranGetSigma FC_FUNC(getsigma, GETSIGMA) -#define fortranGetEpsilon FC_FUNC(getepsilon, GETEPSILON) -#define fortranSetLJCutoff FC_FUNC(setljdefaultcutoff, SETLJDEFAULTCUTOFF) -#define fortranDoLJPair FC_FUNC(do_lj_pair, DO_LJ_PAIR) - - RealType fortranGetSigma(int* atid) { - return OpenMD::LJ::Instance()->getSigma(*atid); - } - RealType fortranGetEpsilon(int* atid) { - return OpenMD::LJ::Instance()->getEpsilon(*atid); - } - void fortranSetLJCutoff(RealType *rcut, int *shiftedPot, int *shiftedFrc) { - return OpenMD::LJ::Instance()->setLJDefaultCutoff(rcut, shiftedPot, - shiftedFrc); - } - void fortranDoLJPair(int *atid1, int *atid2, RealType *d, RealType *rij, - RealType *r2, RealType *rcut, RealType *sw, - RealType *vdwMult, RealType* vpair, RealType* pot, - RealType *f1){ + RealType LJ::getSuggestedCutoffRadius(pair atypes) { + if (!initialized_) initialize(); - return OpenMD::LJ::Instance()->do_lj_pair(atid1, atid2, d, rij, r2, rcut, - sw, vdwMult, vpair, pot, f1); + int atid1 = atypes.first->getIdent(); + int atid2 = atypes.second->getIdent(); + int ljtid1 = LJtids[atid1]; + int ljtid2 = LJtids[atid2]; + + if (ljtid1 == -1 || ljtid2 == -1) return 0.0; + else { + LJInteractionData mixer = MixingMap[ljtid1][ljtid2]; + return 2.5 * mixer.sigma; + } } + }