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 edu.uiowa.jopenmm.OpenMM_Vec3;
41 import ffx.crystal.Crystal;
42 import ffx.openmm.AndersenThermostat;
43 import ffx.openmm.CMMotionRemover;
44 import ffx.openmm.Force;
45 import ffx.openmm.MonteCarloBarostat;
46 import ffx.potential.MolecularAssembly;
47 import ffx.potential.bonded.Angle;
48 import ffx.potential.bonded.Atom;
49 import ffx.potential.bonded.Bond;
50 import ffx.potential.nonbonded.GeneralizedKirkwood;
51 import ffx.potential.nonbonded.ParticleMeshEwald;
52 import ffx.potential.nonbonded.VanDerWaals;
53 import ffx.potential.nonbonded.VanDerWaalsForm;
54 import ffx.potential.nonbonded.implicit.ChandlerCavitation;
55 import ffx.potential.nonbonded.implicit.DispersionRegion;
56 import ffx.potential.nonbonded.pme.AlchemicalParameters;
57 import ffx.potential.parameters.BondType;
58 import ffx.potential.parameters.ForceField;
59 import ffx.utilities.Constants;
60 import org.apache.commons.configuration2.CompositeConfiguration;
61 import org.checkerframework.checker.units.qual.A;
62
63 import javax.annotation.Nullable;
64 import java.util.logging.Logger;
65
66 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_NmPerAngstrom;
67 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_Boolean.OpenMM_True;
68 import static ffx.potential.parameters.VDWType.VDW_TYPE.LENNARD_JONES;
69 import static ffx.utilities.Constants.KCAL_TO_GRAM_ANG2_PER_PS2;
70 import static ffx.utilities.Constants.kB;
71 import static java.lang.String.format;
72 import static org.apache.commons.math3.util.FastMath.cos;
73 import static org.apache.commons.math3.util.FastMath.sqrt;
74 import static org.apache.commons.math3.util.FastMath.toRadians;
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90 public class OpenMMSystem extends ffx.openmm.System {
91
92 private static final Logger logger = Logger.getLogger(OpenMMSystem.class.getName());
93
94
95
96
97 private final OpenMMEnergy openMMEnergy;
98
99
100
101 private final ForceField forceField;
102
103
104
105 private final Atom[] atoms;
106
107
108
109
110 private boolean updateBondedTerms = false;
111
112
113
114 private BondForce bondForce = null;
115
116
117
118 private AngleForce angleForce = null;
119
120
121
122 private InPlaneAngleForce inPlaneAngleForce = null;
123
124
125
126 private StretchBendForce stretchBendForce = null;
127
128
129
130 private UreyBradleyForce ureyBradleyForce = null;
131
132
133
134 private OutOfPlaneBendForce outOfPlaneBendForce = null;
135
136
137
138 private PiOrbitalTorsionForce piOrbitalTorsionForce = null;
139
140
141
142 private TorsionForce torsionForce = null;
143
144
145
146 private ImproperTorsionForce improperTorsionForce = null;
147
148
149
150 private AmoebaTorsionTorsionForce amoebaTorsionTorsionForce = null;
151
152
153
154 private StretchTorsionForce stretchTorsionForce = null;
155
156
157
158 private AngleTorsionForce angleTorsionForce = null;
159
160
161
162 private RestrainTorsionsForce restrainTorsionsForce = null;
163
164
165
166 private RestrainPositionsForce restrainPositionsForce = null;
167
168
169
170 private RestrainGroupsForce restrainGroupsForce = null;
171
172
173
174 private AmoebaVdwForce amoebaVDWForce = null;
175
176
177
178 private AmoebaMultipoleForce amoebaMultipoleForce = null;
179
180
181
182 private AmoebaGeneralizedKirkwoodForce amoebaGeneralizedKirkwoodForce = null;
183
184
185
186 private AmoebaWcaDispersionForce amoebaWcaDispersionForce = null;
187
188
189
190 private AmoebaGKCavitationForce amoebaGKCavitationForce = null;
191
192
193
194 private FixedChargeGBForce fixedChargeGBForce = null;
195
196
197
198 private FixedChargeNonbondedForce fixedChargeNonBondedForce = null;
199
200
201
202 private FixedChargeAlchemicalForces fixedChargeAlchemicalForces = null;
203
204
205
206 private AndersenThermostat andersenThermostat = null;
207
208
209
210 private MonteCarloBarostat monteCarloBarostat = null;
211
212
213
214 private CMMotionRemover cmMotionRemover = null;
215
216
217
218 private boolean softcoreCreated = false;
219
220
221
222
223
224
225 public OpenMMSystem(OpenMMEnergy openMMEnergy) {
226 this.openMMEnergy = openMMEnergy;
227
228 MolecularAssembly molecularAssembly = openMMEnergy.getMolecularAssembly();
229 forceField = molecularAssembly.getForceField();
230 atoms = molecularAssembly.getAtomArray();
231
232
233 try {
234 addAtoms();
235 } catch (Exception e) {
236 logger.severe(" Atom without mass encountered.");
237 }
238
239 logger.info(format("\n OpenMM system created with %d atoms.", atoms.length));
240
241 }
242
243
244
245
246 public void addForces() {
247
248
249
250 boolean rigidHydrogen = forceField.getBoolean("RIGID_HYDROGEN", false);
251 boolean rigidBonds = forceField.getBoolean("RIGID_BONDS", false);
252 boolean rigidHydrogenAngles = forceField.getBoolean("RIGID_HYDROGEN_ANGLES", false);
253 if (rigidHydrogen) {
254 addHydrogenConstraints();
255 }
256 if (rigidBonds) {
257 addUpBondConstraints();
258 }
259 if (rigidHydrogenAngles) {
260 setUpHydrogenAngleConstraints();
261 }
262
263 logger.info("\n Bonded Terms\n");
264
265
266 if (rigidBonds) {
267 logger.info(" Not creating AmoebaBondForce because bonds are constrained.");
268 } else {
269 bondForce = (BondForce) BondForce.constructForce(openMMEnergy);
270 addForce(bondForce);
271 }
272
273
274 angleForce = (AngleForce) AngleForce.constructForce(openMMEnergy);
275 addForce(angleForce);
276
277
278 inPlaneAngleForce = (InPlaneAngleForce) InPlaneAngleForce.constructForce(openMMEnergy);
279 addForce(inPlaneAngleForce);
280
281
282 stretchBendForce = (StretchBendForce) StretchBendForce.constructForce(openMMEnergy);
283 addForce(stretchBendForce);
284
285
286 ureyBradleyForce = (UreyBradleyForce) UreyBradleyForce.constructForce(openMMEnergy);
287 addForce(ureyBradleyForce);
288
289
290 outOfPlaneBendForce = (OutOfPlaneBendForce) OutOfPlaneBendForce.constructForce(openMMEnergy);
291 addForce(outOfPlaneBendForce);
292
293
294 piOrbitalTorsionForce = (PiOrbitalTorsionForce) PiOrbitalTorsionForce.constructForce(openMMEnergy);
295 addForce(piOrbitalTorsionForce);
296
297
298 torsionForce = (TorsionForce) TorsionForce.constructForce(openMMEnergy);
299 addForce(torsionForce);
300
301
302 improperTorsionForce = (ImproperTorsionForce) ImproperTorsionForce.constructForce(openMMEnergy);
303 addForce(improperTorsionForce);
304
305
306 stretchTorsionForce = (StretchTorsionForce) StretchTorsionForce.constructForce(openMMEnergy);
307 addForce(stretchTorsionForce);
308
309
310 angleTorsionForce = (AngleTorsionForce) AngleTorsionForce.constructForce(openMMEnergy);
311 addForce(angleTorsionForce);
312
313
314 amoebaTorsionTorsionForce = (AmoebaTorsionTorsionForce) AmoebaTorsionTorsionForce.constructForce(openMMEnergy);
315 addForce(amoebaTorsionTorsionForce);
316
317
318 restrainPositionsForce = (RestrainPositionsForce) RestrainPositionsForce.constructForce(openMMEnergy);
319 addForce(restrainPositionsForce);
320
321
322 for (BondType.BondFunction function : BondType.BondFunction.values()) {
323 RestrainBondsForce restrainBondsForce = (RestrainBondsForce) RestrainBondsForce.constructForce(function, openMMEnergy);
324 addForce(restrainBondsForce);
325 }
326
327
328 restrainTorsionsForce = (RestrainTorsionsForce) RestrainTorsionsForce.constructForce(openMMEnergy);
329 addForce(restrainTorsionsForce);
330
331
332 restrainGroupsForce = (RestrainGroupsForce) RestrainGroupsForce.constructForce(openMMEnergy);
333 addForce(restrainGroupsForce);
334
335 setDefaultPeriodicBoxVectors();
336
337 VanDerWaals vdW = openMMEnergy.getVdwNode();
338 if (vdW != null) {
339 logger.info("\n Non-Bonded Terms");
340 VanDerWaalsForm vdwForm = vdW.getVDWForm();
341 if (vdwForm.vdwType == LENNARD_JONES) {
342 fixedChargeNonBondedForce = (FixedChargeNonbondedForce) FixedChargeNonbondedForce.constructForce(openMMEnergy);
343 addForce(fixedChargeNonBondedForce);
344 if (fixedChargeNonBondedForce != null) {
345 GeneralizedKirkwood gk = openMMEnergy.getGK();
346 if (gk != null) {
347 fixedChargeGBForce = (FixedChargeGBForce) FixedChargeGBForce.constructForce(openMMEnergy);
348 addForce(fixedChargeGBForce);
349 }
350 }
351 } else {
352
353 amoebaVDWForce = (AmoebaVdwForce) AmoebaVdwForce.constructForce(openMMEnergy);
354 addForce(amoebaVDWForce);
355
356
357 amoebaMultipoleForce = (AmoebaMultipoleForce) AmoebaMultipoleForce.constructForce(openMMEnergy);
358 addForce(amoebaMultipoleForce);
359
360
361 GeneralizedKirkwood gk = openMMEnergy.getGK();
362 if (gk != null) {
363 amoebaGeneralizedKirkwoodForce = (AmoebaGeneralizedKirkwoodForce) AmoebaGeneralizedKirkwoodForce.constructForce(openMMEnergy);
364 addForce(amoebaGeneralizedKirkwoodForce);
365
366
367 DispersionRegion dispersionRegion = gk.getDispersionRegion();
368 if (dispersionRegion != null) {
369 amoebaWcaDispersionForce = (AmoebaWcaDispersionForce) AmoebaWcaDispersionForce.constructForce(openMMEnergy);
370 addForce(amoebaWcaDispersionForce);
371 }
372
373
374 ChandlerCavitation chandlerCavitation = gk.getChandlerCavitation();
375 if (chandlerCavitation != null && chandlerCavitation.getGaussVol() != null) {
376 amoebaGKCavitationForce = (AmoebaGKCavitationForce) AmoebaGKCavitationForce.constructForce(openMMEnergy);
377 addForce(amoebaGKCavitationForce);
378 }
379 }
380 }
381 }
382 }
383
384
385
386
387
388 public void removeForce(Force force) {
389 if (force != null) {
390 removeForce(force.getForceIndex());
391 }
392 }
393
394
395
396
397
398
399 public void addAndersenThermostatForce(double targetTemp) {
400 double collisionFreq = forceField.getDouble("COLLISION_FREQ", 0.1);
401 addAndersenThermostatForce(targetTemp, collisionFreq);
402 }
403
404
405
406
407
408
409
410 public void addAndersenThermostatForce(double targetTemp, double collisionFreq) {
411 if (andersenThermostat != null) {
412 removeForce(andersenThermostat);
413 andersenThermostat = null;
414 }
415 andersenThermostat = new AndersenThermostat(targetTemp, collisionFreq);
416 addForce(andersenThermostat);
417 logger.info("\n Adding an Andersen thermostat");
418 logger.info(format(" Target Temperature: %6.2f (K)", targetTemp));
419 logger.info(format(" Collision Frequency: %6.2f (1/psec)", collisionFreq));
420 }
421
422
423
424
425 public void addCOMMRemoverForce() {
426 if (cmMotionRemover != null) {
427 removeForce(cmMotionRemover);
428 cmMotionRemover = null;
429 }
430 int frequency = 100;
431 cmMotionRemover = new CMMotionRemover(frequency);
432 addForce(cmMotionRemover);
433 logger.info("\n Added a Center of Mass Motion Remover");
434 logger.info(format(" Frequency: %6d", frequency));
435 }
436
437
438
439
440
441
442
443
444 public void addMonteCarloBarostatForce(double targetPressure, double targetTemp, int frequency) {
445 if (monteCarloBarostat != null) {
446 removeForce(monteCarloBarostat);
447 monteCarloBarostat = null;
448 }
449 double pressureInBar = targetPressure * Constants.ATM_TO_BAR;
450 monteCarloBarostat = new MonteCarloBarostat(pressureInBar, targetTemp, frequency);
451 CompositeConfiguration properties = openMMEnergy.getMolecularAssembly().getProperties();
452 if (properties.containsKey("barostat-seed")) {
453 int randomSeed = properties.getInt("barostat-seed", 0);
454 logger.info(format(" Setting random seed %d for Monte Carlo Barostat", randomSeed));
455 monteCarloBarostat.setRandomNumberSeed(randomSeed);
456 }
457 addForce(monteCarloBarostat);
458 logger.info("\n Added a Monte Carlo Barostat");
459 logger.info(format(" Target Pressure: %6.2f (atm)", targetPressure));
460 logger.info(format(" Target Temperature: %6.2f (K)", targetTemp));
461 logger.info(format(" MC Move Frequency: %6d", frequency));
462 }
463
464
465
466
467
468
469 public int calculateDegreesOfFreedom() {
470
471 int dof = openMMEnergy.getNumberOfVariables();
472
473 dof = dof - getNumConstraints();
474
475 if (cmMotionRemover != null) {
476 dof -= 3;
477 }
478 return dof;
479 }
480
481 public double getTemperature(double kineticEnergy) {
482 double dof = calculateDegreesOfFreedom();
483 return 2.0 * kineticEnergy * KCAL_TO_GRAM_ANG2_PER_PS2 / (kB * dof);
484 }
485
486
487
488
489
490
491 public ForceField getForceField() {
492 return forceField;
493 }
494
495
496
497
498
499
500 public Crystal getCrystal() {
501 return openMMEnergy.getCrystal();
502 }
503
504
505
506
507
508
509 public int getNumberOfVariables() {
510 return openMMEnergy.getNumberOfVariables();
511 }
512
513
514
515
516 public void free() {
517 if (getPointer() != null) {
518 logger.fine(" Free OpenMM system.");
519 destroy();
520 logger.fine(" Free OpenMM system completed.");
521 }
522 }
523
524 public void setUpdateBondedTerms(boolean updateBondedTerms) {
525 this.updateBondedTerms = updateBondedTerms;
526 }
527
528
529
530
531
532
533
534
535
536
537
538 private void setDefaultPeriodicBoxVectors() {
539 Crystal crystal = openMMEnergy.getCrystal();
540 if (!crystal.aperiodic()) {
541 OpenMM_Vec3 a = new OpenMM_Vec3();
542 OpenMM_Vec3 b = new OpenMM_Vec3();
543 OpenMM_Vec3 c = new OpenMM_Vec3();
544 double[][] Ai = crystal.Ai;
545 a.x = Ai[0][0] * OpenMM_NmPerAngstrom;
546 a.y = Ai[0][1] * OpenMM_NmPerAngstrom;
547 a.z = Ai[0][2] * OpenMM_NmPerAngstrom;
548 b.x = Ai[1][0] * OpenMM_NmPerAngstrom;
549 b.y = Ai[1][1] * OpenMM_NmPerAngstrom;
550 b.z = Ai[1][2] * OpenMM_NmPerAngstrom;
551 c.x = Ai[2][0] * OpenMM_NmPerAngstrom;
552 c.y = Ai[2][1] * OpenMM_NmPerAngstrom;
553 c.z = Ai[2][2] * OpenMM_NmPerAngstrom;
554 setDefaultPeriodicBoxVectors(a, b, c);
555 }
556 }
557
558
559
560
561
562
563 public void updateParameters(@Nullable Atom[] atoms) {
564 VanDerWaals vanDerWaals = openMMEnergy.getVdwNode();
565 if (vanDerWaals != null) {
566 boolean vdwLambdaTerm = vanDerWaals.getLambdaTerm();
567 if (vdwLambdaTerm) {
568 double lambdaVDW = vanDerWaals.getLambda();
569 if (fixedChargeNonBondedForce != null) {
570 if (!softcoreCreated) {
571 fixedChargeAlchemicalForces = new FixedChargeAlchemicalForces(openMMEnergy, fixedChargeNonBondedForce);
572 addForce(fixedChargeAlchemicalForces.getFixedChargeSoftcoreForce());
573 addForce(fixedChargeAlchemicalForces.getAlchemicalAlchemicalStericsForce());
574 addForce(fixedChargeAlchemicalForces.getNonAlchemicalAlchemicalStericsForce());
575
576 openMMEnergy.getContext().reinitialize(OpenMM_True);
577 softcoreCreated = true;
578 }
579
580 openMMEnergy.getContext().setParameter("vdw_lambda", lambdaVDW);
581 } else if (amoebaVDWForce != null) {
582
583 openMMEnergy.getContext().setParameter("AmoebaVdwLambda", lambdaVDW);
584 if (softcoreCreated) {
585 ParticleMeshEwald pme = openMMEnergy.getPmeNode();
586
587 if (pme == null || !pme.getLambdaTerm()) {
588 return;
589 }
590 } else {
591 softcoreCreated = true;
592 }
593 }
594 }
595 }
596
597
598
599 if (updateBondedTerms) {
600 if (bondForce != null) {
601 bondForce.updateForce(openMMEnergy);
602 }
603 if (angleForce != null) {
604 angleForce.updateForce(openMMEnergy);
605 }
606 if (inPlaneAngleForce != null) {
607 inPlaneAngleForce.updateForce(openMMEnergy);
608 }
609 if (stretchBendForce != null) {
610 stretchBendForce.updateForce(openMMEnergy);
611 }
612 if (ureyBradleyForce != null) {
613 ureyBradleyForce.updateForce(openMMEnergy);
614 }
615 if (outOfPlaneBendForce != null) {
616 outOfPlaneBendForce.updateForce(openMMEnergy);
617 }
618 if (piOrbitalTorsionForce != null) {
619 piOrbitalTorsionForce.updateForce(openMMEnergy);
620 }
621 if (torsionForce != null) {
622 torsionForce.updateForce(openMMEnergy);
623 }
624 if (improperTorsionForce != null) {
625 improperTorsionForce.updateForce(openMMEnergy);
626 }
627 }
628
629 if (restrainTorsionsForce != null) {
630 restrainTorsionsForce.updateForce(openMMEnergy);
631 }
632
633 if (atoms == null || atoms.length == 0) {
634 return;
635 }
636
637
638 if (fixedChargeNonBondedForce != null) {
639 fixedChargeNonBondedForce.updateForce(atoms, openMMEnergy);
640 }
641
642
643 if (fixedChargeGBForce != null) {
644 fixedChargeGBForce.updateForce(atoms, openMMEnergy);
645 }
646
647
648 if (amoebaVDWForce != null) {
649 amoebaVDWForce.updateForce(atoms, openMMEnergy);
650 }
651
652
653 if (amoebaMultipoleForce != null) {
654 amoebaMultipoleForce.updateForce(atoms, openMMEnergy);
655 }
656
657
658 if (amoebaGeneralizedKirkwoodForce != null) {
659 amoebaGeneralizedKirkwoodForce.updateForce(atoms, openMMEnergy);
660 }
661
662
663 if (amoebaWcaDispersionForce != null) {
664 amoebaWcaDispersionForce.updateForce(atoms, openMMEnergy);
665 }
666
667
668 if (amoebaGKCavitationForce != null) {
669 amoebaGKCavitationForce.updateForce(atoms, openMMEnergy);
670 }
671 }
672
673
674
675
676
677 private void addAtoms() throws Exception {
678 for (Atom atom : atoms) {
679 double mass = atom.getMass();
680 if (mass < 0.0) {
681 throw new Exception(" Atom with mass less than 0.");
682 }
683 if (mass == 0.0) {
684 logger.info(format(" Atom %s has zero mass.", atom));
685 }
686 addParticle(mass);
687 }
688 }
689
690
691
692
693 public void updateAtomMass() {
694 int index = 0;
695 for (Atom atom : atoms) {
696 double mass = 0.0;
697 if (atom.isActive()) {
698 mass = atom.getMass();
699 }
700 setParticleMass(index++, mass);
701 }
702 }
703
704 public boolean hasAmoebaCavitationForce() {
705 return amoebaGKCavitationForce != null;
706 }
707
708
709
710
711 private void addUpBondConstraints() {
712 Bond[] bonds = openMMEnergy.getBonds();
713 if (bonds == null || bonds.length < 1) {
714 return;
715 }
716
717 logger.info(" Adding constraints for all bonds.");
718 for (Bond bond : bonds) {
719 Atom atom1 = bond.getAtom(0);
720 Atom atom2 = bond.getAtom(1);
721 int iAtom1 = atom1.getXyzIndex() - 1;
722 int iAtom2 = atom2.getXyzIndex() - 1;
723 addConstraint(iAtom1, iAtom2, bond.bondType.distance * OpenMM_NmPerAngstrom);
724 }
725 }
726
727
728
729
730 private void addHydrogenConstraints() {
731 Bond[] bonds = openMMEnergy.getBonds();
732 if (bonds == null || bonds.length < 1) {
733 return;
734 }
735
736 logger.info(" Adding constraints for hydrogen bonds.");
737 for (Bond bond : bonds) {
738 Atom atom1 = bond.getAtom(0);
739 Atom atom2 = bond.getAtom(1);
740 if (atom1.isHydrogen() || atom2.isHydrogen()) {
741 BondType bondType = bond.bondType;
742 int iAtom1 = atom1.getXyzIndex() - 1;
743 int iAtom2 = atom2.getXyzIndex() - 1;
744 addConstraint(iAtom1, iAtom2, bondType.distance * OpenMM_NmPerAngstrom);
745 }
746 }
747 }
748
749
750
751
752 private void setUpHydrogenAngleConstraints() {
753 Angle[] angles = openMMEnergy.getAngles();
754 if (angles == null || angles.length < 1) {
755 return;
756 }
757
758 logger.info(" Adding hydrogen angle constraints.");
759
760 for (Angle angle : angles) {
761 if (isHydrogenAngle(angle)) {
762 Atom atom1 = angle.getAtom(0);
763 Atom atom3 = angle.getAtom(2);
764
765
766
767 Bond bond1 = angle.getBond(0);
768 double distance1 = bond1.bondType.distance;
769
770 Bond bond2 = angle.getBond(1);
771 double distance2 = bond2.bondType.distance;
772
773
774 double angleVal = angle.angleType.angle[angle.nh];
775
776
777 double falseBondLength = sqrt(
778 distance1 * distance1 + distance2 * distance2 - 2.0 * distance1 * distance2 * cos(
779 toRadians(angleVal)));
780
781 int iAtom1 = atom1.getXyzIndex() - 1;
782 int iAtom3 = atom3.getXyzIndex() - 1;
783 addConstraint(iAtom1, iAtom3, falseBondLength * OpenMM_NmPerAngstrom);
784 }
785 }
786 }
787
788
789
790
791
792
793
794
795 private boolean isHydrogenAngle(Angle angle) {
796 if (angle.containsHydrogen()) {
797
798 double angleVal = angle.angleType.angle[angle.nh];
799
800
801
802 if (angleVal < 160.0) {
803 Atom atom1 = angle.getAtom(0);
804 Atom atom2 = angle.getAtom(1);
805 Atom atom3 = angle.getAtom(2);
806
807
808 return atom1.isHydrogen() && atom3.isHydrogen() && !atom2.isHydrogen();
809 }
810 }
811 return false;
812 }
813
814 }