--- branches/development/src/nonbonded/SC.cpp 2011/05/27 16:45:44 1571 +++ branches/development/src/nonbonded/SC.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 @@ -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(); @@ -134,16 +84,13 @@ namespace OpenMD { return sqrt(alpha1 * alpha2); else 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(); @@ -311,12 +259,10 @@ namespace OpenMD { RealType rcij = mixer.rCut; if ( *(idat.rij) < rcij) { - *(idat.rho_i_at_j) = mixer.phi->getValueAt( *(idat.rij) ); - *(idat.rho_j_at_i) = *(idat.rho_i_at_j); - } else { - *(idat.rho_i_at_j) = 0.0; - *(idat.rho_j_at_i) = 0.0; - } + RealType rho = mixer.phi->getValueAt( *(idat.rij) ); + *(idat.rho1) += rho; + *(idat.rho2) += rho; + } return; } @@ -326,10 +272,16 @@ namespace OpenMD { if (!initialized_) initialize(); SCAtomData data1 = SCMap[sdat.atype]; - - *(sdat.frho) = - data1.c * data1.epsilon * sqrt( *(sdat.rho) ); + + 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; + if (sdat.doParticlePot) { + *(sdat.particlePot) += u; + } + return; } @@ -365,22 +317,24 @@ namespace OpenMD { *(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 - // 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. + 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.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[METALLIC_FAMILY] += pot_temp; + (*(idat.pot))[METALLIC_FAMILY] += pot_temp; } return;