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.CustomCompoundBondForce;
41 import ffx.openmm.DoubleArray;
42 import ffx.openmm.Force;
43 import ffx.openmm.IntArray;
44 import ffx.potential.ForceFieldEnergy;
45 import ffx.potential.bonded.StretchBend;
46 import ffx.potential.terms.StretchBendPotentialEnergy;
47
48 import java.util.logging.Logger;
49
50 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_KJPerKcal;
51 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_NmPerAngstrom;
52 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_RadiansPerDegree;
53 import static java.lang.String.format;
54
55
56
57
58 public class StretchBendForce extends CustomCompoundBondForce {
59
60 private static final Logger logger = Logger.getLogger(StretchBendForce.class.getName());
61
62
63
64
65
66
67 public StretchBendForce(StretchBendPotentialEnergy stretchBendPotentialEnergy) {
68 super(3, StretchBendPotentialEnergy.getStretchBendEnergyString());
69 StretchBend[] stretchBends = stretchBendPotentialEnergy.getStretchBendArray();
70 addPerBondParameter("r12");
71 addPerBondParameter("r23");
72 addPerBondParameter("theta0");
73 addPerBondParameter("k1");
74 addPerBondParameter("k2");
75 setName("AmoebaStretchBend");
76
77 IntArray particles = new IntArray(0);
78 DoubleArray parameters = new DoubleArray(0);
79 for (StretchBend stretchBend : stretchBends) {
80 int i1 = stretchBend.getAtom(0).getArrayIndex();
81 int i2 = stretchBend.getAtom(1).getArrayIndex();
82 int i3 = stretchBend.getAtom(2).getArrayIndex();
83 double r12 = stretchBend.bond0Eq * OpenMM_NmPerAngstrom;
84 double r23 = stretchBend.bond1Eq * OpenMM_NmPerAngstrom;
85 double theta0 = stretchBend.angleEq * OpenMM_RadiansPerDegree;
86 double k1 = stretchBend.force0 * OpenMM_KJPerKcal / OpenMM_NmPerAngstrom;
87 double k2 = stretchBend.force1 * OpenMM_KJPerKcal / OpenMM_NmPerAngstrom;
88 particles.append(i1);
89 particles.append(i2);
90 particles.append(i3);
91 parameters.append(r12);
92 parameters.append(r23);
93 parameters.append(theta0);
94 parameters.append(k1);
95 parameters.append(k2);
96 addBond(particles, parameters);
97 particles.resize(0);
98 parameters.resize(0);
99 }
100 particles.destroy();
101 parameters.destroy();
102
103 int forceGroup = stretchBendPotentialEnergy.getForceGroup();
104 setForceGroup(forceGroup);
105 logger.info(format(" Stretch-Bends: %10d", stretchBends.length));
106 logger.fine(format(" Force Group: %10d", forceGroup));
107 }
108
109
110
111
112
113
114
115
116 public StretchBendForce(StretchBendPotentialEnergy stretchBendPotentialEnergy,
117 int topology, OpenMMDualTopologyEnergy openMMDualTopologyEnergy) {
118 super(3, StretchBendPotentialEnergy.getStretchBendEnergyString());
119 StretchBend[] stretchBends = stretchBendPotentialEnergy.getStretchBendArray();
120 addPerBondParameter("r12");
121 addPerBondParameter("r23");
122 addPerBondParameter("theta0");
123 addPerBondParameter("k1");
124 addPerBondParameter("k2");
125 setName("AmoebaStretchBend");
126
127 double scale = openMMDualTopologyEnergy.getTopologyScale(topology);
128
129 IntArray particles = new IntArray(0);
130 DoubleArray parameters = new DoubleArray(0);
131 for (StretchBend stretchBend : stretchBends) {
132 int i1 = stretchBend.getAtom(0).getArrayIndex();
133 int i2 = stretchBend.getAtom(1).getArrayIndex();
134 int i3 = stretchBend.getAtom(2).getArrayIndex();
135 double r12 = stretchBend.bond0Eq * OpenMM_NmPerAngstrom;
136 double r23 = stretchBend.bond1Eq * OpenMM_NmPerAngstrom;
137 double theta0 = stretchBend.angleEq * OpenMM_RadiansPerDegree;
138 double k1 = stretchBend.force0 * OpenMM_KJPerKcal / OpenMM_NmPerAngstrom;
139 double k2 = stretchBend.force1 * OpenMM_KJPerKcal / OpenMM_NmPerAngstrom;
140
141 if (!stretchBend.applyLambda()) {
142 k1 = k1 * scale;
143 k2 = k2 * scale;
144 }
145 i1 = openMMDualTopologyEnergy.mapToDualTopologyIndex(topology, i1);
146 i2 = openMMDualTopologyEnergy.mapToDualTopologyIndex(topology, i2);
147 i3 = openMMDualTopologyEnergy.mapToDualTopologyIndex(topology, i3);
148 particles.append(i1);
149 particles.append(i2);
150 particles.append(i3);
151 parameters.append(r12);
152 parameters.append(r23);
153 parameters.append(theta0);
154 parameters.append(k1);
155 parameters.append(k2);
156 addBond(particles, parameters);
157 particles.resize(0);
158 parameters.resize(0);
159 }
160 particles.destroy();
161 parameters.destroy();
162
163 int forceGroup = stretchBendPotentialEnergy.getForceGroup();
164 setForceGroup(forceGroup);
165 logger.info(format(" Stretch-Bends: %10d", stretchBends.length));
166 logger.fine(format(" Force Group: %10d", forceGroup));
167 }
168
169
170
171
172
173
174
175 public static Force constructForce(OpenMMEnergy openMMEnergy) {
176 StretchBendPotentialEnergy stretchBendPotentialEnergy = openMMEnergy.getStretchBendPotentialEnergy();
177 if (stretchBendPotentialEnergy == null) {
178 return null;
179 }
180 return new StretchBendForce(stretchBendPotentialEnergy);
181 }
182
183
184
185
186
187
188
189
190 public static Force constructForce(int topology, OpenMMDualTopologyEnergy openMMDualTopologyEnergy) {
191 ForceFieldEnergy forceFieldEnergy = openMMDualTopologyEnergy.getForceFieldEnergy(topology);
192 StretchBendPotentialEnergy stretchBendPotentialEnergy = forceFieldEnergy.getStretchBendPotentialEnergy();
193 if (stretchBendPotentialEnergy == null) {
194 return null;
195 }
196 return new StretchBendForce(stretchBendPotentialEnergy, topology, openMMDualTopologyEnergy);
197 }
198
199
200
201
202
203
204 public void updateForce(OpenMMEnergy openMMEnergy) {
205 StretchBendPotentialEnergy stretchBendPotentialEnergy = openMMEnergy.getStretchBendPotentialEnergy();
206 if (stretchBendPotentialEnergy == null) {
207 return;
208 }
209 StretchBend[] stretchBends = stretchBendPotentialEnergy.getStretchBendArray();
210
211 IntArray particles = new IntArray(0);
212 DoubleArray parameters = new DoubleArray(0);
213 int index = 0;
214 for (StretchBend stretchBend : stretchBends) {
215 int i1 = stretchBend.getAtom(0).getArrayIndex();
216 int i2 = stretchBend.getAtom(1).getArrayIndex();
217 int i3 = stretchBend.getAtom(2).getArrayIndex();
218 double r12 = stretchBend.bond0Eq * OpenMM_NmPerAngstrom;
219 double r23 = stretchBend.bond1Eq * OpenMM_NmPerAngstrom;
220 double theta0 = stretchBend.angleEq * OpenMM_RadiansPerDegree;
221 double k1 = stretchBend.force0 * OpenMM_KJPerKcal / OpenMM_NmPerAngstrom;
222 double k2 = stretchBend.force1 * OpenMM_KJPerKcal / OpenMM_NmPerAngstrom;
223 particles.append(i1);
224 particles.append(i2);
225 particles.append(i3);
226 parameters.append(r12);
227 parameters.append(r23);
228 parameters.append(theta0);
229 parameters.append(k1);
230 parameters.append(k2);
231 setBondParameters(index++, particles, parameters);
232 particles.resize(0);
233 parameters.resize(0);
234 }
235 particles.destroy();
236 parameters.destroy();
237
238 updateParametersInContext(openMMEnergy.getContext());
239 }
240
241
242
243
244
245
246
247 public void updateForce(int topology, OpenMMDualTopologyEnergy openMMDualTopologyEnergy) {
248 ForceFieldEnergy forceFieldEnergy = openMMDualTopologyEnergy.getForceFieldEnergy(topology);
249 StretchBendPotentialEnergy stretchBendPotentialEnergy = forceFieldEnergy.getStretchBendPotentialEnergy();
250 if (stretchBendPotentialEnergy == null) {
251 return;
252 }
253 StretchBend[] stretchBends = stretchBendPotentialEnergy.getStretchBendArray();
254
255 double scale = openMMDualTopologyEnergy.getTopologyScale(topology);
256
257 IntArray particles = new IntArray(0);
258 DoubleArray parameters = new DoubleArray(0);
259 int index = 0;
260 for (StretchBend stretchBend : stretchBends) {
261 int i1 = stretchBend.getAtom(0).getArrayIndex();
262 int i2 = stretchBend.getAtom(1).getArrayIndex();
263 int i3 = stretchBend.getAtom(2).getArrayIndex();
264 double r12 = stretchBend.bond0Eq * OpenMM_NmPerAngstrom;
265 double r23 = stretchBend.bond1Eq * OpenMM_NmPerAngstrom;
266 double theta0 = stretchBend.angleEq * OpenMM_RadiansPerDegree;
267 double k1 = stretchBend.force0 * OpenMM_KJPerKcal / OpenMM_NmPerAngstrom;
268 double k2 = stretchBend.force1 * OpenMM_KJPerKcal / OpenMM_NmPerAngstrom;
269
270 if (!stretchBend.applyLambda()) {
271 k1 = k1 * scale;
272 k2 = k2 * scale;
273 }
274 i1 = openMMDualTopologyEnergy.mapToDualTopologyIndex(topology, i1);
275 i2 = openMMDualTopologyEnergy.mapToDualTopologyIndex(topology, i2);
276 i3 = openMMDualTopologyEnergy.mapToDualTopologyIndex(topology, i3);
277 particles.append(i1);
278 particles.append(i2);
279 particles.append(i3);
280 parameters.append(r12);
281 parameters.append(r23);
282 parameters.append(theta0);
283 parameters.append(k1);
284 parameters.append(k2);
285 setBondParameters(index++, particles, parameters);
286 particles.resize(0);
287 parameters.resize(0);
288 }
289 particles.destroy();
290 parameters.destroy();
291
292 updateParametersInContext(openMMDualTopologyEnergy.getContext());
293 }
294 }