--- branches/development/src/nonbonded/LJ.cpp 2010/07/19 18:59:59 1471 +++ branches/development/src/nonbonded/LJ.cpp 2010/10/02 19:53:32 1502 @@ -46,25 +46,11 @@ #include "nonbonded/LJ.hpp" #include "utils/simError.h" - 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), shiftedPot_(false), + shiftedFrc_(false), forceField_(NULL) {} - LJ* LJ::Instance() { - if (!_instance) { - _instance = new LJ(); - } - return _instance; - } - LJParam LJ::getLJParam(AtomType* atomType) { // Do sanity checking on the AtomType we were passed before @@ -106,28 +92,12 @@ namespace OpenMD { return ljParam.sigma; } - RealType LJ::getSigma(int atid) { - 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") @@ -139,22 +109,6 @@ namespace OpenMD { RealType LJ::getEpsilon(AtomType* atomType) { LJParam ljParam = getLJParam(atomType); return ljParam.epsilon; - } - - RealType LJ::getEpsilon(int atid) { - 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; - - return getEpsilon(atype); } RealType LJ::getEpsilon(AtomType* atomType1, AtomType* atomType2) { @@ -164,27 +118,27 @@ namespace OpenMD { } void LJ::initialize() { - ForceField::AtomTypeContainer atomTypes = forceField_->getAtomTypes(); + ForceField::AtomTypeContainer* atomTypes = forceField_->getAtomTypes(); ForceField::AtomTypeContainer::MapTypeIterator i; AtomType* at; - for (at = atomTypes.beginType(i); at != NULL; - at = atomTypes.nextType(i)) { + for (at = atomTypes->beginType(i); at != NULL; + at = atomTypes->nextType(i)) { if (at->isLennardJones()) addType(at); } - ForceField::NonBondedInteractionTypeContainer nbiTypes = forceField_->getNonBondedInteractionTypes(); + ForceField::NonBondedInteractionTypeContainer* nbiTypes = forceField_->getNonBondedInteractionTypes(); ForceField::NonBondedInteractionTypeContainer::MapTypeIterator j; NonBondedInteractionType* nbt; - for (nbt = nbiTypes.beginType(j); nbt != NULL; - nbt = nbiTypes.nextType(j)) { + for (nbt = nbiTypes->beginType(j); nbt != NULL; + nbt = nbiTypes->nextType(j)) { if (nbt->isLennardJones()) { - std::pair atypes = nbt->getAtomTypes(); + pair atypes = nbt->getAtomTypes(); GenericData* data = nbt->getPropertyByName("LennardJones"); if (data == NULL) { @@ -225,11 +179,12 @@ namespace OpenMD { void LJ::addType(AtomType* atomType){ RealType sigma1 = getSigma(atomType); RealType epsilon1 = getEpsilon(atomType); - + // add it to the map: AtomTypeProperties atp = atomType->getATP(); - std::pair::iterator,bool> ret; - ret = LJMap.insert( std::pair(atp.ident, atomType) ); + + pair::iterator,bool> ret; + ret = LJMap.insert( pair(atp.ident, atomType) ); if (ret.second == false) { sprintf( painCave.errMsg, "LJ already had a previous entry with ident %d\n", @@ -285,10 +240,8 @@ namespace OpenMD { } } - 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(); RealType ros; @@ -298,101 +251,58 @@ namespace OpenMD { RealType myDeriv = 0.0; RealType myDerivC = 0.0; - std::pair key = std::make_pair(at1, at2); + std::pair key = std::make_pair(idat.atype1, + idat.atype2); 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; + rcos = idat.rcut * sigmai; getLJfunc(rcos, myPotC, myDerivC); myDerivC = 0.0; } else if (LJ::shiftedFrc_) { - rcos = rcut * sigmai; + 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; } - RealType pot_temp = vdwMult * epsilon * (myPot - myPotC); - vpair += pot_temp; + RealType pot_temp = idat.vdwMult * epsilon * (myPot - myPotC); + idat.vpair += pot_temp; - RealType dudr = sw * vdwMult * epsilon * (myDeriv - myDerivC)*sigmai; + RealType dudr = idat.sw * idat.vdwMult * epsilon * (myDeriv - + myDerivC)*sigmai; - pot += sw * pot_temp; - f1 = d * dudr / rij; + idat.pot += idat.sw * pot_temp; + idat.f1 = idat.d * dudr / idat.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(); - - AtomType* atype1 = LJMap[*atid1]; - AtomType* atype2 = LJMap[*atid2]; - - 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); - return; } - void LJ::getLJfunc(const RealType r, RealType pot, RealType deriv) { + void LJ::getLJfunc(RealType r, RealType &pot, RealType &deriv) { + RealType ri = 1.0 / r; RealType ri2 = ri * ri; RealType ri6 = ri2 * ri2 * ri2; RealType ri7 = ri6 * ri; RealType ri12 = ri6 * ri6; RealType ri13 = ri12 * ri; - + pot = 4.0 * (ri12 - ri6); deriv = 24.0 * (ri7 - 2.0 * ri13); + 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){ - - return OpenMD::LJ::Instance()->do_lj_pair(atid1, atid2, d, rij, r2, rcut, - sw, vdwMult, vpair, pot, f1); - } -}