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