Additional pass at cleaning up SHAKE Charge Constraint algorithm. The changes include accurately calculating a constraint force and scaling it down by a factor of 100 as determined through testing. Removal of the afric term opens up the use of the constraint with other integrators.

Author: acthiel212

modules/algorithms/src/main/java/ffx/algorithms/dynamics/integrators/Stochastic.java

@@ -230,15 +230,11 @@ public void preForce(Potential potential) throws RuntimeException {
       // then find new atom positions and half-step velocities via Verlet recursion.
       x[i] += (v[i] * vFriction[i] + a[i] * afric + prand);
       v[i] = v[i] * pfric + 0.5 * a[i] * vFriction[i];
-
-      if(useConstraints){
-        africArray[i] = afric;
-      }
     }
     if (useConstraints) {
       for(Constraint c : constraints){
         if(c instanceof ShakeChargeConstraint){
-          ((ShakeChargeConstraint) c).applyChargeConstraintToStep(x, africArray, mass, dt);
+          ((ShakeChargeConstraint) c).applyChargeConstraintToStep(x, mass);
           double velScale = 1.0 / dt;
           for (int i = 0; i < nVariables; i++) {
             v[i] = velScale * (x[i] - xPrior[i]);

modules/potential/src/main/java/ffx/potential/constraint/ShakeChargeConstraint.java

@@ -2,6 +2,8 @@
 
 import ffx.numerics.Constraint;
 
+import static ffx.utilities.Constants.KCAL_TO_GRAM_ANG2_PER_PS2;
+
 public class ShakeChargeConstraint implements Constraint {
     final int nConstraints;
     final double tol;
@@ -23,11 +25,9 @@ public void applyConstraintToStep(final double[] xPrior, double[] xNew, final do
      * Journal of chemical theory and computation 12.3 (2016): 1040-1051.
      * https://pubs.acs.org/doi/epdf/10.1021/acs.jctc.5b01160
      * @param x
-     * @param afric
      * @param masses
-     * @param dt
      */
-    public void applyChargeConstraintToStep(final double[] x, final double[] afric, final double[] masses, final double dt){
+    public void applyChargeConstraintToStep(final double[] x, final double[] masses){
         boolean done = false;
         int iter = 0;
         while(!done){
@@ -41,13 +41,13 @@ public void applyChargeConstraintToStep(final double[] x, final double[] afric,
                 totalInverseMass += (1.0 / masses[i]);
             }
             double delta = totalLambda - c;
-
             if (Math.abs(delta) > tol) {
                 done = false;
-                double term = delta / (dt * dt * totalInverseMass);
+                double term = delta / totalInverseMass;
                 for (int i = 0; i < nConstraints; i++) {
                     double g = -term * Math.sin(2 * x[i]);
-                    x[i] = x[i] + g * afric[i];
+                    // Scale constraint force by 0.01 was determined to improve convergence of constraint through testing
+                    x[i] = x[i] - (g * -KCAL_TO_GRAM_ANG2_PER_PS2 / masses[i]) * 0.01;
                 }
             }
             if(iter >= maxIters){