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 ffx.openmm.CustomAngleForce;
41  import ffx.openmm.DoubleArray;
42  import ffx.openmm.Force;
43  import ffx.potential.ForceFieldEnergy;
44  import ffx.potential.bonded.Angle;
45  import ffx.potential.bonded.Atom;
46  import ffx.potential.parameters.AngleType;
47  import ffx.potential.parameters.ForceField;
48  import ffx.potential.terms.AnglePotentialEnergy;
49  
50  import java.util.logging.Level;
51  import java.util.logging.Logger;
52  
53  import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_KJPerKcal;
54  import static java.lang.String.format;
55  
56  /**
57   * OpenMM Angle Force.
58   */
59  public class AngleForce extends CustomAngleForce {
60  
61    private static final Logger logger = Logger.getLogger(AngleForce.class.getName());
62  
63    private int nAngles = 0;
64    private final boolean manyBodyTitration;
65    private final boolean rigidHydrogenAngles;
66  
67    /**
68     * Create an OpenMM Angle Force.
69     *
70     * @param openMMEnergy The OpenMM Energy instance that contains the angles.
71     */
72    public AngleForce(AnglePotentialEnergy anglePotentialEnergy, OpenMMEnergy openMMEnergy) {
73      super(anglePotentialEnergy.getAngleEnergyString());
74      ForceField forceField = openMMEnergy.getMolecularAssembly().getForceField();
75      manyBodyTitration = forceField.getBoolean("MANYBODY_TITRATION", false);
76      rigidHydrogenAngles = forceField.getBoolean("RIGID_HYDROGEN_ANGLES", false);
77      addPerAngleParameter("theta0");
78      addPerAngleParameter("k");
79      setName("Angle");
80  
81      DoubleArray parameters = new DoubleArray(0);
82      Angle[] angles = anglePotentialEnergy.getAngleArray();
83      for (Angle angle : angles) {
84        AngleType angleType = angle.getAngleType();
85        AngleType.AngleMode angleMode = angleType.angleMode;
86        if (!manyBodyTitration && angleMode == AngleType.AngleMode.IN_PLANE) {
87          // Skip In-Plane angles unless this is ManyBody Titration.
88        } else if (isHydrogenAngle(angle) && rigidHydrogenAngles) {
89          logger.log(Level.INFO, " Constrained angle %s was not added the AngleForce.", angle);
90        } else {
91          int i1 = angle.getAtom(0).getArrayIndex();
92          int i2 = angle.getAtom(1).getArrayIndex();
93          int i3 = angle.getAtom(2).getArrayIndex();
94  
95          double theta0 = angleType.angle[angle.nh];
96          double k = OpenMM_KJPerKcal * angleType.angleUnit * angleType.forceConstant;
97          if (angleMode == AngleType.AngleMode.IN_PLANE) {
98            // This is a place-holder Angle, in case the In-Plane Angle is swtiched to a
99            // Normal Angle during in the udpateAngleForce.
100           k = 0.0;
101         }
102         parameters.append(theta0);
103         parameters.append(k);
104         addAngle(i1, i2, i3, parameters);
105         nAngles++;
106         parameters.resize(0);
107       }
108     }
109     parameters.destroy();
110 
111     if (nAngles > 0) {
112       int forceGroup = anglePotentialEnergy.getForceGroup();
113       setForceGroup(forceGroup);
114       logger.info(format("  Angles:                            %10d", nAngles));
115       logger.fine(format("   Force Group:                      %10d", forceGroup));
116     }
117   }
118 
119   /**
120    * Create an OpenMM Angle Force.
121    *
122    * @param anglePotentialEnergy     The AnglePotentialEnergy instance.
123    * @param topology                 The topology index for the OpenMM System.
124    * @param openMMDualTopologyEnergy The OpenMMDualTopologyEnergy instance.
125    */
126   public AngleForce(AnglePotentialEnergy anglePotentialEnergy,
127                     int topology, OpenMMDualTopologyEnergy openMMDualTopologyEnergy) {
128     super(anglePotentialEnergy.getAngleEnergyString());
129     Angle[] angles = anglePotentialEnergy.getAngleArray();
130     addPerAngleParameter("theta0");
131     addPerAngleParameter("k");
132     setName("Angle");
133 
134     ForceFieldEnergy forceFieldEnergy = openMMDualTopologyEnergy.getForceFieldEnergy(topology);
135     ForceField forceField = forceFieldEnergy.getMolecularAssembly().getForceField();
136     manyBodyTitration = forceField.getBoolean("MANYBODY_TITRATION", false);
137     rigidHydrogenAngles = forceField.getBoolean("RIGID_HYDROGEN_ANGLES", false);
138     if (manyBodyTitration || rigidHydrogenAngles) {
139       logger.severe("Dual Topology does not support rigid hydrogen angles or many body titration.");
140     }
141 
142     double scale = openMMDualTopologyEnergy.getTopologyScale(topology);
143 
144     DoubleArray parameters = new DoubleArray(0);
145     for (Angle angle : angles) {
146       AngleType angleType = angle.angleType;
147       AngleType.AngleMode angleMode = angleType.angleMode;
148       if (angleMode == AngleType.AngleMode.IN_PLANE) {
149         // Skip In-Plane angles.
150       } else {
151         int i1 = angle.getAtom(0).getArrayIndex();
152         int i2 = angle.getAtom(1).getArrayIndex();
153         int i3 = angle.getAtom(2).getArrayIndex();
154 
155         double theta0 = angleType.angle[angle.nh];
156         double k = OpenMM_KJPerKcal * angleType.angleUnit * angleType.forceConstant;
157         // Don't apply lambda scale to alchemical angle
158         if (!angle.applyLambda()) {
159           k = k * scale;
160         }
161         parameters.append(theta0);
162         parameters.append(k);
163         i1 = openMMDualTopologyEnergy.mapToDualTopologyIndex(topology, i1);
164         i2 = openMMDualTopologyEnergy.mapToDualTopologyIndex(topology, i2);
165         i3 = openMMDualTopologyEnergy.mapToDualTopologyIndex(topology, i3);
166         addAngle(i1, i2, i3, parameters);
167         nAngles++;
168         parameters.resize(0);
169       }
170     }
171     parameters.destroy();
172 
173     if (nAngles > 0) {
174       int forceGroup = anglePotentialEnergy.getForceGroup();
175       setForceGroup(forceGroup);
176       logger.info(format("  Angles:                            %10d", nAngles));
177       logger.fine(format("   Force Group:                      %10d", forceGroup));
178     }
179   }
180 
181   /**
182    * Convenience method to construct an OpenMM Angle Force.
183    *
184    * @param openMMEnergy The OpenMM Energy instance that contains the angles.
185    * @return An Angle Force, or null if there are no angles.
186    */
187   public static Force constructForce(OpenMMEnergy openMMEnergy) {
188     AnglePotentialEnergy anglePotentialEnergy = openMMEnergy.getAnglePotentialEnergy();
189     if (anglePotentialEnergy == null) {
190       return null;
191     }
192     AngleForce angleForce = new AngleForce(anglePotentialEnergy, openMMEnergy);
193     if (angleForce.nAngles > 0) {
194       return angleForce;
195     }
196     angleForce.destroy();
197     return null;
198   }
199 
200   /**
201    * Add a bond force to the OpenMM System
202    *
203    * @param topology                 The topology index for the OpenMM System.
204    * @param openMMDualTopologyEnergy The OpenMMDualTopologyEnergy instance.
205    */
206   public static Force constructForce(int topology, OpenMMDualTopologyEnergy openMMDualTopologyEnergy) {
207     ForceFieldEnergy forceFieldEnergy = openMMDualTopologyEnergy.getForceFieldEnergy(topology);
208     AnglePotentialEnergy anglePotentialEnergy = forceFieldEnergy.getAnglePotentialEnergy();
209     if (anglePotentialEnergy == null) {
210       return null;
211     }
212     AngleForce angleForce = new AngleForce(anglePotentialEnergy, topology, openMMDualTopologyEnergy);
213     if (angleForce.nAngles > 0) {
214       return angleForce;
215     }
216     angleForce.destroy();
217     return null;
218   }
219 
220   /**
221    * Update an existing angle force for the OpenMM System.
222    *
223    * @param openMMEnergy The OpenMM Energy instance that contains the angles.
224    */
225   public void updateForce(OpenMMEnergy openMMEnergy) {
226     AnglePotentialEnergy anglePotentialEnergy = openMMEnergy.getAnglePotentialEnergy();
227     if (anglePotentialEnergy == null) {
228       return;
229     }
230     Angle[] angles = anglePotentialEnergy.getAngleArray();
231 
232     DoubleArray parameters = new DoubleArray(0);
233     int index = 0;
234     for (Angle angle : angles) {
235       AngleType.AngleMode angleMode = angle.angleType.angleMode;
236       if (!manyBodyTitration && angleMode == AngleType.AngleMode.IN_PLANE) {
237         // Skip In-Plane angles unless this is ManyBody Titration.
238       } else if (!rigidHydrogenAngles || !isHydrogenAngle(angle)) {
239         // Update angles that do not involve rigid hydrogen atoms.
240         int i1 = angle.getAtom(0).getArrayIndex();
241         int i2 = angle.getAtom(1).getArrayIndex();
242         int i3 = angle.getAtom(2).getArrayIndex();
243         double theta0 = angle.angleType.angle[angle.nh];
244         double k = OpenMM_KJPerKcal * angle.angleType.angleUnit * angle.angleType.forceConstant;
245         if (angleMode == AngleType.AngleMode.IN_PLANE) {
246           // Zero the force constant for In-Plane Angles.
247           k = 0.0;
248         }
249         parameters.append(theta0);
250         parameters.append(k);
251         setAngleParameters(index++, i1, i2, i3, parameters);
252         parameters.resize(0);
253       }
254     }
255     parameters.destroy();
256     updateParametersInContext(openMMEnergy.getContext());
257   }
258 
259   /**
260    * Update an existing angle force for the OpenMM System.
261    *
262    * @param topology                 The topology index for the OpenMM System.
263    * @param openMMDualTopologyEnergy The OpenMMDualTopologyEnergy instance.
264    */
265   public void updateForce(int topology, OpenMMDualTopologyEnergy openMMDualTopologyEnergy) {
266     ForceFieldEnergy forceFieldEnergy = openMMDualTopologyEnergy.getForceFieldEnergy(topology);
267     AnglePotentialEnergy anglePotentialEnergy = forceFieldEnergy.getAnglePotentialEnergy();
268     Angle[] angles = anglePotentialEnergy.getAngleArray();
269     if (angles == null || angles.length < 1) {
270       return;
271     }
272 
273     double scale = openMMDualTopologyEnergy.getTopologyScale(topology);
274 
275     DoubleArray parameters = new DoubleArray(0);
276     int index = 0;
277     for (Angle angle : angles) {
278       AngleType.AngleMode angleMode = angle.angleType.angleMode;
279       if (angleMode == AngleType.AngleMode.IN_PLANE) {
280         // Skip In-Plane angles.
281       }
282       // Update angles.
283       else {
284         int i1 = angle.getAtom(0).getArrayIndex();
285         int i2 = angle.getAtom(1).getArrayIndex();
286         int i3 = angle.getAtom(2).getArrayIndex();
287         double theta0 = angle.angleType.angle[angle.nh];
288         double k = OpenMM_KJPerKcal * angle.angleType.angleUnit * angle.angleType.forceConstant;
289         // Don't apply lambda scale to alchemical angle
290         if (!angle.applyLambda()) {
291           k = k * scale;
292         }
293         if (angleMode == AngleType.AngleMode.IN_PLANE) {
294           // Zero the force constant for In-Plane Angles.
295           k = 0.0;
296         }
297         parameters.append(theta0);
298         parameters.append(k);
299         i1 = openMMDualTopologyEnergy.mapToDualTopologyIndex(topology, i1);
300         i2 = openMMDualTopologyEnergy.mapToDualTopologyIndex(topology, i2);
301         i3 = openMMDualTopologyEnergy.mapToDualTopologyIndex(topology, i3);
302         setAngleParameters(index++, i1, i2, i3, parameters);
303         parameters.resize(0);
304       }
305     }
306     parameters.destroy();
307     updateParametersInContext(openMMDualTopologyEnergy.getContext());
308   }
309 
310   /**
311    * Check to see if an angle is a hydrogen angle. This method only returns true for hydrogen
312    * angles that are less than 160 degrees.
313    *
314    * @param angle Angle to check.
315    * @return boolean indicating whether an angle is a hydrogen angle that is less than 160 degrees.
316    */
317   private boolean isHydrogenAngle(Angle angle) {
318     if (angle.containsHydrogen()) {
319       // Equilibrium angle value in degrees.
320       double angleVal = angle.angleType.angle[angle.nh];
321       // Make sure angle is less than 160 degrees because greater than 160 degrees will not be
322       // constrained
323       // well using the law of cosines.
324       if (angleVal < 160.0) {
325         Atom atom1 = angle.getAtom(0);
326         Atom atom2 = angle.getAtom(1);
327         Atom atom3 = angle.getAtom(2);
328         // Setting constraints only on angles where atom1 or atom3 is a hydrogen while atom2 is
329         // not a hydrogen.
330         return atom1.isHydrogen() && atom3.isHydrogen() && !atom2.isHydrogen();
331       }
332     }
333     return false;
334   }
335 
336 }