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