View Javadoc
1   package ffx.potential.constraint;
2   
3   import ffx.numerics.Constraint;
4   
5   public class ShakeChargeConstraint implements Constraint {
6       final int nConstraints;
7       final double tol;
8       final int c;
9       final double maxIters = 150;
10  
11      public ShakeChargeConstraint(int nConstraints, int c, double tol){
12          this.nConstraints = nConstraints;
13          this.c = c;
14          this.tol = tol;
15      }
16      @Override
17      public void applyConstraintToStep(final double[] xPrior, double[] xNew, final double[] masses, double tol){
18      }
19  
20      /**
21       * This method follows the SHAKE Charge Constraint laid out in Appendix B of
22       * Donnini, Serena, et al. "Charge-neutral constant pH molecular dynamics simulations using a parsimonious proton buffer."
23       * Journal of chemical theory and computation 12.3 (2016): 1040-1051.
24       * https://pubs.acs.org/doi/epdf/10.1021/acs.jctc.5b01160
25       * @param x
26       * @param afric
27       * @param masses
28       * @param dt
29       */
30      public void applyChargeConstraintToStep(final double[] x, final double[] afric, final double[] masses, final double dt){
31          boolean done = false;
32          int iter = 0;
33          while(!done){
34              done = true;
35              iter++;
36              double totalLambda = 0.0;
37              double totalInverseMass = 0.0;
38              for (int i = 0; i < nConstraints; i++) {
39                  double lambda = Math.sin(x[i]) * Math.sin(x[i]);
40                  totalLambda += lambda;
41                  totalInverseMass += (1.0 / masses[i]);
42              }
43              double delta = totalLambda - c;
44  
45              if (Math.abs(delta) > tol) {
46                  done = false;
47                  double term = delta / (dt * dt * totalInverseMass);
48                  for (int i = 0; i < nConstraints; i++) {
49                      double g = -term * Math.sin(2 * x[i]);
50                      x[i] = x[i] + g * afric[i];
51                  }
52              }
53              if(iter >= maxIters){
54                  throw new RuntimeException("SHAKE Charge Constraint -- Warning, Charge Constraint not Satisfied");
55              }
56          }
57      }
58  
59      @Override
60      public void applyConstraintToVelocities(final double[] x, double[] v, final double[] masses, double tol){
61  
62      }
63  
64      @Override
65      public int[] constrainedAtomIndices() {
66          return new int[0];
67      }
68  
69      @Override
70      public boolean constraintSatisfied(double[] x, double tol) {
71          return false;
72      }
73  
74      @Override
75      public boolean constraintSatisfied(double[] x, double[] v, double xTol, double vTol) {
76          return false;
77      }
78  
79      @Override
80      public int getNumDegreesFrozen() {
81          return 0;
82      }
83  }
84  
85