1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38 package ffx.potential.openmm;
39
40 import ffx.openmm.CustomIntegrator;
41
42 import java.util.logging.Logger;
43
44 import static ffx.utilities.Constants.KCAL_TO_KJ;
45 import static ffx.utilities.Constants.R;
46 import static java.lang.Math.exp;
47 import static java.lang.String.format;
48 import static java.util.Arrays.copyOfRange;
49 import static org.apache.commons.math3.util.FastMath.sqrt;
50
51
52
53
54 public class CustomMTSLangevinIntegrator extends CustomIntegrator {
55
56 private static final Logger logger = Logger.getLogger(CustomMTSLangevinIntegrator.class.getName());
57
58
59
60
61
62
63
64
65
66 public CustomMTSLangevinIntegrator(double dt, double temperature,
67 double frictionCoeff, boolean hasAmoebaCavitationForce) {
68 super(dt);
69
70 int n = 4;
71
72
73 int[] forceGroups = {1, 0};
74
75 int[] subSteps = {1, 4};
76 if (hasAmoebaCavitationForce) {
77 n = 8;
78
79
80
81 forceGroups = new int[]{2, 1, 0};
82
83
84
85 subSteps = new int[]{1, 2, 8};
86 }
87
88 addGlobalVariable("a", exp(-frictionCoeff * dt / n));
89 addGlobalVariable("b", sqrt(1.0 - exp(-2.0 * frictionCoeff * dt / n)));
90 addGlobalVariable("kT", R * temperature * KCAL_TO_KJ);
91 addPerDofVariable("x1", 0.0);
92 StringBuilder sb = new StringBuilder(" Update Context State\n");
93 addUpdateContextState();
94 createMTSLangevinSubStep(1, forceGroups, subSteps, sb);
95
96 logger.finest(" Langevin-MTS steps:" + sb);
97 addConstrainVelocities();
98 logger.info(" Custom MTS Langevin Integrator");
99 logger.info(format(" Time step: %6.2f (fsec)", dt * 1000));
100 logger.info(format(" Inner Time step: %6.2f (fsec)", dt / n * 1000));
101 logger.info(format(" Friction Coefficient: %6.2f (1/psec)", frictionCoeff));
102 }
103
104
105
106
107
108
109
110
111 public void createMTSLangevinSubStep(int parentSubsteps, int[] forceGroups, int[] subSteps,
112 StringBuilder sb) {
113 int forceGroup = forceGroups[0];
114 int steps = subSteps[0];
115 int stepsPerParentStep = steps / parentSubsteps;
116 if (stepsPerParentStep < 1 || steps % parentSubsteps != 0) {
117 throw new IllegalArgumentException(
118 "The number for substeps for each group must be a multiple of the number for the previous group");
119 }
120 if (forceGroup < 0 || forceGroup > 31) {
121 throw new IllegalArgumentException("Force group must be between 0 and 31");
122 }
123
124 sb.append(" Force Group: ").append(forceGroup).append(" ForceGroup length: ")
125 .append(forceGroups.length).append(" Steps: ").append(steps)
126 .append(" Step Per Parent Step: ").append(stepsPerParentStep).append(" Parent Sub Steps: ")
127 .append(parentSubsteps).append("\n");
128
129 for (int i = 0; i < stepsPerParentStep; i++) {
130 String step = "v+0.5*(dt/" + steps + ")*f" + forceGroup + "/m";
131 sb.append(" v = ").append(step).append("\n");
132 addComputePerDof("v", step);
133
134 if (forceGroups.length == 1) {
135 step = "x+(dt/" + (2 * steps) + ")*v";
136 sb.append(" x = ").append(step).append("\n");
137 addComputePerDof("x", step);
138 step = "a*v + b*sqrt(kT/m)*gaussian";
139 sb.append(" v = ").append(step).append("\n");
140 addComputePerDof("v", step);
141 step = "x+(dt/" + (2 * steps) + ")*v";
142 sb.append(" x = ").append(step).append("\n");
143 addComputePerDof("x", step);
144 step = "x";
145 sb.append(" x1 = ").append(step).append("\n");
146 addComputePerDof("x1", step);
147 sb.append(" Constrain Positions\n");
148 addConstrainPositions();
149 step = "v+(x-x1)/(dt/" + steps + ")";
150 sb.append(" v = ").append(step).append("\n");
151 addComputePerDof("v", step);
152 sb.append(" Constrain Velocities\n");
153 addConstrainVelocities();
154 } else {
155 createMTSLangevinSubStep(steps, copyOfRange(forceGroups, 1, forceGroups.length),
156 copyOfRange(subSteps, 1, subSteps.length), sb);
157 }
158 step = "v+0.5*(dt/" + steps + ")*f" + forceGroup + "/m";
159 sb.append(" v = ").append(step).append("\n");
160 addComputePerDof("v", step);
161 }
162 }
163 }