--- branches/development/src/nonbonded/Electrostatic.cpp 2010/10/02 19:53:32 1502 +++ trunk/src/nonbonded/Electrostatic.cpp 2013/08/07 15:24:16 1925 @@ -34,27 +34,60 @@ * 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). - * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008). - * [4] Vardeman & Gezelter, in progress (2009). + * [2] Fennell & Gezelter, J. Chem. Phys. 124 234104 (2006). + * [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). */ #include #include #include +#include #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" +#include "math/erfc.hpp" +#include "math/SquareMatrix.hpp" +#include "primitives/Molecule.hpp" +#ifdef IS_MPI +#include +#endif - namespace OpenMD { Electrostatic::Electrostatic(): name_("Electrostatic"), initialized_(false), - forceField_(NULL) {} + forceField_(NULL), info_(NULL), + haveCutoffRadius_(false), + haveDampingAlpha_(false), + haveDielectric_(false), + haveElectroSplines_(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_["TAYLOR_SHIFTED"] = esm_TAYLOR_SHIFTED; + 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 @@ -62,266 +95,637 @@ namespace OpenMD { // Charge-Dipole, assuming charges are measured in electrons, and // dipoles are measured in debyes pre12_ = 69.13373; - // Dipole-Dipole, assuming dipoles are measured in debyes + // Dipole-Dipole, assuming dipoles are measured in Debye pre22_ = 14.39325; // Charge-Quadrupole, assuming charges are measured in electrons, and // quadrupoles are measured in 10^-26 esu cm^2 - // This unit is also known affectionately as an esu centi-barn. + // This unit is also known affectionately as an esu centibarn. pre14_ = 69.13373; - + // Dipole-Quadrupole, assuming dipoles are measured in debyes and + // quadrupoles in esu centibarns: + pre24_ = 14.39325; + // Quadrupole-Quadrupole, assuming esu centibarns: + pre44_ = 14.39325; + // conversions for the simulation box dipole moment chargeToC_ = 1.60217733e-19; angstromToM_ = 1.0e-10; debyeToCm_ = 3.33564095198e-30; - // number of points for electrostatic splines + // Default number of points for electrostatic splines np_ = 100; // 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; - // 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)) { - - if (at->isElectrostatic()) - addType(at); + // 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\",\n" + "\t\"taylor_shifted\", or \"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 (!haveDefaultCutoff_) { + if (!haveCutoffRadius_) { 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_; - 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; + + if (screeningMethod_ == DAMPED || summationMethod_ == esm_EWALD_FULL) { + 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\n" + "\tinput file. A default value of %f (1/ang) will be used for the\n" + "\tcutoff of %f (ang).\n", + dampingAlpha_, cutoffRadius_); + painCave.severity = OPENMD_INFO; + painCave.isFatal = 0; simError(); + } else { + dampingAlpha_ = simParams_->getDampingAlpha(); } + haveDampingAlpha_ = true; + } - alpha2_ = dampingAlpha_ * dampingAlpha_; - alpha4_ = alpha2_ * alpha2_; - alpha6_ = alpha4_ * alpha2_; - alpha8_ = alpha4_ * alpha4_; - - constEXP_ = exp(-alpha2_ * defaultCutoff2_); - invRootPi_ = 0.56418958354775628695; - alphaPi_ = 2.0 * dampingAlpha_ * invRootPi_; - c1c_ = erfc(dampingAlpha_ * defaultCutoff_) * 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_; - c5c_ = 8.0 * alphaPi_ * alpha6_ + 7.0 * c4c_ * rcuti2_; - c6c_ = 16.0 * alphaPi_ * alpha8_ + 9.0 * c5c_ * rcuti2_; + Etypes.clear(); + Etids.clear(); + FQtypes.clear(); + FQtids.clear(); + ElectrostaticMap.clear(); + Jij.clear(); + nElectro_ = 0; + nFlucq_ = 0; + + Etids.resize( forceField_->getNAtomType(), -1); + FQtids.resize( forceField_->getNAtomType(), -1); + + set::iterator at; + for (at = simTypes_.begin(); at != simTypes_.end(); ++at) { + if ((*at)->isElectrostatic()) nElectro_++; + if ((*at)->isFluctuatingCharge()) nFlucq_++; + } + + Jij.resize(nFlucq_); + + for (at = simTypes_.begin(); at != simTypes_.end(); ++at) { + if ((*at)->isElectrostatic()) addType(*at); + } + + if (summationMethod_ == esm_REACTION_FIELD) { + preRF_ = (dielectric_ - 1.0) / + ((2.0 * dielectric_ + 1.0) * pow(cutoffRadius_,3) ); + } + + RealType b0c, b1c, b2c, b3c, b4c, b5c; + RealType db0c_1, db0c_2, db0c_3, db0c_4, db0c_5; + RealType a2, expTerm, invArootPi; + + RealType r = cutoffRadius_; + RealType r2 = r * r; + RealType ric = 1.0 / r; + RealType ric2 = ric * ric; + + if (screeningMethod_ == DAMPED) { + a2 = dampingAlpha_ * dampingAlpha_; + invArootPi = 1.0 / (dampingAlpha_ * sqrt(M_PI)); + expTerm = exp(-a2 * r2); + // values of Smith's B_l functions at the cutoff radius: + b0c = erfc(dampingAlpha_ * r) / r; + b1c = ( b0c + 2.0*a2 * expTerm * invArootPi) / r2; + b2c = (3.0 * b1c + pow(2.0*a2, 2) * expTerm * invArootPi) / r2; + b3c = (5.0 * b2c + pow(2.0*a2, 3) * expTerm * invArootPi) / r2; + b4c = (7.0 * b3c + pow(2.0*a2, 4) * expTerm * invArootPi) / r2; + b5c = (9.0 * b4c + pow(2.0*a2, 5) * expTerm * invArootPi) / r2; + // Half the Smith self piece: + selfMult1_ = - a2 * invArootPi; + selfMult2_ = - 2.0 * a2 * a2 * invArootPi / 3.0; + selfMult4_ = - 4.0 * a2 * a2 * a2 * invArootPi / 5.0; } else { - c1c_ = rcuti_; - c2c_ = c1c_ * rcuti_; - c3c_ = 3.0 * c2c_ * rcuti_; - c4c_ = 5.0 * c3c_ * rcuti2_; - c5c_ = 7.0 * c4c_ * rcuti2_; - c6c_ = 9.0 * c5c_ * rcuti2_; + a2 = 0.0; + b0c = 1.0 / r; + b1c = ( b0c) / r2; + b2c = (3.0 * b1c) / r2; + b3c = (5.0 * b2c) / r2; + b4c = (7.0 * b3c) / r2; + b5c = (9.0 * b4c) / r2; + selfMult1_ = 0.0; + selfMult2_ = 0.0; + selfMult4_ = 0.0; } - - if (summationMethod_ == REACTION_FIELD) { - if (haveDielectric_) { - preRF_ = (dielectric_ - 1.0) / - ((2.0 * dielectric_ + 1.0) * defaultCutoff2_ * defaultCutoff_); - preRF2_ = 2.0 * preRF_; + + // higher derivatives of B_0 at the cutoff radius: + db0c_1 = -r * b1c; + db0c_2 = -b1c + r2 * b2c; + db0c_3 = 3.0*r*b2c - r2*r*b3c; + db0c_4 = 3.0*b2c - 6.0*r2*b3c + r2*r2*b4c; + db0c_5 = -15.0*r*b3c + 10.0*r2*r*b4c - r2*r2*r*b5c; + + if (summationMethod_ != esm_EWALD_FULL) { + selfMult1_ -= b0c; + selfMult2_ += (db0c_2 + 2.0*db0c_1*ric) / 3.0; + selfMult4_ -= (db0c_4 + 4.0*db0c_3*ric) / 15.0; + } + + // working variables for the splines: + RealType ri, ri2; + RealType b0, b1, b2, b3, b4, b5; + RealType db0_1, db0_2, db0_3, db0_4, db0_5; + RealType f, fc, f0; + RealType g, gc, g0, g1, g2, g3, g4; + RealType h, hc, h1, h2, h3, h4; + RealType s, sc, s2, s3, s4; + RealType t, tc, t3, t4; + RealType u, uc, u4; + + // working variables for Taylor expansion: + RealType rmRc, rmRc2, rmRc3, rmRc4; + + // Approximate using splines using a maximum of 0.1 Angstroms + // between points. + int nptest = int((cutoffRadius_ + 2.0) / 0.1); + np_ = (np_ > nptest) ? np_ : nptest; + + // Add a 2 angstrom safety window to deal with cutoffGroups that + // have charged atoms longer than the cutoffRadius away from each + // other. Splining is almost certainly the best choice here. + // Direct calls to erfc would be preferrable if it is a very fast + // implementation. + + RealType dx = (cutoffRadius_ + 2.0) / RealType(np_); + + // Storage vectors for the computed functions + vector rv; + vector v01v; + vector v11v; + vector v21v, v22v; + vector v31v, v32v; + vector v41v, v42v, v43v; + + for (int i = 1; i < np_ + 1; i++) { + r = RealType(i) * dx; + rv.push_back(r); + + ri = 1.0 / r; + ri2 = ri * ri; + + r2 = r * r; + expTerm = exp(-a2 * r2); + + // Taylor expansion factors (no need for factorials this way): + rmRc = r - cutoffRadius_; + rmRc2 = rmRc * rmRc / 2.0; + rmRc3 = rmRc2 * rmRc / 3.0; + rmRc4 = rmRc3 * rmRc / 4.0; + + // values of Smith's B_l functions at r: + if (screeningMethod_ == DAMPED) { + b0 = erfc(dampingAlpha_ * r) * ri; + b1 = ( b0 + 2.0*a2 * expTerm * invArootPi) * ri2; + b2 = (3.0 * b1 + pow(2.0*a2, 2) * expTerm * invArootPi) * ri2; + b3 = (5.0 * b2 + pow(2.0*a2, 3) * expTerm * invArootPi) * ri2; + b4 = (7.0 * b3 + pow(2.0*a2, 4) * expTerm * invArootPi) * ri2; + b5 = (9.0 * b4 + pow(2.0*a2, 5) * expTerm * invArootPi) * ri2; } else { - sprintf( painCave.errMsg, "Electrostatic::initialize has no Dielectric" - " value!\n"); - painCave.severity = OPENMD_ERROR; + b0 = ri; + b1 = ( b0) * ri2; + b2 = (3.0 * b1) * ri2; + b3 = (5.0 * b2) * ri2; + b4 = (7.0 * b3) * ri2; + b5 = (9.0 * b4) * ri2; + } + + // higher derivatives of B_0 at r: + db0_1 = -r * b1; + db0_2 = -b1 + r2 * b2; + db0_3 = 3.0*r*b2 - r2*r*b3; + db0_4 = 3.0*b2 - 6.0*r2*b3 + r2*r2*b4; + db0_5 = -15.0*r*b3 + 10.0*r2*r*b4 - r2*r2*r*b5; + + f = b0; + fc = b0c; + f0 = f - fc - rmRc*db0c_1; + + g = db0_1; + gc = db0c_1; + g0 = g - gc; + g1 = g0 - rmRc *db0c_2; + g2 = g1 - rmRc2*db0c_3; + g3 = g2 - rmRc3*db0c_4; + g4 = g3 - rmRc4*db0c_5; + + h = db0_2; + hc = db0c_2; + h1 = h - hc; + h2 = h1 - rmRc *db0c_3; + h3 = h2 - rmRc2*db0c_4; + h4 = h3 - rmRc3*db0c_5; + + s = db0_3; + sc = db0c_3; + s2 = s - sc; + s3 = s2 - rmRc *db0c_4; + s4 = s3 - rmRc2*db0c_5; + + t = db0_4; + tc = db0c_4; + t3 = t - tc; + t4 = t3 - rmRc *db0c_5; + + u = db0_5; + uc = db0c_5; + u4 = u - uc; + + // in what follows below, the various v functions are used for + // potentials and torques, while the w functions show up in the + // forces. + + switch (summationMethod_) { + case esm_SHIFTED_FORCE: + + v01 = f - fc - rmRc*gc; + v11 = g - gc - rmRc*hc; + v21 = g*ri - gc*ric - rmRc*(hc - gc*ric)*ric; + v22 = h - g*ri - (hc - gc*ric) - rmRc*(sc - (hc - gc*ric)*ric); + v31 = (h-g*ri)*ri - (hc-gc*ric)*ric - rmRc*(sc-2.0*(hc-gc*ric)*ric)*ric; + v32 = (s - 3.0*(h-g*ri)*ri) - (sc - 3.0*(hc-gc*ric)*ric) + - rmRc*(tc - 3.0*(sc-2.0*(hc-gc*ric)*ric)*ric); + v41 = (h - g*ri)*ri2 - (hc - gc*ric)*ric2 + - rmRc*(sc - 3.0*(hc-gc*ric)*ric)*ric2; + v42 = (s-3.0*(h-g*ri)*ri)*ri - (sc-3.0*(hc-gc*ric)*ric)*ric + - rmRc*(tc - (4.0*sc - 9.0*(hc - gc*ric)*ric)*ric)*ric; + + v43 = (t - 3.0*(2.0*s - 5.0*(h - g*ri)*ri)*ri) + - (tc - 3.0*(2.0*sc - 5.0*(hc - gc*ric)*ric)*ric) + - rmRc*(uc-3.0*(2.0*tc - (7.0*sc - 15.0*(hc - gc*ric)*ric)*ric)*ric); + + dv01 = g - gc; + dv11 = h - hc; + dv21 = (h - g*ri)*ri - (hc - gc*ric)*ric; + dv22 = (s - (h - g*ri)*ri) - (sc - (hc - gc*ric)*ric); + dv31 = (s - 2.0*(h-g*ri)*ri)*ri - (sc - 2.0*(hc-gc*ric)*ric)*ric; + dv32 = (t - 3.0*(s-2.0*(h-g*ri)*ri)*ri) + - (tc - 3.0*(sc-2.0*(hc-gc*ric)*ric)*ric); + dv41 = (s - 3.0*(h - g*ri)*ri)*ri2 - (sc - 3.0*(hc - gc*ric)*ric)*ric2; + dv42 = (t - (4.0*s - 9.0*(h-g*ri)*ri)*ri)*ri + - (tc - (4.0*sc - 9.0*(hc-gc*ric)*ric)*ric)*ric; + dv43 = (u - 3.0*(2.0*t - (7.0*s - 15.0*(h - g*ri)*ri)*ri)*ri) + - (uc - 3.0*(2.0*tc - (7.0*sc - 15.0*(hc - gc*ric)*ric)*ric)*ric); + + break; + + case esm_TAYLOR_SHIFTED: + + v01 = f0; + v11 = g1; + v21 = g2 * ri; + v22 = h2 - v21; + v31 = (h3 - g3 * ri) * ri; + v32 = s3 - 3.0*v31; + v41 = (h4 - g4 * ri) * ri2; + v42 = s4 * ri - 3.0*v41; + v43 = t4 - 6.0*v42 - 3.0*v41; + + dv01 = g0; + dv11 = h1; + dv21 = (h2 - g2*ri)*ri; + dv22 = (s2 - (h2 - g2*ri)*ri); + dv31 = (s3 - 2.0*(h3-g3*ri)*ri)*ri; + dv32 = (t3 - 3.0*(s3-2.0*(h3-g3*ri)*ri)*ri); + dv41 = (s4 - 3.0*(h4 - g4*ri)*ri)*ri2; + dv42 = (t4 - (4.0*s4 - 9.0*(h4-g4*ri)*ri)*ri)*ri; + dv43 = (u4 - 3.0*(2.0*t4 - (7.0*s4 - 15.0*(h4 - g4*ri)*ri)*ri)*ri); + + break; + + case esm_SHIFTED_POTENTIAL: + + v01 = f - fc; + v11 = g - gc; + v21 = g*ri - gc*ric; + v22 = h - g*ri - (hc - gc*ric); + v31 = (h-g*ri)*ri - (hc-gc*ric)*ric; + v32 = (s - 3.0*(h-g*ri)*ri) - (sc - 3.0*(hc-gc*ric)*ric); + v41 = (h - g*ri)*ri2 - (hc - gc*ric)*ric2; + v42 = (s-3.0*(h-g*ri)*ri)*ri - (sc-3.0*(hc-gc*ric)*ric)*ric; + v43 = (t - 3.0*(2.0*s - 5.0*(h - g*ri)*ri)*ri) + - (tc - 3.0*(2.0*sc - 5.0*(hc - gc*ric)*ric)*ric); + + dv01 = g; + dv11 = h; + dv21 = (h - g*ri)*ri; + dv22 = (s - (h - g*ri)*ri); + dv31 = (s - 2.0*(h-g*ri)*ri)*ri; + dv32 = (t - 3.0*(s-2.0*(h-g*ri)*ri)*ri); + dv41 = (s - 3.0*(h - g*ri)*ri)*ri2; + dv42 = (t - (4.0*s - 9.0*(h-g*ri)*ri)*ri)*ri; + dv43 = (u - 3.0*(2.0*t - (7.0*s - 15.0*(h - g*ri)*ri)*ri)*ri); + + break; + + case esm_SWITCHING_FUNCTION: + case esm_HARD: + case esm_EWALD_FULL: + + v01 = f; + v11 = g; + v21 = g*ri; + v22 = h - g*ri; + v31 = (h-g*ri)*ri; + v32 = (s - 3.0*(h-g*ri)*ri); + v41 = (h - g*ri)*ri2; + v42 = (s-3.0*(h-g*ri)*ri)*ri; + v43 = (t - 3.0*(2.0*s - 5.0*(h - g*ri)*ri)*ri); + + dv01 = g; + dv11 = h; + dv21 = (h - g*ri)*ri; + dv22 = (s - (h - g*ri)*ri); + dv31 = (s - 2.0*(h-g*ri)*ri)*ri; + dv32 = (t - 3.0*(s-2.0*(h-g*ri)*ri)*ri); + dv41 = (s - 3.0*(h - g*ri)*ri)*ri2; + dv42 = (t - (4.0*s - 9.0*(h-g*ri)*ri)*ri)*ri; + dv43 = (u - 3.0*(2.0*t - (7.0*s - 15.0*(h - g*ri)*ri)*ri)*ri); + + break; + + case esm_REACTION_FIELD: + + // following DL_POLY's lead for shifting the image charge potential: + f = b0 + preRF_ * r2; + fc = b0c + preRF_ * cutoffRadius_ * cutoffRadius_; + + g = db0_1 + preRF_ * 2.0 * r; + gc = db0c_1 + preRF_ * 2.0 * cutoffRadius_; + + h = db0_2 + preRF_ * 2.0; + hc = db0c_2 + preRF_ * 2.0; + + v01 = f - fc; + v11 = g - gc; + v21 = g*ri - gc*ric; + v22 = h - g*ri - (hc - gc*ric); + v31 = (h-g*ri)*ri - (hc-gc*ric)*ric; + v32 = (s - 3.0*(h-g*ri)*ri) - (sc - 3.0*(hc-gc*ric)*ric); + v41 = (h - g*ri)*ri2 - (hc - gc*ric)*ric2; + v42 = (s-3.0*(h-g*ri)*ri)*ri - (sc-3.0*(hc-gc*ric)*ric)*ric; + v43 = (t - 3.0*(2.0*s - 5.0*(h - g*ri)*ri)*ri) + - (tc - 3.0*(2.0*sc - 5.0*(hc - gc*ric)*ric)*ric); + + dv01 = g; + dv11 = h; + dv21 = (h - g*ri)*ri; + dv22 = (s - (h - g*ri)*ri); + dv31 = (s - 2.0*(h-g*ri)*ri)*ri; + dv32 = (t - 3.0*(s-2.0*(h-g*ri)*ri)*ri); + dv41 = (s - 3.0*(h - g*ri)*ri)*ri2; + dv42 = (t - (4.0*s - 9.0*(h-g*ri)*ri)*ri)*ri; + dv43 = (u - 3.0*(2.0*t - (7.0*s - 15.0*(h - g*ri)*ri)*ri)*ri); + + break; + + case esm_EWALD_PME: + case esm_EWALD_SPME: + default : + map::iterator i; + std::string meth; + for (i = summationMap_.begin(); i != summationMap_.end(); ++i) { + if ((*i).second == summationMethod_) meth = (*i).first; + } + sprintf( painCave.errMsg, + "Electrostatic::initialize: electrostaticSummationMethod %s \n" + "\thas not been implemented yet. Please select one of:\n" + "\t\"hard\", \"shifted_potential\", or \"shifted_force\"\n", + meth.c_str() ); painCave.isFatal = 1; simError(); + break; } + + // Add these computed values to the storage vectors for spline creation: + v01v.push_back(v01); + v11v.push_back(v11); + v21v.push_back(v21); + v22v.push_back(v22); + v31v.push_back(v31); + v32v.push_back(v32); + v41v.push_back(v41); + v42v.push_back(v42); + v43v.push_back(v43); } - - RealType dx = defaultCutoff_ / RealType(np_ - 1); - RealType rval; - vector rvals; - vector yvals; - for (int i = 0; i < np_; i++) { - rval = RealType(i) * dx; - rvals.push_back(rval); - yvals.push_back(erfc(dampingAlpha_ * rval)); - } - erfcSpline_ = new CubicSpline(); - erfcSpline_->addPoints(rvals, yvals); - haveElectroSpline_ = true; + // construct the spline structures and fill them with the values we've + // computed: + + v01s = new CubicSpline(); + v01s->addPoints(rv, v01v); + v11s = new CubicSpline(); + v11s->addPoints(rv, v11v); + v21s = new CubicSpline(); + v21s->addPoints(rv, v21v); + v22s = new CubicSpline(); + v22s->addPoints(rv, v22v); + v31s = new CubicSpline(); + v31s->addPoints(rv, v31v); + v32s = new CubicSpline(); + v32s->addPoints(rv, v32v); + v41s = new CubicSpline(); + v41s->addPoints(rv, v41v); + v42s = new CubicSpline(); + v42s->addPoints(rv, v42v); + v43s = new CubicSpline(); + v43s->addPoints(rv, v43v); + + haveElectroSplines_ = true; + initialized_ = true; } void Electrostatic::addType(AtomType* atomType){ - + ElectrostaticAtomData electrostaticAtomData; electrostaticAtomData.is_Charge = false; 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 = ma.getDipole(); } - - 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(); - } - electrostaticAtomData.is_SplitDipole = true; - electrostaticAtomData.split_dipole_distance = doubleData->getData(); - } - - 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(); - } - - 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(); - } + if (ma.isQuadrupole()) { electrostaticAtomData.is_Quadrupole = true; - electrostaticAtomData.quadrupole_moments = v3dData->getData(); + electrostaticAtomData.quadrupole = ma.getQuadrupole(); } } - AtomTypeProperties atp = atomType->getATP(); + FluctuatingChargeAdapter fqa = FluctuatingChargeAdapter(atomType); - pair::iterator,bool> ret; - ret = ElectrostaticList.insert( pair(atp.ident, atomType) ); + if (fqa.isFluctuatingCharge()) { + electrostaticAtomData.is_Fluctuating = true; + electrostaticAtomData.electronegativity = fqa.getElectronegativity(); + electrostaticAtomData.hardness = fqa.getHardness(); + electrostaticAtomData.slaterN = fqa.getSlaterN(); + electrostaticAtomData.slaterZeta = fqa.getSlaterZeta(); + } + + int atid = atomType->getIdent(); + int etid = Etypes.size(); + int fqtid = FQtypes.size(); + + pair::iterator,bool> ret; + ret = Etypes.insert( atid ); if (ret.second == false) { sprintf( painCave.errMsg, "Electrostatic already had a previous entry with ident %d\n", - atp.ident); + atid); painCave.severity = OPENMD_INFO; painCave.isFatal = 0; simError(); } - ElectrostaticMap[atomType] = electrostaticAtomData; + Etids[ atid ] = etid; + ElectrostaticMap.push_back(electrostaticAtomData); + + if (electrostaticAtomData.is_Fluctuating) { + ret = FQtypes.insert( atid ); + if (ret.second == false) { + sprintf( painCave.errMsg, + "Electrostatic already had a previous fluctuating charge entry with ident %d\n", + atid ); + painCave.severity = OPENMD_INFO; + painCave.isFatal = 0; + simError(); + } + FQtids[atid] = fqtid; + Jij[fqtid].resize(nFlucq_); + + // Now, iterate over all known fluctuating and add to the + // coulomb integral map: + + std::set::iterator it; + for( it = FQtypes.begin(); it != FQtypes.end(); ++it) { + int etid2 = Etids[ (*it) ]; + int fqtid2 = FQtids[ (*it) ]; + ElectrostaticAtomData eaData2 = ElectrostaticMap[ etid2 ]; + 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 Jvals; + // don't start at i = 0, as rval = 0 is undefined for the + // slater overlap integrals. + for (int i = 1; i < np_+1; i++) { + rval = RealType(i) * dr; + rvals.push_back(rval); + Jvals.push_back(sSTOCoulInt( a, b, m, n, rval * + PhysicalConstants::angstromToBohr ) * + PhysicalConstants::hartreeToKcal ); + } + + CubicSpline* J = new CubicSpline(); + J->addPoints(rvals, Jvals); + Jij[fqtid][fqtid2] = J; + Jij[fqtid2].resize( nFlucq_ ); + Jij[fqtid2][fqtid] = J; + } + } return; } - void Electrostatic::setElectrostaticCutoffRadius( RealType theECR, - RealType theRSW ) { - defaultCutoff_ = theECR; - rrf_ = defaultCutoff_; - rt_ = theRSW; - haveDefaultCutoff_ = true; + void Electrostatic::setCutoffRadius( RealType rCut ) { + cutoffRadius_ = rCut; + haveCutoffRadius_ = true; } + void Electrostatic::setElectrostaticSummationMethod( ElectrostaticSummationMethod esm ) { summationMethod_ = esm; } @@ -337,642 +741,787 @@ 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: - - RealType q_i, q_j, mu_i, mu_j, d_i, d_j; - RealType qxx_i, qyy_i, qzz_i; - RealType qxx_j, qyy_j, qzz_j; - RealType cx_i, cy_i, cz_i; - RealType cx_j, cy_j, cz_j; - RealType cx2, cy2, cz2; - RealType ct_i, ct_j, ct_ij, a1; - RealType riji, ri, ri2, ri3, ri4; - RealType pref, vterm, epot, dudr; - 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 BigR; - - Vector3d Q_i, Q_j; - Vector3d ux_i, uy_i, uz_i; - Vector3d ux_j, uy_j, uz_j; - Vector3d dudux_i, duduy_i, duduz_i; - Vector3d dudux_j, duduy_j, duduz_j; - Vector3d rhatdot2, rhatc4; - Vector3d dVdr; - - pair res; - if (!initialized_) initialize(); - ElectrostaticAtomData data1 = ElectrostaticMap[idat.atype1]; - ElectrostaticAtomData data2 = ElectrostaticMap[idat.atype2]; + data1 = ElectrostaticMap[Etids[idat.atid1]]; + data2 = ElectrostaticMap[Etids[idat.atid2]]; + + U = 0.0; // Potential + F.zero(); // Force + Ta.zero(); // Torque on site a + Tb.zero(); // Torque on site b + Ea.zero(); // Electric field at site a + Eb.zero(); // Electric field at site b + dUdCa = 0.0; // fluctuating charge force at site a + dUdCb = 0.0; // fluctuating charge force at site a - // some variables we'll need independent of electrostatic type: + // Indirect interactions mediated by the reaction field. + indirect_Pot = 0.0; // Potential + indirect_F.zero(); // Force + indirect_Ta.zero(); // Torque on site a + indirect_Tb.zero(); // Torque on site b - riji = 1.0 / idat.rij; - Vector3d rhat = idat.d * riji; + // Excluded potential that is still computed for fluctuating charges + excluded_Pot= 0.0; - // logicals - bool i_is_Charge = data1.is_Charge; - bool i_is_Dipole = data1.is_Dipole; - bool i_is_SplitDipole = data1.is_SplitDipole; - bool i_is_Quadrupole = data1.is_Quadrupole; + // some variables we'll need independent of electrostatic type: - 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; - - if (i_is_Charge) - q_i = data1.charge; - - if (i_is_Dipole) { - mu_i = data1.dipole_moment; - uz_i = idat.eFrame1.getColumn(2); + ri = 1.0 / *(idat.rij); + rhat = *(idat.d) * ri; - ct_i = dot(uz_i, rhat); + // logicals - if (i_is_SplitDipole) - d_i = data1.split_dipole_distance; - - duduz_i = V3Zero; - } - - if (i_is_Quadrupole) { - Q_i = data1.quadrupole_moments; - qxx_i = Q_i.x(); - 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); + a_is_Charge = data1.is_Charge; + a_is_Dipole = data1.is_Dipole; + a_is_Quadrupole = data1.is_Quadrupole; + a_is_Fluctuating = data1.is_Fluctuating; - cx_i = dot(ux_i, rhat); - cy_i = dot(uy_i, rhat); - cz_i = dot(uz_i, rhat); + b_is_Charge = data2.is_Charge; + b_is_Dipole = data2.is_Dipole; + b_is_Quadrupole = data2.is_Quadrupole; + b_is_Fluctuating = data2.is_Fluctuating; - dudux_i = V3Zero; - duduy_i = V3Zero; - duduz_i = V3Zero; + // Obtain all of the required radial function values from the + // spline structures: + + // needed for fields (and forces): + if (a_is_Charge || b_is_Charge) { + v01s->getValueAndDerivativeAt( *(idat.rij), v01, dv01); } + if (a_is_Dipole || b_is_Dipole) { + v11s->getValueAndDerivativeAt( *(idat.rij), v11, dv11); + v11or = ri * v11; + } + if (a_is_Quadrupole || b_is_Quadrupole || (a_is_Dipole && b_is_Dipole)) { + v21s->getValueAndDerivativeAt( *(idat.rij), v21, dv21); + v22s->getValueAndDerivativeAt( *(idat.rij), v22, dv22); + v22or = ri * v22; + } - if (j_is_Charge) - q_j = data2.charge; + // needed for potentials (and forces and torques): + if ((a_is_Dipole && b_is_Quadrupole) || + (b_is_Dipole && a_is_Quadrupole)) { + v31s->getValueAndDerivativeAt( *(idat.rij), v31, dv31); + v32s->getValueAndDerivativeAt( *(idat.rij), v32, dv32); + v31or = v31 * ri; + v32or = v32 * ri; + } + if (a_is_Quadrupole && b_is_Quadrupole) { + v41s->getValueAndDerivativeAt( *(idat.rij), v41, dv41); + v42s->getValueAndDerivativeAt( *(idat.rij), v42, dv42); + v43s->getValueAndDerivativeAt( *(idat.rij), v43, dv43); + v42or = v42 * ri; + v43or = v43 * ri; + } - if (j_is_Dipole) { - mu_j = data2.dipole_moment; - uz_j = idat.eFrame2.getColumn(2); + // calculate the single-site contributions (fields, etc). + + if (a_is_Charge) { + C_a = data1.fixedCharge; - ct_j = dot(uz_j, rhat); - - if (j_is_SplitDipole) - d_j = data2.split_dipole_distance; + if (a_is_Fluctuating) { + C_a += *(idat.flucQ1); + } - duduz_j = V3Zero; + if (idat.excluded) { + *(idat.skippedCharge2) += C_a; + } else { + // only do the field if we're not excluded: + Eb -= C_a * pre11_ * dv01 * rhat; + } } - if (j_is_Quadrupole) { - Q_j = data2.quadrupole_moments; - qxx_j = Q_j.x(); - 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); - - cx_j = dot(ux_j, rhat); - cy_j = dot(uy_j, rhat); - cz_j = dot(uz_j, rhat); - - dudux_j = V3Zero; - duduy_j = V3Zero; - duduz_j = V3Zero; + if (a_is_Dipole) { + D_a = *(idat.dipole1); + rdDa = dot(rhat, D_a); + rxDa = cross(rhat, D_a); + if (!idat.excluded) + Eb -= pre12_ * ((dv11-v11or) * rdDa * rhat + v11or * D_a); } - epot = 0.0; - dVdr = V3Zero; + if (a_is_Quadrupole) { + Q_a = *(idat.quadrupole1); + trQa = Q_a.trace(); + Qar = Q_a * rhat; + rQa = rhat * Q_a; + rdQar = dot(rhat, Qar); + rxQar = cross(rhat, Qar); + if (!idat.excluded) + Eb -= pre14_ * (trQa * rhat * dv21 + 2.0 * Qar * v22or + + rdQar * rhat * (dv22 - 2.0*v22or)); + } - if (i_is_Charge) { + if (b_is_Charge) { + C_b = data2.fixedCharge; - if (j_is_Charge) { - if (screeningMethod_ == DAMPED) { - // assemble the damping variables - res = erfcSpline_->getValueAndDerivativeAt(idat.rij); - erfcVal = res.first; - derfcVal = res.second; - c1 = erfcVal * riji; - c2 = (-derfcVal + c1) * riji; - } else { - c1 = riji; - c2 = c1 * riji; - } - - preVal = idat.electroMult * pre11_ * q_i * q_j; - - if (summationMethod_ == SHIFTED_POTENTIAL) { - vterm = preVal * (c1 - c1c_); - 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_ == REACTION_FIELD) { - rfVal = idat.electroMult * preRF_ * idat.rij * idat.rij; - vterm = preVal * ( riji + rfVal ); - dudr = idat.sw * preVal * ( 2.0 * rfVal - riji ) * riji; - - } else { - vterm = preVal * riji * erfcVal; - - dudr = - idat.sw * preVal * c2; - - } - - idat.vpair += vterm; - epot += idat.sw * vterm; - - dVdr += dudr * rhat; + if (b_is_Fluctuating) + C_b += *(idat.flucQ2); + + if (idat.excluded) { + *(idat.skippedCharge1) += C_b; + } else { + // only do the field if we're not excluded: + Ea += C_b * pre11_ * dv01 * 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; - - if (summationMethod_ == REACTION_FIELD) { - ri2 = riji * riji; - ri3 = ri2 * riji; + } - vterm = - pref * ct_j * ( ri2 - preRF2_ * idat.rij ); - idat.vpair += vterm; - epot += idat.sw * vterm; + if (b_is_Dipole) { + D_b = *(idat.dipole2); + rdDb = dot(rhat, D_b); + rxDb = cross(rhat, D_b); + if (!idat.excluded) + Ea += pre12_ * ((dv11-v11or) * rdDb * rhat + v11or * D_b); + } + + if (b_is_Quadrupole) { + Q_b = *(idat.quadrupole2); + trQb = Q_b.trace(); + Qbr = Q_b * rhat; + rQb = rhat * Q_b; + rdQbr = dot(rhat, Qbr); + rxQbr = cross(rhat, Qbr); + if (!idat.excluded) + Ea += pre14_ * (trQb * rhat * dv21 + 2.0 * Qbr * v22or + + rdQbr * rhat * (dv22 - 2.0*v22or)); + } + + if ((a_is_Fluctuating || b_is_Fluctuating) && idat.excluded) { + J = Jij[FQtids[idat.atid1]][FQtids[idat.atid2]]; + } + + if (a_is_Charge) { + + if (b_is_Charge) { + pref = pre11_ * *(idat.electroMult); + U += C_a * C_b * pref * v01; + F += C_a * C_b * pref * dv01 * rhat; + + // If this is an excluded pair, there are still indirect + // interactions via the reaction field we must worry about: - dVdr += -preSw * (ri3 * (uz_j - 3.0 * ct_j * rhat) - preRF2_*uz_j); - duduz_j += -preSw * rhat * (ri2 - preRF2_ * idat.rij); + if (summationMethod_ == esm_REACTION_FIELD && idat.excluded) { + rfContrib = preRF_ * pref * C_a * C_b * *(idat.r2); + indirect_Pot += rfContrib; + indirect_F += rfContrib * 2.0 * ri * rhat; + } + + // Fluctuating charge forces are handled via Coulomb integrals + // for excluded pairs (i.e. those connected via bonds) and + // with the standard charge-charge interaction otherwise. + if (idat.excluded) { + if (a_is_Fluctuating || b_is_Fluctuating) { + coulInt = J->getValueAt( *(idat.rij) ); + if (a_is_Fluctuating) dUdCa += coulInt * C_b; + if (b_is_Fluctuating) dUdCb += coulInt * C_a; + excluded_Pot += C_a * C_b * coulInt; + } } 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); - ri = 1.0 / BigR; - scale = idat.rij * ri; - } else { - ri = riji; - scale = 1.0; - } - - sc2 = scale * scale; + if (a_is_Fluctuating) dUdCa += C_b * pref * v01; + if (a_is_Fluctuating) dUdCb += C_a * pref * v01; + } + } - if (screeningMethod_ == DAMPED) { - // assemble the damping variables - res = erfcSpline_->getValueAndDerivativeAt(idat.rij); - erfcVal = res.first; - derfcVal = res.second; - c1 = erfcVal * ri; - c2 = (-derfcVal + c1) * ri; - c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * ri; - } else { - c1 = ri; - c2 = c1 * ri; - c3 = 3.0 * c2 * ri; - } - - c2ri = c2 * ri; + if (b_is_Dipole) { + pref = pre12_ * *(idat.electroMult); + U += C_a * pref * v11 * rdDb; + F += C_a * pref * ((dv11 - v11or) * rdDb * rhat + v11or * D_b); + Tb += C_a * pref * v11 * rxDb; - // calculate the potential - pot_term = scale * c2; - vterm = -pref * ct_j * pot_term; - idat.vpair += vterm; - epot += idat.sw * vterm; - - // calculate derivatives for forces and torques + if (a_is_Fluctuating) dUdCa += pref * v11 * rdDb; - dVdr += -preSw * (uz_j * c2ri - ct_j * rhat * sc2 * c3); - duduz_j += -preSw * pot_term * rhat; + // Even if we excluded this pair from direct interactions, we + // still have the reaction-field-mediated charge-dipole + // interaction: + if (summationMethod_ == esm_REACTION_FIELD && idat.excluded) { + rfContrib = C_a * pref * preRF_ * 2.0 * *(idat.rij); + indirect_Pot += rfContrib * rdDb; + indirect_F += rfContrib * D_b / (*idat.rij); + indirect_Tb += C_a * pref * preRF_ * rxDb; } } - if (j_is_Quadrupole) { - // first precalculate some necessary variables - cx2 = cx_j * cx_j; - cy2 = cy_j * cy_j; - cz2 = cz_j * cz_j; - 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; - c1 = erfcVal * riji; - c2 = (-derfcVal + c1) * riji; - c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * riji; - c4 = -4.0 * derfcVal * alpha4_ + 5.0 * c3 * riji * riji; - } else { - c1 = riji; - c2 = c1 * riji; - c3 = 3.0 * c2 * riji; - c4 = 5.0 * c3 * riji * riji; - } + if (b_is_Quadrupole) { + pref = pre14_ * *(idat.electroMult); + U += C_a * pref * (v21 * trQb + v22 * rdQbr); + F += C_a * pref * (trQb * dv21 * rhat + 2.0 * Qbr * v22or); + F += C_a * pref * rdQbr * rhat * (dv22 - 2.0*v22or); + Tb += C_a * pref * 2.0 * rxQbr * v22; - // precompute variables for convenience - preSw = idat.sw * pref; - c2ri = c2 * riji; - c3ri = c3 * riji; - c4rij = c4 * idat.rij; - rhatdot2 = 2.0 * rhat * c3; - rhatc4 = rhat * c4rij; - - // calculate the potential - pot_term = ( qxx_j * (cx2*c3 - c2ri) + - qyy_j * (cy2*c3 - c2ri) + - qzz_j * (cz2*c3 - c2ri) ); - vterm = pref * pot_term; - idat.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)); - - 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 (a_is_Fluctuating) dUdCa += pref * (v21 * trQb + v22 * rdQbr); } } - - if (i_is_Dipole) { - if (j_is_Charge) { - // variables used by all the methods - pref = idat.electroMult * pre12_ * q_j * mu_i; - preSw = idat.sw * pref; + if (a_is_Dipole) { - if (summationMethod_ == REACTION_FIELD) { + if (b_is_Charge) { + pref = pre12_ * *(idat.electroMult); - ri2 = riji * riji; - ri3 = ri2 * riji; + U -= C_b * pref * v11 * rdDa; + F -= C_b * pref * ((dv11-v11or) * rdDa * rhat + v11or * D_a); + Ta -= C_b * pref * v11 * rxDa; - vterm = pref * ct_i * ( ri2 - preRF2_ * idat.rij ); - idat.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); - - } else { - - // determine inverse r if we are using split dipoles - if (i_is_SplitDipole) { - BigR = sqrt(idat.r2 + 0.25 * d_i * d_i); - ri = 1.0 / BigR; - scale = idat.rij * ri; - } else { - ri = riji; - scale = 1.0; - } - - sc2 = scale * scale; - - if (screeningMethod_ == DAMPED) { - // assemble the damping variables - res = erfcSpline_->getValueAndDerivativeAt(idat.rij); - erfcVal = res.first; - derfcVal = res.second; - c1 = erfcVal * ri; - c2 = (-derfcVal + c1) * ri; - c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * ri; - } else { - c1 = ri; - c2 = c1 * ri; - c3 = 3.0 * c2 * ri; - } - - c2ri = c2 * ri; - - // calculate the potential - pot_term = c2 * scale; - vterm = pref * ct_i * pot_term; - idat.vpair += vterm; - epot += idat.sw * vterm; + if (b_is_Fluctuating) dUdCb -= pref * v11 * rdDa; - // calculate derivatives for the forces and torques - dVdr += preSw * (uz_i * c2ri - ct_i * rhat * sc2 * c3); - duduz_i += preSw * pot_term * rhat; + // Even if we excluded this pair from direct interactions, + // we still have the reaction-field-mediated charge-dipole + // interaction: + if (summationMethod_ == esm_REACTION_FIELD && idat.excluded) { + rfContrib = C_b * pref * preRF_ * 2.0 * *(idat.rij); + indirect_Pot -= rfContrib * rdDa; + indirect_F -= rfContrib * D_a / (*idat.rij); + indirect_Ta -= C_b * pref * preRF_ * rxDa; } } - if (j_is_Dipole) { - // variables used by all methods - ct_ij = dot(uz_i, uz_j); + if (b_is_Dipole) { + pref = pre22_ * *(idat.electroMult); + DadDb = dot(D_a, D_b); + DaxDb = cross(D_a, D_b); - pref = idat.electroMult * pre22_ * mu_i * mu_j; - preSw = idat.sw * pref; + U -= pref * (DadDb * v21 + rdDa * rdDb * v22); + F -= pref * (dv21 * DadDb * rhat + v22or * (rdDb * D_a + rdDa * D_b)); + F -= pref * (rdDa * rdDb) * (dv22 - 2.0*v22or) * rhat; + Ta += pref * ( v21 * DaxDb - v22 * rdDb * rxDa); + Tb += pref * (-v21 * DaxDb - v22 * rdDa * rxDb); - if (summationMethod_ == REACTION_FIELD) { - ri2 = riji * riji; - ri3 = ri2 * riji; - ri4 = ri2 * ri2; + // Even if we excluded this pair from direct interactions, we + // still have the reaction-field-mediated dipole-dipole + // interaction: + if (summationMethod_ == esm_REACTION_FIELD && idat.excluded) { + rfContrib = -pref * preRF_ * 2.0; + indirect_Pot += rfContrib * DadDb; + indirect_Ta += rfContrib * DaxDb; + indirect_Tb -= rfContrib * DaxDb; + } + } - vterm = pref * ( ri3 * (ct_ij - 3.0 * ct_i * ct_j) - - preRF2_ * ct_ij ); - idat.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); + if (b_is_Quadrupole) { + pref = pre24_ * *(idat.electroMult); + DadQb = D_a * Q_b; + DadQbr = dot(D_a, Qbr); + DaxQbr = cross(D_a, Qbr); - 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); + U -= pref * ((trQb*rdDa + 2.0*DadQbr)*v31 + rdDa*rdQbr*v32); + F -= pref * (trQb*D_a + 2.0*DadQb) * v31or; + F -= pref * (trQb*rdDa + 2.0*DadQbr) * (dv31-v31or) * rhat; + F -= pref * (D_a*rdQbr + 2.0*rdDa*rQb) * v32or; + F -= pref * (rdDa * rdQbr * rhat * (dv32-3.0*v32or)); + Ta += pref * ((-trQb*rxDa + 2.0 * DaxQbr)*v31 - rxDa*rdQbr*v32); + Tb += pref * ((2.0*cross(DadQb, rhat) - 2.0*DaxQbr)*v31 + - 2.0*rdDa*rxQbr*v32); + } + } - } 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); - } else { - BigR = sqrt(idat.r2 + 0.25 * d_i * d_i); - } - ri = 1.0 / BigR; - scale = idat.rij * ri; - } else { - if (j_is_SplitDipole) { - BigR = sqrt(idat.r2 + 0.25 * d_j * d_j); - ri = 1.0 / BigR; - scale = idat.rij * ri; - } else { - ri = riji; - scale = 1.0; - } - } - if (screeningMethod_ == DAMPED) { - // assemble damping variables - res = erfcSpline_->getValueAndDerivativeAt(idat.rij); - erfcVal = res.first; - derfcVal = res.second; - c1 = erfcVal * ri; - c2 = (-derfcVal + c1) * ri; - c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * ri; - c4 = -4.0 * derfcVal * alpha4_ + 5.0 * c3 * ri * ri; - } else { - c1 = ri; - c2 = c1 * ri; - c3 = 3.0 * c2 * ri; - c4 = 5.0 * c3 * ri * ri; - } + if (a_is_Quadrupole) { + if (b_is_Charge) { + pref = pre14_ * *(idat.electroMult); + U += C_b * pref * (v21 * trQa + v22 * rdQar); + F += C_b * pref * (trQa * rhat * dv21 + 2.0 * Qar * v22or); + F += C_b * pref * rdQar * rhat * (dv22 - 2.0*v22or); + Ta += C_b * pref * 2.0 * rxQar * v22; - // precompute variables for convenience - sc2 = scale * scale; - cti3 = ct_i * sc2 * c3; - ctj3 = ct_j * sc2 * c3; - ctidotj = ct_i * ct_j * sc2; - preSwSc = preSw * scale; - c2ri = c2 * ri; - c3ri = c3 * ri; - c4rij = c4 * idat.rij; + if (b_is_Fluctuating) dUdCb += pref * (v21 * trQa + v22 * rdQar); + } + if (b_is_Dipole) { + pref = pre24_ * *(idat.electroMult); + DbdQa = D_b * Q_a; + DbdQar = dot(D_b, Qar); + DbxQar = cross(D_b, Qar); - // calculate the potential - pot_term = (ct_ij * c2ri - ctidotj * c3); - vterm = pref * pot_term; - idat.vpair += vterm; - epot += idat.sw * vterm; - - // calculate derivatives for the forces and torques - dVdr += preSwSc * ( ctidotj * rhat * c4rij - - (ct_i*uz_j + ct_j*uz_i + ct_ij*rhat) * c3ri); - - duduz_i += preSw * (uz_j * c2ri - ctj3 * rhat); - duduz_j += preSw * (uz_i * c2ri - cti3 * rhat); - } + U += pref * ((trQa*rdDb + 2.0*DbdQar)*v31 + rdDb*rdQar*v32); + F += pref * (trQa*D_b + 2.0*DbdQa) * v31or; + F += pref * (trQa*rdDb + 2.0*DbdQar) * (dv31-v31or) * rhat; + F += pref * (D_b*rdQar + 2.0*rdDb*rQa) * v32or; + F += pref * (rdDb * rdQar * rhat * (dv32-3.0*v32or)); + Ta += pref * ((-2.0*cross(DbdQa, rhat) + 2.0*DbxQar)*v31 + + 2.0*rdDb*rxQar*v32); + Tb += pref * ((trQa*rxDb - 2.0 * DbxQar)*v31 + rxDb*rdQar*v32); } - } + if (b_is_Quadrupole) { + pref = pre44_ * *(idat.electroMult); // yes + QaQb = Q_a * Q_b; + trQaQb = QaQb.trace(); + rQaQb = rhat * QaQb; + QaQbr = QaQb * rhat; + QaxQb = cross(Q_a, Q_b); + rQaQbr = dot(rQa, Qbr); + rQaxQbr = cross(rQa, Qbr); + + U += pref * (trQa * trQb + 2.0 * trQaQb) * v41; + U += pref * (trQa * rdQbr + trQb * rdQar + 4.0 * rQaQbr) * v42; + U += pref * (rdQar * rdQbr) * v43; - if (i_is_Quadrupole) { - if (j_is_Charge) { - // precompute some necessary variables - cx2 = cx_i * cx_i; - cy2 = cy_i * cy_i; - cz2 = cz_i * cz_i; + F += pref * rhat * (trQa * trQb + 2.0 * trQaQb)*dv41; + F += pref*rhat*(trQa*rdQbr + trQb*rdQar + 4.0*rQaQbr)*(dv42-2.0*v42or); + F += pref * rhat * (rdQar * rdQbr)*(dv43 - 4.0*v43or); - pref = idat.electroMult * pre14_ * q_j * one_third_; + F += pref * 2.0 * (trQb*rQa + trQa*rQb) * v42or; + F += pref * 4.0 * (rQaQb + QaQbr) * v42or; + F += pref * 2.0 * (rQa*rdQbr + rdQar*rQb) * v43or; - if (screeningMethod_ == DAMPED) { - // assemble the damping variables - res = erfcSpline_->getValueAndDerivativeAt(idat.rij); - erfcVal = res.first; - derfcVal = res.second; - c1 = erfcVal * riji; - c2 = (-derfcVal + c1) * riji; - c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * riji; - c4 = -4.0 * derfcVal * alpha4_ + 5.0 * c3 * riji * riji; - } else { - c1 = riji; - c2 = c1 * riji; - c3 = 3.0 * c2 * riji; - c4 = 5.0 * c3 * riji * riji; - } - - // precompute some variables for convenience - preSw = idat.sw * pref; - c2ri = c2 * riji; - c3ri = c3 * riji; - c4rij = c4 * idat.rij; - rhatdot2 = 2.0 * rhat * c3; - rhatc4 = rhat * c4rij; + Ta += pref * (- 4.0 * QaxQb * v41); + Ta += pref * (- 2.0 * trQb * cross(rQa, rhat) + + 4.0 * cross(rhat, QaQbr) + - 4.0 * rQaxQbr) * v42; + Ta += pref * 2.0 * cross(rhat,Qar) * rdQbr * v43; - // calculate the potential - pot_term = ( qxx_i * (cx2 * c3 - c2ri) + - qyy_i * (cy2 * c3 - c2ri) + - qzz_i * (cz2 * c3 - c2ri) ); - - vterm = pref * pot_term; - idat.vpair += vterm; - epot += idat.sw * vterm; - // calculate the derivatives for the forces and torques + Tb += pref * (+ 4.0 * QaxQb * v41); + Tb += pref * (- 2.0 * trQa * cross(rQb, rhat) + - 4.0 * cross(rQaQb, rhat) + + 4.0 * rQaxQbr) * v42; + // Possible replacement for line 2 above: + // + 4.0 * cross(rhat, QbQar) - 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)); + Tb += pref * 2.0 * cross(rhat,Qbr) * rdQar * v43; - dudux_i += preSw * qxx_i * cx_i * rhatdot2; - duduy_i += preSw * qyy_i * cy_i * rhatdot2; - duduz_i += preSw * qzz_i * cz_i * rhatdot2; } } - 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.doElectricField) { + *(idat.eField1) += Ea * *(idat.electroMult); + *(idat.eField2) += Eb * *(idat.electroMult); } - 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 (a_is_Fluctuating) *(idat.dVdFQ1) += dUdCa * *(idat.sw); + if (b_is_Fluctuating) *(idat.dVdFQ2) += dUdCb * *(idat.sw); - return; - } + if (!idat.excluded) { + + *(idat.vpair) += U; + (*(idat.pot))[ELECTROSTATIC_FAMILY] += U * *(idat.sw); + *(idat.f1) += F * *(idat.sw); + + if (a_is_Dipole || a_is_Quadrupole) + *(idat.t1) += Ta * *(idat.sw); - void Electrostatic::calcSkipCorrection(SkipCorrectionData skdat) { + if (b_is_Dipole || b_is_Quadrupole) + *(idat.t2) += Tb * *(idat.sw); + + } else { - 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; + // only accumulate the forces and torques resulting from the + // indirect reaction field terms. - RealType q_i, q_j; + *(idat.vpair) += indirect_Pot; + (*(idat.excludedPot))[ELECTROSTATIC_FAMILY] += excluded_Pot; + (*(idat.pot))[ELECTROSTATIC_FAMILY] += *(idat.sw) * indirect_Pot; + *(idat.f1) += *(idat.sw) * indirect_F; + + if (a_is_Dipole || a_is_Quadrupole) + *(idat.t1) += *(idat.sw) * indirect_Ta; + + if (b_is_Dipole || b_is_Quadrupole) + *(idat.t2) += *(idat.sw) * indirect_Tb; + } + return; + } - // The skippedCharge computation is needed by the real-space cutoff methods - // (i.e. shifted force and shifted potential) + void Electrostatic::calcSelfCorrection(SelfData &sdat) { - if (i_is_Charge) { - q_i = data1.charge; - skdat.skippedCharge2 += q_i; - } + if (!initialized_) initialize(); - if (j_is_Charge) { - q_j = data2.charge; - skdat.skippedCharge1 += q_j; + ElectrostaticAtomData data = ElectrostaticMap[Etids[sdat.atid]]; + + // logicals + bool i_is_Charge = data.is_Charge; + bool i_is_Dipole = data.is_Dipole; + bool i_is_Quadrupole = data.is_Quadrupole; + bool i_is_Fluctuating = data.is_Fluctuating; + RealType C_a = data.fixedCharge; + RealType self(0.0), preVal, DdD, trQ, trQQ; + + if (i_is_Dipole) { + DdD = data.dipole.lengthSquare(); } + + if (i_is_Fluctuating) { + C_a += *(sdat.flucQ); + // dVdFQ is really a force, so this is negative the derivative + *(sdat.dVdFQ) -= *(sdat.flucQ) * data.hardness + data.electronegativity; + (*(sdat.excludedPot))[ELECTROSTATIC_FAMILY] += (*sdat.flucQ) * + (*(sdat.flucQ) * data.hardness * 0.5 + data.electronegativity); + } - // 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: + switch (summationMethod_) { + case esm_REACTION_FIELD: - 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); - } + // Self potential [see Wang and Hermans, "Reaction Field + // Molecular Dynamics Simulation with Friedman’s Image Charge + // Method," J. Phys. Chem. 99, 12001-12007 (1995).] + preVal = pre11_ * preRF_ * C_a * C_a; + (*(sdat.pot))[ELECTROSTATIC_FAMILY] -= 0.5 * preVal / cutoffRadius_; } + 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); - } + (*(sdat.pot))[ELECTROSTATIC_FAMILY] -= pre22_ * preRF_ * DdD; } - // accumulate the forces and torques resulting from the self term - skdat.pot += myPot; - skdat.f1 += dVdr; + break; + case esm_SHIFTED_FORCE: + case esm_SHIFTED_POTENTIAL: + case esm_TAYLOR_SHIFTED: + case esm_EWALD_FULL: + if (i_is_Charge) + self += selfMult1_ * pre11_ * C_a * (C_a + *(sdat.skippedCharge)); if (i_is_Dipole) - *(skdat.t1) -= cross(uz_i, duduz_i); - if (j_is_Dipole) - *(skdat.t2) -= cross(uz_j, duduz_j); + self += selfMult2_ * pre22_ * DdD; + if (i_is_Quadrupole) { + trQ = data.quadrupole.trace(); + trQQ = (data.quadrupole * data.quadrupole).trace(); + self += selfMult4_ * pre44_ * (2.0*trQQ + trQ*trQ); + if (i_is_Charge) + self -= selfMult2_ * pre14_ * 2.0 * C_a * trQ; + } + (*(sdat.pot))[ELECTROSTATIC_FAMILY] += self; + break; + default: + break; } } + + 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; + } + + + void Electrostatic::ReciprocalSpaceSum(RealType& pot) { - void Electrostatic::calcSelfCorrection(SelfCorrectionData scdat) { - RealType mu1, preVal, chg1, self; + RealType kPot = 0.0; + RealType kVir = 0.0; - if (!initialized_) initialize(); + const RealType mPoleConverter = 0.20819434; // converts from the + // internal units of + // Debye (for dipoles) + // or Debye-angstroms + // (for quadrupoles) to + // electron angstroms or + // electron-angstroms^2 - ElectrostaticAtomData data = ElectrostaticMap[scdat.atype]; - - // logicals + const RealType eConverter = 332.0637778; // convert the + // Charge-Charge + // electrostatic + // interactions into kcal / + // mol assuming distances + // are measured in + // angstroms. - bool i_is_Charge = data.is_Charge; - bool i_is_Dipole = data.is_Dipole; + Mat3x3d hmat = info_->getSnapshotManager()->getCurrentSnapshot()->getHmat(); + Vector3d box = hmat.diagonals(); + RealType boxMax = box.max(); + + //int kMax = int(2.0 * M_PI / (pow(dampingAlpha_,2)*cutoffRadius_ * boxMax) ); + int kMax = 7; + int kSqMax = kMax*kMax + 2; + + int kLimit = kMax+1; + int kLim2 = 2*kMax+1; + int kSqLim = kSqMax; + + vector AK(kSqLim+1, 0.0); + RealType xcl = 2.0 * M_PI / box.x(); + RealType ycl = 2.0 * M_PI / box.y(); + RealType zcl = 2.0 * M_PI / box.z(); + RealType rcl = 2.0 * M_PI / boxMax; + RealType rvol = 2.0 * M_PI /(box.x() * box.y() * box.z()); + + if(dampingAlpha_ < 1.0e-12) return; + + RealType ralph = -0.25/pow(dampingAlpha_,2); + + // Calculate and store exponential factors + + vector > elc; + vector > emc; + vector > enc; + vector > els; + vector > ems; + vector > ens; - if (summationMethod_ == REACTION_FIELD) { - if (i_is_Dipole) { - mu1 = data.dipole_moment; - preVal = pre22_ * preRF2_ * mu1 * mu1; - scdat.pot -= 0.5 * preVal; - - // The self-correction term adds into the reaction field vector - Vector3d uz_i = scdat.eFrame->getColumn(2); - Vector3d ei = preVal * uz_i; + + int nMax = info_->getNAtoms(); + + elc.resize(kLimit+1); + emc.resize(kLimit+1); + enc.resize(kLimit+1); + els.resize(kLimit+1); + ems.resize(kLimit+1); + ens.resize(kLimit+1); - // This looks very wrong. A vector crossed with itself is zero. - *(scdat.t) -= cross(uz_i, ei); - } - } else if (summationMethod_ == SHIFTED_FORCE || summationMethod_ == SHIFTED_POTENTIAL) { - if (i_is_Charge) { - chg1 = data.charge; - if (screeningMethod_ == DAMPED) { - self = - 0.5 * (c1c_ + alphaPi_) * chg1 * (chg1 + scdat.skippedCharge) * pre11_; - } else { - self = - 0.5 * rcuti_ * chg1 * (chg1 + scdat.skippedCharge) * pre11_; + for (int j = 0; j < kLimit+1; j++) { + elc[j].resize(nMax); + emc[j].resize(nMax); + enc[j].resize(nMax); + els[j].resize(nMax); + ems[j].resize(nMax); + ens[j].resize(nMax); + } + + Vector3d t( 2.0 * M_PI ); + t.Vdiv(t, box); + + + SimInfo::MoleculeIterator mi; + Molecule::AtomIterator ai; + int i; + Vector3d r; + Vector3d tt; + + for (Molecule* mol = info_->beginMolecule(mi); mol != NULL; + mol = info_->nextMolecule(mi)) { + for(Atom* atom = mol->beginAtom(ai); atom != NULL; + atom = mol->nextAtom(ai)) { + + i = atom->getLocalIndex(); + r = atom->getPos(); + info_->getSnapshotManager()->getCurrentSnapshot()->wrapVector(r); + + tt.Vmul(t, r); + + elc[1][i] = 1.0; + emc[1][i] = 1.0; + enc[1][i] = 1.0; + els[1][i] = 0.0; + ems[1][i] = 0.0; + ens[1][i] = 0.0; + + elc[2][i] = cos(tt.x()); + emc[2][i] = cos(tt.y()); + enc[2][i] = cos(tt.z()); + els[2][i] = sin(tt.x()); + ems[2][i] = sin(tt.y()); + ens[2][i] = sin(tt.z()); + + for(int l = 3; l <= kLimit; l++) { + elc[l][i]=elc[l-1][i]*elc[2][i]-els[l-1][i]*els[2][i]; + emc[l][i]=emc[l-1][i]*emc[2][i]-ems[l-1][i]*ems[2][i]; + enc[l][i]=enc[l-1][i]*enc[2][i]-ens[l-1][i]*ens[2][i]; + els[l][i]=els[l-1][i]*elc[2][i]+elc[l-1][i]*els[2][i]; + ems[l][i]=ems[l-1][i]*emc[2][i]+emc[l-1][i]*ems[2][i]; + ens[l][i]=ens[l-1][i]*enc[2][i]+enc[l-1][i]*ens[2][i]; } - scdat.pot += self; } } + + // Calculate and store AK coefficients: + + RealType eksq = 1.0; + RealType expf = 0.0; + if (ralph < 0.0) expf = exp(ralph*rcl*rcl); + for (i = 1; i <= kSqLim; i++) { + RealType rksq = float(i)*rcl*rcl; + eksq = expf*eksq; + AK[i] = eConverter * eksq/rksq; + } + + /* + * Loop over all k vectors k = 2 pi (ll/Lx, mm/Ly, nn/Lz) + * the values of ll, mm and nn are selected so that the symmetry of + * reciprocal lattice is taken into account i.e. the following + * rules apply. + * + * ll ranges over the values 0 to kMax only. + * + * mm ranges over 0 to kMax when ll=0 and over + * -kMax to kMax otherwise. + * nn ranges over 1 to kMax when ll=mm=0 and over + * -kMax to kMax otherwise. + * + * Hence the result of the summation must be doubled at the end. + */ + + std::vector clm(nMax, 0.0); + std::vector slm(nMax, 0.0); + std::vector ckr(nMax, 0.0); + std::vector skr(nMax, 0.0); + std::vector ckc(nMax, 0.0); + std::vector cks(nMax, 0.0); + std::vector dkc(nMax, 0.0); + std::vector dks(nMax, 0.0); + std::vector qkc(nMax, 0.0); + std::vector qks(nMax, 0.0); + std::vector dxk(nMax, V3Zero); + std::vector qxk(nMax, V3Zero); + RealType rl, rm, rn; + Vector3d kVec; + Vector3d Qk; + Mat3x3d k2; + RealType ckcs, ckss, dkcs, dkss, qkcs, qkss; + int atid; + ElectrostaticAtomData data; + RealType C, dk, qk; + Vector3d D; + Mat3x3d Q; + + int mMin = kLimit; + int nMin = kLimit + 1; + for (int l = 1; l <= kLimit; l++) { + int ll = l - 1; + rl = xcl * float(ll); + for (int mmm = mMin; mmm <= kLim2; mmm++) { + int mm = mmm - kLimit; + int m = abs(mm) + 1; + rm = ycl * float(mm); + // Set temporary products of exponential terms + for (Molecule* mol = info_->beginMolecule(mi); mol != NULL; + mol = info_->nextMolecule(mi)) { + for(Atom* atom = mol->beginAtom(ai); atom != NULL; + atom = mol->nextAtom(ai)) { + + i = atom->getLocalIndex(); + if(mm < 0) { + clm[i]=elc[l][i]*emc[m][i]+els[l][i]*ems[m][i]; + slm[i]=els[l][i]*emc[m][i]-ems[m][i]*elc[l][i]; + } else { + clm[i]=elc[l][i]*emc[m][i]-els[l][i]*ems[m][i]; + slm[i]=els[l][i]*emc[m][i]+ems[m][i]*elc[l][i]; + } + } + } + for (int nnn = nMin; nnn <= kLim2; nnn++) { + int nn = nnn - kLimit; + int n = abs(nn) + 1; + rn = zcl * float(nn); + // Test on magnitude of k vector: + int kk=ll*ll + mm*mm + nn*nn; + if(kk <= kSqLim) { + kVec = Vector3d(rl, rm, rn); + k2 = outProduct(kVec, kVec); + // Calculate exp(ikr) terms + for (Molecule* mol = info_->beginMolecule(mi); mol != NULL; + mol = info_->nextMolecule(mi)) { + for(Atom* atom = mol->beginAtom(ai); atom != NULL; + atom = mol->nextAtom(ai)) { + i = atom->getLocalIndex(); + + if (nn < 0) { + ckr[i]=clm[i]*enc[n][i]+slm[i]*ens[n][i]; + skr[i]=slm[i]*enc[n][i]-clm[i]*ens[n][i]; + + } else { + ckr[i]=clm[i]*enc[n][i]-slm[i]*ens[n][i]; + skr[i]=slm[i]*enc[n][i]+clm[i]*ens[n][i]; + } + } + } + + // Calculate scalar and vector products for each site: + + for (Molecule* mol = info_->beginMolecule(mi); mol != NULL; + mol = info_->nextMolecule(mi)) { + for(Atom* atom = mol->beginAtom(ai); atom != NULL; + atom = mol->nextAtom(ai)) { + i = atom->getLocalIndex(); + int atid = atom->getAtomType()->getIdent(); + data = ElectrostaticMap[Etids[atid]]; + + if (data.is_Charge) { + C = data.fixedCharge; + if (atom->isFluctuatingCharge()) C += atom->getFlucQPos(); + ckc[i] = C * ckr[i]; + cks[i] = C * skr[i]; + } + + if (data.is_Dipole) { + D = atom->getDipole() * mPoleConverter; + dk = dot(D, kVec); + dxk[i] = cross(D, kVec); + dkc[i] = dk * ckr[i]; + dks[i] = dk * skr[i]; + } + if (data.is_Quadrupole) { + Q = atom->getQuadrupole(); + Q *= mPoleConverter; + Qk = Q * kVec; + qk = dot(kVec, Qk); + qxk[i] = cross(kVec, Qk); + qkc[i] = qk * ckr[i]; + qks[i] = qk * skr[i]; + } + } + } + + // calculate vector sums + + ckcs = std::accumulate(ckc.begin(),ckc.end(),0.0); + ckss = std::accumulate(cks.begin(),cks.end(),0.0); + dkcs = std::accumulate(dkc.begin(),dkc.end(),0.0); + dkss = std::accumulate(dks.begin(),dks.end(),0.0); + qkcs = std::accumulate(qkc.begin(),qkc.end(),0.0); + qkss = std::accumulate(qks.begin(),qks.end(),0.0); + +#ifdef IS_MPI + MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &ckcs, 1, MPI::REALTYPE, + MPI::SUM); + MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &ckss, 1, MPI::REALTYPE, + MPI::SUM); + MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &dkcs, 1, MPI::REALTYPE, + MPI::SUM); + MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &dkss, 1, MPI::REALTYPE, + MPI::SUM); + MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &qkcs, 1, MPI::REALTYPE, + MPI::SUM); + MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &qkss, 1, MPI::REALTYPE, + MPI::SUM); +#endif + + // Accumulate potential energy and virial contribution: + + kPot += 2.0 * rvol * AK[kk]*((ckss+dkcs-qkss)*(ckss+dkcs-qkss) + + (ckcs-dkss-qkcs)*(ckcs-dkss-qkcs)); + + kVir += 2.0 * rvol * AK[kk]*(ckcs*ckcs+ckss*ckss + +4.0*(ckss*dkcs-ckcs*dkss) + +3.0*(dkcs*dkcs+dkss*dkss) + -6.0*(ckss*qkss+ckcs*qkcs) + +8.0*(dkss*qkcs-dkcs*qkss) + +5.0*(qkss*qkss+qkcs*qkcs)); + + // Calculate force and torque for each site: + + for (Molecule* mol = info_->beginMolecule(mi); mol != NULL; + mol = info_->nextMolecule(mi)) { + for(Atom* atom = mol->beginAtom(ai); atom != NULL; + atom = mol->nextAtom(ai)) { + + i = atom->getLocalIndex(); + atid = atom->getAtomType()->getIdent(); + data = ElectrostaticMap[Etids[atid]]; + + RealType qfrc = AK[kk]*((cks[i]+dkc[i]-qks[i])*(ckcs-dkss-qkcs) + - (ckc[i]-dks[i]-qkc[i])*(ckss+dkcs-qkss)); + RealType qtrq1 = AK[kk]*(skr[i]*(ckcs-dkss-qkcs) + -ckr[i]*(ckss+dkcs-qkss)); + RealType qtrq2 = 2.0*AK[kk]*(ckr[i]*(ckcs-dkss-qkcs) + +skr[i]*(ckss+dkcs-qkss)); + + atom->addFrc( 4.0 * rvol * qfrc * kVec ); + + if (data.is_Dipole) { + atom->addTrq( 4.0 * rvol * qtrq1 * dxk[i] ); + } + if (data.is_Quadrupole) { + atom->addTrq( 4.0 * rvol * qtrq2 * qxk[i] ); + } + } + } + } + } + nMin = 1; + } + mMin = 1; + } + pot += kPot; } }