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.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
67
68
69
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
85
86 public static final double DEFAULT_TORTOR_UNIT = 1.0;
87
88
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
99
100 public final int[] atomClasses;
101
102
103
104 public final double[] energy;
105
106
107
108 public final int nx;
109
110
111
112 public final int ny;
113
114
115
116 public final double[] tx;
117
118
119
120 public final double[] ty;
121
122
123
124 public final double[] dx;
125
126
127
128 public final double[] dy;
129
130
131
132 public final double[] dxy;
133
134
135
136 private final int[] gridPoints;
137
138
139
140
141
142
143
144
145
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
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
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
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
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
283
284
285
286
287
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
300
301
302
303
304
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
331 torsion1[i] = parseDouble(tokens[0]);
332 torsion2[i] = parseDouble(tokens[1]);
333 energy[i] = parseDouble(tokens[2]);
334
335 torsion1[i + 1] = parseDouble(tokens[3]);
336 torsion2[i + 1] = parseDouble(tokens[4]);
337 energy[i + 1] = parseDouble(tokens[5]);
338
339 i += 1;
340 } else if (tokens.length == 9) {
341
342 torsion1[i] = parseDouble(tokens[0]);
343 torsion2[i] = parseDouble(tokens[1]);
344 energy[i] = parseDouble(tokens[2]);
345
346 torsion1[i + 1] = parseDouble(tokens[3]);
347 torsion2[i + 1] = parseDouble(tokens[4]);
348 energy[i + 1] = parseDouble(tokens[5]);
349
350 torsion1[i + 2] = parseDouble(tokens[6]);
351 torsion2[i + 2] = parseDouble(tokens[7]);
352 energy[i + 2] = parseDouble(tokens[8]);
353
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
371
372
373
374
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
414
415
416
417
418 public static String reverseKey(int[] c) {
419 return c[4] + " " + c[3] + " " + c[2] + " " + c[1] + " " + c[0];
420 }
421
422
423
424
425
426
427
428 public static String sortKey(int[] c) {
429 return c[0] + " " + c[1] + " " + c[2] + " " + c[3] + " " + c[4];
430 }
431
432
433
434
435
436
437
438
439
440
441
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
449 if (cytsyp(n, dm, du, cr) == 1) {
450
451 cytsys(n, dm, du, cr, rs, c);
452 return 1;
453 }
454 return -1;
455 }
456
457
458
459
460
461
462
463
464
465
466
467 private static int cytsyp(int n, double[] dm, double[] du, double[] cr) {
468 double eps = 0.00000001;
469
470 if (n < 3) {
471 return -2;
472 }
473
474
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
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
533
534
535
536
537
538
539
540
541 private static void cytsys(int n, double[] dm, double[] du, double[] cr, double[] rs, double[] c) {
542
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
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
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
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
599
600 @Override
601 public int hashCode() {
602 return Arrays.hashCode(atomClasses);
603 }
604
605
606
607
608
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
619
620
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
647
648
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
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
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
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
683
684
685
686
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
704
705
706
707
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
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
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
729
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
745
746
747
748
749
750
751
752
753
754
755
756
757
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
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
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
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
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
797 arraycopy(g, 0, s1, 0, n + 1);
798
799
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
810
811
812
813
814
815
816
817
818
819
820
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
830 double mean = 0.5 * (fn[0] + fn[n]);
831 fn[0] = mean;
832 fn[n] = mean;
833
834
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
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
859 if (cytsy(n, dm, du, rc, rs, c) != 1) {
860 return;
861 }
862
863
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 }