--- branches/development/src/nonbonded/EAM.cpp 2011/11/22 20:38:56 1665 +++ trunk/src/nonbonded/EAM.cpp 2014/11/01 14:12:16 2033 @@ -35,7 +35,7 @@ * * [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). + * [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). */ @@ -54,112 +54,32 @@ namespace OpenMD { EAM::EAM() : name_("EAM"), initialized_(false), forceField_(NULL), 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::getPhi(AtomType* atomType1, AtomType* atomType2) { + EAMAdapter ea1 = EAMAdapter(atomType1); + EAMAdapter ea2 = EAMAdapter(atomType2); + CubicSpline* z1 = ea1.getZ(); + CubicSpline* z2 = ea2.getZ(); - 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); + // Thise prefactors convert the charge-charge interactions into + // kcal / mol all were computed assuming distances are measured in + // angstroms Charge-Charge, assuming charges are measured in + // electrons. Matches value in Electrostatics.cpp + pre11_ = 332.0637778; // make the r grid: - // 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); @@ -175,7 +95,7 @@ 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]; // only use z(r) if we're inside this atom's cutoff radius, @@ -183,14 +103,13 @@ 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; - + phi = pre11_ * (zi * zj) / r; + phivals.push_back(phi); } - CubicSpline* cs = new CubicSpline(); cs->addPoints(rvals, phivals); return cs; @@ -202,7 +121,6 @@ namespace OpenMD { } void EAM::initialize() { - // set up the mixing method: ForceFieldOptions& fopts = forceField_->getForceFieldOptions(); string EAMMixMeth = fopts.getEAMMixingMethod(); @@ -216,16 +134,24 @@ namespace OpenMD { mixMeth_ = eamUnknown; // find all of the EAM atom Types: - ForceField::AtomTypeContainer* atomTypes = forceField_->getAtomTypes(); - ForceField::AtomTypeContainer::MapTypeIterator i; - AtomType* at; + EAMtypes.clear(); + EAMtids.clear(); + EAMdata.clear(); + MixingMap.clear(); + nEAM_ = 0; + + EAMtids.resize( forceField_->getNAtomType(), -1); - for (at = atomTypes->beginType(i); at != NULL; - at = atomTypes->nextType(i)) { - - if (at->isEAM()) - addType(at); + set::iterator at; + for (at = simTypes_.begin(); at != simTypes_.end(); ++at) { + if ((*at)->isEAM()) nEAM_++; } + EAMdata.resize(nEAM_); + MixingMap.resize(nEAM_); + + for (at = simTypes_.begin(); at != simTypes_.end(); ++at) { + if ((*at)->isEAM()) addType(*at); + } // find all of the explicit EAM interactions (setfl): ForceField::NonBondedInteractionTypeContainer* nbiTypes = forceField_->getNonBondedInteractionTypes(); @@ -278,47 +204,52 @@ 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(); + int atid = atomType->getIdent(); + int eamtid = EAMtypes.size(); - pair::iterator,bool> ret; - ret = EAMlist.insert( pair(atp.ident, atomType) ); + pair::iterator,bool> ret; + ret = EAMtypes.insert( atid ); if (ret.second == false) { sprintf( painCave.errMsg, "EAM already had a previous entry with ident %d\n", - atp.ident); + atid); painCave.severity = OPENMD_INFO; painCave.isFatal = 0; simError(); } - EAMMap[atomType] = eamAtomData; + + EAMtids[atid] = eamtid; + EAMdata[eamtid] = eamAtomData; + MixingMap[eamtid].resize(nEAM_); // Now, iterate over all known types and add to the mixing map: - map::iterator it; - for( it = EAMMap.begin(); it != EAMMap.end(); ++it) { + std::set::iterator it; + for( it = EAMtypes.begin(); it != EAMtypes.end(); ++it) { - AtomType* atype2 = (*it).first; + int eamtid2 = EAMtids[ (*it) ]; + AtomType* atype2 = forceField_->getAtomType( (*it) ); EAMInteractionData mixer; mixer.phi = getPhi(atomType, atype2); + mixer.rcut = mixer.phi->getLimits().second; mixer.explicitlySet = false; - pair key1, key2; - key1 = make_pair(atomType, atype2); - key2 = make_pair(atype2, atomType); + MixingMap[eamtid2].resize( nEAM_ ); - MixingMap[key1] = mixer; - if (key2 != key1) { - MixingMap[key2] = mixer; + MixingMap[eamtid][eamtid2] = mixer; + if (eamtid2 != eamtid) { + MixingMap[eamtid2][eamtid] = mixer; } } return; @@ -340,15 +271,15 @@ namespace OpenMD { cs->addPoints(rVals, phiVals); mixer.phi = cs; + mixer.rcut = mixer.phi->getLimits().second; mixer.explicitlySet = true; - pair key1, key2; - key1 = make_pair(atype1, atype2); - key2 = make_pair(atype2, atype1); + int eamtid1 = EAMtids[ atype1->getIdent() ]; + int eamtid2 = EAMtids[ atype2->getIdent() ]; - MixingMap[key1] = mixer; - if (key2 != key1) { - MixingMap[key2] = mixer; + MixingMap[eamtid1][eamtid2] = mixer; + if (eamtid2 != eamtid1) { + MixingMap[eamtid2][eamtid1] = mixer; } return; } @@ -357,18 +288,19 @@ namespace OpenMD { if (!initialized_) initialize(); - EAMAtomData data1 = EAMMap[idat.atypes.first]; - EAMAtomData data2 = EAMMap[idat.atypes.second]; - + EAMAtomData &data1 = EAMdata[EAMtids[idat.atid1]]; + EAMAtomData &data2 = EAMdata[EAMtids[idat.atid2]]; + if (haveCutoffRadius_) if ( *(idat.rij) > eamRcut_) return; - if ( *(idat.rij) < data1.rcut) - *(idat.rho1) += data1.rho->getValueAt( *(idat.rij)); - + if ( *(idat.rij) < data1.rcut) { + *(idat.rho2) += data1.rho->getValueAt( *(idat.rij)); + } - if ( *(idat.rij) < data2.rcut) - *(idat.rho2) += data2.rho->getValueAt( *(idat.rij)); + if ( *(idat.rij) < data2.rcut) { + *(idat.rho1) += data2.rho->getValueAt( *(idat.rij)); + } return; } @@ -376,17 +308,16 @@ namespace OpenMD { void EAM::calcFunctional(SelfData &sdat) { if (!initialized_) initialize(); + EAMAtomData &data1 = EAMdata[ EAMtids[sdat.atid] ]; + + data1.F->getValueAndDerivativeAt( *(sdat.rho), *(sdat.frho), + *(sdat.dfrhodrho) ); - EAMAtomData data1 = EAMMap[ sdat.atype ]; - - pair result = data1.F->getValueAndDerivativeAt( *(sdat.rho) ); + (*(sdat.pot))[METALLIC_FAMILY] += *(sdat.frho); + if (sdat.doParticlePot) { + *(sdat.particlePot) += *(sdat.frho); + } - *(sdat.frho) = result.first; - *(sdat.dfrhodrho) = result.second; - - (*(sdat.pot))[METALLIC_FAMILY] += result.first; - *(sdat.particlePot) += result.first; - return; } @@ -398,15 +329,17 @@ namespace OpenMD { if (haveCutoffRadius_) if ( *(idat.rij) > eamRcut_) return; - pair res; + + int eamtid1 = EAMtids[idat.atid1]; + int eamtid2 = EAMtids[idat.atid2]; + EAMAtomData &data1 = EAMdata[eamtid1]; + EAMAtomData &data2 = EAMdata[eamtid2]; - 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); @@ -414,49 +347,35 @@ namespace OpenMD { 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; + data1.rho->getValueAndDerivativeAt( *(idat.rij), rha, drha); + CubicSpline* phi = MixingMap[eamtid1][eamtid1].phi; + phi->getValueAndDerivativeAt( *(idat.rij), pha, dpha); } 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; + data2.rho->getValueAndDerivativeAt( *(idat.rij), rhb, drhb ); + CubicSpline* phi = MixingMap[eamtid2][eamtid2].phi; + phi->getValueAndDerivativeAt( *(idat.rij), phb, dphb); } - 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; - + if ( *(idat.rij) < MixingMap[eamtid1][eamtid2].rcut) { + MixingMap[eamtid1][eamtid2].phi->getValueAndDerivativeAt( *(idat.rij), + phab, dvpdr); + } break; case eamUnknown: default: @@ -476,28 +395,27 @@ namespace OpenMD { dudr = drhojdr* *(idat.dfrho1) + drhoidr* *(idat.dfrho2) + dvpdr; *(idat.f1) += *(idat.d) * dudr / *(idat.rij); - - // 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. + + 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. + + *(idat.particlePot1) += data2.F->getValueAt( *(idat.rho2) - rha ) + - *(idat.frho2); + + *(idat.particlePot2) += data1.F->getValueAt( *(idat.rho1) - rhb) + - *(idat.frho1); + } - *(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; - - return; - + (*(idat.pot))[METALLIC_FAMILY] += phab; + *(idat.vpair) += phab; + return; } RealType EAM::getSuggestedCutoffRadius(pair atypes) { @@ -505,21 +423,22 @@ namespace OpenMD { RealType cut = 0.0; - map::iterator it; - - it = EAMMap.find(atypes.first); - if (it != EAMMap.end()) { - EAMAtomData data1 = (*it).second; + int atid1 = atypes.first->getIdent(); + int atid2 = atypes.second->getIdent(); + int eamtid1 = EAMtids[atid1]; + int eamtid2 = EAMtids[atid2]; + + if (eamtid1 != -1) { + EAMAtomData data1 = EAMdata[eamtid1]; cut = data1.rcut; } - it = EAMMap.find(atypes.second); - if (it != EAMMap.end()) { - EAMAtomData data2 = (*it).second; + if (eamtid2 != -1) { + EAMAtomData data2 = EAMdata[eamtid2]; if (data2.rcut > cut) cut = data2.rcut; } - + return cut; } }