--- branches/development/src/nonbonded/SC.cpp 2010/10/02 19:53:32 1502 +++ branches/development/src/nonbonded/SC.cpp 2011/04/08 21:25:19 1545 @@ -302,46 +302,46 @@ 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; + 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 { - ddat.rho_i_at_j = 0.0; - ddat.rho_j_at_i = 0.0; + idat.rho_i_at_j = 0.0; + idat.rho_j_at_i = 0.0; } return; } - void SC::calcFunctional(FunctionalData fdat) { + void SC::calcFunctional(SelfData &sdat) { if (!initialized_) initialize(); - SCAtomData data1 = SCMap[fdat.atype]; + SCAtomData data1 = SCMap[sdat.atype]; - fdat.frho = - data1.c * data1.epsilon * sqrt(fdat.rho); - fdat.dfrhodrho = 0.5 * fdat.frho / fdat.rho; + sdat.frho = - data1.c * data1.epsilon * sqrt(sdat.rho); + sdat.dfrhodrho = 0.5 * sdat.frho / sdat.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; @@ -359,7 +359,7 @@ namespace OpenMD { RealType dvpdr = res.second; RealType pot_temp = vptmp - vcij; - idat.vpair += pot_temp; + idat.vpair[3] += pot_temp; RealType dudr = drhodr * (idat.dfrho1 + idat.dfrho2) + dvpdr; @@ -380,9 +380,22 @@ namespace OpenMD { 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[3] += 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; + } + } }