1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
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
58
59
60
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
73 public SturmMethod() {
74 RELERROR = 1.0e-15;
75 MAXIT = 100;
76 MAX_ITER_SECANT = 20;
77 }
78
79
80
81
82
83
84
85
86
87
88
89
90
91
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
138
139
140
141 public double[][] getr_o() {
142 return xyz_o;
143 }
144
145
146
147
148
149
150
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
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
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
226 nroots = atmin[0] - atmax[0];
227 sturmBisection(np, sseq, min, max, atmin[0], atmax[0], this.roots);
228
229
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
248
249
250
251
252
253
254
255
256
257
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
373
374
375
376
377
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
394
395
396
397
398
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
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
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
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
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
439
440
441
442
443
444
445
446 private int numRoots(int np, Polynomial[] sseq, int[] atneg, int[] atpos) {
447
448 int atposinf = 0;
449 int atneginf = 0;
450
451
452 double lf = sseq[0].coefficients[sseq[0].order];
453
454 for (int i = 1; i <= np; i++)
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
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
491
492
493
494
495
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
515
516
517
518
519
520
521
522
523
524
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
536 if (modifiedRegulaFalsi(sseq[0].order, sseq[0].coefficients, min, max, roots)) {
537 fillArray(roots, this.roots);
538 return;
539 }
540
541
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
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
582
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
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
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
637
638
639
640
641
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
657
658
659
660
661
662
663
664
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
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
718 private static class Polynomial {
719
720 static final int MAX_ORDER = 16;
721
722 public final double[] coefficients;
723
724 public int order;
725
726
727 private Polynomial() {
728 coefficients = new double[MAX_ORDER + 1];
729 }
730 }
731 }