Updates to AmoebaVdwForce, AmoebaGeneralizedKirkwoodForce and AmoebaWcaDispersionForce (#4647)

* Update the AMOEBA OpenMM API for vdW, GK and WCA

* Changes needed for the Corrigan et al Generalized Kirkwood model and minor changes to the vdW force to support CpHMD

* Add casts to real for uses of POW in GK; Pass force by reference within the WCA kernel

* Update swigInputConfig for Amoeba vdW and GK forces

* Update TestAPIUnits.testAmoebaVdwForce

* Set the units for getSolventDielectric and getSoluteDielectric to None

* Update default dispersion offset parameter for the AmoebaWcaDispersionForce

* Remove overloaded getParticleParameters and setParticleParameters from AmoebaGeneralizedKirkwoodForce

* Update the AmoebaWcaDispersionForce TestAPIUnits tests to reflect using the correct units for the C++ parameter default values; Update the alanine-dipeptide-amoeba-forces to reflect the updated GK model

* Move neck descreening constants into AmoebaGeneralizedKirkwoodForceImpl; set the default GK dielecticOffset to 0.09; set the default WCA shctd parameter to 0.82

* Fix Python test cases for WCA and GK

* Load AMOEBA/GK parameters into an array of float4

* Cleaned up the AmoebaGeneralizedKirkwoodForce based on feedback from Peter; the one case where backwards compatibility remains a challenge is application of the dielectric offset parameter - in the prior code this was only applied to the nonpolar cavity term, but not to calculation of Born radii; in this revision the dielectric offset is applied to BOTH the nonpolar cavity term and to calculation of Born radii. At this point I can't think of elegant way to maintain backwards compatibility that isn't confusing, nor does it make sense (at least to me) to only apply the concept of the dieletric offset to one aspect (i.e. only to nonpolar cavity or only to Born radii calculation) but not to both.

* Remove 'using std::vector' from AmoebaGeneralizedKirkwoodForceImpl.h; divide by 10 instead of multiplying by 0.1f in amoebaGk.cc

* Added a parameter called descreenOffset, which is applied during calculation of effective Born radii for GK. The parameter dielectricOffset is only used for the nonpolar cavity term consistent with its prior use. All tests in TestAmoebaGeneralizedKirkwoodForce.h are now backwards compatible with their behavior prior to this PR.

* Change two constants in amoebaGk.cc to single precision; Improved the documentation for getNeckConstants in AmoebaGeneralizedKirkwoodForceImpl.h

* Fix comment for setTanhRescaling in AmoebaGeneralizedKirkwoodForce.h, Fix comment for setTanhParameters in AmoebaReferenceGeneralizedKirkwoodForce.h; set the type of parameter GeneralizedKirkwoodTanhRescaling to bool in AmoebaGeneralizedKirkwoodForceProxy.cpp; In ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel return references of per particle parameters instead of copies; update AmoebaReferenceKernels.h method signatures for per particle parameters to return const vector references

* Minor tweaks to the documentation for the tanh rescaling flag

* Improve the comments for the get and setTanhParameters in AmoebaGeneralizedKirkwoodForce.h and AmoebaReferenceGeneralizedKirkwoodForce.h

plugins/amoeba/openmmapi/include/openmm/AmoebaGeneralizedKirkwoodForce.h

@@ -69,6 +69,10 @@ public:
      * Add the  parameters for a particle.  This should be called once for each particle
      * in the System.  When it is called for the i'th time, it specifies the parameters for the i'th particle.
      *
+     * This method is provided for backwards compatibility. Compared to the alternative five parameter addParticle
+     * method, the descreenRadius parameter is set to base radius value and the neckFactor is set to zero
+     * (no neck descreening).
+     *
      * @param charge         the charge of the particle, measured in units of the proton charge
      * @param radius         the atomic radius of the particle, measured in nm
      * @param scalingFactor  the scaling factor for the particle
@@ -76,25 +80,55 @@ public:
      */
     int addParticle(double charge, double radius, double scalingFactor);
 
+    /**
+     * Add the  parameters for a particle.  This should be called once for each particle
+     * in the System.  When it is called for the i'th time, it specifies the parameters for the i'th particle.
+     * <p>
+     * For generalized Born / generalized Kirkwood methods, the radius of each atom has two roles. The first
+     * is to define the base radius of the atom when computing its effective radius. This base radius is usually
+     * parameterized against solvation free energy differences. The second role is to describe how much continuum
+     * water is displaced when the atom descreens water for the calculation of the Born radii of other atoms.
+     * Separation of the two roles into the "radius" and "descreenRadius" parameters gives model developers more
+     * control over these separate roles.
+     * <p>
+     * For example, the fitting of base "radius" values will usually result in deviation from the force field's
+     * van der Waals definition of Rmin (or sigma). The descreenRadius can be defined separately using force field
+     * van der Waals Rmin values, which maintains consistency of atomic sizes during the HCT pairwise
+     * descreening integral. The "scalingFactor" is applied to the descreenRadius during the HCT pairwise
+     * descreening integral, while the neckFactor (if greater than zero) includes neck contributions to descreening.
+     *
+     * @param charge         the charge of the particle, measured in units of the proton charge
+     * @param radius         the atomic radius of the particle, measured in nm
+     * @param scalingFactor  the scaling factor for the particle (unitless)
+     * @param descreenRadius the atomic radius of the particle for descreening, measure in nm
+     * @param neckFactor     the scaling factor for interstitial neck descreening (unitless)
+     * @return the index of the particle that was added
+    */
+    int addParticle(double charge, double radius, double scalingFactor, double descreenRadius, double neckFactor);
+
     /**
      * Get the force field parameters for a particle.
-     * 
+     *
      * @param      index          the index of the particle for which to get parameters
      * @param[out] charge         the charge of the particle, measured in units of the proton charge
      * @param[out] radius         the atomic radius of the particle, measured in nm
      * @param[out] scalingFactor  the scaling factor for the particle
-     */
-    void getParticleParameters(int index, double& charge, double& radius, double& scalingFactor) const;
+     * @param[out] descreenRadius the atomic radius of the particle for descreening, measure in nm
+     * @param[out] neckFactor     the scaling factor for interstitial neck descreening (unitless)
+    */
+    void getParticleParameters(int index, double& charge, double& radius, double& scalingFactor, double& descreenRadius, double& neckFactor) const;
 
     /**
      * Set the force field parameters for a particle.
-     * 
-     * @param index          the index of the particle for which to set parameters
-     * @param charge         the charge of the particle, measured in units of the proton charge
-     * @param radius         the atomic radius of the particle, measured in nm
-     * @param scalingFactor  the scaling factor for the particle
-     */
-    void setParticleParameters(int index, double charge, double radius, double scalingFactor);
+    *
+    * @param index          the index of the particle for which to set parameters
+    * @param charge         the charge of the particle, measured in units of the proton charge
+    * @param radius         the atomic radius of the particle, measured in nm
+    * @param scalingFactor  the scaling factor for the particle
+    * @param descreenRadius the atomic radius of the particle for descreening, measure in nm
+    * @param neckFactor     the scaling factor for interstitial neck descreening (unitless)
+    */
+    void setParticleParameters(int index, double charge, double radius, double scalingFactor, double descreenRadius, double neckFactor);
 
     /**
      * Get the dielectric constant for the solvent.
@@ -119,38 +153,116 @@ public:
 
     /**
      * Set the dielectric constant for the solute.
+     *
+     * @param dielectric The solute dielectric constant.
      */
     void setSoluteDielectric(double dielectric) {
         soluteDielectric = dielectric;
     }
 
+    /**
+     * Get the flag signaling whether the solute descreening integral is rescaled by a Tanh function
+     * to account for interstitial spaces. If True then tanh rescaling is used.
+    */
+    bool getTanhRescaling() const {
+        return tanhRescaling;
+    }
+
+    /**
+     * Set the flag signaling whether the solute descreening integral is rescaled by a Tanh function
+     * to account for interstitial spaces.
+     *
+     * @param tanhRescale False to turn off Tanh rescaling; true to turn on.
+    */
+    void setTanhRescaling(bool tanhRescale) {
+        this->tanhRescaling = tanhRescale;
+    }
+
+    /**
+     * Get Tanh function parameters b0, b1 and b2.
+     *
+     * @param b0 The first tanh parameter.
+     * @param b1 The second tanh parameter.
+     * @param b2 The third tanh parameter.
+     */
+    void getTanhParameters(double& b0, double& b1, double& b2) const {
+        b0 = this->beta0;
+        b1 = this->beta1;
+        b2 = this->beta2;
+    }
+
+    /**
+     * Set the the Tanh function parameters to account for interstitial spaces.
+     *
+     * @param b0 The first tanh parameter.
+     * @param b1 The second tanh parameter.
+     * @param b2 The third tanh parameter.
+     */
+    void setTanhParameters(double b0, double b1, double b2) {
+        this->beta0 = b0;
+        this->beta1 = b1;
+        this->beta2 = b2;
+    }
+
+    /**
+    * Get the offset added to the atomic radius of each atom that sets the beginning of the
+    * descreening integral when calculating effective Born radii.
+    */
+    double getDescreenOffset() const;
+
+    /**
+     * Get the offset added to the atomic radius of each atom that sets the beginning of the
+     * descreening integral when calculating effective Born radii.
+     *
+     * @param descreenOffet The descreening offset (nm).
+     */
+    void setDescreenOffset(double descreenOffet);
+
     /**
      * Get the flag signaling whether the cavity term should be included
      */
     int getIncludeCavityTerm() const;
 
     /**
-     * Set the flag signaling whether the cavity term should be included
+     * Set the flag signaling whether the cavity term should be included.
+     *
+     * @param includeCavityTerm Zero to turn off the cavity term; one to turn on.
      */
     void setIncludeCavityTerm(int includeCavityTerm);
 
     /**
-     * Get the probe radius (nm) used in SASA contribution
+     * Get the probe radius (nm) used for the cavity contribution.
      */
     double getProbeRadius() const;
 
     /**
-     * Set the probe radius (nm) used in SASA contribution
+     * Set the probe radius (nm) used for the cavity contribution.
+     *
+     * @param probeRadius The probeRadius for the cavity term.
      */
     void setProbeRadius(double probeRadius);
 
+    /**
+     * Get the dielectric offset (nm) used for cavity contribution.
+     */
+    double getDielectricOffset() const;
+
+    /**
+     * Set the dielectric offset (nm) used for cavity contribution.
+     *
+     * @param dielectricOffset The dielectric offset (nm).
+     */
+    void setDielectricOffset(double dielectricOffset);
+
     /**
      * Get the surface area factor kJ/(nm*nm) used in SASA contribution
      */
     double getSurfaceAreaFactor() const;
 
     /**
-     * Set the surface area factor kJ/(nm*nm) used in SASA contribution
+     * Set the surface area factor kJ/(nm*nm) used in SASA contribution.
+     *
+     * @param surfaceAreaFactor The surface area factor in kJ/(nm*nm).
      */
     void setSurfaceAreaFactor(double surfaceAreaFactor);
     /**
@@ -177,8 +289,10 @@ protected:
 private:
     class ParticleInfo;
     int includeCavityTerm;
+    bool tanhRescaling;
     double solventDielectric, soluteDielectric, dielectricOffset,
            probeRadius, surfaceAreaFactor;
+    double beta0, beta1, beta2, descreenOffset;
     std::vector<ParticleInfo> particles;
 };
 
@@ -188,12 +302,12 @@ private:
  */
 class AmoebaGeneralizedKirkwoodForce::ParticleInfo {
 public:
-    double charge, radius, scalingFactor;
+    double charge, radius, scalingFactor, descreenRadius, neckFactor;
     ParticleInfo() {
-        charge = radius = scalingFactor = 0.0;
+        charge = radius = scalingFactor = descreenRadius = neckFactor = 0.0;
     }
-    ParticleInfo(double charge, double radius, double scalingFactor) :
-        charge(charge), radius(radius), scalingFactor(scalingFactor) {
+    ParticleInfo(double charge, double radius, double scalingFactor, double descreenRadius, double neckFactor) :
+        charge(charge), radius(radius), scalingFactor(scalingFactor), descreenRadius(descreenRadius), neckFactor(neckFactor) {
     }
 };
 

plugins/amoeba/openmmapi/include/openmm/AmoebaVdwForce.h

@@ -173,10 +173,11 @@ public:
      * @param reductionFactor the fraction of the distance along the line from the parent particle to this particle
      *                        at which the interaction site should be placed
      * @param isAlchemical    if true, this vdW particle is undergoing an alchemical change.
-     * @param typeIndex       the index of the particle type for this particle
+     * @param typeIndex       the index of the particle type for this particle.
+     * @param scaleFactor     a scale factor to apply to all interactions involving this particle (used for CpHMD).
      */
     void setParticleParameters(int particleIndex, int parentIndex, double sigma, double epsilon, 
-                               double reductionFactor, bool isAlchemical=false, int typeIndex=-1);
+                               double reductionFactor, bool isAlchemical=false, int typeIndex=-1, double scaleFactor = 1.0);
 
     /**
      * Get the force field parameters for a vdw particle.
@@ -189,9 +190,10 @@ public:
      *                             at which the interaction site should be placed
      * @param[out] isAlchemical    if true, this vdW particle is undergoing an alchemical change.
      * @param[out] typeIndex       the index of the particle type for this particle
+     * @param[out] scaleFactor     a scale factor to apply to all interactions involving this particle (used for CpHMD).
      */
     void getParticleParameters(int particleIndex, int& parentIndex, double& sigma, double& epsilon, 
-                               double& reductionFactor, bool& isAlchemical, int& typeIndex) const;
+                               double& reductionFactor, bool& isAlchemical, int& typeIndex, double& scaleFactor) const;
 
     /**
      * Add the force field parameters for a vdw particle.  This version is used when parameters
@@ -203,9 +205,10 @@ public:
      * @param reductionFactor the fraction of the distance along the line from the parent particle to this particle
      *                        at which the interaction site should be placed
      * @param isAlchemical    if true, this vdW particle is undergoing an alchemical change.
+     * @param scaleFactor     a scale factor to apply to all interactions involving this particle (used for CpHMD).
      * @return index of added particle
      */
-    int addParticle(int parentIndex, double sigma, double epsilon, double reductionFactor, bool isAlchemical = false);
+    int addParticle(int parentIndex, double sigma, double epsilon, double reductionFactor, bool isAlchemical = false, double scaleFactor = 1.0);
 
     /**
      * Add the force field parameters for a vdw particle. This version is used when parameters
@@ -216,9 +219,10 @@ public:
      * @param reductionFactor the fraction of the distance along the line from the parent particle to this particle
      *                        at which the interaction site should be placed
      * @param isAlchemical    if true, this vdW particle is undergoing an alchemical change.
+     * @param scaleFactor     a scale factor to apply to all interactions involving this particle (used for CpHMD).
      * @return index of added particle
      */
-    int addParticle(int parentIndex, int typeIndex, double reductionFactor, bool isAlchemical = false);
+    int addParticle(int parentIndex, int typeIndex, double reductionFactor, bool isAlchemical = false, double scaleFactor = 1.0);
 
     /**
      * Add a particle type.
@@ -489,17 +493,20 @@ private:
 class AmoebaVdwForce::VdwInfo {
 public:
     int parentIndex, typeIndex;
-    double reductionFactor, sigma, epsilon, cutoff;
+    double reductionFactor, sigma, epsilon, scaleFactor;
     bool isAlchemical;
     VdwInfo() {
         parentIndex = -1;
+        typeIndex   = -1;
         reductionFactor      = 0.0;
         sigma                = 1.0;
         epsilon              = 0.0;
+        scaleFactor          = 1.0;
         isAlchemical         = false;
     }
-    VdwInfo(int parentIndex, double sigma, double epsilon, int typeIndex, double reductionFactor, bool isAlchemical) :
-        parentIndex(parentIndex), reductionFactor(reductionFactor), sigma(sigma), epsilon(epsilon), typeIndex(typeIndex), isAlchemical(isAlchemical) {
+    VdwInfo(int parentIndex, double sigma, double epsilon, int typeIndex, double reductionFactor, bool isAlchemical, double scaleFactor) :
+        parentIndex(parentIndex), reductionFactor(reductionFactor), sigma(sigma), epsilon(epsilon),
+        typeIndex(typeIndex), isAlchemical(isAlchemical), scaleFactor(scaleFactor) {
     }
 };
 

plugins/amoeba/openmmapi/include/openmm/internal/AmoebaGeneralizedKirkwoodForceImpl.h

@@ -62,6 +62,54 @@ public:
         return kernel;
     }
     void updateParametersInContext(ContextImpl& context);
+
+    /**
+     * Compute a "neck" descreening contribution.
+     * <p>
+     * Using the tabulated Aij and Bij values, compute the neck integral
+     * using Equations 13 and 14 from Aguilar, Shadrach and Onufriev (10.1021/ct100392h).
+     *
+     * @param r separation distance.
+     * @param radius Radius of the current atom.
+     * @param radiusK Radius of the descreening atom.
+     * @param sneck The neck scale factor.
+     * @return The integral value.
+     */
+    static double neckDescreen(double r, double radius, double radiusK, double sneck);
+
+    /**
+     * Get Neck Aij and Bij constants.
+     *
+     * Based on the radii of two input atoms ("radius" for the atom being descreened
+     * and "radiusK" for the descreening atom), the neck parameters Aij and Bij are
+     * interpolated from tabulated values. The definitions of Aij and Bij are defined by
+     * Eq. 11 from Corrigan et al. (10.1063/5.0158914) while Figure 2 describes
+     * a neck descreening region.
+     *
+     * @param radius Radius of the current atom.
+     * @param radiusK Radius of the descreening atom.
+     * @param aij The Aij neck constant.
+     * @param bij the Bij neck constant.
+     */
+    static void getNeckConstants(double radius, double radiusK, double &aij, double &bij);
+
+     /**
+      * The array of neck tabulated radii values.
+      */
+    static const std::vector<float>& getNeckRadii();
+
+    const static int NUM_NECK_RADII = 45;
+
+    /**
+     * The tabulated Aij parameters.
+     */
+    const static float (&getAij())[NUM_NECK_RADII][NUM_NECK_RADII];
+
+     /**
+      * The tabulated Bij parameters.
+      */
+     const static float (&getBij())[NUM_NECK_RADII][NUM_NECK_RADII];
+
 private:
     const AmoebaGeneralizedKirkwoodForce& owner;
     Kernel kernel;

plugins/amoeba/openmmapi/src/AmoebaGeneralizedKirkwoodForce.cpp

@@ -36,35 +36,45 @@
 
 using namespace OpenMM;
 
-AmoebaGeneralizedKirkwoodForce::AmoebaGeneralizedKirkwoodForce() : solventDielectric(78.3), soluteDielectric(1.0), dielectricOffset(0.009), includeCavityTerm(1), probeRadius(0.14) {
-
+AmoebaGeneralizedKirkwoodForce::AmoebaGeneralizedKirkwoodForce() :
+    solventDielectric(78.3), soluteDielectric(1.0), dielectricOffset(0.009), includeCavityTerm(1), probeRadius(0.14),
+    tanhRescaling(false), beta0(0.9563), beta1(0.2578), beta2(0.0810), descreenOffset(0.0) {
      surfaceAreaFactor = -6.0* 3.1415926535*0.0216*1000.0*0.4184;
 }
 
 int AmoebaGeneralizedKirkwoodForce::addParticle(double charge, double radius, double scalingFactor) {
-    particles.push_back(ParticleInfo(charge, radius, scalingFactor));
+    particles.push_back(ParticleInfo(charge, radius, scalingFactor, radius, 0.0));
     return particles.size()-1;
 }
 
-void AmoebaGeneralizedKirkwoodForce::getParticleParameters(int index, double& charge, double& radius, double& scalingFactor) const {
+int AmoebaGeneralizedKirkwoodForce::addParticle(double charge, double radius, double scalingFactor, double descreenRadius, double neckFactor) {
+    particles.push_back(ParticleInfo(charge, radius, scalingFactor, descreenRadius, neckFactor));
+    return particles.size()-1;
+}
+
+void AmoebaGeneralizedKirkwoodForce::getParticleParameters(int index, double& charge, double& radius, double& scalingFactor, double& descreenRadius, double& neckFactor) const {
     charge = particles[index].charge;
     radius = particles[index].radius;
     scalingFactor = particles[index].scalingFactor;
+    descreenRadius = particles[index].descreenRadius;
+    neckFactor = particles[index].neckFactor;
 }
 
-void AmoebaGeneralizedKirkwoodForce::setParticleParameters(int index, double charge, double radius, double scalingFactor) {
+void AmoebaGeneralizedKirkwoodForce::setParticleParameters(int index, double charge, double radius, double scalingFactor, double descreenRadius, double neckFactor) {
     particles[index].charge = charge;
     particles[index].radius = radius;
     particles[index].scalingFactor = scalingFactor;
+    particles[index].descreenRadius = descreenRadius;
+    particles[index].neckFactor = neckFactor;
 }
-/*
+
 double AmoebaGeneralizedKirkwoodForce::getDielectricOffset() const {
     return dielectricOffset;
 }
 
 void AmoebaGeneralizedKirkwoodForce::setDielectricOffset(double inputDielectricOffset) {
     dielectricOffset = inputDielectricOffset;
-} */
+}
 
 int AmoebaGeneralizedKirkwoodForce::getIncludeCavityTerm() const {
     return includeCavityTerm;
@@ -90,6 +100,14 @@ void AmoebaGeneralizedKirkwoodForce::setSurfaceAreaFactor(double inputSurfaceAre
     surfaceAreaFactor = inputSurfaceAreaFactor;
 }
 
+double AmoebaGeneralizedKirkwoodForce::getDescreenOffset() const {
+    return descreenOffset;
+}
+
+void AmoebaGeneralizedKirkwoodForce::setDescreenOffset(double inputDescreenOffet) {
+    descreenOffset = inputDescreenOffet;
+}
+
 ForceImpl* AmoebaGeneralizedKirkwoodForce::createImpl() const {
     return new AmoebaGeneralizedKirkwoodForceImpl(*this);
 }

plugins/amoeba/openmmapi/src/AmoebaVdwForce.cpp

@@ -44,23 +44,23 @@ AmoebaVdwForce::AmoebaVdwForce() : nonbondedMethod(NoCutoff), potentialFunction(
         useTypes(false), alchemicalMethod(None), n(5), alpha(0.7) {
 }
 
-int AmoebaVdwForce::addParticle(int parentIndex, double sigma, double epsilon, double reductionFactor, bool isAlchemical) {
+int AmoebaVdwForce::addParticle(int parentIndex, double sigma, double epsilon, double reductionFactor, bool isAlchemical, double scaleFactor) {
     if (useTypes)
         throw OpenMMException("AmoebaVdwForce: must use the same version of addParticle() for all particles");
-    parameters.push_back(VdwInfo(parentIndex, sigma, epsilon, -1, reductionFactor, isAlchemical));
+    parameters.push_back(VdwInfo(parentIndex, sigma, epsilon, -1, reductionFactor, isAlchemical, scaleFactor));
     return parameters.size()-1;
 }
 
-int AmoebaVdwForce::addParticle(int parentIndex, int typeIndex, double reductionFactor, bool isAlchemical) {
+int AmoebaVdwForce::addParticle(int parentIndex, int typeIndex, double reductionFactor, bool isAlchemical, double scaleFactor) {
     if (parameters.size() > 0 && !useTypes)
         throw OpenMMException("AmoebaVdwForce: must use the same version of addParticle() for all particles");
     useTypes = true;
-    parameters.push_back(VdwInfo(parentIndex, 1.0, 0.0, typeIndex, reductionFactor, isAlchemical));
+    parameters.push_back(VdwInfo(parentIndex, 1.0, 0.0, typeIndex, reductionFactor, isAlchemical, scaleFactor));
     return parameters.size()-1;
 }
 
-void AmoebaVdwForce::getParticleParameters(int particleIndex, int& parentIndex,
-                                           double& sigma, double& epsilon, double& reductionFactor, bool& isAlchemical, int& typeIndex) const {
+void AmoebaVdwForce::getParticleParameters(int particleIndex, int& parentIndex,double& sigma, double& epsilon, double& reductionFactor,
+                                           bool& isAlchemical, int& typeIndex, double& scaleFactor) const {
     ASSERT_VALID_INDEX(particleIndex, parameters);
     parentIndex     = parameters[particleIndex].parentIndex;
     sigma           = parameters[particleIndex].sigma;
@@ -68,10 +68,11 @@ void AmoebaVdwForce::getParticleParameters(int particleIndex, int& parentIndex,
     reductionFactor = parameters[particleIndex].reductionFactor;
     isAlchemical    = parameters[particleIndex].isAlchemical;
     typeIndex       = parameters[particleIndex].typeIndex;
+    scaleFactor     = parameters[particleIndex].scaleFactor;
 }
 
-void AmoebaVdwForce::setParticleParameters(int particleIndex, int parentIndex,
-                                           double sigma, double epsilon, double reductionFactor, bool isAlchemical, int typeIndex) {
+void AmoebaVdwForce::setParticleParameters(int particleIndex, int parentIndex,double sigma, double epsilon, double reductionFactor,
+                                           bool isAlchemical, int typeIndex, double scaleFactor) {
     ASSERT_VALID_INDEX(particleIndex, parameters);
     parameters[particleIndex].parentIndex     = parentIndex;
     parameters[particleIndex].sigma           = sigma;
@@ -79,6 +80,7 @@ void AmoebaVdwForce::setParticleParameters(int particleIndex, int parentIndex,
     parameters[particleIndex].reductionFactor = reductionFactor;
     parameters[particleIndex].isAlchemical    = isAlchemical;
     parameters[particleIndex].typeIndex       = typeIndex;
+    parameters[particleIndex].scaleFactor     = scaleFactor;
 }
 
 int AmoebaVdwForce::addParticleType(double sigma, double epsilon) {

plugins/amoeba/openmmapi/src/AmoebaVdwForceImpl.cpp

@@ -55,9 +55,9 @@ void AmoebaVdwForceImpl::initialize(ContextImpl& context) {
         throw OpenMMException("AmoebaVdwForce must have exactly as many particles as the System it belongs to.");
     for (int i = 0; i < owner.getNumParticles(); i++) {
         int parentIndex, typeIndex;
-        double sigma, epsilon, reductionFactor;
+        double sigma, epsilon, reductionFactor, scaleFactor;
         bool isAlchemical;
-        owner.getParticleParameters(i, parentIndex, sigma, epsilon, reductionFactor, isAlchemical, typeIndex);
+        owner.getParticleParameters(i, parentIndex, sigma, epsilon, reductionFactor, isAlchemical, typeIndex, scaleFactor);
         if (sigma < 0)
             throw OpenMMException("AmoebaVdwForce: sigma for a particle cannot be negative");
         if (owner.getPotentialFunction() == AmoebaVdwForce::Buffered147 && sigma == 0)
@@ -101,11 +101,11 @@ void AmoebaVdwForceImpl::createParameterMatrix(const AmoebaVdwForce& force, vect
     if (force.getUseParticleTypes()) {
         // We get the types directly from the particles.
 
-        double sigma, epsilon, reduction;
+        double sigma, epsilon, reduction, scaleFactor;
         int parent;
         bool isAlchemical;
         for (int i = 0; i < numParticles; i++)
-            force.getParticleParameters(i, parent, sigma, epsilon, reduction, isAlchemical, type[i]);
+            force.getParticleParameters(i, parent, sigma, epsilon, reduction, isAlchemical, type[i], scaleFactor);
         numTypes = force.getNumParticleTypes();
         typeSigma.resize(numTypes);
         typeEpsilon.resize(numTypes);
@@ -117,10 +117,10 @@ void AmoebaVdwForceImpl::createParameterMatrix(const AmoebaVdwForce& force, vect
 
         map<pair<double, double>, int> typeForParams;
         for (int i = 0; i < numParticles; i++) {
-            double sigma, epsilon, reduction;
+            double sigma, epsilon, reduction, scaleFactor;
             int parent, typeIndex;
             bool isAlchemical;
-            force.getParticleParameters(i, parent, sigma, epsilon, reduction, isAlchemical, typeIndex);
+            force.getParticleParameters(i, parent, sigma, epsilon, reduction, isAlchemical, typeIndex, scaleFactor);
             pair<double, double> params = make_pair(sigma, epsilon);
             map<pair<double, double>, int>::iterator entry = typeForParams.find(params);
             if (entry == typeForParams.end()) {

plugins/amoeba/openmmapi/src/AmoebaWcaDispersionForce.cpp

@@ -38,14 +38,22 @@
 using namespace OpenMM;
 
 AmoebaWcaDispersionForce::AmoebaWcaDispersionForce() {
-    epso      = 0.1100;
-    epsh      = 0.0135;
-    rmino     = 1.7025;
-    rminh     = 1.3275;
-    awater    = 0.033428;
-    slevy     = 1.0;
-    shctd     = 0.81;
-    dispoff   = 0.26;
+    // Amoeba Water '03 vdW parameters (Diameters in Angstroms; Well depth in kcal/mole)
+    // vdw           1               3.4050     0.1100
+    // vdw           2               2.6550     0.0135      0.910
+
+    // Convert kcal/mol to kJ/mol
+    epso = 0.1100 * 4.184e0;
+    epsh = 0.0135 * 4.184e0;
+    // Convert A to nm.
+    rmino = 1.7025e-01;
+    rminh = 1.3275e-01;
+    dispoff = 1.056e-01;
+    // Convert water number density from water / A^3 to water / nm^3.
+    awater = 0.033428e03;
+    // No units.
+    slevy = 1.0;
+    shctd = 0.82;
 }
 
 int AmoebaWcaDispersionForce::addParticle(double radius, double epsilon) {

plugins/amoeba/openmmapi/src/AmoebaWcaDispersionForceImpl.cpp

@@ -61,66 +61,66 @@ double AmoebaWcaDispersionForceImpl::calcForcesAndEnergy(ContextImpl& context, b
         return kernel.getAs<CalcAmoebaWcaDispersionForceKernel>().execute(context, includeForces, includeEnergy);
     return 0.0;
 }
-void AmoebaWcaDispersionForceImpl::getMaximumDispersionEnergy(const AmoebaWcaDispersionForce& force, int particleIndex, double& maxDispersionEnergy) {
 
-    const double pi  = 3.1415926535897;
+static double tailCorrection(double ri, double eps, double rmin) {
+    if (ri < rmin) {
+        double r3 = ri * ri * ri;
+        double rmin3 = rmin * rmin * rmin;
+        return -eps * 4.0e0 * M_PI * (rmin3 - r3) / 3.0e0 - eps * M_PI * 18.0e0 * rmin3 / 11.0e0;
+    } else {
+        double ri2 = ri * ri;
+        double ri4 = ri2 * ri2;
+        double ri7 = ri * ri2 * ri4;
+        double ri11 = ri7 * ri4;
+        double rmin2 = rmin * rmin;
+        double rmin4 = rmin2 * rmin2;
+        double rmin7 = rmin * rmin2 * rmin4;
+        return 2.0e0 * M_PI * eps * rmin7 * (2.0e0 * rmin7 - 11.0e0 * ri7) / (11.0e0 * ri11);
+    }
+}
 
-    // from last loop in subroutine knp in ksolv.f
+void AmoebaWcaDispersionForceImpl::getMaximumDispersionEnergy(const AmoebaWcaDispersionForce &force, int particleIndex,
+                                                              double &maxDispersionEnergy) {
 
-    double rdisp, epsi;
-    force.getParticleParameters(particleIndex, rdisp, epsi);
-    if (epsi <= 0.0 || rdisp <= 0.0) {
+    // Atom i parameters
+    double rmini, epsi;
+    force.getParticleParameters(particleIndex, rmini, epsi);
+    if (epsi <= 0.0 || rmini <= 0.0) {
         maxDispersionEnergy = 0.0;
         return;
     }
-    double rmini     = rdisp;
-    rdisp           += force.getDispoff();
-    double epso      = force.getEpso();
-    double emixo     = std::sqrt(epso) + std::sqrt(epsi);
-           emixo     = 4.0*epso*epsi/(emixo*emixo);
-
-    double rmino     = force.getRmino();
-    double rmino2    = rmino*rmino;
-    double rmini2    = rmini*rmini;
-    double rmixo     = 2.0*(rmino2*rmino + rmini2*rmini) / (rmino2 + rmini2);
-
-    double rmixo3    = rmixo*rmixo*rmixo;
-    double rmixo7    = rmixo*rmixo3*rmixo3;
-    double ao        = emixo*rmixo7;
-
-    double epsh      = force.getEpsh();
-    double emixh     = std::sqrt(epsh) + std::sqrt(epsi);
-           emixh     = 4.0*epsh*epsi/(emixh*emixh);
-    double rminh     = force.getRminh();
-    double rminh2    = rminh*rminh;
-    double rmixh     = rminh*rminh + rmini2;
-           rmixh     = 2.0 * (rminh2*rminh + rmini2*rmini) / (rminh2 + rmini2);
-    double rmixh3    = rmixh*rmixh*rmixh;
-    double rmixh7    = rmixh3*rmixh3*rmixh;
-    double ah        = emixh*rmixh7;
-
-    double rdisp3    = rdisp*rdisp*rdisp;
-    double rdisp7    = rdisp*rdisp3*rdisp3;
-    double rdisp11   = rdisp7*rdisp3*rdisp;
-
-    double cdisp;
-    if (rdisp < rmixh) {
-        cdisp = -4.0*pi*emixh*(rmixh3-rdisp3)/3.0 - emixh*18.0/11.0*rmixh3*pi;
-    } else {
-        cdisp = 2.0*pi*(2.0*rmixh7-11.0*rdisp7)*ah/ (11.0*rdisp11);
-    }
-    cdisp *= 2.0;
-    if (rdisp < rmixo) {
-        cdisp -= 4.0*pi*emixo*(rmixo3-rdisp3)/3.0;
-        cdisp -= emixo*18.0/11.0*rmixo3*pi;
-    } else {
-        cdisp += 2.0*pi*(2.0*rmixo7-11.0*rdisp7) * ao/(11.0*rdisp11);
-    }
-    maxDispersionEnergy = force.getSlevy()*force.getAwater()*cdisp;
+    double rmini2 = rmini * rmini;
+    double sepsi = std::sqrt(epsi);
+
+    // Atom i with water oxygen.
+    double epso = force.getEpso();
+    double emixo = std::sqrt(epso) + sepsi;
+    emixo = 4.0 * epso * epsi / (emixo * emixo);
+
+    double rmino = force.getRmino();
+    double rmino2 = rmino * rmino;
+    double rmixo = 2.0 * (rmino2 * rmino + rmini2 * rmini) / (rmino2 + rmini2);
+
+    // Atom i with water hydrogen.
+    double epsh = force.getEpsh();
+    double emixh = std::sqrt(epsh) + sepsi;
+    emixh = 4.0 * epsh * epsi / (emixh * emixh);
+
+    double rminh = force.getRminh();
+    double rminh2 = rminh * rminh;
+    double rmixh = 2.0 * (rminh2 * rminh + rmini2 * rmini) / (rminh2 + rmini2);
+
+    // Compute the tail correction (i.e. the dispersion integral in the absence of any other solute atoms).
+    double dispersionOffest = force.getDispoff();
+    double riO = rmixo / 2.0e0 + dispersionOffest;
+    double cdisp = tailCorrection(riO, emixo, rmixo);
+    double riH = rmixh / 2.0e0 + dispersionOffest;
+    cdisp += 2.0e0 * tailCorrection(riH, emixh, rmixh);
+
+    maxDispersionEnergy = force.getSlevy() * force.getAwater() * cdisp;
 }
 
-double AmoebaWcaDispersionForceImpl::getTotalMaximumDispersionEnergy(const AmoebaWcaDispersionForce& force) {
-
+double AmoebaWcaDispersionForceImpl::getTotalMaximumDispersionEnergy(const AmoebaWcaDispersionForce &force) {
     double totalMaximumDispersionEnergy = 0.0;
     for (int ii = 0; ii < force.getNumParticles(); ii++) {
         double maximumDispersionEnergy;

plugins/amoeba/platforms/common/src/AmoebaCommonKernels.cpp

@@ -1722,10 +1722,10 @@ public:
     ForceInfo(const AmoebaGeneralizedKirkwoodForce& force) : force(force) {
     }
     bool areParticlesIdentical(int particle1, int particle2) {
-        double charge1, charge2, radius1, radius2, scale1, scale2;
-        force.getParticleParameters(particle1, charge1, radius1, scale1);
-        force.getParticleParameters(particle2, charge2, radius2, scale2);
-        return (charge1 == charge2 && radius1 == radius2 && scale1 == scale2);
+        double charge1, charge2, radius1, radius2, scale1, scale2, descreen1, descreen2, neck1, neck2;
+        force.getParticleParameters(particle1, charge1, radius1, scale1, descreen1, neck1);
+        force.getParticleParameters(particle2, charge2, radius2, scale2, descreen2, neck2);
+        return (charge1 == charge2 && radius1 == radius2 && scale1 == scale2 && descreen1 == descreen2 && neck1 == neck2);
     }
 private:
     const AmoebaGeneralizedKirkwoodForce& force;
@@ -1747,7 +1747,16 @@ void CommonCalcAmoebaGeneralizedKirkwoodForceKernel::initialize(const System& sy
     NonbondedUtilities& nb = cc.getNonbondedUtilities();
     int paddedNumAtoms = cc.getPaddedNumAtoms();
     int elementSize = (cc.getUseDoublePrecision() ? sizeof(double) : sizeof(float));
-    params.initialize<mm_float2>(cc, paddedNumAtoms, "amoebaGkParams");
+    params.initialize<mm_float4>(cc, paddedNumAtoms,"amoebaGkParams");
+    vector<float> neckRadiiVector = AmoebaGeneralizedKirkwoodForceImpl::getNeckRadii();
+    unsigned long numNeckRadii = neckRadiiVector.size();
+    unsigned long numNeckRadii2 = numNeckRadii * numNeckRadii;
+    neckRadii.initialize<float>(cc, numNeckRadii, "neckRadii");
+    neckRadii.upload(neckRadiiVector);
+    neckA.initialize<float>( cc, numNeckRadii2, "neckA");
+    neckA.upload(AmoebaGeneralizedKirkwoodForceImpl::getAij());
+    neckB.initialize<float>( cc, numNeckRadii2, "neckB");
+    neckB.upload(AmoebaGeneralizedKirkwoodForceImpl::getBij());
     bornRadii.initialize(cc, paddedNumAtoms, elementSize, "bornRadii");
     field.initialize(cc, 3*paddedNumAtoms, sizeof(long long), "gkField");
     bornSum.initialize<long long>(cc, paddedNumAtoms, "bornSum");
@@ -1762,14 +1771,12 @@ void CommonCalcAmoebaGeneralizedKirkwoodForceKernel::initialize(const System& sy
     cc.addAutoclearBuffer(field);
     cc.addAutoclearBuffer(bornSum);
     cc.addAutoclearBuffer(bornForce);
-    vector<mm_float2> paramsVector(paddedNumAtoms);
+    vector<mm_float4> paramsVector(paddedNumAtoms);
     for (int i = 0; i < force.getNumParticles(); i++) {
-        double charge, radius, scalingFactor;
-        force.getParticleParameters(i, charge, radius, scalingFactor);
-        paramsVector[i] = mm_float2((float) radius, (float) (scalingFactor*radius));
-        
+        double charge, radius, scalingFactor, descreenRadius, neckFactor;
+        force.getParticleParameters(i, charge, radius, scalingFactor, descreenRadius, neckFactor);
+        paramsVector[i] = mm_float4((float) radius,(float) scalingFactor, (float) descreenRadius, (float) neckFactor);
         // Make sure the charge matches the one specified by the AmoebaMultipoleForce.
-        
         double charge2, thole, damping, polarity;
         int axisType, atomX, atomY, atomZ;
         vector<double> dipole, quadrupole;
@@ -1814,11 +1821,28 @@ void CommonCalcAmoebaGeneralizedKirkwoodForceKernel::initialize(const System& sy
         defines["MUTUAL_POLARIZATION"] = "";
     else if (polarizationType == AmoebaMultipoleForce::Extrapolated)
         defines["EXTRAPOLATED_POLARIZATION"] = "";
+    defines["DESCREEN_OFFSET"] = cc.doubleToString(force.getDescreenOffset());
+    tanhRescaling = force.getTanhRescaling();
+    // Max Born radius is 3 nm (30 A).
+    defines["BIG_RADIUS"] = cc.doubleToString(3.0);
+    if (tanhRescaling) {
+        defines["TANH_RESCALING"] = "true";
+        double beta0, beta1, beta2;
+        force.getTanhParameters(beta0, beta1, beta2);
+        defines["BETA0"] = cc.doubleToString(beta0);
+        defines["BETA1"] = cc.doubleToString(beta1);
+        defines["BETA2"] = cc.doubleToString(beta2);
+    } else {
+        defines["TANH_RESCALING"] = "false";
+        defines["BETA0"] = cc.doubleToString(0.0);
+        defines["BETA1"] = cc.doubleToString(0.0);
+        defines["BETA2"] = cc.doubleToString(0.0);
+    }
     includeSurfaceArea = force.getIncludeCavityTerm();
     if (includeSurfaceArea) {
         defines["SURFACE_AREA_FACTOR"] = cc.doubleToString(force.getSurfaceAreaFactor());
         defines["PROBE_RADIUS"] = cc.doubleToString(force.getProbeRadius());
-        defines["DIELECTRIC_OFFSET"] = cc.doubleToString(0.009);
+        defines["DIELECTRIC_OFFSET"] = cc.doubleToString(force.getDielectricOffset());
     }
     cc.addForce(new ForceInfo(force));
 }
@@ -1877,6 +1901,9 @@ void CommonCalcAmoebaGeneralizedKirkwoodForceKernel::computeBornRadii(ComputeArr
         computeBornSumKernel->addArg(cc.getPosq());
         computeBornSumKernel->addArg(params);
         computeBornSumKernel->addArg();
+        computeBornSumKernel->addArg(neckRadii);
+        computeBornSumKernel->addArg(neckA);
+        computeBornSumKernel->addArg(neckB);
         reduceBornSumKernel = program->createKernel("reduceBornSum");
         reduceBornSumKernel->addArg(bornSum);
         reduceBornSumKernel->addArg(params);
@@ -1900,6 +1927,10 @@ void CommonCalcAmoebaGeneralizedKirkwoodForceKernel::computeBornRadii(ComputeArr
         chainRuleKernel->addArg();
         chainRuleKernel->addArg();
         chainRuleKernel->addArg(params);
+        chainRuleKernel->addArg(neckRadii);
+        chainRuleKernel->addArg(neckA);
+        chainRuleKernel->addArg(neckB);
+        chainRuleKernel->addArg(bornSum);
         chainRuleKernel->addArg(bornRadii);
         chainRuleKernel->addArg(bornForce);
         ediffKernel = program->createKernel("computeEDiffForce");
@@ -1968,12 +1999,11 @@ void CommonCalcAmoebaGeneralizedKirkwoodForceKernel::copyParametersToContext(Con
         throw OpenMMException("updateParametersInContext: The number of particles has changed");
     
     // Record the per-particle parameters.
-    
-    vector<mm_float2> paramsVector(cc.getPaddedNumAtoms());
+    vector<mm_float4> paramsVector(cc.getPaddedNumAtoms());
     for (int i = 0; i < force.getNumParticles(); i++) {
-        double charge, radius, scalingFactor;
-        force.getParticleParameters(i, charge, radius, scalingFactor);
-        paramsVector[i] = mm_float2((float) radius, (float) (scalingFactor*radius));
+        double charge, radius, scalingFactor, descreenRadius, neckFactor;
+        force.getParticleParameters(i, charge, radius, scalingFactor, descreenRadius, neckFactor);
+        paramsVector[i] = mm_float4((float) radius, (float) scalingFactor, (float) descreenRadius, (float) neckFactor);
     }
     params.upload(paramsVector);
     cc.invalidateMolecules();
@@ -1989,11 +2019,11 @@ public:
     }
     bool areParticlesIdentical(int particle1, int particle2) {
         int iv1, iv2, type1, type2;
-        double sigma1, sigma2, epsilon1, epsilon2, reduction1, reduction2;
+        double sigma1, sigma2, epsilon1, epsilon2, reduction1, reduction2, scaleFactor1, scaleFactor2;
         bool isAlchemical1, isAlchemical2;
-        force.getParticleParameters(particle1, iv1, sigma1, epsilon1, reduction1, isAlchemical1, type1);
-        force.getParticleParameters(particle2, iv2, sigma2, epsilon2, reduction2, isAlchemical2, type2);
-        return (sigma1 == sigma2 && epsilon1 == epsilon2 && reduction1 == reduction2 && isAlchemical1 == isAlchemical2 && type1 == type2);
+        force.getParticleParameters(particle1, iv1, sigma1, epsilon1, reduction1, isAlchemical1, type1, scaleFactor1);
+        force.getParticleParameters(particle2, iv2, sigma2, epsilon2, reduction2, isAlchemical2, type2, scaleFactor2);
+        return (sigma1 == sigma2 && epsilon1 == epsilon2 && reduction1 == reduction2 && isAlchemical1 == isAlchemical2 && type1 == type2 && scaleFactor1 == scaleFactor2);
     }
 private:
     const AmoebaVdwForce& force;
@@ -2014,6 +2044,7 @@ void CommonCalcAmoebaVdwForceKernel::initialize(const System& system, const Amoe
     int paddedNumAtoms = cc.getPaddedNumAtoms();
     bondReductionAtoms.initialize<int>(cc, paddedNumAtoms, "bondReductionAtoms");
     bondReductionFactors.initialize<float>(cc, paddedNumAtoms, "bondReductionFactors");
+    scaleFactors.initialize<float>(cc, paddedNumAtoms, "scaleFactors");
     tempPosq.initialize(cc, paddedNumAtoms, cc.getUseDoublePrecision() ? sizeof(mm_double4) : sizeof(mm_float4), "tempPosq");
     tempForces.initialize<long long>(cc, 3*paddedNumAtoms, "tempForces");
     
@@ -2035,6 +2066,7 @@ void CommonCalcAmoebaVdwForceKernel::initialize(const System& system, const Amoe
     vector<float> isAlchemicalVec(paddedNumAtoms, 0);
     vector<int> bondReductionAtomsVec(paddedNumAtoms, 0);
     vector<float> bondReductionFactorsVec(paddedNumAtoms, 0);
+    vector<float> scaleFactorsVec(paddedNumAtoms, 0);
     vector<vector<int> > exclusions(cc.getNumAtoms());
 
     // Handle Alchemical parameters.
@@ -2046,17 +2078,19 @@ void CommonCalcAmoebaVdwForceKernel::initialize(const System& system, const Amoe
 
     for (int i = 0; i < force.getNumParticles(); i++) {
         int ivIndex, type;
-        double sigma, epsilon, reductionFactor;
+        double sigma, epsilon, reductionFactor, scaleFactor;
         bool alchemical;
-        force.getParticleParameters(i, ivIndex, sigma, epsilon, reductionFactor, alchemical, type);
+        force.getParticleParameters(i, ivIndex, sigma, epsilon, reductionFactor, alchemical, type, scaleFactor);
         isAlchemicalVec[i] = (alchemical) ? 1.0f : 0.0f;
         bondReductionAtomsVec[i] = ivIndex;
         bondReductionFactorsVec[i] = (float) reductionFactor;
+        scaleFactorsVec[i] = (float) scaleFactor;
         force.getParticleExclusions(i, exclusions[i]);
         exclusions[i].push_back(i);
     }
     bondReductionAtoms.upload(bondReductionAtomsVec);
     bondReductionFactors.upload(bondReductionFactorsVec);
+    scaleFactors.upload(scaleFactorsVec);
     if (force.getUseDispersionCorrection())
         dispersionCoefficient = AmoebaVdwForceImpl::calcDispersionCorrection(system, force);
     else
@@ -2069,7 +2103,6 @@ void CommonCalcAmoebaVdwForceKernel::initialize(const System& system, const Amoe
     nonbonded = cc.createNonbondedUtilities();
     nonbonded->addParameter(ComputeParameterInfo(atomType, "atomType", "int", 1));
     nonbonded->addArgument(ComputeParameterInfo(sigmaEpsilon, "sigmaEpsilon", "float", 2));
-
     if (hasAlchemical) {
        isAlchemical.upload(isAlchemicalVec);
        currentVdwLambda = 1.0f;
@@ -2077,7 +2110,7 @@ void CommonCalcAmoebaVdwForceKernel::initialize(const System& system, const Amoe
        nonbonded->addParameter(ComputeParameterInfo(isAlchemical, "isAlchemical", "float", 1));
        nonbonded->addArgument(ComputeParameterInfo(vdwLambda, "vdwLambda", "float", 1));
     }
-    
+    nonbonded->addParameter(ComputeParameterInfo(scaleFactors, "scaleFactor", "float", 1));
     // Create the interaction kernel.
     
     map<string, string> replacements;
@@ -2170,18 +2203,21 @@ void CommonCalcAmoebaVdwForceKernel::copyParametersToContext(ContextImpl& contex
     vector<float> isAlchemicalVec(cc.getPaddedNumAtoms(), 0);
     vector<int> bondReductionAtomsVec(cc.getPaddedNumAtoms(), 0);
     vector<float> bondReductionFactorsVec(cc.getPaddedNumAtoms(), 0);
+    vector<float> scaleFactorsVec(cc.getPaddedNumAtoms(), 0);
     for (int i = 0; i < force.getNumParticles(); i++) {
         int ivIndex, type;
-        double sigma, epsilon, reductionFactor;
+        double sigma, epsilon, reductionFactor, scaleFactor;
         bool alchemical;
-        force.getParticleParameters(i, ivIndex, sigma, epsilon, reductionFactor, alchemical, type);
+        force.getParticleParameters(i, ivIndex, sigma, epsilon, reductionFactor, alchemical, type, scaleFactor);
         isAlchemicalVec[i] = (alchemical) ? 1.0f : 0.0f;
         bondReductionAtomsVec[i] = ivIndex;
         bondReductionFactorsVec[i] = (float) reductionFactor;
+        scaleFactorsVec[i] = (float) scaleFactor;
     }
     if (hasAlchemical) isAlchemical.upload(isAlchemicalVec);
     bondReductionAtoms.upload(bondReductionAtomsVec);
     bondReductionFactors.upload(bondReductionFactorsVec);
+    scaleFactors.upload(scaleFactorsVec);
     if (force.getUseDispersionCorrection())
         dispersionCoefficient = AmoebaVdwForceImpl::calcDispersionCorrection(system, force);
     else
@@ -2240,6 +2276,7 @@ void CommonCalcAmoebaWcaDispersionForceKernel::initialize(const System& system,
     defines["RMINO"] = cc.doubleToString(force.getRmino());
     defines["RMINH"] = cc.doubleToString(force.getRminh());
     defines["AWATER"] = cc.doubleToString(force.getAwater());
+    defines["DISPOFF"] = cc.doubleToString(force.getDispoff());
     defines["SHCTD"] = cc.doubleToString(force.getShctd());
     defines["M_PI"] = cc.doubleToString(M_PI);
     ComputeProgram program = cc.compileProgram(CommonAmoebaKernelSources::amoebaWcaForce, defines);

plugins/amoeba/platforms/common/src/AmoebaCommonKernels.h

@@ -298,15 +298,19 @@ public:
      * @param force      the AmoebaGeneralizedKirkwoodForce to copy the parameters from
      */
     void copyParametersToContext(ContextImpl& context, const AmoebaGeneralizedKirkwoodForce& force);
+
 private:
     class ForceInfo;
     ComputeContext& cc;
     const System& system;
-    bool includeSurfaceArea, hasInitializedKernels;
+    bool includeSurfaceArea, tanhRescaling, hasInitializedKernels;
     int computeBornSumThreads, gkForceThreads, chainRuleThreads, ediffThreads;
     AmoebaMultipoleForce::PolarizationType polarizationType;
     std::map<std::string, std::string> defines;
     ComputeArray params;
+    ComputeArray neckRadii;
+    ComputeArray neckA;
+    ComputeArray neckB;
     ComputeArray bornSum;
     ComputeArray bornRadii;
     ComputeArray bornForce;
@@ -362,6 +366,8 @@ private:
     float currentVdwLambda;
     // Per particle alchemical flag.
     ComputeArray isAlchemical;
+    // Per particle scale factor.
+    ComputeArray scaleFactors;
 
     double dispersionCoefficient;
     ComputeArray sigmaEpsilon, atomType;

plugins/amoeba/platforms/common/src/kernels/amoebaVdwForce2.cc

@@ -7,7 +7,7 @@
     real tempForce = 0.0f;
     float2 pairSigmaEpsilon = sigmaEpsilon[atomType1+atomType2*NUM_TYPES];
     real sigma = pairSigmaEpsilon.x;
-    real epsilon = pairSigmaEpsilon.y;
+    real epsilon = pairSigmaEpsilon.y * scaleFactor1 * scaleFactor2;
     real softcore = 0.0f;
 #if VDW_ALCHEMICAL_METHOD == 1
     if (isAlchemical1 != isAlchemical2) { 

plugins/amoeba/platforms/common/src/kernels/amoebaWcaForce.cc

@@ -29,10 +29,145 @@ DEVICE void initParticleParameters(float radius, float epsilon, real* rmixo, rea
     *rmixh = 2*(rminh2*RMINH + radius2*radius) / (rminh2+radius2);
 }
 
+DEVICE real integralBeforeRMin(real eps, real r, real r2, real sk2,
+        real lik2, real lik3, real lik4, real uik2, real uik3, real uik4) {
+        return -eps * (4 * M_PI / (48 * r) * (3 * (lik4 - uik4) - 8 * r * (lik3 - uik3) + 6 * (r2 - sk2) * (lik2 - uik2)));
+}
+
+DEVICE real integralBeforeRminDerivative(real ri, real eps, real rmin, real r, real r2,
+                                         real r3, real sk, real sk2, real lik, real lik2,
+                                         real lik3, real uik, real uik2, real uik3) {
+    real dl;
+    if (ri > r - sk) {
+        dl = (-lik2 + 2 * r2 + 2 * sk2) * lik2;
+    } else {
+        dl = (-lik3 + 4 * lik2 * r - 6 * lik * r2 + 2 * lik * sk2 + 4 * r3 - 4 * r * sk2) * lik;
+    }
+    real du;
+    if (r + sk > rmin) {
+        du = -(-uik2 + 2 * r2 + 2 * sk2) * uik2;
+    } else {
+        du = -(-uik3 + 4 * uik2 * r - 6 * uik * r2 + 2 * uik * sk2 + 4 * r3 - 4 * r * sk2) * uik;
+    }
+    return -eps * M_PI * (dl + du) / (4 * r2);
+}
+
+DEVICE real integralAfterRmin(real eps, real rmin7, real r, real r2, real sk2,
+                               real lik, real lik2, real lik3, real lik4, real lik5, real lik10,
+                                real lik11, real lik12, real uik, real uik2, real uik3, real uik4,
+                                real uik5, real uik10, real uik11, real uik12) {
+    real er7 = eps * rmin7;
+    real term = 4 * M_PI / (120 * r * lik5 * uik5)
+            * (15 * uik * lik * r * (uik4 - lik4)
+               - 10 * uik2 * lik2 * (uik3 - lik3)
+               + 6 * (sk2 - r2) * (uik5 - lik5));
+    real term2 = 4 * M_PI / (2640 * r * lik12 * uik12)
+             * (120 * uik * lik * r * (uik11 - lik11)
+                - 66 * uik2 * lik2 * (uik10 - lik10)
+                + 55 * (sk2 - r2) * (uik12 - lik12));
+    real idisp = -2 * er7 * term;
+    real irep = er7 * rmin7 * term2;
+    return irep + idisp;
+}
+
+DEVICE real integralAfterRminDerivative(real ri, real eps, real rmin, real rmin7, real rmax,
+                                        real r, real r2, real r3, real sk, real sk2, real lik,
+                                        real lik2, real lik3, real lik5, real lik6, real lik12,
+                                        real lik13, real uik, real uik2, real uik3, real uik6,
+                                        real uik13) {
+    real er7 = eps * rmin7;
+    real lowerTerm = lik2 * r + r3 - r * sk2;
+    real upperTerm = uik2 * r + r3 - r * sk2;
+
+    real dl;
+    if (ri > r - sk || rmax < rmin) {
+        dl = -(-5 * lik2 + 3 * r2 + 3 * sk2) / lik5;
+    } else {
+        dl = (5 * lik3 - 33 * lik * r2 - 3 * lik * sk2 + 15 * lowerTerm) / lik6;
+    }
+    real du = -(5 * uik3 - 33 * uik * r2 - 3 * uik * sk2 + 15 * upperTerm) / uik6;
+    real de = -2 * M_PI * er7 * (dl + du) / (15 * r2);
+
+    if (ri > r - sk || rmax < rmin) {
+        dl = -(-6 * lik2 + 5 * r2 + 5 * sk2) / lik12;
+    } else {
+        dl = (6 * lik3 - 125 * lik * r2 - 5 * lik * sk2 + 60 * lowerTerm) / lik13;
+    }
+    du = -(6 * uik3 - 125 * uik * r2 - 5 * uik * sk2 + 60 * upperTerm) / uik13;
+    de += M_PI * er7 * rmin7 * (dl + du) / (60 * r2);
+    return de;
+}
+
+DEVICE real interact(real factor, real ri, real sk, real rmix, real emix,
+                     real r, real r2, real r3, real3 *force) {
+    real sum = 0;
+    // Nothing to do if the integral begins beyond r + sk (i.e. atom k does not exclude solvent)
+    if (ri < r + sk) {
+        // Zero out the derivative contribution of atom k.
+        real de = 0;
+        real sk2 = sk * sk;
+        // Compute the maximum of 1) the beginning of the integral and 2) closest edge of atom K.
+        real iStart = ri > r - sk ? ri : r - sk;
+        // Use this as the lower limit for integrating the constant eps value below Rmin.
+        real lik = iStart;
+        // Interaction with water from lik to Rmin; nothing to do if the lower limit is greater than Rmin.
+        if (lik < rmix) {
+            real lik2 = lik * lik;
+            real lik3 = lik2 * lik;
+            real lik4 = lik3 * lik;
+            // Upper limit is the minimum of Rmin and the farthest edge of atom K.
+            real uik = r + sk < rmix ? r + sk : rmix;
+            real uik2 = uik * uik;
+            real uik3 = uik2 * uik;
+            real uik4 = uik3 * uik;
+            sum = integralBeforeRMin(emix, r, r2, sk2, lik2, lik3, lik4, uik2, uik3, uik4);
+            de = integralBeforeRminDerivative(ri, emix, rmix, r, r2, r3, sk, sk2, lik, lik2, lik3, uik, uik2, uik3);
+        }
+        // Upper limit the variable part of Uwca always the farthest edge of atom K.
+        real uik = r + sk;
+        // Interaction with water beyond Rmin, from lik to uik = r + sk.
+        if (uik > rmix) {
+            // Start the integral at the max of 1) iStart and 2) Rmin.
+            lik = iStart > rmix ? iStart : rmix;
+            real lik2 = lik * lik;
+            real lik3 = lik2 * lik;
+            real lik4 = lik3 * lik;
+            real lik5 = lik4 * lik;
+            real lik6 = lik5 * lik;
+            real lik10 = lik5 * lik5;
+            real lik11 = lik10 * lik;
+            real lik12 = lik11 * lik;
+            real uik2 = uik * uik;
+            real uik3 = uik2 * uik;
+            real uik4 = uik3 * uik;
+            real uik5 = uik4 * uik;
+            real uik10 = uik5 * uik5;
+            real uik11 = uik10 * uik;
+            real uik12 = uik11 * uik;
+            real rmix3 = rmix * rmix * rmix;
+            real rmix7 = rmix3 * rmix3 * rmix;
+            sum += integralAfterRmin(emix, rmix7, r, r2, sk2, lik, lik2, lik3, lik4, lik5, lik10, lik11, lik12, uik,
+                    uik2, uik3, uik4, uik5, uik10, uik11, uik12);
+            real lik13 = lik12 * lik;
+            real uik6 = uik5 * uik;
+            real uik13 = uik12 * uik;
+            de += integralAfterRminDerivative(ri, emix, rmix, rmix7, iStart, r, r2, r3, sk, sk2, lik, lik2, lik3,
+                    lik5, lik6, lik12, lik13, uik, uik2, uik3, uik6, uik13);
+        }
+        // Increment the individual dispersion gradient components.
+        de *= factor / r;
+        (*force).x += de;
+        (*force).y += de;
+        (*force).z += de;
+    }
+    return factor * sum;
+}
+
+
 DEVICE void computeOneInteraction(AtomData atom1, AtomData atom2, real rmixo, real rmixh, real emixo, real emixh, real3* force, real* energy) {
     // get deltaR and r between 2 atoms
-    
-    *force = atom2.pos - atom1.pos;
+
+    *force = atom1.pos - atom2.pos;
     real r2 = dot(*force, *force);
     if (r2 <= 0) {
         *force = make_real3(0);
@@ -42,152 +177,28 @@ DEVICE void computeOneInteraction(AtomData atom1, AtomData atom2, real rmixo, re
     real rI = RSQRT(r2);
     real r = RECIP(rI);
 
-    real sk = atom2.radius*SHCTD;
-    real sk2 = sk*sk;
-    if (atom1.radius >= (r+sk)) {
-        *force = make_real3(0);
-        *energy = 0;
-        return;
-    }
-
-    real rmax = atom1.radius > (r - sk) ? atom1.radius : (r - sk);
-    real lik = rmax;
-    real lik2 = lik*lik;
-    real lik3 = lik2*lik;
-    real lik4 = lik2*lik2;
- 
-    real uik = (r+sk) < rmixo ? (r+sk) : rmixo;
-    real uik2 = uik*uik;
-    real uik3 = uik2*uik;
-    real uik4 = uik2*uik2;
-
-    real term = 4*M_PI/(48*r)*(3*(lik4-uik4) - 8*r*(lik3-uik3) + 6*(r2-sk2)*(lik2-uik2));
+    real xr = (*force).x;
+    real yr = (*force).y;
+    real zr = (*force).z;
+    real r3 = r2 * r;
 
-    real r3 = r2*r;
-    real dl1 = lik2*(-lik2 + 2*(r2 + sk2));
-    real dl2 = lik*(-lik3 + 4*lik2*r - 6*lik*r2 + 2*lik*sk2 + 4*r3 - 4*r*sk2);
-    real dl = atom1.radius > (r-sk)? dl1 : dl2;
+    real sk = atom2.radius * SHCTD;
 
-    real du1 = uik2*(-uik2 + 2*(r2 + sk2));
-    real du2 = uik*(-uik3 + 4*uik2*r - 2*uik*(3*r2 - sk2) + 4*r*(r2 - sk2));
-    real du = (r+sk) > rmixo ? -du1 : -du2;
+    // Start of integration of dispersion for atom i with water oxygen.
+    real riO = rmixo * 0.5f + DISPOFF;
+    real nO = 1.0f;
 
-    real mask2 = lik < rmixo ? 1 : 0;
-    real sum = -mask2*(emixo*term);
-    real de = -mask2*emixo*M_PI*(dl+du)/(4*r2);
+    // Start of integration of dispersion for atom i with water hydrogen.
+    real riH = rmixh * 0.5f + DISPOFF;
+    real nH = 2.0f;
 
-    uik = (r+sk) < rmixh ? (r+sk) : rmixh;
-    uik2 = uik*uik;
-    uik3 = uik2*uik;
-    uik4 = uik2*uik2;
+    *force = make_real3(0);
+    *energy = interact(nO, riO, sk, rmixo, emixo, r, r2, r3, force) +
+             interact(nH, riH, sk, rmixh, emixh, r, r2, r3, force);
 
-    term = (M_PI)/ (12*r) * (3*(lik4-uik4) - 8*r*(lik3-uik3) + 6*(r2-sk2)*(lik2-uik2));
-
-    dl1 = lik2*(-lik2 + 2*r2 + 2*sk2);
-    dl2 = lik*(-lik3 + 4*lik2*r - 6*lik*r2 + 2*lik*sk2 + 4*r3 - 4*r*sk2);
-    dl = atom1.radius > (r-sk) ? dl1 : dl2;
-
-    du1 = -uik2*(-uik2 + 2*r2 + 2*sk2);
-    du2 = -uik*(-uik3 + 4*uik2*r - 6*uik*r2 + 2*uik*sk2 + 4*r3 - 4*r*sk2);
-    du = (r+sk) > rmixh ? du1 : du2;
-
-    mask2 = lik < rmixh ? 1 : 0;
-    sum -= mask2*(2*emixh*term);
-    de -= mask2*(2*emixh*M_PI*(dl+du)/(4*r2));
-
-    uik = r + sk;
-    uik2 = uik*uik;
-    uik3 = uik2*uik;
-    uik4 = uik2*uik2;
-    real uik5 = uik4*uik;
-    real uik6 = uik3*uik3;
-    real uik10 = uik5*uik5;
-    real uik11 = uik10*uik;
-    real uik12 = uik6*uik6;
-    real uik13 = uik12*uik;
-
-    lik = rmax > rmixo ? rmax : rmixo;
-    lik2 = lik*lik;
-    lik3 = lik2*lik;
-    lik4 = lik2*lik2;
-    real lik5 = lik4*lik;
-    real lik6 = lik3*lik3;
-    real lik10 = lik5*lik5;
-    real lik11 = lik10*lik;
-    real lik12 = lik6*lik6;
-    real lik13 = lik12*lik;
-
-    term = 4*M_PI/(120*r*lik5*uik5)*(15*uik*lik*r*(uik4-lik4) - 10*uik2*lik2*(uik3-lik3) + 6*(sk2-r2)*(uik5-lik5));
-    dl1 = (-5*lik2 + 3*r2 + 3*sk2)/lik5;
-    dl2 = (5*lik3 - 33*lik*r2 - 3*lik*sk2 + 15*(lik2*r+r3-r*sk2))/lik6;
-    dl = (atom1.radius > (r-sk)) || (rmax < rmixo) ? -dl1 : dl2;
-
-    du = (-5*uik3 + 33*uik*r2 + 3*uik*sk2 - 15*(uik2*r+r3-r*sk2))/uik6;
-
-    real rmixo7 = rmixo*rmixo*rmixo;
-    rmixo7 = rmixo7*rmixo7*rmixo;
-    real ao = emixo*rmixo7;
-
-    real idisp = -2*ao*term;
-    mask2 = uik > rmixo ? 1 : 0;
-
-    de -= mask2*(2*ao*M_PI*(dl + du)/(15*r2));
-
-    term = 4*M_PI/(2640*r*lik12*uik12) * (120*uik*lik*r*(uik11-lik11) - 66*uik2*lik2*(uik10-lik10) + 55*(sk2-r2)*(uik12-lik12));
-
-    dl1 = (6*lik2 - 5*r2 - 5*sk2)/lik12;
-    dl2 = (6*lik3 - 125*lik*r2 - 5*lik*sk2 + 60*(lik2*r+r3-r*sk2))/lik13;
-    dl = (atom1.radius > (r-sk)) || (rmax < rmixo) ? dl1 : dl2;
-
-    du = (-6*uik3 + 125*uik*r2 + 5*uik*sk2 - 60*(uik2*r+r3-r*sk2))/uik13;
-
-    de += mask2*(ao*rmixo7*M_PI*(dl + du)/(60*r2));
-    real irep = ao*rmixo7*term;
-    sum += mask2*(irep + idisp);
-
-    lik = rmax > rmixh ? rmax : rmixh;
-    lik2 = lik*lik;
-    lik3 = lik2*lik;
-    lik4 = lik2*lik2;
-    lik5 = lik4*lik;
-    lik6 = lik3*lik3;
-    lik10 = lik5*lik5;
-    lik11 = lik10*lik;
-    lik12 = lik6*lik6;
-    lik13 = lik12*lik;
-
-    term = 4*M_PI / (120*r*lik5*uik5) * (15*uik*lik*r*(uik4-lik4) - 10*uik2*lik2*(uik3-lik3) + 6*(sk2-r2)*(uik5-lik5));
-
-    dl1 = (-5*lik2 + 3*r2 + 3*sk2)/lik5;
-    dl2 = (5*lik3 - 33*lik*r2 - 3*lik*sk2+ 15*(lik2*r+r3-r*sk2))/lik6;
-    dl = (atom1.radius > (r-sk)) || (rmax < rmixh) ? -dl1 : dl2;
-
-    du = -(5*uik3 - 33*uik*r2 - 3*uik*sk2 + 15*(uik2*r+r3-r*sk2))/uik6;
-
-    real rmixh7 = rmixh*rmixh*rmixh;
-    rmixh7 = rmixh7*rmixh7*rmixh;
-    real ah = emixh*rmixh7;
-
-    idisp = -4*ah*term;
-
-    mask2 = uik > rmixh ? 1 : 0;
-    de -= mask2*(4*ah*M_PI*(dl + du)/(15*r2));
-
-    term = 4*M_PI / (2640*r*lik12*uik12) * (120*uik*lik*r*(uik11-lik11) - 66*uik2*lik2*(uik10-lik10) + 55*(sk2-r2)*(uik12-lik12));
-
-    dl1 = -(-6*lik2 + 5*r2 + 5*sk2)/lik12;
-    dl2 = (6*lik3 - 125*lik*r2 - 5*lik*sk2 + 60*(lik2*r+r3-r*sk2))/lik13;
-    dl = ((atom1.radius > (r-sk)) || (rmax < rmixh)) ? dl1 : dl2;
-
-    du = -(6*uik3 - 125*uik*r2 -5*uik*sk2 + 60*(uik2*r+r3-r*sk2))/uik13;
-    irep = 2*ah*rmixh7*term;
-
-    de += mask2*(ah*rmixh7*M_PI*(dl+du)/(30*r2));
-    sum += mask2*(irep+idisp);
-
-    *energy = sum;
-    de *= -AWATER*rI;
-    *force *= de;
+    (*force).x *= AWATER * xr;
+    (*force).y *= AWATER * yr;
+    (*force).z *= AWATER * zr;
 }
 
 /**
@@ -263,4 +274,4 @@ KERNEL void computeWCAForce(GLOBAL mm_ulong* RESTRICT forceBuffers, GLOBAL mixed
         pos++;
     } while (pos < end);
     energyBuffer[GLOBAL_ID] -= AWATER*energy;
-}
+}
\ No newline at end of file

plugins/amoeba/platforms/reference/src/AmoebaReferenceKernels.cpp

@@ -294,19 +294,18 @@ AmoebaReferenceMultipoleForce* ReferenceCalcAmoebaMultipoleForceKernel::setupAmo
         amoebaReferenceGeneralizedKirkwoodForce->setSurfaceAreaFactor(gkKernel->getSurfaceAreaFactor());
         amoebaReferenceGeneralizedKirkwoodForce->setIncludeCavityTerm(gkKernel->getIncludeCavityTerm());
         amoebaReferenceGeneralizedKirkwoodForce->setDirectPolarization(gkKernel->getDirectPolarization());
-
-        vector<double> parameters; 
-        gkKernel->getAtomicRadii(parameters);
-        amoebaReferenceGeneralizedKirkwoodForce->setAtomicRadii(parameters);
-
-        gkKernel->getScaleFactors(parameters);
-        amoebaReferenceGeneralizedKirkwoodForce->setScaleFactors(parameters);
-
-        gkKernel->getCharges(parameters);
-        amoebaReferenceGeneralizedKirkwoodForce->setCharges(parameters);
+        amoebaReferenceGeneralizedKirkwoodForce->setTanhRescaling(gkKernel->getTanhRescaling());
+        double beta0, beta1, beta2;
+        gkKernel->getTanhParameters(beta0, beta1, beta2);
+        amoebaReferenceGeneralizedKirkwoodForce->setTanhParameters(beta0, beta1, beta2);
+        amoebaReferenceGeneralizedKirkwoodForce->setDescreenOffset(gkKernel->getDescreenOffset());
+        amoebaReferenceGeneralizedKirkwoodForce->setAtomicRadii(gkKernel->getAtomicRadii());
+        amoebaReferenceGeneralizedKirkwoodForce->setScaleFactors(gkKernel->getScaleFactors());
+        amoebaReferenceGeneralizedKirkwoodForce->setCharges(gkKernel->getCharges());
+        amoebaReferenceGeneralizedKirkwoodForce->setDescreenRadii(gkKernel->getDescreenRadii());
+        amoebaReferenceGeneralizedKirkwoodForce->setNeckFactors(gkKernel->getNeckFactors());
 
         // calculate Grycuk Born radii
-
         vector<Vec3>& posData   = extractPositions(context);
         amoebaReferenceGeneralizedKirkwoodForce->calculateGrycukBornRadii(posData);
 
@@ -534,6 +533,20 @@ int ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::getIncludeCavityTerm() co
     return includeCavityTerm;
 }
 
+bool ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::getTanhRescaling() const {
+    return tanhRescaling;
+}
+
+void ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::getTanhParameters(double &b0, double &b1, double &b2) const {
+    b0 = beta0;
+    b1 = beta1;
+    b2 = beta2;
+}
+
+double ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::getDescreenOffset() const {
+    return descreenOffset;
+}
+
 int ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::getDirectPolarization() const {
     return directPolarization;
 }
@@ -558,19 +571,24 @@ double ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::getSurfaceAreaFactor()
     return surfaceAreaFactor;
 }
 
-void ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::getAtomicRadii(vector<double>& outputAtomicRadii) const {
-    outputAtomicRadii.resize(atomicRadii.size());
-    copy(atomicRadii.begin(), atomicRadii.end(), outputAtomicRadii.begin());
+const vector<double>& ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::getAtomicRadii() const {
+    return atomicRadii;
+}
+
+const vector<double>& ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::getCharges() const {
+    return charges;
+}
+
+const vector<double>& ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::getScaleFactors() const {
+    return scaleFactors;
 }
 
-void ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::getScaleFactors(vector<double>& outputScaleFactors) const {
-    outputScaleFactors.resize(scaleFactors.size());
-    copy(scaleFactors.begin(), scaleFactors.end(), outputScaleFactors.begin());
+const vector<double>& ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::getDescreenRadii() const {
+    return descreenRadii;
 }
 
-void ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::getCharges(vector<double>& outputCharges) const {
-    outputCharges.resize(charges.size());
-    copy(charges.begin(), charges.end(), outputCharges.begin());
+const vector<double>& ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::getNeckFactors() const {
+    return neckFactors;
 }
 
 void ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::initialize(const System& system, const AmoebaGeneralizedKirkwoodForce& force) {
@@ -593,13 +611,13 @@ void ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::initialize(const System&
     numParticles = system.getNumParticles();
 
     for (int ii = 0; ii < numParticles; ii++) {
-
-        double particleCharge, particleRadius, scalingFactor;
-        force.getParticleParameters(ii, particleCharge, particleRadius, scalingFactor);
+        double particleCharge, particleRadius, scalingFactor, descreenRadius, neckFactor;
+        force.getParticleParameters(ii, particleCharge, particleRadius, scalingFactor, descreenRadius, neckFactor);
         atomicRadii.push_back(particleRadius);
         scaleFactors.push_back(scalingFactor);
         charges.push_back(particleCharge);
-
+        descreenRadii.push_back(descreenRadius);
+        neckFactors.push_back(neckFactor);
         // Make sure the charge matches the one specified by the AmoebaMultipoleForce.
 
         double charge2, thole, damping, polarity;
@@ -610,13 +628,20 @@ void ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::initialize(const System&
             throw OpenMMException("AmoebaGeneralizedKirkwoodForce and AmoebaMultipoleForce must specify the same charge for every atom.");
         }
 
-    }   
+    }
     includeCavityTerm  = force.getIncludeCavityTerm();
+    tanhRescaling      = force.getTanhRescaling();
+    double b0, b1, b2;
+    force.getTanhParameters(b0, b1, b2);
+    beta0 = b0;
+    beta1 = b1;
+    beta2 = b2;
+    descreenOffset     = force.getDescreenOffset();
     soluteDielectric   = force.getSoluteDielectric();
     solventDielectric  = force.getSolventDielectric();
-    dielectricOffset   = 0.009;
-    probeRadius        = force.getProbeRadius(), 
-    surfaceAreaFactor  = force.getSurfaceAreaFactor(); 
+    dielectricOffset   = force.getDielectricOffset();
+    probeRadius        = force.getProbeRadius();
+    surfaceAreaFactor  = force.getSurfaceAreaFactor();
     directPolarization = amoebaMultipoleForce->getPolarizationType() == AmoebaMultipoleForce::Direct ? 1 : 0;
 }
 
@@ -632,11 +657,13 @@ void ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::copyParametersToContext(
     // Record the values.
 
     for (int i = 0; i < numParticles; ++i) {
-        double particleCharge, particleRadius, scalingFactor;
-        force.getParticleParameters(i, particleCharge, particleRadius, scalingFactor);
+        double particleCharge, particleRadius, scalingFactor, descreenRadius, neckFactor;
+        force.getParticleParameters(i, particleCharge, particleRadius, scalingFactor, descreenRadius, neckFactor);
         atomicRadii[i] = particleRadius;
         scaleFactors[i] = scalingFactor;
         charges[i] = particleCharge;
+        descreenRadii[i] = descreenRadius;
+        neckFactors[i] = neckFactor;
     }
 }
 

plugins/amoeba/platforms/reference/src/AmoebaReferenceKernels.h

@@ -318,6 +318,25 @@ public:
      */
     int getIncludeCavityTerm() const;
 
+    /**
+     *  Get the 'tanh rescaling' flag.
+     *
+     *  @return tanhRescaling
+     */
+    bool getTanhRescaling() const;
+
+    /**
+     * Get Tanh parameters beta0, beta1 and beta2.
+     */
+    void getTanhParameters(double& b0, double& b1, double& b2) const;
+
+    /**
+     *  Get the descreen offset for computing effective Born radii.
+     *
+     *  @return descreenOffset
+     */
+    double getDescreenOffset() const;
+
     /**
      *  Get the number of particles.
      *
@@ -350,10 +369,9 @@ public:
     double getSolventDielectric() const;
 
     /**
-     *  Get the dielectric offset.
+     *  Get the dielectric offset for the cavity term.
      *
      *  @return dielectricOffset
-     *
      */
     double getDielectricOffset() const;
 
@@ -375,27 +393,28 @@ public:
 
     /**
      *  Get the vector of particle radii.
-     *
-     *  @param atomicRadii vector of atomic radii
-     *
      */
-    void getAtomicRadii(std::vector<double>& atomicRadii) const;
+    const vector<double>& getAtomicRadii() const;
 
     /**
      *  Get the vector of scale factors.
-     *
-     *  @param scaleFactors vector of scale factors
-     *
      */
-    void getScaleFactors(std::vector<double>& scaleFactors) const;
+    const vector<double>& getScaleFactors() const;
 
     /**
      *  Get the vector of charges.
-     *
-     *  @param charges vector of charges
-     *
      */
-    void getCharges(std::vector<double>& charges) const;
+    const vector<double>& getCharges() const;
+
+    /**
+     *  Get the vector of descreening radii.
+     */
+    const vector<double>& getDescreenRadii() const;
+
+    /**
+     *  Get the vector of neck scaling factors.
+     */
+    const vector<double>& getNeckFactors() const;
 
     /**
      * Copy changed parameters over to a context.
@@ -411,12 +430,19 @@ private:
     std::vector<double> atomicRadii;
     std::vector<double> scaleFactors;
     std::vector<double> charges;
+    std::vector<double> descreenRadii;
+    std::vector<double> neckFactors;
     double soluteDielectric;
     double solventDielectric;
     double dielectricOffset;
     double probeRadius;
     double surfaceAreaFactor;
     int includeCavityTerm;
+    bool tanhRescaling;
+    double beta0;
+    double beta1;
+    double beta2;
+    double descreenOffset;
     int directPolarization;
     const System& system;
 };

plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceGeneralizedKirkwoodForce.cpp

@@ -22,20 +22,25 @@
  */
 
 #include "AmoebaReferenceGeneralizedKirkwoodForce.h"
+#include "openmm/internal/AmoebaGeneralizedKirkwoodForceImpl.h"
 #include <cmath>
 
 using std::vector;
 using namespace OpenMM;
 
 AmoebaReferenceGeneralizedKirkwoodForce::AmoebaReferenceGeneralizedKirkwoodForce() : _numParticles(0),
-                                                                                      _includeCavityTerm(1),
-                                                                                      _directPolarization(0),
-                                                                                      _soluteDielectric(1.0),
-                                                                                      _solventDielectric(78.3),
-                                                                                      _dielectricOffset(0.009),
-                                                                                      _probeRadius(0.14),
-                                                                                      _surfaceAreaFactor(0.0054) {
-
+                                                                                     _includeCavityTerm(1),
+                                                                                     _directPolarization(0),
+                                                                                     _soluteDielectric(1.0),
+                                                                                     _solventDielectric(78.3),
+                                                                                     _dielectricOffset(0.009),
+                                                                                     _probeRadius(0.14),
+                                                                                     _surfaceAreaFactor(0.0054),
+                                                                                     _tanhRescaling(false),
+                                                                                     _beta0(0.9563),
+                                                                                     _beta1(0.2578),
+                                                                                     _beta2(0.0810),
+                                                                                     _descreenOffset(0.0) {
 }
 
 void AmoebaReferenceGeneralizedKirkwoodForce::setNumParticles(int numParticles) {
@@ -54,6 +59,41 @@ int AmoebaReferenceGeneralizedKirkwoodForce::getIncludeCavityTerm() const {
     return _includeCavityTerm;
 }
 
+void AmoebaReferenceGeneralizedKirkwoodForce::setTanhRescaling(bool tanhRescaling) {
+    _tanhRescaling = tanhRescaling;
+}
+
+bool AmoebaReferenceGeneralizedKirkwoodForce::getTanhRescaling() const {
+    return _tanhRescaling;
+}
+
+/**
+ * Get Tanh parameters beta0, beta1 and beta2.
+ */
+void AmoebaReferenceGeneralizedKirkwoodForce::getTanhParameters(double& b0, double& b1, double& b2) const {
+    b0 = _beta0;
+    b1 = _beta1;
+    b2 = _beta2;
+}
+
+/**
+ * Set the flag signaling whether the solute integral is rescaled by a Tanh function
+ * to account for interstitial spaces.
+*/
+void AmoebaReferenceGeneralizedKirkwoodForce::setTanhParameters(double b0, double b1, double b2) {
+    _beta0 = b0;
+    _beta1 = b1;
+    _beta2 = b2;
+}
+
+void AmoebaReferenceGeneralizedKirkwoodForce::setDescreenOffset(double descreenOffset) {
+    _descreenOffset  = descreenOffset;
+}
+
+double AmoebaReferenceGeneralizedKirkwoodForce::getDescreenOffset() const {
+    return _descreenOffset;
+}
+
 void AmoebaReferenceGeneralizedKirkwoodForce::setDirectPolarization(int directPolarization) {
     _directPolarization = directPolarization;
 }
@@ -127,77 +167,145 @@ void AmoebaReferenceGeneralizedKirkwoodForce::setCharges(const vector<double>& c
     copy(charges.begin(), charges.end(), _charges.begin());
 }
 
+void AmoebaReferenceGeneralizedKirkwoodForce::setDescreenRadii(const vector<double> &descreenRadii) {
+    _descreenRadii.resize(descreenRadii.size());
+    copy(descreenRadii.begin(), descreenRadii.end(), _descreenRadii.begin());
+}
+
+void AmoebaReferenceGeneralizedKirkwoodForce::getDescreenRadii(vector<double> &descreenRadii) const {
+    descreenRadii.resize(_atomicRadii.size());
+    copy(_descreenRadii.begin(), _descreenRadii.end(), descreenRadii.begin());
+}
+
+void AmoebaReferenceGeneralizedKirkwoodForce::setNeckFactors(const vector<double> &neckFactors) {
+    _neckFactors.resize(neckFactors.size());
+    copy(neckFactors.begin(), neckFactors.end(), _neckFactors.begin());
+}
+
+void AmoebaReferenceGeneralizedKirkwoodForce::getNeckFactors(vector<double> &neckFactors) const {
+    neckFactors.resize(_neckFactors.size());
+    copy(_neckFactors.begin(), _neckFactors.end(), neckFactors.begin());
+}
+
 void AmoebaReferenceGeneralizedKirkwoodForce::getGrycukBornRadii(vector<double>& bornRadii) const {
     bornRadii.resize(_bornRadii.size());
     copy(_bornRadii.begin(), _bornRadii.end(), bornRadii.begin());
 }
 
-void AmoebaReferenceGeneralizedKirkwoodForce::calculateGrycukBornRadii(const vector<Vec3>& particlePositions) {
+void AmoebaReferenceGeneralizedKirkwoodForce::getSoluteIntegral(vector<double> &soluteIntegral) const {
+    soluteIntegral.resize(_soluteIntegral.size());
+    copy(_soluteIntegral.begin(), _soluteIntegral.end(), soluteIntegral.begin());
+}
 
-    const double bigRadius = 1000.0;
+void AmoebaReferenceGeneralizedKirkwoodForce::calculateGrycukBornRadii(const vector<Vec3> &particlePositions) {
+
+    // Set the radius to 30 Angstroms (3 nm) if either the base radius is zero, or the
+    // descreening integral is negative.
+    const double MAX_RADIUS = 3.0;
+    const double RECIP_MAX_RADIUS3 = pow(MAX_RADIUS, -3.0);
+    const double PI4_3 = 4.0 / 3.0 * M_PI;
+    const double INVERSE_PI4_3 = 1.0 / PI4_3;
+    const double ONE_THIRD = 1.0 / 3.0;
 
     _bornRadii.resize(_numParticles);
+    _soluteIntegral.resize(_numParticles);
+
     for (unsigned int ii = 0; ii < _numParticles; ii++) {
 
         if (_atomicRadii[ii] <= 0.0) {
-            _bornRadii[ii] = bigRadius;
+            _bornRadii[ii] = MAX_RADIUS;
             continue;
         }
 
+        double integralStartI = max(_atomicRadii[ii], _descreenRadii[ii]) + _descreenOffset;
+
         double bornSum = 0.0;
+        double neckSum = 0.0;
         for (unsigned int jj = 0; jj < _numParticles; jj++) {
 
-            if (ii == jj || _atomicRadii[jj] < 0.0)continue;
-          
-            double xr       = particlePositions[jj][0] - particlePositions[ii][0];
-            double yr       = particlePositions[jj][1] - particlePositions[ii][1];
-            double zr       = particlePositions[jj][2] - particlePositions[ii][2];
+            double sk = _descreenRadii[jj] * _scaleFactors[jj];
+
+            if (ii == jj || integralStartI <= 0.0 || sk <= 0.0) continue;
 
-            double r2       = xr*xr + yr*yr + zr*zr;
-            double r        = sqrt(r2);
+            double xr = particlePositions[jj][0] - particlePositions[ii][0];
+            double yr = particlePositions[jj][1] - particlePositions[ii][1];
+            double zr = particlePositions[jj][2] - particlePositions[ii][2];
 
-            double sk       = _atomicRadii[jj]*_scaleFactors[jj];
+            double r2 = xr * xr + yr * yr + zr * zr;
+            double r = sqrt(r2);
 
             // If atom ii engulfs the descreening atom, then continue.
-            if (_atomicRadii[ii] > r + sk) continue;
-
-            double sk2      = sk*sk;
-
-            if ((_atomicRadii[ii] + r) < sk) {
-                double lik       = _atomicRadii[ii];
-                double uik       = sk - r;  
-                double lik3      = lik*lik*lik;
-                double uik3      = uik*uik*uik;
-                bornSum             -= (1.0/uik3 - 1.0/lik3);
-            }   
-        
-            double uik = r + sk; 
+            if (integralStartI > r + sk) continue;
+
+            double sk2 = sk * sk;
+
+            if ((integralStartI + r) < sk) {
+                double lik = integralStartI;
+                double uik = sk - r;
+                double lik3 = lik * lik * lik;
+                double uik3 = uik * uik * uik;
+                bornSum -= (1.0 / uik3 - 1.0 / lik3);
+            }
+
+            double uik = r + sk;
             double lik;
-            if ((_atomicRadii[ii] + r) < sk) {
-                lik = sk - r;  
-            } else if (r < (_atomicRadii[ii] + sk)) {
-                lik = _atomicRadii[ii];
+            if ((integralStartI + r) < sk) {
+                lik = sk - r;
+            } else if (r < (integralStartI + sk)) {
+                lik = integralStartI;
             } else {
-                lik = r - sk; 
-            }   
-        
-            double l2          = lik*lik; 
-            double l4          = l2*l2;
-            double lr          = lik*r;
-            double l4r         = l4*r;
-        
-            double u2          = uik*uik;
-            double u4          = u2*u2;
-            double ur          = uik*r;
-            double u4r         = u4*r;
-        
-            double term        = (3.0*(r2-sk2) + 6.0*u2 - 8.0*ur)/u4r - (3.0*(r2-sk2) + 6.0*l2 - 8.0*lr)/l4r;
-            bornSum           += term/16.0;
-        
+                lik = r - sk;
+            }
+
+            double l2 = lik * lik;
+            double l4 = l2 * l2;
+            double lr = lik * r;
+            double l4r = l4 * r;
+
+            double u2 = uik * uik;
+            double u4 = u2 * u2;
+            double ur = uik * r;
+            double u4r = u4 * r;
+
+            double term =
+                    (3.0 * (r2 - sk2) + 6.0 * u2 - 8.0 * ur) / u4r
+                    - (3.0 * (r2 - sk2) + 6.0 * l2 - 8.0 * lr) / l4r;
+            bornSum += term / 16.0;
+
+            double mixedNeckScale = 0.5 * (_neckFactors[ii] + _neckFactors[jj]);
+            if (mixedNeckScale > 0.0) {
+                neckSum += AmoebaGeneralizedKirkwoodForceImpl::neckDescreen(r, integralStartI, _descreenRadii[jj], mixedNeckScale);
+            }
+        }
+
+        bornSum = bornSum * PI4_3 + neckSum;
+        _soluteIntegral[ii] = bornSum;
+
+        double baseRadiusI3 = _atomicRadii[ii] * _atomicRadii[ii] * _atomicRadii[ii];
+
+        if (_tanhRescaling) {
+            // Set up tanh function components
+            double rhoi3Psi = baseRadiusI3 * _soluteIntegral[ii];
+            double rhoi6Psi2 = rhoi3Psi * rhoi3Psi;
+            double rhoi9Psi3 = rhoi6Psi2 * rhoi3Psi;
+            // If the output of the tanh function is 1.0, then the Born radius will be MaxBornRadius
+            double tanh_constant = PI4_3 * ((1.0 / baseRadiusI3) - RECIP_MAX_RADIUS3);
+            bornSum = tanh_constant * tanh(_beta0 * rhoi3Psi - _beta1 * rhoi6Psi2 + _beta2 * rhoi9Psi3);
+        }
+
+        bornSum = PI4_3 / baseRadiusI3 - bornSum;
+
+        _bornRadii[ii] = (bornSum <= 0.0) ? MAX_RADIUS : pow(INVERSE_PI4_3 * bornSum, -ONE_THIRD);
+
+        // Born radius must be at least as large as its base radius.
+        if (_bornRadii[ii] < _atomicRadii[ii]) {
+            _bornRadii[ii] = _atomicRadii[ii];
+        }
+
+        // Maximum Born radius is 30.0 Angstroms.
+        if (_bornRadii[ii] > MAX_RADIUS) {
+            _bornRadii[ii] = MAX_RADIUS;
         }
-        bornSum        = 1.0/(_atomicRadii[ii]*_atomicRadii[ii]*_atomicRadii[ii]) - bornSum;
-        _bornRadii[ii] = (bornSum <= 0.0) ? bigRadius : pow(bornSum, -1.0/3.0);
     }
 
-    return;
 }

plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceGeneralizedKirkwoodForce.h

@@ -81,6 +81,53 @@ public:
      */
     void setIncludeCavityTerm(int includeCavityTerm);
 
+    /**
+     * Get the tanhRescaling flag, which controls whether tanh rescaling of the solute integral is performed.
+     *
+     * @return tanhRescaling The value tanhRescaling flag.
+     */
+    bool getTanhRescaling() const;
+
+    /**
+     *  Set the tanhRescaling flag.
+     *
+     *  @param tanhRescaling if true, tanh rescaling of the solute integral is performed.
+     */
+    void setTanhRescaling(bool tanhRescaling);
+
+    /**
+     * Get Tanh parameters b0, b1 and b2.
+     *
+     * @param b0 The first tanh parameter.
+     * @param b1 The second tanh parameter.
+     * @param b2 The third tanh parameter.
+    */
+    void getTanhParameters(double& b0, double& b1, double& b2) const;
+
+    /**
+     * Set the the Tanh function parameters that account for interstitial spaces.
+     *
+     * @param b0 The first tanh parameter.
+     * @param b1 The second tanh parameter.
+     * @param b2 The third tanh parameter.
+     */
+    void setTanhParameters(double b0, double b1, double b2);
+
+    /**
+     *  Get the descreen offset for calculating Born radii.
+     *
+     *  @return descreenOffset
+     */
+    double getDescreenOffset() const;
+
+    /**
+     * Set descreen offset for calculating Born radii.
+     *
+     * @param descreenOffset descreen offset
+     *
+     */
+    void setDescreenOffset(double descreenOffset);
+
     /**
      *  Get directPolarization flag 
      *
@@ -128,7 +175,7 @@ public:
     void setSolventDielectric(double solventDielectric);
 
     /**
-     *  Get dielectric offset
+     *  Get dielectric offset for the cavity term.
      *
      *  @return dielectricOffset
      *
@@ -136,7 +183,7 @@ public:
     double getDielectricOffset() const;
 
     /**
-     * Set dielectric offset
+     * Set dielectric offset for the cavity term.
      *
      * @param dielectricOffset dielectric offset
      *
@@ -215,6 +262,36 @@ public:
      */
     void setCharges(const vector<double>& charges);
 
+    /**
+     * Set descreen radii
+     *
+     * @param descreenRadii input vector of descreen radii
+     *
+    */
+    void setDescreenRadii(const vector<double>& descreenRadii);
+
+    /**
+     * Get descreen radii
+     *
+     * @param descreenRadii output vector of descreen radii
+     *
+     */
+    void getDescreenRadii(vector<double>& descreenRadii) const;
+
+    /**
+     * Set neck scale factors
+     *
+     * @param neckFactors input vector of neck scale factors
+     */
+    void setNeckFactors(const vector<double>& neckFactors);
+
+    /**
+     * Get neck scale factors
+     *
+     * @param neckFactors output vector of neck scale factors
+     */
+    void getNeckFactors(vector<double>& neckFactors) const;
+
     /**
      * Calculate Grycuk Born radii
      *
@@ -229,13 +306,23 @@ public:
      * @param bornRadii vector of Born radii
      *
      */
-    void getGrycukBornRadii(vector<double>& bornRadii) const;     
+    void getGrycukBornRadii(vector<double>& bornRadii) const;
+
+    /**
+     * Get Grycik Solute Integral (must have called calculateGrycukBornRadii())
+     *
+     * @param bornRadii vector of Born radii
+     *
+     */
+    void getSoluteIntegral(vector<double>& soluteIntegral) const;
 
 private:
 
     int _numParticles;
     int _includeCavityTerm;
     int _directPolarization;
+    bool _tanhRescaling;
+    double _descreenOffset;
 
     double _soluteDielectric;
     double _solventDielectric;
@@ -243,10 +330,17 @@ private:
     double _probeRadius;
     double _surfaceAreaFactor;
 
+    double _beta0;
+    double _beta1;
+    double _beta2;
+
     std::vector<double> _atomicRadii;
     std::vector<double> _scaleFactors;
     std::vector<double> _charges;
+    std::vector<double> _descreenRadii;
+    std::vector<double> _neckFactors;
 
+    std::vector<double> _soluteIntegral;
     std::vector<double> _bornRadii;
 
 };

plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceMultipoleForce.cpp

@@ -23,6 +23,7 @@
  */
 
 #include "AmoebaReferenceMultipoleForce.h"
+#include "openmm/internal/AmoebaGeneralizedKirkwoodForceImpl.h"
 #include "SimTKOpenMMRealType.h"
 #include "jama_svd.h"
 #include <algorithm>
@@ -35,6 +36,7 @@
 // make sure that erf() and erfc() are defined.
 #include "openmm/internal/MSVC_erfc.h"
 
+
 using std::vector;
 using namespace OpenMM;
 
@@ -2117,8 +2119,7 @@ AmoebaReferenceMultipoleForce::UpdateInducedDipoleFieldStruct::UpdateInducedDipo
 AmoebaReferenceGeneralizedKirkwoodMultipoleForce::AmoebaReferenceGeneralizedKirkwoodMultipoleForce(AmoebaReferenceGeneralizedKirkwoodForce* amoebaReferenceGeneralizedKirkwoodForce) :
                AmoebaReferenceMultipoleForce(NoCutoff),
                _amoebaReferenceGeneralizedKirkwoodForce(amoebaReferenceGeneralizedKirkwoodForce),
-               _gkc(2.455)
-{
+               _gkc(2.455), _beta0(0.9563), _beta1(0.2578), _beta2(0.0810) {
 
     double solventDielectric  = _amoebaReferenceGeneralizedKirkwoodForce->getSolventDielectric();
 
@@ -2127,17 +2128,19 @@ AmoebaReferenceGeneralizedKirkwoodMultipoleForce::AmoebaReferenceGeneralizedKirk
     _fq                           =  3.0*(1.0 - solventDielectric)/(2.0 + 3.0*solventDielectric);
 
     _amoebaReferenceGeneralizedKirkwoodForce->getGrycukBornRadii(_bornRadii);
+    _amoebaReferenceGeneralizedKirkwoodForce->getSoluteIntegral(_soluteIntegral);
     _amoebaReferenceGeneralizedKirkwoodForce->getAtomicRadii(_atomicRadii);
-    _amoebaReferenceGeneralizedKirkwoodForce->getScaleFactors(_scaledRadii);
-
-    _includeCavityTerm            = _amoebaReferenceGeneralizedKirkwoodForce->getIncludeCavityTerm();
-    _probeRadius                  = _amoebaReferenceGeneralizedKirkwoodForce->getProbeRadius();
-    _surfaceAreaFactor            = _amoebaReferenceGeneralizedKirkwoodForce->getSurfaceAreaFactor();
-    _dielectricOffset             = _amoebaReferenceGeneralizedKirkwoodForce->getDielectricOffset();
-
-    for (unsigned int ii = 0; ii < _scaledRadii.size(); ii++) {
-        _scaledRadii[ii] *= _atomicRadii[ii];
-    }
+    _amoebaReferenceGeneralizedKirkwoodForce->getScaleFactors(_scaleFactors);
+    _amoebaReferenceGeneralizedKirkwoodForce->getDescreenRadii(_descreenRadii);
+    _amoebaReferenceGeneralizedKirkwoodForce->getNeckFactors(_neckFactors);
+
+    _includeCavityTerm = _amoebaReferenceGeneralizedKirkwoodForce->getIncludeCavityTerm();
+    _probeRadius = _amoebaReferenceGeneralizedKirkwoodForce->getProbeRadius();
+    _surfaceAreaFactor = _amoebaReferenceGeneralizedKirkwoodForce->getSurfaceAreaFactor();
+    _dielectricOffset = _amoebaReferenceGeneralizedKirkwoodForce->getDielectricOffset();
+    _tanhRescaling = _amoebaReferenceGeneralizedKirkwoodForce->getTanhRescaling();
+    _amoebaReferenceGeneralizedKirkwoodForce->getTanhParameters(_beta0, _beta1, _beta2);
+    _descreenOffset = _amoebaReferenceGeneralizedKirkwoodForce->getDescreenOffset();
 }
 
 AmoebaReferenceGeneralizedKirkwoodMultipoleForce::~AmoebaReferenceGeneralizedKirkwoodMultipoleForce()
@@ -2165,6 +2168,12 @@ double AmoebaReferenceGeneralizedKirkwoodMultipoleForce::getDielectricOffset() c
      return _dielectricOffset;
 };
 
+double AmoebaReferenceGeneralizedKirkwoodMultipoleForce::getDescreenOffset() const
+{
+    return _descreenOffset;
+};
+
+
 void AmoebaReferenceGeneralizedKirkwoodMultipoleForce::zeroFixedMultipoleFields()
 {
     this->AmoebaReferenceMultipoleForce::zeroFixedMultipoleFields();
@@ -4060,67 +4069,125 @@ double AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateElectrostatic(
     return energy;
 }
 
-void AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateGrycukChainRulePairIxn(const MultipoleParticleData& particleI, const MultipoleParticleData& particleJ,
-                                                                                       const vector<double>& dBorn, vector<Vec3>& forces) const
-{
-    unsigned int iIndex             = particleI.particleIndex;
-    unsigned int jIndex             = particleJ.particleIndex;
+static double neckDescreenDerivative(double r, double radius, double radiusK, double sneck) {
+    // Radius of water in nm.
+    double radiusWater = 0.14;
+
+    if (r > radius + radiusK + 2 * radiusWater) {
+        return 0.0;
+    }
 
-    double pi43                     = (4.0/3.0)*M_PI;
+    // Get Aij and Bij
+    double Aij, Bij;
+    AmoebaGeneralizedKirkwoodForceImpl::getNeckConstants(radius, radiusK, Aij, Bij);
+    // Use Aij and Bij to get neck value using derivative of Equation 13 from Aguilar/Onufriev 2010 paper
+    double rMinusBij = r - Bij;
+    double rMinusBij3 = rMinusBij * rMinusBij * rMinusBij;
+    double rMinusBij4 = rMinusBij3 * rMinusBij;
+    double radiiMinusr = radius + radiusK + 2.0 * radiusWater - r;
+    double radiiMinusr3 = radiiMinusr * radiiMinusr * radiiMinusr;
+    double radiiMinusr4 = radiiMinusr3 * radiiMinusr;
+
+    double PI4_3 = 4.0 * M_PI / 3.0;
+    return 4.0 * PI4_3 * (sneck * Aij * rMinusBij3 * radiiMinusr4 - sneck * Aij * rMinusBij4 * radiiMinusr3);
+}
+
+void AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateGrycukChainRulePairIxn(
+        const MultipoleParticleData &particleI, const MultipoleParticleData &particleJ,
+        const vector<double> &dBorn, vector<Vec3> &forces) const {
+    unsigned int iIndex = particleI.particleIndex;
+    unsigned int jIndex = particleJ.particleIndex;
+
+    if (iIndex == jIndex) return;
+
+    double bornRadiusI = _bornRadii[iIndex];
+
+    const double bigRadius = 3.0;
+    if (bornRadiusI >= bigRadius) return;
+
+    double third = 1.0 / 3.0;
+    double pi43 = 4.0 * third * M_PI;
 
     double lik, uik;
     double lik4, uik4;
-    double factor                   = -pow(M_PI, (1.0/3.0))*pow(6.0,(2.0/3.0))/9.0;
-    double term                     = pi43/(_bornRadii[iIndex]*_bornRadii[iIndex]*_bornRadii[iIndex]);
-          term                      = factor/pow(term, (4.0/3.0));
-
-    Vec3 deltaR                     = particleJ.position - particleI.position;
+    double factor = -pow(M_PI, third) * pow(6.0, 2.0 * third) / 9.0;
+    double term = pi43 / (_bornRadii[iIndex] * _bornRadii[iIndex] * _bornRadii[iIndex]);
+    term = factor / pow(term, 4.0 * third);
+    if (_tanhRescaling) {
+        double rhoi = _atomicRadii[iIndex];
+        double rhoi3 = rhoi * rhoi * rhoi;
+        double rhoi3Psi = rhoi3 * _soluteIntegral[iIndex];
+        double rhoi6Psi2 = rhoi3Psi * rhoi3Psi;
+        double rhoi9Psi3 = rhoi6Psi2 * rhoi3Psi;
+        double rhoi6Psi = rhoi3 * rhoi3 * _soluteIntegral[iIndex];
+        double rhoi9Psi2 = rhoi6Psi2 * rhoi3;
+        double tanhTerm = tanh(_beta0 * rhoi3Psi - _beta1 * rhoi6Psi2 + _beta2 * rhoi9Psi3);
+        double tanh2 = tanhTerm * tanhTerm;
+        double chainRuleTerm = _beta0 * rhoi3 - 2.0 * _beta1 * rhoi6Psi + 3.0 * _beta2 * rhoi9Psi2;
+        double recipBigRad = 1.0 / bigRadius;
+        double recipBigRad3 = recipBigRad * recipBigRad * recipBigRad;
+        double tanh_constant = pi43 * ((1.0 / rhoi3) - recipBigRad3);
+        term = term * tanh_constant * chainRuleTerm * (1.0 - tanh2);
+    }
 
-    double sk                       = _scaledRadii[jIndex];
-    double sk2                      = sk*sk;
-    double r2                       = deltaR.dot(deltaR);
-    double r                        = sqrt(r2);
-    double de                       = 0.0;
+    Vec3 deltaR = particleJ.position - particleI.position;
+    double sk = _scaleFactors[jIndex] * _descreenRadii[jIndex];
+    if (sk <= 0.0) return;
+    double sk2 = sk * sk;
+    double r2 = deltaR.dot(deltaR);
+    double r = sqrt(r2);
+    double de = 0.0;
 
-    // If atom index engulfs the descreening atom, then there is no descreening. 
-    if (_atomicRadii[iIndex] > r + sk) return;
+    double baseRadiusI = max(_atomicRadii[iIndex], _descreenRadii[iIndex]) + _descreenOffset;
+    // If atom index engulfs the descreening atom, then there is no descreening.
+    if (baseRadiusI > r + sk) return;
 
-    if ((_atomicRadii[iIndex] + r) < sk) {
+    if ((baseRadiusI + r) < sk) {
         double uik4;
-        uik                = sk - r;
-        uik4               = uik*uik;
-        uik4               = uik4*uik4;
-        de                 = -4.0*M_PI/uik4;
+        uik = sk - r;
+        uik4 = uik * uik;
+        uik4 = uik4 * uik4;
+        de = -4.0 * M_PI / uik4;
     }
 
-    if ((_atomicRadii[iIndex] + r) < sk) {
-        lik          = sk - r;
-        lik4         = lik*lik;
-        lik4         = lik4*lik4;
-        de          += 0.25*M_PI*(sk2-4.0*sk*r+17.0*r2)/ (r2*lik4);
-    } else if (r < (_atomicRadii[iIndex] +sk)) {
-        lik          = _atomicRadii[iIndex];
-        lik4         = lik*lik;
-        lik4         = lik4*lik4;
-        de          += 0.25*M_PI*(2.0*_atomicRadii[iIndex]*_atomicRadii[iIndex]-sk2-r2)/ (r2*lik4);
+    if ((baseRadiusI + r) < sk) {
+        lik = sk - r;
+        lik4 = lik * lik;
+        lik4 = lik4 * lik4;
+        de += 0.25 * M_PI * (sk2 - 4.0 * sk * r + 17.0 * r2) / (r2 * lik4);
+    } else if (r < (baseRadiusI + sk)) {
+        lik = baseRadiusI;
+        lik4 = lik * lik;
+        lik4 = lik4 * lik4;
+        de += 0.25 * M_PI * (2.0 * baseRadiusI * baseRadiusI - sk2 - r2) / (r2 * lik4);
     } else {
-        lik          = r - sk;
-        lik4         = lik*lik;
-        lik4         = lik4*lik4;
-        de          += 0.25*M_PI*(sk2-4.0*sk*r+r2)/ (r2*lik4);
+        lik = r - sk;
+        lik4 = lik * lik;
+        lik4 = lik4 * lik4;
+        de += 0.25 * M_PI * (sk2 - 4.0 * sk * r + r2) / (r2 * lik4);
+    }
+    uik = r + sk;
+    uik4 = uik * uik;
+    uik4 = uik4 * uik4;
+
+    de -= 0.25 * M_PI * (sk2 + 4.0 * sk * r + r2) / (r2 * uik4);
+
+    double mixedNeckScale =  0.5 * (_neckFactors[iIndex] + _neckFactors[jIndex]);
+    if (mixedNeckScale > 0.0) {
+        de += neckDescreenDerivative(r, baseRadiusI, _descreenRadii[jIndex], mixedNeckScale);
     }
-    uik                = r + sk;
-    uik4               = uik*uik;
-    uik4               = uik4*uik4;
 
-    de                -= 0.25*M_PI*(sk2+4.0*sk*r+r2)/ (r2*uik4);
-    double dbr         = term*de/r;
-          de           = dbr*dBorn[iIndex];
+    double dbr = term * de / r;
+    de = dbr * dBorn[iIndex];
+
+    forces[iIndex] -= deltaR * de;
+    forces[jIndex] += deltaR * de;
 
-    forces[iIndex]    -= deltaR*de;
-    forces[jIndex]    += deltaR*de;
 }
 
+
+
+
 double AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateCavityTermEnergyAndForces(vector<double>& dBorn) const
 {
 

plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceMultipoleForce.h

@@ -1239,10 +1239,16 @@ public:
      * Get dielectric offset.
      *
      * @return dielectric offset
-     *
      */
     double getDielectricOffset() const;
 
+    /**
+     * Get the descreen offset.
+     *
+     * @return descreen offset.
+     */
+    double getDescreenOffset() const;
+
 private:
 
     AmoebaReferenceGeneralizedKirkwoodForce* _amoebaReferenceGeneralizedKirkwoodForce;
@@ -1252,10 +1258,17 @@ private:
     double _fd;
     double _fq;
 
+    double _beta0;
+    double _beta1;
+    double _beta2;
+
     std::vector<double> _atomicRadii;
-    std::vector<double> _scaledRadii;
+    std::vector<double> _scaleFactors;
+    std::vector<double> _descreenRadii;
+    std::vector<double> _neckFactors;
+
     std::vector<double> _bornRadii;
-    std::vector<double> _bornForce;
+    std::vector<double> _soluteIntegral;
 
     std::vector<Vec3> _gkField;
     std::vector<Vec3> _inducedDipoleS;
@@ -1269,6 +1282,8 @@ private:
     double _probeRadius;
     double _surfaceAreaFactor;
     double _dielectricOffset;
+    double _tanhRescaling;
+    double _descreenOffset;
 
     /**
      * Zero fixed multipole fields.