--- branches/development/src/nonbonded/EAM.cpp 2010/07/26 19:50:53 1480 +++ branches/development/src/nonbonded/EAM.cpp 2010/10/03 22:18:59 1505 @@ -50,24 +50,9 @@ namespace OpenMD { namespace OpenMD { - bool EAM::initialized_ = false; - RealType EAM::eamRcut_ = 0.0; - EAMMixingMethod EAM::mixMeth_ = eamJohnson; - 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) {} - EAM* EAM::_instance = NULL; - - 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 @@ -111,7 +96,7 @@ namespace OpenMD { RealType dr = eamParam.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); CubicSpline* cs = new CubicSpline(); cs->addPoints(rvals, eamParam.Z); @@ -129,7 +114,7 @@ namespace OpenMD { RealType dr = eamParam.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); CubicSpline* cs = new CubicSpline(); cs->addPoints(rvals, eamParam.rho); @@ -144,12 +129,12 @@ namespace OpenMD { vector scaledF; for (int i = 0; i < nrho; i++) { - rhovals.push_back(i * drho); + rhovals.push_back(RealType(i) * drho); scaledF.push_back( eamParam.F[i] * 23.06054 ); } CubicSpline* cs = new CubicSpline(); - cs->addPoints(rhovals, eamParam.F); + cs->addPoints(rhovals, scaledF); return cs; } @@ -161,17 +146,24 @@ namespace OpenMD { // make the r grid: - // set rcut to be the smaller of the two atomic rcuts - RealType rcut = eamParam1.rcut < eamParam2.rcut ? - eamParam1.rcut : eamParam2.rcut; + // 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, eamParam2.rcut); + rmax = max(rmax, eamParam2.nr * eamParam2.dr); + // 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); + RealType dr = min(eamParam1.dr, eamParam2.dr); + + int nr = int(rmax/dr + 0.5); + 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: @@ -184,10 +176,17 @@ namespace OpenMD { for (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 <= eamParam1.rcut ? z1->getValueAt(r) : 0.0; + zj = r <= eamParam2.rcut ? z2->getValueAt(r) : 0.0; + phi = 331.999296 * (zi * zj) / r; + phivals.push_back(phi); } @@ -200,7 +199,7 @@ namespace OpenMD { // set up the mixing method: ForceFieldOptions& fopts = forceField_->getForceFieldOptions(); - std::string EAMMixMeth = fopts.getEAMMixingMethod(); + string EAMMixMeth = fopts.getEAMMixingMethod(); toUpper(EAMMixMeth); if (EAMMixMeth == "JOHNSON") @@ -232,7 +231,7 @@ namespace OpenMD { if (nbt->isEAM()) { - std::pair atypes = nbt->getAtomTypes(); + pair atypes = nbt->getAtomTypes(); GenericData* data = nbt->getPropertyByName("EAM"); if (data == NULL) { @@ -283,8 +282,8 @@ namespace OpenMD { // 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(atp.ident, atomType) ); if (ret.second == false) { sprintf( painCave.errMsg, "EAM already had a previous entry with ident %d\n", @@ -298,7 +297,7 @@ 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).first; @@ -307,9 +306,9 @@ namespace OpenMD { 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) { @@ -337,9 +336,9 @@ namespace OpenMD { 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) { @@ -348,49 +347,46 @@ namespace OpenMD { return; } - void EAM::calcDensity(AtomType* at1, AtomType* at2, const RealType rij, - RealType &rho_i_at_j, RealType &rho_j_at_i) { + void EAM::calcDensity(DensityData ddat) { if (!initialized_) initialize(); - EAMAtomData data1 = EAMMap[at1]; - EAMAtomData data2 = EAMMap[at2]; + EAMAtomData data1 = EAMMap[ddat.atype1]; + EAMAtomData data2 = EAMMap[ddat.atype2]; - if (rij < data1.rcut) rho_i_at_j = data1.rho->getValueAt(rij); - if (rij < data2.rcut) rho_j_at_i = data2.rho->getValueAt(rij); + if (ddat.rij < data1.rcut) + ddat.rho_i_at_j = data1.rho->getValueAt(ddat.rij); + + if (ddat.rij < data2.rcut) + ddat.rho_j_at_i = data2.rho->getValueAt(ddat.rij); + return; } - void EAM::calcFunctional(AtomType* at1, RealType rho, RealType &frho, - RealType &dfrhodrho) { + void EAM::calcFunctional(FunctionalData fdat) { if (!initialized_) initialize(); - EAMAtomData data1 = EAMMap[at1]; + EAMAtomData data1 = EAMMap[fdat.atype]; - pair result = data1.F->getValueAndDerivativeAt(rho); + pair result = data1.F->getValueAndDerivativeAt(fdat.rho); - frho = result.first; - dfrhodrho = result.second; + fdat.frho = result.first; + fdat.dfrhodrho = result.second; return; } - void EAM::calcForce(AtomType* at1, AtomType* at2, Vector3d d, - RealType rij, RealType r2, RealType sw, - RealType &vpair, RealType &pot, Vector3d &f1, - RealType rho_i, RealType rho_j, - RealType dfrhodrho_i, RealType dfrhodrho_j, - RealType &fshift_i, RealType &fshift_j) { + void EAM::calcForce(InteractionData idat) { if (!initialized_) initialize(); - + pair res; - if (rij < eamRcut_) { + if (idat.rij < eamRcut_) { - EAMAtomData data1 = EAMMap[at1]; - EAMAtomData data2 = EAMMap[at2]; + EAMAtomData data1 = EAMMap[idat.atype1]; + EAMAtomData data2 = EAMMap[idat.atype2]; // get type-specific cutoff radii @@ -402,22 +398,22 @@ namespace OpenMD { RealType phab, dvpdr; RealType drhoidr, drhojdr, dudr; - if (rij < rci) { - res = data1.rho->getValueAndDerivativeAt(rij); + if (idat.rij < rci) { + res = data1.rho->getValueAndDerivativeAt(idat.rij); rha = res.first; drha = res.second; - res = MixingMap[make_pair(at1, at1)].phi->getValueAndDerivativeAt(rij); + res = MixingMap[make_pair(idat.atype1, idat.atype1)].phi->getValueAndDerivativeAt(idat.rij); pha = res.first; dpha = res.second; } - if (rij < rcj) { - res = data2.rho->getValueAndDerivativeAt(rij); + if (idat.rij < rcj) { + res = data2.rho->getValueAndDerivativeAt(idat.rij); rhb = res.first; drhb = res.second; - res = MixingMap[make_pair(at2, at2)].phi->getValueAndDerivativeAt(rij); + res = MixingMap[make_pair(idat.atype2, idat.atype2)].phi->getValueAndDerivativeAt(idat.rij); phb = res.first; dphb = res.second; } @@ -428,13 +424,13 @@ namespace OpenMD { switch(mixMeth_) { case eamJohnson: - if (rij < rci) { + 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 (rij < rcj) { + 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))); @@ -443,8 +439,7 @@ namespace OpenMD { break; case eamDaw: - - res = MixingMap[make_pair(at1,at2)].phi->getValueAndDerivativeAt(rij); + res = MixingMap[make_pair(idat.atype1,idat.atype2)].phi->getValueAndDerivativeAt(idat.rij); phab = res.first; dvpdr = res.second; @@ -464,9 +459,9 @@ namespace OpenMD { drhoidr = drha; drhojdr = drhb; - dudr = drhojdr*dfrhodrho_i + drhoidr*dfrhodrho_j + dvpdr; + dudr = drhojdr*idat.dfrho1 + drhoidr*idat.dfrho2 + dvpdr; - f1 = d * dudr / rij; + 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 @@ -480,121 +475,39 @@ namespace OpenMD { // 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 ); + idat.fshift1 = data1.F->getValueAt( idat.rho1 - rhb ); + idat.fshift2 = data1.F->getValueAt( idat.rho2 - rha ); - pot += phab; + idat.pot += phab; - vpair += phab; + idat.vpair += phab; } return; } + RealType EAM::getSuggestedCutoffRadius(AtomType* at1, AtomType* at2) { + if (!initialized_) initialize(); - void EAM::calc_eam_prepair_rho(int *atid1, int *atid2, RealType *rij, - RealType* rho_i_at_j, RealType* rho_j_at_i){ + RealType cut = 0.0; - if (!initialized_) initialize(); - - AtomType* atype1 = EAMlist[*atid1]; - AtomType* atype2 = EAMlist[*atid2]; - - calcDensity(atype1, atype2, *rij, *rho_i_at_j, *rho_j_at_i); + map::iterator it; - return; - } + it = EAMMap.find(at1); + if (it != EAMMap.end()) { + EAMAtomData data1 = (*it).second; + cut = data1.rcut; + } - void EAM::calc_eam_preforce_Frho(int *atid1, RealType *rho, RealType *frho, - RealType *dfrhodrho) { + it = EAMMap.find(at2); + if (it != EAMMap.end()) { + EAMAtomData data2 = (*it).second; + if (data2.rcut > cut) + cut = data2.rcut; + } - if (!initialized_) initialize(); - - AtomType* atype1 = EAMlist[*atid1]; - - calcFunctional(atype1, *rho, *frho, *dfrhodrho); - - return; + return cut; } - RealType EAM::getEAMcut(int *atid1) { - - if (!initialized_) initialize(); - - AtomType* atype1 = EAMlist[*atid1]; - - return getRcut(atype1); - } - - 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) { - - if (!initialized_) initialize(); - - AtomType* atype1 = EAMlist[*atid1]; - AtomType* atype2 = EAMlist[*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(); - - return; - } - - void EAM::setCutoffEAM(RealType *thisRcut) { - eamRcut_ = *thisRcut; - } } -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) -#define fortranGetEAMcut FC_FUNC(geteamcut, GETEAMCUT) - - - void fortranCalcDensity(int *atid1, int *atid2, RealType *rij, - RealType *rho_i_at_j, RealType *rho_j_at_i) { - - return OpenMD::EAM::Instance()->calc_eam_prepair_rho(atid1, atid2, rij, - rho_i_at_j, - rho_j_at_i); - } - void fortranCalcFunctional(int *atid1, RealType *rho, RealType *frho, - RealType *dfrhodrho) { - - return OpenMD::EAM::Instance()->calc_eam_preforce_Frho(atid1, rho, frho, - dfrhodrho); - - } - void fortranSetCutoffEAM(RealType *rcut) { - return OpenMD::EAM::Instance()->setCutoffEAM(rcut); - } - void fortranCalcForce(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); - } - RealType fortranGetEAMcut(int* atid) { - return OpenMD::EAM::Instance()->getEAMcut(atid); - } - -}