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 java.io.BufferedReader;
45  import java.io.IOException;
46  import java.util.Arrays;
47  import java.util.Comparator;
48  import java.util.HashMap;
49  import java.util.Map;
50  import java.util.logging.Level;
51  import java.util.logging.Logger;
52  
53  import static ffx.potential.parameters.ForceField.ForceFieldType.TORTORS;
54  import static ffx.utilities.Constants.KCAL_TO_KJ;
55  import static ffx.utilities.PropertyGroup.EnergyUnitConversion;
56  import static ffx.utilities.PropertyGroup.PotentialFunctionParameter;
57  import static java.lang.Double.parseDouble;
58  import static java.lang.Integer.parseInt;
59  import static java.lang.String.format;
60  import static java.lang.System.arraycopy;
61  import static java.util.Arrays.sort;
62  import static org.apache.commons.math3.util.FastMath.abs;
63  
64  /**
65   * The TorsionTorsionType class defines a Torsion-Torsion spline.
66   *
67   * @author Michael J. Schnieders
68   * @since 1.0
69   */
70  @FFXProperty(name = "tortors", clazz = String[].class, propertyGroup = PotentialFunctionParameter, description = """
71      [7 integers, then multiple lines of 2 integers and 1 real]
72      Provides the values for a single torsion-torsion parameter.
73      The first five integer modifiers give the atom class numbers for the atoms involved in the two adjacent torsional angles to be defined.
74      The last two integer modifiers contain the number of data grid points that lie along each axis of the torsion-torsion map.
75      For example, this value will be 13 for a 30 degree torsional angle spacing, i.e., 360/30 = 12, but 13
76      values are required since data values for -180 and +180 degrees must both be supplied.
77      The subsequent lines contain the torsion-torsion map data as the integer values in degrees of each
78      torsional angle and the target energy value in kcal/mole.
79      """)
80  public final class TorsionTorsionType extends BaseType implements Comparator<String> {
81  
82    /**
83     * Default units to convert Torsion-Torsion energy to kcal/mole.
84     */
85    public static final double DEFAULT_TORTOR_UNIT = 1.0;
86    /**
87     * Convert Torsion-Torsion energy to kcal/mole.
88     */
89    @FFXProperty(name = "tortorunit", propertyGroup = EnergyUnitConversion, defaultValue = "1.0", description = """
90        Sets the scale factor needed to convert the energy value computed by the torsion-torsion potential into units of kcal/mole.
91        The correct value is force field dependent and typically provided in the header of the master force field parameter file.
92        """)
93    public double torTorUnit = DEFAULT_TORTOR_UNIT;
94  
95    private static final Logger logger = Logger.getLogger(TorsionTorsionType.class.getName());
96    /**
97     * Atom classes that form this Torsion-Torsion type.
98     */
99    public final int[] atomClasses;
100   /**
101    * Energy values.
102    */
103   public final double[] energy;
104   /**
105    * Number of points along x.
106    */
107   public final int nx;
108   /**
109    * Number of point along y.
110    */
111   public final int ny;
112   /**
113    * Torsion values along x.
114    */
115   public final double[] tx;
116   /**
117    * Torsion values along y.
118    */
119   public final double[] ty;
120   /**
121    * First derivative along x.
122    */
123   public final double[] dx;
124   /**
125    * First derivative along y.
126    */
127   public final double[] dy;
128   /**
129    * Second derivatives.
130    */
131   public final double[] dxy;
132   /**
133    * Grid points.
134    */
135   private final int[] gridPoints;
136 
137   /**
138    * Constructor for TorsionTorsionType.
139    *
140    * @param atomClasses an array of int.
141    * @param gridPoints  an array of int.
142    * @param torsion1    an array of double.
143    * @param torsion2    an array of double.
144    * @param energy      an array of double.
145    */
146   public TorsionTorsionType(int[] atomClasses, int[] gridPoints, double[] torsion1,
147                             double[] torsion2, double[] energy) {
148     super(TORTORS, sortKey(atomClasses));
149     this.atomClasses = atomClasses;
150     nx = gridPoints[0];
151     ny = gridPoints[1];
152     if (nx != ny) {
153       logger.severe("Untested TORTOR parameters: nx != ny: " + nx + ", " + ny);
154     }
155 
156     this.energy = energy;
157     this.gridPoints = gridPoints;
158     tx = new double[nx];
159     ty = new double[ny];
160     dx = new double[nx * ny];
161     dy = new double[nx * ny];
162     dxy = new double[nx * ny];
163     sort(torsion1);
164     sort(torsion2);
165     tx[0] = torsion1[0];
166     ty[0] = torsion2[0];
167     int j1 = 1;
168     int j2 = 1;
169     for (int i = 1; i < nx; i++) {
170       while (torsion1[j1] == tx[i - 1]) {
171         j1++;
172       }
173       while (torsion2[j2] == ty[i - 1]) {
174         j2++;
175       }
176       tx[i] = torsion1[j1];
177       ty[i] = torsion2[j2];
178     }
179 
180     // Check for cyclic energy.
181     boolean isCyclic = true;
182     double eps = 0.0001;
183     if (abs(tx[0] - tx[nx - 1]) - 360.0 > eps) {
184       isCyclic = false;
185       if (logger.isLoggable(Level.FINEST)) {
186         logger.finest(" tortor is aperiodic: " + tx[0] + ", " + tx[nx - 1]);
187       }
188     }
189     if (isCyclic) {
190       for (int i = 0; i < ny; i++) {
191         int k = i * nx;
192         if (abs(energy[k] - energy[k + nx - 1]) > eps) {
193           isCyclic = false;
194           if (logger.isLoggable(Level.FINEST)) {
195             logger.finest(" tortor is apreriodic: " + k + ", " + (k + nx - 1) + ": " + abs(
196                 energy[k] - energy[k + nx - 1]));
197           }
198           break;
199         }
200       }
201     }
202     if (isCyclic) {
203       int k = (ny - 1) * nx;
204       for (int i = 0; i < nx; i++) {
205         if (abs(energy[i] - energy[i + k]) > eps) {
206           if (logger.isLoggable(Level.FINEST)) {
207             logger.fine(
208                 " tortor is aperiodic: " + i + ", " + i + k + ": " + abs(energy[i] - energy[i + k]));
209           }
210           isCyclic = false;
211           break;
212         }
213       }
214     }
215 
216     boolean cyclic = isCyclic;
217     double[] tmp1 = new double[nx];
218     double[] tmp2 = new double[nx];
219     double[] tmp3 = new double[nx];
220     double[] tmp4 = new double[nx];
221     double[] tmp5 = new double[nx];
222     double[] tmp6 = new double[nx];
223     double[] tmp7 = new double[nx];
224     double[] bs = new double[nx];
225     double[] cs = new double[nx];
226     double[] ds = new double[nx];
227 
228     // Spline fit the derivatives about the first torsion.
229     arraycopy(tx, 0, tmp1, 0, nx);
230 
231     int m = 0;
232     for (int j = 0; j < ny; j++) {
233       arraycopy(energy, m, tmp2, 0, nx);
234       if (cyclic) {
235         cspline(nx - 1, tmp1, tmp2, bs, cs, ds, tmp3, tmp4, tmp5, tmp6, tmp7);
236       } else {
237         nspline(nx - 1, tmp1, tmp2, 0.0, 0.0, bs, cs, tmp3, tmp4, tmp5, tmp6, tmp7);
238       }
239       arraycopy(bs, 0, dx, m, nx);
240       m = m + nx;
241     }
242 
243     // Spline fit the derivatives about the second torsion.
244     arraycopy(ty, 0, tmp1, 0, ny);
245 
246     m = 0;
247     for (int j = 0; j < nx; j++) {
248       for (int k = 0; k < ny; k++) {
249         tmp2[k] = energy[m + k * nx];
250       }
251       if (cyclic) {
252         cspline(ny - 1, tmp1, tmp2, bs, cs, ds, tmp3, tmp4, tmp5, tmp6, tmp7);
253       } else {
254         nspline(ny - 1, tmp1, tmp2, 0.0, 0.0, bs, cs, tmp3, tmp4, tmp5, tmp6, tmp7);
255       }
256       for (int k = 0; k < ny; k++) {
257         dy[m + k * nx] = bs[k];
258       }
259       m = m + 1;
260     }
261 
262     // Spline fit the cross derivatives about both torsions.
263     m = 0;
264     for (int j = 0; j < nx; j++) {
265       for (int k = 0; k < ny; k++) {
266         tmp2[k] = dx[m + k * nx];
267       }
268       if (cyclic) {
269         cspline(ny - 1, tmp1, tmp2, bs, cs, ds, tmp3, tmp4, tmp5, tmp6, tmp7);
270       } else {
271         nspline(ny - 1, tmp1, tmp2, 0.0, 0.0, bs, cs, tmp3, tmp4, tmp5, tmp6, tmp7);
272       }
273       for (int k = 0; k < ny; k++) {
274         dxy[m + k * nx] = bs[k];
275       }
276       m = m + 1;
277     }
278   }
279 
280   /**
281    * average.
282    *
283    * @param torsionTorsionType1 a {@link ffx.potential.parameters.TorsionTorsionType} object.
284    * @param torsionTorsionType2 a {@link ffx.potential.parameters.TorsionTorsionType} object.
285    * @param atomClasses         an array of {@link int} objects.
286    * @return a {@link ffx.potential.parameters.TorsionTorsionType} object.
287    */
288   public static TorsionTorsionType average(TorsionTorsionType torsionTorsionType1,
289                                            TorsionTorsionType torsionTorsionType2, int[] atomClasses) {
290     if (torsionTorsionType1 == null || torsionTorsionType2 == null || atomClasses == null) {
291       return null;
292     }
293     return null;
294   }
295 
296   /**
297    * Construct a TorsionTorsionType from multiple input lines.
298    *
299    * @param input  The overall input String.
300    * @param tokens The input String tokenized.
301    * @param br     a BufferedReader instance.
302    * @return a TorsionTorsionType instance.
303    */
304   public static TorsionTorsionType parse(String input, String[] tokens, BufferedReader br) {
305     if (tokens.length < 8) {
306       logger.log(Level.WARNING, "Invalid TORTORS type:\n{0}", input);
307     } else {
308       try {
309         int[] atomClasses = new int[5];
310         for (int i = 0; i < 5; i++) {
311           atomClasses[i] = parseInt(tokens[i + 1]);
312         }
313         int[] gridPoints = new int[2];
314         gridPoints[0] = parseInt(tokens[6]);
315         gridPoints[1] = parseInt(tokens[7]);
316         int points = gridPoints[0] * gridPoints[1];
317         double[] torsion1 = new double[points];
318         double[] torsion2 = new double[points];
319         double[] energy = new double[points];
320         for (int i = 0; i < points; i++) {
321           input = br.readLine();
322           tokens = input.trim().split(" +");
323           if (tokens.length == 3) {
324             torsion1[i] = parseDouble(tokens[0]);
325             torsion2[i] = parseDouble(tokens[1]);
326             energy[i] = parseDouble(tokens[2]);
327           } else if (tokens.length == 6) {
328             // First entry.
329             torsion1[i] = parseDouble(tokens[0]);
330             torsion2[i] = parseDouble(tokens[1]);
331             energy[i] = parseDouble(tokens[2]);
332             // Second entry.
333             torsion1[i + 1] = parseDouble(tokens[3]);
334             torsion2[i + 1] = parseDouble(tokens[4]);
335             energy[i + 1] = parseDouble(tokens[5]);
336             // Increment i by 1 since we read 2 values.
337             i += 1;
338           } else if (tokens.length == 9) {
339             // First entry.
340             torsion1[i] = parseDouble(tokens[0]);
341             torsion2[i] = parseDouble(tokens[1]);
342             energy[i] = parseDouble(tokens[2]);
343             // Second entry.
344             torsion1[i + 1] = parseDouble(tokens[3]);
345             torsion2[i + 1] = parseDouble(tokens[4]);
346             energy[i + 1] = parseDouble(tokens[5]);
347             // Third entry.
348             torsion1[i + 2] = parseDouble(tokens[6]);
349             torsion2[i + 2] = parseDouble(tokens[7]);
350             energy[i + 2] = parseDouble(tokens[8]);
351             // Increment i by 2 since we read 3 values.
352             i += 2;
353           } else {
354             logger.log(Level.WARNING, "Invalid TORTORS type:\n{0}", input);
355             return null;
356           }
357         }
358         return new TorsionTorsionType(atomClasses, gridPoints, torsion1, torsion2, energy);
359       } catch (NumberFormatException | IOException e) {
360         String message = "Exception parsing TORTORS type:\n" + input + "\n";
361         logger.log(Level.SEVERE, message, e);
362       }
363     }
364     return null;
365   }
366 
367   /**
368    * Construct a TorsionTorsionType from a single input line.
369    *
370    * @param input  The overall input String.
371    * @param tokens The input String tokenized.
372    * @return a TorsionTorsionType instance.
373    */
374   public static TorsionTorsionType parse(String input, String[] tokens) {
375     if (tokens.length < 8) {
376       logger.log(Level.WARNING, "Invalid TORTORS type:\n{0}", input);
377     } else {
378       try {
379         int[] atomClasses = new int[5];
380         for (int i = 0; i < 5; i++) {
381           atomClasses[i] = parseInt(tokens[i + 1]);
382         }
383         int[] gridPoints = new int[2];
384         gridPoints[0] = parseInt(tokens[6]);
385         gridPoints[1] = parseInt(tokens[7]);
386         int points = gridPoints[0] * gridPoints[1];
387         int numTokens = points * 3 + 8;
388         if (tokens.length < numTokens) {
389           logger.log(Level.WARNING, "Invalid TORTORS type:\n{0}", input);
390           return null;
391         }
392         double[] torsion1 = new double[points];
393         double[] torsion2 = new double[points];
394         double[] energy = new double[points];
395         int index = 8;
396         for (int i = 0; i < points; i++) {
397           torsion1[i] = parseDouble(tokens[index++]);
398           torsion2[i] = parseDouble(tokens[index++]);
399           energy[i] = parseDouble(tokens[index++]);
400         }
401         return new TorsionTorsionType(atomClasses, gridPoints, torsion1, torsion2, energy);
402       } catch (NumberFormatException e) {
403         String message = "Exception parsing TORTORS type:\n" + input + "\n";
404         logger.log(Level.SEVERE, message, e);
405       }
406     }
407     return null;
408   }
409 
410   /**
411    * Reversed key for the Torsion-Torsion lookup.
412    *
413    * @param c atomClasses
414    * @return lookup key
415    */
416   public static String reverseKey(int[] c) {
417     return c[4] + " " + c[3] + " " + c[2] + " " + c[1] + " " + c[0];
418   }
419 
420   /**
421    * No sorting is done for the Torsion-Torsion lookup.
422    *
423    * @param c atomClasses
424    * @return lookup key
425    */
426   public static String sortKey(int[] c) {
427     return c[0] + " " + c[1] + " " + c[2] + " " + c[3] + " " + c[4];
428   }
429 
430   /**
431    * Solves a system of linear equations for a cyclically tridiagonal, symmetric, positive definite
432    * matrix.
433    *
434    * @param n
435    * @param dm
436    * @param du
437    * @param rs
438    * @param c
439    * @return positive or negative 1.
440    */
441   private static int cytsy(int n, double[] dm, double[] du, double[] cr, double[] rs, double[] c) {
442     if (n < 3) {
443       return -2;
444     }
445 
446     // Factorization of the input matrix.
447     if (cytsyp(n, dm, du, cr) == 1) {
448       // Update and back-substitute as necessary.
449       cytsys(n, dm, du, cr, rs, c);
450       return 1;
451     }
452     return -1;
453   }
454 
455   /**
456    * Finds the Cholesky factors of a cyclically tridiagonal symmetric, positive definite matrix given
457    * by two vectors.
458    *
459    * @param n
460    * @param dm
461    * @param du
462    * @param cr
463    * @return an integer.
464    */
465   private static int cytsyp(int n, double[] dm, double[] du, double[] cr) {
466     double eps = 0.00000001;
467     // Check for n < 3.
468     if (n < 3) {
469       return -2;
470     }
471 
472     // Check if the matrix is positive definite.
473     double row = abs(dm[1]) + abs(du[1]) + abs(du[n]);
474     if (row == 0.0) {
475       return 0;
476     }
477     double d = 1.0 / row;
478     if (dm[1] < 0.0) {
479       return -1;
480     } else if (abs(dm[1]) * d < eps) {
481       return 0;
482     }
483 
484     // Factoring while checking for a positive definite and strong non-singular matrix.
485     double temp1 = du[1];
486     double temp2 = 0.0;
487     du[1] = du[1] / dm[1];
488     cr[1] = du[n] / dm[1];
489     for (int i = 2; i < n; i++) {
490       row = abs(dm[i]) + abs(du[i]) + abs(temp1);
491       if (row == 0.0) {
492         return 0;
493       }
494       d = 1.0 / row;
495       dm[i] = dm[i] - temp1 * du[i - 1];
496       if (dm[i] < 0.0) {
497         return -1;
498       } else if (abs(dm[i]) * d < eps) {
499         return 0;
500       }
501       if (i < (n - 1)) {
502         cr[i] = -temp1 * cr[i - 1] / dm[i];
503         temp1 = du[i];
504         du[i] = du[i] / dm[i];
505       } else {
506         temp2 = du[i];
507         du[i] = (du[i] - temp1 * cr[i - 1]) / dm[i];
508       }
509     }
510     row = abs(du[n]) + abs(dm[n]) + abs(temp2);
511     if (row == 0.0) {
512       return 0;
513     }
514     d = 1.0 / row;
515     dm[n] = dm[n] - dm[n - 1] * du[n - 1] * du[n - 1];
516     temp1 = 0.0;
517     for (int i = 1; i < n - 1; i++) {
518       temp1 = temp1 + dm[i] * cr[i] * cr[i];
519     }
520     dm[n] = dm[n] - temp1;
521     if (dm[n] < 0.0) {
522       return -1;
523     } else if (abs(dm[n]) * d < eps) {
524       return 0;
525     }
526     return 1;
527   }
528 
529   /**
530    * Solves a cyclically tridiagonal linear system given the Cholesky factors.
531    *
532    * @param n
533    * @param dm
534    * @param du
535    * @param cr
536    * @param rs
537    * @param c
538    */
539   private static void cytsys(int n, double[] dm, double[] du, double[] cr, double[] rs, double[] c) {
540     // Updating phase.
541     double temp = rs[1];
542     rs[1] = temp / dm[1];
543     double sum = cr[1] * temp;
544     for (int i = 2; i < n; i++) {
545       temp = rs[i] - du[i - 1] * temp;
546       rs[i] = temp / dm[i];
547       if (i != (n - 1)) {
548         sum = sum + cr[i] * temp;
549       }
550     }
551     temp = rs[n] - du[n - 1] * temp;
552     temp = temp - sum;
553     rs[n] = temp / dm[n];
554 
555     // Back-substitution phase.
556     c[n] = rs[n];
557     c[n - 1] = rs[n - 1] - du[n - 1] * c[n];
558     for (int i = n - 2; i >= 1; i--) {
559       c[i] = rs[i] - du[i] * c[i + 1] - cr[i] * c[n];
560     }
561   }
562 
563   /**
564    * {@inheritDoc}
565    */
566   @Override
567   public int compare(String key1, String key2) {
568     String[] keys1 = key1.split(" ");
569     String[] keys2 = key2.split(" ");
570     int c1 = parseInt(keys1[2]);
571     int c2 = parseInt(keys2[2]);
572     if (c1 < c2) {
573       return -1;
574     } else if (c1 > c2) {
575       return 1;
576     }
577     return 0;
578   }
579 
580   /**
581    * {@inheritDoc}
582    */
583   @Override
584   public boolean equals(Object o) {
585     if (this == o) {
586       return true;
587     }
588     if (o == null || getClass() != o.getClass()) {
589       return false;
590     }
591     TorsionTorsionType torsionTorsionType = (TorsionTorsionType) o;
592     return Arrays.equals(atomClasses, torsionTorsionType.atomClasses);
593   }
594 
595   /**
596    * {@inheritDoc}
597    */
598   @Override
599   public int hashCode() {
600     return Arrays.hashCode(atomClasses);
601   }
602 
603   /**
604    * Increment the atom classes by a value.
605    *
606    * @param increment The increment to add to the atom classes.
607    */
608   public void incrementClasses(int increment) {
609     for (int i = 0; i < atomClasses.length; i++) {
610       atomClasses[i] += increment;
611     }
612     setKey(sortKey(atomClasses));
613   }
614 
615   /**
616    * Remap new atom classes to known internal ones.
617    *
618    * @param typeMap a lookup between new atom types and known atom types.
619    */
620   public void patchClasses(HashMap<AtomType, AtomType> typeMap) {
621     int count = 0;
622     for (AtomType newType : typeMap.keySet()) {
623 
624       for (int atomClass : atomClasses) {
625         if (atomClass == newType.atomClass) {
626           count++;
627         }
628       }
629     }
630     if (count > 0 && count < atomClasses.length) {
631       for (AtomType newType : typeMap.keySet()) {
632         for (int i = 0; i < atomClasses.length; i++) {
633           if (atomClasses[i] == newType.atomClass) {
634             AtomType knownType = typeMap.get(newType);
635             atomClasses[i] = knownType.atomClass;
636           }
637         }
638       }
639       setKey(sortKey(atomClasses));
640     }
641   }
642 
643   /**
644    * {@inheritDoc}
645    *
646    * <p>Nicely formatted torsion-torsion type.
647    */
648   @Override
649   public String toString() {
650     StringBuilder tortorBuffer = new StringBuilder("tortors");
651     for (int i : atomClasses) {
652       tortorBuffer.append(format("  %5d", i));
653     }
654     tortorBuffer.append(format("  %2d  %2d", gridPoints[0], gridPoints[1]));
655     for (int i = 0; i < energy.length; i++) {
656       // Emit one entry.
657       int nxi = i % nx;
658       int nyi = i / ny;
659       tortorBuffer.append(format(" \\\n  % 6.1f  % 6.1f  % 8.5f", tx[nxi], ty[nyi], energy[i]));
660       // Emit a second entry if available.
661       if (i < energy.length - 1) {
662         i++;
663         nxi = i % nx;
664         nyi = i / ny;
665         tortorBuffer.append(format("  % 6.1f  % 6.1f  % 8.5f", tx[nxi], ty[nyi], energy[i]));
666       }
667       // Emit a third entry if available.
668       if (i < energy.length - 1) {
669         i++;
670         nxi = i % nx;
671         nyi = i / ny;
672         tortorBuffer.append(format("  % 6.1f  % 6.1f  % 8.5f", tx[nxi], ty[nyi], energy[i]));
673       }
674     }
675     return tortorBuffer.toString();
676   }
677 
678 
679   /**
680    * Create an AmoebaTorsionTorsionForce Element.
681    *
682    * @param doc        the Document instance.
683    * @param forceField the ForceField instance to grab constants from.
684    * @return the AmoebaTorsionTorsionForce Element.
685    */
686   public static Element getXMLForce(Document doc, ForceField forceField) {
687     Map<String, TorsionTorsionType> types = forceField.getTorsionTorsionTypes();
688     if (!types.values().isEmpty()) {
689       Element node = doc.createElement("AmoebaTorsionTorsionForce");
690       int i = 0;
691       for (TorsionTorsionType torsionTorsionType : types.values()) {
692         torsionTorsionType.toXML(doc, node, i);
693         i++;
694       }
695       return node;
696     }
697     return null;
698   }
699 
700   /**
701    * Write TorsionTorsionType to OpenMM XML format.
702    *
703    * @param doc        the Document instance.
704    * @param torTorNode the Element to attach TorsionTorsion and TorsionTorsionGrid to.
705    * @param label      the label used to match TorsionTorsion and TorsionTorsionGrid.
706    */
707   public void toXML(Document doc, Element torTorNode, int label) {
708     Element tortors = doc.createElement("TorsionTorsion");
709     Element ttGrid = doc.createElement("TorsionTorsionGrid");
710 
711     // Set TorsionTorsion node
712     tortors.setAttribute("class1", format("%d", atomClasses[0]));
713     tortors.setAttribute("class2", format("%d", atomClasses[1]));
714     tortors.setAttribute("class3", format("%d", atomClasses[2]));
715     tortors.setAttribute("class4", format("%d", atomClasses[3]));
716     tortors.setAttribute("class5", format("%d", atomClasses[4]));
717     tortors.setAttribute("nx", format("%d", gridPoints[0]));
718     tortors.setAttribute("ny", format("%d", gridPoints[1]));
719 
720     // Set TorsionTorsionGrid
721     ttGrid.setAttribute("nx", format("%d", gridPoints[0]));
722     ttGrid.setAttribute("ny", format("%d", gridPoints[1]));
723     for (int i = 0; i < energy.length; i++) {
724       int nxi = i % nx;
725       int nyi = i / ny;
726       // OpenMM appears to offer some support for storing dE/dX, dE/dY and d2E/dXdY
727       // instead of computing these values on the fly.
728       Element grid = doc.createElement("Grid");
729       grid.setAttribute("angle1", format("%f", tx[nxi]));
730       grid.setAttribute("angle2", format("%f", ty[nyi]));
731       grid.setAttribute("f", format("%f", energy[i] * KCAL_TO_KJ));
732       ttGrid.appendChild(grid);
733     }
734 
735     tortors.setAttribute("grid", String.valueOf(label));
736     ttGrid.setAttribute("grid", String.valueOf(label));
737     torTorNode.appendChild(tortors);
738     torTorNode.appendChild(ttGrid);
739   }
740 
741   /**
742    * Computes the coefficients for an aperiodic interpolating cubic spline.
743    *
744    * @param n
745    * @param x0
746    * @param y0
747    * @param y21
748    * @param y2n
749    * @param s1
750    * @param s2
751    * @param h
752    * @param g
753    * @param dy
754    * @param dla
755    * @param dmu
756    */
757   private void nspline(int n, double[] x0, double[] y0, double y21, double y2n, double[] s1,
758                        double[] s2, double[] h, double[] g, double[] dy, double[] dla, double[] dmu) {
759 
760     // Calculate the intervals.
761     for (int i = 0; i < n; i++) {
762       h[i] = x0[i + 1] - x0[i];
763       dy[i] = (y0[i + 1] - y0[i]) / h[i];
764     }
765 
766     // Get the coeffcients.
767     for (int i = 1; i < n; i++) {
768       dla[i] = h[i] / (h[i] + h[i - 1]);
769       dmu[i] = 1.0 - dla[i];
770       g[i] = 3.0 * (dla[i] * dy[i - 1] + dmu[i] * dy[i]);
771     }
772 
773     // Set the initial value via natural boundary condition.
774     dla[n] = 1.0;
775     dla[0] = 0.0;
776     dmu[n] = 0.0;
777     dmu[0] = 1.0;
778     g[0] = 3.0 * dy[0] - 0.5 * h[0] * y21;
779     g[n] = 3.0 * dy[n - 1] + 0.5 * h[n - 1] * y2n;
780 
781     // Solve the triagonal system of linear equations.
782     dmu[0] = 0.5 * dmu[0];
783     g[0] = 0.5 * g[0];
784     for (int i = 1; i <= n; i++) {
785       double t = 2.0 - dmu[i - 1] * dla[i];
786       dmu[i] = dmu[i] / t;
787       g[i] = (g[i] - g[i - 1] * dla[i]) / t;
788     }
789 
790     for (int i = n - 1; i >= 0; i--) {
791       g[i] = g[i] - dmu[i] * g[i + 1];
792     }
793 
794     // Get the first derivative at each grid point.
795     arraycopy(g, 0, s1, 0, n + 1);
796 
797     // Get the second derivative at each grid point.
798     s2[0] = y21;
799     s2[n] = y2n;
800     for (int i = 1; i < n; i++) {
801       s2[i] =
802           6.0 * (y0[i + 1] - y0[i]) / (h[i] * h[i]) - 4.0 * s1[i] / h[i] - 2.0 * s1[i + 1] / h[i];
803     }
804   }
805 
806   /**
807    * Computes the coefficients for a periodic interpolating cubic spline.
808    *
809    * @param n
810    * @param xn
811    * @param fn
812    * @param b
813    * @param d
814    * @param h
815    * @param du
816    * @param dm
817    * @param rc
818    * @param rs
819    */
820   private void cspline(int n, double[] xn, double[] fn, double[] b, double[] c, double[] d,
821                        double[] h, double[] du, double[] dm, double[] rc, double[] rs) {
822     double eps = 0.000001;
823     if (abs(fn[n] - fn[0]) > eps) {
824       logger.severe("TORTOR values are not periodic.");
825     }
826 
827     // Make the spline exactly periodic making the ends equal to their mean.
828     double mean = 0.5 * (fn[0] + fn[n]);
829     fn[0] = mean;
830     fn[n] = mean;
831 
832     // Set up auxiliary variables and matrix elements on first call.
833     for (int i = 0; i < n; i++) {
834       h[i] = xn[i + 1] - xn[i];
835     }
836     h[n] = h[0];
837 
838     if (n - 1 >= 0) {
839       arraycopy(h, 1, du, 1, n - 1);
840     }
841 
842     du[n] = h[0];
843     for (int i = 1; i <= n; i++) {
844       dm[i] = 2.0 * (h[i - 1] + h[i]);
845     }
846 
847     // Compute the RHS.
848     double temp1 = (fn[1] - fn[0]) / h[0];
849     for (int i = 1; i < n; i++) {
850       double temp2 = (fn[i + 1] - fn[i]) / h[i];
851       rs[i] = 3.0 * (temp2 - temp1);
852       temp1 = temp2;
853     }
854     rs[n] = 3.0 * ((fn[1] - fn[0]) / h[0] - temp1);
855 
856     // Solve the linear system with factorization.
857     if (cytsy(n, dm, du, rc, rs, c) != 1) {
858       return;
859     }
860 
861     // Compute remaining spline coefficients.
862     c[0] = c[n];
863     for (int i = 0; i < n; i++) {
864       b[i] = (fn[i + 1] - fn[i]) / h[i] - h[i] / 3.0 * (c[i + 1] + 2.0 * c[i]);
865       d[i] = (c[i + 1] - c[i]) / (3.0 * h[i]);
866     }
867     b[n] = (fn[1] - fn[n]) / h[n] - h[n] / 3.0 * (c[1] + 2.0 * c[n]);
868   }
869 }