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