--- branches/development/src/nonbonded/SC.cpp 2010/08/10 18:34:59 1489 +++ branches/development/src/nonbonded/SC.cpp 2012/05/19 02:58:35 1711 @@ -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 @@ -49,95 +50,31 @@ namespace OpenMD { namespace OpenMD { - bool SC::initialized_ = false; - RealType SC::scRcut_ = 0.0; - int SC::np_ = 3000; - ForceField* SC::forceField_ = NULL; - map SC::SClist; - map SC::SCMap; - map, SCInteractionData> SC::MixingMap; - - SC* SC::_instance = NULL; - SC* SC::Instance() { - if (!_instance) { - _instance = new SC(); - } - return _instance; - } + 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(); @@ -149,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); } @@ -168,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; @@ -178,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(); @@ -315,191 +250,106 @@ namespace OpenMD { return; } - void SC::calcDensity(AtomType* at1, AtomType* at2, const RealType rij, - RealType &rho_i_at_j, RealType &rho_j_at_i) { + void SC::calcDensity(InteractionData &idat) { if (!initialized_) initialize(); - SCInteractionData mixer = MixingMap[make_pair(at1, at2)]; + SCInteractionData mixer = MixingMap[ idat.atypes ]; - rho_i_at_j = mixer.phi->getValueAt(rij); - rho_j_at_i = rho_i_at_j; + RealType rcij = mixer.rCut; + if ( *(idat.rij) < rcij) { + RealType rho = mixer.phi->getValueAt( *(idat.rij) ); + *(idat.rho1) += rho; + *(idat.rho2) += rho; + } + return; } - void SC::calcFunctional(AtomType* at1, RealType rho, RealType &frho, - RealType &dfrhodrho) { + void SC::calcFunctional(SelfData &sdat) { if (!initialized_) initialize(); - SCAtomData data1 = SCMap[at1]; - - frho = - data1.c * data1.epsilon * sqrt(rho); - dfrhodrho = 0.5 * frho / rho; - + 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; + if (sdat.doParticlePot) { + *(sdat.particlePot) += u; + } + return; } - + - void SC::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 SC::calcForce(InteractionData &idat) { if (!initialized_) initialize(); - SCAtomData data1 = SCMap[at1]; - SCAtomData data2 = SCMap[at1]; + SCAtomData data1 = SCMap[idat.atypes.first]; + SCAtomData data2 = SCMap[idat.atypes.second]; - SCInteractionData mixer = MixingMap[make_pair(at1, at2)]; + SCInteractionData mixer = MixingMap[idat.atypes]; RealType rcij = mixer.rCut; - RealType vcij = mixer.vCut; - pair res; - - res = mixer.phi->getValueAndDerivativeAt(rij); - RealType rhtmp = res.first; - RealType drhodr = res.second; - - res = mixer.V->getValueAndDerivativeAt(rij); - RealType vptmp = res.first; - RealType dvpdr = res.second; - - RealType pot_temp = vptmp - vcij; - vpair += pot_temp; - - RealType dudr = drhodr * (dfrhodrho_i + dfrhodrho_j) + dvpdr; - - f1 += d * dudr / rij; + if ( *(idat.rij) < rcij) { + RealType vcij = mixer.vCut; + + pair res; + + res = mixer.phi->getValueAndDerivativeAt( *(idat.rij) ); + RealType rhtmp = res.first; + RealType drhodr = res.second; + + res = mixer.V->getValueAndDerivativeAt( *(idat.rij) ); + RealType vptmp = res.first; + RealType dvpdr = res.second; + + RealType pot_temp = vptmp - vcij; + *(idat.vpair) += pot_temp; + + RealType dudr = drhodr * ( *(idat.dfrho1) + *(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 - // 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.c * data1.epsilon * sqrt(rho_i - rhtmp); - fshift_j = - data2.c * data2.epsilon * sqrt(rho_j - rhtmp); - - pot += pot_temp; + 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.pot))[METALLIC_FAMILY] += pot_temp; + } + return; } + RealType SC::getSuggestedCutoffRadius(pair atypes) { + if (!initialized_) initialize(); - void SC::calc_sc_prepair_rho(int *atid1, int *atid2, RealType *rij, - RealType* rho_i_at_j, RealType* rho_j_at_i){ - - if (!initialized_) initialize(); - - AtomType* atype1 = SClist[*atid1]; - AtomType* atype2 = SClist[*atid2]; - - calcDensity(atype1, atype2, *rij, *rho_i_at_j, *rho_j_at_i); - - return; + map, SCInteractionData>::iterator it; + it = MixingMap.find(atypes); + if (it == MixingMap.end()) + return 0.0; + else { + SCInteractionData mixer = (*it).second; + return mixer.rCut; + } } - - void SC::calc_sc_preforce_Frho(int *atid1, RealType *rho, RealType *frho, - RealType *dfrhodrho) { - - if (!initialized_) initialize(); - - AtomType* atype1 = SClist[*atid1]; - - calcFunctional(atype1, *rho, *frho, *dfrhodrho); - - return; - } - - RealType SC::getSCcut(int *atid1) { - - if (!initialized_) initialize(); - - AtomType* atype1 = SClist[*atid1]; - - return 2.0 * getAlpha(atype1); - } - - void SC::do_sc_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 = SClist[*atid1]; - AtomType* atype2 = SClist[*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 SC::setCutoffSC(RealType *thisRcut) { - scRcut_ = *thisRcut; - } } - -extern "C" { - -#define fortranCalcDensity FC_FUNC(calc_sc_prepair_rho, CALC_SC_PREPAIR_RHO) -#define fortranCalcFunctional FC_FUNC(calc_sc_preforce_frho, CALC_SC_PREFORCE_FRHO) -#define fortranCalcForce FC_FUNC(do_sc_pair, DO_SC_PAIR) -#define fortranSetCutoffSC FC_FUNC(setcutoffsc, SETCUTOFFSC) -#define fortranGetSCcut FC_FUNC(getsccut, GETSCCUT) - - - void fortranCalcDensity(int *atid1, int *atid2, RealType *rij, - RealType *rho_i_at_j, RealType *rho_j_at_i) { - - return OpenMD::SC::Instance()->calc_sc_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::SC::Instance()->calc_sc_preforce_Frho(atid1, rho, frho, - dfrhodrho); - - } - void fortranSetCutoffSC(RealType *rcut) { - return OpenMD::SC::Instance()->setCutoffSC(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::SC::Instance()->do_sc_pair(atid1, atid2, d, rij, - r2, sw, vpair, - pot, f1, rho1, - rho2, dfrho1, dfrho2, - fshift1, fshift2); - } - RealType fortranGetSCcut(int* atid) { - return OpenMD::SC::Instance()->getSCcut(atid); - } - -}