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.openmm.DoubleArray;
42 import ffx.openmm.Force;
43 import ffx.openmm.IntArray;
44 import ffx.openmm.amoeba.MultipoleForce;
45 import ffx.potential.bonded.Atom;
46 import ffx.potential.nonbonded.ParticleMeshEwald;
47 import ffx.potential.nonbonded.ReciprocalSpace;
48 import ffx.potential.nonbonded.pme.AlchemicalParameters;
49 import ffx.potential.nonbonded.pme.Polarization;
50 import ffx.potential.nonbonded.pme.SCFAlgorithm;
51 import ffx.potential.parameters.ForceField;
52 import ffx.potential.parameters.MultipoleType;
53 import ffx.potential.parameters.PolarizeType;
54
55 import java.util.logging.Level;
56 import java.util.logging.Logger;
57
58 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaMultipoleForce_CovalentType.OpenMM_AmoebaMultipoleForce_Covalent12;
59 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaMultipoleForce_CovalentType.OpenMM_AmoebaMultipoleForce_Covalent13;
60 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaMultipoleForce_CovalentType.OpenMM_AmoebaMultipoleForce_Covalent14;
61 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaMultipoleForce_CovalentType.OpenMM_AmoebaMultipoleForce_Covalent15;
62 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaMultipoleForce_CovalentType.OpenMM_AmoebaMultipoleForce_PolarizationCovalent11;
63 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaMultipoleForce_MultipoleAxisTypes.OpenMM_AmoebaMultipoleForce_Bisector;
64 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaMultipoleForce_MultipoleAxisTypes.OpenMM_AmoebaMultipoleForce_NoAxisType;
65 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaMultipoleForce_MultipoleAxisTypes.OpenMM_AmoebaMultipoleForce_ThreeFold;
66 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaMultipoleForce_MultipoleAxisTypes.OpenMM_AmoebaMultipoleForce_ZBisect;
67 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaMultipoleForce_MultipoleAxisTypes.OpenMM_AmoebaMultipoleForce_ZOnly;
68 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaMultipoleForce_MultipoleAxisTypes.OpenMM_AmoebaMultipoleForce_ZThenX;
69 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaMultipoleForce_NonbondedMethod.OpenMM_AmoebaMultipoleForce_NoCutoff;
70 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaMultipoleForce_NonbondedMethod.OpenMM_AmoebaMultipoleForce_PME;
71 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaMultipoleForce_PolarizationType.OpenMM_AmoebaMultipoleForce_Direct;
72 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaMultipoleForce_PolarizationType.OpenMM_AmoebaMultipoleForce_Extrapolated;
73 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaMultipoleForce_PolarizationType.OpenMM_AmoebaMultipoleForce_Mutual;
74 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_NmPerAngstrom;
75 import static java.lang.String.format;
76 import static org.apache.commons.math3.util.FastMath.sqrt;
77
78
79
80
81 public class AmoebaMultipoleForce extends MultipoleForce {
82
83 private static final Logger logger = Logger.getLogger(AmoebaMultipoleForce.class.getName());
84
85 public AmoebaMultipoleForce(OpenMMEnergy openMMEnergy) {
86 ParticleMeshEwald pme = openMMEnergy.getPmeNode();
87 if (pme == null) {
88 destroy();
89 return;
90 }
91
92 int[][] axisAtom = pme.getAxisAtoms();
93 double quadrupoleConversion = OpenMM_NmPerAngstrom * OpenMM_NmPerAngstrom;
94 double polarityConversion = OpenMM_NmPerAngstrom * OpenMM_NmPerAngstrom * OpenMM_NmPerAngstrom;
95 double dampingFactorConversion = sqrt(OpenMM_NmPerAngstrom);
96
97 ForceField forceField = openMMEnergy.getMolecularAssembly().getForceField();
98
99 Polarization polarization = pme.getPolarizationType();
100 double doPolarization = 1.0;
101 SCFAlgorithm scfAlgorithm = null;
102 if (polarization != Polarization.MUTUAL) {
103 setPolarizationType(OpenMM_AmoebaMultipoleForce_Direct);
104 if (pme.getPolarizationType() == Polarization.NONE) {
105 doPolarization = 0.0;
106 }
107 } else {
108 String algorithm = forceField.getString("SCF_ALGORITHM", "CG");
109 try {
110 algorithm = algorithm.replaceAll("-", "_").toUpperCase();
111 scfAlgorithm = SCFAlgorithm.valueOf(algorithm);
112 } catch (Exception e) {
113 scfAlgorithm = SCFAlgorithm.CG;
114 }
115
116 if (scfAlgorithm == SCFAlgorithm.EPT) {
117
118
119
120
121
122 setPolarizationType(OpenMM_AmoebaMultipoleForce_Extrapolated);
123 DoubleArray exptCoefficients = new DoubleArray(4);
124 exptCoefficients.set(0, -0.154);
125 exptCoefficients.set(1, 0.017);
126 exptCoefficients.set(2, 0.657);
127 exptCoefficients.set(3, 0.475);
128 setExtrapolationCoefficients(exptCoefficients);
129 exptCoefficients.destroy();
130 } else {
131 setPolarizationType(OpenMM_AmoebaMultipoleForce_Mutual);
132 }
133 }
134
135 DoubleArray dipoles = new DoubleArray(3);
136 DoubleArray quadrupoles = new DoubleArray(9);
137
138 AlchemicalParameters alchemicalParameters = pme.getAlchemicalParameters();
139 boolean lambdaTerm = pme.getLambdaTerm();
140 double permLambda = alchemicalParameters.permLambda;
141 double polarLambda = alchemicalParameters.polLambda;
142
143 Atom[] atoms = openMMEnergy.getMolecularAssembly().getAtomArray();
144 int nAtoms = atoms.length;
145 for (int i = 0; i < nAtoms; i++) {
146 Atom atom = atoms[i];
147 MultipoleType multipoleType = pme.getMultipoleType(i);
148 PolarizeType polarType = pme.getPolarizeType(i);
149
150
151 int axisType = switch (multipoleType.frameDefinition) {
152 case NONE -> OpenMM_AmoebaMultipoleForce_NoAxisType;
153 case ZONLY -> OpenMM_AmoebaMultipoleForce_ZOnly;
154 case ZTHENX -> OpenMM_AmoebaMultipoleForce_ZThenX;
155 case BISECTOR -> OpenMM_AmoebaMultipoleForce_Bisector;
156 case ZTHENBISECTOR -> OpenMM_AmoebaMultipoleForce_ZBisect;
157 case THREEFOLD -> OpenMM_AmoebaMultipoleForce_ThreeFold;
158 };
159
160 double useFactor = 1.0;
161 if (!atoms[i].getUse() || !atoms[i].getElectrostatics()) {
162 useFactor = 0.0;
163 }
164
165 double permScale = useFactor;
166 double polarScale = doPolarization * useFactor;
167 if (lambdaTerm && atom.applyLambda()) {
168 permScale *= permLambda;
169 polarScale *= polarLambda;
170 }
171
172
173 for (int j = 0; j < 3; j++) {
174 dipoles.set(j, multipoleType.dipole[j] * OpenMM_NmPerAngstrom * permScale);
175 }
176 int l = 0;
177 for (int j = 0; j < 3; j++) {
178 for (int k = 0; k < 3; k++) {
179 quadrupoles.set(l++, multipoleType.quadrupole[j][k] * quadrupoleConversion * permScale / 3.0);
180 }
181 }
182
183 int zaxis = -1;
184 int xaxis = -1;
185 int yaxis = -1;
186 int[] refAtoms = axisAtom[i];
187 if (refAtoms != null) {
188 zaxis = refAtoms[0];
189 if (refAtoms.length > 1) {
190 xaxis = refAtoms[1];
191 if (refAtoms.length > 2) {
192 yaxis = refAtoms[2];
193 }
194 }
195 } else {
196 axisType = OpenMM_AmoebaMultipoleForce_NoAxisType;
197 }
198
199 double charge = multipoleType.charge * permScale;
200
201
202 addMultipole(charge, dipoles, quadrupoles, axisType, zaxis, xaxis, yaxis, polarType.thole,
203 polarType.pdamp * dampingFactorConversion, polarType.polarizability * polarityConversion * polarScale);
204 }
205 dipoles.destroy();
206 quadrupoles.destroy();
207
208 Crystal crystal = openMMEnergy.getCrystal();
209 double cutoff = pme.getEwaldCutoff();
210 double aewald = pme.getEwaldCoefficient();
211 if (!crystal.aperiodic()) {
212 setNonbondedMethod(OpenMM_AmoebaMultipoleForce_PME);
213 setCutoffDistance(cutoff * OpenMM_NmPerAngstrom);
214 setAEwald(aewald / OpenMM_NmPerAngstrom);
215
216 double ewaldTolerance = 1.0e-04;
217 setEwaldErrorTolerance(ewaldTolerance);
218
219 IntArray gridDimensions = new IntArray(3);
220 ReciprocalSpace recip = pme.getReciprocalSpace();
221 gridDimensions.set(0, recip.getXDim());
222 gridDimensions.set(1, recip.getYDim());
223 gridDimensions.set(2, recip.getZDim());
224 setPmeGridDimensions(gridDimensions);
225 gridDimensions.destroy();
226 } else {
227 setNonbondedMethod(OpenMM_AmoebaMultipoleForce_NoCutoff);
228 }
229
230 setMutualInducedMaxIterations(500);
231 double poleps = pme.getPolarEps();
232 setMutualInducedTargetEpsilon(poleps);
233
234 int[][] ip11 = pme.getPolarization11();
235
236 IntArray covalentMap = new IntArray(0);
237 for (int i = 0; i < nAtoms; i++) {
238 Atom ai = atoms[i];
239
240
241 covalentMap.resize(0);
242 for (Atom ak : ai.get12List()) {
243 covalentMap.append(ak.getIndex() - 1);
244 }
245 setCovalentMap(i, OpenMM_AmoebaMultipoleForce_Covalent12, covalentMap);
246
247
248 covalentMap.resize(0);
249 for (Atom ak : ai.get13List()) {
250 covalentMap.append(ak.getIndex() - 1);
251 }
252 setCovalentMap(i, OpenMM_AmoebaMultipoleForce_Covalent13, covalentMap);
253
254
255 covalentMap.resize(0);
256 for (Atom ak : ai.get14List()) {
257 covalentMap.append(ak.getIndex() - 1);
258 }
259 setCovalentMap(i, OpenMM_AmoebaMultipoleForce_Covalent14, covalentMap);
260
261
262 covalentMap.resize(0);
263 for (Atom ak : ai.get15List()) {
264 covalentMap.append(ak.getIndex() - 1);
265 }
266 setCovalentMap(i, OpenMM_AmoebaMultipoleForce_Covalent15, covalentMap);
267
268
269 covalentMap.resize(0);
270 for (int j = 0; j < ip11[i].length; j++) {
271 covalentMap.append(ip11[i][j]);
272 }
273 setCovalentMap(i, OpenMM_AmoebaMultipoleForce_PolarizationCovalent11, covalentMap);
274
275
276 }
277 covalentMap.destroy();
278
279 int forceGroup = forceField.getInteger("PME_FORCE_GROUP", 1);
280 setForceGroup(forceGroup);
281 if (logger.isLoggable(Level.INFO)) {
282 StringBuilder sb = new StringBuilder();
283 sb.append(" Electrostatics\n");
284 sb.append(format(" Polarization: %8s\n", polarization.toString()));
285 if (polarization == Polarization.MUTUAL) {
286 sb.append(format(" SCF Convergence Criteria: %8.3e\n", poleps));
287 } else if (scfAlgorithm == SCFAlgorithm.EPT) {
288 sb.append(format(" SCF Algorithm: %8s\n", scfAlgorithm));
289 }
290 if (aewald > 0.0) {
291 sb.append(" Particle-mesh Ewald\n");
292 sb.append(format(" Ewald Coefficient: %8.3f\n", aewald));
293 sb.append(format(" Particle Cut-Off: %8.3f (A)", cutoff));
294 } else if (cutoff != Double.POSITIVE_INFINITY) {
295 sb.append(format(" Electrostatics Cut-Off: %8.3f (A)", cutoff));
296 } else {
297 sb.append(" Electrostatics Cut-Off: NONE");
298 }
299 logger.info(sb.toString());
300
301 if (lambdaTerm) {
302 sb = new StringBuilder(" Alchemical Parameters\n");
303 double permLambdaStart = alchemicalParameters.permLambdaStart;
304 double permLambdaEnd = alchemicalParameters.permLambdaEnd;
305 sb.append(format(" Permanent Multipole Range: %5.3f-%5.3f\n", permLambdaStart, permLambdaEnd));
306 double permLambdaAlpha = alchemicalParameters.permLambdaAlpha;
307 if (permLambdaAlpha != 0.0) {
308 logger.severe(" Permanent multipole softcore not supported for OpenMM.");
309 }
310 double permLambdaExponent = alchemicalParameters.permLambdaExponent;
311 sb.append(format(" Permanent Multipole Lambda Exponent: %5.3f\n", permLambdaExponent));
312 if (polarization != Polarization.NONE) {
313 double polLambdaExponent = alchemicalParameters.polLambdaExponent;
314 double polLambdaStart = alchemicalParameters.polLambdaStart;
315 double polLambdaEnd = alchemicalParameters.polLambdaEnd;
316 sb.append(format(" Polarization Lambda Exponent: %5.3f\n", polLambdaExponent));
317 sb.append(format(" Polarization Range: %5.3f-%5.3f\n", polLambdaStart, polLambdaEnd));
318 boolean doNoLigandCondensedSCF = alchemicalParameters.doNoLigandCondensedSCF;
319 if (doNoLigandCondensedSCF) {
320 logger.severe(" Condensed SCF without a ligand is not supported for OpenMM.");
321 }
322 }
323 boolean doLigandGKElec = alchemicalParameters.doLigandGKElec;
324 if (doLigandGKElec) {
325 logger.severe(" Isolated ligand electrostatics are not supported for OpenMM.");
326 }
327 logger.info(sb.toString());
328 }
329
330 logger.log(Level.FINE, format(" Force group:\t\t%d\n", forceGroup));
331 }
332 }
333
334
335
336
337
338
339
340 public static Force constructForce(OpenMMEnergy openMMEnergy) {
341 ParticleMeshEwald pme = openMMEnergy.getPmeNode();
342 if (pme == null) {
343 return null;
344 }
345 return new AmoebaMultipoleForce(openMMEnergy);
346 }
347
348 public void updateForce(Atom[] atoms, OpenMMEnergy openMMEnergy) {
349 ParticleMeshEwald pme = openMMEnergy.getPmeNode();
350 double quadrupoleConversion = OpenMM_NmPerAngstrom * OpenMM_NmPerAngstrom;
351 double polarityConversion = OpenMM_NmPerAngstrom * OpenMM_NmPerAngstrom * OpenMM_NmPerAngstrom;
352 double dampingFactorConversion = sqrt(OpenMM_NmPerAngstrom);
353
354 double doPolarization = 1.0;
355 if (pme.getPolarizationType() == Polarization.NONE) {
356 doPolarization = 0.0;
357 }
358
359 DoubleArray dipoles = new DoubleArray(3);
360 DoubleArray quadrupoles = new DoubleArray(9);
361
362 AlchemicalParameters alchemicalParameters = pme.getAlchemicalParameters();
363 boolean lambdaTerm = pme.getLambdaTerm();
364 if (lambdaTerm) {
365 AlchemicalParameters.AlchemicalMode alchemicalMode = alchemicalParameters.mode;
366
367 if (alchemicalMode != AlchemicalParameters.AlchemicalMode.SCALE) {
368 logger.severe(format(" Alchemical mode %s not supported for OpenMM.", alchemicalMode));
369 }
370
371 if (alchemicalParameters.permLambdaAlpha != 0.0) {
372 logger.severe(" Permanent multipole softcore not supported for OpenMM.");
373 }
374
375 if (alchemicalParameters.doLigandGKElec || alchemicalParameters.doLigandVaporElec) {
376 logger.severe(" Isolated ligand electrostatics are not supported for OpenMM.");
377 }
378
379 if (alchemicalParameters.doNoLigandCondensedSCF) {
380 logger.severe(" Condensed SCF without a ligand is not supported for OpenMM.");
381 }
382 }
383
384 double permLambda = alchemicalParameters.permLambda;
385 double polarLambda = alchemicalParameters.polLambda;
386
387 for (Atom atom : atoms) {
388 int index = atom.getXyzIndex() - 1;
389 MultipoleType multipoleType = pme.getMultipoleType(index);
390 PolarizeType polarizeType = pme.getPolarizeType(index);
391 int[] axisAtoms = atom.getAxisAtomIndices();
392
393 double permScale = 1.0;
394 double polarScale = doPolarization;
395
396 if (!atom.getUse() || !atom.getElectrostatics()) {
397 permScale = 0.0;
398 polarScale = 0.0;
399 }
400
401 if (atom.applyLambda()) {
402 permScale *= permLambda;
403 polarScale *= polarLambda;
404 }
405
406
407 int axisType = switch (multipoleType.frameDefinition) {
408 case NONE -> OpenMM_AmoebaMultipoleForce_NoAxisType;
409 case ZONLY -> OpenMM_AmoebaMultipoleForce_ZOnly;
410 case ZTHENX -> OpenMM_AmoebaMultipoleForce_ZThenX;
411 case BISECTOR -> OpenMM_AmoebaMultipoleForce_Bisector;
412 case ZTHENBISECTOR -> OpenMM_AmoebaMultipoleForce_ZBisect;
413 case THREEFOLD -> OpenMM_AmoebaMultipoleForce_ThreeFold;
414 };
415
416
417 for (int j = 0; j < 3; j++) {
418 dipoles.set(j, multipoleType.dipole[j] * OpenMM_NmPerAngstrom * permScale);
419 }
420 int l = 0;
421 for (int j = 0; j < 3; j++) {
422 for (int k = 0; k < 3; k++) {
423 quadrupoles.set(l++, multipoleType.quadrupole[j][k] * quadrupoleConversion / 3.0 * permScale);
424 }
425 }
426
427 int zaxis = 1;
428 int xaxis = 1;
429 int yaxis = 1;
430
431 if (axisAtoms != null) {
432 zaxis = axisAtoms[0];
433 if (axisAtoms.length > 1) {
434 xaxis = axisAtoms[1];
435 if (axisAtoms.length > 2) {
436 yaxis = axisAtoms[2];
437 }
438 }
439 } else {
440 axisType = OpenMM_AmoebaMultipoleForce_NoAxisType;
441 }
442
443
444 setMultipoleParameters(index, multipoleType.charge * permScale,
445 dipoles, quadrupoles, axisType, zaxis, xaxis, yaxis,
446 polarizeType.thole, polarizeType.pdamp * dampingFactorConversion,
447 polarizeType.polarizability * polarityConversion * polarScale);
448 }
449
450 dipoles.destroy();
451 quadrupoles.destroy();
452 updateParametersInContext(openMMEnergy.getContext());
453 }
454
455 }