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.CustomAngleForce;
41 import ffx.openmm.DoubleArray;
42 import ffx.openmm.Force;
43 import ffx.potential.ForceFieldEnergy;
44 import ffx.potential.bonded.Angle;
45 import ffx.potential.bonded.Atom;
46 import ffx.potential.parameters.AngleType;
47 import ffx.potential.parameters.ForceField;
48 import ffx.potential.terms.AnglePotentialEnergy;
49
50 import java.util.logging.Level;
51 import java.util.logging.Logger;
52
53 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_KJPerKcal;
54 import static java.lang.String.format;
55
56
57
58
59 public class AngleForce extends CustomAngleForce {
60
61 private static final Logger logger = Logger.getLogger(AngleForce.class.getName());
62
63 private int nAngles = 0;
64 private final boolean manyBodyTitration;
65 private final boolean rigidHydrogenAngles;
66
67
68
69
70
71
72 public AngleForce(AnglePotentialEnergy anglePotentialEnergy, OpenMMEnergy openMMEnergy) {
73 super(anglePotentialEnergy.getAngleEnergyString());
74 ForceField forceField = openMMEnergy.getMolecularAssembly().getForceField();
75 manyBodyTitration = forceField.getBoolean("MANYBODY_TITRATION", false);
76 rigidHydrogenAngles = forceField.getBoolean("RIGID_HYDROGEN_ANGLES", false);
77 addPerAngleParameter("theta0");
78 addPerAngleParameter("k");
79 setName("Angle");
80
81 DoubleArray parameters = new DoubleArray(0);
82 Angle[] angles = anglePotentialEnergy.getAngleArray();
83 for (Angle angle : angles) {
84 AngleType angleType = angle.getAngleType();
85 AngleType.AngleMode angleMode = angleType.angleMode;
86 if (!manyBodyTitration && angleMode == AngleType.AngleMode.IN_PLANE) {
87
88 } else if (isHydrogenAngle(angle) && rigidHydrogenAngles) {
89 logger.log(Level.INFO, " Constrained angle %s was not added the AngleForce.", angle);
90 } else {
91 int i1 = angle.getAtom(0).getArrayIndex();
92 int i2 = angle.getAtom(1).getArrayIndex();
93 int i3 = angle.getAtom(2).getArrayIndex();
94
95 double theta0 = angleType.angle[angle.nh];
96 double k = OpenMM_KJPerKcal * angleType.angleUnit * angleType.forceConstant;
97 if (angleMode == AngleType.AngleMode.IN_PLANE) {
98
99
100 k = 0.0;
101 }
102 parameters.append(theta0);
103 parameters.append(k);
104 addAngle(i1, i2, i3, parameters);
105 nAngles++;
106 parameters.resize(0);
107 }
108 }
109 parameters.destroy();
110
111 if (nAngles > 0) {
112 int forceGroup = anglePotentialEnergy.getForceGroup();
113 setForceGroup(forceGroup);
114 logger.info(format(" Angles: %10d", nAngles));
115 logger.fine(format(" Force Group: %10d", forceGroup));
116 }
117 }
118
119
120
121
122
123
124
125
126 public AngleForce(AnglePotentialEnergy anglePotentialEnergy,
127 int topology, OpenMMDualTopologyEnergy openMMDualTopologyEnergy) {
128 super(anglePotentialEnergy.getAngleEnergyString());
129 Angle[] angles = anglePotentialEnergy.getAngleArray();
130 addPerAngleParameter("theta0");
131 addPerAngleParameter("k");
132 setName("Angle");
133
134 ForceFieldEnergy forceFieldEnergy = openMMDualTopologyEnergy.getForceFieldEnergy(topology);
135 ForceField forceField = forceFieldEnergy.getMolecularAssembly().getForceField();
136 manyBodyTitration = forceField.getBoolean("MANYBODY_TITRATION", false);
137 rigidHydrogenAngles = forceField.getBoolean("RIGID_HYDROGEN_ANGLES", false);
138 if (manyBodyTitration || rigidHydrogenAngles) {
139 logger.severe("Dual Topology does not support rigid hydrogen angles or many body titration.");
140 }
141
142 double scale = openMMDualTopologyEnergy.getTopologyScale(topology);
143
144 DoubleArray parameters = new DoubleArray(0);
145 for (Angle angle : angles) {
146 AngleType angleType = angle.angleType;
147 AngleType.AngleMode angleMode = angleType.angleMode;
148 if (angleMode == AngleType.AngleMode.IN_PLANE) {
149
150 } else {
151 int i1 = angle.getAtom(0).getArrayIndex();
152 int i2 = angle.getAtom(1).getArrayIndex();
153 int i3 = angle.getAtom(2).getArrayIndex();
154
155 double theta0 = angleType.angle[angle.nh];
156 double k = OpenMM_KJPerKcal * angleType.angleUnit * angleType.forceConstant;
157
158 if (!angle.applyLambda()) {
159 k = k * scale;
160 }
161 parameters.append(theta0);
162 parameters.append(k);
163 i1 = openMMDualTopologyEnergy.mapToDualTopologyIndex(topology, i1);
164 i2 = openMMDualTopologyEnergy.mapToDualTopologyIndex(topology, i2);
165 i3 = openMMDualTopologyEnergy.mapToDualTopologyIndex(topology, i3);
166 addAngle(i1, i2, i3, parameters);
167 nAngles++;
168 parameters.resize(0);
169 }
170 }
171 parameters.destroy();
172
173 if (nAngles > 0) {
174 int forceGroup = anglePotentialEnergy.getForceGroup();
175 setForceGroup(forceGroup);
176 logger.info(format(" Angles: %10d", nAngles));
177 logger.fine(format(" Force Group: %10d", forceGroup));
178 }
179 }
180
181
182
183
184
185
186
187 public static Force constructForce(OpenMMEnergy openMMEnergy) {
188 AnglePotentialEnergy anglePotentialEnergy = openMMEnergy.getAnglePotentialEnergy();
189 if (anglePotentialEnergy == null) {
190 return null;
191 }
192 AngleForce angleForce = new AngleForce(anglePotentialEnergy, openMMEnergy);
193 if (angleForce.nAngles > 0) {
194 return angleForce;
195 }
196 angleForce.destroy();
197 return null;
198 }
199
200
201
202
203
204
205
206 public static Force constructForce(int topology, OpenMMDualTopologyEnergy openMMDualTopologyEnergy) {
207 ForceFieldEnergy forceFieldEnergy = openMMDualTopologyEnergy.getForceFieldEnergy(topology);
208 AnglePotentialEnergy anglePotentialEnergy = forceFieldEnergy.getAnglePotentialEnergy();
209 if (anglePotentialEnergy == null) {
210 return null;
211 }
212 AngleForce angleForce = new AngleForce(anglePotentialEnergy, topology, openMMDualTopologyEnergy);
213 if (angleForce.nAngles > 0) {
214 return angleForce;
215 }
216 angleForce.destroy();
217 return null;
218 }
219
220
221
222
223
224
225 public void updateForce(OpenMMEnergy openMMEnergy) {
226 AnglePotentialEnergy anglePotentialEnergy = openMMEnergy.getAnglePotentialEnergy();
227 if (anglePotentialEnergy == null) {
228 return;
229 }
230 Angle[] angles = anglePotentialEnergy.getAngleArray();
231
232 DoubleArray parameters = new DoubleArray(0);
233 int index = 0;
234 for (Angle angle : angles) {
235 AngleType.AngleMode angleMode = angle.angleType.angleMode;
236 if (!manyBodyTitration && angleMode == AngleType.AngleMode.IN_PLANE) {
237
238 } else if (!rigidHydrogenAngles || !isHydrogenAngle(angle)) {
239
240 int i1 = angle.getAtom(0).getArrayIndex();
241 int i2 = angle.getAtom(1).getArrayIndex();
242 int i3 = angle.getAtom(2).getArrayIndex();
243 double theta0 = angle.angleType.angle[angle.nh];
244 double k = OpenMM_KJPerKcal * angle.angleType.angleUnit * angle.angleType.forceConstant;
245 if (angleMode == AngleType.AngleMode.IN_PLANE) {
246
247 k = 0.0;
248 }
249 parameters.append(theta0);
250 parameters.append(k);
251 setAngleParameters(index++, i1, i2, i3, parameters);
252 parameters.resize(0);
253 }
254 }
255 parameters.destroy();
256 updateParametersInContext(openMMEnergy.getContext());
257 }
258
259
260
261
262
263
264
265 public void updateForce(int topology, OpenMMDualTopologyEnergy openMMDualTopologyEnergy) {
266 ForceFieldEnergy forceFieldEnergy = openMMDualTopologyEnergy.getForceFieldEnergy(topology);
267 AnglePotentialEnergy anglePotentialEnergy = forceFieldEnergy.getAnglePotentialEnergy();
268 Angle[] angles = anglePotentialEnergy.getAngleArray();
269 if (angles == null || angles.length < 1) {
270 return;
271 }
272
273 double scale = openMMDualTopologyEnergy.getTopologyScale(topology);
274
275 DoubleArray parameters = new DoubleArray(0);
276 int index = 0;
277 for (Angle angle : angles) {
278 AngleType.AngleMode angleMode = angle.angleType.angleMode;
279 if (angleMode == AngleType.AngleMode.IN_PLANE) {
280
281 }
282
283 else {
284 int i1 = angle.getAtom(0).getArrayIndex();
285 int i2 = angle.getAtom(1).getArrayIndex();
286 int i3 = angle.getAtom(2).getArrayIndex();
287 double theta0 = angle.angleType.angle[angle.nh];
288 double k = OpenMM_KJPerKcal * angle.angleType.angleUnit * angle.angleType.forceConstant;
289
290 if (!angle.applyLambda()) {
291 k = k * scale;
292 }
293 if (angleMode == AngleType.AngleMode.IN_PLANE) {
294
295 k = 0.0;
296 }
297 parameters.append(theta0);
298 parameters.append(k);
299 i1 = openMMDualTopologyEnergy.mapToDualTopologyIndex(topology, i1);
300 i2 = openMMDualTopologyEnergy.mapToDualTopologyIndex(topology, i2);
301 i3 = openMMDualTopologyEnergy.mapToDualTopologyIndex(topology, i3);
302 setAngleParameters(index++, i1, i2, i3, parameters);
303 parameters.resize(0);
304 }
305 }
306 parameters.destroy();
307 updateParametersInContext(openMMDualTopologyEnergy.getContext());
308 }
309
310
311
312
313
314
315
316
317 private boolean isHydrogenAngle(Angle angle) {
318 if (angle.containsHydrogen()) {
319
320 double angleVal = angle.angleType.angle[angle.nh];
321
322
323
324 if (angleVal < 160.0) {
325 Atom atom1 = angle.getAtom(0);
326 Atom atom2 = angle.getAtom(1);
327 Atom atom3 = angle.getAtom(2);
328
329
330 return atom1.isHydrogen() && atom3.isHydrogen() && !atom2.isHydrogen();
331 }
332 }
333 return false;
334 }
335
336 }