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.ForceFieldEnergy;
46 import ffx.potential.bonded.Atom;
47 import ffx.potential.nonbonded.ParticleMeshEwald;
48 import ffx.potential.nonbonded.ReciprocalSpace;
49 import ffx.potential.nonbonded.pme.AlchemicalParameters;
50 import ffx.potential.nonbonded.pme.Polarization;
51 import ffx.potential.nonbonded.pme.SCFAlgorithm;
52 import ffx.potential.parameters.ForceField;
53 import ffx.potential.parameters.MultipoleType;
54 import ffx.potential.parameters.PolarizeType;
55
56 import java.util.HashSet;
57 import java.util.Set;
58 import java.util.logging.Level;
59 import java.util.logging.Logger;
60
61 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaMultipoleForce_CovalentType.OpenMM_AmoebaMultipoleForce_Covalent12;
62 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaMultipoleForce_CovalentType.OpenMM_AmoebaMultipoleForce_Covalent13;
63 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaMultipoleForce_CovalentType.OpenMM_AmoebaMultipoleForce_Covalent14;
64 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaMultipoleForce_CovalentType.OpenMM_AmoebaMultipoleForce_Covalent15;
65 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaMultipoleForce_CovalentType.OpenMM_AmoebaMultipoleForce_PolarizationCovalent11;
66 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaMultipoleForce_MultipoleAxisTypes.OpenMM_AmoebaMultipoleForce_Bisector;
67 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaMultipoleForce_MultipoleAxisTypes.OpenMM_AmoebaMultipoleForce_NoAxisType;
68 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaMultipoleForce_MultipoleAxisTypes.OpenMM_AmoebaMultipoleForce_ThreeFold;
69 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaMultipoleForce_MultipoleAxisTypes.OpenMM_AmoebaMultipoleForce_ZBisect;
70 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaMultipoleForce_MultipoleAxisTypes.OpenMM_AmoebaMultipoleForce_ZOnly;
71 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaMultipoleForce_MultipoleAxisTypes.OpenMM_AmoebaMultipoleForce_ZThenX;
72 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaMultipoleForce_NonbondedMethod.OpenMM_AmoebaMultipoleForce_NoCutoff;
73 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaMultipoleForce_NonbondedMethod.OpenMM_AmoebaMultipoleForce_PME;
74 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaMultipoleForce_PolarizationType.OpenMM_AmoebaMultipoleForce_Direct;
75 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaMultipoleForce_PolarizationType.OpenMM_AmoebaMultipoleForce_Extrapolated;
76 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaMultipoleForce_PolarizationType.OpenMM_AmoebaMultipoleForce_Mutual;
77 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_NmPerAngstrom;
78 import static java.lang.String.format;
79 import static org.apache.commons.math3.util.FastMath.sqrt;
80
81
82
83
84 public class AmoebaMultipoleForce extends MultipoleForce {
85
86 private static final Logger logger = Logger.getLogger(AmoebaMultipoleForce.class.getName());
87
88
89
90
91
92
93 public AmoebaMultipoleForce(OpenMMEnergy openMMEnergy) {
94 ParticleMeshEwald pme = openMMEnergy.getPmeNode();
95 if (pme == null) {
96 destroy();
97 return;
98 }
99
100 double doPolarization = configureForce(openMMEnergy);
101
102 int[][] axisAtom = pme.getAxisAtoms();
103 double quadrupoleConversion = OpenMM_NmPerAngstrom * OpenMM_NmPerAngstrom;
104 double polarityConversion = OpenMM_NmPerAngstrom * OpenMM_NmPerAngstrom * OpenMM_NmPerAngstrom;
105 double dampingFactorConversion = sqrt(OpenMM_NmPerAngstrom);
106
107 boolean lambdaTerm = pme.getLambdaTerm();
108 AlchemicalParameters alchemicalParameters = pme.getAlchemicalParameters();
109 double permLambda = alchemicalParameters.permLambda;
110 double polarLambda = alchemicalParameters.polLambda;
111 DoubleArray dipoles = new DoubleArray(3);
112 DoubleArray quadrupoles = new DoubleArray(9);
113
114 Atom[] atoms = openMMEnergy.getMolecularAssembly().getAtomArray();
115 int nAtoms = atoms.length;
116 for (int i = 0; i < nAtoms; i++) {
117 Atom atom = atoms[i];
118 MultipoleType multipoleType = pme.getMultipoleType(i);
119 PolarizeType polarType = pme.getPolarizeType(i);
120
121
122 int axisType = switch (multipoleType.frameDefinition) {
123 case NONE -> OpenMM_AmoebaMultipoleForce_NoAxisType;
124 case ZONLY -> OpenMM_AmoebaMultipoleForce_ZOnly;
125 case ZTHENX -> OpenMM_AmoebaMultipoleForce_ZThenX;
126 case BISECTOR -> OpenMM_AmoebaMultipoleForce_Bisector;
127 case ZTHENBISECTOR -> OpenMM_AmoebaMultipoleForce_ZBisect;
128 case THREEFOLD -> OpenMM_AmoebaMultipoleForce_ThreeFold;
129 };
130
131 double useFactor = 1.0;
132 if (!atom.getUse() || !atom.getElectrostatics()) {
133 useFactor = 0.0;
134 }
135
136 double permScale = useFactor;
137 double polarScale = doPolarization * useFactor;
138 if (lambdaTerm && atom.applyLambda()) {
139 permScale *= permLambda;
140 polarScale *= polarLambda;
141 }
142
143
144 for (int j = 0; j < 3; j++) {
145 dipoles.set(j, multipoleType.dipole[j] * OpenMM_NmPerAngstrom * permScale);
146 }
147 int l = 0;
148 for (int j = 0; j < 3; j++) {
149 for (int k = 0; k < 3; k++) {
150 quadrupoles.set(l++, multipoleType.quadrupole[j][k] * quadrupoleConversion * permScale / 3.0);
151 }
152 }
153
154 int zaxis = -1;
155 int xaxis = -1;
156 int yaxis = -1;
157 int[] refAtoms = axisAtom[i];
158 if (refAtoms != null) {
159 zaxis = refAtoms[0];
160 if (refAtoms.length > 1) {
161 xaxis = refAtoms[1];
162 if (refAtoms.length > 2) {
163 yaxis = refAtoms[2];
164 }
165 }
166 } else {
167 axisType = OpenMM_AmoebaMultipoleForce_NoAxisType;
168 }
169
170 double charge = multipoleType.charge * permScale;
171
172
173 addMultipole(charge, dipoles, quadrupoles, axisType, zaxis, xaxis, yaxis, polarType.thole,
174 polarType.pdamp * dampingFactorConversion, polarType.polarizability * polarityConversion * polarScale);
175 }
176 dipoles.destroy();
177 quadrupoles.destroy();
178
179 int[][] ip11 = pme.getPolarization11();
180 IntArray covalentMap = new IntArray(0);
181 for (int i = 0; i < nAtoms; i++) {
182 Atom ai = atoms[i];
183
184
185 covalentMap.resize(0);
186 for (Atom ak : ai.get12List()) {
187 covalentMap.append(ak.getIndex() - 1);
188 }
189 setCovalentMap(i, OpenMM_AmoebaMultipoleForce_Covalent12, covalentMap);
190
191
192 covalentMap.resize(0);
193 for (Atom ak : ai.get13List()) {
194 covalentMap.append(ak.getIndex() - 1);
195 }
196 setCovalentMap(i, OpenMM_AmoebaMultipoleForce_Covalent13, covalentMap);
197
198
199 covalentMap.resize(0);
200 for (Atom ak : ai.get14List()) {
201 covalentMap.append(ak.getIndex() - 1);
202 }
203 setCovalentMap(i, OpenMM_AmoebaMultipoleForce_Covalent14, covalentMap);
204
205
206 covalentMap.resize(0);
207 for (Atom ak : ai.get15List()) {
208 covalentMap.append(ak.getIndex() - 1);
209 }
210 setCovalentMap(i, OpenMM_AmoebaMultipoleForce_Covalent15, covalentMap);
211
212
213 covalentMap.resize(0);
214 for (int j = 0; j < ip11[i].length; j++) {
215 covalentMap.append(ip11[i][j]);
216 }
217 setCovalentMap(i, OpenMM_AmoebaMultipoleForce_PolarizationCovalent11, covalentMap);
218
219
220 }
221 covalentMap.destroy();
222 }
223
224
225
226
227
228
229
230 public AmoebaMultipoleForce(int topology, OpenMMDualTopologyEnergy openMMDualTopologyEnergy) {
231
232 int otherTopology = 1 - topology;
233
234 ForceFieldEnergy forceFieldEnergy = openMMDualTopologyEnergy.getForceFieldEnergy(topology);
235 ForceFieldEnergy otherForceFieldEnergy = openMMDualTopologyEnergy.getForceFieldEnergy(otherTopology);
236
237 ParticleMeshEwald pme = forceFieldEnergy.getPmeNode();
238 ParticleMeshEwald otherPME = otherForceFieldEnergy.getPmeNode();
239 if (pme == null || otherPME == null) {
240 destroy();
241 return;
242 }
243
244 double doPolarization = configureForce(forceFieldEnergy);
245
246
247
248
249 double quadrupoleConversion = OpenMM_NmPerAngstrom * OpenMM_NmPerAngstrom;
250 double polarityConversion = OpenMM_NmPerAngstrom * OpenMM_NmPerAngstrom * OpenMM_NmPerAngstrom;
251 double dampingFactorConversion = sqrt(OpenMM_NmPerAngstrom);
252
253 boolean lambdaTerm = pme.getLambdaTerm();
254 AlchemicalParameters alchemicalParameters = pme.getAlchemicalParameters();
255 double permLambda = alchemicalParameters.permLambda;
256 double polarLambda = alchemicalParameters.polLambda;
257 DoubleArray dipoles = new DoubleArray(3);
258 DoubleArray quadrupoles = new DoubleArray(9);
259
260 int nAtoms = openMMDualTopologyEnergy.getNumberOfAtoms();
261
262
263 for (int i = 0; i < nAtoms; i++) {
264 Atom atom = openMMDualTopologyEnergy.getDualTopologyAtom(topology, i);
265 if (atom.getTopologyIndex() == topology) {
266 int index = atom.getArrayIndex();
267 MultipoleType multipoleType = pme.getMultipoleType(index);
268 PolarizeType polarType = pme.getPolarizeType(index);
269
270
271 int axisType = switch (multipoleType.frameDefinition) {
272 case NONE -> OpenMM_AmoebaMultipoleForce_NoAxisType;
273 case ZONLY -> OpenMM_AmoebaMultipoleForce_ZOnly;
274 case ZTHENX -> OpenMM_AmoebaMultipoleForce_ZThenX;
275 case BISECTOR -> OpenMM_AmoebaMultipoleForce_Bisector;
276 case ZTHENBISECTOR -> OpenMM_AmoebaMultipoleForce_ZBisect;
277 case THREEFOLD -> OpenMM_AmoebaMultipoleForce_ThreeFold;
278 };
279
280
281
282 double useFactor = 1.0;
283
284
285 if (!atom.getUse() || !atom.getElectrostatics()) {
286
287 useFactor = 0.0;
288 }
289
290
291 double permScale = useFactor;
292 double polarScale = doPolarization * useFactor;
293 if (lambdaTerm && atom.applyLambda()) {
294 permScale *= permLambda;
295 polarScale *= polarLambda;
296 }
297
298
299 double charge = multipoleType.charge * permScale;
300 for (int j = 0; j < 3; j++) {
301 dipoles.set(j, multipoleType.dipole[j] * OpenMM_NmPerAngstrom * permScale);
302 }
303 int l = 0;
304 for (int j = 0; j < 3; j++) {
305 for (int k = 0; k < 3; k++) {
306 quadrupoles.set(l++, multipoleType.quadrupole[j][k] * quadrupoleConversion * permScale / 3.0);
307 }
308 }
309
310 int zaxis = -1;
311 int xaxis = -1;
312 int yaxis = -1;
313 int[] refAtoms = atom.getAxisAtomIndices();
314 if (refAtoms != null) {
315 zaxis = refAtoms[0];
316 zaxis = openMMDualTopologyEnergy.mapToDualTopologyIndex(topology, zaxis);
317 if (refAtoms.length > 1) {
318 xaxis = refAtoms[1];
319 xaxis = openMMDualTopologyEnergy.mapToDualTopologyIndex(topology, xaxis);
320 if (refAtoms.length > 2) {
321 yaxis = refAtoms[2];
322 yaxis = openMMDualTopologyEnergy.mapToDualTopologyIndex(topology, yaxis);
323 }
324 }
325 } else {
326 axisType = OpenMM_AmoebaMultipoleForce_NoAxisType;
327 }
328
329
330 addMultipole(charge, dipoles, quadrupoles, axisType, zaxis, xaxis, yaxis,
331 polarType.thole, polarType.pdamp * dampingFactorConversion,
332 polarType.polarizability * polarityConversion * polarScale);
333 } else {
334
335
336 int axisType = OpenMM_AmoebaMultipoleForce_NoAxisType;
337
338 double charge = 0.0;
339 for (int j = 0; j < 3; j++) {
340 dipoles.set(j, 0.0);
341 }
342 int l = 0;
343 for (int j = 0; j < 3; j++) {
344 for (int k = 0; k < 3; k++) {
345 quadrupoles.set(l++, 0.0);
346 }
347 }
348
349 int zaxis = -1;
350 int xaxis = -1;
351 int yaxis = -1;
352
353 double thole = 0.39;
354 double pdamp = 1.0;
355
356 addMultipole(charge, dipoles, quadrupoles, axisType, zaxis, xaxis, yaxis,
357 thole, pdamp, 0.0);
358 }
359 }
360 dipoles.destroy();
361 quadrupoles.destroy();
362
363 int[][] ip11 = pme.getPolarization11();
364 int[][] ip11Other = otherPME.getPolarization11();
365
366 IntArray covalentMap = new IntArray(0);
367
368 Set<Integer> covalentSet = new HashSet<>();
369
370 for (int i = 0; i < nAtoms; i++) {
371 Atom atom = openMMDualTopologyEnergy.getDualTopologyAtom(topology, i);
372 Atom otherAtom = openMMDualTopologyEnergy.getDualTopologyAtom(otherTopology, i);
373 int index = atom.getArrayIndex();
374 int otherIndex = otherAtom.getArrayIndex();
375
376
377 covalentSet.clear();
378 for (Atom ak : atom.get12List()) {
379 covalentSet.add(ak.getTopologyAtomIndex());
380 }
381 for (Atom ak : otherAtom.get12List()) {
382 covalentSet.add(ak.getTopologyAtomIndex());
383 }
384 covalentMap.resize(0);
385 for (int ak : covalentSet) {
386 covalentMap.append(ak);
387 }
388 setCovalentMap(i, OpenMM_AmoebaMultipoleForce_Covalent12, covalentMap);
389
390
391 covalentSet.clear();
392 for (Atom ak : atom.get13List()) {
393 covalentSet.add(ak.getTopologyAtomIndex());
394 }
395 for (Atom ak : otherAtom.get13List()) {
396 covalentSet.add(ak.getTopologyAtomIndex());
397 }
398 covalentMap.resize(0);
399 for (int ak : covalentSet) {
400 covalentMap.append(ak);
401 }
402 setCovalentMap(i, OpenMM_AmoebaMultipoleForce_Covalent13, covalentMap);
403
404
405 covalentSet.clear();
406 for (Atom ak : atom.get14List()) {
407 covalentSet.add(ak.getTopologyAtomIndex());
408 }
409 for (Atom ak : otherAtom.get14List()) {
410 covalentSet.add(ak.getTopologyAtomIndex());
411 }
412 covalentMap.resize(0);
413 for (int ak : covalentSet) {
414 covalentMap.append(ak);
415 }
416 setCovalentMap(i, OpenMM_AmoebaMultipoleForce_Covalent14, covalentMap);
417
418
419 covalentSet.clear();
420 for (Atom ak : atom.get15List()) {
421 covalentSet.add(ak.getTopologyAtomIndex());
422 }
423 for (Atom ak : otherAtom.get15List()) {
424 covalentSet.add(ak.getTopologyAtomIndex());
425 }
426 covalentMap.resize(0);
427 for (int ak : covalentSet) {
428 covalentMap.append(ak);
429 }
430 setCovalentMap(i, OpenMM_AmoebaMultipoleForce_Covalent15, covalentMap);
431
432
433 covalentSet.clear();
434 for (int k : ip11[index]) {
435 int value = openMMDualTopologyEnergy.mapToDualTopologyIndex(topology, k);
436 covalentSet.add(value);
437 }
438 for (int k : ip11Other[otherIndex]) {
439 int value = openMMDualTopologyEnergy.mapToDualTopologyIndex(otherTopology, k);
440 covalentSet.add(value);
441 }
442 covalentMap.resize(0);
443 for (int k : covalentSet) {
444 covalentMap.append(k);
445 }
446 setCovalentMap(i, OpenMM_AmoebaMultipoleForce_PolarizationCovalent11, covalentMap);
447
448 }
449 covalentMap.destroy();
450 }
451
452
453
454
455
456
457
458 private double configureForce(ForceFieldEnergy forceFieldEnergy) {
459 ParticleMeshEwald pme = forceFieldEnergy.getPmeNode();
460 ForceField forceField = forceFieldEnergy.getMolecularAssembly().getForceField();
461
462 Polarization polarization = pme.getPolarizationType();
463 double doPolarization = 1.0;
464 SCFAlgorithm scfAlgorithm = null;
465 if (polarization != Polarization.MUTUAL) {
466 setPolarizationType(OpenMM_AmoebaMultipoleForce_Direct);
467 if (pme.getPolarizationType() == Polarization.NONE) {
468 doPolarization = 0.0;
469 }
470 } else {
471 String algorithm = forceField.getString("SCF_ALGORITHM", "CG");
472 try {
473 algorithm = algorithm.replaceAll("-", "_").toUpperCase();
474 scfAlgorithm = SCFAlgorithm.valueOf(algorithm);
475 } catch (Exception e) {
476 scfAlgorithm = SCFAlgorithm.CG;
477 }
478 if (scfAlgorithm == SCFAlgorithm.EPT) {
479
480
481
482
483
484 setPolarizationType(OpenMM_AmoebaMultipoleForce_Extrapolated);
485 DoubleArray exptCoefficients = new DoubleArray(4);
486 exptCoefficients.set(0, -0.154);
487 exptCoefficients.set(1, 0.017);
488 exptCoefficients.set(2, 0.657);
489 exptCoefficients.set(3, 0.475);
490 setExtrapolationCoefficients(exptCoefficients);
491 exptCoefficients.destroy();
492 } else {
493 setPolarizationType(OpenMM_AmoebaMultipoleForce_Mutual);
494 }
495 }
496
497 Crystal crystal = forceFieldEnergy.getCrystal();
498 double cutoff = pme.getEwaldCutoff();
499 double aewald = pme.getEwaldCoefficient();
500 if (!crystal.aperiodic()) {
501 setNonbondedMethod(OpenMM_AmoebaMultipoleForce_PME);
502 setCutoffDistance(cutoff * OpenMM_NmPerAngstrom);
503 double ewaldTolerance = 1.0e-04;
504 setEwaldErrorTolerance(ewaldTolerance);
505 ReciprocalSpace recip = pme.getReciprocalSpace();
506 int nx = recip.getXDim();
507 int ny = recip.getYDim();
508 int nz = recip.getZDim();
509 setPMEParameters(aewald / OpenMM_NmPerAngstrom, nx, ny, nz);
510 } else {
511 setNonbondedMethod(OpenMM_AmoebaMultipoleForce_NoCutoff);
512 }
513
514 setMutualInducedMaxIterations(500);
515 double poleps = pme.getPolarEps();
516 setMutualInducedTargetEpsilon(poleps);
517
518 AlchemicalParameters alchemicalParameters = pme.getAlchemicalParameters();
519 boolean lambdaTerm = pme.getLambdaTerm();
520
521 int forceGroup = forceField.getInteger("PME_FORCE_GROUP", 1);
522 setForceGroup(forceGroup);
523 if (logger.isLoggable(Level.INFO)) {
524 StringBuilder sb = new StringBuilder();
525 sb.append("\n Electrostatics\n");
526 sb.append(format(" Polarization: %8s\n", polarization.toString()));
527 if (polarization == Polarization.MUTUAL) {
528 sb.append(format(" SCF Convergence Criteria: %8.3e\n", poleps));
529 } else if (scfAlgorithm == SCFAlgorithm.EPT) {
530 sb.append(format(" SCF Algorithm: %8s\n", scfAlgorithm));
531 }
532 if (aewald > 0.0) {
533 sb.append(" Particle-mesh Ewald\n");
534 sb.append(format(" Ewald Coefficient: %8.3f\n", aewald));
535 sb.append(format(" Particle Cut-Off: %8.3f (A)", cutoff));
536 } else if (cutoff != Double.POSITIVE_INFINITY) {
537 sb.append(format(" Electrostatics Cut-Off: %8.3f (A)", cutoff));
538 } else {
539 sb.append(" Electrostatics Cut-Off: NONE");
540 }
541 logger.info(sb.toString());
542 if (lambdaTerm) {
543 sb = new StringBuilder(" Alchemical Parameters\n");
544 double permLambdaStart = alchemicalParameters.permLambdaStart;
545 double permLambdaEnd = alchemicalParameters.permLambdaEnd;
546 sb.append(format(" Permanent Multipole Range: %5.3f-%5.3f\n", permLambdaStart, permLambdaEnd));
547 double permLambdaAlpha = alchemicalParameters.permLambdaAlpha;
548 if (permLambdaAlpha != 0.0) {
549 logger.severe(" Permanent multipole softcore not supported for OpenMM.");
550 }
551 double permLambdaExponent = alchemicalParameters.permLambdaExponent;
552 sb.append(format(" Permanent Multipole Lambda Exponent: %5.3f\n", permLambdaExponent));
553 if (polarization != Polarization.NONE) {
554 double polLambdaExponent = alchemicalParameters.polLambdaExponent;
555 double polLambdaStart = alchemicalParameters.polLambdaStart;
556 double polLambdaEnd = alchemicalParameters.polLambdaEnd;
557 sb.append(format(" Polarization Lambda Exponent: %5.3f\n", polLambdaExponent));
558 sb.append(format(" Polarization Range: %5.3f-%5.3f\n", polLambdaStart, polLambdaEnd));
559 boolean doNoLigandCondensedSCF = alchemicalParameters.doNoLigandCondensedSCF;
560 if (doNoLigandCondensedSCF) {
561 logger.severe(" Condensed SCF without a ligand is not supported for OpenMM.");
562 }
563 }
564 boolean doLigandGKElec = alchemicalParameters.doLigandGKElec;
565 if (doLigandGKElec) {
566 logger.severe(" Isolated ligand electrostatics are not supported for OpenMM.");
567 }
568 logger.info(sb.toString());
569 }
570 logger.log(Level.FINE, format(" Force group:\t\t%d\n", forceGroup));
571 }
572 return doPolarization;
573 }
574
575
576
577
578
579
580
581 public static Force constructForce(OpenMMEnergy openMMEnergy) {
582 ParticleMeshEwald pme = openMMEnergy.getPmeNode();
583 if (pme == null) {
584 return null;
585 }
586 return new AmoebaMultipoleForce(openMMEnergy);
587 }
588
589
590
591
592
593
594
595
596 public static Force constructForce(int topology, OpenMMDualTopologyEnergy openMMDualTopologyEnergy) {
597 ForceFieldEnergy forceFieldEnergy = openMMDualTopologyEnergy.getForceFieldEnergy(topology);
598 ParticleMeshEwald pme = forceFieldEnergy.getPmeNode();
599 if (pme == null) {
600 return null;
601 }
602 return new AmoebaMultipoleForce(topology, openMMDualTopologyEnergy);
603 }
604
605
606
607
608
609
610
611 public void updateForce(Atom[] atoms, OpenMMEnergy openMMEnergy) {
612 ParticleMeshEwald pme = openMMEnergy.getPmeNode();
613 double quadrupoleConversion = OpenMM_NmPerAngstrom * OpenMM_NmPerAngstrom;
614 double polarityConversion = OpenMM_NmPerAngstrom * OpenMM_NmPerAngstrom * OpenMM_NmPerAngstrom;
615 double dampingFactorConversion = sqrt(OpenMM_NmPerAngstrom);
616
617 double doPolarization = 1.0;
618 if (pme.getPolarizationType() == Polarization.NONE) {
619 doPolarization = 0.0;
620 }
621
622 DoubleArray dipoles = new DoubleArray(3);
623 DoubleArray quadrupoles = new DoubleArray(9);
624
625 AlchemicalParameters alchemicalParameters = pme.getAlchemicalParameters();
626 boolean lambdaTerm = pme.getLambdaTerm();
627 if (lambdaTerm) {
628 AlchemicalParameters.AlchemicalMode alchemicalMode = alchemicalParameters.mode;
629
630 if (alchemicalMode != AlchemicalParameters.AlchemicalMode.SCALE) {
631 logger.severe(format(" Alchemical mode %s not supported for OpenMM.", alchemicalMode));
632 }
633
634 if (alchemicalParameters.permLambdaAlpha != 0.0) {
635 logger.severe(" Permanent multipole softcore not supported for OpenMM.");
636 }
637
638 if (alchemicalParameters.doLigandGKElec || alchemicalParameters.doLigandVaporElec) {
639 logger.severe(" Isolated ligand electrostatics are not supported for OpenMM.");
640 }
641
642 if (alchemicalParameters.doNoLigandCondensedSCF) {
643 logger.severe(" Condensed SCF without a ligand is not supported for OpenMM.");
644 }
645 }
646
647 double permLambda = alchemicalParameters.permLambda;
648 double polarLambda = alchemicalParameters.polLambda;
649
650 for (Atom atom : atoms) {
651 int index = atom.getXyzIndex() - 1;
652 MultipoleType multipoleType = pme.getMultipoleType(index);
653 PolarizeType polarizeType = pme.getPolarizeType(index);
654 int[] axisAtoms = atom.getAxisAtomIndices();
655
656 double permScale = 1.0;
657 double polarScale = doPolarization;
658
659 if (!atom.getUse() || !atom.getElectrostatics()) {
660 permScale = 0.0;
661 polarScale = 0.0;
662 }
663
664 if (atom.applyLambda()) {
665 permScale *= permLambda;
666 polarScale *= polarLambda;
667 }
668
669
670 int axisType = switch (multipoleType.frameDefinition) {
671 case NONE -> OpenMM_AmoebaMultipoleForce_NoAxisType;
672 case ZONLY -> OpenMM_AmoebaMultipoleForce_ZOnly;
673 case ZTHENX -> OpenMM_AmoebaMultipoleForce_ZThenX;
674 case BISECTOR -> OpenMM_AmoebaMultipoleForce_Bisector;
675 case ZTHENBISECTOR -> OpenMM_AmoebaMultipoleForce_ZBisect;
676 case THREEFOLD -> OpenMM_AmoebaMultipoleForce_ThreeFold;
677 };
678
679
680 for (int j = 0; j < 3; j++) {
681 dipoles.set(j, multipoleType.dipole[j] * OpenMM_NmPerAngstrom * permScale);
682 }
683 int l = 0;
684 for (int j = 0; j < 3; j++) {
685 for (int k = 0; k < 3; k++) {
686 quadrupoles.set(l++, multipoleType.quadrupole[j][k] * quadrupoleConversion / 3.0 * permScale);
687 }
688 }
689
690 int zaxis = 1;
691 int xaxis = 1;
692 int yaxis = 1;
693
694 if (axisAtoms != null) {
695 zaxis = axisAtoms[0];
696 if (axisAtoms.length > 1) {
697 xaxis = axisAtoms[1];
698 if (axisAtoms.length > 2) {
699 yaxis = axisAtoms[2];
700 }
701 }
702 } else {
703 axisType = OpenMM_AmoebaMultipoleForce_NoAxisType;
704 }
705
706
707 setMultipoleParameters(index, multipoleType.charge * permScale,
708 dipoles, quadrupoles, axisType, zaxis, xaxis, yaxis,
709 polarizeType.thole, polarizeType.pdamp * dampingFactorConversion,
710 polarizeType.polarizability * polarityConversion * polarScale);
711 }
712
713 dipoles.destroy();
714 quadrupoles.destroy();
715 updateParametersInContext(openMMEnergy.getContext());
716 }
717
718
719
720
721
722
723
724
725 public void updateForce(Atom[] atoms, int topology, OpenMMDualTopologyEnergy openMMDualTopologyEnergy) {
726 ForceFieldEnergy forceFieldEnergy = openMMDualTopologyEnergy.getForceFieldEnergy(topology);
727 ParticleMeshEwald pme = forceFieldEnergy.getPmeNode();
728 double quadrupoleConversion = OpenMM_NmPerAngstrom * OpenMM_NmPerAngstrom;
729 double polarityConversion = OpenMM_NmPerAngstrom * OpenMM_NmPerAngstrom * OpenMM_NmPerAngstrom;
730 double dampingFactorConversion = sqrt(OpenMM_NmPerAngstrom);
731
732 double doPolarization = 1.0;
733 if (pme.getPolarizationType() == Polarization.NONE) {
734 doPolarization = 0.0;
735 }
736
737
738 double scaleDT = Math.sqrt(openMMDualTopologyEnergy.getTopologyScale(topology));
739
740 DoubleArray dipoles = new DoubleArray(3);
741 DoubleArray quadrupoles = new DoubleArray(9);
742
743 AlchemicalParameters alchemicalParameters = pme.getAlchemicalParameters();
744 boolean lambdaTerm = pme.getLambdaTerm();
745 if (lambdaTerm) {
746 AlchemicalParameters.AlchemicalMode alchemicalMode = alchemicalParameters.mode;
747
748 if (alchemicalMode != AlchemicalParameters.AlchemicalMode.SCALE) {
749 logger.severe(format(" Alchemical mode %s not supported for OpenMM.", alchemicalMode));
750 }
751
752 if (alchemicalParameters.permLambdaAlpha != 0.0) {
753 logger.severe(" Permanent multipole softcore not supported for OpenMM.");
754 }
755
756 if (alchemicalParameters.doLigandGKElec || alchemicalParameters.doLigandVaporElec) {
757 logger.severe(" Isolated ligand electrostatics are not supported for OpenMM.");
758 }
759
760 if (alchemicalParameters.doNoLigandCondensedSCF) {
761 logger.severe(" Condensed SCF without a ligand is not supported for OpenMM.");
762 }
763 }
764
765 double permLambda = alchemicalParameters.permLambda;
766 double polarLambda = alchemicalParameters.polLambda;
767
768 for (Atom atom : atoms) {
769 if (atom.getTopologyIndex() != topology) {
770
771 continue;
772 }
773
774 int index = atom.getArrayIndex();
775 MultipoleType multipoleType = pme.getMultipoleType(index);
776 PolarizeType polarizeType = pme.getPolarizeType(index);
777 int[] axisAtoms = atom.getAxisAtomIndices();
778
779 double permScale = scaleDT;
780 double polarScale = doPolarization;
781 if (!atom.getUse() || !atom.getElectrostatics()) {
782 permScale = 0.0;
783 polarScale = 0.0;
784 }
785
786 if (atom.applyLambda()) {
787 permScale *= permLambda;
788 polarScale *= polarLambda;
789 }
790
791
792 int axisType = switch (multipoleType.frameDefinition) {
793 case NONE -> OpenMM_AmoebaMultipoleForce_NoAxisType;
794 case ZONLY -> OpenMM_AmoebaMultipoleForce_ZOnly;
795 case ZTHENX -> OpenMM_AmoebaMultipoleForce_ZThenX;
796 case BISECTOR -> OpenMM_AmoebaMultipoleForce_Bisector;
797 case ZTHENBISECTOR -> OpenMM_AmoebaMultipoleForce_ZBisect;
798 case THREEFOLD -> OpenMM_AmoebaMultipoleForce_ThreeFold;
799 };
800
801
802 for (int j = 0; j < 3; j++) {
803 dipoles.set(j, multipoleType.dipole[j] * OpenMM_NmPerAngstrom * permScale);
804 }
805 int l = 0;
806 for (int j = 0; j < 3; j++) {
807 for (int k = 0; k < 3; k++) {
808 quadrupoles.set(l++, multipoleType.quadrupole[j][k] * quadrupoleConversion / 3.0 * permScale);
809 }
810 }
811
812 int zaxis = 1;
813 int xaxis = 1;
814 int yaxis = 1;
815
816 if (axisAtoms != null) {
817 zaxis = axisAtoms[0];
818 zaxis = openMMDualTopologyEnergy.mapToDualTopologyIndex(topology, zaxis);
819 if (axisAtoms.length > 1) {
820 xaxis = axisAtoms[1];
821 xaxis = openMMDualTopologyEnergy.mapToDualTopologyIndex(topology, xaxis);
822 if (axisAtoms.length > 2) {
823 yaxis = axisAtoms[2];
824 yaxis = openMMDualTopologyEnergy.mapToDualTopologyIndex(topology, yaxis);
825 }
826 }
827 } else {
828 axisType = OpenMM_AmoebaMultipoleForce_NoAxisType;
829 }
830
831
832 int dualTopologyIndex = openMMDualTopologyEnergy.mapToDualTopologyIndex(topology, index);
833 setMultipoleParameters(dualTopologyIndex, multipoleType.charge * permScale,
834 dipoles, quadrupoles, axisType, zaxis, xaxis, yaxis,
835 polarizeType.thole, polarizeType.pdamp * dampingFactorConversion,
836 polarizeType.polarizability * polarityConversion * polarScale);
837 }
838
839 dipoles.destroy();
840 quadrupoles.destroy();
841 updateParametersInContext(openMMDualTopologyEnergy.getContext());
842 }
843
844 }