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-2024.
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.bonded;
39  
40  import static java.lang.String.format;
41  import static java.lang.System.arraycopy;
42  import static org.apache.commons.math3.util.FastMath.abs;
43  import static org.apache.commons.math3.util.FastMath.exp;
44  import static org.apache.commons.math3.util.FastMath.max;
45  
46  import ffx.potential.MolecularAssembly;
47  import ffx.potential.parsers.PDBFilter;
48  import java.io.File;
49  import java.util.ArrayList;
50  import java.util.Arrays;
51  import java.util.List;
52  import java.util.logging.Level;
53  import java.util.logging.Logger;
54  import org.apache.commons.io.FilenameUtils;
55  
56  /**
57   * SturmMethod class.
58   *
59   * @author Mallory R. Tollefson
60   * @since 1.0
61   */
62  public class SturmMethod {
63  
64    private static final Logger logger = Logger.getLogger(SturmMethod.class.getName());
65    final int MAXPOW = 32;
66    private final double RELERROR;
67    private final int MAXIT;
68    private final int MAX_ITER_SECANT;
69    private final double[][] xyz_o = new double[5][3];
70    private double[] roots;
71  
72    /** SturmMethod constructor with termination criteria for polynomial solver. */
73    public SturmMethod() {
74      RELERROR = 1.0e-15;
75      MAXIT = 100;
76      MAX_ITER_SECANT = 20;
77    }
78  
79    /**
80     * Calculates the modulus of u(x)/v(x) leaving it in r, it returns 0 if r(x) is constant. This
81     * function assumes the leading coefficient of v is 1 or -1.
82     *
83     * <p>Modp was originally a static int and returned r.ord as an int. The buildSturm function
84     * requires that a boolean is returned from the modp function, so modp is set to return a boolean.
85     * This new boolean return could be problematic elsewhere in the code if modp is used to return a
86     * value. This should be checked.
87     *
88     * @param u
89     * @param v
90     * @param r
91     * @return true if r.ord > 0.
92     */
93    private static boolean modp(Polynomial u, Polynomial v, Polynomial r) {
94      double[] nr = r.coefficients;
95      int end = u.order;
96      double[] uc = u.coefficients;
97  
98      if (end + 1 >= 0) arraycopy(uc, 0, nr, 0, end + 1);
99  
100     if (v.coefficients[v.order] < 0.0) {
101       for (int k = u.order - v.order - 1; k >= 0; k -= 2) {
102         r.coefficients[k] = -r.coefficients[k];
103       }
104 
105       for (int k = u.order - v.order; k >= 0; k--) {
106         for (int j = v.order + k - 1; j >= k; j--) {
107           r.coefficients[j] =
108               -r.coefficients[j] - r.coefficients[v.order + k] * v.coefficients[j - k];
109         }
110       }
111     } else {
112       for (int k = u.order - v.order; k >= 0; k--) {
113         for (int j = v.order + k - 1; j >= k; j--) {
114           r.coefficients[j] -= r.coefficients[v.order + k] * v.coefficients[j - k];
115         }
116       }
117     }
118 
119     int k = v.order - 1;
120     while (k >= 0 && abs(r.coefficients[k]) < 1.0e-18) {
121       r.coefficients[k] = 0.0;
122       k--;
123     }
124 
125     r.order = max(k, 0);
126 
127     if (r.order > 0) {
128       return true;
129     } else if (r.order == 0) {
130       return false;
131     } else {
132       return false;
133     }
134   }
135 
136   /**
137    * Used only in JUnit testing.
138    *
139    * @return oxygen coordinates.
140    */
141   public double[][] getr_o() {
142     return xyz_o;
143   }
144 
145   /**
146    * Solve using the Sturm method.
147    *
148    * @param order the order of the polynomial.
149    * @param poly_coeffs an array of {@link double} objects.
150    * @param roots an array of {@link double} objects.
151    */
152   public int solveSturm(int order, double[] poly_coeffs, double[] roots) {
153     Polynomial[] sseq = new Polynomial[Polynomial.MAX_ORDER * 2];
154     double min, max;
155     int nroots, nchanges, np;
156     int[] atmin = new int[1];
157     int[] atmax = new int[1];
158     this.roots = roots;
159 
160     for (int i = 0; i < Polynomial.MAX_ORDER * 2; i++) {
161       sseq[i] = new Polynomial();
162     }
163 
164     if (order + 1 >= 0) arraycopy(poly_coeffs, 0, sseq[0].coefficients, 0, order + 1);
165 
166     if (logger.isLoggable(Level.FINE)) {
167       StringBuilder string = new StringBuilder();
168       for (int i = order; i >= 0; i--) {
169         string.append(" Coefficients in Sturm solver\n");
170         string.append(format("%d %f\n", i, sseq[0].coefficients[i]));
171       }
172       logger.fine(string.toString());
173     }
174 
175     np = buildSturm(order, sseq);
176 
177     if (logger.isLoggable(Level.FINE)) {
178       StringBuilder string1 = new StringBuilder();
179       string1.append(" Sturm sequence for:\n");
180       for (int i = order; i >= 0; i--) {
181         string1.append(format("%f ", sseq[0].coefficients[i]));
182         string1.append("\n");
183       }
184       for (int i = 0; i <= np; i++) {
185         for (int j = sseq[i].order; j >= 0; j--) {
186           string1.append(format("%f ", sseq[i].coefficients[j]));
187           string1.append("\n");
188         }
189       }
190       logger.fine(string1.toString());
191     }
192 
193     // get the number of real roots
194     nroots = numRoots(np, sseq, atmin, atmax);
195     if (logger.isLoggable(Level.FINE)) {
196       logger.fine(format(" Number of real roots: %d\n", nroots));
197     }
198 
199     // calculate the bracket that the roots live in
200     min = -1.0;
201     nchanges = numChanges(np, sseq, min);
202 
203     for (int i = 0; nchanges != atmin[0] && i != MAXPOW; i++) {
204       min *= 10.0;
205       nchanges = numChanges(np, sseq, min);
206     }
207 
208     if (nchanges != atmin[0]) {
209       logger.fine(" Solve: unable to bracket all negative roots\n");
210       atmin[0] = nchanges;
211     }
212 
213     max = 1.0;
214     nchanges = numChanges(np, sseq, max);
215 
216     for (int i = 0; nchanges != atmax[0] && i != MAXPOW; i++) {
217       max *= 10.0;
218       nchanges = numChanges(np, sseq, max);
219     }
220     if (nchanges != atmax[0]) {
221       logger.fine(" Solve: unable to bracket all positive roots\n");
222       atmax[0] = nchanges;
223     }
224 
225     // Perform the bisection
226     nroots = atmin[0] - atmax[0];
227     sturmBisection(np, sseq, min, max, atmin[0], atmax[0], this.roots);
228 
229     // write out the roots
230     if (logger.isLoggable(Level.FINE)) {
231       if (nroots == 1) {
232         logger.fine(format("\n One distinct real root at x = %f\n", this.roots[0]));
233       } else {
234         StringBuilder string2 = new StringBuilder();
235         string2.append(format("\n %d distinct real roots for x: \n", nroots));
236         for (int i = 0; i != nroots; i++) {
237           string2.append(format("%f\n", this.roots[i]));
238         }
239         logger.fine(string2.toString());
240       }
241     }
242 
243     return nroots;
244   }
245 
246   /**
247    * Write out loop coordinates and determine oxygen placement.
248    *
249    * @param r_n an array of {@link double} objects.
250    * @param r_a an array of {@link double} objects.
251    * @param r_c an array of {@link double} objects.
252    * @param stt_res a int.
253    * @param end_res a int.
254    * @param molAss a {@link ffx.potential.MolecularAssembly} object.
255    * @param counter a int.
256    * @param writeFile a boolean.
257    * @return the File.
258    */
259   public File writePDBBackbone(
260       double[][] r_n,
261       double[][] r_a,
262       double[][] r_c,
263       int stt_res,
264       int end_res,
265       MolecularAssembly molAss,
266       int counter,
267       boolean writeFile) {
268 
269     Polymer[] newChain = molAss.getChains();
270     List<Atom> backBoneAtoms;
271     double[] xyz_n = new double[3];
272     double[] xyz_a = new double[3];
273     double[] xyz_c = new double[3];
274 
275     xyz_o[0][0] = 0.0;
276     xyz_o[0][1] = 0.0;
277     xyz_o[0][2] = 0.0;
278     xyz_o[4][0] = 0.0;
279     xyz_o[4][1] = 0.0;
280     xyz_o[4][2] = 0.0;
281 
282     List<Atom> OAtoms = new ArrayList<>();
283 
284     for (int i = stt_res + 1; i < end_res; i++) {
285       Residue newResidue = newChain[0].getResidue(i);
286       backBoneAtoms = newResidue.getBackboneAtoms();
287       for (Atom backBoneAtom : backBoneAtoms) {
288         switch (backBoneAtom.getAtomType().name) {
289           case "C":
290             xyz_c[0] = r_c[i - stt_res][0];
291             xyz_c[1] = r_c[i - stt_res][1];
292             xyz_c[2] = r_c[i - stt_res][2];
293             backBoneAtom.moveTo(xyz_c);
294             break;
295           case "N":
296             xyz_n[0] = r_n[i - stt_res][0];
297             xyz_n[1] = r_n[i - stt_res][1];
298             xyz_n[2] = r_n[i - stt_res][2];
299             backBoneAtom.moveTo(xyz_n);
300             break;
301           case "CA":
302             xyz_a[0] = r_a[i - stt_res][0];
303             xyz_a[1] = r_a[i - stt_res][1];
304             xyz_a[2] = r_a[i - stt_res][2];
305             backBoneAtom.moveTo(xyz_a);
306             break;
307           case "O":
308             OAtoms.add(backBoneAtom);
309             break;
310           default:
311             newResidue.deleteAtom(backBoneAtom);
312             break;
313         }
314       }
315 
316       List<Atom> sideChainAtoms = newResidue.getSideChainAtoms();
317       for (Atom sideChainAtom : sideChainAtoms) {
318         newResidue.deleteAtom(sideChainAtom);
319       }
320     }
321 
322     int oCount = 0;
323     for (int i = stt_res + 1; i < end_res; i++) {
324       Residue newResidue = newChain[0].getResidue(i);
325       backBoneAtoms = newResidue.getBackboneAtoms();
326       Atom CA = new Atom("CA");
327       Atom N = new Atom("N");
328       Atom C = new Atom("C");
329       Atom O = OAtoms.get(oCount);
330 
331       for (Atom backBoneAtom : backBoneAtoms) {
332         switch (backBoneAtom.getAtomType().name) {
333           case "C":
334             C = backBoneAtom;
335             break;
336           case "N":
337             N = backBoneAtom;
338             break;
339           case "CA":
340             CA = backBoneAtom;
341             break;
342           default:
343             break;
344         }
345       }
346       BondedUtils.intxyz(O, C, 1.2255, CA, 122.4, N, 180, 0);
347       xyz_o[i - stt_res][0] = O.getX();
348       xyz_o[i - stt_res][1] = O.getY();
349       xyz_o[i - stt_res][2] = O.getZ();
350 
351       oCount++;
352     }
353 
354     File file = molAss.getFile();
355 
356     String filename = FilenameUtils.removeExtension(file.getAbsolutePath());
357     if (!filename.contains("_loop")) {
358       filename = filename + "_loop";
359     }
360 
361     File modifiedFile = new File(filename + ".pdb_" + counter);
362     PDBFilter modFilter = new PDBFilter(modifiedFile, molAss, null, null);
363 
364     if (writeFile) {
365       modFilter.writeFile(modifiedFile, true);
366     }
367 
368     return (modifiedFile);
369   }
370 
371   /**
372    * The hyperTan method.
373    *
374    * @param a
375    * @param x
376    * @return 1.0 if (ax > 100), -1.0 if (ax < -100) and (exp(ax)-exp(-ax))/(exp(ax)+exp(-ax))
377    *     otherwise.
378    */
379   private double hyperTan(double a, double x) {
380     double ax = a * x;
381     if (ax > 100.0) {
382       return 1.0;
383     } else if (ax < -100.0) {
384       return -1.0;
385     } else {
386       double exp_x1 = exp(ax);
387       double exp_x2 = exp(-ax);
388       return (exp_x1 - exp_x2) / (exp_x1 + exp_x2);
389     }
390   }
391 
392   /**
393    * Build up a sturm sequence for a polynomial, and return the number of polynomials in the
394    * sequence.
395    *
396    * @param ord
397    * @param sseq
398    * @return the number of polynomials in the sequence.
399    */
400   private int buildSturm(int ord, Polynomial[] sseq) {
401     double f;
402     double[] fp;
403     double[] fc;
404 
405     sseq[0].order = ord;
406     sseq[1].order = ord - 1;
407 
408     // Calculate the derivative and normalise the leading coefficient
409     f = abs(sseq[0].coefficients[ord] * ord);
410     fp = sseq[1].coefficients;
411     fc = Arrays.copyOfRange(sseq[0].coefficients, 1, sseq[0].coefficients.length);
412 
413     int j = 0;
414     for (int i = 1; i <= ord; i++) {
415       fp[j] = fc[j] * i / f;
416       j++;
417     }
418 
419     // Construct the rest of the Sturm sequence. (Double check this... )
420     int i;
421     for (i = 0;
422         i < sseq[0].coefficients.length - 2 && modp(sseq[i], sseq[i + 1], sseq[i + 2]);
423         i++) {
424       // reverse the sign and normalise
425       f = -Math.abs(sseq[i + 2].coefficients[sseq[i + 2].order]);
426       for (j = sseq[i + 2].order; j >= 0; j--) {
427         sseq[i + 2].coefficients[j] /= f;
428       }
429     }
430 
431     // Reverse the sign.
432     sseq[i + 2].coefficients[0] = -sseq[i + 2].coefficients[0];
433 
434     return (sseq[0].order - sseq[i + 2].order);
435   }
436 
437   /**
438    * Return the number of distinct real roots of the polynomial described in sseq.
439    *
440    * @param np
441    * @param sseq
442    * @param atneg
443    * @param atpos
444    * @return the number of distinct real roots of sseq.
445    */
446   private int numRoots(int np, Polynomial[] sseq, int[] atneg, int[] atpos) {
447 
448     int atposinf = 0;
449     int atneginf = 0;
450 
451     // Changes at positive infinity.
452     double lf = sseq[0].coefficients[sseq[0].order];
453 
454     for (int i = 1; i <= np; i++) // for (s = sseq + 1; s <= sseq + np; s++)
455     {
456       double f = sseq[i].coefficients[sseq[i].order];
457       if (lf == 0.0 || lf * f < 0) {
458         atposinf++;
459       }
460       lf = f;
461     }
462 
463     // Changes at negative infinity.
464     if ((sseq[0].order & 1) != 0) {
465       lf = -sseq[0].coefficients[sseq[0].order];
466     } else {
467       lf = sseq[0].coefficients[sseq[0].order];
468     }
469 
470     for (int i = 1; i <= np; i++) {
471       double f;
472       if ((sseq[i].order & 1) != 0) {
473         f = -sseq[i].coefficients[sseq[i].order];
474       } else {
475         f = sseq[i].coefficients[sseq[i].order];
476       }
477       if (lf == 0.0 || lf * f < 0) {
478         atneginf++;
479       }
480       lf = f;
481     }
482 
483     atneg[0] = atneginf;
484     atpos[0] = atposinf;
485 
486     return (atneginf - atposinf);
487   }
488 
489   /**
490    * Return the number of sign changes in the Sturm sequence in sseq at the value a.
491    *
492    * @param np
493    * @param sseq
494    * @param a
495    * @return the number of sign changes.
496    */
497   private int numChanges(int np, Polynomial[] sseq, double a) {
498 
499     int changes = 0;
500     double lf = evalPoly(sseq[0].order, sseq[0].coefficients, a);
501     for (int i = 1; i <= np; i++) {
502       double f = evalPoly(sseq[i].order, sseq[i].coefficients, a);
503       if (lf == 0.0 || lf * f < 0) {
504         changes++;
505       }
506 
507       lf = f;
508     }
509 
510     return changes;
511   }
512 
513   /**
514    * Uses a bisection based on the Sturm sequence for the polynomial described in sseq to isolate
515    * intervals in which roots occur, the roots are returned in the roots array in order of
516    * magnitude.
517    *
518    * @param np
519    * @param sseq
520    * @param min
521    * @param max
522    * @param atmin
523    * @param atmax
524    * @param roots
525    */
526   private void sturmBisection(
527       int np, Polynomial[] sseq, double min, double max, int atmin, int atmax, double[] roots) {
528     double mid = 0;
529     int n1 = 0;
530     int n2 = 0;
531     int its, atmid, nroot;
532 
533     nroot = atmin - atmax;
534     if (nroot == 1) {
535       // First try a less expensive technique
536       if (modifiedRegulaFalsi(sseq[0].order, sseq[0].coefficients, min, max, roots)) {
537         fillArray(roots, this.roots);
538         return;
539       }
540 
541       // When code reaches this point, the root must be evaluated using the Sturm sequence.
542       for (its = 0; its < MAXIT; its++) {
543         mid = (min + max) / 2;
544         atmid = numChanges(np, sseq, mid);
545         if (abs(mid) > RELERROR) {
546           if (abs((max - min) / mid) < RELERROR) {
547             roots[0] = mid;
548             fillArray(roots, this.roots);
549             return;
550           }
551         } else if (abs(max - min) < RELERROR) {
552           roots[0] = mid;
553           fillArray(roots, this.roots);
554           return;
555         }
556         if ((atmin - atmid) == 0) {
557           min = mid;
558         } else {
559           max = mid;
560         }
561       }
562       if (its == MAXIT) {
563         logger.info(
564             format(
565                 " sbisect: overflow min %f max %f diff %f nroot %d n1 %d n2 %d\n",
566                 min, max, max - min, nroot, n1, n2));
567         roots[0] = mid;
568       }
569       fillArray(roots, this.roots);
570       return;
571     }
572 
573     // More than one root in the interval, must bisect.
574     for (its = 0; its < MAXIT; its++) {
575       mid = (min + max) / 2;
576       atmid = numChanges(np, sseq, mid);
577       n1 = atmin - atmid;
578       n2 = atmid - atmax;
579       if (n1 != 0 && n2 != 0) {
580         sturmBisection(np, sseq, min, mid, atmin, atmid, roots);
581         // Keep a reference to this array so that its values can be copied into the
582         // roots array later.
583         final double[] rootsSection = Arrays.copyOfRange(roots, n1, roots.length);
584         sturmBisection(np, sseq, mid, max, atmid, atmax, rootsSection);
585         fillArray(rootsSection, roots);
586         break;
587       }
588       if (n1 == 0) {
589         min = mid;
590       } else {
591         max = mid;
592       }
593     }
594 
595     if (its == MAXIT) {
596       for (n1 = atmax; n1 < atmin; n1++) {
597         roots[n1 - atmax] = mid;
598       }
599     }
600   }
601 
602   /**
603    * This method fills the target array with values from the values array. The context of this
604    * method makes the most sense when the values array is smaller than the target array, resulting
605    * in the data from the values array simply being added to the target array at the first point
606    * that would entail the final value being copied into the target array's final index. See below
607    * for an example input and output.
608    *
609    * <pre>
610    * Inputs:
611    * valuesArray: [ 10, 20, 30,  0 ]
612    * targetArray: [  1,  2,  3,  4,  0,  0,  0,  0,  0,  0];
613    *
614    * Output (via array reference with later use of target array):
615    * [  1,  2,  3,  4,  0,  0, 10, 20, 30, 0 ]
616    * </pre>
617    *
618    * @param valuesArray The array that will hold the values to copy into the target array.
619    * @param targetArray The array that will receive the values in the values array.
620    */
621   private void fillArray(final double[] valuesArray, final double[] targetArray) {
622     final int tarLen = targetArray.length;
623     final int valLen = valuesArray.length;
624     if (tarLen >= valLen) {
625       for (int tarIndex = tarLen - valLen, valIndex = 0;
626           tarIndex < tarLen;
627           tarIndex++, valIndex++) {
628         targetArray[tarIndex] = valuesArray[valIndex];
629       }
630     } else {
631       logger.info(" sbisect: length of values array is too long for the target array; not filling");
632     }
633   }
634 
635   /**
636    * Evaluate polynomial defined in coef returning its value.
637    *
638    * @param ord
639    * @param coef
640    * @param x
641    * @return value of the polynomial at x.
642    */
643   private double evalPoly(int ord, double[] coef, double x) {
644 
645     double[] fp = coef;
646     double f = fp[ord];
647 
648     int i = ord;
649     for (i--; i >= 0; i--) {
650       f = x * f + fp[i];
651     }
652     return f;
653   }
654 
655   /**
656    * Uses the modified regula-falsi method to evaluate the root in interval [a,b] of the polynomial
657    * described in coef. The root is returned in val. The routine returns zero if it can't converge.
658    *
659    * @param ord
660    * @param coef
661    * @param a
662    * @param b
663    * @param val
664    * @return 0 if the method does not converge.
665    */
666   private boolean modifiedRegulaFalsi(int ord, double[] coef, double a, double b, double[] val) {
667     double fa = coef[ord];
668     double fb = fa;
669 
670     for (int i = ord - 1; i >= 0; i--) {
671       fa = a * fa + coef[i];
672       fb = b * fb + coef[i];
673     }
674     if (fa * fb > 0.0) {
675       return false;
676     }
677     double lfx = fa;
678     for (int its = 0; its < MAX_ITER_SECANT; its++) {
679       double x = (fb * a - fa * b) / (fb - fa);
680       // constrain that x stays in the bounds
681       if (x < a || x > b) {
682         x = 0.5 * (a + b);
683       }
684       double fx = coef[ord];
685       for (int i = ord - 1; i >= 0; i--) {
686         fx = x * fx + coef[i];
687       }
688       if (abs(x) > RELERROR) {
689         if (abs(fx / x) < RELERROR) {
690           val[0] = x;
691           return true;
692         }
693       } else if (abs(fx) < RELERROR) {
694         val[0] = x;
695         return true;
696       }
697       if ((fa * fx) < 0) {
698         b = x;
699         fb = fx;
700         if ((lfx * fx) > 0) {
701           fa /= 2;
702         }
703       } else {
704         a = x;
705         fa = fx;
706         if ((lfx * fx) > 0) {
707           fb /= 2;
708         }
709       }
710 
711       lfx = fx;
712     }
713 
714     return false;
715   }
716 
717   /** The Polynomial class describes a single polynomial of up to 16th degree. */
718   private static class Polynomial {
719 
720     static final int MAX_ORDER = 16;
721     /** The coefficients of the polynomial. */
722     public final double[] coefficients;
723     /** The order of the polynomial. */
724     public int order;
725 
726     /** Private constructor. */
727     private Polynomial() {
728       coefficients = new double[MAX_ORDER + 1];
729     }
730   }
731 }