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