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