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.IMPROPER;
41  import static ffx.potential.parameters.ForceField.ForceFieldType.TORSION;
42  import static ffx.utilities.KeywordGroup.EnergyUnitConversion;
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  import static java.util.Arrays.copyOf;
48  import static org.apache.commons.math3.util.FastMath.cos;
49  import static org.apache.commons.math3.util.FastMath.sin;
50  import static org.apache.commons.math3.util.FastMath.toRadians;
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 TorsionType class defines a torsional angle.
61   *
62   * @author Michael J. Schnieders
63   * @since 1.0
64   */
65  @FFXKeyword(name = "improper", clazz = String.class, keywordGroup = PotentialFunctionParameter,
66      description = "[4 integers and 2 reals]"
67          + "Provides the values for a single CHARMM-style improper dihedral angle parameter. "
68          + "The integer modifiers give the atom class numbers for the four kinds of atoms involved in the torsion which is to be defined. "
69          + "The real number modifiers give the force constant value for the deviation from the target improper torsional angle, and the target value for the torsional angle, respectively. "
70          + "The default units for the improper force constant are kcal/mole/radian^2, but this can be controlled via the impropunit keyword.")
71  @FFXKeyword(name = "torsion", clazz = String.class, keywordGroup = PotentialFunctionParameter,
72      description = "[4 integers and up to 6 real/real/integer triples] "
73          + "Provides the values for a single torsional angle parameter. "
74          + "The first four integer modifiers give the atom class numbers for the atoms involved in the torsional angle to be defined. "
75          + "Each of the remaining triples of real/real/integer modifiers give the amplitude, phase offset in degrees and periodicity of a particular torsional function term, respectively. "
76          + "Periodicities through 6-fold are allowed for torsional parameters.")
77  public final class TorsionType extends BaseType implements Comparator<String> {
78  
79    private static final Logger logger = Logger.getLogger(TorsionType.class.getName());
80    /** Atom classes that for this Torsion angle. */
81    public final int[] atomClasses;
82    /** Number of terms in the Fourier series. */
83    public final int terms;
84    /** Amplitudes of the Fourier series. */
85    public final double[] amplitude;
86    /** Phases of the Fourier series in degrees. */
87    public final double[] phase;
88    /** Cosine of the phase angle. */
89    public final double[] cosine;
90    /** Sine of the phase angle. */
91    public final double[] sine;
92    /** Periodicity of the Fourier series. */
93    public final int[] periodicity;
94    /** The torsion mode in use. */
95    private final TorsionMode torsionMode;
96  
97    public static final double DEFAULT_TORSION_UNIT = 1.0;
98    /** Unit conversion. */
99    @FFXKeyword(name = "torsionunit", keywordGroup = EnergyUnitConversion, defaultValue = "1.0",
100       description = "Sets the scale factor needed to convert the energy value computed by the torsional angle potential into units of kcal/mole. "
101           + "The correct value is force field dependent and typically provided in the header of the master force field parameter file.")
102   public double torsionUnit = DEFAULT_TORSION_UNIT;
103 
104   /**
105    * TorsionType Constructor.
106    *
107    * @param atomClasses Atom classes.
108    * @param amplitude Amplitudes of the Fourier series.
109    * @param phase Phases of the Fourier series in degrees.
110    * @param periodicity Periodicity of the Fourier series.
111    */
112   public TorsionType(int[] atomClasses, double[] amplitude, double[] phase, int[] periodicity) {
113     this(atomClasses, amplitude, phase, periodicity, TorsionMode.NORMAL);
114   }
115 
116   /**
117    * TorsionType Constructor.
118    *
119    * @param atomClasses Atom classes.
120    * @param amplitude Amplitudes of the Fourier series.
121    * @param phase Phases of the Fourier series in degrees.
122    * @param periodicity Periodicity of the Fourier series.
123    * @param torsionMode Define the TorsionMode for this TorsionType.
124    */
125   public TorsionType(
126       int[] atomClasses,
127       double[] amplitude,
128       double[] phase,
129       int[] periodicity,
130       TorsionMode torsionMode) {
131     super(TORSION, sortKey(atomClasses));
132     this.atomClasses = atomClasses;
133     int max = 1;
134     for (int i1 : periodicity) {
135       if (i1 > max) {
136         max = i1;
137       }
138     }
139     terms = max;
140     this.amplitude = new double[terms];
141     this.phase = new double[terms];
142     this.periodicity = new int[terms];
143     for (int i = 0; i < terms; i++) {
144       this.periodicity[i] = i + 1;
145     }
146     for (int i = 0; i < amplitude.length; i++) {
147       int j = periodicity[i] - 1;
148       if (j >= 0 && j < terms) {
149         this.amplitude[j] = amplitude[i];
150         this.phase[j] = phase[i];
151       }
152     }
153 
154     cosine = new double[terms];
155     sine = new double[terms];
156     for (int i = 0; i < terms; i++) {
157       double angle = toRadians(this.phase[i]);
158       cosine[i] = cos(angle);
159       sine[i] = sin(angle);
160     }
161 
162     this.torsionMode = torsionMode;
163     if (torsionMode == TorsionMode.IMPROPER) {
164       forceFieldType = IMPROPER;
165     }
166   }
167 
168   /**
169    * average.
170    *
171    * @param torsionType1 a {@link ffx.potential.parameters.TorsionType} object.
172    * @param torsionType2 a {@link ffx.potential.parameters.TorsionType} object.
173    * @param atomClasses an array of {@link int} objects.
174    * @return a {@link ffx.potential.parameters.TorsionType} object.
175    */
176   public static TorsionType average(
177       TorsionType torsionType1, TorsionType torsionType2, int[] atomClasses) {
178     if (torsionType1 == null || torsionType2 == null || atomClasses == null) {
179       return null;
180     }
181     int len = torsionType1.amplitude.length;
182     if (len != torsionType2.amplitude.length) {
183       return null;
184     }
185     double[] amplitude = new double[len];
186     double[] phase = new double[len];
187     int[] periodicity = new int[len];
188     for (int i = 0; i < len; i++) {
189       amplitude[i] = (torsionType1.amplitude[i] + torsionType2.amplitude[i]) / 2.0;
190       phase[i] = (torsionType1.phase[i] + torsionType2.phase[i]) / 2.0;
191       periodicity[i] = (torsionType1.periodicity[i] + torsionType2.periodicity[i]) / 2;
192     }
193 
194     return new TorsionType(atomClasses, amplitude, phase, periodicity);
195   }
196 
197   /**
198    * Construct a TorsionType from an input string.
199    *
200    * @param input The overall input String.
201    * @param tokens The input String tokenized.
202    * @return a TorsionType instance.
203    */
204   public static TorsionType parse(String input, String[] tokens) {
205     if (tokens.length < 5) {
206       logger.log(Level.WARNING, "Invalid TORSION type:\n{0}", input);
207     } else {
208       try {
209         int[] atomClasses = new int[4];
210         atomClasses[0] = parseInt(tokens[1]);
211         atomClasses[1] = parseInt(tokens[2]);
212         atomClasses[2] = parseInt(tokens[3]);
213         atomClasses[3] = parseInt(tokens[4]);
214         int terms = (tokens.length - 5) / 3;
215         double[] amplitude = new double[terms];
216         double[] phase = new double[terms];
217         int[] periodicity = new int[terms];
218         int index = 5;
219         for (int i = 0; i < terms; i++) {
220           amplitude[i] = parseDouble(tokens[index++]);
221           phase[i] = parseDouble(tokens[index++]);
222           periodicity[i] = parseInt(tokens[index++]);
223         }
224         return new TorsionType(atomClasses, amplitude, phase, periodicity);
225       } catch (NumberFormatException e) {
226         String message = "NumberFormatException parsing TORSION type:\n" + input + "\n";
227         logger.log(Level.SEVERE, message, e);
228       } catch (Exception e) {
229         String message = "Exception parsing TORSION type:\n" + input + "\n";
230         logger.log(Level.SEVERE, message, e);
231       }
232     }
233     return null;
234   }
235 
236   /**
237    * Construct a TorsionType with TorsionMode.IMPROPER from an input string.
238    *
239    * @param input The overall input String.
240    * @param tokens The input String tokenized.
241    * @return a TorsionType instance.
242    */
243   public static TorsionType parseImproper(String input, String[] tokens) {
244     if (tokens.length < 5) {
245       logger.log(Level.WARNING, "Invalid IMPROPER type:\n{0}", input);
246     } else {
247       try {
248         int[] atomClasses = new int[4];
249         atomClasses[0] = parseInt(tokens[1]);
250         atomClasses[1] = parseInt(tokens[2]);
251         atomClasses[2] = parseInt(tokens[3]);
252         atomClasses[3] = parseInt(tokens[4]);
253         double[] amplitude = new double[1];
254         double[] phase = new double[1];
255         int[] periodicity = new int[1];
256         int index = 5;
257         amplitude[0] = parseDouble(tokens[index++]);
258         phase[0] = parseDouble(tokens[index]);
259         periodicity[0] = 1;
260         return new TorsionType(atomClasses, amplitude, phase, periodicity, TorsionMode.IMPROPER);
261       } catch (NumberFormatException e) {
262         String message = "Exception parsing IMPROPER type:\n" + input + "\n";
263         logger.log(Level.SEVERE, message, e);
264       }
265     }
266     return null;
267   }
268 
269   /**
270    * This method sorts the atom classes for the torsion.
271    *
272    * @param c atomClasses
273    * @return lookup key
274    * @since 1.0
275    */
276   public static String sortKey(int[] c) {
277     if (c == null || c.length != 4) {
278       return null;
279     }
280     if (c[1] < c[2]) {
281       // Do nothing.
282     } else if (c[2] < c[1]) {
283       // Reverse the order.
284       int temp = c[0];
285       c[0] = c[3];
286       c[3] = temp;
287       temp = c[1];
288       c[1] = c[2];
289       c[2] = temp;
290     } else if (c[1] == c[2]) {
291       if (c[0] > c[3]) {
292         // Reverse the order.
293         int temp = c[0];
294         c[0] = c[3];
295         c[3] = temp;
296       }
297     } else if (c[0] <= c[3]) {
298       // Do nothing.
299     } else {
300       // Reverse the order.
301       int temp = c[0];
302       c[0] = c[3];
303       c[3] = temp;
304       temp = c[1];
305       c[1] = c[2];
306       c[2] = temp;
307     }
308 
309     return c[0] + " " + c[1] + " " + c[2] + " " + c[3];
310   }
311 
312   /**
313    * {@inheritDoc}
314    *
315    * @since 1.0
316    */
317   @Override
318   public int compare(String s1, String s2) {
319     String[] keys1 = s1.split(" ");
320     String[] keys2 = s2.split(" ");
321     int[] c1 = new int[4];
322     int[] c2 = new int[4];
323 
324     for (int i = 0; i < 4; i++) {
325       c1[i] = parseInt(keys1[i]);
326       c2[i] = parseInt(keys2[i]);
327     }
328 
329     if (c1[1] < c2[1]) {
330       return -1;
331     } else if (c1[1] > c2[1]) {
332       return 1;
333     } else if (c1[2] < c2[2]) {
334       return -1;
335     } else if (c1[2] > c2[2]) {
336       return 1;
337     } else if (c1[0] < c2[0]) {
338       return -1;
339     } else if (c1[0] > c2[0]) {
340       return 1;
341     } else if (c1[3] < c2[3]) {
342       return -1;
343     } else if (c1[3] > c2[3]) {
344       return 1;
345     }
346 
347     return 0;
348   }
349 
350   /**
351    * {@inheritDoc}
352    *
353    * <p>Override the default <code>equals</code> method.
354    *
355    * @since 1.0
356    */
357   @Override
358   public boolean equals(Object o) {
359     if (this == o) {
360       return true;
361     }
362     if (o == null || getClass() != o.getClass()) {
363       return false;
364     }
365     TorsionType torsionType = (TorsionType) o;
366     return Arrays.equals(atomClasses, torsionType.atomClasses);
367   }
368 
369   /**
370    * {@inheritDoc}
371    *
372    * <p>Implementation of the <code>hashCode</code> method.
373    *
374    * @since 1.0
375    */
376   @Override
377   public int hashCode() {
378     return Arrays.hashCode(atomClasses);
379   }
380 
381   /**
382    * incrementClasses
383    *
384    * @param increment a int.
385    */
386   public void incrementClasses(int increment) {
387     for (int i = 0; i < atomClasses.length; i++) {
388       if (atomClasses[i] != 0) {
389         atomClasses[i] += increment;
390       }
391     }
392     setKey(sortKey(atomClasses));
393   }
394 
395   /**
396    * Remap new atom classes to known internal ones.
397    *
398    * @param typeMap a lookup between new atom types and known atom types.
399    * @return a {@link ffx.potential.parameters.TorsionType} object.
400    */
401   public TorsionType patchClasses(HashMap<AtomType, AtomType> typeMap) {
402     int count = 0;
403     int len = atomClasses.length;
404 
405     // Look for new TorsionTypes that contain 1 to 3 mapped atom classes.
406     for (AtomType newType : typeMap.keySet()) {
407       for (int atomClass : atomClasses) {
408         if (atomClass == newType.atomClass) {
409           count++;
410         }
411       }
412     }
413 
414     // If found, create a new TorsionType that bridges to known classes.
415     if (count == 1 || count == 2 || count == 3) {
416       int[] newClasses = copyOf(atomClasses, len);
417       for (AtomType newType : typeMap.keySet()) {
418         for (int i = 0; i < len; i++) {
419           if (atomClasses[i] == newType.atomClass) {
420             AtomType knownType = typeMap.get(newType);
421             newClasses[i] = knownType.atomClass;
422           }
423         }
424       }
425       return new TorsionType(newClasses, amplitude, phase, periodicity);
426     }
427     return null;
428   }
429 
430   /**
431    * {@inheritDoc}
432    *
433    * <p>Nicely formatted Torsion angle.
434    *
435    * @since 1.0
436    */
437   @Override
438   public String toString() {
439     StringBuilder torsionBuffer;
440     if (torsionMode == TorsionMode.IMPROPER) {
441       torsionBuffer = new StringBuilder("improper");
442     } else {
443       torsionBuffer = new StringBuilder("torsion");
444     }
445     for (int i : atomClasses) {
446       torsionBuffer.append(format(" %5d", i));
447     }
448 
449     boolean nonZero = false;
450     for (double v : amplitude) {
451       if (v != 0.0) {
452         nonZero = true;
453         break;
454       }
455     }
456 
457     for (int i = 0; i < amplitude.length; i++) {
458       if (amplitude[i] == 0.0 && (i > 0 || nonZero)) {
459         continue;
460       }
461       if (torsionMode == TorsionMode.NORMAL) {
462         torsionBuffer.append(format(" %7.3f %7.3f %1d", amplitude[i], phase[i], periodicity[i]));
463       } else {
464         torsionBuffer.append(format(" %7.3f %7.3f", amplitude[i], phase[i]));
465       }
466     }
467     return torsionBuffer.toString();
468   }
469 
470   /**
471    * setScaleFactor.
472    *
473    * @param scale a double.
474    */
475   void setScaleFactor(double scale) {
476     for (int i = 0; i < amplitude.length; i++) {
477       amplitude[i] *= scale;
478     }
479   }
480 
481   /** Torsion modes include Normal or In-Plane */
482   public enum TorsionMode {
483     NORMAL,
484     IMPROPER
485   }
486 }