--- branches/development/src/nonbonded/Electrostatic.cpp 2010/10/02 19:53:32 1502 +++ branches/development/src/nonbonded/Electrostatic.cpp 2011/11/22 20:38:56 1665 @@ -34,9 +34,10 @@ * work. Good starting points are: * * [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005). - * [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006). + * [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 @@ -47,14 +48,34 @@ #include "utils/simError.h" #include "types/NonBondedInteractionType.hpp" #include "types/DirectionalAtomType.hpp" +#include "io/Globals.hpp" - namespace OpenMD { Electrostatic::Electrostatic(): name_("Electrostatic"), initialized_(false), - forceField_(NULL) {} + forceField_(NULL), info_(NULL), + haveCutoffRadius_(false), + haveDampingAlpha_(false), + haveDielectric_(false), + haveElectroSpline_(false) + {} void Electrostatic::initialize() { + + Globals* simParams_ = info_->getSimParams(); + + summationMap_["HARD"] = esm_HARD; + summationMap_["NONE"] = esm_HARD; + summationMap_["SWITCHING_FUNCTION"] = esm_SWITCHING_FUNCTION; + summationMap_["SHIFTED_POTENTIAL"] = esm_SHIFTED_POTENTIAL; + summationMap_["SHIFTED_FORCE"] = esm_SHIFTED_FORCE; + summationMap_["REACTION_FIELD"] = esm_REACTION_FIELD; + summationMap_["EWALD_FULL"] = esm_EWALD_FULL; + summationMap_["EWALD_PME"] = esm_EWALD_PME; + summationMap_["EWALD_SPME"] = esm_EWALD_SPME; + screeningMap_["DAMPED"] = DAMPED; + screeningMap_["UNDAMPED"] = UNDAMPED; + // these prefactors convert the multipole interactions into kcal / mol // all were computed assuming distances are measured in angstroms // Charge-Charge, assuming charges are measured in electrons @@ -79,20 +100,111 @@ namespace OpenMD { // variables to handle different summation methods for long-range // electrostatics: - summationMethod_ = NONE; + summationMethod_ = esm_HARD; screeningMethod_ = UNDAMPED; dielectric_ = 1.0; one_third_ = 1.0 / 3.0; - haveDefaultCutoff_ = false; - haveDampingAlpha_ = false; - haveDielectric_ = false; - haveElectroSpline_ = false; + // check the summation method: + if (simParams_->haveElectrostaticSummationMethod()) { + string myMethod = simParams_->getElectrostaticSummationMethod(); + toUpper(myMethod); + map::iterator i; + i = summationMap_.find(myMethod); + if ( i != summationMap_.end() ) { + summationMethod_ = (*i).second; + } else { + // throw error + sprintf( painCave.errMsg, + "Electrostatic::initialize: Unknown electrostaticSummationMethod.\n" + "\t(Input file specified %s .)\n" + "\telectrostaticSummationMethod must be one of: \"hard\",\n" + "\t\"shifted_potential\", \"shifted_force\", or \n" + "\t\"reaction_field\".\n", myMethod.c_str() ); + painCave.isFatal = 1; + simError(); + } + } else { + // set ElectrostaticSummationMethod to the cutoffMethod: + if (simParams_->haveCutoffMethod()){ + string myMethod = simParams_->getCutoffMethod(); + toUpper(myMethod); + map::iterator i; + i = summationMap_.find(myMethod); + if ( i != summationMap_.end() ) { + summationMethod_ = (*i).second; + } + } + } + + if (summationMethod_ == esm_REACTION_FIELD) { + if (!simParams_->haveDielectric()) { + // throw warning + sprintf( painCave.errMsg, + "SimInfo warning: dielectric was not specified in the input file\n\tfor the reaction field correction method.\n" + "\tA default value of %f will be used for the dielectric.\n", dielectric_); + painCave.isFatal = 0; + painCave.severity = OPENMD_INFO; + simError(); + } else { + dielectric_ = simParams_->getDielectric(); + } + haveDielectric_ = true; + } + + if (simParams_->haveElectrostaticScreeningMethod()) { + string myScreen = simParams_->getElectrostaticScreeningMethod(); + toUpper(myScreen); + map::iterator i; + i = screeningMap_.find(myScreen); + if ( i != screeningMap_.end()) { + screeningMethod_ = (*i).second; + } else { + sprintf( painCave.errMsg, + "SimInfo error: Unknown electrostaticScreeningMethod.\n" + "\t(Input file specified %s .)\n" + "\telectrostaticScreeningMethod must be one of: \"undamped\"\n" + "or \"damped\".\n", myScreen.c_str() ); + painCave.isFatal = 1; + simError(); + } + } + + // check to make sure a cutoff value has been set: + if (!haveCutoffRadius_) { + sprintf( painCave.errMsg, "Electrostatic::initialize has no Default " + "Cutoff value!\n"); + painCave.severity = OPENMD_ERROR; + painCave.isFatal = 1; + simError(); + } + + if (screeningMethod_ == DAMPED) { + if (!simParams_->haveDampingAlpha()) { + // first set a cutoff dependent alpha value + // we assume alpha depends linearly with rcut from 0 to 20.5 ang + dampingAlpha_ = 0.425 - cutoffRadius_* 0.02; + if (dampingAlpha_ < 0.0) dampingAlpha_ = 0.0; + + // throw warning + sprintf( painCave.errMsg, + "Electrostatic::initialize: dampingAlpha was not specified in the input file.\n" + "\tA default value of %f (1/ang) will be used for the cutoff of\n\t%f (ang).\n", + dampingAlpha_, cutoffRadius_); + painCave.severity = OPENMD_INFO; + painCave.isFatal = 0; + simError(); + } else { + dampingAlpha_ = simParams_->getDampingAlpha(); + } + haveDampingAlpha_ = true; + } + // find all of the Electrostatic atom Types: ForceField::AtomTypeContainer* atomTypes = forceField_->getAtomTypes(); ForceField::AtomTypeContainer::MapTypeIterator i; AtomType* at; - + for (at = atomTypes->beginType(i); at != NULL; at = atomTypes->nextType(i)) { @@ -100,40 +212,25 @@ namespace OpenMD { addType(at); } - // check to make sure a cutoff value has been set: - if (!haveDefaultCutoff_) { - sprintf( painCave.errMsg, "Electrostatic::initialize has no Default " - "Cutoff value!\n"); - painCave.severity = OPENMD_ERROR; - painCave.isFatal = 1; - simError(); - } - defaultCutoff2_ = defaultCutoff_ * defaultCutoff_; - rcuti_ = 1.0 / defaultCutoff_; + cutoffRadius2_ = cutoffRadius_ * cutoffRadius_; + rcuti_ = 1.0 / cutoffRadius_; rcuti2_ = rcuti_ * rcuti_; rcuti3_ = rcuti2_ * rcuti_; rcuti4_ = rcuti2_ * rcuti2_; if (screeningMethod_ == DAMPED) { - if (!haveDampingAlpha_) { - sprintf( painCave.errMsg, "Electrostatic::initialize has no " - "DampingAlpha value!\n"); - painCave.severity = OPENMD_ERROR; - painCave.isFatal = 1; - simError(); - } - + alpha2_ = dampingAlpha_ * dampingAlpha_; alpha4_ = alpha2_ * alpha2_; alpha6_ = alpha4_ * alpha2_; alpha8_ = alpha4_ * alpha4_; - constEXP_ = exp(-alpha2_ * defaultCutoff2_); + constEXP_ = exp(-alpha2_ * cutoffRadius2_); invRootPi_ = 0.56418958354775628695; alphaPi_ = 2.0 * dampingAlpha_ * invRootPi_; - c1c_ = erfc(dampingAlpha_ * defaultCutoff_) * rcuti_; + c1c_ = erfc(dampingAlpha_ * cutoffRadius_) * rcuti_; c2c_ = alphaPi_ * constEXP_ * rcuti_ + c1c_ * rcuti_; c3c_ = 2.0 * alphaPi_ * alpha2_ + 3.0 * c2c_ * rcuti_; c4c_ = 4.0 * alphaPi_ * alpha4_ + 5.0 * c3c_ * rcuti2_; @@ -148,21 +245,18 @@ namespace OpenMD { c6c_ = 9.0 * c5c_ * rcuti2_; } - if (summationMethod_ == REACTION_FIELD) { - if (haveDielectric_) { - preRF_ = (dielectric_ - 1.0) / - ((2.0 * dielectric_ + 1.0) * defaultCutoff2_ * defaultCutoff_); - preRF2_ = 2.0 * preRF_; - } else { - sprintf( painCave.errMsg, "Electrostatic::initialize has no Dielectric" - " value!\n"); - painCave.severity = OPENMD_ERROR; - painCave.isFatal = 1; - simError(); - } + if (summationMethod_ == esm_REACTION_FIELD) { + preRF_ = (dielectric_ - 1.0) / + ((2.0 * dielectric_ + 1.0) * cutoffRadius2_ * cutoffRadius_); + preRF2_ = 2.0 * preRF_; } - - RealType dx = defaultCutoff_ / RealType(np_ - 1); + + // Add a 2 angstrom safety window to deal with cutoffGroups that + // have charged atoms longer than the cutoffRadius away from each + // other. Splining may not be the best choice here. Direct calls + // to erfc might be preferrable. + + RealType dx = (cutoffRadius_ + 2.0) / RealType(np_ - 1); RealType rval; vector rvals; vector yvals; @@ -283,6 +377,12 @@ namespace OpenMD { simError(); } + // Quadrupoles in OpenMD are set as the diagonal elements + // of the diagonalized traceless quadrupole moment tensor. + // The column vectors of the unitary matrix that diagonalizes + // the quadrupole moment tensor become the eFrame (or the + // electrostatic version of the body-fixed frame. + Vector3dGenericData* v3dData = dynamic_cast(data); if (v3dData == NULL) { sprintf( painCave.errMsg, @@ -315,13 +415,15 @@ namespace OpenMD { return; } - void Electrostatic::setElectrostaticCutoffRadius( RealType theECR, - RealType theRSW ) { - defaultCutoff_ = theECR; - rrf_ = defaultCutoff_; - rt_ = theRSW; - haveDefaultCutoff_ = true; + void Electrostatic::setCutoffRadius( RealType rCut ) { + cutoffRadius_ = rCut; + rrf_ = cutoffRadius_; + haveCutoffRadius_ = true; } + + void Electrostatic::setSwitchingRadius( RealType rSwitch ) { + rt_ = rSwitch; + } void Electrostatic::setElectrostaticSummationMethod( ElectrostaticSummationMethod esm ) { summationMethod_ = esm; } @@ -337,7 +439,7 @@ namespace OpenMD { haveDielectric_ = true; } - void Electrostatic::calcForce(InteractionData idat) { + void Electrostatic::calcForce(InteractionData &idat) { // utility variables. Should clean these up and use the Vector3d and // Mat3x3d to replace as many as we can in future versions: @@ -351,12 +453,13 @@ namespace OpenMD { RealType ct_i, ct_j, ct_ij, a1; RealType riji, ri, ri2, ri3, ri4; RealType pref, vterm, epot, dudr; + RealType vpair(0.0); RealType scale, sc2; RealType pot_term, preVal, rfVal; RealType c2ri, c3ri, c4rij, cti3, ctj3, ctidotj; RealType preSw, preSwSc; RealType c1, c2, c3, c4; - RealType erfcVal, derfcVal; + RealType erfcVal(1.0), derfcVal(0.0); RealType BigR; Vector3d Q_i, Q_j; @@ -367,17 +470,23 @@ namespace OpenMD { Vector3d rhatdot2, rhatc4; Vector3d dVdr; + // variables for indirect (reaction field) interactions for excluded pairs: + RealType indirect_Pot(0.0); + RealType indirect_vpair(0.0); + Vector3d indirect_dVdr(V3Zero); + Vector3d indirect_duduz_i(V3Zero), indirect_duduz_j(V3Zero); + pair res; if (!initialized_) initialize(); - ElectrostaticAtomData data1 = ElectrostaticMap[idat.atype1]; - ElectrostaticAtomData data2 = ElectrostaticMap[idat.atype2]; + ElectrostaticAtomData data1 = ElectrostaticMap[idat.atypes.first]; + ElectrostaticAtomData data2 = ElectrostaticMap[idat.atypes.second]; // some variables we'll need independent of electrostatic type: - riji = 1.0 / idat.rij; - Vector3d rhat = idat.d * riji; + riji = 1.0 / *(idat.rij) ; + Vector3d rhat = *(idat.d) * riji; // logicals @@ -391,12 +500,16 @@ namespace OpenMD { bool j_is_SplitDipole = data2.is_SplitDipole; bool j_is_Quadrupole = data2.is_Quadrupole; - if (i_is_Charge) + if (i_is_Charge) { q_i = data1.charge; + if (idat.excluded) { + *(idat.skippedCharge2) += q_i; + } + } if (i_is_Dipole) { mu_i = data1.dipole_moment; - uz_i = idat.eFrame1.getColumn(2); + uz_i = idat.eFrame1->getColumn(2); ct_i = dot(uz_i, rhat); @@ -412,9 +525,9 @@ namespace OpenMD { qyy_i = Q_i.y(); qzz_i = Q_i.z(); - ux_i = idat.eFrame1.getColumn(0); - uy_i = idat.eFrame1.getColumn(1); - uz_i = idat.eFrame1.getColumn(2); + ux_i = idat.eFrame1->getColumn(0); + uy_i = idat.eFrame1->getColumn(1); + uz_i = idat.eFrame1->getColumn(2); cx_i = dot(ux_i, rhat); cy_i = dot(uy_i, rhat); @@ -425,12 +538,17 @@ namespace OpenMD { duduz_i = V3Zero; } - if (j_is_Charge) + if (j_is_Charge) { q_j = data2.charge; + if (idat.excluded) { + *(idat.skippedCharge1) += q_j; + } + } + if (j_is_Dipole) { mu_j = data2.dipole_moment; - uz_j = idat.eFrame2.getColumn(2); + uz_j = idat.eFrame2->getColumn(2); ct_j = dot(uz_j, rhat); @@ -446,9 +564,9 @@ namespace OpenMD { qyy_j = Q_j.y(); qzz_j = Q_j.z(); - ux_j = idat.eFrame2.getColumn(0); - uy_j = idat.eFrame2.getColumn(1); - uz_j = idat.eFrame2.getColumn(2); + ux_j = idat.eFrame2->getColumn(0); + uy_j = idat.eFrame2->getColumn(1); + uz_j = idat.eFrame2->getColumn(2); cx_j = dot(ux_j, rhat); cy_j = dot(uy_j, rhat); @@ -467,9 +585,13 @@ namespace OpenMD { if (j_is_Charge) { if (screeningMethod_ == DAMPED) { // assemble the damping variables - res = erfcSpline_->getValueAndDerivativeAt(idat.rij); - erfcVal = res.first; - derfcVal = res.second; + //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) ); + //erfcVal = res.first; + //derfcVal = res.second; + + erfcVal = erfc(dampingAlpha_ * *(idat.rij)); + derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2)); + c1 = erfcVal * riji; c2 = (-derfcVal + c1) * riji; } else { @@ -477,56 +599,76 @@ namespace OpenMD { c2 = c1 * riji; } - preVal = idat.electroMult * pre11_ * q_i * q_j; + preVal = *(idat.electroMult) * pre11_ * q_i * q_j; - if (summationMethod_ == SHIFTED_POTENTIAL) { + if (summationMethod_ == esm_SHIFTED_POTENTIAL) { vterm = preVal * (c1 - c1c_); - dudr = -idat.sw * preVal * c2; + dudr = - *(idat.sw) * preVal * c2; - } else if (summationMethod_ == SHIFTED_FORCE) { - vterm = preVal * ( c1 - c1c_ + c2c_*(idat.rij - defaultCutoff_) ); - dudr = idat.sw * preVal * (c2c_ - c2); + } else if (summationMethod_ == esm_SHIFTED_FORCE) { + vterm = preVal * ( c1 - c1c_ + c2c_*( *(idat.rij) - cutoffRadius_) ); + dudr = *(idat.sw) * preVal * (c2c_ - c2); - } else if (summationMethod_ == REACTION_FIELD) { - rfVal = idat.electroMult * preRF_ * idat.rij * idat.rij; + } else if (summationMethod_ == esm_REACTION_FIELD) { + rfVal = preRF_ * *(idat.rij) * *(idat.rij); + vterm = preVal * ( riji + rfVal ); - dudr = idat.sw * preVal * ( 2.0 * rfVal - riji ) * riji; + dudr = *(idat.sw) * preVal * ( 2.0 * rfVal - riji ) * riji; + + // if this is an excluded pair, there are still indirect + // interactions via the reaction field we must worry about: + if (idat.excluded) { + indirect_vpair += preVal * rfVal; + indirect_Pot += *(idat.sw) * preVal * rfVal; + indirect_dVdr += *(idat.sw) * preVal * 2.0 * rfVal * riji * rhat; + } + } else { - vterm = preVal * riji * erfcVal; - dudr = - idat.sw * preVal * c2; + vterm = preVal * riji * erfcVal; + dudr = - *(idat.sw) * preVal * c2; } - - idat.vpair += vterm; - epot += idat.sw * vterm; - dVdr += dudr * rhat; + vpair += vterm; + epot += *(idat.sw) * vterm; + dVdr += dudr * rhat; } if (j_is_Dipole) { // pref is used by all the possible methods - pref = idat.electroMult * pre12_ * q_i * mu_j; - preSw = idat.sw * pref; + pref = *(idat.electroMult) * pre12_ * q_i * mu_j; + preSw = *(idat.sw) * pref; - if (summationMethod_ == REACTION_FIELD) { + if (summationMethod_ == esm_REACTION_FIELD) { ri2 = riji * riji; ri3 = ri2 * riji; - vterm = - pref * ct_j * ( ri2 - preRF2_ * idat.rij ); - idat.vpair += vterm; - epot += idat.sw * vterm; + vterm = - pref * ct_j * ( ri2 - preRF2_ * *(idat.rij) ); + vpair += vterm; + epot += *(idat.sw) * vterm; dVdr += -preSw * (ri3 * (uz_j - 3.0 * ct_j * rhat) - preRF2_*uz_j); - duduz_j += -preSw * rhat * (ri2 - preRF2_ * idat.rij); + duduz_j += -preSw * rhat * (ri2 - preRF2_ * *(idat.rij) ); + // Even if we excluded this pair from direct interactions, + // we still have the reaction-field-mediated charge-dipole + // interaction: + + if (idat.excluded) { + indirect_vpair += pref * ct_j * preRF2_ * *(idat.rij); + indirect_Pot += preSw * ct_j * preRF2_ * *(idat.rij); + indirect_dVdr += preSw * preRF2_ * uz_j; + indirect_duduz_j += preSw * rhat * preRF2_ * *(idat.rij); + } + } else { // determine the inverse r used if we have split dipoles if (j_is_SplitDipole) { - BigR = sqrt(idat.r2 + 0.25 * d_j * d_j); + BigR = sqrt( *(idat.r2) + 0.25 * d_j * d_j); ri = 1.0 / BigR; - scale = idat.rij * ri; + scale = *(idat.rij) * ri; } else { ri = riji; scale = 1.0; @@ -536,9 +678,11 @@ namespace OpenMD { if (screeningMethod_ == DAMPED) { // assemble the damping variables - res = erfcSpline_->getValueAndDerivativeAt(idat.rij); - erfcVal = res.first; - derfcVal = res.second; + //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) ); + //erfcVal = res.first; + //derfcVal = res.second; + erfcVal = erfc(dampingAlpha_ * *(idat.rij)); + derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2)); c1 = erfcVal * ri; c2 = (-derfcVal + c1) * ri; c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * ri; @@ -553,8 +697,8 @@ namespace OpenMD { // calculate the potential pot_term = scale * c2; vterm = -pref * ct_j * pot_term; - idat.vpair += vterm; - epot += idat.sw * vterm; + vpair += vterm; + epot += *(idat.sw) * vterm; // calculate derivatives for forces and torques @@ -569,13 +713,15 @@ namespace OpenMD { cx2 = cx_j * cx_j; cy2 = cy_j * cy_j; cz2 = cz_j * cz_j; - pref = idat.electroMult * pre14_ * q_i * one_third_; + pref = *(idat.electroMult) * pre14_ * q_i * one_third_; if (screeningMethod_ == DAMPED) { // assemble the damping variables - res = erfcSpline_->getValueAndDerivativeAt(idat.rij); - erfcVal = res.first; - derfcVal = res.second; + //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) ); + //erfcVal = res.first; + //derfcVal = res.second; + erfcVal = erfc(dampingAlpha_ * *(idat.rij)); + derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2)); c1 = erfcVal * riji; c2 = (-derfcVal + c1) * riji; c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * riji; @@ -588,10 +734,10 @@ namespace OpenMD { } // precompute variables for convenience - preSw = idat.sw * pref; + preSw = *(idat.sw) * pref; c2ri = c2 * riji; c3ri = c3 * riji; - c4rij = c4 * idat.rij; + c4rij = c4 * *(idat.rij) ; rhatdot2 = 2.0 * rhat * c3; rhatc4 = rhat * c4rij; @@ -600,8 +746,8 @@ namespace OpenMD { qyy_j * (cy2*c3 - c2ri) + qzz_j * (cz2*c3 - c2ri) ); vterm = pref * pot_term; - idat.vpair += vterm; - epot += idat.sw * vterm; + vpair += vterm; + epot += *(idat.sw) * vterm; // calculate derivatives for the forces and torques @@ -619,29 +765,40 @@ namespace OpenMD { if (j_is_Charge) { // variables used by all the methods - pref = idat.electroMult * pre12_ * q_j * mu_i; - preSw = idat.sw * pref; + pref = *(idat.electroMult) * pre12_ * q_j * mu_i; + preSw = *(idat.sw) * pref; - if (summationMethod_ == REACTION_FIELD) { + if (summationMethod_ == esm_REACTION_FIELD) { ri2 = riji * riji; ri3 = ri2 * riji; - vterm = pref * ct_i * ( ri2 - preRF2_ * idat.rij ); - idat.vpair += vterm; - epot += idat.sw * vterm; + vterm = pref * ct_i * ( ri2 - preRF2_ * *(idat.rij) ); + vpair += vterm; + epot += *(idat.sw) * vterm; dVdr += preSw * (ri3 * (uz_i - 3.0 * ct_i * rhat) - preRF2_ * uz_i); - duduz_i += preSw * rhat * (ri2 - preRF2_ * idat.rij); + duduz_i += preSw * rhat * (ri2 - preRF2_ * *(idat.rij) ); + + // Even if we excluded this pair from direct interactions, + // we still have the reaction-field-mediated charge-dipole + // interaction: + + if (idat.excluded) { + indirect_vpair += -pref * ct_i * preRF2_ * *(idat.rij); + indirect_Pot += -preSw * ct_i * preRF2_ * *(idat.rij); + indirect_dVdr += -preSw * preRF2_ * uz_i; + indirect_duduz_i += -preSw * rhat * preRF2_ * *(idat.rij); + } } else { // determine inverse r if we are using split dipoles if (i_is_SplitDipole) { - BigR = sqrt(idat.r2 + 0.25 * d_i * d_i); + BigR = sqrt( *(idat.r2) + 0.25 * d_i * d_i); ri = 1.0 / BigR; - scale = idat.rij * ri; + scale = *(idat.rij) * ri; } else { ri = riji; scale = 1.0; @@ -651,9 +808,11 @@ namespace OpenMD { if (screeningMethod_ == DAMPED) { // assemble the damping variables - res = erfcSpline_->getValueAndDerivativeAt(idat.rij); - erfcVal = res.first; - derfcVal = res.second; + //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) ); + //erfcVal = res.first; + //derfcVal = res.second; + erfcVal = erfc(dampingAlpha_ * *(idat.rij)); + derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2)); c1 = erfcVal * ri; c2 = (-derfcVal + c1) * ri; c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * ri; @@ -668,8 +827,8 @@ namespace OpenMD { // calculate the potential pot_term = c2 * scale; vterm = pref * ct_i * pot_term; - idat.vpair += vterm; - epot += idat.sw * vterm; + vpair += vterm; + epot += *(idat.sw) * vterm; // calculate derivatives for the forces and torques dVdr += preSw * (uz_i * c2ri - ct_i * rhat * sc2 * c3); @@ -681,18 +840,18 @@ namespace OpenMD { // variables used by all methods ct_ij = dot(uz_i, uz_j); - pref = idat.electroMult * pre22_ * mu_i * mu_j; - preSw = idat.sw * pref; + pref = *(idat.electroMult) * pre22_ * mu_i * mu_j; + preSw = *(idat.sw) * pref; - if (summationMethod_ == REACTION_FIELD) { + if (summationMethod_ == esm_REACTION_FIELD) { ri2 = riji * riji; ri3 = ri2 * riji; ri4 = ri2 * ri2; vterm = pref * ( ri3 * (ct_ij - 3.0 * ct_i * ct_j) - preRF2_ * ct_ij ); - idat.vpair += vterm; - epot += idat.sw * vterm; + vpair += vterm; + epot += *(idat.sw) * vterm; a1 = 5.0 * ct_i * ct_j - ct_ij; @@ -701,21 +860,28 @@ namespace OpenMD { duduz_i += preSw * (ri3 * (uz_j - 3.0 * ct_j * rhat) - preRF2_*uz_j); duduz_j += preSw * (ri3 * (uz_i - 3.0 * ct_i * rhat) - preRF2_*uz_i); + if (idat.excluded) { + indirect_vpair += - pref * preRF2_ * ct_ij; + indirect_Pot += - preSw * preRF2_ * ct_ij; + indirect_duduz_i += -preSw * preRF2_ * uz_j; + indirect_duduz_j += -preSw * preRF2_ * uz_i; + } + } else { if (i_is_SplitDipole) { if (j_is_SplitDipole) { - BigR = sqrt(idat.r2 + 0.25 * d_i * d_i + 0.25 * d_j * d_j); + BigR = sqrt( *(idat.r2) + 0.25 * d_i * d_i + 0.25 * d_j * d_j); } else { - BigR = sqrt(idat.r2 + 0.25 * d_i * d_i); + BigR = sqrt( *(idat.r2) + 0.25 * d_i * d_i); } ri = 1.0 / BigR; - scale = idat.rij * ri; + scale = *(idat.rij) * ri; } else { if (j_is_SplitDipole) { - BigR = sqrt(idat.r2 + 0.25 * d_j * d_j); + BigR = sqrt( *(idat.r2) + 0.25 * d_j * d_j); ri = 1.0 / BigR; - scale = idat.rij * ri; + scale = *(idat.rij) * ri; } else { ri = riji; scale = 1.0; @@ -723,9 +889,11 @@ namespace OpenMD { } if (screeningMethod_ == DAMPED) { // assemble damping variables - res = erfcSpline_->getValueAndDerivativeAt(idat.rij); - erfcVal = res.first; - derfcVal = res.second; + //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) ); + //erfcVal = res.first; + //derfcVal = res.second; + erfcVal = erfc(dampingAlpha_ * *(idat.rij)); + derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2)); c1 = erfcVal * ri; c2 = (-derfcVal + c1) * ri; c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * ri; @@ -745,13 +913,13 @@ namespace OpenMD { preSwSc = preSw * scale; c2ri = c2 * ri; c3ri = c3 * ri; - c4rij = c4 * idat.rij; + c4rij = c4 * *(idat.rij) ; // calculate the potential pot_term = (ct_ij * c2ri - ctidotj * c3); vterm = pref * pot_term; - idat.vpair += vterm; - epot += idat.sw * vterm; + vpair += vterm; + epot += *(idat.sw) * vterm; // calculate derivatives for the forces and torques dVdr += preSwSc * ( ctidotj * rhat * c4rij - @@ -770,13 +938,15 @@ namespace OpenMD { cy2 = cy_i * cy_i; cz2 = cz_i * cz_i; - pref = idat.electroMult * pre14_ * q_j * one_third_; + pref = *(idat.electroMult) * pre14_ * q_j * one_third_; if (screeningMethod_ == DAMPED) { // assemble the damping variables - res = erfcSpline_->getValueAndDerivativeAt(idat.rij); - erfcVal = res.first; - derfcVal = res.second; + //res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) ); + //erfcVal = res.first; + //derfcVal = res.second; + erfcVal = erfc(dampingAlpha_ * *(idat.rij)); + derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2)); c1 = erfcVal * riji; c2 = (-derfcVal + c1) * riji; c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * riji; @@ -789,10 +959,10 @@ namespace OpenMD { } // precompute some variables for convenience - preSw = idat.sw * pref; + preSw = *(idat.sw) * pref; c2ri = c2 * riji; c3ri = c3 * riji; - c4rij = c4 * idat.rij; + c4rij = c4 * *(idat.rij) ; rhatdot2 = 2.0 * rhat * c3; rhatc4 = rhat * c4rij; @@ -802,8 +972,8 @@ namespace OpenMD { qzz_i * (cz2 * c3 - c2ri) ); vterm = pref * pot_term; - idat.vpair += vterm; - epot += idat.sw * vterm; + vpair += vterm; + epot += *(idat.sw) * vterm; // calculate the derivatives for the forces and torques @@ -817,162 +987,87 @@ namespace OpenMD { } } - idat.pot += epot; - idat.f1 += dVdr; - if (i_is_Dipole || i_is_Quadrupole) - idat.t1 -= cross(uz_i, duduz_i); - if (i_is_Quadrupole) { - idat.t1 -= cross(ux_i, dudux_i); - idat.t1 -= cross(uy_i, duduy_i); - } + if (!idat.excluded) { + *(idat.vpair) += vpair; + (*(idat.pot))[ELECTROSTATIC_FAMILY] += epot; + *(idat.f1) += dVdr; + + if (i_is_Dipole || i_is_Quadrupole) + *(idat.t1) -= cross(uz_i, duduz_i); + if (i_is_Quadrupole) { + *(idat.t1) -= cross(ux_i, dudux_i); + *(idat.t1) -= cross(uy_i, duduy_i); + } + + if (j_is_Dipole || j_is_Quadrupole) + *(idat.t2) -= cross(uz_j, duduz_j); + if (j_is_Quadrupole) { + *(idat.t2) -= cross(uz_j, dudux_j); + *(idat.t2) -= cross(uz_j, duduy_j); + } + + } else { + + // only accumulate the forces and torques resulting from the + // indirect reaction field terms. - if (j_is_Dipole || j_is_Quadrupole) - idat.t2 -= cross(uz_j, duduz_j); - if (j_is_Quadrupole) { - idat.t2 -= cross(uz_j, dudux_j); - idat.t2 -= cross(uz_j, duduy_j); + *(idat.vpair) += indirect_vpair; + (*(idat.pot))[ELECTROSTATIC_FAMILY] += indirect_Pot; + *(idat.f1) += indirect_dVdr; + + if (i_is_Dipole) + *(idat.t1) -= cross(uz_i, indirect_duduz_i); + if (j_is_Dipole) + *(idat.t2) -= cross(uz_j, indirect_duduz_j); } + return; } - - void Electrostatic::calcSkipCorrection(SkipCorrectionData skdat) { - - if (!initialized_) initialize(); - ElectrostaticAtomData data1 = ElectrostaticMap[skdat.atype1]; - ElectrostaticAtomData data2 = ElectrostaticMap[skdat.atype2]; - - // logicals - - bool i_is_Charge = data1.is_Charge; - bool i_is_Dipole = data1.is_Dipole; - - bool j_is_Charge = data2.is_Charge; - bool j_is_Dipole = data2.is_Dipole; - - RealType q_i, q_j; - - // The skippedCharge computation is needed by the real-space cutoff methods - // (i.e. shifted force and shifted potential) - - if (i_is_Charge) { - q_i = data1.charge; - skdat.skippedCharge2 += q_i; - } - - if (j_is_Charge) { - q_j = data2.charge; - skdat.skippedCharge1 += q_j; - } - - // the rest of this function should only be necessary for reaction field. - - if (summationMethod_ == REACTION_FIELD) { - RealType riji, ri2, ri3; - RealType q_i, mu_i, ct_i; - RealType q_j, mu_j, ct_j; - RealType preVal, rfVal, vterm, dudr, pref, myPot; - Vector3d dVdr, uz_i, uz_j, duduz_i, duduz_j, rhat; - - // some variables we'll need independent of electrostatic type: - - riji = 1.0 / skdat.rij; - rhat = skdat.d * riji; - - if (i_is_Dipole) { - mu_i = data1.dipole_moment; - uz_i = skdat.eFrame1->getColumn(2); - ct_i = dot(uz_i, rhat); - duduz_i = V3Zero; - } - - if (j_is_Dipole) { - mu_j = data2.dipole_moment; - uz_j = skdat.eFrame2->getColumn(2); - ct_j = dot(uz_j, rhat); - duduz_j = V3Zero; - } - - if (i_is_Charge) { - if (j_is_Charge) { - preVal = skdat.electroMult * pre11_ * q_i * q_j; - rfVal = preRF_ * skdat.rij * skdat.rij; - vterm = preVal * rfVal; - myPot += skdat.sw * vterm; - dudr = skdat.sw * preVal * 2.0 * rfVal * riji; - dVdr += dudr * rhat; - } - - if (j_is_Dipole) { - ri2 = riji * riji; - ri3 = ri2 * riji; - pref = skdat.electroMult * pre12_ * q_i * mu_j; - vterm = - pref * ct_j * ( ri2 - preRF2_ * skdat.rij ); - myPot += skdat.sw * vterm; - dVdr += -skdat.sw * pref * ( ri3 * ( uz_j - 3.0 * ct_j * rhat) - preRF2_ * uz_j); - duduz_j += -skdat.sw * pref * rhat * (ri2 - preRF2_ * skdat.rij); - } - } - if (i_is_Dipole) { - if (j_is_Charge) { - ri2 = riji * riji; - ri3 = ri2 * riji; - pref = skdat.electroMult * pre12_ * q_j * mu_i; - vterm = - pref * ct_i * ( ri2 - preRF2_ * skdat.rij ); - myPot += skdat.sw * vterm; - dVdr += skdat.sw * pref * ( ri3 * ( uz_i - 3.0 * ct_i * rhat) - preRF2_ * uz_i); - duduz_i += skdat.sw * pref * rhat * (ri2 - preRF2_ * skdat.rij); - } - } - - // accumulate the forces and torques resulting from the self term - skdat.pot += myPot; - skdat.f1 += dVdr; - - if (i_is_Dipole) - *(skdat.t1) -= cross(uz_i, duduz_i); - if (j_is_Dipole) - *(skdat.t2) -= cross(uz_j, duduz_j); - } - } - - void Electrostatic::calcSelfCorrection(SelfCorrectionData scdat) { + void Electrostatic::calcSelfCorrection(SelfData &sdat) { RealType mu1, preVal, chg1, self; if (!initialized_) initialize(); - - ElectrostaticAtomData data = ElectrostaticMap[scdat.atype]; + + ElectrostaticAtomData data = ElectrostaticMap[sdat.atype]; // logicals - bool i_is_Charge = data.is_Charge; bool i_is_Dipole = data.is_Dipole; - if (summationMethod_ == REACTION_FIELD) { + if (summationMethod_ == esm_REACTION_FIELD) { if (i_is_Dipole) { mu1 = data.dipole_moment; preVal = pre22_ * preRF2_ * mu1 * mu1; - scdat.pot -= 0.5 * preVal; + (*(sdat.pot))[ELECTROSTATIC_FAMILY] -= 0.5 * preVal; // The self-correction term adds into the reaction field vector - Vector3d uz_i = scdat.eFrame->getColumn(2); + Vector3d uz_i = sdat.eFrame->getColumn(2); Vector3d ei = preVal * uz_i; // This looks very wrong. A vector crossed with itself is zero. - *(scdat.t) -= cross(uz_i, ei); + *(sdat.t) -= cross(uz_i, ei); } - } else if (summationMethod_ == SHIFTED_FORCE || summationMethod_ == SHIFTED_POTENTIAL) { + } else if (summationMethod_ == esm_SHIFTED_FORCE || summationMethod_ == esm_SHIFTED_POTENTIAL) { if (i_is_Charge) { chg1 = data.charge; if (screeningMethod_ == DAMPED) { - self = - 0.5 * (c1c_ + alphaPi_) * chg1 * (chg1 + scdat.skippedCharge) * pre11_; + self = - 0.5 * (c1c_ + alphaPi_) * chg1 * (chg1 + *(sdat.skippedCharge)) * pre11_; } else { - self = - 0.5 * rcuti_ * chg1 * (chg1 + scdat.skippedCharge) * pre11_; + self = - 0.5 * rcuti_ * chg1 * (chg1 + *(sdat.skippedCharge)) * pre11_; } - scdat.pot += self; + (*(sdat.pot))[ELECTROSTATIC_FAMILY] += self; } } } + + RealType Electrostatic::getSuggestedCutoffRadius(pair atypes) { + // This seems to work moderately well as a default. There's no + // inherent scale for 1/r interactions that we can standardize. + // 12 angstroms seems to be a reasonably good guess for most + // cases. + return 12.0; + } }