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