--- branches/development/src/nonbonded/Electrostatic.cpp 2010/10/02 19:53:32 1502 +++ branches/development/src/nonbonded/Electrostatic.cpp 2012/08/29 18:13:11 1787 @@ -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,15 +47,41 @@ #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" - 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_["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,13 +89,18 @@ 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; @@ -79,102 +111,453 @@ namespace OpenMD { // variables to handle different summation methods for long-range // electrostatics: - summationMethod_ = NONE; + summationMethod_ = esm_HARD; screeningMethod_ = UNDAMPED; dielectric_ = 1.0; - one_third_ = 1.0 / 3.0; - haveDefaultCutoff_ = false; - haveDampingAlpha_ = false; - haveDielectric_ = false; - haveElectroSpline_ = false; - // 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\", or \n" + "\t\"reaction_field\".\n", myMethod.c_str() ); + painCave.isFatal = 1; + simError(); + } + } else { + // set ElectrostaticSummationMethod to the cutoffMethod: + if (simParams_->haveCutoffMethod()){ + string myMethod = simParams_->getCutoffMethod(); + toUpper(myMethod); + map::iterator i; + i = summationMap_.find(myMethod); + if ( i != summationMap_.end() ) { + summationMethod_ = (*i).second; + } + } } + if (summationMethod_ == esm_REACTION_FIELD) { + if (!simParams_->haveDielectric()) { + // throw warning + sprintf( painCave.errMsg, + "SimInfo warning: dielectric was not specified in the input file\n\tfor the reaction field correction method.\n" + "\tA default value of %f will be used for the dielectric.\n", dielectric_); + painCave.isFatal = 0; + painCave.severity = OPENMD_INFO; + simError(); + } else { + dielectric_ = simParams_->getDielectric(); + } + haveDielectric_ = true; + } + + if (simParams_->haveElectrostaticScreeningMethod()) { + string myScreen = simParams_->getElectrostaticScreeningMethod(); + toUpper(myScreen); + map::iterator i; + i = screeningMap_.find(myScreen); + if ( i != screeningMap_.end()) { + screeningMethod_ = (*i).second; + } else { + sprintf( painCave.errMsg, + "SimInfo error: Unknown electrostaticScreeningMethod.\n" + "\t(Input file specified %s .)\n" + "\telectrostaticScreeningMethod must be one of: \"undamped\"\n" + "or \"damped\".\n", myScreen.c_str() ); + painCave.isFatal = 1; + simError(); + } + } + // check to make sure a cutoff value has been set: - if (!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) { + 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_; + // 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)) { - constEXP_ = exp(-alpha2_ * defaultCutoff2_); - invRootPi_ = 0.56418958354775628695; - alphaPi_ = 2.0 * dampingAlpha_ * invRootPi_; + 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; - 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_; + 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; + selfMult_ = b0c + 2.0 * a2 * invArootPi; } 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_; - } - - if (summationMethod_ == REACTION_FIELD) { - if (haveDielectric_) { - preRF_ = (dielectric_ - 1.0) / - ((2.0 * dielectric_ + 1.0) * defaultCutoff2_ * defaultCutoff_); - preRF2_ = 2.0 * preRF_; + 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; + selfMult_ = b0c; + } + + // 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; + + // 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 f0; + RealType g0, g1, g2, g3, g4; + RealType h1, h2, h3, h4; + RealType s2, s3, s4; + RealType t3, t4; + RealType u4; + + // working variables for Taylor expansion: + RealType rmRc, rmRc2, rmRc3, rmRc4; + + // 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, v02v; + vector v11v, v12v, v13v; + vector v21v, v22v, v23v, v24v; + vector v31v, v32v, v33v, v34v, v35v; + vector v41v, v42v, v43v, v44v, v45v, v46v; + + 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; + + + switch (summationMethod_) { + case esm_SHIFTED_FORCE: + f0 = b0 - b0c - rmRc*db0c_1; + + g0 = db0_1 - db0c_1; + g1 = g0 - rmRc *db0c_2; + g2 = g1 - rmRc2*db0c_3; + g3 = g2 - rmRc3*db0c_4; + g4 = g3 - rmRc4*db0c_5; + + h1 = db0_2 - db0c_2; + h2 = h1 - rmRc *db0c_3; + h3 = h2 - rmRc2*db0c_4; + h4 = h3 - rmRc3*db0c_5; + + s2 = db0_3 - db0c_3; + s3 = s2 - rmRc *db0c_4; + s4 = s3 - rmRc2*db0c_5; + + t3 = db0_4 - db0c_4; + t4 = t3 - rmRc *db0c_5; + + u4 = db0_5 - db0c_5; + break; + + case esm_SHIFTED_POTENTIAL: + f0 = b0 - b0c; + + g0 = db0_1; + g1 = db0_1 - db0c_1; + g2 = g1 - rmRc *db0c_2; + g3 = g2 - rmRc2*db0c_3; + g4 = g3 - rmRc3*db0c_4; + + h1 = db0_2; + h2 = db0_2 - db0c_2; + h3 = h2 - rmRc *db0c_3; + h4 = h3 - rmRc2*db0c_4; + + s2 = db0_3; + s3 = db0_3 - db0c_3; + s4 = s3 - rmRc *db0c_4; + + t3 = db0_4; + t4 = db0_4 - db0c_4; + + u4 = db0_5; + break; + + case esm_SWITCHING_FUNCTION: + case esm_HARD: + f0 = b0; + + g0 = db0_1; + g1 = g0; + g2 = g1; + g3 = g2; + g4 = g3; + + h1 = db0_2; + h2 = h1; + h3 = h2; + h4 = h3; + + s2 = db0_3; + s3 = s2; + s4 = s3; + + t3 = db0_4; + t4 = t3; + + u4 = db0_5; + break; + + case esm_REACTION_FIELD: + + // following DL_POLY's lead for shifting the image charge potential: + f0 = b0 + preRF_ * r2 + - (b0c + preRF_ * cutoffRadius_ * cutoffRadius_); + + g0 = db0_1 + preRF_ * 2.0 * r; + g1 = g0; + g2 = g1; + g3 = g2; + g4 = g3; + + h1 = db0_2 + preRF_ * 2.0; + h2 = h1; + h3 = h2; + h4 = h3; + + s2 = db0_3; + s3 = s2; + s4 = s3; + + t3 = db0_4; + t4 = t3; + + u4 = db0_5; + break; + + case esm_EWALD_FULL: + 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; } - } - - 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)); + + v01 = f0; + v02 = g0; + + v11 = g1; + v12 = g1 * ri; + v13 = h1 - v12; + + v21 = g2 * ri; + v22 = h2 - v21; + v23 = v22 * ri; + v24 = s2 - 3.0*v23; + + v31 = (h3 - g3 * ri) * ri; + v32 = s3 - 3.0*v31; + v33 = v31 * ri; + v34 = v32 * ri; + v35 = t3 - 6.0*v34 - 3.0*v33; + + v41 = (h4 - g4 * ri) * ri2; + v42 = s4 * ri - 3.0*v41; + v43 = t4 - 6.0*v42 - 3.0*v41; + v44 = v42 * ri; + v45 = v43 * ri; + v46 = u4 - 10.0*v45 - 15.0*v44; + + // Add these computed values to the storage vectors for spline creation: + v01v.push_back(v01); + v02v.push_back(v02); + + v11v.push_back(v11); + v12v.push_back(v12); + v13v.push_back(v13); + + v21v.push_back(v21); + v22v.push_back(v22); + v23v.push_back(v23); + v24v.push_back(v24); + + v31v.push_back(v31); + v32v.push_back(v32); + v33v.push_back(v33); + v34v.push_back(v34); + v35v.push_back(v35); + + v41v.push_back(v41); + v42v.push_back(v42); + v43v.push_back(v43); + v44v.push_back(v44); + v45v.push_back(v45); + v46v.push_back(v46); } - 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); + v02s = new CubicSpline(); + v02s->addPoints(rv, v02v); + + v11s = new CubicSpline(); + v11s->addPoints(rv, v11v); + v12s = new CubicSpline(); + v12s->addPoints(rv, v12v); + v13s = new CubicSpline(); + v13s->addPoints(rv, v13v); + + v21s = new CubicSpline(); + v21s->addPoints(rv, v21v); + v22s = new CubicSpline(); + v22s->addPoints(rv, v22v); + v23s = new CubicSpline(); + v23s->addPoints(rv, v23v); + v24s = new CubicSpline(); + v24s->addPoints(rv, v24v); + + v31s = new CubicSpline(); + v31s->addPoints(rv, v31v); + v32s = new CubicSpline(); + v32s->addPoints(rv, v32v); + v33s = new CubicSpline(); + v33s->addPoints(rv, v33v); + v34s = new CubicSpline(); + v34s->addPoints(rv, v34v); + v35s = new CubicSpline(); + v35s->addPoints(rv, v35v); + + v41s = new CubicSpline(); + v41s->addPoints(rv, v41v); + v42s = new CubicSpline(); + v42s->addPoints(rv, v42v); + v43s = new CubicSpline(); + v43s->addPoints(rv, v43v); + v44s = new CubicSpline(); + v44s->addPoints(rv, v44v); + v45s = new CubicSpline(); + v45s->addPoints(rv, v45v); + v46s = new CubicSpline(); + v46s->addPoints(rv, v46v); + + haveElectroSplines_ = true; + initialized_ = true; } @@ -183,145 +566,104 @@ namespace OpenMD { 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); + 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 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); + + pair key1, key2; + key1 = make_pair(atomType, atype2); + key2 = make_pair(atype2, atomType); + + Jij[key1] = J; + Jij[key2] = 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 +679,454 @@ 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 C_a, C_b; // Charges + Vector3d D_a, D_b; // Dipoles (space-fixed) + Mat3x3d Q_a, Q_b; // Quadrupoles (space-fixed) - 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; + RealType ri, ri2, ri3, ri4; // Distance utility scalars + RealType rdDa, rdDb; // Dipole utility scalars + Vector3d rxDa, rxDb; // Dipole utility vectors + RealType rdQar, rdQbr, trQa, trQb; // Quadrupole utility scalars + Vector3d Qar, Qbr, rQa, rQb, rxQar, rxQbr; // Quadrupole utility vectors + RealType pref; - 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; + RealType DadDb, trQaQb, DadQbr, DbdQar; // Cross-interaction scalars + Vector3d DaxDb, DadQb, DbdQa, DaxQbr, DbxQar; // Cross-interaction vectors + Vector3d rQaQb, QaQbr, QaxQb; + Mat3x3d QaQb; // Cross-interaction matrices - pair res; + RealType U(0.0); // Potential + Vector3d F(0.0); // Force + Vector3d Ta(0.0); // Torque on site a + Vector3d Tb(0.0); // Torque on site b + Vector3d Ea(0.0); // Electric field at site a + Vector3d Eb(0.0); // Electric field at site b + RealType dUdCa(0.0); // fluctuating charge force at site a + RealType dUdCb(0.0); // fluctuating charge force at site a + // Indirect interactions mediated by the reaction field. + RealType indirect_Pot(0.0); // Potential + Vector3d indirect_F(0.0); // Force + Vector3d indirect_Ta(0.0); // Torque on site a + Vector3d indirect_Tb(0.0); // Torque on site b + + // Excluded potential that is still computed for fluctuating charges + RealType excluded_Pot(0.0); + + RealType rfContrib, coulInt; + + // spline for coulomb integral + CubicSpline* J; + if (!initialized_) initialize(); - ElectrostaticAtomData data1 = ElectrostaticMap[idat.atype1]; - ElectrostaticAtomData data2 = ElectrostaticMap[idat.atype2]; + ElectrostaticAtomData data1 = ElectrostaticMap[idat.atypes.first]; + ElectrostaticAtomData data2 = ElectrostaticMap[idat.atypes.second]; // some variables we'll need independent of electrostatic type: - riji = 1.0 / idat.rij; - Vector3d rhat = idat.d * riji; - + ri = 1.0 / *(idat.rij); + Vector3d rhat = *(idat.d) * ri; + ri2 = ri * ri; + // 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; + bool a_is_Charge = data1.is_Charge; + bool a_is_Dipole = data1.is_Dipole; + bool a_is_Quadrupole = data1.is_Quadrupole; + bool a_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 b_is_Charge = data2.is_Charge; + bool b_is_Dipole = data2.is_Dipole; + bool b_is_Quadrupole = data2.is_Quadrupole; + bool b_is_Fluctuating = data2.is_Fluctuating; + + // Obtain all of the required radial function values from the + // spline structures: - if (i_is_Charge) - q_i = data1.charge; + if (a_is_Charge && b_is_Charge) { + v01 = v01s->getValueAt( *(idat.rij) ); + v02 = v02s->getValueAt( *(idat.rij) ); + } + if ((a_is_Charge && b_is_Dipole) || (b_is_Charge && a_is_Dipole)) { + v11 = v11s->getValueAt( *(idat.rij) ); + v12 = v12s->getValueAt( *(idat.rij) ); + v13 = v13s->getValueAt( *(idat.rij) ); + } + if ((a_is_Charge && b_is_Quadrupole) || + (b_is_Charge && a_is_Quadrupole) || + (a_is_Dipole && b_is_Dipole)) { + v21 = v21s->getValueAt( *(idat.rij) ); + v22 = v22s->getValueAt( *(idat.rij) ); + v23 = v23s->getValueAt( *(idat.rij) ); + v24 = v24s->getValueAt( *(idat.rij) ); + } + if ((a_is_Dipole && b_is_Quadrupole) || + (b_is_Dipole && a_is_Quadrupole)) { + v31 = v31s->getValueAt( *(idat.rij) ); + v32 = v32s->getValueAt( *(idat.rij) ); + v33 = v33s->getValueAt( *(idat.rij) ); + v34 = v34s->getValueAt( *(idat.rij) ); + v35 = v35s->getValueAt( *(idat.rij) ); + } + if (a_is_Quadrupole && b_is_Quadrupole) { + v41 = v41s->getValueAt( *(idat.rij) ); + v42 = v42s->getValueAt( *(idat.rij) ); + v43 = v43s->getValueAt( *(idat.rij) ); + v44 = v44s->getValueAt( *(idat.rij) ); + v45 = v45s->getValueAt( *(idat.rij) ); + v46 = v46s->getValueAt( *(idat.rij) ); + } - if (i_is_Dipole) { - mu_i = data1.dipole_moment; - uz_i = idat.eFrame1.getColumn(2); - - ct_i = dot(uz_i, rhat); - if (i_is_SplitDipole) - d_i = data1.split_dipole_distance; + // calculate the single-site contributions (fields, etc). + + if (a_is_Charge) { + C_a = data1.fixedCharge; - duduz_i = V3Zero; + if (a_is_Fluctuating) { + C_a += *(idat.flucQ1); + } + + if (idat.excluded) { + *(idat.skippedCharge2) += C_a; + } + Eb -= C_a * pre11_ * v02 * rhat; } - 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); - - cx_i = dot(ux_i, rhat); - cy_i = dot(uy_i, rhat); - cz_i = dot(uz_i, rhat); - - dudux_i = V3Zero; - duduy_i = V3Zero; - duduz_i = V3Zero; + if (a_is_Dipole) { + D_a = *(idat.dipole1); + rdDa = dot(rhat, D_a); + rxDa = cross(rhat, D_a); + Eb -= pre12_ * (v13 * rdDa * rhat + v12 * D_a); } - - if (j_is_Charge) - q_j = data2.charge; - - if (j_is_Dipole) { - mu_j = data2.dipole_moment; - uz_j = idat.eFrame2.getColumn(2); + + 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); + Eb -= pre14_ * ((trQa * rhat + 2.0 * Qar) * v23 + rdQar * rhat * v24); + } + + if (b_is_Charge) { + C_b = data2.fixedCharge; - ct_j = dot(uz_j, rhat); - - if (j_is_SplitDipole) - d_j = data2.split_dipole_distance; + if (b_is_Fluctuating) + C_b += *(idat.flucQ2); - duduz_j = V3Zero; + if (idat.excluded) { + *(idat.skippedCharge1) += C_b; + } + Ea += C_b * pre11_ * v02 * 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 (b_is_Dipole) { + D_b = *(idat.dipole2); + rdDb = dot(rhat, D_b); + rxDb = cross(rhat, D_b); + Ea += pre12_ * (v13 * rdDb * rhat + v12 * D_b); } - epot = 0.0; - dVdr = V3Zero; + 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); + Ea += pre14_ * ((trQb * rhat + 2.0 * Qbr) * v23 + rdQbr * rhat * v24); + } - if (i_is_Charge) { + if ((a_is_Fluctuating || b_is_Fluctuating) && idat.excluded) { + J = Jij[idat.atypes]; + } + + if (a_is_Charge) { - 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; - } + if (b_is_Charge) { + pref = pre11_ * *(idat.electroMult); + U += C_a * C_b * pref * v01; + F += C_a * C_b * pref * v02 * rhat; + + // If this is an excluded pair, there are still indirect + // interactions via the reaction field we must worry about: - preVal = idat.electroMult * pre11_ * q_i * q_j; + 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; + } - if (summationMethod_ == SHIFTED_POTENTIAL) { - vterm = preVal * (c1 - c1c_); - dudr = -idat.sw * preVal * c2; + // 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. - } else if (summationMethod_ == SHIFTED_FORCE) { - vterm = preVal * ( c1 - c1c_ + c2c_*(idat.rij - defaultCutoff_) ); - dudr = idat.sw * preVal * (c2c_ - c2); + 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 { + if (a_is_Fluctuating) dUdCa += C_b * pref * v01; + if (a_is_Fluctuating) dUdCb += C_a * pref * v01; + } + } - } 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; + if (b_is_Dipole) { + pref = pre12_ * *(idat.electroMult); + U += C_a * pref * v11 * rdDb; + F += C_a * pref * (v13 * rdDb * rhat + v12 * D_b); + Tb += C_a * pref * v11 * rxDb; - } else { - vterm = preVal * riji * erfcVal; + if (a_is_Fluctuating) dUdCa += pref * v11 * rdDb; - dudr = - idat.sw * preVal * c2; + // 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; } - - idat.vpair += vterm; - epot += idat.sw * vterm; - - dVdr += dudr * rhat; } - if (j_is_Dipole) { - // pref is used by all the possible methods - pref = idat.electroMult * pre12_ * q_i * mu_j; - preSw = idat.sw * pref; + if (b_is_Quadrupole) { + pref = pre14_ * *(idat.electroMult); + U += C_a * pref * (v21 * trQb + v22 * rdQbr); + F += C_a * pref * (trQb * rhat + 2.0 * Qbr) * v23; + F += C_a * pref * rdQbr * rhat * v24; + Tb += C_a * pref * 2.0 * rxQbr * v22; - 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 (a_is_Fluctuating) dUdCa += pref * (v21 * trQb + v22 * rdQbr); + } + } - dVdr += -preSw * (ri3 * (uz_j - 3.0 * ct_j * rhat) - preRF2_*uz_j); - duduz_j += -preSw * rhat * (ri2 - preRF2_ * idat.rij); + if (a_is_Dipole) { - } 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 (b_is_Charge) { + pref = pre12_ * *(idat.electroMult); - 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; + U -= C_b * pref * v11 * rdDa; + F -= C_b * pref * (v13 * rdDa * rhat + v12 * D_a); + Ta -= C_b * pref * v11 * rxDa; - // 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 (b_is_Fluctuating) dUdCb -= pref * v11 * rdDa; - 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_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_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_Dipole) { + pref = pre22_ * *(idat.electroMult); + DadDb = dot(D_a, D_b); + DaxDb = cross(D_a, D_b); + + U -= pref * (DadDb * v21 + rdDa * rdDb * v22); + F -= pref * (DadDb * rhat + rdDb * D_a + rdDa * D_b)*v23; + F -= pref * (rdDa * rdDb) * v24 * rhat; + Ta += pref * ( v21 * DaxDb - v22 * rdDb * rxDa); + Tb += pref * (-v21 * DaxDb + v22 * rdDa * rxDb); + + // 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; } - // 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 + if (b_is_Quadrupole) { + pref = pre24_ * *(idat.electroMult); + DadQb = D_a * Q_b; + DadQbr = dot(D_a, Qbr); + DaxQbr = cross(D_a, Qbr); - 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; + U -= pref * ((trQb*rdDa + 2.0*DadQbr)*v31 + rdDa*rdQbr*v32); + F -= pref * (trQb*D_a + 2.0*DadQb) * v33; + F -= pref * (trQb*rdDa*rhat + 2.0*DadQbr*rhat + D_a*rdQbr + + 2.0*rdDa*rQb)*v34; + F -= pref * (rdDa * rdQbr * rhat * v35); + 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); } } - - 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_Quadrupole) { + if (b_is_Charge) { + pref = pre14_ * *(idat.electroMult); + U += C_b * pref * (v21 * trQa + v22 * rdQar); + F += C_b * pref * (trQa * rhat + 2.0 * Qar) * v23; + F += C_b * pref * rdQar * rhat * v24; + Ta += C_b * pref * 2.0 * rxQar * v22; - if (summationMethod_ == REACTION_FIELD) { + 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); - ri2 = riji * riji; - ri3 = ri2 * riji; + U += pref * ((trQa*rdDb + 2.0*DbdQar)*v31 + rdDb*rdQar*v32); + F += pref * (trQa*D_b + 2.0*DbdQa) * v33; + F += pref * (trQa*rdDb*rhat + 2.0*DbdQar*rhat + D_b*rdQar + + 2.0*rdDb*rQa)*v34; + F += pref * (rdDb * rdQar * rhat * v35); + 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); + QaQb = Q_a * Q_b; + trQaQb = QaQb.trace(); + rQaQb = rhat * QaQb; + QaQbr = QaQb * rhat; + QaxQb = cross(Q_a, Q_b); - 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; + U += pref * (trQa * trQb + 2.0*trQaQb) * v41; + U += pref * (trQa*rdQbr + trQb*rdQar + 4.0*dot(rQa, Qbr)) * v42; + U += pref * (rdQar * rdQbr) * v43; - // calculate derivatives for the forces and torques - dVdr += preSw * (uz_i * c2ri - ct_i * rhat * sc2 * c3); - duduz_i += preSw * pot_term * rhat; - } - } + F += pref * (trQa*trQb*rhat + 2.0*trQaQb*rhat)*v44; + F += pref * (2.0*trQb*rQa + 2.0*trQa*rQb)*v44; + F += pref * (4.0* QaQb * rhat + 4.0 * rhat * QaQb)*v44; + F += pref * (trQa*rdQbr*rhat + trQb*rdQar*rhat + + 4.0*dot(rQa, Qbr)*rhat)*v45; + F += pref * (2.0*rQa*rdQbr + 2.0*rdQar*rQb)*v45; + F += pref * (rdQar*rdQbr*rhat) * v46; - if (j_is_Dipole) { - // variables used by all methods - ct_ij = dot(uz_i, uz_j); + Ta += pref * (-4.0 * QaxQb * v41); + Ta += pref * (-2.0*trQb*cross(rQa, rhat) + + 4.0*cross(rhat, QaQbr) + - 4.0*cross(rQa, Qbr)) * v42; + Ta += pref * 2.0 * cross(rhat,Qar) * rdQbr * v43; - pref = idat.electroMult * pre22_ * mu_i * mu_j; - preSw = idat.sw * pref; + Tb += pref * (4.0 * QaxQb * v41); + Tb += pref * (-2.0*trQa*cross(rQb, rhat) + - 4.0*cross(rQaQb, rhat) + + 4.0*cross(rQa, Qbr))*v42; + Tb += pref * 2.0 * cross(rhat,Qbr) * rdQar * v43; - if (summationMethod_ == REACTION_FIELD) { - ri2 = riji * riji; - ri3 = ri2 * riji; - ri4 = ri2 * ri2; - - vterm = pref * ( ri3 * (ct_ij - 3.0 * ct_i * ct_j) - - preRF2_ * ct_ij ); - idat.vpair += vterm; - epot += idat.sw * vterm; - - a1 = 5.0 * ct_i * ct_j - ct_ij; - - dVdr += preSw * 3.0 * 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); - - } 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; - } - - // 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; - - // 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); - } } } - 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; + if (idat.doElectricField) { + *(idat.eField1) += Ea * *(idat.electroMult); + *(idat.eField2) += Eb * *(idat.electroMult); + } - pref = idat.electroMult * pre14_ * q_j * one_third_; + if (a_is_Fluctuating) *(idat.dVdFQ1) += dUdCa * *(idat.sw); + if (b_is_Fluctuating) *(idat.dVdFQ2) += dUdCb * *(idat.sw); - 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; + 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); - // 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; + if (b_is_Dipole || b_is_Quadrupole) + *(idat.t2) += Tb * *(idat.sw); + + } else { - // calculate the derivatives for the forces and torques + // only accumulate the forces and torques resulting from the + // indirect reaction field terms. - 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)); - - 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.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; } - - 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 (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); - } - return; } + + void Electrostatic::calcSelfCorrection(SelfData &sdat) { - void Electrostatic::calcSkipCorrection(SkipCorrectionData skdat) { - if (!initialized_) initialize(); + + ElectrostaticAtomData data = ElectrostaticMap[sdat.atype]; - ElectrostaticAtomData data1 = ElectrostaticMap[skdat.atype1]; - ElectrostaticAtomData data2 = ElectrostaticMap[skdat.atype2]; - // logicals - - bool i_is_Charge = data1.is_Charge; - bool i_is_Dipole = data1.is_Dipole; - - bool j_is_Charge = data2.is_Charge; - bool j_is_Dipole = data2.is_Dipole; - - RealType q_i, q_j; + bool i_is_Charge = data.is_Charge; + bool i_is_Dipole = data.is_Dipole; + bool i_is_Fluctuating = data.is_Fluctuating; + RealType C_a = data.fixedCharge; + RealType self, preVal, DadDa; - // The skippedCharge computation is needed by the real-space cutoff methods - // (i.e. shifted force and shifted potential) - - if (i_is_Charge) { - q_i = data1.charge; - skdat.skippedCharge2 += q_i; + if (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); } - if (j_is_Charge) { - q_j = data2.charge; - skdat.skippedCharge1 += q_j; - } - - // the rest of this function should only be necessary for reaction field. - - if (summationMethod_ == REACTION_FIELD) { - RealType riji, ri2, ri3; - RealType q_i, mu_i, ct_i; - RealType q_j, mu_j, ct_j; - RealType preVal, rfVal, vterm, dudr, pref, myPot; - Vector3d dVdr, uz_i, uz_j, duduz_i, duduz_j, rhat; - - // some variables we'll need independent of electrostatic type: + 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); - } + DadDa = data.dipole.lengthSquare(); + (*(sdat.pot))[ELECTROSTATIC_FAMILY] -= pre22_ * preRF_ * DadDa; } - // accumulate the forces and torques resulting from the self term - skdat.pot += myPot; - skdat.f1 += dVdr; + break; - if (i_is_Dipole) - *(skdat.t1) -= cross(uz_i, duduz_i); - if (j_is_Dipole) - *(skdat.t2) -= cross(uz_j, duduz_j); - } - } - - void Electrostatic::calcSelfCorrection(SelfCorrectionData scdat) { - RealType mu1, preVal, chg1, self; - - if (!initialized_) initialize(); - - ElectrostaticAtomData data = ElectrostaticMap[scdat.atype]; - - // logicals - - bool i_is_Charge = data.is_Charge; - bool i_is_Dipole = data.is_Dipole; - - 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; - - // 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) { + case esm_SHIFTED_FORCE: + case esm_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_; - } - scdat.pot += self; + self = -0.5 * selfMult_ * C_a * (C_a + *(sdat.skippedCharge)) * pre11_; + (*(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; + } }