--- branches/development/src/nonbonded/EAM.cpp 2011/06/03 21:39:49 1575 +++ branches/development/src/nonbonded/EAM.cpp 2012/05/19 02:58:35 1711 @@ -36,7 +36,8 @@ * [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). + * [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). + * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). */ #include @@ -51,98 +52,13 @@ namespace OpenMD { namespace OpenMD { EAM::EAM() : name_("EAM"), initialized_(false), forceField_(NULL), - mixMeth_(eamJohnson), eamRcut_(0.0) {} + mixMeth_(eamJohnson), eamRcut_(0.0), haveCutoffRadius_(false) {} - EAMParam EAM::getEAMParam(AtomType* atomType) { - - // Do sanity checking on the AtomType we were passed before - // building any data structures: - if (!atomType->isEAM()) { - sprintf( painCave.errMsg, - "EAM::getEAMParam was passed an atomType (%s) that does not\n" - "\tappear to be an embedded atom method (EAM) atom.\n", - atomType->getName().c_str()); - painCave.severity = OPENMD_ERROR; - painCave.isFatal = 1; - simError(); - } - - GenericData* data = atomType->getPropertyByName("EAM"); - if (data == NULL) { - sprintf( painCave.errMsg, "EAM::getEAMParam could not find EAM\n" - "\tparameters for atomType %s.\n", - atomType->getName().c_str()); - painCave.severity = OPENMD_ERROR; - painCave.isFatal = 1; - simError(); - } - - EAMParamGenericData* eamData = dynamic_cast(data); - if (eamData == NULL) { - sprintf( painCave.errMsg, - "EAM::getEAMParam could not convert GenericData to EAMParam for\n" - "\tatom type %s\n", atomType->getName().c_str()); - painCave.severity = OPENMD_ERROR; - painCave.isFatal = 1; - simError(); - } - - return eamData->getData(); - } - - CubicSpline* EAM::getZ(AtomType* atomType) { - EAMParam eamParam = getEAMParam(atomType); - int nr = eamParam.nr; - RealType dr = eamParam.dr; - vector rvals; - - for (int i = 0; i < nr; i++) rvals.push_back(RealType(i) * dr); - - CubicSpline* cs = new CubicSpline(); - cs->addPoints(rvals, eamParam.Z); - return cs; - } - - RealType EAM::getRcut(AtomType* atomType) { - EAMParam eamParam = getEAMParam(atomType); - return eamParam.rcut; - } - - CubicSpline* EAM::getRho(AtomType* atomType) { - EAMParam eamParam = getEAMParam(atomType); - int nr = eamParam.nr; - RealType dr = eamParam.dr; - vector rvals; - - for (int i = 0; i < nr; i++) rvals.push_back(RealType(i) * dr); - - CubicSpline* cs = new CubicSpline(); - cs->addPoints(rvals, eamParam.rho); - return cs; - } - - CubicSpline* EAM::getF(AtomType* atomType) { - EAMParam eamParam = getEAMParam(atomType); - int nrho = eamParam.nrho; - RealType drho = eamParam.drho; - vector rhovals; - vector scaledF; - - for (int i = 0; i < nrho; i++) { - rhovals.push_back(RealType(i) * drho); - scaledF.push_back( eamParam.F[i] * 23.06054 ); - } - - CubicSpline* cs = new CubicSpline(); - cs->addPoints(rhovals, scaledF); - return cs; - } - - CubicSpline* EAM::getPhi(AtomType* atomType1, AtomType* atomType2) { - EAMParam eamParam1 = getEAMParam(atomType1); - EAMParam eamParam2 = getEAMParam(atomType2); - CubicSpline* z1 = getZ(atomType1); - CubicSpline* z2 = getZ(atomType2); + CubicSpline* EAM::getPhi(AtomType* atomType1, AtomType* atomType2) { + EAMAdapter ea1 = EAMAdapter(atomType1); + EAMAdapter ea2 = EAMAdapter(atomType2); + CubicSpline* z1 = ea1.getZ(); + CubicSpline* z2 = ea2.getZ(); // make the r grid: @@ -150,15 +66,15 @@ namespace OpenMD { // we need phi out to the largest value we'll encounter in the radial space; RealType rmax = 0.0; - rmax = max(rmax, eamParam1.rcut); - rmax = max(rmax, eamParam1.nr * eamParam1.dr); + rmax = max(rmax, ea1.getRcut()); + rmax = max(rmax, ea1.getNr() * ea1.getDr()); - rmax = max(rmax, eamParam2.rcut); - rmax = max(rmax, eamParam2.nr * eamParam2.dr); + rmax = max(rmax, ea2.getRcut()); + rmax = max(rmax, ea2.getNr() * ea2.getDr()); // use the smallest dr (finest grid) to build our grid: - RealType dr = min(eamParam1.dr, eamParam2.dr); + RealType dr = min(ea1.getDr(), ea2.getDr()); int nr = int(rmax/dr + 0.5); @@ -182,8 +98,8 @@ namespace OpenMD { // means that our phi grid goes out beyond the cutoff of the // pair potential - zi = r <= eamParam1.rcut ? z1->getValueAt(r) : 0.0; - zj = r <= eamParam2.rcut ? z2->getValueAt(r) : 0.0; + zi = r <= ea1.getRcut() ? z1->getValueAt(r) : 0.0; + zj = r <= ea2.getRcut() ? z2->getValueAt(r) : 0.0; phi = 331.999296 * (zi * zj) / r; @@ -195,6 +111,11 @@ namespace OpenMD { return cs; } + void EAM::setCutoffRadius( RealType rCut ) { + eamRcut_ = rCut; + haveCutoffRadius_ = true; + } + void EAM::initialize() { // set up the mixing method: @@ -272,22 +193,22 @@ namespace OpenMD { void EAM::addType(AtomType* atomType){ + EAMAdapter ea = EAMAdapter(atomType); EAMAtomData eamAtomData; - - eamAtomData.rho = getRho(atomType); - eamAtomData.F = getF(atomType); - eamAtomData.Z = getZ(atomType); - eamAtomData.rcut = getRcut(atomType); + eamAtomData.rho = ea.getRho(); + eamAtomData.F = ea.getF(); + eamAtomData.Z = ea.getZ(); + eamAtomData.rcut = ea.getRcut(); + // add it to the map: - AtomTypeProperties atp = atomType->getATP(); pair::iterator,bool> ret; - ret = EAMlist.insert( pair(atp.ident, atomType) ); + ret = EAMlist.insert( pair(atomType->getIdent(), atomType) ); if (ret.second == false) { sprintf( painCave.errMsg, "EAM already had a previous entry with ident %d\n", - atp.ident); + atomType->getIdent()); painCave.severity = OPENMD_INFO; painCave.isFatal = 0; simError(); @@ -353,19 +274,22 @@ namespace OpenMD { EAMAtomData data1 = EAMMap[idat.atypes.first]; EAMAtomData data2 = EAMMap[idat.atypes.second]; - + + if (haveCutoffRadius_) + if ( *(idat.rij) > eamRcut_) return; + if ( *(idat.rij) < data1.rcut) *(idat.rho1) += data1.rho->getValueAt( *(idat.rij)); - - + + if ( *(idat.rij) < data2.rcut) *(idat.rho2) += data2.rho->getValueAt( *(idat.rij)); - - return; + + return; } - + void EAM::calcFunctional(SelfData &sdat) { - + if (!initialized_) initialize(); EAMAtomData data1 = EAMMap[ sdat.atype ]; @@ -375,8 +299,10 @@ namespace OpenMD { *(sdat.frho) = result.first; *(sdat.dfrhodrho) = result.second; - sdat.pot[METALLIC_FAMILY] += result.first; - *(sdat.particlePot) += result.first; + (*(sdat.pot))[METALLIC_FAMILY] += result.first; + if (sdat.doParticlePot) { + *(sdat.particlePot) += result.first; + } return; } @@ -386,87 +312,89 @@ namespace OpenMD { if (!initialized_) initialize(); + if (haveCutoffRadius_) + if ( *(idat.rij) > eamRcut_) return; + pair res; - if ( *(idat.rij) < eamRcut_) { - - EAMAtomData data1 = EAMMap[idat.atypes.first]; - EAMAtomData data2 = EAMMap[idat.atypes.second]; - - // get type-specific cutoff radii - - RealType rci = data1.rcut; - RealType rcj = data2.rcut; + EAMAtomData data1 = EAMMap[idat.atypes.first]; + EAMAtomData data2 = EAMMap[idat.atypes.second]; + + // get type-specific cutoff radii + + RealType rci = data1.rcut; + RealType rcj = data2.rcut; + + RealType rha(0.0), drha(0.0), rhb(0.0), drhb(0.0); + RealType pha(0.0), dpha(0.0), phb(0.0), dphb(0.0); + RealType phab(0.0), dvpdr(0.0); + RealType drhoidr, drhojdr, dudr; + + if ( *(idat.rij) < rci) { + res = data1.rho->getValueAndDerivativeAt( *(idat.rij)); + rha = res.first; + drha = res.second; - RealType rha(0.0), drha(0.0), rhb(0.0), drhb(0.0); - RealType pha(0.0), dpha(0.0), phb(0.0), dphb(0.0); - RealType phab(0.0), dvpdr(0.0); - RealType drhoidr, drhojdr, dudr; + res = MixingMap[make_pair(idat.atypes.first, idat.atypes.first)].phi->getValueAndDerivativeAt( *(idat.rij) ); + pha = res.first; + dpha = res.second; + } + + if ( *(idat.rij) < rcj) { + res = data2.rho->getValueAndDerivativeAt( *(idat.rij) ); + rhb = res.first; + drhb = res.second; - if ( *(idat.rij) < rci) { - res = data1.rho->getValueAndDerivativeAt( *(idat.rij)); - rha = res.first; - drha = res.second; + res = MixingMap[make_pair(idat.atypes.second, idat.atypes.second)].phi->getValueAndDerivativeAt( *(idat.rij) ); + phb = res.first; + dphb = res.second; + } - res = MixingMap[make_pair(idat.atypes.first, idat.atypes.first)].phi->getValueAndDerivativeAt( *(idat.rij) ); - pha = res.first; - dpha = res.second; + switch(mixMeth_) { + case eamJohnson: + + if ( *(idat.rij) < rci) { + phab = phab + 0.5 * (rhb / rha) * pha; + dvpdr = dvpdr + 0.5*((rhb/rha)*dpha + + pha*((drhb/rha) - (rhb*drha/rha/rha))); } - + + + if ( *(idat.rij) < rcj) { - res = data2.rho->getValueAndDerivativeAt( *(idat.rij) ); - rhb = res.first; - drhb = res.second; - - res = MixingMap[make_pair(idat.atypes.second, idat.atypes.second)].phi->getValueAndDerivativeAt( *(idat.rij) ); - phb = res.first; - dphb = res.second; + phab = phab + 0.5 * (rha / rhb) * phb; + dvpdr = dvpdr + 0.5 * ((rha/rhb)*dphb + + phb*((drha/rhb) - (rha*drhb/rhb/rhb))); } - - switch(mixMeth_) { - case eamJohnson: - - if ( *(idat.rij) < rci) { - phab = phab + 0.5 * (rhb / rha) * pha; - dvpdr = dvpdr + 0.5*((rhb/rha)*dpha + - pha*((drhb/rha) - (rhb*drha/rha/rha))); - } - - - - if ( *(idat.rij) < rcj) { - phab = phab + 0.5 * (rha / rhb) * phb; - dvpdr = dvpdr + 0.5 * ((rha/rhb)*dphb + - phb*((drha/rhb) - (rha*drhb/rhb/rhb))); - } - - break; - - case eamDaw: - res = MixingMap[idat.atypes].phi->getValueAndDerivativeAt( *(idat.rij)); - phab = res.first; - dvpdr = res.second; - - break; - case eamUnknown: - default: - - sprintf(painCave.errMsg, - "EAM::calcForce hit a mixing method it doesn't know about!\n" - ); - painCave.severity = OPENMD_ERROR; - painCave.isFatal = 1; - simError(); - - } - drhoidr = drha; - drhojdr = drhb; - - dudr = drhojdr* *(idat.dfrho1) + drhoidr* *(idat.dfrho2) + dvpdr; - - *(idat.f1) = *(idat.d) * dudr / *(idat.rij); + break; + + case eamDaw: + res = MixingMap[idat.atypes].phi->getValueAndDerivativeAt( *(idat.rij)); + phab = res.first; + dvpdr = res.second; + + break; + case eamUnknown: + default: + + sprintf(painCave.errMsg, + "EAM::calcForce hit a mixing method it doesn't know about!\n" + ); + painCave.severity = OPENMD_ERROR; + painCave.isFatal = 1; + simError(); + + } + + drhoidr = drha; + drhojdr = drhb; + + dudr = drhojdr* *(idat.dfrho1) + drhoidr* *(idat.dfrho2) + dvpdr; + + *(idat.f1) += *(idat.d) * dudr / *(idat.rij); + if (idat.doParticlePot) { // particlePot is the difference between the full potential and // the full potential without the presence of a particular // particle (atom1). @@ -475,18 +403,18 @@ namespace OpenMD { // need to recompute the density at atom2 assuming atom1 didn't // contribute. This then requires recomputing the density // functional for atom2 as well. - + *(idat.particlePot1) += data2.F->getValueAt( *(idat.rho2) - rha ) - *(idat.frho2); - + *(idat.particlePot2) += data1.F->getValueAt( *(idat.rho1) - rhb) - *(idat.frho1); - - idat.pot[METALLIC_FAMILY] += phab; - - *(idat.vpair) += phab; } - + + (*(idat.pot))[METALLIC_FAMILY] += phab; + + *(idat.vpair) += phab; + return; }