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.parameters;
39  
40  import ffx.utilities.FFXProperty;
41  import org.w3c.dom.Document;
42  import org.w3c.dom.Element;
43  
44  import javax.annotation.Nullable;
45  import java.util.Arrays;
46  import java.util.Comparator;
47  import java.util.HashMap;
48  import java.util.Map;
49  import java.util.logging.Level;
50  import java.util.logging.Logger;
51  
52  import static ffx.potential.parameters.ForceField.ForceFieldType.ANGLE;
53  import static ffx.potential.parameters.ForceField.ForceFieldType.ANGLEP;
54  import static ffx.utilities.Constants.DEGREES_PER_RADIAN;
55  import static ffx.utilities.Constants.KCAL_TO_KJ;
56  import static ffx.utilities.PropertyGroup.EnergyUnitConversion;
57  import static ffx.utilities.PropertyGroup.LocalGeometryFunctionalForm;
58  import static ffx.utilities.PropertyGroup.PotentialFunctionParameter;
59  import static java.lang.Double.parseDouble;
60  import static java.lang.Integer.parseInt;
61  import static java.lang.String.format;
62  import static java.lang.String.valueOf;
63  import static java.lang.System.arraycopy;
64  import static org.apache.commons.math3.util.FastMath.PI;
65  import static org.apache.commons.math3.util.FastMath.pow;
66  
67  /**
68   * The AngleType class defines one harmonic angle bend energy term.
69   *
70   * @author Michael J. Schnieders
71   * @since 1.0
72   */
73  @FFXProperty(name = "angle", clazz = String.class, propertyGroup = PotentialFunctionParameter, description = """
74      [3 integers and 4 reals]
75      Provides the values for a single bond angle bending parameter.
76      The integer modifiers give the atom class numbers for the three kinds of atoms involved in the angle which is to be defined.
77      The real number modifiers give the force constant value for the angle and up to three ideal bond angles in degrees.
78      In most cases only one ideal bond angle is given, and that value is used for all occurrences of the specified bond angle.
79      If all three ideal angles are given, the values apply when the central atom of the angle is attached to 0, 1 or 2 additional hydrogen atoms, respectively.
80      This "hydrogen environment" option is provided to implement the corresponding feature of the AMOEBA force field.
81      The default units for the force constant are kcal/mole/radian^2, but this can be controlled via the angleunit keyword.
82      """)
83  @FFXProperty(name = "anglep", clazz = String.class, propertyGroup = PotentialFunctionParameter, description = """
84      [3 integers and 3 reals]
85      Provides the values for a single projected in-plane bond angle bending parameter.
86      The integer modifiers give the atom class numbers for the three kinds of atoms involved in the angle which is to be defined.
87      The real number modifiers give the force constant value for the angle and up to two ideal bond angles in degrees.
88      In most cases only one ideal bond angle is given, and that value is used for all occurrences of the specified bond angle.
89      If all two ideal angles are given, the values apply when the central atom of the angle is attached to 0 or 1 additional hydrogen atoms, respectively.
90      This "hydrogen environment" option is provided to implement the corresponding feature of the AMOEBA force field.
91      The default units for the force constant are kcal/mole/radian^2, but this can be controlled via the angleunit keyword.
92      """)
93  public final class AngleType extends BaseType implements Comparator<String> {
94  
95    /**
96     * Default convert angle bending energy to kcal/mole.
97     */
98    public static final double DEFAULT_ANGLE_UNIT = pow(PI / 180.0, 2.0);
99    /**
100    * Default cubic coefficient in angle bending potential.
101    */
102   public static final double DEFAULT_ANGLE_CUBIC = 0.0;
103   /**
104    * Default quartic coefficient in angle bending potential.
105    */
106   public static final double DEFAULT_ANGLE_QUARTIC = 0.0;
107   /**
108    * Default pentic coefficient in angle bending potential.
109    */
110   public static final double DEFAULT_ANGLE_PENTIC = 0.0;
111   /**
112    * Default quintic coefficient in angle bending potential.
113    */
114   public static final double DEFAULT_ANGLE_SEXTIC = 0.0;
115 
116   /**
117    * Convert angle bending energy to kcal/mole.
118    */
119   @FFXProperty(name = "angleunit", propertyGroup = EnergyUnitConversion, defaultValue = "(Pi/180)^2", description = """
120       Sets the scale factor needed to convert the energy value computed by the bond angle bending potential into units of kcal/mole.
121       The correct value is force field dependent and typically provided in the header of the master force field parameter file.
122       """)
123   public double angleUnit = DEFAULT_ANGLE_UNIT;
124 
125   /**
126    * Cubic coefficient in angle bending potential.
127    */
128   @FFXProperty(name = "angle-cubic", propertyGroup = LocalGeometryFunctionalForm, defaultValue = "0.0", description = """
129       Sets the value of the cubic term in the Taylor series expansion form of the bond angle bending potential energy.
130       The real number modifier gives the value of the coefficient as a multiple of the quadratic coefficient.
131       This term multiplied by the angle bending energy unit conversion factor, the force constant,
132       and the cube of the deviation of the bond angle from its ideal value gives the cubic contribution to the angle bending energy.
133       The default value in the absence of the angle-cubic keyword is zero; i.e., the cubic angle bending term is omitted.
134       """)
135   public double cubic = DEFAULT_ANGLE_CUBIC;
136 
137   /**
138    * Quartic coefficient in angle bending potential.
139    */
140   @FFXProperty(name = "angle-quartic", propertyGroup = LocalGeometryFunctionalForm, defaultValue = "0.0", description = """
141       Sets the value of the quartic term in the Taylor series expansion form of the bond angle bending potential energy.
142       The real number modifier gives the value of the coefficient as a multiple of the quadratic coefficient.
143       This term multiplied by the angle bending energy unit conversion factor, the force constant,
144       and the forth power of the deviation of the bond angle from its ideal value gives the quartic contribution to the angle bending energy.
145       The default value in the absence of the angle-quartic keyword is zero; i.e., the quartic angle bending term is omitted.
146       """)
147   public double quartic = DEFAULT_ANGLE_QUARTIC;
148 
149   /**
150    * Pentic coefficient in angle bending potential.
151    */
152   @FFXProperty(name = "angle-pentic", propertyGroup = LocalGeometryFunctionalForm, defaultValue = "0.0", description = """
153       Sets the value of the fifth power term in the Taylor series expansion form of the bond angle bending potential energy.
154       The real number modifier gives the value of the coefficient as a multiple of the quadratic coefficient.
155       This term multiplied by the angle bending energy unit conversion factor, the force constant,
156       and the fifth power of the deviation of the bond angle from its ideal value gives the pentic contribution to the angle bending energy.
157       The default value in the absence of the angle-pentic keyword is zero; i.e., the pentic angle bending term is omitted.
158       """)
159   public double pentic = DEFAULT_ANGLE_PENTIC;
160 
161   /**
162    * Sextic coefficient in angle bending potential.
163    */
164   @FFXProperty(name = "angle-sextic", propertyGroup = LocalGeometryFunctionalForm, defaultValue = "0.0", description = """
165       Sets the value of the sixth power term in the Taylor series expansion form of the bond angle bending potential energy.
166       The real number modifier gives the value of the coefficient as a multiple of the quadratic coefficient.
167       This term multiplied by the angle bending energy unit conversion factor, the force constant,
168       and the sixth power of the deviation of the bond angle from its ideal value gives the sextic contribution to the angle bending energy.
169       The default value in the absence of the angle-sextic keyword is zero; i.e., the sextic angle bending term is omitted.
170       """)
171   public double sextic = DEFAULT_ANGLE_SEXTIC;
172 
173   /**
174    * A Logger for the AngleType class.
175    */
176   private static final Logger logger = Logger.getLogger(AngleType.class.getName());
177   /**
178    * Atom classes that for this Angle type.
179    */
180   public final int[] atomClasses;
181   /**
182    * Force constant (Kcal/mole/radian^2).
183    */
184   public final double forceConstant;
185   /**
186    * Equilibrium angle (degrees). There can be up to three equilibrium angles, depending on the
187    * number of attached hydrogens (0, 1, or 2).
188    */
189   public final double[] angle;
190   /**
191    * The angle mode in use.
192    */
193   public final AngleMode angleMode;
194   /**
195    * The angle function in use.
196    */
197   public AngleFunction angleFunction;
198 
199   /**
200    * The default AngleType constructor defines use of the Sextic AngleFunction.
201    *
202    * @param atomClasses   an array of int.
203    * @param forceConstant a double.
204    * @param angle         an array of double.
205    */
206   public AngleType(int[] atomClasses, double forceConstant, double[] angle) {
207     this(atomClasses, forceConstant, angle, AngleMode.NORMAL);
208   }
209 
210   /**
211    * Constructor for In-Plane AngleType.
212    *
213    * @param atomClasses   an array of int.
214    * @param forceConstant a double.
215    * @param angle         an array of double.
216    * @param angleMode     the AngleMode to apply.
217    */
218   public AngleType(int[] atomClasses, double forceConstant, double[] angle, AngleMode angleMode) {
219     this(atomClasses, forceConstant, angle, angleMode, AngleFunction.SEXTIC);
220   }
221 
222   /**
223    * Constructor for In-Plane AngleType.
224    *
225    * @param atomClasses   an array of int.
226    * @param forceConstant a double.
227    * @param angle         an array of double.
228    * @param angleMode     the AngleMode to apply.
229    * @param angleFunction the AngleFunction to use.
230    */
231   public AngleType(int[] atomClasses, double forceConstant, double[] angle, AngleMode angleMode,
232                    AngleFunction angleFunction) {
233     super(ANGLE, sortKey(atomClasses));
234     this.atomClasses = atomClasses;
235     this.forceConstant = forceConstant;
236     this.angle = angle;
237     this.angleMode = angleMode;
238     this.angleFunction = angleFunction;
239     if (angleMode == AngleMode.IN_PLANE) {
240       forceFieldType = ANGLEP;
241     }
242   }
243 
244   /**
245    * Average two AngleType instances. The atom classes that define the new type must be supplied.
246    * If the angle modes or angle functions differ, null is returned.
247    *
248    * @param angleType1  a {@link ffx.potential.parameters.AngleType} object.
249    * @param angleType2  a {@link ffx.potential.parameters.AngleType} object.
250    * @param atomClasses the atom classes for the new type.
251    * @return a {@link ffx.potential.parameters.AngleType} object or null if averaging is not possible.
252    */
253   public static AngleType average(@Nullable AngleType angleType1,
254                                   @Nullable AngleType angleType2,
255                                   @Nullable int[] atomClasses) {
256     if (angleType1 == null || angleType2 == null || atomClasses == null) {
257       return null;
258     }
259     AngleMode angleMode = angleType1.angleMode;
260     if (angleMode != angleType2.angleMode) {
261       return null;
262     }
263     AngleFunction angleFunction = angleType1.angleFunction;
264     if (angleFunction != angleType2.angleFunction) {
265       return null;
266     }
267     int len = angleType1.angle.length;
268     if (len != angleType2.angle.length) {
269       return null;
270     }
271     double forceConstant = (angleType1.forceConstant + angleType2.forceConstant) / 2.0;
272     double[] angle = new double[len];
273     for (int i = 0; i < len; i++) {
274       angle[i] = (angleType1.angle[i] + angleType2.angle[i]) / 2.0;
275     }
276 
277     return new AngleType(atomClasses, forceConstant, angle, angleMode, angleFunction);
278   }
279 
280   /**
281    * Construct an AngleType from an input string.
282    *
283    * @param input  The overall input String.
284    * @param tokens The input String tokenized.
285    * @return an AngleType instance.
286    */
287   public static AngleType parse(String input, String[] tokens) {
288     if (tokens.length < 6) {
289       logger.log(Level.WARNING, "Invalid ANGLE type:\n{0}", input);
290     } else {
291       int[] atomClasses = new int[3];
292       double forceConstant;
293       int angles;
294       double[] bondAngle;
295       try {
296         atomClasses[0] = parseInt(tokens[1]);
297         atomClasses[1] = parseInt(tokens[2]);
298         atomClasses[2] = parseInt(tokens[3]);
299         forceConstant = parseDouble(tokens[4]);
300         angles = tokens.length - 5;
301         bondAngle = new double[angles];
302         for (int i = 0; i < angles; i++) {
303           bondAngle[i] = parseDouble(tokens[5 + i]);
304         }
305       } catch (NumberFormatException e) {
306         String message = "Exception parsing ANGLE type:\n" + input + "\n";
307         logger.log(Level.SEVERE, message, e);
308         return null;
309       }
310       double[] newBondAngle = new double[angles];
311       arraycopy(bondAngle, 0, newBondAngle, 0, angles);
312       return new AngleType(atomClasses, forceConstant, newBondAngle);
313     }
314     return null;
315   }
316 
317   /**
318    * Construct an In-Plane AngleType from an input string.
319    *
320    * @param input  The overall input String.
321    * @param tokens The input String tokenized.
322    * @return an AngleType instance.
323    */
324   public static AngleType parseInPlane(String input, String[] tokens) {
325     if (tokens.length < 6) {
326       logger.log(Level.WARNING, "Invalid ANGLEP type:\n{0}", input);
327     } else {
328       int[] atomClasses = new int[3];
329       double forceConstant;
330       int angles;
331       double[] bondAngle;
332       try {
333         atomClasses[0] = parseInt(tokens[1]);
334         atomClasses[1] = parseInt(tokens[2]);
335         atomClasses[2] = parseInt(tokens[3]);
336         forceConstant = parseDouble(tokens[4]);
337         angles = tokens.length - 5;
338         bondAngle = new double[angles];
339         for (int i = 0; i < angles; i++) {
340           bondAngle[i] = parseDouble(tokens[5 + i]);
341         }
342       } catch (NumberFormatException e) {
343         String message = "Exception parsing ANGLEP type:\n" + input + "\n";
344         logger.log(Level.SEVERE, message, e);
345         return null;
346       }
347       double[] newBondAngle = new double[angles];
348       arraycopy(bondAngle, 0, newBondAngle, 0, angles);
349       return new AngleType(atomClasses, forceConstant, newBondAngle, AngleMode.IN_PLANE);
350     }
351     return null;
352   }
353 
354   /**
355    * This method sorts the atom classes as: min, c[1], max
356    *
357    * @param c atomClasses
358    * @return lookup key
359    */
360   public static String sortKey(int[] c) {
361     if (c == null || c.length != 3) {
362       return null;
363     }
364     if (c[0] > c[2]) {
365       int temp = c[0];
366       c[0] = c[2];
367       c[2] = temp;
368     }
369     return c[0] + " " + c[1] + " " + c[2];
370   }
371 
372   /**
373    * {@inheritDoc}
374    */
375   @Override
376   public int compare(String key1, String key2) {
377     String[] keys1 = key1.split(" ");
378     String[] keys2 = key2.split(" ");
379     int[] c1 = new int[3];
380     int[] c2 = new int[3];
381     for (int i = 0; i < 3; i++) {
382       c1[i] = Integer.parseInt(keys1[i]);
383       c2[i] = Integer.parseInt(keys2[i]);
384     }
385 
386     if (c1[1] < c2[1]) {
387       return -1;
388     } else if (c1[1] > c2[1]) {
389       return 1;
390     } else if (c1[0] < c2[0]) {
391       return -1;
392     } else if (c1[0] > c2[0]) {
393       return 1;
394     } else if (c1[2] < c2[2]) {
395       return -1;
396     } else if (c1[2] > c2[2]) {
397       return 1;
398     }
399 
400     return 0;
401   }
402 
403   /**
404    * {@inheritDoc}
405    */
406   @Override
407   public boolean equals(Object o) {
408     if (this == o) {
409       return true;
410     }
411     if (o == null || getClass() != o.getClass()) {
412       return false;
413     }
414     AngleType angleType = (AngleType) o;
415     return Arrays.equals(atomClasses, angleType.atomClasses);
416   }
417 
418   /**
419    * {@inheritDoc}
420    */
421   @Override
422   public int hashCode() {
423     return Arrays.hashCode(atomClasses);
424   }
425 
426   /**
427    * incrementClasses
428    *
429    * @param increment the value to add to each atom class.
430    */
431   public void incrementClasses(int increment) {
432     for (int i = 0; i < atomClasses.length; i++) {
433       atomClasses[i] += increment;
434     }
435     setKey(sortKey(atomClasses));
436   }
437 
438   /**
439    * Remap new atom classes to known internal ones.
440    *
441    * @param typeMap a lookup between new atom types and known atom types.
442    * @return a {@link ffx.potential.parameters.AngleType} object.
443    */
444   public AngleType patchClasses(HashMap<AtomType, AtomType> typeMap) {
445     int count = 0;
446     int len = atomClasses.length;
447 
448     // Look for new AngleTypes that contain 1 or 2 mapped atom classes.
449     for (AtomType newType : typeMap.keySet()) {
450       for (int atomClass : atomClasses) {
451         if (atomClass == newType.atomClass) {
452           count++;
453         }
454       }
455     }
456 
457     // If found, create a new AngleType that bridges to known classes.
458     if (count == 1 || count == 2) {
459       int[] newClasses = Arrays.copyOf(atomClasses, len);
460       for (AtomType newType : typeMap.keySet()) {
461         for (int i = 0; i < len; i++) {
462           if (atomClasses[i] == newType.atomClass) {
463             AtomType knownType = typeMap.get(newType);
464             newClasses[i] = knownType.atomClass;
465           }
466         }
467       }
468       return new AngleType(newClasses, forceConstant, angle, angleMode, angleFunction);
469     }
470 
471     return null;
472   }
473 
474   /**
475    * Set the AngleFunction.
476    *
477    * @param angleFunction the AngleFunction.
478    */
479   public void setAngleFunction(AngleFunction angleFunction) {
480     this.angleFunction = angleFunction;
481   }
482 
483   /**
484    * {@inheritDoc}
485    *
486    * <p>Nicely formatted Angle bending string.
487    */
488   @Override
489   public String toString() {
490     StringBuilder angleString;
491     if (angleMode == AngleMode.NORMAL) {
492       angleString = new StringBuilder(
493           format("angle  %5d  %5d  %5d  %6.2f", atomClasses[0], atomClasses[1], atomClasses[2],
494               forceConstant));
495     } else {
496       angleString = new StringBuilder(
497           format("anglep  %5d  %5d  %5d  %6.2f", atomClasses[0], atomClasses[1], atomClasses[2],
498               forceConstant));
499     }
500     for (double eq : angle) {
501       angleString.append(format("  %6.2f", eq));
502     }
503     return angleString.toString();
504   }
505 
506   /**
507    * Create an AmoebaAngleForce
508    *
509    * @param doc        the Document instance.
510    * @param forceField the ForceField to define constants from.
511    * @return the AmoebaAngleForce Element.
512    */
513   public static Element getXMLForce(Document doc, ForceField forceField) {
514     Map<String, AngleType> angMap = forceField.getAngleTypes();
515     angMap.putAll(forceField.getAnglepTypes());
516     if (!angMap.values().isEmpty()) {
517       Element node = doc.createElement("AmoebaAngleForce");
518       node.setAttribute("angle-cubic", valueOf(forceField.getDouble("angle-cubic", DEFAULT_ANGLE_CUBIC)));
519       node.setAttribute("angle-quartic", valueOf(forceField.getDouble("angle-quartic", DEFAULT_ANGLE_QUARTIC)));
520       node.setAttribute("angle-pentic", valueOf(forceField.getDouble("angle-pentic", DEFAULT_ANGLE_PENTIC)));
521       node.setAttribute("angle-sextic", valueOf(forceField.getDouble("angle-sextic", DEFAULT_ANGLE_SEXTIC)));
522       for (AngleType angleType : angMap.values()) {
523         node.appendChild(angleType.toXML(doc));
524       }
525       return node;
526     }
527     return null;
528   }
529 
530   /**
531    * Write AngleType to OpenMM XML format.
532    *
533    * @param doc the Document instance.
534    * @return the Angle element.
535    */
536   public Element toXML(Document doc) {
537     Element node = doc.createElement("Angle");
538     node.setAttribute("class1", format("%d", atomClasses[0]));
539     node.setAttribute("class2", format("%d", atomClasses[1]));
540     node.setAttribute("class3", format("%d", atomClasses[2]));
541     // Convert Kcal/mol/radian^2 to KJ/mol/deg^2
542     node.setAttribute("k", format("%.17f", forceConstant * KCAL_TO_KJ / (DEGREES_PER_RADIAN * DEGREES_PER_RADIAN)));
543     int i = 1;
544     for (double eq : angle) {
545       node.setAttribute(format("angle%d", i), format("%f", eq));
546       i++;
547     }
548     if (angleMode == AngleMode.NORMAL) {
549       node.setAttribute("inPlane", "False");
550     } else {
551       node.setAttribute("inPlane", "True");
552     }
553     return node;
554   }
555 
556   /**
557    * Angle function types include harmonic or sextic.
558    */
559   public enum AngleFunction {
560     HARMONIC, SEXTIC
561   }
562 
563   /**
564    * Angle modes include Normal or In-Plane
565    */
566   public enum AngleMode {
567     NORMAL, IN_PLANE
568   }
569 }