View Javadoc
1   // ******************************************************************************
2   //
3   // Title:       Force Field X.
4   // Description: Force Field X - Software for Molecular Biophysics.
5   // Copyright:   Copyright (c) Michael J. Schnieders 2001-2025.
6   //
7   // This file is part of Force Field X.
8   //
9   // Force Field X is free software; you can redistribute it and/or modify it
10  // under the terms of the GNU General Public License version 3 as published by
11  // the Free Software Foundation.
12  //
13  // Force Field X is distributed in the hope that it will be useful, but WITHOUT
14  // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
15  // FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
16  // details.
17  //
18  // You should have received a copy of the GNU General Public License along with
19  // Force Field X; if not, write to the Free Software Foundation, Inc., 59 Temple
20  // Place, Suite 330, Boston, MA 02111-1307 USA
21  //
22  // Linking this library statically or dynamically with other modules is making a
23  // combined work based on this library. Thus, the terms and conditions of the
24  // GNU General Public License cover the whole combination.
25  //
26  // As a special exception, the copyright holders of this library give you
27  // permission to link this library with independent modules to produce an
28  // executable, regardless of the license terms of these independent modules, and
29  // to copy and distribute the resulting executable under terms of your choice,
30  // provided that you also meet, for each linked independent module, the terms
31  // and conditions of the license of that module. An independent module is a
32  // module which is not derived from or based on this library. If you modify this
33  // library, you may extend this exception to your version of the library, but
34  // you are not obligated to do so. If you do not wish to do so, delete this
35  // exception statement from your version.
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   * Create and manage an OpenMM System.
78   *
79   * <p>The definition of a System involves four elements:
80   *
81   * <p>The particles and constraints are defined directly by the System object, while forces are
82   * defined by objects that extend the Force class. After creating a System, call addParticle() once
83   * for each particle, addConstraint() for each constraint, and addForce() for each Force.
84   *
85   * <p>In addition, particles may be designated as "virtual sites". These are particles whose
86   * positions are computed automatically based on the positions of other particles. To define a
87   * virtual site, call setVirtualSite(), passing in a VirtualSite object that defines the rules for
88   * computing its position.
89   */
90  public class OpenMMSystem extends ffx.openmm.System {
91  
92    private static final Logger logger = Logger.getLogger(OpenMMSystem.class.getName());
93  
94    /**
95     * The ForceFieldEnergyOpenMM instance.
96     */
97    private final OpenMMEnergy openMMEnergy;
98    /**
99     * The Force Field in use.
100    */
101   private final ForceField forceField;
102   /**
103    * Array of atoms in the system.
104    */
105   private final Atom[] atoms;
106   /**
107    * This flag indicates bonded force constants and equilibria are updated (e.g. during ManyBody
108    * titration).
109    */
110   private boolean updateBondedTerms = false;
111   /**
112    * OpenMM Custom Bond Force
113    */
114   private BondForce bondForce = null;
115   /**
116    * OpenMM Custom Angle Force
117    */
118   private AngleForce angleForce = null;
119   /**
120    * OpenMM Custom In-Plane Angle Force
121    */
122   private InPlaneAngleForce inPlaneAngleForce = null;
123   /**
124    * OpenMM Custom Stretch-Bend Force
125    */
126   private StretchBendForce stretchBendForce = null;
127   /**
128    * OpenMM Custom Urey-Bradley Force
129    */
130   private UreyBradleyForce ureyBradleyForce = null;
131   /**
132    * OpenMM Custom Out-of-Plane Bend Force
133    */
134   private OutOfPlaneBendForce outOfPlaneBendForce = null;
135   /**
136    * OpenMM Custom Pi-Torsion Force
137    */
138   private PiOrbitalTorsionForce piOrbitalTorsionForce = null;
139   /**
140    * OpenMM AMOEBA Torsion Force.
141    */
142   private TorsionForce torsionForce = null;
143   /**
144    * OpenMM Improper Torsion Force.
145    */
146   private ImproperTorsionForce improperTorsionForce = null;
147   /**
148    * OpenMM Torsion-Torsion Force.
149    */
150   private AmoebaTorsionTorsionForce amoebaTorsionTorsionForce = null;
151   /**
152    * OpenMM Stretch-Torsion Force.
153    */
154   private StretchTorsionForce stretchTorsionForce = null;
155   /**
156    * OpenMM Angle-Torsion Force.
157    */
158   private AngleTorsionForce angleTorsionForce = null;
159   /**
160    * OpenMM Restraint-Torsion Force.
161    */
162   private RestrainTorsionsForce restrainTorsionsForce = null;
163   /**
164    * OpenMM Restrain-Position Force.
165    */
166   private RestrainPositionsForce restrainPositionsForce = null;
167   /**
168    * OpenMM Restrain-Groups Force.
169    */
170   private RestrainGroupsForce restrainGroupsForce = null;
171   /**
172    * OpenMM AMOEBA van der Waals Force.
173    */
174   private AmoebaVdwForce amoebaVDWForce = null;
175   /**
176    * OpenMM AMOEBA Multipole Force.
177    */
178   private AmoebaMultipoleForce amoebaMultipoleForce = null;
179   /**
180    * OpenMM Generalized Kirkwood Force.
181    */
182   private AmoebaGeneralizedKirkwoodForce amoebaGeneralizedKirkwoodForce = null;
183   /**
184    * OpenMM AMOEBA WCA Dispersion Force.
185    */
186   private AmoebaWcaDispersionForce amoebaWcaDispersionForce = null;
187   /**
188    * OpenMM AMOEBA WCA Cavitation Force.
189    */
190   private AmoebaGKCavitationForce amoebaGKCavitationForce = null;
191   /**
192    * OpenMM Custom GB Force.
193    */
194   private FixedChargeGBForce fixedChargeGBForce = null;
195   /**
196    * OpenMM Fixed Charge Non-Bonded Force.
197    */
198   private FixedChargeNonbondedForce fixedChargeNonBondedForce = null;
199   /**
200    * Custom forces to handle alchemical transformations for fixed charge systems.
201    */
202   private FixedChargeAlchemicalForces fixedChargeAlchemicalForces = null;
203   /**
204    * OpenMM thermostat. Currently, an Andersen thermostat is supported.
205    */
206   private AndersenThermostat andersenThermostat = null;
207   /**
208    * Barostat to be added if NPT (isothermal-isobaric) dynamics is requested.
209    */
210   private MonteCarloBarostat monteCarloBarostat = null;
211   /**
212    * OpenMM center-of-mass motion remover.
213    */
214   private CMMotionRemover cmMotionRemover = null;
215   /**
216    * Fixed charge softcore vdW force boolean.
217    */
218   private boolean softcoreCreated = false;
219 
220   /**
221    * OpenMMSystem constructor.
222    *
223    * @param openMMEnergy ForceFieldEnergyOpenMM instance.
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     // Load atoms.
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    * Add forces to the system.
245    */
246   public void addForces() {
247 
248     // Set up rigid constraints.
249     // These flags need to be set before bonds and angles are created below.
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     // Add Bond Force.
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     // Add Angle Force.
274     angleForce = (AngleForce) AngleForce.constructForce(openMMEnergy);
275     addForce(angleForce);
276 
277     // Add In-Plane Angle Force.
278     inPlaneAngleForce = (InPlaneAngleForce) InPlaneAngleForce.constructForce(openMMEnergy);
279     addForce(inPlaneAngleForce);
280 
281     // Add Stretch-Bend Force.
282     stretchBendForce = (StretchBendForce) StretchBendForce.constructForce(openMMEnergy);
283     addForce(stretchBendForce);
284 
285     // Add Urey-Bradley Force.
286     ureyBradleyForce = (UreyBradleyForce) UreyBradleyForce.constructForce(openMMEnergy);
287     addForce(ureyBradleyForce);
288 
289     // Out-of Plane Bend Force.
290     outOfPlaneBendForce = (OutOfPlaneBendForce) OutOfPlaneBendForce.constructForce(openMMEnergy);
291     addForce(outOfPlaneBendForce);
292 
293     // Add Pi-Torsion Force.
294     piOrbitalTorsionForce = (PiOrbitalTorsionForce) PiOrbitalTorsionForce.constructForce(openMMEnergy);
295     addForce(piOrbitalTorsionForce);
296 
297     // Add Torsion Force.
298     torsionForce = (TorsionForce) TorsionForce.constructForce(openMMEnergy);
299     addForce(torsionForce);
300 
301     // Add Improper Torsion Force.
302     improperTorsionForce = (ImproperTorsionForce) ImproperTorsionForce.constructForce(openMMEnergy);
303     addForce(improperTorsionForce);
304 
305     // Add Stretch-Torsion coupling terms.
306     stretchTorsionForce = (StretchTorsionForce) StretchTorsionForce.constructForce(openMMEnergy);
307     addForce(stretchTorsionForce);
308 
309     // Add Angle-Torsion coupling terms.
310     angleTorsionForce = (AngleTorsionForce) AngleTorsionForce.constructForce(openMMEnergy);
311     addForce(angleTorsionForce);
312 
313     // Add Torsion-Torsion Force.
314     amoebaTorsionTorsionForce = (AmoebaTorsionTorsionForce) AmoebaTorsionTorsionForce.constructForce(openMMEnergy);
315     addForce(amoebaTorsionTorsionForce);
316 
317     // Add Restrain-Position force.
318     restrainPositionsForce = (RestrainPositionsForce) RestrainPositionsForce.constructForce(openMMEnergy);
319     addForce(restrainPositionsForce);
320 
321     // Add a Restrain-Bond force for each functional form.
322     for (BondType.BondFunction function : BondType.BondFunction.values()) {
323       RestrainBondsForce restrainBondsForce = (RestrainBondsForce) RestrainBondsForce.constructForce(function, openMMEnergy);
324       addForce(restrainBondsForce);
325     }
326 
327     // Add Restraint-Torsions
328     restrainTorsionsForce = (RestrainTorsionsForce) RestrainTorsionsForce.constructForce(openMMEnergy);
329     addForce(restrainTorsionsForce);
330 
331     // Add Restrain-Groups force.
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         // Add vdW Force.
353         amoebaVDWForce = (AmoebaVdwForce) AmoebaVdwForce.constructForce(openMMEnergy);
354         addForce(amoebaVDWForce);
355 
356         // Add Multipole Force.
357         amoebaMultipoleForce = (AmoebaMultipoleForce) AmoebaMultipoleForce.constructForce(openMMEnergy);
358         addForce(amoebaMultipoleForce);
359 
360         // Add Generalized Kirkwood Force.
361         GeneralizedKirkwood gk = openMMEnergy.getGK();
362         if (gk != null) {
363           amoebaGeneralizedKirkwoodForce = (AmoebaGeneralizedKirkwoodForce) AmoebaGeneralizedKirkwoodForce.constructForce(openMMEnergy);
364           addForce(amoebaGeneralizedKirkwoodForce);
365 
366           // Add WCA Dispersion Force.
367           DispersionRegion dispersionRegion = gk.getDispersionRegion();
368           if (dispersionRegion != null) {
369             amoebaWcaDispersionForce = (AmoebaWcaDispersionForce) AmoebaWcaDispersionForce.constructForce(openMMEnergy);
370             addForce(amoebaWcaDispersionForce);
371           }
372 
373           // Add a GaussVol Cavitation Force.
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    * Remove a force from the OpenMM System.
386    * The OpenMM memory associated with the removed Force object is deleted.
387    */
388   public void removeForce(Force force) {
389     if (force != null) {
390       removeForce(force.getForceIndex());
391     }
392   }
393 
394   /**
395    * Add an Andersen thermostat to the system.
396    *
397    * @param targetTemp Target temperature in Kelvins.
398    */
399   public void addAndersenThermostatForce(double targetTemp) {
400     double collisionFreq = forceField.getDouble("COLLISION_FREQ", 0.1);
401     addAndersenThermostatForce(targetTemp, collisionFreq);
402   }
403 
404   /**
405    * Add an Andersen thermostat to the system.
406    *
407    * @param targetTemp    Target temperature in Kelvins.
408    * @param collisionFreq Collision frequency in 1/psec.
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    * Adds a force that removes center-of-mass motion.
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    * Add a Monte Carlo Barostat to the system.
439    *
440    * @param targetPressure The target pressure (in atm).
441    * @param targetTemp     The target temperature.
442    * @param frequency      The frequency to apply the barostat.
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    * Calculate the number of degrees of freedom.
466    *
467    * @return Number of degrees of freedom.
468    */
469   public int calculateDegreesOfFreedom() {
470     // Begin from the 3 times the number of active atoms.
471     int dof = openMMEnergy.getNumberOfVariables();
472     // Remove OpenMM constraints.
473     dof = dof - getNumConstraints();
474     // Remove center of mass motion.
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    * Get the ForceField in use.
488    *
489    * @return ForceField instance.
490    */
491   public ForceField getForceField() {
492     return forceField;
493   }
494 
495   /**
496    * Get the Crystal instance.
497    *
498    * @return the Crystal instance.
499    */
500   public Crystal getCrystal() {
501     return openMMEnergy.getCrystal();
502   }
503 
504   /**
505    * Get the number of variables.
506    *
507    * @return the number of variables.
508    */
509   public int getNumberOfVariables() {
510     return openMMEnergy.getNumberOfVariables();
511   }
512 
513   /**
514    * Destroy the system.
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    * Set the default values of the vectors defining the axes of the periodic box (measured in nm).
530    *
531    * <p>Any newly created Context will have its box vectors set to these. They will affect any
532    * Force added to the System that uses periodic boundary conditions.
533    *
534    * <p>Triclinic boxes are supported, but the vectors must satisfy certain requirements. In
535    * particular, a must point in the x direction, b must point "mostly" in the y direction, and c
536    * must point "mostly" in the z direction. See the documentation for details.
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    * Update parameters if the Use flags and/or Lambda value has changed.
560    *
561    * @param atoms Atoms in this list are considered.
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             // Re-initialize the context.
576             openMMEnergy.getContext().reinitialize(OpenMM_True);
577             softcoreCreated = true;
578           }
579           // Update the lambda value.
580           openMMEnergy.getContext().setParameter("vdw_lambda", lambdaVDW);
581         } else if (amoebaVDWForce != null) {
582           // Update the lambda value.
583           openMMEnergy.getContext().setParameter("AmoebaVdwLambda", lambdaVDW);
584           if (softcoreCreated) {
585             ParticleMeshEwald pme = openMMEnergy.getPmeNode();
586             // Avoid any updateParametersInContext calls if vdwLambdaTerm is true, but not other alchemical terms.
587             if (pme == null || !pme.getLambdaTerm()) {
588               return;
589             }
590           } else {
591             softcoreCreated = true;
592           }
593         }
594       }
595     }
596 
597     // Note Stretch-Torsion and Angle-Torsion terms (for nucleic acids)
598     // and Torsion-Torsion terms (for protein backbones) are not updated yet.
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     // Update fixed charge non-bonded parameters.
638     if (fixedChargeNonBondedForce != null) {
639       fixedChargeNonBondedForce.updateForce(atoms, openMMEnergy);
640     }
641 
642     // Update fixed charge GB parameters.
643     if (fixedChargeGBForce != null) {
644       fixedChargeGBForce.updateForce(atoms, openMMEnergy);
645     }
646 
647     // Update AMOEBA vdW parameters.
648     if (amoebaVDWForce != null) {
649       amoebaVDWForce.updateForce(atoms, openMMEnergy);
650     }
651 
652     // Update AMOEBA polarizable multipole parameters.
653     if (amoebaMultipoleForce != null) {
654       amoebaMultipoleForce.updateForce(atoms, openMMEnergy);
655     }
656 
657     // Update GK force.
658     if (amoebaGeneralizedKirkwoodForce != null) {
659       amoebaGeneralizedKirkwoodForce.updateForce(atoms, openMMEnergy);
660     }
661 
662     // Update WCA Force.
663     if (amoebaWcaDispersionForce != null) {
664       amoebaWcaDispersionForce.updateForce(atoms, openMMEnergy);
665     }
666 
667     // Update WCA Force.
668     if (amoebaGKCavitationForce != null) {
669       amoebaGKCavitationForce.updateForce(atoms, openMMEnergy);
670     }
671   }
672 
673   /**
674    * Adds atoms from the molecular assembly to the OpenMM System and reports to the user the number
675    * of particles added.
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    * This method sets the mass of inactive atoms to zero.
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    * Add a constraint to every bond.
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    * Add a constraint to every bond that includes a hydrogen atom.
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    * Add a constraint to every angle that includes two hydrogen atoms.
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         // Calculate a "false bond" length between atoms 1 and 3 to constrain the angle using the
766         // law of cosines.
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         // Equilibrium angle value in degrees.
774         double angleVal = angle.angleType.angle[angle.nh];
775 
776         // Law of cosines.
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    * Check to see if an angle is a hydrogen angle. This method only returns true for hydrogen
790    * angles that are less than 160 degrees.
791    *
792    * @param angle Angle to check.
793    * @return boolean indicating whether an angle is a hydrogen angle that is less than 160 degrees.
794    */
795   private boolean isHydrogenAngle(Angle angle) {
796     if (angle.containsHydrogen()) {
797       // Equilibrium angle value in degrees.
798       double angleVal = angle.angleType.angle[angle.nh];
799       // Make sure angle is less than 160 degrees because greater than 160 degrees will not be
800       // constrained
801       // well using the law of cosines.
802       if (angleVal < 160.0) {
803         Atom atom1 = angle.getAtom(0);
804         Atom atom2 = angle.getAtom(1);
805         Atom atom3 = angle.getAtom(2);
806         // Setting constraints only on angles where atom1 or atom3 is a hydrogen while atom2 is
807         // not a hydrogen.
808         return atom1.isHydrogen() && atom3.isHydrogen() && !atom2.isHydrogen();
809       }
810     }
811     return false;
812   }
813 
814 }