--- branches/development/src/nonbonded/EAM.cpp 2010/07/23 20:45:40 1478 +++ branches/development/src/nonbonded/EAM.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,124 +46,40 @@ #include #include "nonbonded/EAM.hpp" #include "utils/simError.h" +#include "types/NonBondedInteractionType.hpp" namespace OpenMD { - bool EAM::initialized_ = false; - ForceField* EAM::forceField_ = NULL; - std::map EAM::EAMlist; - std::map EAM::EAMMap; - std::map, EAMInteractionData> EAM::MixingMap; + EAM::EAM() : name_("EAM"), initialized_(false), forceField_(NULL), + mixMeth_(eamJohnson), eamRcut_(0.0), haveCutoffRadius_(false) {} - EAM* EAM::_instance = NULL; + CubicSpline* EAM::getPhi(AtomType* atomType1, AtomType* atomType2) { + EAMAdapter ea1 = EAMAdapter(atomType1); + EAMAdapter ea2 = EAMAdapter(atomType2); + CubicSpline* z1 = ea1.getZ(); + CubicSpline* z2 = ea2.getZ(); - EAM* EAM::Instance() { - if (!_instance) { - _instance = new EAM(); - } - return _instance; - } - - 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(); - } + // make the r grid: - 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(i * dr); - - CubicSpline* cs = new CubicSpline(); - cs->addPoints(rvals, eamParam.Z); - return cs; - } - CubicSpline* EAM::getRho(AtomType* atomType) { - EAMParam eamParam = getEAMParam(atomType); - int nr = eamParam.nr; - RealType dr = eamParam.dr; - vector rvals; + // we need phi out to the largest value we'll encounter in the radial space; - for (int i = 0; i < nr; i++) rvals.push_back(i * dr); - - CubicSpline* cs = new CubicSpline(); - cs->addPoints(rvals, eamParam.rho); - return cs; - } + RealType rmax = 0.0; + rmax = max(rmax, ea1.getRcut()); + rmax = max(rmax, ea1.getNr() * ea1.getDr()); - 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(i * drho); - scaledF.push_back( eamParam.F[i] * 23.06054 ); - } - - CubicSpline* cs = new CubicSpline(); - cs->addPoints(rhovals, eamParam.F); - 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); + rmax = max(rmax, ea2.getRcut()); + rmax = max(rmax, ea2.getNr() * ea2.getDr()); - // make the r grid: + // use the smallest dr (finest grid) to build our grid: - // set rcut to be the smaller of the two atomic rcuts + RealType dr = min(ea1.getDr(), ea2.getDr()); - RealType rcut = eamParam1.rcut < eamParam2.rcut ? - eamParam1.rcut : eamParam2.rcut; + int nr = int(rmax/dr + 0.5); - // use the smallest dr (finest grid) to build our grid: - - RealType dr = eamParam1.dr < eamParam2.dr ? eamParam1.dr : eamParam2.dr; - int nr = int(rcut/dr); vector rvals; - for (int i = 0; i < nr; i++) rvals.push_back(i*dr); + for (int i = 0; i < nr; i++) rvals.push_back(RealType(i*dr)); // construct the pair potential: @@ -173,12 +90,19 @@ namespace OpenMD { phivals.push_back(0.0); - for (int i = 1; i < rvals.size(); i++ ) { + for (unsigned int i = 1; i < rvals.size(); i++ ) { r = rvals[i]; - zi = z1->getValueAt(r); - zj = z2->getValueAt(r); + // only use z(r) if we're inside this atom's cutoff radius, + // otherwise, we'll use zero for the charge. This effectively + // means that our phi grid goes out beyond the cutoff of the + // pair potential + + zi = r <= ea1.getRcut() ? z1->getValueAt(r) : 0.0; + zj = r <= ea2.getRcut() ? z2->getValueAt(r) : 0.0; + phi = 331.999296 * (zi * zj) / r; + phivals.push_back(phi); } @@ -187,12 +111,18 @@ namespace OpenMD { return cs; } + void EAM::setCutoffRadius( RealType rCut ) { + eamRcut_ = rCut; + haveCutoffRadius_ = true; + } + void EAM::initialize() { // set up the mixing method: - ForceFieldOptions ffo = forceField_->getForceFieldOptions(); - string EAMMixMeth = toUpperCopy(ffo.getEAMMixingMethod()); - + ForceFieldOptions& fopts = forceField_->getForceFieldOptions(); + string EAMMixMeth = fopts.getEAMMixingMethod(); + toUpper(EAMMixMeth); + if (EAMMixMeth == "JOHNSON") mixMeth_ = eamJohnson; else if (EAMMixMeth == "DAW") @@ -222,7 +152,7 @@ namespace OpenMD { if (nbt->isEAM()) { - std::pair atypes = nbt->getAtomTypes(); + pair atypes = nbt->getAtomTypes(); GenericData* data = nbt->getPropertyByName("EAM"); if (data == NULL) { @@ -247,9 +177,9 @@ namespace OpenMD { simError(); } - EAMMix eamParam = eamData->getData(); + EAMMixingParam eamParam = eamData->getData(); - vector phiAB = eamParam.phiAB; + vector phiAB = eamParam.phi; RealType dr = eamParam.dr; int nr = eamParam.nr; @@ -263,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(); - std::pair::iterator,bool> ret; - ret = EAMlist.insert( std::pair(atp.ident, atomType) ); + pair::iterator,bool> ret; + 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(); @@ -288,18 +218,18 @@ namespace OpenMD { // Now, iterate over all known types and add to the mixing map: - std::map::iterator it; + map::iterator it; for( it = EAMMap.begin(); it != EAMMap.end(); ++it) { - AtomType* atype2 = (*it).second; + AtomType* atype2 = (*it).first; EAMInteractionData mixer; mixer.phi = getPhi(atomType, atype2); mixer.explicitlySet = false; - std::pair key1, key2; - key1 = std::make_pair(atomType, atype2); - key2 = std::make_pair(atype2, atomType); + pair key1, key2; + key1 = make_pair(atomType, atype2); + key2 = make_pair(atype2, atomType); MixingMap[key1] = mixer; if (key2 != key1) { @@ -319,17 +249,17 @@ namespace OpenMD { EAMInteractionData mixer; CubicSpline* cs = new CubicSpline(); - vector rvals; + vector rVals; - for (int i = 0; i < nr; i++) rvals.push_back(i * dr); + for (int i = 0; i < nr; i++) rVals.push_back(i * dr); cs->addPoints(rVals, phiVals); mixer.phi = cs; mixer.explicitlySet = true; - std::pair key1, key2; - key1 = std::make_pair(atype1, atype2); - key2 = std::make_pair(atype2, atype1); + pair key1, key2; + key1 = make_pair(atype1, atype2); + key2 = make_pair(atype2, atype1); MixingMap[key1] = mixer; if (key2 != key1) { @@ -338,245 +268,178 @@ namespace OpenMD { return; } - void EAM::calcDensity(AtomType* at1, AtomType* at2, Vector3d d, - RealType rij, RealType r2, RealType rho_i_at_j, - RealType rho_j_at_i) { - + void EAM::calcDensity(InteractionData &idat) { + if (!initialized_) initialize(); - - EAMAtomData data1 = EAMMap[at1]; - EAMAtomData data2 = EAMMap[at2]; - - if (rij < data1.rcut) rho_i_at_j = data1.rho->getValueAt(rij); - if (rij < data2.rcut) rho_j_at_i = data2.rho->getValueAt(rij); - return; + + 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; } - - void EAM::calcFunctional(AtomType* at1, RealType rho, RealType frho, - RealType dfrhodrho) { - + + void EAM::calcFunctional(SelfData &sdat) { + if (!initialized_) initialize(); - EAMAtomData data1 = EAMMap[at1]; + EAMAtomData data1 = EAMMap[ sdat.atype ]; - pair result = data1.F->getValueAndDerivativeAt(rho); + pair result = data1.F->getValueAndDerivativeAt( *(sdat.rho) ); - frho = result.first; - dfrhodrho = result.second; + *(sdat.frho) = result.first; + *(sdat.dfrhodrho) = result.second; + + (*(sdat.pot))[METALLIC_FAMILY] += result.first; + if (sdat.doParticlePot) { + *(sdat.particlePot) += result.first; + } + return; } - void EAM::calcForce(AtomType* at1, AtomType* at2, Vector3d d, - RealType rij, RealType r2, RealType sw, - RealType &vpair, RealType &pot, Vector3d &f1, - RealType rho1, RealType rho2, RealType dfrho1, - RealType dfrho2, RealType fshift1, RealType fshift2) { + void EAM::calcForce(InteractionData &idat) { if (!initialized_) initialize(); - + + if (haveCutoffRadius_) + if ( *(idat.rij) > eamRcut_) return; + pair res; - if (rij < eamRcut_) { - - EAMAtomData data1 = EAMMap[at1]; - EAMAtomData data2 = EAMMap[at2]; - - // get type-specific cutoff radii - - RealType rci = data1.rcut; - RealType rcj = data2.rcut; - - RealType rha, drha, rhb, drhb; - RealType pha, dpha, phb, dphb; - RealType phab, dvpdr; - RealType drhoidr, drhojdr, dudr; - - if (rij < rci) { - res = data1.rho->getValueAndDerivativeAt(rij); - rha = res.first; - drha = res.second; - - res = MixingMap[make_pair(at1, at1)].phi->getValueAndDerivativeAt(rij); - pha = res.first; - dpha = res.second; - } + 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; + + 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; + + res = MixingMap[make_pair(idat.atypes.second, idat.atypes.second)].phi->getValueAndDerivativeAt( *(idat.rij) ); + phb = res.first; + dphb = res.second; + } - if (rij < rcj) { - res = data2.rho->getValueAndDerivativeAt(rij); - rhb = res.first; - drhb = res.second; - - res = MixingMap[make_pair(at2, at2)].phi->getValueAndDerivativeAt(rij); - phb = res.first; - dphb = 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))); } - - phab = 0.0; - dvpdr = 0.0; - - switch(mixMeth_) { - case eamJohnson: - - if (rij < rci) { - phab = phab + 0.5 * (rhb / rha) * pha; - dvpdr = dvpdr + 0.5*((rhb/rha)*dpha + - pha*((drhb/rha) - (rhb*drha/rha/rha))); - } - - if (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[make_pair(at1,at2)].phi->getValueAndDerivativeAt(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(); - + + + + 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))); } - drhoidr = drha; - drhojdr = drhb; - - dudr = drhojdr*dfrhodrho_i + drhoidr*dfrhodrho_j + dvpdr; - - f1 = d * dudr / 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); - // particle_pot is the difference between the full potential - // and the full potential without the presence of a particular + if (idat.doParticlePot) { + // particlePot is the difference between the full potential and + // the full potential without the presence of a particular // particle (atom1). // - // This reduces the density at other particle locations, so - // we need to recompute the density at atom2 assuming atom1 - // didn't contribute. This then requires recomputing the - // density functional for atom2 as well. - // - // Most of the particle_pot heavy lifting comes from the - // pair interaction, and will be handled by vpair. - - fshift_i = data1.F->getValueAt( rho_i - rhb ); - fshift_j = data1.F->getValueAt( rho_j - rha ); - - pot += phab; - - vpair += phab; + // This reduces the density at other particle locations, so we + // 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); } - - return; - } - - - void EAM::calc_eam_prepair_rho(int *atid1, int *atid2, RealType *d, - RealType *rij, RealType *r2, - RealType* rho_i_at_j, RealType* rho_j_at_i){ - if (!initialized_) initialize(); - - AtomType* atype1 = EAMlist[*atid1]; - AtomType* atype2 = EAMlist[*atid2]; + (*(idat.pot))[METALLIC_FAMILY] += phab; - Vector3d disp(d[0], d[1], d[2]); - - calcDensity(atype1, atype2, disp, *rij, *r2, *rho_i_at_j, *rho_j_at_i); - - return; - } - - void EAM::calc_eam_preforce_Frho(int *atid1, RealType *rho, RealType *frho, - RealType *dfrhodrho) { - - if (!initialized_) initialize(); - - AtomType* atype1 = EAMlist[*atid1]; - - calcFunctional(atype1, *rho, *frho, *dfrhodrho); + *(idat.vpair) += phab; + + return; - return; } - void EAM::do_eam_pair(int *atid1, int *atid2, RealType *d, RealType *rij, - RealType *r2, RealType *sw, RealType *vpair, - RealType *pot, RealType *f1, RealType *rho1, - RealType *rho2, RealType *dfrho1, RealType *dfrho2, - RealType *fshift1, RealType *fshift2) { + RealType EAM::getSuggestedCutoffRadius(pair atypes) { + if (!initialized_) initialize(); - if (!initialized_) initialize(); - - AtomType* atype1 = EAMMap[*atid1]; - AtomType* atype2 = EAMMap[*atid2]; - - Vector3d disp(d[0], d[1], d[2]); - Vector3d frc(f1[0], f1[1], f1[2]); - - calcForce(atype1, atype2, disp, *rij, *r2, *sw, *vpair, *pot, frc, - *rho1, *rho2, *dfrho1, *dfrho2, *fshift1, *fshift2); - - f1[0] = frc.x(); - f1[1] = frc.y(); - f1[2] = frc.z(); + RealType cut = 0.0; - return; - } - - void EAM::setCutoffEAM(RealType *thisRcut) { - eamRcut_ = thisRcut; - } -} + map::iterator it; -extern "C" { - -#define fortranCalcDensity FC_FUNC(calc_eam_prepair_rho, CALC_EAM_PREPAIR_RHO) -#define fortranCalcFunctional FC_FUNC(calc_eam_preforce_frho, CALC_EAM_PREFORCE_FRHO) -#define fortranCalcForce FC_FUNC(do_eam_pair, DO_EAM_PAIR) -#define fortranSetCutoffEAM FC_FUNC(setcutoffeam, SETCUTOFFEAM) - - RealType fortranCalcDensity(int *atid1, int *atid2, RealType *d, - RealType *rij, RealType *r2, - RealType *rho_i_at_j, RealType *rho_j_at_i) { + it = EAMMap.find(atypes.first); + if (it != EAMMap.end()) { + EAMAtomData data1 = (*it).second; + cut = data1.rcut; + } - return OpenMD::EAM::Instance()->calc_eam_prepair_rho(*atid1, *atid2, *d, - *rij, *r2, - *rho_i_at_j, - *rho_j_at_i); - } - RealType fortranCalcFunctional(int *atid1, RealType *rho, RealType *frho, - RealType *dfrhodrho) { + it = EAMMap.find(atypes.second); + if (it != EAMMap.end()) { + EAMAtomData data2 = (*it).second; + if (data2.rcut > cut) + cut = data2.rcut; + } - return OpenMD::EAM::Instance()->calc_eam_preforce_Frho(*atid1, - *rho, - *frho, - *dfrhodrho); - + return cut; } - void fortranSetEAMCutoff(RealType *rcut) { - return OpenMD::EAM::Instance()->setCutoffEAM(rcut); - } - void fortranDoEAMPair(int *atid1, int *atid2, RealType *d, RealType *rij, - RealType *r2, RealType *sw, RealType *vpair, - RealType *pot, RealType *f1, RealType *rho1, - RealType *rho2, RealType *dfrho1, RealType *dfrho2, - RealType *fshift1, RealType *fshift2){ - - return OpenMD::EAM::Instance()->do_eam_pair(*atid1, *atid2, *d, *rij, - *r2, *sw, *vpair, - *pot, *f1, *rho1, - *rho2, *dfrho1, *dfrho2, - *fshift1, *fshift2); - } } +