--- branches/development/src/nonbonded/LJ.cpp 2010/10/02 19:53:32 1502 +++ branches/development/src/nonbonded/LJ.cpp 2011/04/27 18:38:15 1549 @@ -240,53 +240,54 @@ namespace OpenMD { } } - void LJ::calcForce(InteractionData idat) { + void LJ::calcForce(InteractionData &idat) { if (!initialized_) initialize(); - 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(idat.atype1, - idat.atype2); - LJInteractionData mixer = MixingMap[key]; - - RealType sigmai = mixer.sigmai; - RealType epsilon = mixer.epsilon; + map, LJInteractionData>::iterator it; + it = MixingMap.find(idat.atypes); - - ros = idat.rij * sigmai; - - getLJfunc(ros, myPot, myDeriv); + if (it != MixingMap.end()) { + + LJInteractionData mixer = (*it).second; + + 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; + + ros = idat.rij * sigmai; + + getLJfunc(ros, myPot, myDeriv); + + if (shiftedPot_) { + rcos = idat.rcut * sigmai; + getLJfunc(rcos, myPotC, myDerivC); + myDerivC = 0.0; + } else if (LJ::shiftedFrc_) { + rcos = idat.rcut * sigmai; + getLJfunc(rcos, myPotC, myDerivC); + myPotC = myPotC + myDerivC * (idat.rij - idat.rcut) * sigmai; + } else { + myPotC = 0.0; + myDerivC = 0.0; + } - if (shiftedPot_) { - rcos = idat.rcut * sigmai; - getLJfunc(rcos, myPotC, myDerivC); - myDerivC = 0.0; - } else if (LJ::shiftedFrc_) { - rcos = idat.rcut * sigmai; - getLJfunc(rcos, 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[0] += idat.sw * pot_temp; + idat.f1 = idat.d * dudr / idat.rij; } - - RealType pot_temp = idat.vdwMult * epsilon * (myPot - myPotC); - idat.vpair += pot_temp; - - RealType dudr = idat.sw * idat.vdwMult * epsilon * (myDeriv - - myDerivC)*sigmai; - - idat.pot += idat.sw * pot_temp; - idat.f1 = idat.d * dudr / idat.rij; - return; - } void LJ::getLJfunc(RealType r, RealType &pot, RealType &deriv) { @@ -304,5 +305,16 @@ namespace OpenMD { return; } + RealType LJ::getSuggestedCutoffRadius(pair atypes) { + if (!initialized_) initialize(); + map, LJInteractionData>::iterator it; + it = MixingMap.find(atypes); + if (it == MixingMap.end()) + return 0.0; + else { + LJInteractionData mixer = (*it).second; + return 2.5 * mixer.sigma; + } + } }