--- branches/development/src/nonbonded/SC.cpp 2011/11/22 20:38:56 1665 +++ trunk/src/nonbonded/SC.cpp 2013/07/01 21:09:37 1895 @@ -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,78 +54,36 @@ 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(); - } + SC::~SC() { + initialized_ = false; - RealType SC::getC(AtomType* atomType) { - SCParam scParam = getSCParam(atomType); - return scParam.c; + MixingMap.clear(); + SCtypes.clear(); + SCdata.clear(); + SCtids.clear(); } - - 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(); @@ -137,28 +95,33 @@ 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); } - void SC::initialize() { + void SC::initialize() { // find all of the SC atom Types: - ForceField::AtomTypeContainer* atomTypes = forceField_->getAtomTypes(); - ForceField::AtomTypeContainer::MapTypeIterator i; - AtomType* at; + SCtypes.clear(); + SCtids.clear(); + SCdata.clear(); + MixingMap.clear(); + nSC_ = 0; - for (at = atomTypes->beginType(i); at != NULL; - at = atomTypes->nextType(i)) { - if (at->isSC()) - addType(at); - } + SCtids.resize( forceField_->getNAtomType(), -1); + + set::iterator at; + for (at = simTypes_.begin(); at != simTypes_.end(); ++at) { + if ((*at)->isSC()) nSC_++; + } + SCdata.resize(nSC_); + MixingMap.resize(nSC_); + for (at = simTypes_.begin(); at != simTypes_.end(); ++at) { + if ((*at)->isSC()) addType((*at)); + } initialized_ = true; } @@ -166,37 +129,42 @@ 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(); + + // add it to the map: + int atid = atomType->getIdent(); + int sctid = SCtypes.size(); - pair::iterator,bool> ret; - ret = SClist.insert( pair(atp.ident, atomType) ); + pair::iterator,bool> ret; + ret = SCtypes.insert( atid ); if (ret.second == false) { sprintf( painCave.errMsg, "SC already had a previous entry with ident %d\n", - atp.ident); + atid ); painCave.severity = OPENMD_INFO; painCave.isFatal = 0; simError(); } - - SCMap[atomType] = scAtomData; + SCtids[atid] = sctid; + SCdata[sctid] = scAtomData; + MixingMap[sctid].resize(nSC_); + // Now, iterate over all known types and add to the mixing map: - map::iterator it; - for( it = SCMap.begin(); it != SCMap.end(); ++it) { + std::set::iterator it; + for( it = SCtypes.begin(); it != SCtypes.end(); ++it) { - AtomType* atype2 = (*it).first; + int sctid2 = SCtids[ (*it) ]; + AtomType* atype2 = forceField_->getAtomType( (*it) ); SCInteractionData mixer; @@ -235,13 +203,11 @@ namespace OpenMD { mixer.explicitlySet = false; - pair key1, key2; - key1 = make_pair(atomType, atype2); - key2 = make_pair(atype2, atomType); + MixingMap[sctid2].resize( nSC_ ); - MixingMap[key1] = mixer; - if (key2 != key1) { - MixingMap[key2] = mixer; + MixingMap[sctid][sctid2] = mixer; + if (sctid2 != sctid) { + MixingMap[sctid2][sctid] = mixer; } } return; @@ -292,13 +258,12 @@ namespace OpenMD { mixer.explicitlySet = true; - pair key1, key2; - key1 = make_pair(atype1, atype2); - key2 = make_pair(atype2, atype1); - - MixingMap[key1] = mixer; - if (key2 != key1) { - MixingMap[key2] = mixer; + int sctid1 = SCtids[ atype1->getIdent() ]; + int sctid2 = SCtids[ atype2->getIdent() ]; + + MixingMap[sctid1][sctid2] = mixer; + if (sctid2 != sctid1) { + MixingMap[sctid2][sctid1] = mixer; } return; } @@ -306,8 +271,10 @@ namespace OpenMD { void SC::calcDensity(InteractionData &idat) { if (!initialized_) initialize(); + int sctid1 = SCtids[idat.atid1]; + int sctid2 = SCtids[idat.atid2]; - SCInteractionData mixer = MixingMap[ idat.atypes ]; + SCInteractionData &mixer = MixingMap[sctid1][sctid2]; RealType rcij = mixer.rCut; @@ -324,15 +291,17 @@ namespace OpenMD { if (!initialized_) initialize(); - SCAtomData data1 = SCMap[sdat.atype]; + SCAtomData &data1 = SCdata[SCtids[sdat.atid]]; 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; - + if (sdat.doParticlePot) { + *(sdat.particlePot) += u; + } + return; } @@ -341,26 +310,23 @@ namespace OpenMD { if (!initialized_) initialize(); - SCAtomData data1 = SCMap[idat.atypes.first]; - SCAtomData data2 = SCMap[idat.atypes.second]; + int &sctid1 = SCtids[idat.atid1]; + int &sctid2 = SCtids[idat.atid2]; - SCInteractionData mixer = MixingMap[idat.atypes]; + SCAtomData &data1 = SCdata[sctid1]; + SCAtomData &data2 = SCdata[sctid2]; + SCInteractionData &mixer = MixingMap[sctid1][sctid2]; + RealType rcij = mixer.rCut; if ( *(idat.rij) < rcij) { - RealType vcij = mixer.vCut; + RealType vcij = mixer.vCut; + RealType rhtmp, drhodr, vptmp, dvpdr; - pair res; + mixer.phi->getValueAndDerivativeAt( *(idat.rij), rhtmp, drhodr ); + mixer.V->getValueAndDerivativeAt( *(idat.rij), vptmp, dvpdr); - 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; @@ -368,20 +334,22 @@ namespace OpenMD { *(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. - - *(idat.particlePot1) -= data2.c * data2.epsilon * - sqrt( *(idat.rho2) - rhtmp) + *(idat.frho2); + 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.particlePot2) -= data1.c * data1.epsilon * + sqrt( *(idat.rho1) - rhtmp) + *(idat.frho1); + } (*(idat.pot))[METALLIC_FAMILY] += pot_temp; } @@ -392,13 +360,15 @@ namespace OpenMD { RealType SC::getSuggestedCutoffRadius(pair atypes) { if (!initialized_) initialize(); - map, SCInteractionData>::iterator it; - it = MixingMap.find(atypes); - if (it == MixingMap.end()) + int atid1 = atypes.first->getIdent(); + int atid2 = atypes.second->getIdent(); + int &sctid1 = SCtids[atid1]; + int &sctid2 = SCtids[atid2]; + + if (sctid1 == -1 || sctid2 == -1) { return 0.0; - else { - SCInteractionData mixer = (*it).second; - return mixer.rCut; + } else { + return MixingMap[sctid1][sctid2].rCut; } } }