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.UREYBRAD;
53  import static ffx.utilities.Constants.ANG_TO_NM;
54  import static ffx.utilities.Constants.KCAL_TO_KJ;
55  import static ffx.utilities.PropertyGroup.EnergyUnitConversion;
56  import static ffx.utilities.PropertyGroup.LocalGeometryFunctionalForm;
57  import static ffx.utilities.PropertyGroup.PotentialFunctionParameter;
58  import static java.lang.Double.parseDouble;
59  import static java.lang.Integer.parseInt;
60  import static java.lang.String.format;
61  import static java.lang.String.valueOf;
62  
63  /**
64   * The UreyBradleyType class defines one harmonic UreyBradley cross term.
65   *
66   * @author Michael J. Schnieders
67   * @since 1.0
68   */
69  @FFXProperty(name = "ureybrad", clazz = String.class, propertyGroup = PotentialFunctionParameter, description = """
70      [3 integers and 2 reals]
71      Provides the values for a single Urey-Bradley cross term potential parameter.
72      The integer modifiers give the atom class numbers for the three kinds of atoms
73      involved in the angle for which a Urey-Bradley term is to be defined.
74      The real number modifiers give the force constant value for the term and the target value for the 1-3 distance in Angstroms.
75      The default units for the force constant are kcal/mole/Ang^2, but this can be controlled via the ureyunit keyword
76      """)
77  public final class UreyBradleyType extends BaseType implements Comparator<String> {
78  
79    /**
80     * Default conversion Urey-Bradley stretch energy to kcal/mole.
81     */
82    public static final double DEFAULT_UREY_UNIT = 1.0;
83    /**
84     * Default cubic coefficient in Urey-Bradley stretch potential.
85     */
86    public static final double DEFAULT_UREY_CUBIC = 0.0;
87    /**
88     * Default quartic coefficient in Urey-Bradley stretch potential.
89     */
90    public static final double DEFAULT_UREY_QUARTIC = 0.0;
91  
92    /**
93     * Convert Urey-Bradley stretch energy to kcal/mole.
94     */
95    @FFXProperty(name = "ureyunit", propertyGroup = EnergyUnitConversion, defaultValue = "1.0", description = """
96        Sets the scale factor needed to convert the energy value computed by the Urey-Bradley potential into units of kcal/mole.
97        The correct value is force field dependent and typically provided in the header of the master force field parameter file.
98        """)
99    public double ureyUnit = DEFAULT_UREY_UNIT;
100 
101   /**
102    * Cubic coefficient in bond stretch potential.
103    */
104   @FFXProperty(name = "urey-cubic", propertyGroup = LocalGeometryFunctionalForm, defaultValue = "0.0", description = """
105       Sets the value of the cubic term in the Taylor series expansion form of the Urey-Bradley potential energy.
106       The real number modifier gives the value of the coefficient as a multiple of the quadratic coefficient.
107       The default value in the absence of the urey-cubic keyword is zero; i.e., the cubic Urey-Bradley term is omitted.
108       """)
109   public double cubic = DEFAULT_UREY_CUBIC;
110 
111   /**
112    * Quartic coefficient in bond stretch potential.
113    */
114   @FFXProperty(name = "urey-quartic", propertyGroup = LocalGeometryFunctionalForm, defaultValue = "0.0", description = """
115       Sets the value of the quartic term in the Taylor series expansion form of the Urey-Bradley potential energy.
116       The real number modifier gives the value of the coefficient as a multiple of the quadratic coefficient.
117       The default value in the absence of the urey-quartic keyword is zero; i.e., the quartic Urey-Bradley term is omitted.
118       """)
119   public double quartic = DEFAULT_UREY_QUARTIC;
120 
121   /**
122    * A Logger for the UreyBradleyType class.
123    */
124   private static final Logger logger = Logger.getLogger(UreyBradleyType.class.getName());
125   /**
126    * Atom classes that form this Urey-Bradley cross term.
127    */
128   public final int[] atomClasses;
129   /**
130    * Force constant (Kcal/mole/angstroms^2).
131    */
132   public final double forceConstant;
133   /**
134    * Equilibrium 1-3 separation (Angstroms).
135    */
136   public final double distance;
137 
138   /**
139    * UreyBradleyType constructor.
140    *
141    * @param atomClasses   Atom classes.
142    * @param forceConstant Force constant (Kcal/mole/angstroms^2).
143    * @param distance      Equilibrium 1-3 separation (Angstroms).
144    */
145   public UreyBradleyType(int[] atomClasses, double forceConstant, double distance) {
146     super(UREYBRAD, sortKey(atomClasses));
147     this.atomClasses = atomClasses;
148     this.forceConstant = forceConstant;
149     this.distance = distance;
150   }
151 
152   /**
153    * Average two UreyBradleyType objects.
154    *
155    * @param ureyBradleyType1 First UreyBradleyType.
156    * @param ureyBradleyType2 Second UreyBradleyType.
157    * @param atomClasses       Atom classes to use for the average.
158    * @return A new UreyBradleyType with averaged properties, or null if inputs are invalid.
159    */
160   public static UreyBradleyType average(@Nullable UreyBradleyType ureyBradleyType1,
161                                         @Nullable UreyBradleyType ureyBradleyType2,
162                                         @Nullable int[] atomClasses) {
163     if (ureyBradleyType1 == null || ureyBradleyType2 == null || atomClasses == null) {
164       return null;
165     }
166 
167     double forceConstant = (ureyBradleyType1.forceConstant + ureyBradleyType2.forceConstant) / 2.0;
168     double distance = (ureyBradleyType1.distance + ureyBradleyType2.distance) / 2.0;
169 
170     return new UreyBradleyType(atomClasses, forceConstant, distance);
171   }
172 
173   /**
174    * Construct a UreyBradleyType from an input string.
175    *
176    * @param input  The overall input String.
177    * @param tokens The input String tokenized.
178    * @return a UreyBradleyType instance.
179    */
180   public static UreyBradleyType parse(String input, String[] tokens) {
181     if (tokens.length < 5) {
182       logger.log(Level.WARNING, "Invalid UREYBRAD type:\n{0}", input);
183     } else {
184       try {
185         int[] atomClasses = new int[3];
186         atomClasses[0] = parseInt(tokens[1]);
187         atomClasses[1] = parseInt(tokens[2]);
188         atomClasses[2] = parseInt(tokens[3]);
189         double forceConstant = parseDouble(tokens[4]);
190         double distance = parseDouble(tokens[5]);
191         return new UreyBradleyType(atomClasses, forceConstant, distance);
192       } catch (NumberFormatException e) {
193         String message = "Exception parsing UREYBRAD type:\n" + input + "\n";
194         logger.log(Level.SEVERE, message, e);
195       }
196     }
197     return null;
198   }
199 
200   /**
201    * This method sorts the atom classes as: min, c[1], max
202    *
203    * @param c atomClasses
204    * @return lookup key
205    */
206   public static String sortKey(int[] c) {
207     if (c == null || c.length != 3) {
208       return null;
209     }
210     if (c[0] > c[2]) {
211       int temp = c[0];
212       c[0] = c[2];
213       c[2] = temp;
214     }
215 
216     return c[0] + " " + c[1] + " " + c[2];
217   }
218 
219   /**
220    * {@inheritDoc}
221    */
222   @Override
223   public int compare(String key1, String key2) {
224     String[] keys1 = key1.split(" ");
225     String[] keys2 = key2.split(" ");
226     int[] c1 = new int[3];
227     int[] c2 = new int[3];
228     for (int i = 0; i < 3; i++) {
229       c1[i] = parseInt(keys1[i]);
230       c2[i] = parseInt(keys2[i]);
231     }
232     if (c1[1] < c2[1]) {
233       return -1;
234     } else if (c1[1] > c2[1]) {
235       return 1;
236     } else if (c1[0] < c2[0]) {
237       return -1;
238     } else if (c1[0] > c2[0]) {
239       return 1;
240     } else if (c1[2] < c2[2]) {
241       return -1;
242     } else if (c1[2] > c2[2]) {
243       return 1;
244     }
245     return 0;
246   }
247 
248   /**
249    * {@inheritDoc}
250    */
251   @Override
252   public boolean equals(Object o) {
253     if (this == o) {
254       return true;
255     }
256     if (o == null || getClass() != o.getClass()) {
257       return false;
258     }
259     UreyBradleyType ureyBradleyType = (UreyBradleyType) o;
260     return Arrays.equals(atomClasses, ureyBradleyType.atomClasses);
261   }
262 
263   /**
264    * {@inheritDoc}
265    */
266   @Override
267   public int hashCode() {
268     return Arrays.hashCode(atomClasses);
269   }
270 
271   /**
272    * Increment the atom classes by a specified amount.
273    *
274    * @param increment The increment to apply to the atom classes.
275    */
276   public void incrementClasses(int increment) {
277     for (int i = 0; i < atomClasses.length; i++) {
278       atomClasses[i] += increment;
279     }
280     setKey(sortKey(atomClasses));
281   }
282 
283   /**
284    * Remap new atom classes to known internal ones.
285    *
286    * @param typeMap a lookup between new atom types and known atom types.
287    */
288   public void patchClasses(HashMap<AtomType, AtomType> typeMap) {
289     int count = 0;
290     for (AtomType newType : typeMap.keySet()) {
291       for (int atomClass : atomClasses) {
292         if (atomClass == newType.atomClass) {
293           count++;
294         }
295       }
296     }
297     if (count > 0 && count < atomClasses.length) {
298       for (AtomType newType : typeMap.keySet()) {
299         for (int i = 0; i < atomClasses.length; i++) {
300           if (atomClasses[i] == newType.atomClass) {
301             AtomType knownType = typeMap.get(newType);
302             atomClasses[i] = knownType.atomClass;
303           }
304         }
305       }
306       setKey(sortKey(atomClasses));
307     }
308   }
309 
310   /**
311    * {@inheritDoc}
312    *
313    * <p>Nicely formatted Urey-Bradley string.
314    */
315   @Override
316   public String toString() {
317     return format("ureybrad  %5d  %5d  %5d  %6.2f  %7.4f", atomClasses[0], atomClasses[1],
318         atomClasses[2], forceConstant, distance);
319   }
320 
321   /**
322    * Create an AmoebaUreyBradleyForce element.
323    *
324    * @param doc        the Document instance.
325    * @param forceField the ForceField to grab constants from.
326    * @return the AmoebaUreyBradleyForce element.
327    */
328   public static Element getXMLForce(Document doc, ForceField forceField) {
329     Map<String, UreyBradleyType> types = forceField.getUreyBradleyTypes();
330     if (!types.values().isEmpty()) {
331       Element node = doc.createElement("AmoebaUreyBradleyForce");
332       node.setAttribute("cubic", valueOf(forceField.getDouble("urey-cubic", DEFAULT_UREY_CUBIC)));
333       node.setAttribute("quartic", valueOf(forceField.getDouble("urey-quartic", DEFAULT_UREY_QUARTIC)));
334       for (UreyBradleyType ureyBradleyType : types.values()) {
335         node.appendChild(ureyBradleyType.toXML(doc));
336       }
337       return node;
338     }
339     return null;
340   }
341 
342   /**
343    * Write UreyBradleyType to OpenMM XML format.
344    */
345   public Element toXML(Document doc) {
346     Element node = doc.createElement("UreyBradley");
347     node.setAttribute("class1", format("%d", atomClasses[0]));
348     node.setAttribute("class2", format("%d", atomClasses[1]));
349     node.setAttribute("class3", format("%d", atomClasses[2]));
350     // Convert from kcal/mol/A^2 to KJ/mol/nm^2
351     node.setAttribute("k", format("%f", forceConstant * KCAL_TO_KJ / (ANG_TO_NM * ANG_TO_NM)));
352     // Convert from Angstroms to nm
353     node.setAttribute("d", format("%f", distance * ANG_TO_NM));
354     return node;
355   }
356 }