--- branches/development/src/nonbonded/SC.cpp 2010/10/02 19:53:32 1502 +++ branches/development/src/nonbonded/SC.cpp 2012/05/18 21:44:02 1710 @@ -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 @@ -53,78 +54,27 @@ namespace OpenMD { SC::SC() : name_("SC"), initialized_(false), forceField_(NULL), scRcut_(0.0), np_(3000) {} - SCParam SC::getSCParam(AtomType* atomType) { - - // Do sanity checking on the AtomType we were passed before - // building any data structures: - if (!atomType->isSC()) { - sprintf( painCave.errMsg, - "SC::getSCParam was passed an atomType (%s) that does not\n" - "\tappear to be a Sutton-Chen (SC) atom.\n", - atomType->getName().c_str()); - painCave.severity = OPENMD_ERROR; - painCave.isFatal = 1; - simError(); - } - - GenericData* data = atomType->getPropertyByName("SC"); - if (data == NULL) { - sprintf( painCave.errMsg, "SC::getSCParam could not find SC\n" - "\tparameters for atomType %s.\n", - atomType->getName().c_str()); - painCave.severity = OPENMD_ERROR; - painCave.isFatal = 1; - simError(); - } - - SCParamGenericData* scData = dynamic_cast(data); - if (scData == NULL) { - sprintf( painCave.errMsg, - "SC::getSCParam could not convert GenericData to SCParamGenericData\n" - "\tfor atom type %s\n", atomType->getName().c_str()); - painCave.severity = OPENMD_ERROR; - painCave.isFatal = 1; - simError(); - } - - return scData->getData(); - } - - RealType SC::getC(AtomType* atomType) { - SCParam scParam = getSCParam(atomType); - return scParam.c; - } - - RealType SC::getM(AtomType* atomType) { - SCParam scParam = getSCParam(atomType); - return scParam.m; - } - RealType SC::getM(AtomType* atomType1, AtomType* atomType2) { - RealType m1 = getM(atomType1); - RealType m2 = getM(atomType2); + SuttonChenAdapter sca1 = SuttonChenAdapter(atomType1); + SuttonChenAdapter sca2 = SuttonChenAdapter(atomType2); + RealType m1 = sca1.getM(); + RealType m2 = sca2.getM(); return 0.5 * (m1 + m2); } - RealType SC::getN(AtomType* atomType) { - SCParam scParam = getSCParam(atomType); - return scParam.n; - } - RealType SC::getN(AtomType* atomType1, AtomType* atomType2) { - RealType n1 = getN(atomType1); - RealType n2 = getN(atomType2); + SuttonChenAdapter sca1 = SuttonChenAdapter(atomType1); + SuttonChenAdapter sca2 = SuttonChenAdapter(atomType2); + RealType n1 = sca1.getN(); + RealType n2 = sca2.getN(); return 0.5 * (n1 + n2); } - RealType SC::getAlpha(AtomType* atomType) { - SCParam scParam = getSCParam(atomType); - return scParam.alpha; - } - RealType SC::getAlpha(AtomType* atomType1, AtomType* atomType2) { - RealType alpha1 = getAlpha(atomType1); - RealType alpha2 = getAlpha(atomType2); + SuttonChenAdapter sca1 = SuttonChenAdapter(atomType1); + SuttonChenAdapter sca2 = SuttonChenAdapter(atomType2); + RealType alpha1 = sca1.getAlpha(); + RealType alpha2 = sca2.getAlpha(); ForceFieldOptions& fopts = forceField_->getForceFieldOptions(); std::string DistanceMix = fopts.getDistanceMixingRule(); @@ -136,14 +86,11 @@ namespace OpenMD { return 0.5 * (alpha1 + alpha2); } - RealType SC::getEpsilon(AtomType* atomType) { - SCParam scParam = getSCParam(atomType); - return scParam.epsilon; - } - - RealType SC::getEpsilon(AtomType* atomType1, AtomType* atomType2) { - RealType epsilon1 = getEpsilon(atomType1); - RealType epsilon2 = getEpsilon(atomType2); + RealType SC::getEpsilon(AtomType* atomType1, AtomType* atomType2) { + SuttonChenAdapter sca1 = SuttonChenAdapter(atomType1); + SuttonChenAdapter sca2 = SuttonChenAdapter(atomType2); + RealType epsilon1 = sca1.getEpsilon(); + RealType epsilon2 = sca2.getEpsilon(); return sqrt(epsilon1 * epsilon2); } @@ -155,7 +102,8 @@ namespace OpenMD { for (at = atomTypes->beginType(i); at != NULL; at = atomTypes->nextType(i)) { - if (at->isSC()) + SuttonChenAdapter sca = SuttonChenAdapter(at); + if (sca.isSuttonChen()) addType(at); } initialized_ = true; @@ -165,24 +113,24 @@ namespace OpenMD { void SC::addType(AtomType* atomType){ + SuttonChenAdapter sca = SuttonChenAdapter(atomType); SCAtomData scAtomData; - scAtomData.c = getC(atomType); - scAtomData.m = getM(atomType); - scAtomData.n = getN(atomType); - scAtomData.alpha = getAlpha(atomType); - scAtomData.epsilon = getEpsilon(atomType); + scAtomData.c = sca.getC(); + scAtomData.m = sca.getM(); + scAtomData.n = sca.getN(); + scAtomData.alpha = sca.getAlpha(); + scAtomData.epsilon = sca.getEpsilon(); scAtomData.rCut = 2.0 * scAtomData.alpha; // add it to the map: - AtomTypeProperties atp = atomType->getATP(); pair::iterator,bool> ret; - ret = SClist.insert( pair(atp.ident, atomType) ); + ret = SClist.insert( pair(atomType->getIdent(), atomType) ); if (ret.second == false) { sprintf( painCave.errMsg, "SC already had a previous entry with ident %d\n", - atp.ident); + atomType->getIdent() ); painCave.severity = OPENMD_INFO; painCave.isFatal = 0; simError(); @@ -302,87 +250,102 @@ namespace OpenMD { return; } - void SC::calcDensity(DensityData ddat) { + void SC::calcDensity(InteractionData &idat) { if (!initialized_) initialize(); - SCInteractionData mixer = MixingMap[make_pair(ddat.atype1, ddat.atype2)]; + SCInteractionData mixer = MixingMap[ idat.atypes ]; RealType rcij = mixer.rCut; - if (ddat.rij < rcij) { - ddat.rho_i_at_j = mixer.phi->getValueAt(ddat.rij); - ddat.rho_j_at_i = ddat.rho_i_at_j; - } else { - ddat.rho_i_at_j = 0.0; - ddat.rho_j_at_i = 0.0; - } - + if ( *(idat.rij) < rcij) { + RealType rho = mixer.phi->getValueAt( *(idat.rij) ); + *(idat.rho1) += rho; + *(idat.rho2) += rho; + } + return; } - void SC::calcFunctional(FunctionalData fdat) { + void SC::calcFunctional(SelfData &sdat) { if (!initialized_) initialize(); - SCAtomData data1 = SCMap[fdat.atype]; + SCAtomData data1 = SCMap[sdat.atype]; + + RealType u = - data1.c * data1.epsilon * sqrt( *(sdat.rho) ); + *(sdat.frho) = u; + *(sdat.dfrhodrho) = 0.5 * *(sdat.frho) / *(sdat.rho); + + (*(sdat.pot))[METALLIC_FAMILY] += u; + *(sdat.particlePot) += u; - fdat.frho = - data1.c * data1.epsilon * sqrt(fdat.rho); - fdat.dfrhodrho = 0.5 * fdat.frho / fdat.rho; - return; } - void SC::calcForce(InteractionData idat) { + void SC::calcForce(InteractionData &idat) { if (!initialized_) initialize(); - SCAtomData data1 = SCMap[idat.atype1]; - SCAtomData data2 = SCMap[idat.atype2]; + SCAtomData data1 = SCMap[idat.atypes.first]; + SCAtomData data2 = SCMap[idat.atypes.second]; - SCInteractionData mixer = MixingMap[make_pair(idat.atype1, idat.atype2)]; + SCInteractionData mixer = MixingMap[idat.atypes]; RealType rcij = mixer.rCut; - if (idat.rij < rcij) { + if ( *(idat.rij) < rcij) { RealType vcij = mixer.vCut; pair res; - res = mixer.phi->getValueAndDerivativeAt(idat.rij); + res = mixer.phi->getValueAndDerivativeAt( *(idat.rij) ); RealType rhtmp = res.first; RealType drhodr = res.second; - res = mixer.V->getValueAndDerivativeAt(idat.rij); + res = mixer.V->getValueAndDerivativeAt( *(idat.rij) ); RealType vptmp = res.first; RealType dvpdr = res.second; RealType pot_temp = vptmp - vcij; - idat.vpair += pot_temp; + *(idat.vpair) += pot_temp; - RealType dudr = drhodr * (idat.dfrho1 + idat.dfrho2) + dvpdr; + RealType dudr = drhodr * ( *(idat.dfrho1) + *(idat.dfrho2) ) + dvpdr; - idat.f1 += idat.d * dudr / idat.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 + // 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. + // 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.c * data2.epsilon * + sqrt( *(idat.rho2) - rhtmp) + *(idat.frho2); + + *(idat.particlePot2) -= data1.c * data1.epsilon * + sqrt( *(idat.rho1) - rhtmp) + *(idat.frho1); - idat.fshift1 = - data1.c * data1.epsilon * sqrt(idat.rho1 - rhtmp); - idat.fshift2 = - data2.c * data2.epsilon * sqrt(idat.rho2 - rhtmp); - - idat.pot += pot_temp; + (*(idat.pot))[METALLIC_FAMILY] += pot_temp; } return; } + + RealType SC::getSuggestedCutoffRadius(pair atypes) { + if (!initialized_) initialize(); + + map, SCInteractionData>::iterator it; + it = MixingMap.find(atypes); + if (it == MixingMap.end()) + return 0.0; + else { + SCInteractionData mixer = (*it).second; + return mixer.rCut; + } + } }