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.crystal.Crystal;
41 import ffx.numerics.switching.UnivariateSwitchingFunction;
42 import ffx.potential.DualTopologyEnergy;
43 import ffx.potential.ForceFieldEnergy;
44 import ffx.potential.MolecularAssembly;
45 import ffx.potential.Platform;
46 import ffx.potential.bonded.Atom;
47 import ffx.potential.parameters.ForceField;
48 import ffx.potential.utils.EnergyException;
49
50 import javax.annotation.Nullable;
51 import java.util.logging.Level;
52 import java.util.logging.Logger;
53
54 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_Boolean.OpenMM_True;
55 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_State_DataType.OpenMM_State_Energy;
56 import static java.lang.Double.isFinite;
57 import static java.lang.String.format;
58
59 public class OpenMMDualTopologyEnergy extends DualTopologyEnergy implements OpenMMPotential {
60
61 private static final Logger logger = Logger.getLogger(OpenMMDualTopologyEnergy.class.getName());
62
63
64
65
66 private final OpenMMContext openMMContext;
67
68
69
70 private final OpenMMDualTopologySystem openMMDualTopologySystem;
71
72
73
74 private final Atom[] atoms;
75
76
77
78 private final MolecularAssembly molecularAssembly1;
79
80
81
82 private final MolecularAssembly molecularAssembly2;
83
84
85
86
87
88
89
90
91
92 public OpenMMDualTopologyEnergy(MolecularAssembly topology1, MolecularAssembly topology2,
93 UnivariateSwitchingFunction switchFunction, Platform requestedPlatform) {
94 super(topology1, topology2, switchFunction);
95
96 logger.info("\n Initializing an OpenMM Dual Topology System.");
97
98 molecularAssembly1 = topology1;
99 molecularAssembly2 = topology2;
100
101
102 Crystal crystal1 = molecularAssembly1.getCrystal();
103 int symOps1 = crystal1.spaceGroup.getNumberOfSymOps();
104 Crystal crystal2 = molecularAssembly2.getCrystal();
105 int symOps2 = crystal2.spaceGroup.getNumberOfSymOps();
106 if (symOps1 > 1 || symOps2 > 1) {
107 logger.severe(" OpenMM does not support symmetry operators.");
108 }
109
110
111 ForceField forceField = topology1.getForceField();
112 ffx.openmm.Platform openMMPlatform = OpenMMContext.loadPlatform(requestedPlatform, forceField);
113
114
115 openMMDualTopologySystem = new OpenMMDualTopologySystem(this);
116 openMMDualTopologySystem.addForces();
117
118 atoms = getDualTopologyAtoms(0);
119
120
121 openMMContext = new OpenMMContext(openMMPlatform, openMMDualTopologySystem, atoms);
122 }
123
124
125
126
127
128 @Override
129 public double energy(double[] x) {
130 return energy(x, false);
131 }
132
133
134
135
136 @Override
137 public double energy(double[] x, boolean verbose) {
138
139 openMMContext.update();
140
141 updateParameters(atoms);
142
143
144 unscaleCoordinates(x);
145 setCoordinates(x);
146
147 OpenMMState openMMState = openMMContext.getOpenMMState(OpenMM_State_Energy);
148 double e = openMMState.potentialEnergy;
149 openMMState.destroy();
150
151 if (!isFinite(e)) {
152 String message = format(" OpenMMDualTopologyEnergy was a non-finite %8g", e);
153 logger.warning(message);
154 throw new EnergyException(message);
155 }
156
157 if (verbose) {
158 logger.log(Level.INFO, format("\n OpenMM Energy: %14.10g", e));
159 }
160
161 return e;
162 }
163
164
165
166
167
168
169
170 public double energyFFX(double[] x) {
171 return super.energy(x, false);
172 }
173
174
175
176
177
178
179
180
181 public double energyFFX(double[] x, boolean verbose) {
182 return super.energy(x, verbose);
183 }
184
185
186
187
188
189
190
191 @Override
192 public void setCoordinates(double[] x) {
193
194 super.setCoordinates(x);
195
196 int n = atoms.length * 3;
197 double[] xall = new double[n];
198 int i = 0;
199 for (Atom atom : atoms) {
200 xall[i] = atom.getX();
201 xall[i + 1] = atom.getY();
202 xall[i + 2] = atom.getZ();
203 i += 3;
204 }
205
206
207 openMMContext.setPositions(xall);
208 }
209
210
211
212
213
214
215 @Override
216 public void setVelocity(double[] v) {
217
218 super.setVelocity(v);
219
220 int n = atoms.length * 3;
221 double[] vall = new double[n];
222 double[] v3 = new double[3];
223 int i = 0;
224 for (Atom atom : atoms) {
225 atom.getVelocity(v3);
226
227 if (!atom.isActive()) {
228 v3[0] = 0.0;
229 v3[1] = 0.0;
230 v3[2] = 0.0;
231 atom.setVelocity(v3);
232 }
233 vall[i] = v3[0];
234 vall[i + 1] = v3[1];
235 vall[i + 2] = v3[2];
236 i += 3;
237 }
238
239
240 openMMContext.setVelocities(vall);
241 }
242
243
244
245
246
247
248
249
250 public MolecularAssembly getMolecularAssembly(int topology) {
251 if (topology == 0) {
252 return molecularAssembly1;
253 } else if (topology == 1) {
254 return molecularAssembly2;
255 } else {
256 throw new IllegalArgumentException(" Invalid topology index: " + topology);
257 }
258 }
259
260
261
262
263
264
265
266 public ForceFieldEnergy getForceFieldEnergy(int topology) {
267 if (topology == 0) {
268 return getForceFieldEnergy1();
269 } else if (topology == 1) {
270 return getForceFieldEnergy2();
271 } else {
272 throw new IllegalArgumentException(" Invalid topology index: " + topology);
273 }
274 }
275
276
277
278
279
280
281 @Override
282 public void updateParameters(@Nullable Atom[] atoms) {
283 if (atoms == null) {
284 atoms = this.atoms;
285 }
286 if (openMMDualTopologySystem != null) {
287 openMMDualTopologySystem.updateParameters(atoms);
288 }
289 }
290
291
292
293
294
295
296 @Override
297 public OpenMMContext getContext() {
298 return openMMContext;
299 }
300
301
302
303
304
305
306
307
308
309
310
311
312 @Override
313 public void updateContext(String integratorName, double timeStep, double temperature, boolean forceCreation) {
314 openMMContext.update(integratorName, timeStep, temperature, forceCreation);
315 }
316
317
318
319
320
321
322
323
324
325 @Override
326 public OpenMMState getOpenMMState(int mask) {
327 return openMMContext.getOpenMMState(mask);
328 }
329
330
331
332
333
334
335 @Override
336 public OpenMMSystem getSystem() {
337 return openMMDualTopologySystem;
338 }
339
340
341
342
343 @Override
344 public boolean setActiveAtoms() {
345 return openMMDualTopologySystem.updateAtomMass();
346 }
347
348 }