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.algorithms.optimize.manybody;
39
40 import ffx.algorithms.AlgorithmListener;
41 import ffx.algorithms.optimize.RotamerOptimization;
42 import ffx.numerics.Potential;
43 import ffx.potential.MolecularAssembly;
44 import ffx.potential.bonded.Atom;
45 import ffx.potential.bonded.Residue;
46 import ffx.potential.bonded.Rotamer;
47 import ffx.potential.openmm.OpenMMEnergy;
48 import ffx.potential.utils.EnergyException;
49 import org.apache.commons.configuration2.CompositeConfiguration;
50
51 import java.io.File;
52 import java.io.IOException;
53 import java.nio.charset.StandardCharsets;
54 import java.nio.file.Files;
55 import java.nio.file.Path;
56 import java.nio.file.Paths;
57 import java.util.ArrayList;
58 import java.util.HashMap;
59 import java.util.List;
60 import java.util.Map;
61 import java.util.Set;
62 import java.util.logging.Level;
63 import java.util.logging.Logger;
64
65 import static ffx.potential.bonded.RotamerLibrary.applyRotamer;
66 import static java.lang.String.format;
67
68 public class EnergyExpansion {
69
70 private static final Logger logger = Logger.getLogger(EnergyExpansion.class.getName());
71
72
73
74 private static final double DEFAULT_OMM_RECALCULATE_THRESHOLD = -200;
75
76
77
78 private static final double DEFAULT_SINGULARITY_THRESHOLD = -1000;
79
80
81
82 private final Map<Integer, Integer[]> selfEnergyMap = new HashMap<>();
83
84
85
86 private final Map<Integer, Integer[]> twoBodyEnergyMap = new HashMap<>();
87
88
89
90 private final Map<Integer, Integer[]> threeBodyEnergyMap = new HashMap<>();
91
92
93
94 private final Map<Integer, Integer[]> fourBodyEnergyMap = new HashMap<>();
95
96
97
98 private final boolean master;
99
100
101
102
103 private final double ommRecalculateThreshold;
104
105
106
107
108
109
110 private final double singularityThreshold;
111
112
113
114 private final boolean potentialIsOpenMM;
115 private final RotamerOptimization rO;
116 private final DistanceMatrix dM;
117 private final EliminatedRotamers eR;
118
119
120
121 private final MolecularAssembly molecularAssembly;
122
123
124
125 private final AlgorithmListener algorithmListener;
126
127
128
129
130 private final List<Residue> allResiduesList;
131
132
133
134 private final int[][] resNeighbors;
135
136
137
138 private final boolean threeBodyTerm;
139
140
141
142 private final boolean decomposeOriginal;
143
144
145
146 private final boolean usingBoxOptimization;
147
148
149
150 private final boolean verbose;
151
152
153
154 private final boolean pruneClashes;
155
156
157
158 private final boolean prunePairClashes;
159
160
161
162 private final int max4BodyCount;
163
164
165
166 private double backboneEnergy;
167
168
169
170 private double[][] selfEnergy;
171
172
173
174
175 private double[][][][] twoBodyEnergy;
176
177
178
179
180 private double[][][][][][] threeBodyEnergy;
181
182 public EnergyExpansion(RotamerOptimization rO, DistanceMatrix dM, EliminatedRotamers eR,
183 MolecularAssembly molecularAssembly, Potential potential, AlgorithmListener algorithmListener,
184 List<Residue> allResiduesList, int[][] resNeighbors, boolean threeBodyTerm,
185 boolean decomposeOriginal, boolean usingBoxOptimization, boolean verbose,
186 boolean pruneClashes, boolean prunePairClashes, boolean master) {
187 this.rO = rO;
188 this.dM = dM;
189 this.eR = eR;
190 this.molecularAssembly = molecularAssembly;
191 this.algorithmListener = algorithmListener;
192 this.allResiduesList = allResiduesList;
193 this.resNeighbors = resNeighbors;
194 this.threeBodyTerm = threeBodyTerm;
195 this.decomposeOriginal = decomposeOriginal;
196 this.usingBoxOptimization = usingBoxOptimization;
197 this.verbose = verbose;
198 this.pruneClashes = pruneClashes;
199 this.prunePairClashes = prunePairClashes;
200 this.master = master;
201
202 CompositeConfiguration properties = molecularAssembly.getProperties();
203 max4BodyCount = properties.getInt("ro-max4BodyCount", Integer.MAX_VALUE);
204 if (max4BodyCount != Integer.MAX_VALUE) {
205 logger.info(format(" Max 4Body Count: %d", max4BodyCount));
206 }
207 singularityThreshold = properties.getDouble("ro-singularityThreshold", DEFAULT_SINGULARITY_THRESHOLD);
208 potentialIsOpenMM = potential instanceof OpenMMEnergy;
209 if (potentialIsOpenMM) {
210 ommRecalculateThreshold = properties.getDouble("ro-ommRecalculateThreshold", DEFAULT_OMM_RECALCULATE_THRESHOLD);
211 } else {
212 ommRecalculateThreshold = -1E200;
213 }
214 }
215
216
217
218
219
220
221 private static void turnOffAtoms(Residue residue) {
222 switch (residue.getResidueType()) {
223 case NA:
224 case AA:
225 List<Atom> atomList = residue.getVariableAtoms();
226 for (Atom atom : atomList) {
227 atom.setUse(false);
228 }
229 break;
230 default:
231 atomList = residue.getAtomList();
232 for (Atom atom : atomList) {
233 atom.setUse(false);
234 }
235 break;
236 }
237 }
238
239
240
241
242
243
244 private static void turnOnAtoms(Residue residue) {
245 switch (residue.getResidueType()) {
246 case NA:
247 case AA:
248 List<Atom> atomList = residue.getVariableAtoms();
249 for (Atom atom : atomList) {
250 atom.setUse(true);
251 }
252 break;
253 default:
254 atomList = residue.getAtomList();
255 for (Atom atom : atomList) {
256 atom.setUse(true);
257 }
258 break;
259 }
260 }
261
262 public HashMap<String, Integer> allocate2BodyJobMap(Residue[] residues, int nResidues,
263 boolean reverseMap) {
264 twoBodyEnergyMap.clear();
265
266 HashMap<String, Integer> reverseJobMapPairs = new HashMap<>();
267 int pairJobIndex = 0;
268 twoBodyEnergy = new double[nResidues][][][];
269 for (int i = 0; i < nResidues; i++) {
270 Residue resi = residues[i];
271 int indexI = allResiduesList.indexOf(resi);
272 Rotamer[] roti = resi.getRotamers();
273 int[] nI = resNeighbors[i];
274 int lenNI = nI.length;
275 twoBodyEnergy[i] = new double[roti.length][lenNI][];
276
277 for (int ri = 0; ri < roti.length; ri++) {
278 if (eR.check(i, ri)) {
279 continue;
280 }
281
282 for (int indJ = 0; indJ < lenNI; indJ++) {
283 int j = nI[indJ];
284 if (rO.checkNeighboringPair(i, j)) {
285 Residue resj = residues[j];
286 int indexJ = allResiduesList.indexOf(resj);
287 Rotamer[] rotj = resj.getRotamers();
288 twoBodyEnergy[i][ri][indJ] = new double[rotj.length];
289 for (int rj = 0; rj < rotj.length; rj++) {
290 if (eR.checkToJ(i, ri, j, rj)) {
291 continue;
292 }
293
294
295 if (dM.checkPairDistThreshold(indexI, ri, indexJ, rj)) {
296 continue;
297 }
298
299 Integer[] pairJob = {i, ri, j, rj};
300 if (decomposeOriginal && (ri != 0 || rj != 0)) {
301 continue;
302 }
303 twoBodyEnergyMap.put(pairJobIndex, pairJob);
304 if (reverseMap) {
305 String revKey = format("%d %d %d %d", i, ri, j, rj);
306 reverseJobMapPairs.put(revKey, pairJobIndex);
307 }
308 pairJobIndex++;
309 }
310 }
311 }
312 }
313 }
314 return reverseJobMapPairs;
315 }
316
317 public HashMap<String, Integer> allocate3BodyJobMap(Residue[] residues, int nResidues,
318 boolean reverseMap) {
319 HashMap<String, Integer> reverseJobMapTrimers = new HashMap<>();
320 threeBodyEnergyMap.clear();
321
322 threeBodyEnergy = new double[nResidues][][][][][];
323 int trimerJobIndex = 0;
324 for (int i = 0; i < nResidues; i++) {
325 Residue resi = residues[i];
326 int indexI = allResiduesList.indexOf(resi);
327 Rotamer[] roti = resi.getRotamers();
328 int lenri = roti.length;
329 int[] nI = resNeighbors[i];
330 int lenNI = nI.length;
331 threeBodyEnergy[i] = new double[lenri][lenNI][][][];
332
333 for (int ri = 0; ri < lenri; ri++) {
334 if (eR.check(i, ri)) {
335 continue;
336 }
337 for (int indJ = 0; indJ < lenNI; indJ++) {
338
339 int j = nI[indJ];
340 Residue resj = residues[j];
341 int indexJ = allResiduesList.indexOf(resj);
342 Rotamer[] rotj = resj.getRotamers();
343 int lenrj = rotj.length;
344 int[] nJ = resNeighbors[j];
345 int lenNJ = nJ.length;
346 threeBodyEnergy[i][ri][indJ] = new double[lenrj][lenNJ][];
347
348 for (int rj = 0; rj < lenrj; rj++) {
349 if (eR.checkToJ(i, ri, j, rj)) {
350 continue;
351 }
352
353 for (int indK = 0; indK < lenNJ; indK++) {
354 int k = nJ[indK];
355 Residue resk = residues[k];
356 int indexK = allResiduesList.indexOf(resk);
357 Rotamer[] rotk = resk.getRotamers();
358 int lenrk = rotk.length;
359 threeBodyEnergy[i][ri][indJ][rj][indK] = new double[lenrk];
360
361 for (int rk = 0; rk < lenrk; rk++) {
362 if (eR.checkToK(i, ri, j, rj, k, rk)) {
363 continue;
364 }
365 if (dM.checkTriDistThreshold(indexI, ri, indexJ, rj, indexK, rk)) {
366 continue;
367 }
368 Integer[] trimerJob = {i, ri, j, rj, k, rk};
369 if (decomposeOriginal && (ri != 0 || rj != 0 || rk != 0)) {
370 continue;
371 }
372 threeBodyEnergyMap.put(trimerJobIndex, trimerJob);
373 if (reverseMap) {
374 String revKey = format("%d %d %d %d %d %d", i, ri, j, rj, k, rk);
375 reverseJobMapTrimers.put(revKey, trimerJobIndex);
376 }
377 trimerJobIndex++;
378 }
379 }
380 }
381 }
382 }
383 }
384 return reverseJobMapTrimers;
385 }
386
387 public void allocate4BodyJobMap(Residue[] residues, int nResidues) {
388 logger.info(" Creating 4-Body jobs...");
389 fourBodyEnergyMap.clear();
390 boolean maxedOut = false;
391
392 int quadJobIndex = 0;
393 for (int i = 0; i < nResidues; i++) {
394 Residue resi = residues[i];
395 Rotamer[] roti = resi.getRotamers();
396 for (int ri = 0; ri < roti.length; ri++) {
397 if (eR.check(i, ri)) {
398 continue;
399 }
400 for (int j = i + 1; j < nResidues; j++) {
401 Residue resj = residues[j];
402 Rotamer[] rotj = resj.getRotamers();
403 for (int rj = 0; rj < rotj.length; rj++) {
404
405
406
407 if (eR.checkToJ(i, ri, j, rj)) {
408 continue;
409 }
410 for (int k = j + 1; k < nResidues; k++) {
411 Residue resk = residues[k];
412 Rotamer[] rotk = resk.getRotamers();
413 for (int rk = 0; rk < rotk.length; rk++) {
414
415
416
417 if (eR.checkToK(i, ri, j, rj, k, rk)) {
418 continue;
419 }
420 for (int l = k + 1; l < nResidues; l++) {
421 Residue resl = residues[l];
422 Rotamer[] rotl = resl.getRotamers();
423 for (int rl = 0; rl < rotl.length; rl++) {
424 if (eR.checkToL(i, ri, j, rj, k, rk, l, rl)) {
425 continue;
426 }
427 Integer[] quadJob = {i, ri, j, rj, k, rk, l, rl};
428 if (decomposeOriginal && (ri != 0 || rj != 0 || rk != 0 || rl != 0)) {
429 continue;
430 }
431 fourBodyEnergyMap.put(quadJobIndex++, quadJob);
432 if (fourBodyEnergyMap.size() > max4BodyCount) {
433 maxedOut = true;
434 break;
435 }
436 }
437 if (maxedOut) {
438 break;
439 }
440 }
441 if (maxedOut) {
442 break;
443 }
444 }
445 if (maxedOut) {
446 break;
447 }
448 }
449 if (maxedOut) {
450 break;
451 }
452 }
453 if (maxedOut) {
454 break;
455 }
456 }
457 if (maxedOut) {
458 break;
459 }
460 }
461 if (maxedOut) {
462 break;
463 }
464 }
465 }
466
467 public HashMap<String, Integer> allocateSelfJobMap(Residue[] residues, int nResidues,
468 boolean reverseMap) {
469 selfEnergyMap.clear();
470
471 HashMap<String, Integer> reverseJobMapSingles = new HashMap<>();
472 int singleJobIndex = 0;
473 selfEnergy = new double[nResidues][];
474 for (int i = 0; i < nResidues; i++) {
475 Residue resi = residues[i];
476 Rotamer[] roti = resi.getRotamers();
477 selfEnergy[i] = new double[roti.length];
478 for (int ri = 0; ri < roti.length; ri++) {
479 if (!eR.check(i, ri)) {
480 Integer[] selfJob = {i, ri};
481 if (decomposeOriginal && ri != 0) {
482 continue;
483 }
484 selfEnergyMap.put(singleJobIndex, selfJob);
485 if (reverseMap) {
486 String revKey = format("%d %d", i, ri);
487 reverseJobMapSingles.put(revKey, singleJobIndex);
488 }
489 singleJobIndex++;
490 }
491 }
492 }
493 return reverseJobMapSingles;
494 }
495
496
497
498
499
500
501
502
503
504
505
506
507 public double compute2BodyEnergy(Residue[] residues, int i, int ri, int j, int rj) {
508 rO.turnOffAllResidues(residues);
509 turnOnResidue(residues[i], ri);
510 turnOnResidue(residues[j], rj);
511 double energy;
512 try {
513 if (algorithmListener != null) {
514 algorithmListener.algorithmUpdate(molecularAssembly);
515 }
516 Rotamer[] rot_i = residues[i].getRotamers();
517 Rotamer[] rot_j = residues[j].getRotamers();
518 double subtract =
519 -backboneEnergy - getSelf(i, ri, rot_i[ri], true) - getSelf(j, rj, rot_j[rj], true);
520
521
522 energy = rO.currentEnergy(residues) + subtract;
523 if (potentialIsOpenMM && energy < ommRecalculateThreshold) {
524 logger.warning(
525 format(" Experimental: re-computing pair energy %s-%d %s-%d using Force Field X",
526 residues[i], ri, residues[j], rj));
527 energy = rO.currentFFXPE() + subtract;
528 }
529 if (energy < singularityThreshold) {
530 String message = format(
531 " Rejecting pair energy for %s-%d %s-%d is %10.5g << %10f, " + "likely in error.",
532 residues[i], ri, residues[j], rj, energy, singularityThreshold);
533 logger.warning(message);
534 throw new EnergyException(message, false, energy);
535 }
536 } finally {
537
538 turnOffResidue(residues[i]);
539 turnOffResidue(residues[j]);
540 }
541 return energy;
542 }
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558 public double compute3BodyEnergy(Residue[] residues, int i, int ri, int j, int rj, int k, int rk) {
559 turnOffAllResidues(residues);
560 turnOnResidue(residues[i], ri);
561 turnOnResidue(residues[j], rj);
562 turnOnResidue(residues[k], rk);
563 Rotamer[] rot_i = residues[i].getRotamers();
564 Rotamer[] rot_j = residues[j].getRotamers();
565 Rotamer[] rot_k = residues[k].getRotamers();
566 double energy;
567 try {
568 if (algorithmListener != null) {
569 algorithmListener.algorithmUpdate(molecularAssembly);
570 }
571 double subtract =
572 -backboneEnergy - getSelf(i, ri, rot_i[ri], true) - getSelf(j, rj, rot_j[rj], true)
573 - getSelf(k, rk, rot_k[rk], true) - get2Body(i, ri, j, rj) - get2Body(i, ri, k, rk)
574 - get2Body(j, rj, k, rk);
575 energy = rO.currentEnergy(residues) + subtract;
576 if (potentialIsOpenMM && energy < ommRecalculateThreshold) {
577 logger.warning(
578 format(" Experimental: re-computing triple energy %s-%d %s-%d %s-%d using Force Field X",
579 residues[i], ri, residues[j], rj, residues[k], rk));
580 energy = rO.currentFFXPE() + subtract;
581 }
582 if (energy < singularityThreshold) {
583 String message = format(" Rejecting triple energy for %s-%d %s-%d %s-%d is %10.5g << %10f, "
584 + "likely in error.", residues[i], ri, residues[j], rj, residues[k], rk, energy,
585 singularityThreshold);
586 logger.warning(message);
587 throw new EnergyException(message);
588 }
589 } finally {
590
591 turnOffResidue(residues[i]);
592 turnOffResidue(residues[j]);
593 turnOffResidue(residues[k]);
594 }
595 return energy;
596 }
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613 public double compute4BodyEnergy(Residue[] residues, int i, int ri, int j, int rj, int k, int rk,
614 int l, int rl) {
615 turnOffAllResidues(residues);
616 turnOnResidue(residues[i], ri);
617 turnOnResidue(residues[j], rj);
618 turnOnResidue(residues[k], rk);
619 turnOnResidue(residues[l], rl);
620 double energy;
621 try {
622 if (algorithmListener != null) {
623 algorithmListener.algorithmUpdate(molecularAssembly);
624 }
625 double subtract =
626 -backboneEnergy - getSelf(i, ri) - getSelf(j, rj) - getSelf(k, rk) - getSelf(l, rl)
627 - get2Body(i, ri, j, rj) - get2Body(i, ri, k, rk) - get2Body(i, ri, l, rl)
628 - get2Body(j, rj, k, rk) - get2Body(j, rj, l, rl) - get2Body(k, rk, l, rl)
629 - get3Body(residues, i, ri, j, rj, k, rk) - get3Body(residues, i, ri, j, rj, l, rl)
630 - get3Body(residues, i, ri, k, rk, l, rl) - get3Body(residues, j, rj, k, rk, l, rl);
631 energy = rO.currentEnergy(residues) + subtract;
632
633 if (potentialIsOpenMM && energy < ommRecalculateThreshold) {
634 logger.warning(format(
635 " Experimental: re-computing quad energy %s-%d %s-%d %s-%d %s-%d using Force Field X",
636 residues[i], ri, residues[j], rj, residues[k], rk, residues[l], rl));
637 energy = rO.currentFFXPE() + subtract;
638 }
639 if (energy < singularityThreshold) {
640 String message = format(
641 " Rejecting quad energy for %s-%d %s-%d %s-%d %s-%d is %10.5g << %10f, "
642 + "likely in error.", residues[i], ri, residues[j], rj, residues[k], rk, residues[l],
643 rl, energy, singularityThreshold);
644 logger.warning(message);
645 throw new EnergyException(message);
646 }
647
648 } finally {
649
650 turnOffResidue(residues[i]);
651 turnOffResidue(residues[j]);
652 turnOffResidue(residues[k]);
653 turnOffResidue(residues[l]);
654 }
655 return energy;
656 }
657
658
659
660
661
662
663
664
665
666
667
668
669
670 public double computeSelfEnergy(Residue[] residues, int i, int ri) {
671 rO.turnOffAllResidues(residues);
672 rO.turnOnResidue(residues[i], ri);
673 double energy;
674 try {
675 if (algorithmListener != null) {
676 algorithmListener.algorithmUpdate(molecularAssembly);
677 }
678 energy = rO.currentEnergy(residues) - backboneEnergy;
679 if (potentialIsOpenMM && energy < ommRecalculateThreshold) {
680 logger.warning(
681 format(" Experimental: re-computing self energy %s-%d using Force Field X", residues[i],
682 ri));
683 energy = rO.currentFFXPE() - backboneEnergy;
684 }
685 if (energy < singularityThreshold) {
686 String message = format(
687 " Rejecting self energy for %s-%d is %10.5g << %10f, " + "likely in error.", residues[i],
688 ri, energy, singularityThreshold);
689 logger.warning(message);
690 throw new EnergyException(message);
691 }
692 } finally {
693 rO.turnOffResidue(residues[i]);
694 }
695
696 Rotamer[] rotamers = residues[i].getRotamers();
697
698 if (rotamers[ri].isTitrating) {
699 double bias = rotamers[ri].getRotamerPhBias();
700 if (logger.isLoggable(Level.FINE)) {
701 logger.fine(format(" %s Self-Energy %16.8f = FF %16.8f - BB %16.8f + Ph Bias %16.8f",
702 rotamers[ri].getName(), energy + bias, energy + backboneEnergy, backboneEnergy, bias));
703 }
704 energy += bias;
705 }
706
707 return energy;
708 }
709
710
711
712
713
714
715
716
717 public static double getTotalRotamerPhBias(List<Residue> residues, int[] rotamers) {
718 double total = 0.0;
719 int n = residues.size();
720 for (int i = 0; i < n; i++) {
721 Rotamer[] rot = residues.get(i).getRotamers();
722 int ri = rotamers[i];
723 if (rot[ri].isTitrating) {
724 total += rot[ri].getRotamerPhBias();
725 }
726 }
727
728 return total;
729 }
730
731
732
733
734
735
736
737
738
739
740
741 public double get2Body(int i, int ri, int j, int rj) {
742
743 if (j < i) {
744 int ii = i;
745 int iri = ri;
746 i = j;
747 ri = rj;
748 j = ii;
749 rj = iri;
750 }
751 try {
752
753 int[] nI = resNeighbors[i];
754 int indJ = -1;
755 for (int l = 0; l < nI.length; l++) {
756 if (nI[l] == j) {
757 indJ = l;
758 break;
759 }
760 }
761 if (indJ == -1) {
762 logger.fine(format(" Residue %d not found in neighbors of %d; assumed past cutoff.", j, i));
763 return 0;
764 } else {
765 return twoBodyEnergy[i][ri][indJ][rj];
766 }
767 } catch (NullPointerException npe) {
768 logger.info(format(" NPE for 2-body energy (%3d,%2d) (%3d,%2d).", i, ri, j, rj));
769 throw npe;
770 }
771 }
772
773
774
775
776
777
778
779
780
781
782
783
784
785 public double get3Body(Residue[] residues, int i, int ri, int j, int rj, int k, int rk) {
786 if (!threeBodyTerm) {
787 return 0.0;
788 }
789
790 if (j < i) {
791 int ii = i;
792 int iri = ri;
793 i = j;
794 ri = rj;
795 j = ii;
796 rj = iri;
797 }
798 if (k < i) {
799 int ii = i;
800 int iri = ri;
801 i = k;
802 ri = rk;
803 k = ii;
804 rk = iri;
805 }
806 if (k < j) {
807 int jj = j;
808 int jrj = rj;
809 j = k;
810 rj = rk;
811 k = jj;
812 rk = jrj;
813 }
814
815
816 int[] nI = resNeighbors[i];
817 int indJ = -1;
818 for (int l = 0; l < nI.length; l++) {
819 if (nI[l] == j) {
820 indJ = l;
821 break;
822 }
823 }
824
825 int[] nJ = resNeighbors[j];
826 int indK = -1;
827 for (int l = 0; l < nJ.length; l++) {
828 if (nJ[l] == k) {
829 indK = l;
830 break;
831 }
832 }
833
834
835
836
837 int indexI = allResiduesList.indexOf(residues[i]);
838 int indexJ = allResiduesList.indexOf(residues[j]);
839 int indexK = allResiduesList.indexOf(residues[k]);
840 if (dM.checkTriDistThreshold(indexI, ri, indexJ, rj, indexK, rk)) {
841 return 0;
842 } else {
843 try {
844 return threeBodyEnergy[i][ri][indJ][rj][indK][rk];
845 } catch (NullPointerException | ArrayIndexOutOfBoundsException ex) {
846 String message = format(
847 " Could not find an energy for 3-body energy (%3d,%2d) (%3d,%2d) (%3d,%2d)", i, ri, j,
848 rj, k, rk);
849 logger.warning(message);
850 throw ex;
851 }
852 }
853 }
854
855 public double getBackboneEnergy() {
856 return backboneEnergy;
857 }
858
859 public void setBackboneEnergy(double backboneEnergy) {
860 this.backboneEnergy = backboneEnergy;
861 }
862
863 public Map<Integer, Integer[]> getFourBodyEnergyMap() {
864 return fourBodyEnergyMap;
865 }
866
867
868
869
870
871
872
873
874 public double getSelf(int i, int ri) {
875 try {
876 return selfEnergy[i][ri];
877 } catch (NullPointerException npe) {
878 logger.info(format(" NPE for self energy (%3d,%2d).", i, ri));
879 throw npe;
880 }
881 }
882
883
884
885
886
887
888
889
890 public double getSelf(int i, int ri, Rotamer rot, boolean excludeFMod) {
891 try {
892 double totalSelf;
893 if (rot.isTitrating && excludeFMod) {
894 totalSelf = selfEnergy[i][ri] - rot.getRotamerPhBias();
895 } else {
896 totalSelf = selfEnergy[i][ri];
897 }
898 return totalSelf;
899 } catch (NullPointerException npe) {
900 logger.info(format(" NPE for self energy (%3d,%2d).", i, ri));
901 throw npe;
902 }
903 }
904
905
906 public Map<Integer, Integer[]> getSelfEnergyMap() {
907 return selfEnergyMap;
908 }
909
910 public Map<Integer, Integer[]> getThreeBodyEnergyMap() {
911 return threeBodyEnergyMap;
912 }
913
914 public Map<Integer, Integer[]> getTwoBodyEnergyMap() {
915 return twoBodyEnergyMap;
916 }
917
918 public int loadEnergyRestart(File restartFile, Residue[] residues) {
919 return loadEnergyRestart(restartFile, residues, -1, null);
920 }
921
922 public int loadEnergyRestart(File restartFile, Residue[] residues, int boxIteration,
923 int[] cellIndices) {
924 try {
925 int nResidues = residues.length;
926 Path path = Paths.get(restartFile.getCanonicalPath());
927 List<String> lines = Files.readAllLines(path, StandardCharsets.UTF_8);
928 List<String> linesThisBox = new ArrayList<>();
929
930 try {
931 backboneEnergy = rO.computeBackboneEnergy(residues);
932 } catch (ArithmeticException ex) {
933 logger.severe(
934 format(" Exception %s in calculating backbone energy; FFX shutting down.", ex));
935 }
936 rO.logIfRank0(format("\n Backbone energy: %s\n", rO.formatEnergy(backboneEnergy)));
937
938 if (usingBoxOptimization && boxIteration >= 0) {
939 boolean foundBox = false;
940 for (int i = 0; i < lines.size(); i++) {
941 String line = lines.get(i);
942 if (line.startsWith("Box")) {
943 String[] tok = line.replaceAll("Box", "").replaceAll(":", ",").replaceAll(" ", "")
944 .split(",");
945 int readIteration = Integer.parseInt(tok[0]);
946 int readCellIndexX = Integer.parseInt(tok[1]);
947 int readCellIndexY = Integer.parseInt(tok[2]);
948 int readCellIndexZ = Integer.parseInt(tok[3]);
949 if (readIteration == boxIteration && readCellIndexX == cellIndices[0]
950 && readCellIndexY == cellIndices[1] && readCellIndexZ == cellIndices[2]) {
951 foundBox = true;
952 for (int j = i + 1; j < lines.size(); j++) {
953 String l = lines.get(j);
954 if (l.startsWith("Box")) {
955 break;
956 }
957 linesThisBox.add(l);
958 }
959 break;
960 }
961 }
962 }
963 if (!foundBox) {
964 rO.logIfRank0(format(" Didn't find restart energies for Box %d: %d,%d,%d", boxIteration,
965 cellIndices[0], cellIndices[1], cellIndices[2]));
966 return 0;
967 } else if (linesThisBox.isEmpty()) {
968 return 0;
969 } else {
970 lines = linesThisBox;
971 }
972 }
973
974 List<String> singleLines = new ArrayList<>();
975 List<String> pairLines = new ArrayList<>();
976 List<String> tripleLines = new ArrayList<>();
977 for (String line : lines) {
978 String[] tok = line.split("\\s");
979 if (tok[0].startsWith("Self")) {
980 singleLines.add(line);
981 } else if (tok[0].startsWith("Pair")) {
982 pairLines.add(line);
983 } else if (tok[0].startsWith("Triple")) {
984 tripleLines.add(line);
985 }
986 }
987 int loaded = 0;
988 if (!tripleLines.isEmpty()) {
989 loaded = 3;
990 } else if (!pairLines.isEmpty()) {
991 loaded = 2;
992 } else if (!singleLines.isEmpty()) {
993 loaded = 1;
994 } else {
995 logger.warning(
996 format(" Empty or unreadable energy restart file: %s.", restartFile.getCanonicalPath()));
997 }
998 if (loaded >= 1) {
999 boolean reverseMap = true;
1000 HashMap<String, Integer> reverseJobMapSingles = allocateSelfJobMap(residues, nResidues,
1001 reverseMap);
1002
1003 for (String line : singleLines) {
1004 try {
1005 String[] tok = line.replace(",", "").replace(":", "").split("\\s+");
1006 int i;
1007 if (tok[1].contains("-")) {
1008 i = nameToNumber(tok[1], residues);
1009 } else {
1010 i = Integer.parseInt(tok[1]);
1011 }
1012 int ri = Integer.parseInt(tok[2]);
1013 double energy = Double.parseDouble(tok[3]);
1014 try {
1015 setSelf(i, ri, energy);
1016 if (verbose) {
1017 rO.logIfRank0(format(" From restart file: Self energy %3d (%8s,%2d): %s", i,
1018 residues[i].toFormattedString(false, true), ri, rO.formatEnergy(energy)));
1019 }
1020 } catch (Exception e) {
1021 if (verbose) {
1022 rO.logIfRank0(format(" Restart file out-of-bounds index: %s", line));
1023 }
1024 }
1025
1026 String revKey = format("%d %d", i, ri);
1027 selfEnergyMap.remove(reverseJobMapSingles.get(revKey));
1028 } catch (NumberFormatException ex) {
1029 logger.log(Level.WARNING, format(" Unparsable line in energy restart file: \n%s", line),
1030 ex);
1031 }
1032 }
1033 rO.logIfRank0(" Loaded self energies from restart file.");
1034
1035
1036 eR.prePruneSelves(residues);
1037
1038
1039 if (pruneClashes) {
1040 eR.pruneSingleClashes(residues);
1041 }
1042 }
1043
1044
1045 condenseEnergyMap(selfEnergyMap);
1046
1047 if (loaded >= 2) {
1048 if (!selfEnergyMap.isEmpty()) {
1049 rO.logIfRank0(
1050 " Double-check that parameters match original run due to missing self-energies.");
1051 }
1052 boolean reverseMap = true;
1053 HashMap<String, Integer> reverseJobMapPairs = allocate2BodyJobMap(residues, nResidues,
1054 reverseMap);
1055
1056
1057 for (String line : pairLines) {
1058 try {
1059 String[] tok = line.replace(",", "").replace(":", "").split("\\s+");
1060 int i;
1061 if (tok[1].contains("-")) {
1062 i = nameToNumber(tok[1], residues);
1063 } else {
1064 i = Integer.parseInt(tok[1]);
1065 }
1066 int ri = Integer.parseInt(tok[2]);
1067 int j;
1068 if (tok[3].contains("-")) {
1069 j = nameToNumber(tok[3], residues);
1070 } else {
1071 j = Integer.parseInt(tok[3]);
1072 }
1073 int rj = Integer.parseInt(tok[4]);
1074 double energy = Double.parseDouble(tok[5]);
1075 try {
1076
1077
1078
1079
1080 if (rO.checkNeighboringPair(i, j)) {
1081
1082
1083 Residue residueI = residues[i];
1084 Residue residueJ = residues[j];
1085 int indexI = allResiduesList.indexOf(residueI);
1086 int indexJ = allResiduesList.indexOf(residueJ);
1087 if (!dM.checkPairDistThreshold(indexI, ri, indexJ, rj)) {
1088 set2Body(i, ri, j, rj, energy);
1089
1090 double resDist = dM.getResidueDistance(indexI, ri, indexJ, rj);
1091 String resDistString = "large";
1092 if (resDist < Double.MAX_VALUE) {
1093 resDistString = format("%5.3f", resDist);
1094 }
1095
1096 double dist = dM.checkDistMatrix(indexI, ri, indexJ, rj);
1097 String distString = " large";
1098 if (dist < Double.MAX_VALUE) {
1099 distString = format("%10.3f", dist);
1100 }
1101
1102 logger.fine(format(" Pair %8s %-2d, %8s %-2d: %s at %s Ang (%s Ang by residue).",
1103 residueI.toFormattedString(false, true), ri,
1104 residueJ.toFormattedString(false, true), rj,
1105 rO.formatEnergy(get2Body(i, ri, j, rj)), distString, resDistString));
1106 }
1107 } else {
1108 logger.fine(format(
1109 "Ignoring a pair-energy from outside the cutoff: 2-energy [(%8s,%2d),(%8s,%2d)]: %12.4f",
1110 residues[i].toFormattedString(false, true), ri,
1111 residues[j].toFormattedString(false, true), rj, energy));
1112 }
1113
1114 if (verbose) {
1115 rO.logIfRank0(
1116 format(" From restart file: Pair energy [(%8s,%2d),(%8s,%2d)]: %12.4f",
1117 residues[i].toFormattedString(false, true), ri,
1118 residues[j].toFormattedString(false, true), rj, energy));
1119 }
1120 } catch (Exception e) {
1121 if (verbose) {
1122 rO.logIfRank0(format(" Restart file out-of-bounds index: %s", line));
1123 }
1124 }
1125
1126 String revKey = format("%d %d %d %d", i, ri, j, rj);
1127 twoBodyEnergyMap.remove(reverseJobMapPairs.get(revKey));
1128 } catch (NumberFormatException ex) {
1129 logger.log(Level.WARNING, format("Unparsable line in energy restart file: \n%s", line),
1130 ex);
1131 }
1132 }
1133 rO.logIfRank0(" Loaded 2-body energies from restart file.");
1134
1135
1136 eR.prePrunePairs(residues);
1137
1138
1139 if (prunePairClashes) {
1140 eR.prunePairClashes(residues);
1141 }
1142 }
1143
1144
1145 condenseEnergyMap(twoBodyEnergyMap);
1146
1147 if (loaded >= 3) {
1148 if (!twoBodyEnergyMap.isEmpty()) {
1149 if (master) {
1150 logger.warning(
1151 "Double-check that parameters match original run! Found trimers in restart file, but pairs job queue is non-empty.");
1152 }
1153 }
1154 boolean reverseMap = true;
1155 HashMap<String, Integer> reverseJobMapTrimers = allocate3BodyJobMap(residues, nResidues,
1156 reverseMap);
1157
1158
1159
1160 for (String line : tripleLines) {
1161 try {
1162 String[] tok = line.replace(",", "").replace(":", "").split("\\s+");
1163 int i;
1164 if (tok[1].contains("-")) {
1165 i = nameToNumber(tok[1], residues);
1166 } else {
1167 i = Integer.parseInt(tok[1]);
1168 }
1169 int ri = Integer.parseInt(tok[2]);
1170 int j;
1171 if (tok[3].contains("-")) {
1172 j = nameToNumber(tok[3], residues);
1173 } else {
1174 j = Integer.parseInt(tok[3]);
1175 }
1176 int rj = Integer.parseInt(tok[4]);
1177 int k;
1178 if (tok[5].contains("-")) {
1179 k = nameToNumber(tok[5], residues);
1180 } else {
1181 k = Integer.parseInt(tok[5]);
1182 }
1183 int rk = Integer.parseInt(tok[6]);
1184 double energy = Double.parseDouble(tok[7]);
1185
1186 try {
1187
1188
1189
1190
1191
1192
1193 if (rO.checkNeighboringTriple(i, j, k)) {
1194
1195
1196 Residue residueI = residues[i];
1197 Residue residueJ = residues[j];
1198 Residue residueK = residues[k];
1199 int indexI = allResiduesList.indexOf(residueI);
1200 int indexJ = allResiduesList.indexOf(residueJ);
1201 int indexK = allResiduesList.indexOf(residueK);
1202 if (!dM.checkTriDistThreshold(indexI, ri, indexJ, rj, indexK, rk)) {
1203 set3Body(residues, i, ri, j, rj, k, rk, energy);
1204
1205 double resDist = dM.get3BodyResidueDistance(indexI, ri, indexJ, rj, indexK, rk);
1206 String resDistString = " large";
1207 if (resDist < Double.MAX_VALUE) {
1208 resDistString = format("%5.3f", resDist);
1209 }
1210
1211 double rawDist = dM.getRawNBodyDistance(indexI, ri, indexJ, rj, indexK, rk);
1212 String distString = " large";
1213 if (rawDist < Double.MAX_VALUE) {
1214 distString = format("%10.3f", rawDist);
1215 }
1216
1217 logger.fine(format(
1218 " 3-Body %8s %-2d, %8s %-2d, %8s %-2d: %s at %s Ang (%s Ang by residue).",
1219 residueI.toFormattedString(false, true), ri,
1220 residueJ.toFormattedString(false, true), rj,
1221 residueK.toFormattedString(false, true), rk,
1222 rO.formatEnergy(get3Body(residues, i, ri, j, rj, k, rk)), distString,
1223 resDistString));
1224 }
1225 } else {
1226 logger.fine(format(
1227 "Ignoring a triple-energy from outside the cutoff: 3-Body %8s %-2d, %8s %-2d, %8s %-2d: %s",
1228 residues[i].toFormattedString(false, true), ri,
1229 residues[j].toFormattedString(false, true), rj,
1230 residues[k].toFormattedString(false, true), rk,
1231 rO.formatEnergy(get3Body(residues, i, ri, j, rj, k, rk))));
1232 }
1233 } catch (ArrayIndexOutOfBoundsException ex) {
1234 if (verbose) {
1235 rO.logIfRank0(format(" Restart file out-of-bounds index: %s", line));
1236 }
1237 } catch (NullPointerException npe) {
1238 if (verbose) {
1239 rO.logIfRank0(format(" NPE in loading 3-body energies: pruning "
1240 + "likely changed! 3-body %s-%d %s-%d %s-%d",
1241 residues[i].toFormattedString(false, true), ri, residues[j], rj, residues[k],
1242 rk));
1243 }
1244 }
1245 if (verbose) {
1246 rO.logIfRank0(
1247 format(" From restart file: Trimer energy %3d %-2d, %3d %-2d, %3d %-2d: %s", i, ri,
1248 j, rj, k, rk, rO.formatEnergy(energy)));
1249 }
1250
1251 String revKey = format("%d %d %d %d %d %d", i, ri, j, rj, k, rk);
1252 threeBodyEnergyMap.remove(reverseJobMapTrimers.get(revKey));
1253 } catch (NumberFormatException ex) {
1254 logger.log(Level.WARNING, format("Unparsable line in energy restart file: \n%s", line),
1255 ex);
1256 }
1257 }
1258 rO.logIfRank0(" Loaded trimer energies from restart file.");
1259 }
1260
1261
1262 condenseEnergyMap(threeBodyEnergyMap);
1263
1264 return loaded;
1265 } catch (IOException ex) {
1266 logger.log(Level.WARNING, "Exception while loading energy restart file.", ex);
1267 }
1268
1269 return 0;
1270 }
1271
1272
1273
1274
1275
1276
1277
1278
1279
1280
1281 public double lowestPairEnergy(Residue[] residues, int i, int ri, int j) {
1282 if (residues == null) {
1283 return 0.0;
1284 }
1285 int n = residues.length;
1286 if (i < 0 || i >= n) {
1287 return 0.0;
1288 }
1289 if (j < 0 || j >= n) {
1290 return 0.0;
1291 }
1292
1293 Rotamer[] rotamers = residues[j].getRotamers();
1294 int nr = rotamers.length;
1295 double energy = Double.MAX_VALUE;
1296 for (int jr = 0; jr < nr; jr++) {
1297 try {
1298 double e = get2Body(i, ri, j, jr);
1299 if (e < energy) {
1300 energy = e;
1301 }
1302 } catch (Exception e) {
1303
1304 }
1305 }
1306 return energy;
1307 }
1308
1309
1310
1311
1312
1313
1314
1315
1316 public double lowestSelfEnergy(Residue[] residues, int i) {
1317 if (residues == null) {
1318 return 0.0;
1319 }
1320 int n = residues.length;
1321 if (i < 0 || i >= n) {
1322 return 0.0;
1323 }
1324 Rotamer[] rotamers = residues[i].getRotamers();
1325 int nr = rotamers.length;
1326 double energy = Double.MAX_VALUE;
1327 for (int ni = 0; ni < nr; ni++) {
1328 try {
1329 double e = getSelf(i, ni);
1330 if (e < energy) {
1331 energy = e;
1332 }
1333 } catch (Exception e) {
1334
1335 }
1336 }
1337 return energy;
1338 }
1339
1340
1341
1342
1343
1344
1345
1346
1347
1348
1349
1350
1351
1352 public boolean minMaxE2(Residue[] residues, double[] minMax, int i, int ri, int j, int rj)
1353 throws IllegalArgumentException {
1354 Residue resi = residues[i];
1355 Residue resj = residues[j];
1356 if (eR.check(i, ri) || eR.check(j, rj) || eR.check(i, ri, j, rj)) {
1357 throw new IllegalArgumentException(
1358 format(" Called for minMaxE2 on an eliminated pair %s-%d %s-%d",
1359 resi.toFormattedString(false, true), ri, resj.toFormattedString(false, true), rj));
1360 }
1361
1362
1363 minMax[0] = 0;
1364
1365 minMax[1] = 0;
1366
1367 int nRes = residues.length;
1368 for (int k = 0; k < nRes; k++) {
1369 if (k == i || k == j) {
1370 continue;
1371 }
1372 Residue resk = residues[k];
1373 Rotamer[] rotsk = resk.getRotamers();
1374 int lenrk = rotsk.length;
1375 double[] minMaxK = new double[2];
1376 minMaxK[0] = Double.MAX_VALUE;
1377 minMaxK[1] = Double.MIN_VALUE;
1378
1379 for (int rk = 0; rk < lenrk; rk++) {
1380 if (eR.check(k, rk)) {
1381
1382 continue;
1383 }
1384 if (eR.check(i, ri, k, rk) || eR.check(j, rj, k, rk)) {
1385
1386
1387
1388
1389 minMaxK[1] = Double.NaN;
1390 } else {
1391
1392
1393
1394 double currentMin = get2Body(i, ri, k, rk) + get2Body(j, rj, k, rk);
1395 double currentMax = currentMin;
1396 if (threeBodyTerm) {
1397
1398 currentMin += get3Body(residues, i, ri, j, rj, k, rk);
1399 currentMax = currentMin;
1400
1401
1402 double[] minMaxTriple = new double[2];
1403 if (minMaxE3(residues, minMaxTriple, i, ri, j, rj, k, rk)) {
1404
1405 assert (Double.isFinite(minMaxTriple[0]) && minMaxTriple[0] != Double.MAX_VALUE);
1406
1407
1408 currentMin += minMaxTriple[0];
1409
1410 if (Double.isFinite(currentMax) && Double.isFinite(minMaxTriple[1])) {
1411 currentMax += minMaxTriple[1];
1412 } else {
1413 currentMax = Double.NaN;
1414 }
1415 } else {
1416
1417 currentMin = Double.NaN;
1418 currentMax = Double.NaN;
1419 }
1420 }
1421
1422 assert (threeBodyTerm || currentMax == currentMin);
1423
1424
1425 if (Double.isFinite(currentMin) && currentMin < minMaxK[0]) {
1426
1427 minMaxK[0] = currentMin;
1428 }
1429
1430 if (Double.isFinite(currentMax) && Double.isFinite(minMaxK[1])) {
1431
1432 minMaxK[1] = Math.max(currentMax, minMaxK[1]);
1433 } else {
1434
1435 minMaxK[1] = Double.NaN;
1436 }
1437 }
1438 }
1439
1440 if (Double.isFinite(minMaxK[0])) {
1441
1442 minMax[0] += minMaxK[0];
1443 } else {
1444
1445 minMax[0] = Double.NaN;
1446 minMax[1] = Double.NaN;
1447 return false;
1448 }
1449 if (Double.isFinite(minMaxK[1]) && Double.isFinite(minMax[1])) {
1450
1451 minMax[1] += minMaxK[1];
1452 } else {
1453
1454 minMax[1] = Double.NaN;
1455 }
1456 }
1457
1458 return Double.isFinite(minMax[0]);
1459 }
1460
1461
1462
1463
1464
1465
1466
1467
1468
1469
1470
1471
1472
1473
1474 public boolean minMaxPairEnergy(Residue[] residues, double[] minMax, int i, int ri, int j) {
1475 Residue residuej = residues[j];
1476 Rotamer[] rotamersj = residuej.getRotamers();
1477 int lenrj = rotamersj.length;
1478 boolean valid = false;
1479 minMax[0] = Double.MAX_VALUE;
1480 minMax[1] = Double.MIN_VALUE;
1481
1482
1483 for (int rj = 0; rj < lenrj; rj++) {
1484
1485 if (eR.check(i, ri) || eR.check(j, rj) || eR.check(i, ri, j, rj)) {
1486 continue;
1487 }
1488
1489 double currMax = get2Body(i, ri, j, rj);
1490 double currMin = currMax;
1491
1492 if (threeBodyTerm) {
1493 double[] minMaxTriple = new double[2];
1494
1495 boolean validPair = minMax2BodySum(residues, minMaxTriple, i, ri, j, rj);
1496 if (!validPair) {
1497
1498 Residue residuei = residues[i];
1499 rO.logIfRank0(format(" Inconsistent Pair: %8s %2d, %8s %2d.",
1500 residuei.toFormattedString(false, true), ri, residuej.toFormattedString(false, true),
1501 rj), Level.INFO);
1502 continue;
1503 }
1504
1505 if (Double.isFinite(currMin) && Double.isFinite(minMaxTriple[0])) {
1506 currMin += minMaxTriple[0];
1507 } else {
1508 currMin = Double.NaN;
1509 }
1510
1511 if (Double.isFinite(currMax) && Double.isFinite(minMaxTriple[1])) {
1512 currMax += minMaxTriple[1];
1513 } else {
1514 currMax = Double.NaN;
1515 }
1516 }
1517
1518 valid = true;
1519 if (Double.isFinite(currMin) && currMin < minMax[0]) {
1520 minMax[0] = currMin;
1521 }
1522
1523 if (Double.isFinite(currMax) && Double.isFinite(minMax[1])) {
1524 if (currMax > minMax[1]) {
1525
1526 minMax[1] = currMax;
1527 }
1528 } else {
1529
1530 minMax[1] = Double.NaN;
1531 }
1532 }
1533
1534
1535
1536 minMax[0] = (minMax[0] == Double.MAX_VALUE) ? Double.NaN : minMax[0];
1537
1538 return valid;
1539 }
1540
1541 public void set2Body(int i, int ri, int j, int rj, double e) {
1542 set2Body(i, ri, j, rj, e, false);
1543 }
1544
1545
1546
1547
1548
1549
1550
1551
1552
1553
1554
1555 public void set2Body(int i, int ri, int j, int rj, double e, boolean quiet) {
1556
1557 if (j < i) {
1558 int ii = i;
1559 int iri = ri;
1560 i = j;
1561 ri = rj;
1562 j = ii;
1563 rj = iri;
1564 }
1565 try {
1566
1567 int[] nI = resNeighbors[i];
1568 int indJ = -1;
1569 for (int l = 0; l < nI.length; l++) {
1570 if (nI[l] == j) {
1571 indJ = l;
1572 break;
1573 }
1574 }
1575 if (indJ == -1) {
1576 throw new IllegalArgumentException(
1577 format(" Residue %d not found in neighbors of %d; assumed past cutoff.", j, i));
1578 } else {
1579 twoBodyEnergy[i][ri][indJ][rj] = e;
1580 }
1581 } catch (NullPointerException npe) {
1582 if (!quiet) {
1583 logger.info(format(" NPE for 2-body energy (%3d,%2d) (%3d,%2d).", i, ri, j, rj));
1584 }
1585 throw npe;
1586 }
1587 }
1588
1589
1590
1591
1592
1593
1594
1595
1596
1597
1598
1599
1600
1601 public void set3Body(Residue[] residues, int i, int ri, int j, int rj, int k, int rk, double e) {
1602 set3Body(residues, i, ri, j, rj, k, rk, e, false);
1603 }
1604
1605
1606
1607
1608
1609
1610
1611
1612
1613
1614
1615
1616
1617
1618
1619 public void set3Body(Residue[] residues, int i, int ri, int j, int rj, int k, int rk, double e,
1620 boolean quiet) throws IllegalStateException {
1621 if (!threeBodyTerm) {
1622 throw new IllegalStateException(
1623 " Attempting to set a 3-body energy when threeBodyTerm is false!");
1624 }
1625
1626 if (j < i) {
1627 int ii = i;
1628 int iri = ri;
1629 i = j;
1630 ri = rj;
1631 j = ii;
1632 rj = iri;
1633 }
1634 if (k < i) {
1635 int ii = i;
1636 int iri = ri;
1637 i = k;
1638 ri = rk;
1639 k = ii;
1640 rk = iri;
1641 }
1642 if (k < j) {
1643 int jj = j;
1644 int jrj = rj;
1645 j = k;
1646 rj = rk;
1647 k = jj;
1648 rk = jrj;
1649 }
1650
1651
1652 int[] nI = resNeighbors[i];
1653 int indJ = -1;
1654 for (int l = 0; l < nI.length; l++) {
1655 if (nI[l] == j) {
1656 indJ = l;
1657 break;
1658 }
1659 }
1660
1661 int[] nJ = resNeighbors[j];
1662 int indK = -1;
1663 for (int l = 0; l < nJ.length; l++) {
1664 if (nJ[l] == k) {
1665 indK = l;
1666 break;
1667 }
1668 }
1669
1670
1671
1672
1673 int indexI = allResiduesList.indexOf(residues[i]);
1674 int indexJ = allResiduesList.indexOf(residues[j]);
1675 int indexK = allResiduesList.indexOf(residues[k]);
1676 if (dM.checkTriDistThreshold(indexI, ri, indexJ, rj, indexK, rk)) {
1677 throw new IllegalArgumentException(
1678 format(" Residue %d not found in neighbors of %d; assumed past cutoff.", j, i));
1679 } else {
1680 try {
1681 threeBodyEnergy[i][ri][indJ][rj][indK][rk] = e;
1682 } catch (NullPointerException | ArrayIndexOutOfBoundsException ex) {
1683 if (!quiet) {
1684 String message = format(
1685 " Could not access threeBodyEnergy array for 3-body energy (%3d,%2d) (%3d,%2d) (%3d,%2d)",
1686 i, ri, j, rj, k, rk);
1687 logger.warning(message);
1688 }
1689 throw ex;
1690 }
1691 }
1692 }
1693
1694 public void setSelf(int i, int ri, double e) {
1695 setSelf(i, ri, e, false);
1696 }
1697
1698
1699
1700
1701
1702
1703
1704
1705
1706 public void setSelf(int i, int ri, double e, boolean quiet) {
1707 try {
1708 selfEnergy[i][ri] = e;
1709 } catch (NullPointerException | ArrayIndexOutOfBoundsException ex) {
1710 if (!quiet) {
1711 logger.warning(format(" NPE or array index error for (%3d,%2d)", i, ri));
1712 }
1713 throw ex;
1714 }
1715 }
1716
1717 public void turnOffAllResidues(Residue[] residues) {
1718 if (residues == null) {
1719 return;
1720 }
1721 for (Residue residue : residues) {
1722 turnOffResidue(residue);
1723 }
1724 }
1725
1726 public void turnOffResidue(Residue residue) {
1727 turnOffAtoms(residue);
1728 applyDefaultRotamer(residue);
1729 }
1730
1731 public void turnOnAllResidues(Residue[] residues) {
1732 if (residues == null) {
1733 return;
1734 }
1735 for (Residue residue : residues) {
1736 turnOnAtoms(residue);
1737 }
1738 }
1739
1740 public void turnOnResidue(Residue residue, int ri) {
1741 turnOnAtoms(residue);
1742 Rotamer[] rotamers = residue.getRotamers();
1743 applyRotamer(residue, rotamers[ri]);
1744 }
1745
1746
1747
1748
1749
1750
1751
1752
1753
1754
1755
1756
1757 private boolean minMax2BodySum(Residue[] residues, double[] minMax, int i, int ri, int j, int rj) {
1758 int nres = residues.length;
1759 double minSum = 0.0;
1760 double maxSum = 0.0;
1761 for (int k = 0; k < nres; k++) {
1762 if (k == i || k == j) {
1763 continue;
1764 }
1765 Residue residuek = residues[k];
1766 Rotamer[] romatersk = residuek.getRotamers();
1767 int lenrk = romatersk.length;
1768 double currentMin = Double.MAX_VALUE;
1769 double currentMax = Double.MIN_VALUE;
1770 for (int rk = 0; rk < lenrk; rk++) {
1771 if (eR.check(k, rk)) {
1772
1773 continue;
1774 }
1775 if (eR.check(i, ri, k, rk) || eR.check(j, rj, k, rk)) {
1776
1777
1778 currentMax = Double.NaN;
1779 } else {
1780 double current = get3Body(residues, i, ri, j, rj, k, rk);
1781 if (Double.isFinite(current) && current < currentMin) {
1782 currentMin = current;
1783 }
1784 if (Double.isFinite(current) && Double.isFinite(currentMax)) {
1785 if (current > currentMax) {
1786 currentMax = current;
1787 }
1788 } else {
1789
1790 currentMax = Double.NaN;
1791 }
1792 }
1793 }
1794 if (currentMin == Double.MAX_VALUE || !Double.isFinite(minSum)) {
1795
1796
1797 minMax[0] = Double.NaN;
1798 minMax[1] = Double.NaN;
1799 return false;
1800 } else {
1801
1802 minSum += currentMin;
1803 }
1804 if (Double.isFinite(maxSum) && Double.isFinite(currentMax)) {
1805 maxSum += currentMax;
1806 } else {
1807 maxSum = Double.NaN;
1808 }
1809 }
1810 minMax[0] = minSum;
1811 minMax[1] = maxSum;
1812 return true;
1813 }
1814
1815
1816
1817
1818
1819
1820
1821
1822
1823
1824
1825
1826
1827
1828
1829
1830 private boolean minMaxE3(Residue[] residues, double[] minMax, int i, int ri, int j, int rj, int k,
1831 int rk) throws IllegalArgumentException {
1832 Residue resi = residues[i];
1833 Residue resj = residues[j];
1834 Residue resk = residues[k];
1835 if (eR.check(i, ri) || eR.check(j, rj) || eR.check(k, rk) || eR.check(i, ri, j, rj)
1836 || eR.check(i, ri, k, rk) || eR.check(j, rj, k, rk)) {
1837
1838 throw new IllegalArgumentException(
1839 format(" Called for minMaxE2 on an eliminated triple %s-%d %s-%d %s-%d",
1840 resi.toFormattedString(false, true), ri, resj.toFormattedString(false, true), rj,
1841 resk.toFormattedString(false, true), rk));
1842 }
1843
1844
1845 minMax[0] = 0;
1846 minMax[1] = 0;
1847 int nRes = residues.length;
1848 for (int l = 0; l < nRes; l++) {
1849 if (l == i || l == j || l == k) {
1850 continue;
1851 }
1852 Residue resl = residues[l];
1853 Rotamer[] rotsl = resl.getRotamers();
1854 int lenrl = rotsl.length;
1855
1856
1857 double currentMax = Double.MIN_VALUE;
1858 double currentMin = Double.MAX_VALUE;
1859
1860 for (int rl = 0; rl < lenrl; rl++) {
1861 if (eR.check(l, rl) || eR.check(k, rk, l, rl)) {
1862
1863 continue;
1864 }
1865
1866 double current;
1867 if (eR.check(i, ri, l, rl) || eR.check(j, rj, l, rl)) {
1868
1869 current = Double.NaN;
1870 } else {
1871
1872 current =
1873 get3Body(residues, i, ri, k, rk, l, rl) + get3Body(residues, j, rj, k, rk, l, rl);
1874 }
1875
1876
1877
1878
1879
1880 if (Double.isFinite(current) && current < currentMin) {
1881
1882 currentMin = current;
1883 }
1884
1885 if (Double.isFinite(current) && Double.isFinite(currentMax)) {
1886 if (current > currentMax) {
1887 currentMax = current;
1888 }
1889 } else {
1890 currentMax = Double.NaN;
1891 }
1892 }
1893
1894 if (Double.isFinite(currentMin)) {
1895 minMax[0] += currentMin;
1896 } else {
1897
1898 minMax[0] = Double.NaN;
1899 minMax[1] = Double.NaN;
1900 return false;
1901 }
1902
1903 if (Double.isFinite(currentMax) && Double.isFinite(minMax[1])) {
1904 minMax[1] += currentMax;
1905 } else {
1906 minMax[1] = Double.NaN;
1907 }
1908
1909 }
1910 return Double.isFinite(minMax[0]);
1911 }
1912
1913
1914
1915
1916
1917
1918 private void applyDefaultRotamer(Residue residue) {
1919 applyRotamer(residue, residue.getRotamers()[0]);
1920 }
1921
1922 private int nameToNumber(String residueString, Residue[] residues) throws NumberFormatException {
1923 int ret = -1;
1924 for (int x = 0; x < residues.length; x++) {
1925 if (residueString.equals(residues[x].toString())) {
1926 ret = x;
1927 break;
1928 }
1929 }
1930 if (ret == -1) {
1931 throw new NumberFormatException();
1932 }
1933 return ret;
1934 }
1935
1936 private void condenseEnergyMap(Map<Integer, Integer[]> energyMap) {
1937 Set<Integer> keys = energyMap.keySet();
1938 HashMap<Integer, Integer[]> tempMap = new HashMap<>();
1939 int count = 0;
1940 for (int key : keys) {
1941 tempMap.put(count, energyMap.get(key));
1942 count++;
1943 }
1944 energyMap.clear();
1945 energyMap.putAll(tempMap);
1946 }
1947 }