--- branches/development/src/nonbonded/Electrostatic.cpp 2011/05/27 16:45:44 1571 +++ branches/development/src/nonbonded/Electrostatic.cpp 2012/05/18 21:44:02 1710 @@ -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 @@ -46,19 +47,26 @@ #include "nonbonded/Electrostatic.hpp" #include "utils/simError.h" #include "types/NonBondedInteractionType.hpp" -#include "types/DirectionalAtomType.hpp" +#include "types/FixedChargeAdapter.hpp" +#include "types/MultipoleAdapter.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_; + + 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; @@ -97,10 +105,6 @@ namespace OpenMD { screeningMethod_ = UNDAMPED; dielectric_ = 1.0; one_third_ = 1.0 / 3.0; - haveCutoffRadius_ = false; - haveDampingAlpha_ = false; - haveDielectric_ = false; - haveElectroSpline_ = false; // check the summation method: if (simParams_->haveElectrostaticSummationMethod()) { @@ -115,7 +119,7 @@ namespace OpenMD { sprintf( painCave.errMsg, "Electrostatic::initialize: Unknown electrostaticSummationMethod.\n" "\t(Input file specified %s .)\n" - "\telectrostaticSummationMethod must be one of: \"none\",\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; @@ -248,7 +252,12 @@ namespace OpenMD { preRF2_ = 2.0 * preRF_; } - RealType dx = cutoffRadius_ / 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; @@ -272,132 +281,42 @@ namespace OpenMD { electrostaticAtomData.is_SplitDipole = false; electrostaticAtomData.is_Quadrupole = false; - if (atomType->isCharge()) { - GenericData* data = atomType->getPropertyByName("Charge"); + FixedChargeAdapter fca = FixedChargeAdapter(atomType); - if (data == NULL) { - sprintf( painCave.errMsg, "Electrostatic::addType could not find " - "Charge\n" - "\tparameters for atomType %s.\n", - atomType->getName().c_str()); - painCave.severity = OPENMD_ERROR; - painCave.isFatal = 1; - simError(); - } - - DoubleGenericData* doubleData = dynamic_cast(data); - if (doubleData == NULL) { - sprintf( painCave.errMsg, - "Electrostatic::addType could not convert GenericData to " - "Charge for\n" - "\tatom type %s\n", atomType->getName().c_str()); - painCave.severity = OPENMD_ERROR; - painCave.isFatal = 1; - simError(); - } + if (fca.isFixedCharge()) { electrostaticAtomData.is_Charge = true; - electrostaticAtomData.charge = doubleData->getData(); + electrostaticAtomData.charge = fca.getCharge(); } - if (atomType->isDirectional()) { - DirectionalAtomType* daType = dynamic_cast(atomType); - - if (daType->isDipole()) { - GenericData* data = daType->getPropertyByName("Dipole"); - - if (data == NULL) { - sprintf( painCave.errMsg, - "Electrostatic::addType could not find Dipole\n" - "\tparameters for atomType %s.\n", - daType->getName().c_str()); - painCave.severity = OPENMD_ERROR; - painCave.isFatal = 1; - simError(); - } - - DoubleGenericData* doubleData = dynamic_cast(data); - if (doubleData == NULL) { - sprintf( painCave.errMsg, - "Electrostatic::addType could not convert GenericData to " - "Dipole Moment\n" - "\tfor atom type %s\n", daType->getName().c_str()); - painCave.severity = OPENMD_ERROR; - painCave.isFatal = 1; - simError(); - } + MultipoleAdapter ma = MultipoleAdapter(atomType); + if (ma.isMultipole()) { + if (ma.isDipole()) { electrostaticAtomData.is_Dipole = true; - electrostaticAtomData.dipole_moment = doubleData->getData(); + electrostaticAtomData.dipole_moment = ma.getDipoleMoment(); } - - if (daType->isSplitDipole()) { - GenericData* data = daType->getPropertyByName("SplitDipoleDistance"); - - if (data == NULL) { - sprintf(painCave.errMsg, - "Electrostatic::addType could not find SplitDipoleDistance\n" - "\tparameter for atomType %s.\n", - daType->getName().c_str()); - painCave.severity = OPENMD_ERROR; - painCave.isFatal = 1; - simError(); - } - - DoubleGenericData* doubleData = dynamic_cast(data); - if (doubleData == NULL) { - sprintf( painCave.errMsg, - "Electrostatic::addType could not convert GenericData to " - "SplitDipoleDistance for\n" - "\tatom type %s\n", daType->getName().c_str()); - painCave.severity = OPENMD_ERROR; - painCave.isFatal = 1; - simError(); - } + if (ma.isSplitDipole()) { electrostaticAtomData.is_SplitDipole = true; - electrostaticAtomData.split_dipole_distance = doubleData->getData(); + electrostaticAtomData.split_dipole_distance = ma.getSplitDipoleDistance(); } - - if (daType->isQuadrupole()) { - GenericData* data = daType->getPropertyByName("QuadrupoleMoments"); - - if (data == NULL) { - sprintf( painCave.errMsg, - "Electrostatic::addType could not find QuadrupoleMoments\n" - "\tparameter for atomType %s.\n", - daType->getName().c_str()); - painCave.severity = OPENMD_ERROR; - painCave.isFatal = 1; - simError(); - } - + if (ma.isQuadrupole()) { // 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, - "Electrostatic::addType could not convert GenericData to " - "Quadrupole Moments for\n" - "\tatom type %s\n", daType->getName().c_str()); - painCave.severity = OPENMD_ERROR; - painCave.isFatal = 1; - simError(); - } electrostaticAtomData.is_Quadrupole = true; - electrostaticAtomData.quadrupole_moments = v3dData->getData(); + electrostaticAtomData.quadrupole_moments = ma.getQuadrupoleMoments(); } } - AtomTypeProperties atp = atomType->getATP(); pair::iterator,bool> ret; - ret = ElectrostaticList.insert( pair(atp.ident, atomType) ); + ret = ElectrostaticList.insert( pair(atomType->getIdent(), + atomType) ); if (ret.second == false) { sprintf( painCave.errMsg, "Electrostatic already had a previous entry with ident %d\n", - atp.ident); + atomType->getIdent() ); painCave.severity = OPENMD_INFO; painCave.isFatal = 0; simError(); @@ -407,13 +326,15 @@ namespace OpenMD { return; } - void Electrostatic::setElectrostaticCutoffRadius( RealType theECR, - RealType theRSW ) { - cutoffRadius_ = theECR; + void Electrostatic::setCutoffRadius( RealType rCut ) { + cutoffRadius_ = rCut; rrf_ = cutoffRadius_; - rt_ = theRSW; haveCutoffRadius_ = true; } + + void Electrostatic::setSwitchingRadius( RealType rSwitch ) { + rt_ = rSwitch; + } void Electrostatic::setElectrostaticSummationMethod( ElectrostaticSummationMethod esm ) { summationMethod_ = esm; } @@ -443,13 +364,15 @@ 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; + RealType two(2.0), three(3.0); Vector3d Q_i, Q_j; Vector3d ux_i, uy_i, uz_i; @@ -459,6 +382,12 @@ 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(); @@ -483,8 +412,12 @@ 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; @@ -517,9 +450,14 @@ 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); @@ -559,9 +497,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 { @@ -580,21 +522,30 @@ namespace OpenMD { dudr = *(idat.sw) * preVal * (c2c_ - c2); } else if (summationMethod_ == esm_REACTION_FIELD) { - rfVal = *(idat.electroMult) * preRF_ * *(idat.rij) * *(idat.rij) ; + rfVal = preRF_ * *(idat.rij) * *(idat.rij); + vterm = preVal * ( riji + rfVal ); 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 * two * rfVal * riji * rhat; + } + } else { - vterm = preVal * riji * erfcVal; + 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) { @@ -607,12 +558,23 @@ namespace OpenMD { ri3 = ri2 * riji; vterm = - pref * ct_j * ( ri2 - preRF2_ * *(idat.rij) ); - *(idat.vpair) += vterm; + vpair += vterm; epot += *(idat.sw) * vterm; - dVdr += -preSw * (ri3 * (uz_j - 3.0 * ct_j * rhat) - preRF2_*uz_j); + dVdr += -preSw * (ri3 * (uz_j - three * ct_j * rhat) - preRF2_*uz_j); 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) { @@ -628,9 +590,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; @@ -645,7 +609,7 @@ namespace OpenMD { // calculate the potential pot_term = scale * c2; vterm = -pref * ct_j * pot_term; - *(idat.vpair) += vterm; + vpair += vterm; epot += *(idat.sw) * vterm; // calculate derivatives for forces and torques @@ -665,9 +629,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 * riji; c2 = (-derfcVal + c1) * riji; c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * riji; @@ -684,7 +650,7 @@ namespace OpenMD { c2ri = c2 * riji; c3ri = c3 * riji; c4rij = c4 * *(idat.rij) ; - rhatdot2 = 2.0 * rhat * c3; + rhatdot2 = two * rhat * c3; rhatc4 = rhat * c4rij; // calculate the potential @@ -692,14 +658,14 @@ namespace OpenMD { qyy_j * (cy2*c3 - c2ri) + qzz_j * (cz2*c3 - c2ri) ); vterm = pref * pot_term; - *(idat.vpair) += vterm; + vpair += vterm; epot += *(idat.sw) * vterm; // calculate derivatives for the forces and torques - dVdr += -preSw * ( qxx_j* (cx2*rhatc4 - (2.0*cx_j*ux_j + rhat)*c3ri) + - qyy_j* (cy2*rhatc4 - (2.0*cy_j*uy_j + rhat)*c3ri) + - qzz_j* (cz2*rhatc4 - (2.0*cz_j*uz_j + rhat)*c3ri)); + dVdr += -preSw * ( qxx_j* (cx2*rhatc4 - (two*cx_j*ux_j + rhat)*c3ri) + + qyy_j* (cy2*rhatc4 - (two*cy_j*uy_j + rhat)*c3ri) + + qzz_j* (cz2*rhatc4 - (two*cz_j*uz_j + rhat)*c3ri)); dudux_j += preSw * qxx_j * cx_j * rhatdot2; duduy_j += preSw * qyy_j * cy_j * rhatdot2; @@ -720,12 +686,23 @@ namespace OpenMD { ri3 = ri2 * riji; vterm = pref * ct_i * ( ri2 - preRF2_ * *(idat.rij) ); - *(idat.vpair) += vterm; + vpair += vterm; epot += *(idat.sw) * vterm; - dVdr += preSw * (ri3 * (uz_i - 3.0 * ct_i * rhat) - preRF2_ * uz_i); + dVdr += preSw * (ri3 * (uz_i - three * ct_i * rhat) - preRF2_ * uz_i); 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 { @@ -743,9 +720,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; @@ -760,7 +739,7 @@ namespace OpenMD { // calculate the potential pot_term = c2 * scale; vterm = pref * ct_i * pot_term; - *(idat.vpair) += vterm; + vpair += vterm; epot += *(idat.sw) * vterm; // calculate derivatives for the forces and torques @@ -783,16 +762,23 @@ namespace OpenMD { vterm = pref * ( ri3 * (ct_ij - 3.0 * ct_i * ct_j) - preRF2_ * ct_ij ); - *(idat.vpair) += vterm; + vpair += vterm; epot += *(idat.sw) * vterm; a1 = 5.0 * ct_i * ct_j - ct_ij; - dVdr += preSw * 3.0 * ri4 * (a1 * rhat - ct_i * uz_j - ct_j * uz_i); + dVdr += preSw * three * ri4 * (a1 * rhat - ct_i * uz_j - ct_j * uz_i); - 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); + duduz_i += preSw * (ri3 * (uz_j - three * ct_j * rhat) - preRF2_*uz_j); + duduz_j += preSw * (ri3 * (uz_i - three * 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) { @@ -815,9 +801,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; @@ -842,7 +830,7 @@ namespace OpenMD { // calculate the potential pot_term = (ct_ij * c2ri - ctidotj * c3); vterm = pref * pot_term; - *(idat.vpair) += vterm; + vpair += vterm; epot += *(idat.sw) * vterm; // calculate derivatives for the forces and torques @@ -866,9 +854,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 * riji; c2 = (-derfcVal + c1) * riji; c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * riji; @@ -885,7 +875,7 @@ namespace OpenMD { c2ri = c2 * riji; c3ri = c3 * riji; c4rij = c4 * *(idat.rij) ; - rhatdot2 = 2.0 * rhat * c3; + rhatdot2 = two * rhat * c3; rhatc4 = rhat * c4rij; // calculate the potential @@ -894,14 +884,14 @@ namespace OpenMD { qzz_i * (cz2 * c3 - c2ri) ); vterm = pref * pot_term; - *(idat.vpair) += vterm; + vpair += vterm; epot += *(idat.sw) * vterm; // calculate the derivatives for the forces and torques - dVdr += -preSw * (qxx_i* (cx2*rhatc4 - (2.0*cx_i*ux_i + rhat)*c3ri) + - qyy_i* (cy2*rhatc4 - (2.0*cy_i*uy_i + rhat)*c3ri) + - qzz_i* (cz2*rhatc4 - (2.0*cz_i*uz_i + rhat)*c3ri)); + dVdr += -preSw * (qxx_i* (cx2*rhatc4 - (two*cx_i*ux_i + rhat)*c3ri) + + qyy_i* (cy2*rhatc4 - (two*cy_i*uy_i + rhat)*c3ri) + + qzz_i* (cz2*rhatc4 - (two*cz_i*uz_i + rhat)*c3ri)); dudux_i += preSw * qxx_i * cx_i * rhatdot2; duduy_i += preSw * qyy_i * cy_i * rhatdot2; @@ -909,136 +899,53 @@ namespace OpenMD { } } - 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); - } + 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. - return; - } - - void Electrostatic::calcSkipCorrection(InteractionData &idat) { - - if (!initialized_) initialize(); - - ElectrostaticAtomData data1 = ElectrostaticMap[idat.atypes.first]; - ElectrostaticAtomData data2 = ElectrostaticMap[idat.atypes.second]; - - // 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; - *(idat.skippedCharge2) += q_i; - } - - if (j_is_Charge) { - q_j = data2.charge; - *(idat.skippedCharge1) += q_j; - } - - // the rest of this function should only be necessary for reaction field. - - if (summationMethod_ == esm_REACTION_FIELD) { - RealType riji, ri2, ri3; - RealType mu_i, ct_i; - RealType mu_j, ct_j; - RealType preVal, rfVal, vterm, dudr, pref, myPot(0.0); - Vector3d dVdr, uz_i, uz_j, duduz_i, duduz_j, rhat; - - // some variables we'll need independent of electrostatic type: + *(idat.vpair) += indirect_vpair; + (*(idat.pot))[ELECTROSTATIC_FAMILY] += indirect_Pot; + *(idat.f1) += indirect_dVdr; - riji = 1.0 / *(idat.rij) ; - rhat = *(idat.d) * riji; - - if (i_is_Dipole) { - mu_i = data1.dipole_moment; - uz_i = idat.eFrame1->getColumn(2); - ct_i = dot(uz_i, rhat); - duduz_i = V3Zero; - } - - if (j_is_Dipole) { - mu_j = data2.dipole_moment; - uz_j = idat.eFrame2->getColumn(2); - ct_j = dot(uz_j, rhat); - duduz_j = V3Zero; - } - - if (i_is_Charge) { - if (j_is_Charge) { - preVal = *(idat.electroMult) * pre11_ * q_i * q_j; - rfVal = preRF_ * *(idat.rij) * *(idat.rij) ; - vterm = preVal * rfVal; - myPot += *(idat.sw) * vterm; - dudr = *(idat.sw) * preVal * 2.0 * rfVal * riji; - dVdr += dudr * rhat; - } - - if (j_is_Dipole) { - ri2 = riji * riji; - ri3 = ri2 * riji; - pref = *(idat.electroMult) * pre12_ * q_i * mu_j; - vterm = - pref * ct_j * ( ri2 - preRF2_ * *(idat.rij) ); - myPot += *(idat.sw) * vterm; - dVdr += - *(idat.sw) * pref * ( ri3 * ( uz_j - 3.0 * ct_j * rhat) - preRF2_ * uz_j); - duduz_j += - *(idat.sw) * pref * rhat * (ri2 - preRF2_ * *(idat.rij) ); - } - } - if (i_is_Dipole) { - if (j_is_Charge) { - ri2 = riji * riji; - ri3 = ri2 * riji; - pref = *(idat.electroMult) * pre12_ * q_j * mu_i; - vterm = - pref * ct_i * ( ri2 - preRF2_ * *(idat.rij) ); - myPot += *(idat.sw) * vterm; - dVdr += *(idat.sw) * pref * ( ri3 * ( uz_i - 3.0 * ct_i * rhat) - preRF2_ * uz_i); - duduz_i += *(idat.sw) * pref * rhat * (ri2 - preRF2_ * *(idat.rij)); - } - } - - // accumulate the forces and torques resulting from the self term - idat.pot[ELECTROSTATIC_FAMILY] += myPot; - *(idat.f1) += dVdr; - if (i_is_Dipole) - *(idat.t1) -= cross(uz_i, duduz_i); + *(idat.t1) -= cross(uz_i, indirect_duduz_i); if (j_is_Dipole) - *(idat.t2) -= cross(uz_j, duduz_j); + *(idat.t2) -= cross(uz_j, indirect_duduz_j); } - } + + + return; + } void Electrostatic::calcSelfCorrection(SelfData &sdat) { RealType mu1, preVal, chg1, self; if (!initialized_) initialize(); - + ElectrostaticAtomData data = ElectrostaticMap[sdat.atype]; // logicals - bool i_is_Charge = data.is_Charge; bool i_is_Dipole = data.is_Dipole; @@ -1046,7 +953,7 @@ namespace OpenMD { if (i_is_Dipole) { mu1 = data.dipole_moment; preVal = pre22_ * preRF2_ * mu1 * mu1; - sdat.pot[2] -= 0.5 * preVal; + (*(sdat.pot))[ELECTROSTATIC_FAMILY] -= 0.5 * preVal; // The self-correction term adds into the reaction field vector Vector3d uz_i = sdat.eFrame->getColumn(2); @@ -1063,7 +970,7 @@ namespace OpenMD { } else { self = - 0.5 * rcuti_ * chg1 * (chg1 + *(sdat.skippedCharge)) * pre11_; } - sdat.pot[ELECTROSTATIC_FAMILY] += self; + (*(sdat.pot))[ELECTROSTATIC_FAMILY] += self; } } }