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