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-2023.
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.BOND;
41  import static ffx.utilities.KeywordGroup.EnergyUnitConversion;
42  import static ffx.utilities.KeywordGroup.LocalGeometryFunctionalForm;
43  import static ffx.utilities.KeywordGroup.PotentialFunctionParameter;
44  import static java.lang.Double.parseDouble;
45  import static java.lang.Integer.parseInt;
46  import static java.lang.String.format;
47  
48  import ffx.utilities.FFXKeyword;
49  
50  import java.util.Arrays;
51  import java.util.Comparator;
52  import java.util.HashMap;
53  import java.util.logging.Level;
54  import java.util.logging.Logger;
55  
56  /**
57   * The BondType class defines one harmonic bond stretch energy term.
58   *
59   * @author Michael J. Schnieders
60   * @since 1.0
61   */
62  @FFXKeyword(name = "bond", clazz = String.class, keywordGroup = PotentialFunctionParameter, description =
63      "[2 integers and 2 reals] " + "Provides the values for a single bond stretching parameter. "
64          + "The integer modifiers give the atom class numbers for the two kinds of atoms involved in the bond which is to be defined. "
65          + "The real number modifiers give the force constant value for the bond and the ideal bond length in Angstroms. "
66          + "The default value of 1.0 is used, if the bondunit keyword is not given in the force field parameter file or the keyfile.")
67  public final class BondType extends BaseType implements Comparator<String> {
68  
69    public static final double DEFAULT_BOND_UNIT = 1.0;
70  
71    /**
72     * Default cubic coefficient in bond stretch potential.
73     */
74    public static final double DEFAULT_BOND_CUBIC = 0.0;
75    /**
76     * Default quartic coefficient in bond stretch potential.
77     */
78    public static final double DEFAULT_BOND_QUARTIC = 0.0;
79  
80    /**
81     * Convert bond stretch energy to kcal/mole.
82     */
83    @FFXKeyword(name = "bondunit", keywordGroup = EnergyUnitConversion, defaultValue = "1.0", description =
84        "Sets the scale factor needed to convert the energy value computed by the bond stretching potential into units of kcal/mole. "
85            + "The correct value is force field dependent and typically provided in the header of the master force field parameter file.")
86    public double bondUnit = DEFAULT_BOND_UNIT;
87  
88    /**
89     * Cubic coefficient in bond stretch potential.
90     */
91    @FFXKeyword(name = "bond-cubic", keywordGroup = LocalGeometryFunctionalForm, defaultValue = "0.0", description =
92        "Sets the value of the cubic term in the Taylor series expansion form of the bond stretching potential energy. "
93            + "The real number modifier gives the value of the coefficient as a multiple of the quadratic coefficient. "
94            + "This term multiplied by the bond stretching energy unit conversion factor, the force constant, "
95            + "and the cube of the deviation of the bond length from its ideal value gives the cubic contribution to the bond stretching energy. "
96            + "The default value in the absence of the bond-cubic keyword is zero; i.e., the cubic bond stretching term is omitted.")
97    public double cubic = DEFAULT_BOND_CUBIC;
98  
99    /**
100    * Quartic coefficient in bond stretch potential.
101    */
102   @FFXKeyword(name = "bond-quartic", keywordGroup = LocalGeometryFunctionalForm, defaultValue = "0.0", description =
103       "Sets the value of the quartic term in the Taylor series expansion form of the bond stretching potential energy. "
104           + "The real number modifier gives the value of the coefficient as a multiple of the quadratic coefficient. "
105           + "This term multiplied by the bond stretching energy unit conversion factor, the force constant, "
106           + "and the forth power of the deviation of the bond length from its ideal value gives the quartic contribution to the bond stretching energy. "
107           + "The default value in the absence of the bond-quartic keyword is zero; i.e., the quartic bond stretching term is omitted.")
108   public double quartic = DEFAULT_BOND_QUARTIC;
109 
110   /**
111    * A Logger for the BondType class.
112    */
113   private static final Logger logger = Logger.getLogger(BondType.class.getName());
114   /**
115    * Atom classes that form this bond stretch.
116    */
117   public final int[] atomClasses;
118   /**
119    * Force constant (Kcal/mol).
120    */
121   public final double forceConstant;
122   /**
123    * Equilibrium separation (Angstroms).
124    */
125   public final double distance;
126   /**
127    * Radius of a flat bottom where energy and force is 0; typically used for restraints. Will almost
128    * always be 0.
129    */
130   public final double flatBottomRadius;
131   /**
132    * The function used by the bond: harmonic or quartic with flat-bottom variants.
133    */
134   public BondFunction bondFunction;
135 
136   /**
137    * The default BondType constructor defines use of the Quartic BondFunction.
138    *
139    * @param atomClasses   int[]
140    * @param forceConstant double
141    * @param distance      double
142    */
143   public BondType(int[] atomClasses, double forceConstant, double distance) {
144     this(atomClasses, forceConstant, distance, BondFunction.QUARTIC, 0.0);
145   }
146 
147   /**
148    * BondType constructor.
149    *
150    * @param atomClasses   int[]
151    * @param forceConstant double
152    * @param distance      double
153    * @param bondFunction  the BondFunction type to apply.
154    */
155   public BondType(int[] atomClasses, double forceConstant, double distance,
156                   BondFunction bondFunction) {
157     this(atomClasses, forceConstant, distance, bondFunction, 0.0);
158   }
159 
160   /**
161    * BondType constructor.
162    *
163    * @param atomClasses      int[]
164    * @param forceConstant    double
165    * @param distance         double
166    * @param bondFunction     the BondFunction type to apply.
167    * @param flatBottomRadius a double.
168    */
169   public BondType(int[] atomClasses, double forceConstant, double distance, BondFunction bondFunction,
170                   double flatBottomRadius) {
171     super(BOND, sortKey(atomClasses));
172     this.atomClasses = atomClasses;
173     this.forceConstant = forceConstant;
174     this.distance = distance;
175     this.bondFunction = bondFunction;
176     this.flatBottomRadius = flatBottomRadius;
177     assert (flatBottomRadius == 0 || bondFunction == BondFunction.FLAT_BOTTOM_HARMONIC);
178   }
179 
180   /**
181    * Average two BondType instances. The atom classes that define the new type must be supplied.
182    *
183    * @param bondType1   a {@link ffx.potential.parameters.BondType} object.
184    * @param bondType2   a {@link ffx.potential.parameters.BondType} object.
185    * @param atomClasses an array of {@link int} objects.
186    * @return a {@link ffx.potential.parameters.BondType} object.
187    */
188   public static BondType average(BondType bondType1, BondType bondType2, int[] atomClasses) {
189     if (bondType1 == null || bondType2 == null || atomClasses == null) {
190       return null;
191     }
192     BondType.BondFunction bondFunction = bondType1.bondFunction;
193     if (bondFunction != bondType2.bondFunction) {
194       return null;
195     }
196 
197     double forceConstant = (bondType1.forceConstant + bondType2.forceConstant) / 2.0;
198     double distance = (bondType1.distance + bondType2.distance) / 2.0;
199 
200     return new BondType(atomClasses, forceConstant, distance, bondFunction);
201   }
202 
203   /**
204    * Construct a BondType from an input string.
205    *
206    * @param input  The overall input String.
207    * @param tokens The input String tokenized.
208    * @return a BondType instance.
209    */
210   public static BondType parse(String input, String[] tokens) {
211     if (tokens.length < 5) {
212       logger.log(Level.WARNING, "Invalid BOND type:\n{0}", input);
213     } else {
214       try {
215         int[] atomClasses = new int[2];
216         atomClasses[0] = parseInt(tokens[1]);
217         atomClasses[1] = parseInt(tokens[2]);
218         double forceConstant = parseDouble(tokens[3]);
219         double distance = parseDouble(tokens[4]);
220         return new BondType(atomClasses, forceConstant, distance);
221       } catch (NumberFormatException e) {
222         String message = "Exception parsing BOND type:\n" + input + "\n";
223         logger.log(Level.SEVERE, message, e);
224       }
225     }
226     return null;
227   }
228 
229   /**
230    * This method sorts the atom classes as: min, max
231    *
232    * @param c atomClasses
233    * @return lookup key
234    */
235   public static String sortKey(int[] c) {
236     if (c == null || c.length != 2) {
237       return null;
238     }
239 
240     int temp;
241     if (c[1] <= c[0]) {
242       temp = c[1];
243       c[1] = c[0];
244       c[0] = temp;
245     }
246 
247     return c[0] + " " + c[1];
248   }
249 
250   /**
251    * {@inheritDoc}
252    */
253   @Override
254   public int compare(String key1, String key2) {
255     String[] keys1 = key1.split(" ");
256     String[] keys2 = key2.split(" ");
257     int[] c1 = new int[2];
258     int[] c2 = new int[2];
259     for (int i = 0; i < 2; i++) {
260       c1[i] = parseInt(keys1[i]);
261       c2[i] = parseInt(keys2[i]);
262     }
263 
264     if (c1[0] < c2[0]) {
265       return -1;
266     } else if (c1[0] > c2[0]) {
267       return 1;
268     } else if (c1[1] < c2[1]) {
269       return -1;
270     } else if (c1[1] > c2[1]) {
271       return 1;
272     }
273 
274     return 0;
275   }
276 
277   /**
278    * {@inheritDoc}
279    */
280   @Override
281   public boolean equals(Object o) {
282     if (this == o) {
283       return true;
284     }
285     if (o == null || getClass() != o.getClass()) {
286       return false;
287     }
288     BondType bondType = (BondType) o;
289     return Arrays.equals(atomClasses, bondType.atomClasses);
290   }
291 
292   /**
293    * {@inheritDoc}
294    */
295   @Override
296   public int hashCode() {
297     return Arrays.hashCode(atomClasses);
298   }
299 
300   /**
301    * incrementClasses
302    *
303    * @param increment The increment to apply to the atom classes.
304    */
305   public void incrementClasses(int increment) {
306     for (int i = 0; i < atomClasses.length; i++) {
307       atomClasses[i] += increment;
308     }
309     setKey(sortKey(atomClasses));
310   }
311 
312   /**
313    * Remap new atom classes to known internal ones.
314    *
315    * @param typeMap a lookup between new atom types and known atom types.
316    * @return a {@link ffx.potential.parameters.BondType} object.
317    */
318   public BondType patchClasses(HashMap<AtomType, AtomType> typeMap) {
319 
320     int count = 0;
321     int len = atomClasses.length;
322 
323     // Look for new BondTypes that contain one mapped atom class.
324     for (AtomType newType : typeMap.keySet()) {
325 
326       for (int atomClass : atomClasses) {
327         if (atomClass == newType.atomClass) {
328           count++;
329         }
330       }
331     }
332 
333     // If found, create a new BondType that bridges to known classes.
334     if (count == 1) {
335       int[] newClasses = Arrays.copyOf(atomClasses, len);
336       for (AtomType newType : typeMap.keySet()) {
337         for (int i = 0; i < len; i++) {
338           if (atomClasses[i] == newType.atomClass) {
339             AtomType knownType = typeMap.get(newType);
340             newClasses[i] = knownType.atomClass;
341           }
342         }
343       }
344 
345       return new BondType(newClasses, forceConstant, distance, bondFunction);
346     }
347     return null;
348   }
349 
350   public void setBondFunction(BondFunction bondFunction) {
351     this.bondFunction = bondFunction;
352   }
353 
354   /**
355    * {@inheritDoc}
356    *
357    * <p>Nicely formatted bond stretch string.
358    */
359   @Override
360   public String toString() {
361     return format("bond  %5d  %5d  %7.2f  %7.4f", atomClasses[0], atomClasses[1], forceConstant,
362         distance);
363   }
364 
365   /**
366    * Describes the function used by the bond. Currently supported: harmonic and quartic with
367    * flat-bottom variants.
368    */
369   public enum BondFunction {
370 
371     // Harmonic bond function.
372     HARMONIC("0.5*k*(r-r0)^2", false),
373 
374     // Quartic bond function.
375     QUARTIC("0.5*k*dv^2*((1+cubic)*dv+(1+quartic)*dv^2);dv=r-r0", false),
376 
377     // Flat bottom harmonic bond function for restraints.
378     FLAT_BOTTOM_HARMONIC(
379         "0.5*k*dv^2;dv=step(dv)*step(dv-fb)*(dv-fb)" + "+step(-dv)*step(-dv-fb)*(-dv-fb);dv=r-r0",
380         true),
381 
382     // Flat bottom Quartic bond function for restraints.
383     FLAT_BOTTOM_QUARTIC("0.5*k*dv^2*((1+cubic)*dv+(1+quartic)*dv^2);"
384         + "dv=step(dv)*step(dv-fb)*(dv-fb)+step(-dv)*step(-dv-fb)*(-dv-fb);dv=r-r0", true);
385 
386     /**
387      * String representation of the mathematical form.
388      */
389     private final String mathematicalForm;
390 
391     /**
392      * Flag to indicate if the bond function has a flat bottom.
393      */
394     private final boolean hasFlatBottom;
395 
396     /**
397      * BondFunction constructor.
398      *
399      * @param mathForm   String representation of the mathematical form.
400      * @param flatBottom Flag to indicate if the bond function has a flat bottom.
401      */
402     BondFunction(String mathForm, boolean flatBottom) {
403       this.mathematicalForm = mathForm;
404       this.hasFlatBottom = flatBottom;
405     }
406 
407     /**
408      * Returns whether this BondFunction has a flat bottom.
409      *
410      * @return Flat bottom.
411      */
412     public boolean hasFlatBottom() {
413       return hasFlatBottom;
414     }
415 
416     /**
417      * Returns the form of this bond as a mathematical expression parsable by OpenMM.
418      *
419      * @return A string describing mathematical form.
420      */
421     public String toMathematicalForm() {
422       return mathematicalForm;
423     }
424   }
425 }