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