--- branches/development/src/nonbonded/Electrostatic.cpp 2011/08/04 20:04:35 1601 +++ branches/development/src/nonbonded/Electrostatic.cpp 2012/06/05 17:49:36 1734 @@ -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 @@ -46,9 +47,14 @@ #include "nonbonded/Electrostatic.hpp" #include "utils/simError.h" #include "types/NonBondedInteractionType.hpp" -#include "types/DirectionalAtomType.hpp" +#include "types/FixedChargeAdapter.hpp" +#include "types/FluctuatingChargeAdapter.hpp" +#include "types/MultipoleAdapter.hpp" #include "io/Globals.hpp" +#include "nonbonded/SlaterIntegrals.hpp" +#include "utils/PhysicalConstants.hpp" + namespace OpenMD { Electrostatic::Electrostatic(): name_("Electrostatic"), initialized_(false), @@ -64,6 +70,7 @@ namespace OpenMD { 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; @@ -116,7 +123,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; @@ -210,7 +217,6 @@ namespace OpenMD { addType(at); } - cutoffRadius2_ = cutoffRadius_ * cutoffRadius_; rcuti_ = 1.0 / cutoffRadius_; rcuti2_ = rcuti_ * rcuti_; @@ -249,7 +255,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,139 +283,105 @@ namespace OpenMD { electrostaticAtomData.is_Dipole = false; electrostaticAtomData.is_SplitDipole = false; electrostaticAtomData.is_Quadrupole = false; + electrostaticAtomData.is_Fluctuating = 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.fixedCharge = 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(); + FluctuatingChargeAdapter fqa = FluctuatingChargeAdapter(atomType); + if (fqa.isFluctuatingCharge()) { + electrostaticAtomData.is_Fluctuating = true; + electrostaticAtomData.electronegativity = fqa.getElectronegativity(); + electrostaticAtomData.hardness = fqa.getHardness(); + electrostaticAtomData.slaterN = fqa.getSlaterN(); + electrostaticAtomData.slaterZeta = fqa.getSlaterZeta(); + } + 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(); } - ElectrostaticMap[atomType] = electrostaticAtomData; + ElectrostaticMap[atomType] = electrostaticAtomData; + + // Now, iterate over all known types and add to the mixing map: + + map::iterator it; + for( it = ElectrostaticMap.begin(); it != ElectrostaticMap.end(); ++it) { + AtomType* atype2 = (*it).first; + ElectrostaticAtomData eaData2 = (*it).second; + if (eaData2.is_Fluctuating && electrostaticAtomData.is_Fluctuating) { + + RealType a = electrostaticAtomData.slaterZeta; + RealType b = eaData2.slaterZeta; + int m = electrostaticAtomData.slaterN; + int n = eaData2.slaterN; + + // Create the spline of the coulombic integral for s-type + // Slater orbitals. Add a 2 angstrom safety window to deal + // with cutoffGroups that have charged atoms longer than the + // cutoffRadius away from each other. + + RealType rval; + RealType dr = (cutoffRadius_ + 2.0) / RealType(np_ - 1); + vector rvals; + vector J1vals; + vector J2vals; + for (int i = 0; i < np_; i++) { + rval = RealType(i) * dr; + rvals.push_back(rval); + J1vals.push_back(electrostaticAtomData.hardness * sSTOCoulInt( a, b, m, n, rval * PhysicalConstants::angstromsToBohr ) ); + // may not be necessary if Slater coulomb integral is symmetric + J2vals.push_back(eaData2.hardness * sSTOCoulInt( b, a, n, m, rval * PhysicalConstants::angstromsToBohr ) ); + } + + CubicSpline* J1 = new CubicSpline(); + J1->addPoints(rvals, J1vals); + CubicSpline* J2 = new CubicSpline(); + J2->addPoints(rvals, J2vals); + + pair key1, key2; + key1 = make_pair(atomType, atype2); + key2 = make_pair(atype2, atomType); + + Jij[key1] = J1; + Jij[key2] = J2; + } + } + return; } @@ -454,6 +431,7 @@ namespace OpenMD { RealType c1, c2, c3, c4; 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; @@ -469,8 +447,13 @@ namespace OpenMD { Vector3d indirect_dVdr(V3Zero); Vector3d indirect_duduz_i(V3Zero), indirect_duduz_j(V3Zero); + RealType coulInt, vFluc1(0.0), vFluc2(0.0); pair res; + // splines for coulomb integrals + CubicSpline* J1; + CubicSpline* J2; + if (!initialized_) initialize(); ElectrostaticAtomData data1 = ElectrostaticMap[idat.atypes.first]; @@ -487,14 +470,21 @@ namespace OpenMD { bool i_is_Dipole = data1.is_Dipole; bool i_is_SplitDipole = data1.is_SplitDipole; bool i_is_Quadrupole = data1.is_Quadrupole; + bool i_is_Fluctuating = data1.is_Fluctuating; bool j_is_Charge = data2.is_Charge; bool j_is_Dipole = data2.is_Dipole; bool j_is_SplitDipole = data2.is_SplitDipole; bool j_is_Quadrupole = data2.is_Quadrupole; + bool j_is_Fluctuating = data2.is_Fluctuating; if (i_is_Charge) { - q_i = data1.charge; + q_i = data1.fixedCharge; + + if (i_is_Fluctuating) { + q_i += *(idat.flucQ1); + } + if (idat.excluded) { *(idat.skippedCharge2) += q_i; } @@ -532,7 +522,11 @@ namespace OpenMD { } if (j_is_Charge) { - q_j = data2.charge; + q_j = data2.fixedCharge; + + if (j_is_Fluctuating) + q_j += *(idat.flucQ2); + if (idat.excluded) { *(idat.skippedCharge1) += q_j; } @@ -570,6 +564,11 @@ namespace OpenMD { duduz_j = V3Zero; } + if (i_is_Fluctuating && j_is_Fluctuating) { + J1 = Jij[idat.atypes]; + J2 = Jij[make_pair(idat.atypes.second, idat.atypes.first)]; + } + epot = 0.0; dVdr = V3Zero; @@ -578,9 +577,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 { @@ -588,7 +591,7 @@ namespace OpenMD { c2 = c1 * riji; } - preVal = *(idat.electroMult) * pre11_ * q_i * q_j; + preVal = *(idat.electroMult) * pre11_; if (summationMethod_ == esm_SHIFTED_POTENTIAL) { vterm = preVal * (c1 - c1c_); @@ -610,19 +613,45 @@ namespace OpenMD { if (idat.excluded) { indirect_vpair += preVal * rfVal; indirect_Pot += *(idat.sw) * preVal * rfVal; - indirect_dVdr += *(idat.sw) * preVal * 2.0 * rfVal * riji * rhat; + indirect_dVdr += *(idat.sw) * preVal * two * rfVal * riji * rhat; } } else { vterm = preVal * riji * erfcVal; dudr = - *(idat.sw) * preVal * c2; + + } + + vpair += vterm * q_i * q_j; + epot += *(idat.sw) * vterm * q_i * q_j; + dVdr += dudr * rhat * q_i * q_j; + if (i_is_Fluctuating) { + if (idat.excluded) { + // vFluc1 is the difference between the direct coulomb integral + // and the normal 1/r-like interaction between point charges. + coulInt = J1->getValueAt( *(idat.rij) ); + vFluc1 = coulInt - (*(idat.sw) * vterm); + } else { + vFluc1 = 0.0; + } + *(idat.dVdFQ1) += ( *(idat.sw) * vterm + vFluc1 ) * q_j; } - vpair += vterm; - epot += *(idat.sw) * vterm; - dVdr += dudr * rhat; + if (j_is_Fluctuating) { + if (idat.excluded) { + // vFluc2 is the difference between the direct coulomb integral + // and the normal 1/r-like interaction between point charges. + coulInt = J2->getValueAt( *(idat.rij) ); + vFluc2 = coulInt - (*(idat.sw) * vterm); + } else { + vFluc2 = 0.0; + } + *(idat.dVdFQ2) += ( *(idat.sw) * vterm + vFluc2 ) * q_i; + } + + } if (j_is_Dipole) { @@ -638,7 +667,7 @@ namespace OpenMD { 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, @@ -667,9 +696,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; @@ -693,6 +724,9 @@ namespace OpenMD { duduz_j += -preSw * pot_term * rhat; } + if (i_is_Fluctuating) { + *(idat.dVdFQ1) += ( *(idat.sw) * vterm ) / q_i; + } } if (j_is_Quadrupole) { @@ -704,9 +738,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; @@ -723,7 +759,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 @@ -736,13 +772,17 @@ namespace OpenMD { // 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; duduz_j += preSw * qzz_j * cz_j * rhatdot2; + if (i_is_Fluctuating) { + *(idat.dVdFQ1) += ( *(idat.sw) * vterm ) / q_i; + } + } } @@ -762,7 +802,7 @@ namespace OpenMD { 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) ); @@ -793,9 +833,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; @@ -817,6 +859,11 @@ namespace OpenMD { dVdr += preSw * (uz_i * c2ri - ct_i * rhat * sc2 * c3); duduz_i += preSw * pot_term * rhat; } + + if (j_is_Fluctuating) { + *(idat.dVdFQ2) += ( *(idat.sw) * vterm ) / q_j; + } + } if (j_is_Dipole) { @@ -838,10 +885,10 @@ namespace OpenMD { 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; @@ -872,9 +919,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; @@ -923,9 +972,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; @@ -942,7 +993,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 @@ -956,13 +1007,18 @@ namespace OpenMD { // 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; duduz_i += preSw * qzz_i * cz_i * rhatdot2; + + if (j_is_Fluctuating) { + *(idat.dVdFQ2) += ( *(idat.sw) * vterm ) / q_j; + } + } } @@ -990,6 +1046,7 @@ namespace OpenMD { // only accumulate the forces and torques resulting from the // indirect reaction field terms. + *(idat.vpair) += indirect_vpair; (*(idat.pot))[ELECTROSTATIC_FAMILY] += indirect_Pot; *(idat.f1) += indirect_dVdr; @@ -1000,13 +1057,11 @@ namespace OpenMD { *(idat.t2) -= cross(uz_j, indirect_duduz_j); } - return; } void Electrostatic::calcSelfCorrection(SelfData &sdat) { - RealType mu1, preVal, chg1, self; - + RealType mu1, preVal, self; if (!initialized_) initialize(); ElectrostaticAtomData data = ElectrostaticMap[sdat.atype]; @@ -1014,6 +1069,14 @@ namespace OpenMD { // logicals bool i_is_Charge = data.is_Charge; bool i_is_Dipole = data.is_Dipole; + bool i_is_Fluctuating = data.is_Fluctuating; + RealType chg1 = data.fixedCharge; + + if (i_is_Fluctuating) { + chg1 += *(sdat.flucQ); + // dVdFQ is really a force, so this is negative the derivative + *(sdat.dVdFQ) -= *(sdat.flucQ) * data.hardness + data.electronegativity; + } if (summationMethod_ == esm_REACTION_FIELD) { if (i_is_Dipole) { @@ -1030,7 +1093,6 @@ namespace OpenMD { } } 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 + *(sdat.skippedCharge)) * pre11_; } else {