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;
39
40 import edu.rit.pj.Comm;
41 import edu.rit.pj.ParallelTeam;
42 import edu.rit.pj.WorkerTeam;
43 import ffx.algorithms.AlgorithmListener;
44 import ffx.algorithms.Terminatable;
45 import ffx.algorithms.mc.MCMove;
46 import ffx.algorithms.optimize.manybody.ManyBodyCell;
47 import ffx.algorithms.optimize.manybody.DistanceMatrix;
48 import ffx.algorithms.optimize.manybody.EliminatedRotamers;
49 import ffx.algorithms.optimize.manybody.EnergyExpansion;
50 import ffx.algorithms.optimize.manybody.EnergyRegion;
51 import ffx.algorithms.optimize.manybody.FourBodyEnergyRegion;
52 import ffx.algorithms.optimize.manybody.GoldsteinPairRegion;
53 import ffx.algorithms.optimize.manybody.RotamerMatrixMC;
54 import ffx.algorithms.optimize.manybody.RotamerMatrixMove;
55 import ffx.algorithms.optimize.manybody.SelfEnergyRegion;
56 import ffx.algorithms.optimize.manybody.ThreeBodyEnergyRegion;
57 import ffx.algorithms.optimize.manybody.TwoBodyEnergyRegion;
58 import ffx.crystal.Crystal;
59 import ffx.crystal.SymOp;
60 import ffx.numerics.Potential;
61 import ffx.potential.DualTopologyEnergy;
62 import ffx.potential.ForceFieldEnergy;
63 import ffx.potential.MolecularAssembly;
64 import ffx.potential.QuadTopologyEnergy;
65 import ffx.potential.bonded.Atom;
66 import ffx.potential.bonded.Polymer;
67 import ffx.potential.bonded.Residue;
68 import ffx.potential.bonded.ResidueState;
69 import ffx.potential.bonded.Rotamer;
70 import ffx.potential.bonded.RotamerLibrary;
71 import ffx.potential.nonbonded.NonbondedCutoff;
72 import ffx.potential.nonbonded.VanDerWaals;
73 import ffx.potential.openmm.OpenMMEnergy;
74 import ffx.potential.parsers.PDBFilter;
75 import ffx.utilities.ObjectPair;
76 import ffx.utilities.Resources;
77 import org.apache.commons.configuration2.CompositeConfiguration;
78 import org.apache.commons.io.FilenameUtils;
79
80 import java.io.BufferedWriter;
81 import java.io.File;
82 import java.io.FileWriter;
83 import java.io.IOException;
84 import java.nio.file.Path;
85 import java.nio.file.Paths;
86 import java.util.ArrayList;
87 import java.util.Arrays;
88 import java.util.Collections;
89 import java.util.Comparator;
90 import java.util.HashSet;
91 import java.util.List;
92 import java.util.Objects;
93 import java.util.Random;
94 import java.util.Set;
95 import java.util.function.BiFunction;
96 import java.util.function.ToDoubleFunction;
97 import java.util.logging.Level;
98 import java.util.logging.Logger;
99 import java.util.stream.Collectors;
100 import java.util.stream.IntStream;
101
102 import static ffx.potential.bonded.Residue.ResidueType.NA;
103 import static ffx.potential.bonded.RotamerLibrary.applyRotamer;
104 import static java.lang.Boolean.parseBoolean;
105 import static java.lang.Double.parseDouble;
106 import static java.lang.Integer.parseInt;
107 import static java.lang.String.format;
108 import static java.lang.System.arraycopy;
109 import static java.util.Arrays.copyOf;
110 import static org.apache.commons.math3.util.FastMath.abs;
111
112
113
114
115
116
117
118
119
120
121 public class RotamerOptimization implements Terminatable {
122
123
124
125
126 private static final Logger logger = Logger.getLogger(RotamerOptimization.class.getName());
127
128
129
130 private static final double FALLBACK_TWO_BODY_CUTOFF = 0;
131
132
133
134 protected final MolecularAssembly molecularAssembly;
135
136
137
138 protected final Potential potential;
139
140
141
142 protected final AlgorithmListener algorithmListener;
143
144
145
146 private final Comm world;
147
148
149
150 private final int numProc;
151
152
153
154 private final int rank;
155
156
157
158 private final boolean rank0;
159
160
161
162 private final boolean print = false;
163
164
165
166 private final boolean verboseEnergies = true;
167
168
169
170 private final boolean writeEnergyRestart = true;
171
172
173
174 private final BoxOptimization boxOpt;
175
176
177
178
179
180
181 private final BiFunction<List<Residue>, List<Rotamer>, File> dirSupplier;
182
183
184
185
186 private final ToDoubleFunction<File> eFunction;
187
188
189
190 private final boolean verbose;
191
192
193
194 private DistanceMatrix dM;
195
196
197
198 private EnergyExpansion eE;
199
200
201
202 private EliminatedRotamers eR;
203
204
205
206 protected RotamerLibrary library = RotamerLibrary.getDefaultLibrary();
207
208
209
210 private GoldsteinPairRegion goldsteinPairRegion;
211
212
213
214 private EnergyRegion energyRegion;
215
216
217
218 private boolean terminate = false;
219
220
221
222 private boolean done = true;
223
224
225
226 private double twoBodyCutoffDist;
227
228
229
230 private boolean threeBodyTerm = false;
231
232
233
234 private double threeBodyCutoffDist;
235
236
237
238 private int[][] resNeighbors;
239
240
241
242 private int[][] bidiResNeighbors;
243
244
245
246 private boolean pruneClashes = true;
247
248
249
250 private boolean prunePairClashes = true;
251
252
253
254 private int evaluatedPermutations = 0;
255
256
257
258 private int evaluatedPermutationsPrint = 0;
259
260
261
262 private List<Residue> residueList;
263
264
265
266 private int[] optimum;
267
268
269
270
271 private int windowSize = 7;
272
273
274
275 private int increment = 3;
276
277
278
279 private boolean revert;
280
281
282
283 private Direction direction = Direction.FORWARD;
284
285
286
287
288 private double distance = 2.0;
289
290
291
292 private DistanceMethod distanceMethod = DistanceMethod.RESIDUE;
293
294
295
296 private Algorithm algorithm = null;
297
298
299
300
301
302 private boolean useGoldstein = true;
303
304
305
306 private int ensembleNumber = 1;
307
308
309
310 private double ensembleBuffer = 0.0;
311
312
313
314 private double ensembleBufferStep = 0.5;
315
316
317
318 private double ensembleEnergy = 0.0;
319
320
321
322 private File ensembleFile;
323
324
325
326 private PDBFilter ensembleFilter;
327
328
329
330
331 private boolean lazyMatrix = false;
332
333
334
335 private double clashThreshold = 25.0;
336
337
338
339
340 private double multiResClashThreshold = 80.0;
341
342
343
344 private double pairClashThreshold = 25.0;
345
346
347
348 private double multiResPairClashAddn = 80.0;
349
350
351
352 private double[] x = null;
353
354
355
356
357 private boolean useForceFieldEnergy = false;
358
359
360
361
362 private double nucleicCorrectionThreshold = 0;
363
364
365
366 private double approximateEnergy = 0;
367
368
369
370
371 private int minNumberAcceptedNARotamers = 10;
372
373
374
375 private double nucleicPruningFactor = 10.0;
376
377
378
379 private double nucleicPairsPruningFactor = ((1.0 + nucleicPruningFactor) / 2);
380
381
382
383 private List<Residue> allResiduesList = null;
384
385
386
387 private Residue[] allResiduesArray = null;
388
389
390
391 private int nAllResidues = 0;
392
393
394
395
396 private int[] optimumSubset;
397
398
399
400 private boolean loadEnergyRestart = false;
401
402
403
404 private File energyRestartFile;
405
406
407
408 private ParallelTeam parallelTeam;
409
410
411
412
413 private boolean compute4BodyEnergy = false;
414
415
416
417 private boolean usingBoxOptimization = false;
418
419
420
421
422 private double superpositionThreshold = 0.25;
423
424
425
426 private boolean decomposeOriginal = false;
427
428
429
430 private boolean monteCarlo = false;
431
432
433
434 private int nMCSteps = 1000000;
435
436
437
438 private double mcTemp = 298.15;
439
440
441
442 private boolean mcUseAll = false;
443
444
445
446 private boolean mcNoEnum = false;
447
448
449
450
451 private boolean printFiles = true;
452
453
454
455 private List<ObjectPair<ResidueState[], Double>> ensembleStates;
456
457
458
459 private int maxRotCheckDepth;
460
461
462
463 private BufferedWriter energyWriter;
464
465
466
467
468 private boolean testing = false;
469
470
471
472 private boolean monteCarloTesting = false;
473
474
475
476 private boolean testSelfEnergyEliminations = false;
477
478
479
480
481
482 private int testPairEnergyEliminations = -1;
483
484
485
486
487
488 private int testTripleEnergyEliminations1 = -1;
489
490
491
492
493
494 private int testTripleEnergyEliminations2 = -1;
495
496
497
498 private boolean selfEliminationOn = true;
499
500
501
502 private boolean pairEliminationOn = true;
503
504
505
506
507
508
509
510
511 public RotamerOptimization(MolecularAssembly molecularAssembly, Potential potential,
512 AlgorithmListener algorithmListener) {
513
514 this.molecularAssembly = molecularAssembly;
515 this.potential = potential;
516 if (potential instanceof ForceFieldEnergy) {
517 ((ForceFieldEnergy) potential).setPrintOnFailure(false, true);
518 } else if (potential instanceof DualTopologyEnergy) {
519 ((DualTopologyEnergy) potential).setPrintOnFailure(false, true);
520 } else if (potential instanceof QuadTopologyEnergy) {
521 ((QuadTopologyEnergy) potential).setPrintOnFailure(false, true);
522 }
523 this.algorithmListener = algorithmListener;
524
525 eFunction = this::currentPE;
526 dirSupplier = (List<Residue> resList, List<Rotamer> rotList) -> null;
527
528 world = Comm.world();
529 numProc = world.size();
530 rank = world.rank();
531
532 rank0 = rank == 0;
533
534 boxOpt = new BoxOptimization();
535
536 CompositeConfiguration properties = molecularAssembly.getProperties();
537 verbose = properties.getBoolean("verbose", false);
538
539
540 ForceFieldEnergy forceFieldEnergy = molecularAssembly.getPotentialEnergy();
541 VanDerWaals vdW = forceFieldEnergy.getVdwNode();
542 if (vdW != null) {
543 NonbondedCutoff nonBondedCutoff = vdW.getNonbondedCutoff();
544 twoBodyCutoffDist = nonBondedCutoff.off;
545 } else {
546 twoBodyCutoffDist = FALLBACK_TWO_BODY_CUTOFF;
547 }
548
549
550
551 boolean testing = properties.getBoolean("manybody-testing", false);
552 if (testing) {
553 turnRotamerSingleEliminationOff();
554 turnRotamerPairEliminationOff();
555 }
556 boolean monteCarloTesting = properties.getBoolean("manybody-testing-mc", false);
557 if (monteCarloTesting) {
558 setMonteCarloTesting(true);
559 }
560
561
562 String computeQuads = properties.getString("ro-compute4BodyEnergy");
563 if (computeQuads != null) {
564 this.compute4BodyEnergy = parseBoolean(computeQuads);
565 logger.info(format(" (KEY) compute4BodyEnergy: %b", this.compute4BodyEnergy));
566 }
567
568
569 String direction = properties.getString("ro-direction");
570 String boxDimensions = properties.getString("ro-boxDimensions");
571 if (direction != null) {
572 this.direction = Direction.valueOf(direction);
573 logger.info(format(" (KEY) direction: %s", this.direction));
574 }
575 if (boxDimensions != null) {
576 boxOpt.update(boxDimensions);
577 }
578
579
580 String ensembleNumber = properties.getString("ro-ensembleNumber");
581 String ensembleEnergy = properties.getString("ro-ensembleEnergy");
582 String ensembleBuffer = properties.getString("ro-ensembleBuffer");
583 if (ensembleNumber != null) {
584 this.ensembleNumber = parseInt(ensembleNumber);
585 this.ensembleBuffer = 5.0;
586 this.ensembleBufferStep = 0.1 * this.ensembleBuffer;
587 logger.info(format(" (KEY) ensembleNumber: %d", this.ensembleNumber));
588 }
589 if (ensembleBuffer != null) {
590 this.ensembleBuffer = parseDouble(ensembleBuffer);
591 this.ensembleBufferStep = 0.1 * this.ensembleBuffer;
592 logger.info(format(" (KEY) ensembleBuffer: %.2f", this.ensembleBuffer));
593 }
594 if (ensembleEnergy != null) {
595 this.ensembleEnergy = parseDouble(ensembleEnergy);
596 logger.info(format(" (KEY) ensembleEnergy: %.2f", this.ensembleEnergy));
597 }
598
599
600 String nucleicPruningFactor = properties.getString("ro-nucleicPruningFactor");
601 String nucleicCorrectionThreshold = properties.getString("ro-nucleicCorrectionThreshold");
602 String minimumNumberAcceptedNARotamers = properties.getString(
603 "ro-minimumNumberAcceptedNARotamers");
604 if (nucleicPruningFactor != null) {
605 double value = parseDouble(nucleicPruningFactor);
606 this.nucleicPruningFactor = (value >= 0 ? value : 1.0);
607 this.nucleicPairsPruningFactor = (1.0 + value) / 2;
608 logger.info(format(" (KEY) nucleicPruningFactor: %.2f", this.nucleicPruningFactor));
609 }
610 if (nucleicCorrectionThreshold != null) {
611 double value = parseDouble(nucleicCorrectionThreshold);
612 this.nucleicCorrectionThreshold = (value >= 0 ? value : 0);
613 logger.info(
614 format(" (KEY) nucleicCorrectionThreshold: %.2f", this.nucleicCorrectionThreshold));
615 }
616 if (minimumNumberAcceptedNARotamers != null) {
617 int value = parseInt(minimumNumberAcceptedNARotamers);
618 this.minNumberAcceptedNARotamers = (value > 0 ? value : 10);
619 logger.info(
620 format(" (KEY) minimumNumberAcceptedNARotamers: %d", this.minNumberAcceptedNARotamers));
621 }
622
623
624 String superpositionThreshold = properties.getString("ro-superpositionThreshold");
625 if (superpositionThreshold != null) {
626 this.superpositionThreshold = parseDouble(superpositionThreshold);
627 logger.info(format(" (KEY) superpositionThreshold: %.2f", this.superpositionThreshold));
628 }
629
630
631 String multiResClashThreshold = properties.getString("ro-multiResClashThreshold");
632 String multiResPairClashAddition = properties.getString("ro-multiResPairClashAddition");
633 if (multiResClashThreshold != null) {
634 this.multiResClashThreshold = parseDouble(multiResClashThreshold);
635 logger.info(format(" (KEY) multiResClashThreshold: %.2f", this.multiResClashThreshold));
636 }
637 if (multiResPairClashAddition != null) {
638 this.multiResPairClashAddn = parseDouble(multiResPairClashAddition);
639 logger.info(format(" (KEY) multiResPairClashAddition: %.2f", this.multiResPairClashAddn));
640 }
641
642
643 String mcTemp = properties.getString("ro-mcTemp");
644 String mcUseAll = properties.getString("ro-mcUseAll");
645 String mcNoEnum = properties.getString("ro-debug-mcNoEnum");
646 if (mcTemp != null) {
647 this.mcTemp = parseDouble(mcTemp);
648 logIfRank0(format(" (KEY) mcTemp: %10.6f", this.mcTemp));
649 }
650 if (mcUseAll != null) {
651 this.mcUseAll = parseBoolean(mcUseAll);
652 logIfRank0(format(" (KEY) mcUseAll: %b", this.mcUseAll));
653 }
654 if (mcNoEnum != null) {
655 this.mcNoEnum = parseBoolean(mcNoEnum);
656 logIfRank0(format(" (KEY) debug-mcNoEnum: %b", this.mcNoEnum));
657 }
658
659 String lazyMatrix = properties.getString("ro-lazyMatrix");
660 if (lazyMatrix != null) {
661 this.lazyMatrix = parseBoolean(lazyMatrix);
662 logger.info(format(" (KEY) lazyMatrix: %b", lazyMatrix));
663 }
664
665 String propStr = properties.getString("ro-maxRotCheckDepth");
666 int defaultMaxRotCheckDepth = 1;
667 if (propStr != null) {
668 try {
669 maxRotCheckDepth = parseInt(propStr);
670 if (maxRotCheckDepth > 3 || maxRotCheckDepth < 0) {
671 throw new IllegalArgumentException(" ro-maxRotCheckDepth must be between 0-3 inclusive!");
672 }
673 } catch (Exception ex) {
674 maxRotCheckDepth = defaultMaxRotCheckDepth;
675 logger.warning(format(" Could not parse %s value %s as valid integer; defaulting to %d",
676 "ro-maxRotCheckDepth", propStr, maxRotCheckDepth));
677 logger.warning(format(" Exception: %s", ex));
678 }
679 } else {
680 maxRotCheckDepth = defaultMaxRotCheckDepth;
681 }
682
683 setUpRestart();
684 }
685
686
687
688
689
690
691 private static void turnOffAtoms(Residue residue) {
692 switch (residue.getResidueType()) {
693 case NA:
694 case AA:
695 List<Atom> atomList = residue.getVariableAtoms();
696 for (Atom atom : atomList) {
697 atom.setUse(false);
698 }
699 break;
700 default:
701 atomList = residue.getAtomList();
702 for (Atom atom : atomList) {
703 atom.setUse(false);
704 }
705 break;
706 }
707 }
708
709
710
711
712
713
714 private static void turnOnAtoms(Residue residue) {
715 switch (residue.getResidueType()) {
716 case NA:
717 case AA:
718 List<Atom> atomList = residue.getVariableAtoms();
719 for (Atom atom : atomList) {
720 atom.setUse(true);
721 }
722 break;
723 default:
724 atomList = residue.getAtomList();
725 for (Atom atom : atomList) {
726 atom.setUse(true);
727 }
728 break;
729 }
730 }
731
732
733
734
735
736
737
738
739
740 public boolean checkNeighboringPair(int i, int j) {
741 assert i != j;
742 final int first;
743 final int second;
744 if (i > j) {
745 first = j;
746 second = i;
747 } else {
748 first = i;
749 second = j;
750 }
751 return Arrays.stream(resNeighbors[first]).anyMatch(l -> l == second);
752 }
753
754
755
756
757
758
759
760
761
762
763
764 public boolean checkNeighboringTriple(int i, int j, int k) {
765 assert i != j && i != k && j != k;
766 int[] vals = new int[]{i, j, k};
767 Arrays.sort(vals);
768 final int first = vals[0];
769 final int second = vals[1];
770 final int third = vals[2];
771 if (!checkNeighboringPair(first, second)) {
772 return false;
773 }
774 return (checkNeighboringPair(first, third) || checkNeighboringPair(second, third));
775 }
776
777
778
779
780
781
782
783
784
785 public boolean checkValidMove(int i, int ri, int[] currentRots) {
786 int nRes = currentRots.length;
787 for (int j = 0; j < nRes; j++) {
788 if (j == i) {
789 continue;
790 }
791 if (eR.check(j, currentRots[j], i, ri)) {
792 return false;
793 }
794 }
795 return true;
796 }
797
798
799
800
801
802
803
804
805
806 public double computeBackboneEnergy(Residue[] residues) throws ArithmeticException {
807
808 Atom[] atoms = molecularAssembly.getAtomArray();
809 for (Atom atom : atoms) {
810 atom.setUse(true);
811 }
812
813
814 turnOffAllResidues(residues);
815
816
817 return currentEnergy(residues);
818 }
819
820
821
822
823
824
825
826
827
828
829 public double computeEnergy(Residue[] residues, int[] rotamers, boolean print) {
830 double selfSum;
831 double pairSum;
832 double threeBodySum;
833 try {
834 if (parallelTeam == null) {
835 parallelTeam = new ParallelTeam();
836 }
837 if (energyRegion == null) {
838 energyRegion = new EnergyRegion(parallelTeam.getThreadCount());
839 }
840 energyRegion.init(eE, residues, rotamers, threeBodyTerm);
841 parallelTeam.execute(energyRegion);
842 selfSum = energyRegion.getSelf();
843 pairSum = energyRegion.getTwoBody();
844 threeBodySum = energyRegion.getThreeBody();
845 } catch (Exception e) {
846 throw new IllegalArgumentException(e);
847 }
848
849 approximateEnergy = eE.getBackboneEnergy() + selfSum + pairSum + threeBodySum;
850
851 if (print) {
852 logger.info(format(" Backbone %s", formatEnergy(eE.getBackboneEnergy())));
853 logger.info(format(" Self Energy %s", formatEnergy(selfSum)));
854 logger.info(format(" Pair Energy %s", formatEnergy(pairSum)));
855 if (!threeBodyTerm) {
856 logger.info(format(" Total Energy up to 2-Body %s", formatEnergy(approximateEnergy)));
857 } else {
858 logger.info(format(" 3-Body Energy %s", formatEnergy(threeBodySum)));
859 logger.info(format(" Total Energy up to 3-Body %s", formatEnergy(approximateEnergy)));
860 }
861 }
862 return approximateEnergy;
863 }
864
865
866
867
868
869
870
871 public double currentEnergy(Residue[] resArray) throws ArithmeticException {
872 return currentEnergy(Arrays.asList(resArray));
873 }
874
875
876
877
878
879
880
881
882
883 public double currentEnergyWrapper(List<Residue> resList) throws ArithmeticException {
884 return currentEnergy(resList);
885 }
886
887
888
889
890
891
892 public double currentFFXPE() {
893 if (x == null) {
894 int n = potential.getNumberOfVariables();
895 x = new double[n];
896 }
897 potential.getCoordinates(x);
898 return ((OpenMMEnergy) potential).energyFFX(x, false);
899 }
900
901
902
903
904
905
906
907 public String formatEnergy(double energy) {
908 if (abs(energy) < 1.0e6) {
909 return format("%16.8f", energy);
910 } else {
911 return format("*%15.4e", energy);
912 }
913 }
914
915 public double getApproximate() {
916 return approximateEnergy;
917 }
918
919 public double getBackboneEnergy() {
920 return eE.getBackboneEnergy();
921 }
922
923 public EliminatedRotamers getEliminatedRotamers() {
924 return eR;
925 }
926
927 public EnergyExpansion getEnergyExpansion() {
928 return eE;
929 }
930
931
932
933
934
935
936 public List<ResidueState[]> getEnsemble() {
937 if (ensembleStates == null) {
938 return null;
939 } else {
940 List<ResidueState[]> states = new ArrayList<>(ensembleStates.size());
941 ensembleStates.forEach((es) -> states.add(es.val()));
942 return states;
943 }
944 }
945
946
947
948
949
950
951 public void setEnsemble(int ensemble) {
952 setEnsemble(ensemble, 5.0);
953 }
954
955
956
957
958
959
960 public List<Residue> getResidues() {
961 return residueList;
962 }
963
964
965
966
967
968
969 public void setResidues(List<Residue> residueList) {
970 this.residueList = residueList;
971 }
972
973
974
975
976
977
978 public File getRestartFile() {
979 return energyRestartFile;
980 }
981
982 public double goldsteinPairSumOverK(Residue[] residues, int lb, int ub, int i, int riA, int riB,
983 int j, int rjC, int rjD, List<Residue> blockedResidues, int[] possK) {
984 double sumOverK = 0.0;
985
986 for (int indK = lb; indK <= ub; indK++) {
987 int k = possK[indK];
988 double minForResK = Double.MAX_VALUE;
989 Residue resk = residues[k];
990 Rotamer[] rotk = resk.getRotamers();
991 int nrk = rotk.length;
992 int rkEvals = 0;
993
994 for (int rk = 0; rk < nrk; rk++) {
995 if (eR.check(k, rk)) {
996 continue;
997 }
998
999 if (eR.check(i, riA, k, rk) || eR.check(j, rjC, k, rk)) {
1000
1001 continue;
1002 }
1003
1004 if (eR.check(i, riB, k, rk) || eR.check(j, rjD, k, rk)) {
1005 blockedResidues.add(resk);
1006 return Double.NaN;
1007 }
1008
1009 rkEvals++;
1010 double currentResK =
1011 eE.get2Body(i, riA, k, rk) - eE.get2Body(i, riB, k, rk) + eE.get2Body(j, rjC, k, rk)
1012 - eE.get2Body(j, rjD, k, rk);
1013
1014 if (threeBodyTerm) {
1015 double sumOverL = (eE.get3Body(residues, i, riA, j, rjC, k, rk) - eE.get3Body(residues, i,
1016 riB, j, rjD, k, rk));
1017
1018 int[] nK = bidiResNeighbors[k];
1019 IntStream lStream = IntStream.concat(Arrays.stream(possK), Arrays.stream(nK));
1020 int[] possL = lStream.filter(l -> (l != i && l != j && l != k)).sorted().distinct()
1021 .toArray();
1022
1023 for (int l : possL) {
1024 if (l == k || l == i || l == j) {
1025 continue;
1026 }
1027 Residue residuel = residues[l];
1028 Rotamer[] rotamersl = residuel.getRotamers();
1029 int nrl = rotamersl.length;
1030 int rlEvaluations = 0;
1031 double minForResL = Double.MAX_VALUE;
1032
1033 for (int rl = 0; rl < nrl; rl++) {
1034
1035 if (eR.check(l, rl) || eR.check(k, rk, l, rl) || eR.check(i, riA, l, rl) || eR.check(j,
1036 rjC, l, rl)) {
1037
1038
1039 continue;
1040 }
1041 if (eR.check(i, riB, l, rl) || eR.check(j, rjD, l, rl)) {
1042
1043
1044 blockedResidues.add(residuel);
1045 return Double.NaN;
1046 }
1047 rlEvaluations++;
1048 double e =
1049 eE.get3Body(residues, i, riA, k, rk, l, rl) - eE.get3Body(residues, i, riB, k, rk,
1050 l, rl) + eE.get3Body(residues, j, rjC, k, rk, l, rl) - eE.get3Body(residues, j,
1051 rjD, k, rk, l, rl);
1052 if (e < minForResL) {
1053 minForResL = e;
1054 }
1055 }
1056 if (rlEvaluations == 0) {
1057 minForResL = 0.0;
1058 }
1059 sumOverL += minForResL;
1060 }
1061 currentResK += sumOverL;
1062 }
1063 if (currentResK < minForResK) {
1064 minForResK = currentResK;
1065 }
1066 }
1067 if (rkEvals == 0) {
1068 minForResK = 0.0;
1069 blockedResidues.add(resk);
1070 }
1071 sumOverK += minForResK;
1072 }
1073 return sumOverK;
1074 }
1075
1076 public void logIfRank0(String msg) {
1077 if (rank0) {
1078 logger.info(msg);
1079 }
1080 }
1081
1082 public void logIfRank0(String msg, Level level) {
1083 if (rank0) {
1084 logger.log(level, msg);
1085 }
1086 }
1087
1088
1089
1090
1091
1092
1093
1094 public double optimize(Algorithm algorithm) {
1095 this.algorithm = algorithm;
1096 return optimize();
1097 }
1098
1099
1100
1101
1102
1103
1104 public double optimize() {
1105 double e = 0;
1106 try {
1107
1108 CompositeConfiguration properties = molecularAssembly.getProperties();
1109 boolean ignoreNA = properties.getBoolean("ignoreNA", false);
1110 if (ignoreNA) {
1111 logger.info(" Ignoring nucleic acids.");
1112 }
1113
1114 logger.info(format("\n Rotamer Library: %s", library.getLibrary()));
1115 logger.info(format(" Algorithm: %s", algorithm));
1116 logger.info(format(" Goldstein Criteria: %b", useGoldstein));
1117 logger.info(format(" Three-Body Energies: %b\n", threeBodyTerm));
1118
1119
1120
1121
1122 allResiduesList = new ArrayList<>();
1123
1124 Polymer[] polymers = molecularAssembly.getChains();
1125 for (Polymer polymer : polymers) {
1126 List<Residue> current = polymer.getResidues();
1127 for (Residue residuej : current) {
1128 residuej.setRotamers(library);
1129 if (residuej.getRotamers() != null) {
1130 if (!(ignoreNA && residuej.getResidueType() == Residue.ResidueType.NA)) {
1131 allResiduesList.add(residuej);
1132 }
1133 }
1134 }
1135 }
1136
1137 sortResidues(allResiduesList);
1138 sortResidues(residueList);
1139
1140
1141 if (ignoreNA) {
1142 List<Residue> onlyAA = new ArrayList<>();
1143 for (Residue res : residueList) {
1144 if (res.getResidueType() != Residue.ResidueType.NA) {
1145 onlyAA.add(res);
1146 }
1147 }
1148 residueList = onlyAA;
1149 }
1150
1151 RotamerLibrary.initializeDefaultAtomicCoordinates(
1152 molecularAssembly.getChains());
1153
1154 nAllResidues = allResiduesList.size();
1155 allResiduesArray = allResiduesList.toArray(new Residue[nAllResidues]);
1156
1157
1158
1159
1160
1161
1162
1163 if (distance > 0) {
1164 dM = new DistanceMatrix(molecularAssembly, algorithmListener, allResiduesArray,
1165 allResiduesList, distanceMethod, distance, twoBodyCutoffDist, threeBodyCutoffDist,
1166 lazyMatrix);
1167 }
1168
1169 if (residueList != null) {
1170
1171
1172 optimum = new int[residueList.size()];
1173
1174 done = false;
1175 terminate = false;
1176
1177 switch (algorithm) {
1178 case INDEPENDENT -> e = independent(residueList);
1179 case BRUTE_FORCE -> e = bruteForce(residueList);
1180 case ALL -> {
1181 e = globalOptimization(residueList);
1182 arraycopy(optimumSubset, 0, optimum, 0, residueList.size());
1183 }
1184 case WINDOW -> e = slidingWindowOptimization(residueList, windowSize, increment, revert, distance,
1185 direction);
1186 case BOX -> e = boxOpt.boxOptimization(residueList);
1187 default -> {
1188 }
1189 }
1190 terminate = false;
1191 done = true;
1192 }
1193 } finally {
1194 try {
1195 if (energyWriter != null) {
1196 energyWriter.close();
1197 }
1198 } catch (IOException ex) {
1199 logger.severe(" Exception in closing buffered energy writer.");
1200 }
1201 }
1202
1203 return e;
1204 }
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216 public double rotamerOptimization(MolecularAssembly molecularAssembly, Residue[] residues, int i,
1217 double lowEnergy, int[] optimum) {
1218
1219
1220 if (i == 0) {
1221 evaluatedPermutations = 0;
1222 }
1223
1224 int nResidues = residues.length;
1225 Residue current = residues[i];
1226 Rotamer[] rotamers = current.getRotamers();
1227 int lenri = rotamers.length;
1228 double currentEnergy = Double.MAX_VALUE;
1229 List<Residue> resList = Arrays.asList(residues);
1230 if (i < nResidues - 1) {
1231
1232
1233 int minRot = -1;
1234 for (int ri = 0; ri < lenri; ri++) {
1235 applyRotamer(current, rotamers[ri]);
1236 double rotEnergy = rotamerOptimization(molecularAssembly, residues, i + 1, lowEnergy,
1237 optimum);
1238 if (rotEnergy < currentEnergy) {
1239 currentEnergy = rotEnergy;
1240 }
1241 if (rotEnergy < lowEnergy) {
1242 minRot = ri;
1243 lowEnergy = rotEnergy;
1244 }
1245 }
1246 if (minRot > -1) {
1247 optimum[i] = minRot;
1248 }
1249 } else {
1250
1251
1252
1253
1254
1255 int[] rotArray = copyOf(optimum, nResidues);
1256 for (int ri = 0; ri < lenri; ri++) {
1257 applyRotamer(current, rotamers[ri]);
1258 double rotEnergy = Double.NaN;
1259 try {
1260 rotArray[nResidues - 1] = ri;
1261 rotEnergy = currentEnergy(resList) + eE.getTotalRotamerPhBias(resList, rotArray);
1262 logger.info(format(" %d Energy: %s", ++evaluatedPermutations, formatEnergy(rotEnergy)));
1263 } catch (ArithmeticException ex) {
1264 logger.info(
1265 format(" %d Energy set to NaN (unreasonable conformation)", ++evaluatedPermutations));
1266 }
1267 if (algorithmListener != null) {
1268 algorithmListener.algorithmUpdate(molecularAssembly);
1269 }
1270 if (rotEnergy < currentEnergy) {
1271 currentEnergy = rotEnergy;
1272 }
1273 if (rotEnergy < lowEnergy) {
1274 lowEnergy = rotEnergy;
1275 optimum[nResidues - 1] = ri;
1276 }
1277 }
1278 }
1279 return currentEnergy;
1280 }
1281
1282
1283
1284
1285
1286
1287
1288 public void setApproxBoxLength(double approxBoxLength) {
1289 boxOpt.approxBoxLength = approxBoxLength;
1290 }
1291
1292
1293
1294
1295
1296
1297 public void setBoxBorderSize(double boxBorderSize) {
1298 boxOpt.cellBorderSize = boxBorderSize;
1299 }
1300
1301
1302
1303
1304
1305
1306 public void setBoxEnd(int boxEnd) {
1307
1308 boxOpt.cellEnd = boxEnd;
1309 }
1310
1311
1312
1313
1314
1315
1316
1317 public void setBoxInclusionCriterion(int boxInclusionCriterion) {
1318 boxOpt.boxInclusionCriterion = boxInclusionCriterion;
1319 }
1320
1321
1322
1323
1324
1325
1326 public void setBoxStart(int boxStart) {
1327 boxOpt.cellStart = boxStart;
1328 }
1329
1330
1331
1332
1333
1334
1335 public void setCoordinatesToEnsemble(int ensnum) {
1336 if (ensembleStates != null && !ensembleStates.isEmpty()) {
1337 ensnum %= ensembleStates.size();
1338 ResidueState.revertAllCoordinates(residueList, ensembleStates.get(ensnum).val());
1339 } else {
1340 throw new IllegalArgumentException(" Ensemble states not initialized!");
1341 }
1342 }
1343
1344
1345
1346
1347
1348
1349 public void setDecomposeOriginal(boolean decomposeOriginal) {
1350 this.decomposeOriginal = decomposeOriginal;
1351 }
1352
1353
1354
1355
1356
1357
1358 public void setDirection(Direction direction) {
1359 this.direction = direction;
1360 }
1361
1362
1363
1364
1365
1366
1367 public void setDistanceCutoff(double distance) {
1368 this.distance = distance;
1369 }
1370
1371
1372
1373
1374
1375
1376 public void setEnergyRestartFile(File file) {
1377 loadEnergyRestart = true;
1378 energyRestartFile = file;
1379 }
1380
1381
1382
1383
1384
1385
1386
1387 public void setEnsemble(int ensemble, double ensembleBuffer) {
1388 this.ensembleNumber = ensemble;
1389 this.ensembleBuffer = ensembleBuffer;
1390 this.ensembleBufferStep = 0.1 * ensembleBuffer;
1391 if (ensemble > 1) {
1392 setPruning(0);
1393 }
1394 }
1395
1396
1397
1398
1399
1400
1401 public void setIncrement(int increment) {
1402 this.increment = increment;
1403 }
1404
1405
1406
1407
1408
1409
1410 public void setMaxRotCheckDepth(int maxRotCheckDepth) {
1411 this.maxRotCheckDepth = maxRotCheckDepth;
1412 }
1413
1414
1415
1416
1417
1418
1419 public void setMinimumNumberAcceptedNARotamers(int minNumberAcceptedNARotamers) {
1420 if (minNumberAcceptedNARotamers > 0) {
1421 this.minNumberAcceptedNARotamers = minNumberAcceptedNARotamers;
1422 } else {
1423 logger.warning(
1424 "\n Minimum number of accepted NA rotamers must be a positive integer.\n Setting to default value 10.\n");
1425 this.minNumberAcceptedNARotamers = 10;
1426 }
1427 }
1428
1429
1430
1431
1432
1433
1434
1435 public void setMonteCarlo(boolean monteCarlo, int nMCsteps) {
1436 this.monteCarlo = monteCarlo;
1437 this.nMCSteps = nMCsteps;
1438 }
1439
1440
1441
1442
1443
1444
1445
1446
1447
1448 public void setMonteCarloTesting(boolean bool) {
1449 this.monteCarloTesting = bool;
1450 }
1451
1452
1453
1454
1455
1456
1457 public void setNucleicCorrectionThreshold(double nucleicCorrectionThreshold) {
1458 if (nucleicCorrectionThreshold >= 0) {
1459 this.nucleicCorrectionThreshold = nucleicCorrectionThreshold;
1460 } else {
1461 logger.warning(
1462 "\n Correction threshold must be >= 0. Setting to default of 0 (threshold inactive).\n");
1463 this.nucleicCorrectionThreshold = 0;
1464 }
1465 }
1466
1467
1468
1469
1470
1471
1472 public void setNucleicPruningFactor(double nucleicPruningFactor) {
1473 this.nucleicPruningFactor = nucleicPruningFactor;
1474 this.nucleicPairsPruningFactor = ((1.0 + nucleicPruningFactor) / 2);
1475 }
1476
1477
1478
1479
1480
1481
1482 public void setNumXYZBoxes(int[] numXYZBoxes) {
1483 arraycopy(numXYZBoxes, 0, boxOpt.numXYZCells, 0, boxOpt.numXYZCells.length);
1484 }
1485
1486
1487
1488
1489
1490
1491 public void setPairClashThreshold(double pairClashThreshold) {
1492 this.pairClashThreshold = pairClashThreshold;
1493 }
1494
1495
1496
1497
1498
1499
1500
1501 public void setPrintFiles(boolean printFiles) {
1502 this.printFiles = printFiles;
1503 }
1504
1505
1506
1507
1508
1509
1510 public void setPruning(int set) {
1511 if (set == 0) {
1512 this.pruneClashes = false;
1513 this.prunePairClashes = false;
1514 } else if (set == 1) {
1515 this.pruneClashes = true;
1516 this.prunePairClashes = false;
1517 } else if (set == 2) {
1518 this.pruneClashes = true;
1519 this.prunePairClashes = true;
1520 }
1521 }
1522
1523
1524
1525
1526
1527
1528 public void setResiduesIgnoreNull(List<Residue> residues) {
1529 residueList = new ArrayList<>();
1530 logger.fine(" Optimizing these residues: ");
1531 for (Residue r : residues) {
1532 if (r.getRotamers() != null) {
1533 residueList.add(r);
1534 logger.fine(format("\t%s", r));
1535 } else {
1536 logger.fine(format(" not \t%s", r));
1537 }
1538 }
1539 }
1540
1541
1542
1543
1544
1545
1546 public void setRevert(boolean revert) {
1547 this.revert = revert;
1548 }
1549
1550 public void setRotamerLibrary(RotamerLibrary lib) {
1551 library = lib;
1552 }
1553
1554
1555
1556
1557
1558
1559 public void setSingletonClashThreshold(double singletonClashThreshold) {
1560 this.clashThreshold = singletonClashThreshold;
1561 }
1562
1563
1564
1565
1566
1567
1568 public void setSuperpositionThreshold(double superpositionThreshold) {
1569 this.superpositionThreshold = superpositionThreshold;
1570 }
1571
1572
1573
1574
1575
1576
1577
1578
1579 public void setThreeBodyCutoff(double threeBodyCutoffDist) {
1580 this.threeBodyCutoffDist = threeBodyCutoffDist;
1581 if (this.threeBodyCutoffDist < 0) {
1582 logger.info("Warning: threeBodyCutoffDist should not be less than 0.");
1583 }
1584 }
1585
1586
1587
1588
1589
1590
1591
1592 public void setThreeBodyEnergy(boolean threeBodyTerm) {
1593 this.threeBodyTerm = threeBodyTerm;
1594 }
1595
1596
1597
1598
1599
1600
1601
1602
1603 public void setTwoBodyCutoff(double twoBodyCutoffDist) {
1604 this.twoBodyCutoffDist = twoBodyCutoffDist;
1605 if (this.twoBodyCutoffDist < 0) {
1606 logger.info("Warning: threeBodyCutoffDist should not be less than 0.");
1607 }
1608 }
1609
1610
1611
1612
1613
1614
1615 public void setUseGoldstein(boolean useGoldstein) {
1616 this.useGoldstein = useGoldstein;
1617 }
1618
1619
1620
1621
1622
1623
1624 public void setWindowSize(int windowSize) {
1625 this.windowSize = windowSize;
1626 if (this.increment > windowSize) {
1627 logger.info(format(" Decreasing increment to match window size %d", windowSize));
1628 this.increment = windowSize;
1629 }
1630 }
1631
1632
1633
1634
1635 @Override
1636 public void terminate() {
1637 terminate = true;
1638 while (!done) {
1639 synchronized (this) {
1640 try {
1641 wait(1);
1642 } catch (InterruptedException e) {
1643 logger.log(Level.WARNING, " Exception terminating rotamer optimization.\n", e);
1644 }
1645 }
1646 }
1647 }
1648
1649
1650
1651
1652 @Override
1653 public String toString() {
1654 if (eR != null) {
1655 return eR.toString();
1656 } else {
1657 return null;
1658 }
1659 }
1660
1661 public void turnOffAllResidues(Residue[] residues) {
1662 if (residues == null) {
1663 return;
1664 }
1665 for (Residue residue : residues) {
1666 turnOffResidue(residue);
1667 }
1668 }
1669
1670 public void turnOffResidue(Residue residue) {
1671 turnOffAtoms(residue);
1672 applyDefaultRotamer(residue);
1673 }
1674
1675 public void turnOnAllResidues(Residue[] residues) {
1676 if (residues == null) {
1677 return;
1678 }
1679 for (Residue residue : residues) {
1680 turnOnAtoms(residue);
1681 }
1682 }
1683
1684 public void turnOnResidue(Residue residue, int ri) {
1685 turnOnAtoms(residue);
1686 Rotamer[] rotamers = residue.getRotamers();
1687 applyRotamer(residue, rotamers[ri]);
1688 }
1689
1690
1691
1692
1693 public void turnRotamerPairEliminationOff() {
1694 logger.info(" Turning off pair eliminations.");
1695 pairEliminationOn = false;
1696 }
1697
1698
1699
1700
1701 public void turnRotamerSingleEliminationOff() {
1702 logger.info(" Turning off single eliminations.");
1703 selfEliminationOn = false;
1704 }
1705
1706
1707
1708
1709
1710
1711
1712
1713
1714
1715
1716
1717 private double decomposedRotamerOptimization(Residue[] residues, int i, double lowEnergy,
1718 int[] optimum, int[] currentRotamers) {
1719
1720
1721 if (i == 0) {
1722 evaluatedPermutations = 0;
1723 }
1724
1725 int nResidues = residues.length;
1726 Residue current = residues[i];
1727 Rotamer[] rotamers = current.getRotamers();
1728 int lenri = rotamers.length;
1729 double currentEnergy = Double.MAX_VALUE;
1730 if (i < nResidues - 1) {
1731
1732
1733 for (int ri = 0; ri < lenri; ri++) {
1734 if (eR.check(i, ri)) {
1735 continue;
1736 }
1737 currentRotamers[i] = ri;
1738 double rotEnergy = decomposedRotamerOptimization(residues, i + 1, lowEnergy, optimum,
1739 currentRotamers);
1740 if (rotEnergy < lowEnergy) {
1741 lowEnergy = rotEnergy;
1742 }
1743 if (rotEnergy < currentEnergy) {
1744 currentEnergy = rotEnergy;
1745 }
1746 }
1747 } else {
1748
1749
1750 for (int ri = 0; ri < lenri; ri++) {
1751 if (eR.check(i, ri)) {
1752 continue;
1753 }
1754 currentRotamers[i] = ri;
1755 double rotEnergy = computeEnergy(residues, currentRotamers, false);
1756 ++evaluatedPermutations;
1757 if (rotEnergy < currentEnergy) {
1758 currentEnergy = rotEnergy;
1759 }
1760 if (rotEnergy < lowEnergy) {
1761
1762
1763
1764
1765
1766 arraycopy(currentRotamers, 0, optimum, 0, optimum.length);
1767 if (evaluatedPermutations > 1) {
1768 logger.info(
1769 format(" Minimum energy update: %f < %f, permutation %d", rotEnergy, lowEnergy,
1770 evaluatedPermutations));
1771 String permutation = " Rotamer permutation: " + optimum[0];
1772 for (int j = 1; j < nResidues; j++) {
1773 permutation = permutation.concat(", " + optimum[j]);
1774 }
1775 logger.info(permutation);
1776 } else {
1777 logger.info(format(" First minimum energy (permutation 1): %f", rotEnergy));
1778 }
1779 lowEnergy = rotEnergy;
1780 }
1781 }
1782 }
1783 return currentEnergy;
1784 }
1785
1786
1787
1788
1789
1790
1791
1792
1793
1794
1795
1796
1797
1798
1799 private double rotamerOptimizationMC(Residue[] residues, int[] optimum, int[] initialRots,
1800 int maxIters, boolean randomizeRots, boolean useAllElims) {
1801
1802 long initTime = -System.nanoTime();
1803 if (randomizeRots) {
1804 randomizeRotamers(initialRots, residues, true);
1805 }
1806
1807 int nRes = residues.length;
1808 arraycopy(initialRots, 0, optimum, 0, nRes);
1809 assert optimum.length == nRes;
1810 assert initialRots.length == nRes;
1811
1812 RotamerMatrixMC rmc = new RotamerMatrixMC(initialRots, residues, useForceFieldEnergy, this);
1813 rmc.setTemperature(mcTemp);
1814 RotamerMatrixMove rmove = new RotamerMatrixMove(useAllElims, initialRots, residues, this, eR,
1815 monteCarloTesting);
1816 List<MCMove> rmList = new ArrayList<>(1);
1817 rmList.add(rmove);
1818
1819 double initialEnergy = computeEnergy(residues, initialRots, false);
1820 double optimumEnergy = initialEnergy;
1821 double currentEnergy = initialEnergy;
1822
1823 int nAccept = 0;
1824 logIfRank0(
1825 format(" Beginning %d iterations of Monte Carlo search " + "starting from energy %10.6f",
1826 maxIters, initialEnergy));
1827
1828 for (int i = 0; i < maxIters; i++) {
1829 if (rmc.mcStep(rmList, currentEnergy)) {
1830 currentEnergy = rmc.lastEnergy();
1831 ++nAccept;
1832 if (currentEnergy < optimumEnergy) {
1833 optimumEnergy = currentEnergy;
1834 arraycopy(initialRots, 0, optimum, 0, nRes);
1835 }
1836 } else {
1837 currentEnergy = rmc.lastEnergy();
1838 }
1839 }
1840
1841 initTime += System.nanoTime();
1842 double fractAccept = ((double) nAccept) / ((double) maxIters);
1843 logIfRank0(
1844 format(" %d steps of DEE-MC completed in %10.6f seconds", maxIters, (initTime * 1.0E-9)));
1845 logIfRank0(format(" Number of steps accepted: %d for %10.6f of total", nAccept, fractAccept));
1846 logIfRank0(format(" Lowest energy found: %10.6f kcal/mol", optimumEnergy));
1847 logIfRank0(format(" Final energy found: %10.6f kcal/mol", currentEnergy));
1848
1849 return optimumEnergy;
1850 }
1851
1852
1853
1854
1855
1856
1857
1858
1859 private void randomizeRotamers(int[] rotamers, Residue[] residues, boolean useAllElims) {
1860 int nRes = rotamers.length;
1861 for (int i = 0; i < nRes; i++) {
1862 Rotamer[] rotsi = residues[i].getRotamers();
1863 int lenri = rotsi.length;
1864 List<Integer> allowedRots = new ArrayList<>(lenri);
1865
1866 for (int ri = 0; ri < lenri; ri++) {
1867 if (!eR.check(i, ri)) {
1868 allowedRots.add(ri);
1869 }
1870 }
1871
1872 Random rand = new Random();
1873 int nRots = allowedRots.size();
1874 if (monteCarloTesting) {
1875 rand.setSeed(nRots);
1876 }
1877 if (nRots > 1) {
1878 boolean validMove = !useAllElims;
1879 int indexRI;
1880 do {
1881 int ri = rand.nextInt(nRots);
1882 indexRI = allowedRots.get(ri);
1883 if (useAllElims) {
1884 validMove = checkValidMove(i, indexRI, rotamers);
1885 }
1886 } while (!validMove);
1887 rotamers[i] = indexRI;
1888 }
1889 }
1890 }
1891
1892
1893
1894
1895
1896
1897
1898
1899
1900
1901
1902
1903
1904
1905 private double rotamerOptimizationDEE(MolecularAssembly molecularAssembly, Residue[] residues,
1906 int i, int[] currentRotamers, double lowEnergy, int[] optimum, double[] permutationEnergies) {
1907
1908 if (i == 0) {
1909 evaluatedPermutations = 0;
1910 }
1911
1912 int nResidues = residues.length;
1913 Residue residuei = residues[i];
1914 Rotamer[] rotamersi = residuei.getRotamers();
1915 int lenri = rotamersi.length;
1916 double currentEnergy = Double.MAX_VALUE;
1917 List<Residue> resList = Arrays.asList(residues);
1918
1919
1920
1921 if (i < nResidues - 1) {
1922
1923 for (int ri = 0; ri < lenri; ri++) {
1924
1925 if (eR.check(i, ri)) {
1926 continue;
1927 }
1928
1929
1930 boolean deadEnd = false;
1931 for (int j = 0; j < i; j++) {
1932 int rj = currentRotamers[j];
1933 deadEnd = eR.check(j, rj, i, ri);
1934 if (deadEnd) {
1935 break;
1936 }
1937 }
1938 if (deadEnd) {
1939 continue;
1940 }
1941 applyRotamer(residuei, rotamersi[ri]);
1942 currentRotamers[i] = ri;
1943 double rotEnergy = rotamerOptimizationDEE(molecularAssembly, residues, i + 1,
1944 currentRotamers, lowEnergy, optimum, permutationEnergies);
1945 if (rotEnergy < currentEnergy) {
1946 currentEnergy = rotEnergy;
1947 }
1948 if (rotEnergy < lowEnergy) {
1949 optimum[i] = ri;
1950 lowEnergy = rotEnergy;
1951 }
1952 }
1953 } else {
1954 if (ensembleStates == null) {
1955 ensembleStates = new ArrayList<>();
1956 }
1957
1958
1959
1960
1961
1962
1963 for (int ri = 0; ri < lenri; ri++) {
1964
1965 if (eR.check(i, ri)) {
1966 continue;
1967 }
1968 currentRotamers[i] = ri;
1969
1970
1971 boolean deadEnd = false;
1972 for (int j = 0; j < i; j++) {
1973 int rj = currentRotamers[j];
1974 deadEnd = eR.check(j, rj, i, ri);
1975 if (deadEnd) {
1976 break;
1977 }
1978 }
1979 if (deadEnd) {
1980 continue;
1981 }
1982 applyRotamer(residuei, rotamersi[ri]);
1983
1984 approximateEnergy = computeEnergy(residues, currentRotamers, false);
1985 double comparisonEnergy = approximateEnergy;
1986 evaluatedPermutations++;
1987
1988 if (useForceFieldEnergy) {
1989 double amoebaEnergy = Double.NaN;
1990 try {
1991
1992 amoebaEnergy =
1993 currentEnergy(resList) + eE.getTotalRotamerPhBias(resList, currentRotamers);
1994 } catch (ArithmeticException ex) {
1995 logger.warning(
1996 format(" Exception %s in calculating full AMOEBA energy for permutation %d", ex,
1997 evaluatedPermutations));
1998 }
1999 comparisonEnergy = amoebaEnergy;
2000 }
2001 if (permutationEnergies != null) {
2002 permutationEnergies[evaluatedPermutations - 1] = comparisonEnergy;
2003 }
2004 if (algorithmListener != null) {
2005 algorithmListener.algorithmUpdate(molecularAssembly);
2006 }
2007 if (ensembleNumber > 1) {
2008 if (rank0 && printFiles) {
2009 try {
2010 FileWriter fw = new FileWriter(ensembleFile, true);
2011 BufferedWriter bw = new BufferedWriter(fw);
2012 bw.write(format("MODEL %d", evaluatedPermutations));
2013 for (int j = 0; j < 75; j++) {
2014 bw.write(" ");
2015 }
2016 bw.newLine();
2017 bw.flush();
2018 ensembleFilter.writeFile(ensembleFile, true);
2019 bw.write("ENDMDL");
2020 for (int j = 0; j < 64; j++) {
2021 bw.write(" ");
2022 }
2023 bw.newLine();
2024 bw.close();
2025 } catch (IOException e) {
2026 logger.warning(format("Exception writing to file: %s", ensembleFile.getName()));
2027 }
2028 }
2029 ResidueState[] states = ResidueState.storeAllCoordinates(residues);
2030 ensembleStates.add(new ObjectPair<>(states, comparisonEnergy));
2031 }
2032
2033 if (comparisonEnergy < currentEnergy) {
2034 currentEnergy = comparisonEnergy;
2035 }
2036
2037 if (comparisonEnergy < lowEnergy) {
2038 lowEnergy = comparisonEnergy;
2039 optimum[i] = ri;
2040 }
2041
2042 if (useForceFieldEnergy) {
2043
2044 logIfRank0(format(" %6e AMOEBA: %12.4f 3-Body: %12.4f Neglected: %12.4f (%12.4f)",
2045 (double) evaluatedPermutations, comparisonEnergy, approximateEnergy,
2046 comparisonEnergy - approximateEnergy, lowEnergy));
2047 } else {
2048 logIfRank0(
2049 format(" %12s %25f %25f", evaluatedPermutations, approximateEnergy, lowEnergy));
2050 }
2051 }
2052
2053 ensembleStates.sort(null);
2054 }
2055 return currentEnergy;
2056 }
2057
2058
2059
2060
2061
2062
2063
2064
2065
2066
2067 private double dryRun(Residue[] residues, int i, int[] currentRotamers) {
2068
2069 if (i == 0) {
2070 evaluatedPermutations = 0;
2071 evaluatedPermutationsPrint = 1000;
2072 }
2073
2074 if (evaluatedPermutations >= evaluatedPermutationsPrint) {
2075 if (evaluatedPermutations % evaluatedPermutationsPrint == 0) {
2076 logIfRank0(
2077 format(" The permutations have reached %10.4e.", (double) evaluatedPermutationsPrint));
2078 evaluatedPermutationsPrint *= 10;
2079 }
2080 }
2081
2082 int nResidues = residues.length;
2083 Residue residuei = residues[i];
2084 Rotamer[] rotamersi = residuei.getRotamers();
2085 int lenri = rotamersi.length;
2086 if (i < nResidues - 1) {
2087 for (int ri = 0; ri < lenri; ri++) {
2088 if (eR.check(i, ri)) {
2089 continue;
2090 }
2091 boolean deadEnd = false;
2092 for (int j = 0; j < i; j++) {
2093 int rj = currentRotamers[j];
2094 deadEnd = eR.check(j, rj, i, ri);
2095 if (deadEnd) {
2096 break;
2097 }
2098 }
2099 if (deadEnd) {
2100 continue;
2101 }
2102 currentRotamers[i] = ri;
2103 dryRun(residues, i + 1, currentRotamers);
2104 }
2105 } else {
2106
2107 for (int ri = 0; ri < lenri; ri++) {
2108 if (eR.check(i, ri)) {
2109 continue;
2110 }
2111 currentRotamers[i] = ri;
2112 boolean deadEnd = false;
2113 for (int j = 0; j < i; j++) {
2114 int rj = currentRotamers[j];
2115 deadEnd = eR.check(j, rj, i, ri);
2116 if (deadEnd) {
2117 break;
2118 }
2119 }
2120 if (!deadEnd) {
2121 evaluatedPermutations++;
2122 }
2123 }
2124 }
2125
2126 return 0.0;
2127 }
2128
2129
2130
2131
2132
2133
2134
2135
2136
2137
2138
2139
2140 private double dryRunForEnsemble(Residue[] residues, int i, int[] currentRotamers,
2141 double gmecEnergy, double[] permutationEnergies, int[][] permutations) {
2142
2143 if (i == 0) {
2144 evaluatedPermutations = 0;
2145 }
2146 int nResidues = residues.length;
2147 Residue residuei = residues[i];
2148 Rotamer[] rotamersi = residuei.getRotamers();
2149 int lenri = rotamersi.length;
2150 if (i < nResidues - 1) {
2151 for (int ri = 0; ri < lenri; ri++) {
2152 if (eR.check(i, ri)) {
2153 continue;
2154 }
2155 boolean deadEnd = false;
2156 for (int j = 0; j < i; j++) {
2157 int rj = currentRotamers[j];
2158 deadEnd = eR.check(j, rj, i, ri);
2159 if (deadEnd) {
2160 break;
2161 }
2162 }
2163 if (deadEnd) {
2164 continue;
2165 }
2166 currentRotamers[i] = ri;
2167 dryRunForEnsemble(residues, i + 1, currentRotamers, gmecEnergy, permutationEnergies,
2168 permutations);
2169 }
2170 } else {
2171
2172 for (int ri = 0; ri < lenri; ri++) {
2173 if (eR.check(i, ri)) {
2174 continue;
2175 }
2176 currentRotamers[i] = ri;
2177 boolean deadEnd = false;
2178 for (int j = 0; j < i; j++) {
2179 int rj = currentRotamers[j];
2180 deadEnd = eR.check(j, rj, i, ri);
2181 if (deadEnd) {
2182 break;
2183 }
2184 }
2185 if (deadEnd) {
2186 continue;
2187 }
2188 if (permutationEnergies[evaluatedPermutations] - gmecEnergy < ensembleEnergy) {
2189 permutations[evaluatedPermutations] = new int[nResidues];
2190 arraycopy(currentRotamers, 0, permutations[evaluatedPermutations], 0, nResidues);
2191 }
2192 evaluatedPermutations++;
2193 }
2194 }
2195 return 0.0;
2196 }
2197
2198
2199
2200
2201
2202
2203 public int[] getOptimumRotamers() {
2204 return optimum;
2205 }
2206
2207
2208
2209
2210
2211
2212
2213 private double independent(List<Residue> residues) {
2214 double e = 0.0;
2215 List<Residue> singletonResidue = new ArrayList<>(Collections.nCopies(1, null));
2216 for (int i = 0; i < residues.size(); i++) {
2217 Residue residue = residues.get(i);
2218 singletonResidue.set(0, residue);
2219 logger.info(format(" Optimizing %s side-chain.", residue));
2220 Rotamer[] rotamers = residue.getRotamers();
2221 e = Double.MAX_VALUE;
2222 int bestRotamer = 0;
2223 double startingEnergy = 0.0;
2224 for (int j = 0; j < rotamers.length; j++) {
2225 Rotamer rotamer = rotamers[j];
2226 RotamerLibrary.applyRotamer(residue, rotamer);
2227 if (algorithmListener != null) {
2228 algorithmListener.algorithmUpdate(molecularAssembly);
2229 }
2230 double newE = Double.NaN;
2231 try {
2232 if (rotamer.isTitrating) {
2233 newE = currentEnergy(singletonResidue) + rotamer.getRotamerPhBias();
2234 } else {
2235 newE = currentEnergy(singletonResidue);
2236 }
2237
2238 if (j == 0) {
2239 startingEnergy = newE;
2240 newE = 0.0;
2241 } else {
2242 newE -= startingEnergy;
2243 }
2244 logger.info(format(" Energy %8s %-2d: %s", residue.toString(rotamers[j]), j, formatEnergy(newE)));
2245 double singularityThreshold = -100000;
2246 if (newE < singularityThreshold) {
2247 String message = format(" Rejecting as energy (%s << %s) is likely an error.", formatEnergy(newE), formatEnergy(singularityThreshold));
2248 logger.info(message);
2249 newE = Double.MAX_VALUE;
2250 }
2251 } catch (ArithmeticException ex) {
2252 logger.info(format(" Exception %s in energy calculations during independent for %s-%d", ex, residue, j));
2253 }
2254 if (newE < e) {
2255 e = newE;
2256 bestRotamer = j;
2257 }
2258 }
2259 Rotamer rotamer = rotamers[bestRotamer];
2260 RotamerLibrary.applyRotamer(residue, rotamer);
2261 optimum[i] = bestRotamer;
2262 logger.info(format(" Best Energy %8s %-2d: %s", residue.toString(rotamer), bestRotamer, formatEnergy(e)));
2263
2264 if (algorithmListener != null) {
2265 algorithmListener.algorithmUpdate(molecularAssembly);
2266 }
2267 if (terminate) {
2268 logger.info(format("\n Terminating after residue %s.\n", residue));
2269 break;
2270 }
2271 }
2272 return e;
2273 }
2274
2275
2276
2277
2278
2279
2280
2281
2282
2283 private boolean firstValidPerm(Residue[] residues, int i, int[] currentRotamers) {
2284 int nResidues = residues.length;
2285 Residue residuei = residues[i];
2286 Rotamer[] rotamersi = residuei.getRotamers();
2287 int lenri = rotamersi.length;
2288 if (i < nResidues - 1) {
2289 for (int ri = 0; ri < lenri; ri++) {
2290 if (eR.check(i, ri)) {
2291 continue;
2292 }
2293 boolean deadEnd = false;
2294 for (int j = 0; j < i; j++) {
2295 int rj = currentRotamers[j];
2296 deadEnd = eR.check(j, rj, i, ri);
2297 if (deadEnd) {
2298 break;
2299 }
2300 }
2301 if (deadEnd) {
2302 continue;
2303 }
2304 currentRotamers[i] = ri;
2305 if (firstValidPerm(residues, i + 1, currentRotamers)) {
2306 return true;
2307 }
2308 }
2309 } else {
2310
2311 for (int ri = 0; ri < lenri; ri++) {
2312 if (eR.check(i, ri)) {
2313 continue;
2314 }
2315 currentRotamers[i] = ri;
2316 boolean deadEnd = false;
2317 for (int j = 0; j < i; j++) {
2318 int rj = currentRotamers[j];
2319 deadEnd = eR.check(j, rj, i, ri);
2320 if (deadEnd) {
2321 break;
2322 }
2323 }
2324 if (!deadEnd) {
2325 break;
2326 }
2327 }
2328 return true;
2329 }
2330 return false;
2331 }
2332
2333
2334
2335
2336
2337
2338
2339 private double globalOptimization(List<Residue> residueList) {
2340 int currentEnsemble = Integer.MAX_VALUE;
2341 Residue[] residues = residueList.toArray(new Residue[0]);
2342 int nResidues = residues.length;
2343 int[] currentRotamers = new int[nResidues];
2344 optimumSubset = new int[nResidues];
2345
2346 int iterations = 0;
2347 boolean finalTry = false;
2348 int bestEnsembleTargetDiffThusFar = Integer.MAX_VALUE;
2349 double bestBufferThusFar = ensembleBuffer;
2350 double startingBuffer = ensembleBuffer;
2351
2352 if (ensembleEnergy > 0.0) {
2353 ensembleBuffer = ensembleEnergy;
2354 applyEliminationCriteria(residues, true, true);
2355
2356
2357 double permutations = 1;
2358 double singletonPermutations = 1;
2359 for (int i = 0; i < nResidues; i++) {
2360 Residue residue = residues[i];
2361 Rotamer[] rotamers = residue.getRotamers();
2362 int nr = rotamers.length;
2363 if (nr > 1) {
2364 int nrot = 0;
2365 for (int ri = 0; ri < nr; ri++) {
2366 if (!eR.eliminatedSingles[i][ri]) {
2367 nrot++;
2368 }
2369 }
2370 permutations *= rotamers.length;
2371 if (nrot > 1) {
2372 singletonPermutations *= nrot;
2373 }
2374 }
2375 }
2376 dryRun(residues, 0, currentRotamers);
2377 double pairTotalElimination = singletonPermutations - (double) evaluatedPermutations;
2378 double afterPairElim = singletonPermutations - pairTotalElimination;
2379 if (evaluatedPermutations == 0) {
2380 logger.severe(
2381 " No valid path through rotamer space found; try recomputing without pruning or using ensemble.");
2382 }
2383 if (rank0 && printFiles && ensembleFile == null) {
2384 File file = molecularAssembly.getFile();
2385 String filename = FilenameUtils.removeExtension(file.getAbsolutePath());
2386 ensembleFile = new File(filename + ".ens");
2387 if (ensembleFile.exists()) {
2388 for (int i = 2; i < 1000; i++) {
2389 ensembleFile = new File(filename + ".ens_" + i);
2390 if (!ensembleFile.exists()) {
2391 break;
2392 }
2393 }
2394 if (ensembleFile.exists()) {
2395 logger.warning(
2396 format(" Versioning failed: appending to end of file %s", ensembleFile.getName()));
2397 }
2398 }
2399 ensembleFilter = new PDBFilter(new File(ensembleFile.getName()), molecularAssembly, null,
2400 null);
2401 logger.info(format(" Ensemble file: %s", ensembleFile.getName()));
2402 }
2403 logIfRank0(format("%30s %35s %35s", "Condition", "Number of Permutations Left",
2404 "Number of Permutations Removed"));
2405 logIfRank0(format("%30s %35s %35s", "No Eliminations", permutations, ""));
2406 logIfRank0(format("%30s %35s %35s", "Single Eliminations", singletonPermutations,
2407 permutations - singletonPermutations));
2408 logIfRank0(
2409 format("%30s %35s %35s", "Pair Eliminations", afterPairElim, pairTotalElimination));
2410 logIfRank0(
2411 format("%30s %35s %35s", "Single and Pair Eliminations", (double) evaluatedPermutations,
2412 pairTotalElimination + (permutations - singletonPermutations)));
2413
2414 logIfRank0("\n Energy of permutations:");
2415 logIfRank0(format(" %12s %25s %25s", "Permutation", "Energy", "Lowest Possible Energy"));
2416
2417 double e;
2418 if (useMonteCarlo()) {
2419 firstValidPerm(residues, 0, currentRotamers);
2420 arraycopy(currentRotamers, 0, optimumSubset, 0, nResidues);
2421 rotamerOptimizationMC(residues, optimumSubset, currentRotamers, nMCSteps, false, mcUseAll);
2422
2423 logIfRank0(" Ensembles not currently compatible with Monte Carlo search");
2424
2425 } else {
2426 double[] permutationEnergies = new double[evaluatedPermutations];
2427 ensembleStates = new ArrayList<>();
2428
2429 e = rotamerOptimizationDEE(molecularAssembly, residues, 0, currentRotamers, Double.MAX_VALUE,
2430 optimumSubset, permutationEnergies);
2431 int[][] acceptedPermutations = new int[evaluatedPermutations][];
2432 logIfRank0(format(
2433 "\n Checking permutations for distance < %5.3f kcal/mol from GMEC energy %10.8f kcal/mol",
2434 ensembleEnergy, e));
2435 dryRunForEnsemble(residues, 0, currentRotamers, e, permutationEnergies,
2436 acceptedPermutations);
2437 int numAcceptedPermutations = 0;
2438
2439 for (int i = 0; i < acceptedPermutations.length; i++) {
2440 if (acceptedPermutations[i] != null) {
2441 ++numAcceptedPermutations;
2442 logIfRank0(
2443 format(" Accepting permutation %d at %8.6f < %8.6f", i, permutationEnergies[i] - e,
2444 ensembleEnergy));
2445 for (int j = 0; j < nResidues; j++) {
2446 Residue residuej = residues[j];
2447 Rotamer[] rotamersj = residuej.getRotamers();
2448 RotamerLibrary.applyRotamer(residuej, rotamersj[acceptedPermutations[i][j]]);
2449 }
2450
2451 ResidueState[] states = ResidueState.storeAllCoordinates(residues);
2452 ensembleStates.add(new ObjectPair<>(states, permutationEnergies[i]));
2453
2454 if (printFiles && rank0) {
2455 try {
2456 FileWriter fw = new FileWriter(ensembleFile, true);
2457 BufferedWriter bw = new BufferedWriter(fw);
2458 bw.write(format("MODEL %d", numAcceptedPermutations));
2459 for (int j = 0; j < 75; j++) {
2460 bw.write(" ");
2461 }
2462 bw.newLine();
2463 bw.flush();
2464 ensembleFilter.writeFile(ensembleFile, true);
2465 bw.write("ENDMDL");
2466 for (int j = 0; j < 64; j++) {
2467 bw.write(" ");
2468 }
2469 bw.newLine();
2470 bw.close();
2471 } catch (IOException ex) {
2472 logger.warning(format(" Exception writing to file: %s", ensembleFile.getName()));
2473 }
2474 }
2475 }
2476 }
2477 logIfRank0(format(" Number of permutations within %5.3f kcal/mol of GMEC energy: %6.4e",
2478 ensembleEnergy, (double) numAcceptedPermutations));
2479 ensembleStates.sort(null);
2480 }
2481
2482 logIfRank0(" Final rotamers:");
2483 logIfRank0(format("%17s %10s %11s %12s %11s", "Residue", "Chi 1", "Chi 2", "Chi 3", "Chi 4"));
2484 for (int i = 0; i < nResidues; i++) {
2485 Residue residue = residues[i];
2486 Rotamer[] rotamers = residue.getRotamers();
2487 int ri = optimumSubset[i];
2488 Rotamer rotamer = rotamers[ri];
2489 logIfRank0(format(" %c (%7s,%2d) %s", residue.getChainID(), residue.toString(rotamer), ri,
2490 rotamer.toAngleString()));
2491 RotamerLibrary.applyRotamer(residue, rotamer);
2492 }
2493 logIfRank0("\n");
2494
2495 double sumSelfEnergy = 0;
2496 double sumPairEnergy = 0;
2497 double sumTrimerEnergy = 0;
2498 for (int i = 0; i < nResidues; i++) {
2499 Residue residue = residues[i];
2500 Rotamer[] rotamers = residue.getRotamers();
2501 int ri = optimumSubset[i];
2502 sumSelfEnergy += eE.getSelf(i, ri);
2503 logIfRank0(
2504 format(" Final self Energy (%8s,%2d): %12.4f", residue.toString(rotamers[ri]), ri,
2505 eE.getSelf(i, ri)));
2506 }
2507 for (int i = 0; i < nResidues - 1; i++) {
2508 Residue residueI = residues[i];
2509 Rotamer[] rotI = residueI.getRotamers();
2510 int ri = optimumSubset[i];
2511 for (int j = i + 1; j < nResidues; j++) {
2512 Residue residueJ = residues[j];
2513 Rotamer[] rotJ = residueJ.getRotamers();
2514 int rj = optimumSubset[j];
2515 sumPairEnergy += eE.get2Body(i, ri, j, rj);
2516 if (eE.get2Body(i, ri, j, rj) > 10.0) {
2517 logIfRank0(format(" Large Final Pair Energy (%8s,%2d) (%8s,%2d): %12.4f",
2518 residueI.toString(rotI[ri]), ri, residueJ.toString(rotJ[rj]), rj,
2519 eE.get2Body(i, ri, j, rj)));
2520 }
2521 }
2522 }
2523
2524 try {
2525
2526 e = currentEnergy(residueList) + eE.getTotalRotamerPhBias(residueList, optimumSubset);
2527 } catch (ArithmeticException ex) {
2528 e = Double.NaN;
2529 logger.severe(
2530 format(" Exception %s in calculating current energy at the end of triples", ex));
2531 }
2532
2533 logIfRank0(format(" %12s %25s %25s", "Type", "Energy", "Lowest Possible Energy"));
2534 logIfRank0(format(" %12s %25f %25s", "Self:", sumSelfEnergy, ""));
2535 logIfRank0(format(" %12s %25f %25s", "Pair:", sumPairEnergy, ""));
2536
2537 approximateEnergy = eE.getBackboneEnergy() + sumSelfEnergy + sumPairEnergy;
2538
2539 if (threeBodyTerm) {
2540 for (int i = 0; i < nResidues - 2; i++) {
2541 int ri = optimumSubset[i];
2542 for (int j = i + 1; j < nResidues - 1; j++) {
2543 int rj = optimumSubset[j];
2544 for (int k = j + 1; k < nResidues; k++) {
2545 int rk = optimumSubset[k];
2546 try {
2547 sumTrimerEnergy += eE.get3Body(residues, i, ri, j, rj, k, rk);
2548 } catch (Exception ex) {
2549 logger.warning(ex.toString());
2550 }
2551 }
2552 }
2553 }
2554 approximateEnergy += sumTrimerEnergy;
2555 double higherOrderEnergy =
2556 e - sumSelfEnergy - sumPairEnergy - sumTrimerEnergy - eE.getBackboneEnergy();
2557 logIfRank0(format(" %12s %25f %25s", "Trimer:", sumTrimerEnergy, ""));
2558 logIfRank0(format(" %12s %25f %25s", "Neglected:", higherOrderEnergy, ""));
2559 } else {
2560 double higherOrderEnergy = e - sumSelfEnergy - sumPairEnergy - eE.getBackboneEnergy();
2561 logIfRank0(format(" %12s %25f %25s", "Neglected:", higherOrderEnergy, ""));
2562 }
2563 logIfRank0(format(" %12s %25f %25s", "Approximate:", approximateEnergy, ""));
2564 return e;
2565 }
2566
2567
2568
2569 int nPerms = 1;
2570 for (Residue residue : residues) {
2571 Rotamer[] rotamers = residue.getRotamers();
2572 int nr = rotamers.length;
2573 if (nr > 1) {
2574 nPerms *= rotamers.length;
2575 }
2576 if (nPerms > ensembleNumber) {
2577 break;
2578 }
2579 }
2580
2581 if (nPerms < ensembleNumber) {
2582 logger.info(format(
2583 " Requested an ensemble of %d, but only %d permutations exist; returning full ensemble",
2584 ensembleNumber, nPerms));
2585 ensembleNumber = nPerms;
2586 }
2587
2588 while (currentEnsemble != ensembleNumber) {
2589 if (monteCarlo) {
2590 logIfRank0(" Ensemble search not currently compatible with Monte Carlo");
2591 ensembleNumber = 1;
2592 }
2593
2594 if (iterations == 0) {
2595 applyEliminationCriteria(residues, true, true);
2596 } else {
2597 applyEliminationCriteria(residues, false, false);
2598 }
2599
2600
2601
2602 double permutations = 1;
2603 double singletonPermutations = 1;
2604 for (int i = 0; i < nResidues; i++) {
2605 Residue residue = residues[i];
2606 Rotamer[] rotamers = residue.getRotamers();
2607 int nr = rotamers.length;
2608 if (nr > 1) {
2609 int nrot = 0;
2610 for (int ri = 0; ri < nr; ri++) {
2611 if (!eR.eliminatedSingles[i][ri]) {
2612 nrot++;
2613 }
2614 }
2615 permutations *= rotamers.length;
2616 if (nrot > 1) {
2617 singletonPermutations *= nrot;
2618 }
2619 }
2620 }
2621
2622 logIfRank0(" Collecting Permutations:");
2623 dryRun(residues, 0, currentRotamers);
2624
2625 double pairTotalElimination = singletonPermutations - (double) evaluatedPermutations;
2626 double afterPairElim = singletonPermutations - pairTotalElimination;
2627 currentEnsemble = evaluatedPermutations;
2628 if (ensembleNumber == 1 && currentEnsemble == 0) {
2629 logger.severe(
2630 " No valid path through rotamer space found; try recomputing without pruning or using ensemble.");
2631 }
2632 if (ensembleNumber > 1) {
2633 if (rank0 && printFiles && ensembleFile == null) {
2634 File file = molecularAssembly.getFile();
2635 String filename = FilenameUtils.removeExtension(file.getAbsolutePath());
2636 ensembleFile = new File(filename + ".ens");
2637 if (ensembleFile.exists()) {
2638 for (int i = 2; i < 1000; i++) {
2639 ensembleFile = new File(filename + ".ens_" + i);
2640 if (!ensembleFile.exists()) {
2641 break;
2642 }
2643 }
2644 if (ensembleFile.exists()) {
2645 logger.warning(
2646 format(" Versioning failed: appending to end of file %s", ensembleFile.getName()));
2647 }
2648 }
2649 ensembleFilter = new PDBFilter(new File(ensembleFile.getName()), molecularAssembly, null,
2650 null);
2651 logger.info(format(" Ensemble file: %s", ensembleFile.getName()));
2652 }
2653 logIfRank0(format(" Ensemble Search Stats: (buffer: %5.3f, current: %d, target: %d)",
2654 ensembleBuffer, currentEnsemble, ensembleNumber));
2655 }
2656 if (ensembleNumber == 1 || finalTry) {
2657 logIfRank0(format("%30s %35s %35s", "Condition", "Number of Permutations Left",
2658 "Number of Permutations Removed"));
2659 logIfRank0(format("%30s %35s %35s", "No Eliminations", permutations, ""));
2660 logIfRank0(format("%30s %35s %35s", "Single Eliminations", singletonPermutations,
2661 permutations - singletonPermutations));
2662 logIfRank0(
2663 format("%30s %35s %35s", "Pair Eliminations", afterPairElim, pairTotalElimination));
2664 logIfRank0(
2665 format("%30s %35s %35s", "Single and Pair Eliminations", (double) evaluatedPermutations,
2666 pairTotalElimination + (permutations - singletonPermutations)));
2667
2668 logIfRank0("\n Energy of permutations:");
2669 logIfRank0(format(" %12s %25s %25s", "Permutation", "Energy", "Lowest Possible Energy"));
2670
2671 break;
2672 }
2673 if (abs(currentEnsemble - ensembleNumber) < bestEnsembleTargetDiffThusFar) {
2674 bestEnsembleTargetDiffThusFar = Math.abs(currentEnsemble - ensembleNumber);
2675 bestBufferThusFar = ensembleBuffer;
2676 }
2677 if (currentEnsemble > ensembleNumber) {
2678 ensembleBuffer -= ensembleBufferStep;
2679 ensembleBufferStep -= (ensembleBufferStep * 0.01);
2680 iterations++;
2681 } else if (currentEnsemble < ensembleNumber) {
2682 ensembleBuffer += ensembleBufferStep;
2683 ensembleBufferStep -= (ensembleBufferStep * 0.01);
2684 iterations++;
2685 }
2686 if (iterations > 100) {
2687 if (currentEnsemble == 0) {
2688
2689
2690 logIfRank0(" Ensemble still empty; increasing buffer energy.");
2691 startingBuffer = 3 * startingBuffer;
2692 setEnsemble(10, startingBuffer);
2693 iterations = 0;
2694 } else {
2695 ensembleBuffer = bestBufferThusFar;
2696 finalTry = true;
2697 }
2698 }
2699 }
2700
2701 if (currentEnsemble == 0) {
2702 logger.warning(
2703 " No valid rotamer permutations found; results will be unreliable. Try increasing the starting ensemble buffer.");
2704 }
2705 double[] permutationEnergyStub = null;
2706
2707 if (useMonteCarlo()) {
2708 firstValidPerm(residues, 0, currentRotamers);
2709 rotamerOptimizationMC(residues, optimumSubset, currentRotamers, nMCSteps, false, mcUseAll);
2710 } else {
2711 rotamerOptimizationDEE(molecularAssembly, residues, 0, currentRotamers, Double.MAX_VALUE,
2712 optimumSubset, permutationEnergyStub);
2713 }
2714
2715 double[] residueEnergy = new double[nResidues];
2716
2717 double sumSelfEnergy = 0;
2718 double sumLowSelfEnergy = 0;
2719
2720 logIfRank0("\n Energy contributions:");
2721 logIfRank0(format(" %15s %25s %25s", "Type", "Energy", "Lowest Possible Energy"));
2722
2723 for (int i = 0; i < nResidues; i++) {
2724 int ri = optimumSubset[i];
2725 Residue residue = residues[i];
2726 Rotamer[] rotamers = residue.getRotamers();
2727 turnOnAtoms(residue);
2728 RotamerLibrary.applyRotamer(residue, rotamers[ri]);
2729 double self = eE.getSelf(i, ri);
2730 residueEnergy[i] = self;
2731 sumSelfEnergy += self;
2732 double lowest = eE.lowestSelfEnergy(residues, i);
2733 sumLowSelfEnergy += lowest;
2734 if (self - lowest > 10.0) {
2735 logIfRank0(
2736 format(" %15s %25f %25f", "Self (" + residues[i] + "," + ri + "):", self, lowest));
2737 }
2738 }
2739
2740 double sumPairEnergy = 0.0;
2741 double sumLowPairEnergy = 0.0;
2742 double[] resPairEnergy = new double[nResidues];
2743 double[] lowPairEnergy = new double[nResidues];
2744 for (int i = 0; i < nResidues - 1; i++) {
2745 StringBuilder sb = new StringBuilder();
2746 int ri = optimumSubset[i];
2747 double sumPairEnergyI = 0;
2748 double sumLowPairEnergyI = 0;
2749 for (int j = i + 1; j < nResidues; j++) {
2750 int rj = optimumSubset[j];
2751 double pair = eE.get2Body(i, ri, j, rj);
2752 residueEnergy[i] += 0.5 * pair;
2753 residueEnergy[j] += 0.5 * pair;
2754 sumPairEnergy += pair;
2755 sumPairEnergyI += pair;
2756 double lowest = eE.lowestPairEnergy(residues, i, ri, j);
2757 sumLowPairEnergy += lowest;
2758 sumLowPairEnergyI += lowest;
2759 resPairEnergy[i] = 0.5 * pair;
2760 resPairEnergy[j] = 0.5 * pair;
2761 lowPairEnergy[i] = 0.5 * lowest;
2762 lowPairEnergy[j] = 0.5 * lowest;
2763 if (resPairEnergy[i] - lowPairEnergy[i] > 10.0) {
2764 sb.append(format(" Pair Energy (%8s,%2d) (%8s,%2d): %12.4f (Lowest: %12.4f).\n",
2765 residues[i].toFormattedString(false, true), ri,
2766 residues[j].toFormattedString(false, true), rj, pair, lowest));
2767 }
2768 }
2769 if (sumPairEnergyI - sumLowPairEnergyI > 10.0) {
2770 logIfRank0(
2771 format(" %15s %25f %25f", "Self (" + residues[i] + "," + ri + ")", sumPairEnergyI,
2772 sumLowPairEnergyI));
2773 sb.trimToSize();
2774 if (!sb.toString().isEmpty()) {
2775 logIfRank0(sb.toString());
2776 }
2777 }
2778 }
2779
2780 double e = Double.NaN;
2781 try {
2782
2783 e = currentEnergy(residueList) + eE.getTotalRotamerPhBias(residueList, optimumSubset);
2784 } catch (ArithmeticException ex) {
2785 logger.severe(
2786 format(" Exception %s in calculating current energy at the end of self and pairs", ex));
2787 }
2788 logIfRank0(format(" %15s %25f %25s", "Backbone:", eE.getBackboneEnergy(), ""));
2789 logIfRank0(format(" %15s %25f %25f", "Self:", sumSelfEnergy, sumLowSelfEnergy));
2790 logIfRank0(format(" %15s %25f %25f", "Pair:", sumPairEnergy, sumLowPairEnergy));
2791
2792 approximateEnergy = eE.getBackboneEnergy() + sumSelfEnergy + sumPairEnergy;
2793
2794 double sumTrimerEnergy = 0;
2795 if (threeBodyTerm) {
2796 for (int i = 0; i < nResidues - 2; i++) {
2797 int ri = optimumSubset[i];
2798 for (int j = i + 1; j < nResidues - 1; j++) {
2799 int rj = optimumSubset[j];
2800 for (int k = j + 1; k < nResidues; k++) {
2801 int rk = optimumSubset[k];
2802 try {
2803 double triple = eE.get3Body(residues, i, ri, j, rj, k, rk);
2804 double thirdTrip = triple / 3.0;
2805 residueEnergy[i] += thirdTrip;
2806 residueEnergy[j] += thirdTrip;
2807 residueEnergy[k] += thirdTrip;
2808 sumTrimerEnergy += triple;
2809 } catch (Exception ex) {
2810 logger.warning(ex.toString());
2811 }
2812 }
2813 }
2814 }
2815 approximateEnergy += sumTrimerEnergy;
2816 double higherOrderEnergy =
2817 e - sumSelfEnergy - sumPairEnergy - sumTrimerEnergy - eE.getBackboneEnergy();
2818 logIfRank0(format(" %15s %25f %25s", "Trimer:", sumTrimerEnergy, ""));
2819 logIfRank0(format(" %15s %25f %25s", "Neglected:", higherOrderEnergy, ""));
2820 } else {
2821 double higherOrderEnergy = e - sumSelfEnergy - sumPairEnergy - eE.getBackboneEnergy();
2822 logIfRank0(format(" %15s %25f %25s", "Neglected:", higherOrderEnergy, ""));
2823 }
2824
2825 logIfRank0(format(" %15s %25f %25s", "Approximate:", approximateEnergy, ""));
2826
2827 logIfRank0("\n Final rotamers:");
2828 logIfRank0(
2829 format("%17s %10s %11s %12s %11s %14s", "Residue", "Chi 1", "Chi 2", "Chi 3", "Chi 4",
2830 "Energy"));
2831 for (int i = 0; i < nResidues; i++) {
2832 Residue residue = residues[i];
2833 Rotamer[] rotamers = residue.getRotamers();
2834 int ri = optimumSubset[i];
2835 Rotamer rotamer = rotamers[ri];
2836 logIfRank0(format(" %3d %c (%7s,%2d) %s %12.4f ", i + 1, residue.getChainID(),
2837 residue.toString(rotamers[ri]), ri, rotamer.toAngleString(), residueEnergy[i]));
2838 RotamerLibrary.applyRotamer(residue, rotamer);
2839 }
2840 logIfRank0("\n");
2841 return e;
2842 }
2843
2844
2845
2846
2847
2848
2849
2850 private boolean useMonteCarlo() {
2851 return monteCarlo && (mcNoEnum || (nMCSteps < evaluatedPermutations));
2852 }
2853
2854
2855
2856
2857
2858
2859
2860 private double bruteForce(List<Residue> residueList) {
2861
2862 Residue[] residues = residueList.toArray(new Residue[0]);
2863 int nResidues = residues.length;
2864
2865
2866 double permutations = 1;
2867 for (Residue residue : residues) {
2868 Rotamer[] rotamers = residue.getRotamers();
2869 permutations *= rotamers.length;
2870 }
2871
2872 logger.info(format(" Number of permutations: %16.8e.", permutations));
2873
2874 double e;
2875 useForceFieldEnergy = molecularAssembly.getProperties()
2876 .getBoolean("ro-use-force-field-energy", false);
2877 if (!useForceFieldEnergy) {
2878
2879 setPruning(0);
2880 rotamerEnergies(residues);
2881 int[] currentRotamers = new int[nResidues];
2882 e = decomposedRotamerOptimization(residues, 0, Double.MAX_VALUE, optimum, currentRotamers);
2883
2884
2885 for (int i = 0; i < nResidues; i++) {
2886 Residue residue = residues[i];
2887 Rotamer[] rotamers = residue.getRotamers();
2888 RotamerLibrary.applyRotamer(residue, rotamers[optimum[i]]);
2889 turnOnAtoms(residue);
2890 }
2891
2892
2893 double fullEnergy = 0;
2894 try {
2895
2896 fullEnergy = currentEnergy(residueList) + eE.getTotalRotamerPhBias(residueList, optimum);
2897 } catch (Exception ex) {
2898 logger.severe(format(" Exception %s in calculating full energy; FFX shutting down", ex));
2899 }
2900
2901 logger.info(format(" Final summation of energies: %16.5f", e));
2902 logger.info(format(" Final energy of optimized structure: %16.5f", fullEnergy));
2903 logger.info(format(" Neglected: %16.5f", fullEnergy - e));
2904 } else {
2905 e = rotamerOptimization(molecularAssembly, residues, 0, Double.MAX_VALUE, optimum);
2906 }
2907
2908 for (int i = 0; i < nResidues; i++) {
2909 Residue residue = residues[i];
2910 Rotamer[] rotamers = residue.getRotamers();
2911 int ri = optimum[i];
2912 Rotamer rotamer = rotamers[ri];
2913 logger.info(format(" %s %s (%d)", residue.getResidueNumber(), rotamer.toString(), ri));
2914 RotamerLibrary.applyRotamer(residue, rotamer);
2915 if (useForceFieldEnergy) {
2916 try {
2917 e = currentEnergy(residueList) + eE.getTotalRotamerPhBias(residueList, optimum);
2918 } catch (ArithmeticException ex) {
2919 logger.fine(
2920 format(" Exception %s in calculating full AMOEBA energy at the end of brute force",
2921 ex));
2922 }
2923 }
2924 }
2925 return e;
2926 }
2927
2928 @SuppressWarnings("fallthrough")
2929 private double slidingWindowOptimization(List<Residue> residueList, int windowSize, int increment,
2930 boolean revert, double distance, Direction direction) {
2931
2932 long beginTime = -System.nanoTime();
2933 boolean incrementTruncated = false;
2934 boolean firstWindowSaved = false;
2935 int counter = 1;
2936 int windowEnd;
2937
2938 int nOptimize = residueList.size();
2939 if (nOptimize < windowSize) {
2940 windowSize = nOptimize;
2941 logger.info(" Window size too small for given residue range; truncating window size.");
2942 }
2943 switch (direction) {
2944 case BACKWARD:
2945 List<Residue> temp = new ArrayList<>();
2946 for (int i = nOptimize - 1; i >= 0; i--) {
2947 temp.add(residueList.get(i));
2948 }
2949 residueList = temp;
2950
2951 case FORWARD:
2952 for (int windowStart = 0; windowStart + (windowSize - 1) < nOptimize;
2953 windowStart += increment) {
2954 long windowTime = -System.nanoTime();
2955
2956 windowEnd = windowStart + (windowSize - 1);
2957
2958 logIfRank0(format("\n Iteration %d of the sliding window.\n", counter++));
2959
2960 Residue firstResidue = residueList.get(windowStart);
2961 Residue lastResidue = residueList.get(windowEnd);
2962 if (firstResidue != lastResidue) {
2963 logIfRank0(
2964 format(" Residues %s ... %s", firstResidue.toString(), lastResidue.toString()));
2965 } else {
2966 logIfRank0(format(" Residue %s", firstResidue.toString()));
2967 }
2968 List<Residue> currentWindow = new ArrayList<>();
2969 for (int i = windowStart; i <= windowEnd; i++) {
2970 Residue residue = residueList.get(i);
2971 currentWindow.add(residue);
2972 }
2973
2974 if (distance > 0) {
2975 for (int i = windowStart; i <= windowEnd; i++) {
2976 Residue residuei = residueList.get(i);
2977 int indexI = allResiduesList.indexOf(residuei);
2978 int lengthRi;
2979 lengthRi = residuei.getRotamers().length;
2980 for (int ri = 0; ri < lengthRi; ri++) {
2981 for (int j = 0; j < nAllResidues; j++) {
2982 Residue residuej = allResiduesArray[j];
2983 Rotamer[] rotamersj = residuej.getRotamers();
2984 if (currentWindow.contains(residuej) || rotamersj == null) {
2985 continue;
2986 }
2987 int lengthRj = rotamersj.length;
2988 for (int rj = 0; rj < lengthRj; rj++) {
2989 double rotamerSeparation = dM.get2BodyDistance(indexI, ri, j, rj);
2990 if (rotamerSeparation <= distance) {
2991 if (!currentWindow.contains(residuej)) {
2992 logIfRank0(format(" Adding residue %s at distance %16.8f Ang from %s %d.",
2993 residuej.toFormattedString(false, true), rotamerSeparation,
2994 residuei.toFormattedString(false, true), ri));
2995 currentWindow.add(residuej);
2996 }
2997 break;
2998 }
2999 }
3000 }
3001 }
3002 }
3003 }
3004
3005
3006
3007
3008
3009
3010
3011
3012
3013
3014
3015
3016
3017
3018
3019
3020
3021
3022
3023
3024
3025 if (counter > 2 && windowSize <= increment && firstResidue.getResidueType() == NA) {
3026 Residue prevResidue = firstResidue.getPreviousResidue();
3027 if (prevResidue != null && prevResidue.getResidueType() == NA
3028 && !currentWindow.contains(prevResidue) && prevResidue.getRotamers() != null) {
3029 logIfRank0(format(
3030 " Adding nucleic acid residue 5' of window start %s to give it flexibility about its sugar pucker.",
3031 prevResidue));
3032 currentWindow.add(prevResidue);
3033 }
3034 }
3035 sortResidues(currentWindow);
3036
3037 if (revert) {
3038 ResidueState[] coordinates = ResidueState.storeAllCoordinates(currentWindow);
3039 double startingEnergy = Double.NaN;
3040 try {
3041 startingEnergy = currentEnergy(currentWindow);
3042 } catch (ArithmeticException ex) {
3043 logger.severe(format(
3044 " Exception %s in calculating starting energy of a window; FFX shutting down",
3045 ex));
3046 }
3047 globalOptimization(currentWindow);
3048 double finalEnergy = Double.NaN;
3049 try {
3050 finalEnergy = currentEnergy(currentWindow);
3051 } catch (ArithmeticException ex) {
3052 logger.severe(
3053 format(" Exception %s in calculating final energy of a window; FFX shutting down",
3054 ex));
3055 }
3056 if (startingEnergy <= finalEnergy) {
3057 logger.info(
3058 "Optimization did not yield a better energy. Reverting to original coordinates.");
3059 ResidueState.revertAllCoordinates(currentWindow, coordinates);
3060 } else {
3061
3062 int i = 0;
3063 for (Residue residue : currentWindow) {
3064 int index = residueList.indexOf(residue);
3065 optimum[index] = optimumSubset[i++];
3066 }
3067 }
3068 } else {
3069 globalOptimization(currentWindow);
3070
3071 int i = 0;
3072 for (Residue residue : currentWindow) {
3073 int index = residueList.indexOf(residue);
3074 optimum[index] = optimumSubset[i++];
3075 }
3076 }
3077
3078 if (!incrementTruncated) {
3079 if (windowStart + (windowSize - 1) + increment > nOptimize - 1) {
3080 increment = nOptimize - windowStart - windowSize;
3081 if (increment == 0) {
3082 break;
3083 }
3084 logger.warning(" Increment truncated in order to optimize entire residue range.");
3085 incrementTruncated = true;
3086 }
3087 }
3088
3089 if (rank0 && printFiles) {
3090 File file = molecularAssembly.getFile();
3091 if (firstWindowSaved) {
3092 file.delete();
3093 }
3094
3095 if (windowStart + windowSize == nOptimize) {
3096 continue;
3097 }
3098 PDBFilter windowFilter = new PDBFilter(file, molecularAssembly, null, null);
3099
3100
3101 try {
3102 windowFilter.writeFile(file, false);
3103 if (firstResidue != lastResidue) {
3104 logger.info(
3105 format(" File with residues %s ... %s in window written to.", firstResidue,
3106 lastResidue));
3107 } else {
3108 logger.info(format(" File with residue %s in window written to.", firstResidue));
3109 }
3110 } catch (Exception e) {
3111 logger.warning(format("Exception writing to file: %s", file.getName()));
3112 }
3113 firstWindowSaved = true;
3114 }
3115 long currentTime = System.nanoTime();
3116 windowTime += currentTime;
3117 logIfRank0(format(" Time elapsed for this iteration: %11.3f sec", windowTime * 1.0E-9));
3118 logIfRank0(
3119 format(" Overall time elapsed: %11.3f sec", (currentTime + beginTime) * 1.0E-9));
3120 }
3121 break;
3122
3123 default:
3124
3125 break;
3126 }
3127 return 0.0;
3128 }
3129
3130
3131
3132
3133
3134
3135 private void sortResidues(List<Residue> residues) {
3136 Comparator<Residue> comparator = Comparator.comparing(Residue::getChainID)
3137 .thenComparingInt(Residue::getResidueNumber);
3138 residues.sort(comparator);
3139 }
3140
3141
3142
3143
3144
3145
3146
3147
3148 private void applyEliminationCriteria(Residue[] residues, boolean getEnergies, boolean verbose) {
3149 Level prevLevel = logger.getLevel();
3150 if (!verbose) {
3151 logger.setLevel(Level.WARNING);
3152 }
3153 if (getEnergies) {
3154 applyEliminationCriteria(residues);
3155 } else {
3156 int i = 0;
3157 boolean pairEliminated;
3158 do {
3159 pairEliminated = false;
3160 if (useGoldstein) {
3161 if (selfEliminationOn) {
3162 i++;
3163 logIfRank0(format("\n Iteration %d: Applying Single Goldstein DEE conditions ", i));
3164
3165 while (goldsteinDriver(residues)) {
3166 i++;
3167 logIfRank0(this.toString());
3168 logIfRank0(
3169 format("\n Iteration %d: Applying Single Rotamer Goldstein DEE conditions ", i));
3170 }
3171 }
3172 if (pairEliminationOn) {
3173 i++;
3174 logIfRank0(
3175 format("\n Iteration %d: Applying Rotamer Pair Goldstein DEE conditions ", i));
3176
3177 pairEliminated = false;
3178 while (goldsteinPairDriver(residues)) {
3179 pairEliminated = true;
3180 i++;
3181 logIfRank0(this.toString());
3182 logIfRank0(
3183 format("\n Iteration %d: Applying Rotamer Pair Goldstein DEE conditions ", i));
3184 }
3185 }
3186 } else {
3187 if (selfEliminationOn) {
3188 i++;
3189 logIfRank0(format("\n Iteration %d: Applying Single DEE conditions ", i));
3190
3191 while (deeRotamerElimination(residues)) {
3192 i++;
3193 logIfRank0(toString());
3194 logIfRank0(format("\n Iteration %d: Applying Single Rotamer DEE conditions ", i));
3195 }
3196 }
3197 if (pairEliminationOn) {
3198 i++;
3199 logIfRank0(format("\n Iteration %d: Applying Rotamer Pair DEE conditions ", i));
3200
3201 pairEliminated = false;
3202 while (deeRotamerPairElimination(residues)) {
3203 pairEliminated = true;
3204 i++;
3205 logIfRank0(toString());
3206 logIfRank0(format("\n Iteration %d: Applying Rotamer Pair DEE conditions ", i));
3207 }
3208 }
3209 }
3210 eR.validateDEE(residues);
3211 logIfRank0(toString());
3212 } while (pairEliminated);
3213 logIfRank0(" Self-consistent DEE rotamer elimination achieved.\n");
3214 }
3215 if (!verbose) {
3216 logger.setLevel(prevLevel);
3217 }
3218 }
3219
3220 private void applyEliminationCriteria(Residue[] residues) {
3221 if (verboseEnergies) {
3222 try {
3223 logIfRank0(format("\n Beginning Energy %s", formatEnergy(currentEnergy(residues))));
3224 } catch (ArithmeticException ex) {
3225 logger.severe(format(" Exception %s in calculating beginning energy; FFX shutting down.", ex));
3226 }
3227 }
3228
3229 rotamerEnergies(residues);
3230
3231 if (testing) {
3232 int nres = residues.length;
3233 eR.onlyPrunedSingles = new boolean[nres][];
3234 eR.onlyPrunedPairs = new boolean[nres][][][];
3235 for (int i = 0; i < nres; i++) {
3236 Residue residuei = residues[i];
3237 Rotamer[] rotamersi = residuei.getRotamers();
3238 int lenri = rotamersi.length;
3239 eR.onlyPrunedSingles[i] = new boolean[lenri];
3240 eR.onlyPrunedSingles[i] = copyOf(eR.eliminatedSingles[i], eR.eliminatedSingles[i].length);
3241 eR.onlyPrunedPairs[i] = new boolean[lenri][][];
3242
3243 for (int ri = 0; ri < lenri; ri++) {
3244 eR.onlyPrunedPairs[i][ri] = new boolean[nres][];
3245 for (int j = i + 1; j < nres; j++) {
3246 Residue residuej = residues[j];
3247 Rotamer[] rotamersj = residuej.getRotamers();
3248 int lenrj = rotamersj.length;
3249 eR.onlyPrunedPairs[i][ri][j] = new boolean[lenrj];
3250 eR.onlyPrunedPairs[i][ri][j] = copyOf(eR.eliminatedPairs[i][ri][j],
3251 eR.eliminatedPairs[i][ri][j].length);
3252 }
3253 }
3254 }
3255 }
3256
3257 if (testSelfEnergyEliminations) {
3258 testSelfEnergyElimination(residues);
3259 } else if (testPairEnergyEliminations > -1) {
3260 testPairEnergyElimination(residues, testPairEnergyEliminations);
3261 } else if (testTripleEnergyEliminations1 > -1 && testTripleEnergyEliminations2 > -1) {
3262 testTripleEnergyElimination(residues, testTripleEnergyEliminations1,
3263 testTripleEnergyEliminations2);
3264 }
3265
3266 if (pruneClashes) {
3267 eR.validateDEE(residues);
3268 }
3269
3270 int i = 0;
3271 boolean pairEliminated;
3272 do {
3273 pairEliminated = false;
3274 if (useGoldstein) {
3275 if (selfEliminationOn) {
3276 i++;
3277 logIfRank0(format("\n Iteration %d: Applying Single Goldstein DEE conditions ", i));
3278
3279 while (goldsteinDriver(residues)) {
3280 i++;
3281 logIfRank0(this.toString());
3282 logIfRank0(
3283 format("\n Iteration %d: Applying Single Rotamer Goldstein DEE conditions ", i));
3284 }
3285 }
3286 if (pairEliminationOn) {
3287 i++;
3288 logIfRank0(format("\n Iteration %d: Applying Rotamer Pair Goldstein DEE conditions ", i));
3289
3290 while (goldsteinPairDriver(residues)) {
3291 pairEliminated = true;
3292 i++;
3293 logIfRank0(this.toString());
3294 logIfRank0(
3295 format("\n Iteration %d: Applying Rotamer Pair Goldstein DEE conditions ", i));
3296 }
3297 }
3298 } else {
3299 if (selfEliminationOn) {
3300 i++;
3301 logIfRank0(format("\n Iteration %d: Applying Single DEE conditions ", i));
3302
3303 while (deeRotamerElimination(residues)) {
3304 i++;
3305 logIfRank0(toString());
3306 logIfRank0(format("\n Iteration %d: Applying Single Rotamer DEE conditions ", i));
3307 }
3308 }
3309 if (pairEliminationOn) {
3310 i++;
3311 logIfRank0(format("\n Iteration %d: Applying Rotamer Pair DEE conditions ", i));
3312
3313 while (deeRotamerPairElimination(residues)) {
3314 pairEliminated = true;
3315 i++;
3316 logIfRank0(toString());
3317 logIfRank0(format("\n Iteration %d: Applying Rotamer Pair DEE conditions ", i));
3318 }
3319 }
3320 }
3321 eR.validateDEE(residues);
3322 logIfRank0(toString());
3323 } while (pairEliminated);
3324
3325 logIfRank0(" Self-consistent DEE rotamer elimination achieved.\n");
3326 }
3327
3328
3329
3330
3331
3332
3333
3334 private double currentEnergy(List<Residue> resList) throws ArithmeticException {
3335 List<Rotamer> rots = resList.stream().filter(Objects::nonNull).map(Residue::getRotamer)
3336 .collect(Collectors.toList());
3337 File energyDir = dirSupplier.apply(resList, rots);
3338 return eFunction.applyAsDouble(energyDir);
3339 }
3340
3341
3342
3343
3344
3345
3346
3347 private double currentPE(File dir) {
3348 if (x == null) {
3349 int nVar = potential.getNumberOfVariables();
3350 x = new double[nVar];
3351 }
3352 potential.getCoordinates(x);
3353 return potential.energy(x);
3354 }
3355
3356
3357
3358
3359
3360
3361 private void generateResidueNeighbors(Residue[] residues) {
3362 int nRes = residues.length;
3363 resNeighbors = new int[nRes][];
3364 bidiResNeighbors = new int[nRes][];
3365 for (int i = 0; i < nRes; i++) {
3366 Set<Integer> nearby = new HashSet<>();
3367 Residue resi = residues[i];
3368 Rotamer[] rotsi = resi.getRotamers();
3369 int lenri = rotsi.length;
3370 int indexI = allResiduesList.indexOf(resi);
3371
3372 for (int j = 0; j < nRes; j++) {
3373 if (i == j) {
3374 continue;
3375 }
3376 Residue resj = residues[j];
3377 Rotamer[] rotsj = resj.getRotamers();
3378 int lenrj = rotsj.length;
3379 int indexJ = allResiduesList.indexOf(resj);
3380
3381 boolean foundClose = false;
3382 for (int ri = 0; ri < lenri; ri++) {
3383 for (int rj = 0; rj < lenrj; rj++) {
3384 if (!dM.checkPairDistThreshold(indexI, ri, indexJ, rj)) {
3385 foundClose = true;
3386 break;
3387 }
3388 }
3389 if (foundClose) {
3390 break;
3391 }
3392 }
3393 if (foundClose) {
3394 nearby.add(j);
3395 }
3396 }
3397
3398
3399 int[] nI = nearby.stream().mapToInt(Integer::intValue).toArray();
3400 bidiResNeighbors[i] = nI;
3401
3402
3403 final int fi = i;
3404 nI = nearby.stream().mapToInt(Integer::intValue).filter(j -> j > fi).toArray();
3405 resNeighbors[i] = nI;
3406 }
3407 }
3408
3409 private void setUpRestart() {
3410 File restartFile;
3411 if (loadEnergyRestart) {
3412 restartFile = energyRestartFile;
3413 } else {
3414 File file = molecularAssembly.getFile();
3415 String filename = FilenameUtils.removeExtension(file.getAbsolutePath());
3416 Path restartPath = Paths.get(filename + ".restart");
3417 restartFile = restartPath.toFile();
3418 energyRestartFile = restartFile;
3419 }
3420 try {
3421 energyWriter = new BufferedWriter(new FileWriter(restartFile, true));
3422 } catch (IOException ex) {
3423 logger.log(Level.SEVERE, "Couldn't open energy restart file.", ex);
3424 }
3425 logger.info(format("\n Energy restart file: %s", restartFile.getName()));
3426 }
3427
3428
3429
3430
3431
3432
3433
3434
3435 private double rotamerEnergies(Residue[] residues) {
3436
3437 if (residues == null) {
3438 logger.warning(" Attempt to compute rotamer energies for an empty array of residues.");
3439 return 0.0;
3440 }
3441
3442 int nResidues = residues.length;
3443 Atom[] atoms = molecularAssembly.getAtomArray();
3444 generateResidueNeighbors(residues);
3445
3446 eR = new EliminatedRotamers(this, dM, allResiduesList, maxRotCheckDepth, clashThreshold,
3447 pairClashThreshold, multiResClashThreshold, nucleicPruningFactor, nucleicPairsPruningFactor,
3448 multiResPairClashAddn, pruneClashes, prunePairClashes, print, residues);
3449
3450 if (decomposeOriginal) {
3451 assert library.getUsingOrigCoordsRotamer();
3452 for (int i = 0; i < nResidues; i++) {
3453 Residue resi = residues[i];
3454 Rotamer[] rotsi = resi.getRotamers();
3455 int lenri = rotsi.length;
3456
3457
3458 for (int ri = 1; ri < lenri; ri++) {
3459 eR.eliminateRotamer(residues, i, ri, false);
3460 }
3461 }
3462 }
3463
3464
3465 for (Atom atom : atoms) {
3466 atom.setUse(true);
3467 }
3468
3469 eE = new EnergyExpansion(this, dM, eR, molecularAssembly, potential, algorithmListener,
3470 allResiduesList, resNeighbors, threeBodyTerm, decomposeOriginal, usingBoxOptimization,
3471 verbose, pruneClashes, prunePairClashes, rank0);
3472
3473
3474 eR.setEnergyExpansion(eE);
3475
3476 int loaded = 0;
3477 if (loadEnergyRestart) {
3478 if (usingBoxOptimization) {
3479 loaded = eE.loadEnergyRestart(energyRestartFile, residues, boxOpt.boxLoadIndex, boxOpt.boxLoadCellIndices);
3480 } else {
3481 loaded = eE.loadEnergyRestart(energyRestartFile, residues);
3482 }
3483 }
3484
3485 long energyStartTime = System.nanoTime();
3486 WorkerTeam energyWorkerTeam = new WorkerTeam(world);
3487
3488 try {
3489 if (loaded < 1) {
3490 eE.allocateSelfJobMap(residues, nResidues, false);
3491 }
3492
3493 SelfEnergyRegion selfEnergyRegion = new SelfEnergyRegion(this, eE, eR, residues, energyWriter,
3494 world, numProc, pruneClashes, rank0, rank, verbose, writeEnergyRestart, printFiles);
3495 energyWorkerTeam.execute(selfEnergyRegion);
3496 long singlesTime = System.nanoTime() - energyStartTime;
3497 logIfRank0(format(" Time for single energies: %12.4g", (singlesTime * 1.0E-9)));
3498 if (logger.isLoggable(Level.FINE)) {
3499 Resources.logResources();
3500 }
3501
3502 if (loaded < 2) {
3503 eE.allocate2BodyJobMap(residues, nResidues, false);
3504 }
3505
3506 TwoBodyEnergyRegion twoBodyEnergyRegion = new TwoBodyEnergyRegion(this, dM, eE, eR, residues,
3507 allResiduesList, energyWriter, world, numProc, prunePairClashes, superpositionThreshold,
3508 rank0, rank, verbose, writeEnergyRestart, printFiles);
3509 energyWorkerTeam.execute(twoBodyEnergyRegion);
3510 long pairsTime = System.nanoTime() - (singlesTime + energyStartTime);
3511
3512 long triplesTime = 0;
3513 long quadsTime = 0;
3514 logIfRank0(format(" Time for 2-body energies: %12.4g", (pairsTime * 1.0E-9)));
3515 if (logger.isLoggable(Level.FINE)) {
3516 Resources.logResources();
3517 }
3518
3519 if (threeBodyTerm) {
3520 if (loaded < 3) {
3521 eE.allocate3BodyJobMap(residues, nResidues, false);
3522 }
3523 ThreeBodyEnergyRegion threeBodyEnergyRegion = new ThreeBodyEnergyRegion(this, dM, eE, eR,
3524 residues, allResiduesList, energyWriter, world, numProc, superpositionThreshold, rank0,
3525 rank, verbose, writeEnergyRestart, printFiles);
3526 energyWorkerTeam.execute(threeBodyEnergyRegion);
3527 triplesTime = System.nanoTime() - (pairsTime + singlesTime + energyStartTime);
3528 logIfRank0(format(" Time for 3-Body energies: %12.4g", (triplesTime * 1.0E-9)));
3529 if (logger.isLoggable(Level.FINE)) {
3530 Resources.logResources();
3531 }
3532 }
3533
3534 if (compute4BodyEnergy) {
3535 eE.allocate4BodyJobMap(residues, nResidues);
3536
3537 FourBodyEnergyRegion fourBodyEnergyRegion = new FourBodyEnergyRegion(this, dM, eE, eR,
3538 residues, allResiduesList, superpositionThreshold);
3539 energyWorkerTeam.execute(fourBodyEnergyRegion);
3540 quadsTime = System.nanoTime() - (triplesTime + pairsTime + singlesTime + energyStartTime);
3541 logIfRank0(format(" Time for 4-Body energies: %12.4g", quadsTime * 1.0E-9));
3542 }
3543
3544 long allTime = singlesTime + pairsTime + triplesTime + quadsTime;
3545 logIfRank0(format(" Time for all energies: %12.4g", allTime * 1.0E-9));
3546 } catch (Exception ex) {
3547 String message = " Exception computing rotamer energies in parallel.";
3548 logger.log(Level.SEVERE, message, ex);
3549 }
3550
3551
3552 for (Atom atom : atoms) {
3553 atom.setUse(true);
3554 }
3555
3556 if (verboseEnergies && rank0) {
3557 try {
3558 double defaultEnergy = currentEnergy(residues);
3559 logger.info(format(" Energy of the system with rotamers in their default conformation: %s",
3560 formatEnergy(defaultEnergy)));
3561 } catch (ArithmeticException ex) {
3562 logger.severe(format(" Exception %s in calculating default energy; FFX shutting down", ex));
3563 }
3564 }
3565 return eE.getBackboneEnergy();
3566 }
3567
3568
3569
3570
3571
3572
3573 private void applyDefaultRotamer(Residue residue) {
3574 RotamerLibrary.applyRotamer(residue, residue.getRotamers()[0]);
3575 }
3576
3577
3578
3579
3580
3581
3582
3583 private boolean deeRotamerElimination(Residue[] residues) {
3584 int nres = residues.length;
3585
3586 boolean eliminated = false;
3587
3588 double[] minMax = new double[2];
3589 double[] minEnergySingles = null;
3590 double[] maxEnergySingles = null;
3591 for (int i = 0; i < nres; i++) {
3592 Residue residuei = residues[i];
3593 Rotamer[] rotamersi = residuei.getRotamers();
3594 int lenri = rotamersi.length;
3595 if (minEnergySingles == null || minEnergySingles.length < lenri) {
3596 minEnergySingles = new double[lenri];
3597 maxEnergySingles = new double[lenri];
3598 }
3599
3600 for (int ri = 0; ri < lenri; ri++) {
3601
3602 if (eR.check(i, ri)) {
3603 continue;
3604 }
3605
3606 minEnergySingles[ri] = eE.getSelf(i, ri);
3607 maxEnergySingles[ri] = minEnergySingles[ri];
3608 for (int j = 0; j < nres; j++) {
3609 if (j == i) {
3610 continue;
3611 }
3612 if (eE.minMaxPairEnergy(residues, minMax, i, ri, j)) {
3613 if (Double.isFinite(minMax[0]) && Double.isFinite(minEnergySingles[ri])) {
3614 minEnergySingles[ri] += minMax[0];
3615 } else {
3616
3617
3618
3619 minEnergySingles[ri] = Double.NaN;
3620 }
3621 if (Double.isFinite(minMax[0]) && Double.isFinite(maxEnergySingles[ri])) {
3622 maxEnergySingles[ri] += minMax[1];
3623 } else {
3624
3625
3626 maxEnergySingles[ri] = Double.NaN;
3627 }
3628 } else {
3629 Residue residuej = residues[j];
3630 logger.info(
3631 format(" Inconsistent Pair: %8s %2d, %8s.", residuei.toFormattedString(false, true),
3632 ri, residuej.toFormattedString(false, true)));
3633
3634 }
3635 }
3636 }
3637
3638
3639
3640
3641
3642 double eliminationEnergy = Double.MAX_VALUE;
3643 int eliminatingRotamer = 0;
3644 for (int ri = 0; ri < lenri; ri++) {
3645 if (eR.check(i, ri)) {
3646 continue;
3647 }
3648 if (Double.isFinite(maxEnergySingles[ri]) && maxEnergySingles[ri] < eliminationEnergy) {
3649 eliminationEnergy = maxEnergySingles[ri];
3650 eliminatingRotamer = ri;
3651 }
3652 }
3653
3654 if (eliminationEnergy == Double.MAX_VALUE) {
3655
3656
3657 logIfRank0(" Could not eliminate any i,ri because eliminationEnergy was never set!",
3658 Level.FINE);
3659 } else {
3660
3661
3662
3663
3664 for (int ri = 0; ri < lenri; ri++) {
3665 if (eR.check(i, ri)) {
3666 continue;
3667 }
3668
3669
3670 if (!Double.isFinite(minEnergySingles[ri])) {
3671 if (eR.eliminateRotamer(residues, i, ri, print)) {
3672 logIfRank0(format(" Rotamer elimination of (%8s,%2d) that always clashes.",
3673 residuei.toFormattedString(false, true), ri));
3674 eliminated = true;
3675 }
3676 }
3677
3678
3679 if (minEnergySingles[ri] > eliminationEnergy + ensembleBuffer) {
3680 if (eR.eliminateRotamer(residues, i, ri, print)) {
3681 logIfRank0(format(" Rotamer elimination of (%8s,%2d) by (%8s,%2d): %12.4f > %6.4f.",
3682 residuei.toFormattedString(false, true), ri,
3683 residuei.toFormattedString(false, true), eliminatingRotamer, minEnergySingles[ri],
3684 eliminationEnergy + ensembleBuffer));
3685 eliminated = true;
3686 }
3687 }
3688 }
3689 }
3690 }
3691 return eliminated;
3692 }
3693
3694
3695
3696
3697
3698
3699
3700
3701 private boolean deeRotamerPairElimination(Residue[] residues) {
3702 int nres = residues.length;
3703 boolean eliminated = false;
3704
3705 for (int i = 0; i < (nres - 1); i++) {
3706 Residue residuei = residues[i];
3707 Rotamer[] rotamersi = residuei.getRotamers();
3708 int lenri = rotamersi.length;
3709
3710
3711 double[][] minPairEnergies = new double[lenri][];
3712 double[][] maxPairEnergies = new double[lenri][];
3713
3714 for (int j = i + 1; j < nres; j++) {
3715 Residue residuej = residues[j];
3716 Rotamer[] rotamersj = residuej.getRotamers();
3717 int lenrj = rotamersj.length;
3718
3719 for (int ri = 0; ri < lenri; ri++) {
3720 if (eR.check(i, ri)) {
3721 continue;
3722 }
3723 minPairEnergies[ri] = new double[lenrj];
3724 maxPairEnergies[ri] = new double[lenrj];
3725
3726 for (int rj = 0; rj < lenrj; rj++) {
3727 if (eR.check(j, rj) || eR.check(i, ri, j, rj)) {
3728 continue;
3729 }
3730 minPairEnergies[ri][rj] =
3731 eE.getSelf(i, ri) + eE.getSelf(j, rj) + eE.get2Body(i, ri, j, rj);
3732 maxPairEnergies[ri][rj] = minPairEnergies[ri][rj];
3733
3734
3735 double[] minMax = new double[2];
3736
3737
3738 if (eE.minMaxE2(residues, minMax, i, ri, j, rj)) {
3739 if (Double.isFinite(minPairEnergies[ri][rj]) && Double.isFinite(minMax[0])) {
3740 minPairEnergies[ri][rj] += minMax[0];
3741 } else {
3742 logger.severe(
3743 format(" An ri-rj pair %s-%d %s-%d with NaN minimum was caught incorrectly!",
3744 residuei.toFormattedString(false, true), ri,
3745 residuej.toFormattedString(false, true), rj));
3746 }
3747 if (Double.isFinite(maxPairEnergies[ri][rj]) && Double.isFinite(minMax[1])) {
3748 maxPairEnergies[ri][rj] += minMax[1];
3749 } else {
3750
3751 maxPairEnergies[ri][rj] = Double.NaN;
3752 }
3753 } else {
3754
3755 minPairEnergies[ri][rj] = Double.NaN;
3756 logger.info(format(" Eliminating pair %s-%d %s-%d that always clashes.",
3757 residuei.toFormattedString(false, true), ri,
3758 residuej.toFormattedString(false, true), rj));
3759 eR.eliminateRotamerPair(residues, i, ri, j, rj, print);
3760 eliminated = true;
3761 }
3762 }
3763 }
3764
3765 double pairEliminationEnergy = Double.MAX_VALUE;
3766 for (int ri = 0; ri < lenri; ri++) {
3767 if (eR.check(i, ri)) {
3768 continue;
3769 }
3770 for (int rj = 0; rj < lenrj; rj++) {
3771 if (eR.check(j, rj) || eR.check(i, ri, j, rj)) {
3772 continue;
3773 }
3774 if (Double.isFinite(maxPairEnergies[ri][rj])
3775 && maxPairEnergies[ri][rj] < pairEliminationEnergy) {
3776 pairEliminationEnergy = maxPairEnergies[ri][rj];
3777 }
3778 }
3779 }
3780
3781 if (pairEliminationEnergy == Double.MAX_VALUE) {
3782 logIfRank0(format(
3783 " All rotamer pairs for residues %s and %s have possible conflicts; cannot perform any eliminations!",
3784 residuei.toFormattedString(false, true), residuej), Level.FINE);
3785 } else {
3786 double comparisonEnergy = pairEliminationEnergy + ensembleBuffer;
3787 for (int ri = 0; ri < lenri; ri++) {
3788 if (eR.check(i, ri)) {
3789 continue;
3790 }
3791 for (int rj = 0; rj < lenrj; rj++) {
3792 if (eR.check(j, rj) || eR.check(i, ri, j, rj)) {
3793 continue;
3794 }
3795 if (minPairEnergies[ri][rj] > comparisonEnergy) {
3796 if (eR.eliminateRotamerPair(residues, i, ri, j, rj, print)) {
3797 eliminated = true;
3798 logIfRank0(format(" Eliminating rotamer pair: %s %d, %s %d (%s > %s + %6.6f)",
3799 residuei.toFormattedString(false, true), ri,
3800 residuej.toFormattedString(false, true), rj,
3801 formatEnergy(minPairEnergies[ri][rj]), formatEnergy(pairEliminationEnergy),
3802 ensembleBuffer), Level.INFO);
3803 } else {
3804
3805 logIfRank0(
3806 format(" Already eliminated rotamer pair! %s %d, %s %d (%s > %1s + %6.6f)",
3807 residuei.toFormattedString(false, true), ri,
3808 residuej.toFormattedString(false, true), rj,
3809 formatEnergy(minPairEnergies[ri][rj]), formatEnergy(pairEliminationEnergy),
3810 ensembleBuffer), Level.WARNING);
3811 }
3812 }
3813 }
3814 }
3815 }
3816
3817 if (eR.pairsToSingleElimination(residues, i, j)) {
3818 eliminated = true;
3819 }
3820 }
3821 }
3822
3823 return eliminated;
3824 }
3825
3826 private boolean goldsteinDriver(Residue[] residues) {
3827 int nres = residues.length;
3828
3829 boolean eliminated = false;
3830
3831 for (int i = 0; i < nres; i++) {
3832 Residue resi = residues[i];
3833 Rotamer[] roti = resi.getRotamers();
3834 int nri = roti.length;
3835
3836 for (int riA = 0; riA < nri; riA++) {
3837
3838 if (eR.check(i, riA)) {
3839 continue;
3840 }
3841 for (int riB = 0; riB < nri; riB++) {
3842
3843 if (riA == riB || eR.check(i, riB)) {
3844 continue;
3845 }
3846 if (goldsteinElimination(residues, i, riA, riB)) {
3847 eliminated = true;
3848 break;
3849 }
3850 }
3851 }
3852 }
3853 if (!eliminated) {
3854 logIfRank0(" No more single rotamers to eliminate.");
3855 }
3856 return eliminated;
3857 }
3858
3859
3860
3861
3862
3863
3864
3865
3866
3867
3868 private boolean goldsteinElimination(Residue[] residues, int i, int riA, int riB) {
3869 Residue resi = residues[i];
3870
3871
3872 double selfDiff = eE.getSelf(i, riA) - eE.getSelf(i, riB);
3873 double goldsteinEnergy = selfDiff;
3874
3875 double sumPairDiff = 0.0;
3876 double sumTripleDiff = 0.0;
3877
3878
3879 for (int j : bidiResNeighbors[i]) {
3880 Residue resj = residues[j];
3881 Rotamer[] rotj = resj.getRotamers();
3882 int nrj = rotj.length;
3883 double minForResJ = Double.MAX_VALUE;
3884 double minPairDiff = 0.0;
3885 double minTripleDiff = 0.0;
3886 int rjEvals = 0;
3887
3888
3889 for (int rj = 0; rj < nrj; rj++) {
3890 if (eR.check(j, rj)) {
3891 continue;
3892 }
3893
3894 if (eR.check(i, riA, j, rj)) {
3895 continue;
3896 }
3897 if (eR.check(i, riB, j, rj)) {
3898
3899
3900
3901
3902
3903 return false;
3904 }
3905
3906 double pairI = eE.get2Body(i, riA, j, rj);
3907 double pairJ = eE.get2Body(i, riB, j, rj);
3908 double pairDiff = pairI - pairJ;
3909
3910 rjEvals++;
3911
3912
3913 double tripleDiff = 0.0;
3914 if (threeBodyTerm) {
3915 IntStream kStream = IntStream.concat(Arrays.stream(bidiResNeighbors[i]),
3916 Arrays.stream(bidiResNeighbors[j]));
3917 int[] possKs = kStream.distinct().sorted().toArray();
3918 for (int k : possKs) {
3919 if (k == i || k == j) {
3920 continue;
3921 }
3922 Residue resk = residues[k];
3923 Rotamer[] rotk = resk.getRotamers();
3924 int nrk = rotk.length;
3925 int rkEvals = 0;
3926 double minForResK = Double.MAX_VALUE;
3927 for (int rk = 0; rk < nrk; rk++) {
3928
3929
3930
3931
3932 if (eR.check(k, rk) || eR.check(j, rj, k, rk) || eR.check(i, riA, k, rk)) {
3933
3934
3935 continue;
3936 }
3937
3938
3939
3940
3941 if (eR.check(i, riB, k, rk)) {
3942
3943 return false;
3944 }
3945
3946 rkEvals++;
3947 double e =
3948 eE.get3Body(residues, i, riA, j, rj, k, rk) - eE.get3Body(residues, i, riB, j, rj,
3949 k, rk);
3950 if (e < minForResK) {
3951 minForResK = e;
3952 }
3953 }
3954
3955 if (rkEvals == 0) {
3956 minForResK = 0.0;
3957 }
3958 tripleDiff += minForResK;
3959 }
3960 }
3961
3962 double sumOverK = pairDiff + tripleDiff;
3963 if (sumOverK < minForResJ) {
3964 minForResJ = sumOverK;
3965 minPairDiff = pairDiff;
3966 minTripleDiff = tripleDiff;
3967 }
3968 }
3969
3970 if (rjEvals == 0) {
3971 minForResJ = 0.0;
3972 minPairDiff = 0.0;
3973 minTripleDiff = 0.0;
3974 }
3975 sumPairDiff += minPairDiff;
3976 sumTripleDiff += minTripleDiff;
3977 goldsteinEnergy += minForResJ;
3978 }
3979
3980 if (goldsteinEnergy > ensembleBuffer) {
3981 if (eR.eliminateRotamer(residues, i, riA, print)) {
3982 logIfRank0(format(" Rotamer elimination of (%8s,%2d) by (%8s,%2d): %12.4f > %6.4f.",
3983 resi.toFormattedString(false, true), riA, resi.toFormattedString(false, true), riB,
3984 goldsteinEnergy, ensembleBuffer));
3985 logIfRank0(format(" Self: %12.4f, Pairs: %12.4f, Triples: %12.4f.", selfDiff, sumPairDiff,
3986 sumTripleDiff));
3987 return true;
3988 }
3989 }
3990 return false;
3991 }
3992
3993
3994
3995
3996
3997
3998
3999 private boolean goldsteinPairDriver(Residue[] residues) {
4000 int nRes = residues.length;
4001 boolean eliminated = false;
4002
4003
4004 for (int i = 0; i < nRes; i++) {
4005 Residue resi = residues[i];
4006 Rotamer[] rotsi = resi.getRotamers();
4007 int lenri = rotsi.length;
4008
4009 for (int riA = 0; riA < lenri; riA++) {
4010
4011 if (eR.check(i, riA)) {
4012 continue;
4013 }
4014
4015
4016 for (int j = 0; j < nRes; j++) {
4017
4018 if (i == j) {
4019 continue;
4020 }
4021 Residue resj = residues[j];
4022 Rotamer[] rotsj = resj.getRotamers();
4023 int lenrj = rotsj.length;
4024 for (int rjC = 0; rjC < lenrj; rjC++) {
4025
4026 if (eR.check(j, rjC) || eR.check(i, riA, j, rjC)) {
4027 continue;
4028 }
4029 boolean breakOut = false;
4030
4031
4032
4033 for (int riB = 0; riB < lenri; riB++) {
4034 if (breakOut) {
4035 break;
4036 }
4037 if (eR.check(i, riB)) {
4038 continue;
4039 }
4040 for (int rjD = 0; rjD < lenrj; rjD++) {
4041 if (breakOut) {
4042 break;
4043 }
4044
4045 if (eR.check(j, rjD) || eR.check(i, riB, j, rjD)) {
4046 continue;
4047 }
4048
4049 if (riA == riB && rjC == rjD) {
4050 continue;
4051 }
4052 if (goldsteinPairElimination(residues, i, riA, riB, j, rjC, rjD)) {
4053 breakOut = true;
4054 eliminated = true;
4055 }
4056 }
4057 }
4058 }
4059 if (eR.pairsToSingleElimination(residues, i, j)) {
4060 eliminated = true;
4061 }
4062 }
4063 }
4064 }
4065 return eliminated;
4066 }
4067
4068
4069
4070
4071
4072
4073
4074
4075
4076
4077
4078
4079
4080 private boolean goldsteinPairElimination(Residue[] residues, int i, int riA, int riB, int j,
4081 int rjC, int rjD) {
4082
4083 List<Residue> missedResidues = null;
4084
4085
4086 double goldsteinEnergy =
4087 eE.getSelf(i, riA) + eE.getSelf(j, rjC) + eE.get2Body(i, riA, j, rjC) - eE.getSelf(i, riB)
4088 - eE.getSelf(j, rjD) - eE.get2Body(i, riB, j, rjD);
4089
4090 try {
4091 if (parallelTeam == null) {
4092 parallelTeam = new ParallelTeam();
4093 }
4094 if (goldsteinPairRegion == null) {
4095 goldsteinPairRegion = new GoldsteinPairRegion(parallelTeam.getThreadCount());
4096 }
4097 goldsteinPairRegion.init(residues, i, riA, riB, j, rjC, rjD, bidiResNeighbors, this);
4098 parallelTeam.execute(goldsteinPairRegion);
4099 goldsteinEnergy += goldsteinPairRegion.getSumOverK();
4100 missedResidues = goldsteinPairRegion.getMissedResidues();
4101 } catch (Exception e) {
4102 logger.log(Level.WARNING, " Exception in GoldsteinPairRegion.", e);
4103 }
4104 if (missedResidues != null && !missedResidues.isEmpty()) {
4105 logIfRank0(format(
4106 " Skipping energy comparison due to a missed residue: i %d riA %d riB %d j %d rjC %d rjD %d",
4107 i, riA, riB, j, rjC, rjD), Level.FINE);
4108 return false;
4109 }
4110
4111 if (goldsteinEnergy > ensembleBuffer) {
4112 if (missedResidues.isEmpty()) {
4113 if (eR.eliminateRotamerPair(residues, i, riA, j, rjC, print)) {
4114 logIfRank0(format(
4115 " Pair elimination of [(%8s,%2d),(%8s,%2d)] by [(%8s,%2d),(%8s,%2d)]: %12.4f > %6.4f",
4116 residues[i].toFormattedString(false, true), riA,
4117 residues[j].toFormattedString(false, true), rjC,
4118 residues[i].toFormattedString(false, true), riB,
4119 residues[j].toFormattedString(false, true), rjD, goldsteinEnergy, ensembleBuffer));
4120 return true;
4121 }
4122 } else {
4123 logIfRank0(format(
4124 " No Pair elimination of [(%8s,%2d),(%8s,%2d)] by [(%8s,%2d),(%8s,%2d)]: %12.4f > %6.4f",
4125 residues[i].toFormattedString(false, true), riA,
4126 residues[j].toFormattedString(false, true), rjC,
4127 residues[i].toFormattedString(false, true), riB,
4128 residues[j].toFormattedString(false, true), rjD, goldsteinEnergy, ensembleBuffer));
4129 StringBuilder sb = new StringBuilder();
4130 for (Residue residueM : missedResidues) {
4131 sb.append(residueM);
4132 }
4133 logIfRank0(format(" due to %s.", sb));
4134 }
4135 }
4136 return false;
4137 }
4138
4139
4140
4141
4142
4143
4144 void setTestOverallOpt(boolean testing) {
4145 this.testing = testing;
4146 distanceMethod = DistanceMethod.ROTAMER;
4147 setTwoBodyCutoff(Double.MAX_VALUE);
4148 setThreeBodyCutoff(9.0);
4149 }
4150
4151
4152
4153
4154
4155
4156
4157
4158 void setTestSelfEnergyEliminations(boolean testSelfEnergyEliminations) {
4159 this.testing = true;
4160 this.testSelfEnergyEliminations = testSelfEnergyEliminations;
4161 testPairEnergyEliminations = -1;
4162 testTripleEnergyEliminations1 = -1;
4163 testTripleEnergyEliminations2 = -1;
4164 distanceMethod = DistanceMethod.ROTAMER;
4165 setTwoBodyCutoff(Double.MAX_VALUE);
4166 setThreeBodyCutoff(9.0);
4167 }
4168
4169
4170
4171
4172
4173
4174
4175
4176 void setTestPairEnergyEliminations(int testPairEnergyEliminations) {
4177 this.testing = true;
4178 this.testPairEnergyEliminations = testPairEnergyEliminations;
4179 testSelfEnergyEliminations = false;
4180 testTripleEnergyEliminations1 = -1;
4181 testTripleEnergyEliminations2 = -1;
4182 distanceMethod = DistanceMethod.ROTAMER;
4183 setTwoBodyCutoff(Double.MAX_VALUE);
4184 setThreeBodyCutoff(9.0);
4185 }
4186
4187
4188
4189
4190
4191
4192
4193
4194
4195
4196 void setTestTripleEnergyEliminations(int testTripleEnergyEliminations1,
4197 int testTripleEnergyEliminations2) {
4198 this.testing = true;
4199 this.testTripleEnergyEliminations1 = testTripleEnergyEliminations1;
4200 this.testTripleEnergyEliminations2 = testTripleEnergyEliminations2;
4201 testSelfEnergyEliminations = false;
4202 testPairEnergyEliminations = -1;
4203 distanceMethod = DistanceMethod.ROTAMER;
4204 setTwoBodyCutoff(Double.MAX_VALUE);
4205 setThreeBodyCutoff(9.0);
4206 }
4207
4208
4209
4210
4211
4212
4213 private void testSelfEnergyElimination(Residue[] residues) {
4214 int nRes = residues.length;
4215 for (int i = 0; i < nRes; i++) {
4216 Residue resI = residues[i];
4217 Rotamer[] rotI = resI.getRotamers();
4218 int nI = rotI.length;
4219 for (int ri = 0; ri < nI; ri++) {
4220 for (int j = i + 1; j < nRes; j++) {
4221 Residue resJ = residues[j];
4222 Rotamer[] rotJ = resJ.getRotamers();
4223 int nJ = rotJ.length;
4224 for (int rj = 0; rj < nJ; rj++) {
4225 try {
4226 eE.set2Body(i, ri, j, rj, 0, true);
4227 } catch (Exception e) {
4228
4229 }
4230 if (threeBodyTerm) {
4231 for (int k = j + 1; k < nRes; k++) {
4232 Residue resK = residues[k];
4233 Rotamer[] rotK = resK.getRotamers();
4234 int nK = rotK.length;
4235 for (int rk = 0; rk < nK; rk++) {
4236 try {
4237 eE.set3Body(residues, i, ri, j, rj, k, rk, 0, true);
4238 } catch (Exception e) {
4239
4240 }
4241 }
4242 }
4243 }
4244 }
4245 }
4246 }
4247 }
4248 }
4249
4250
4251
4252
4253
4254
4255
4256 private void testPairEnergyElimination(Residue[] residues, int resID) {
4257 int nRes = residues.length;
4258
4259 if (resID >= nRes) {
4260 return;
4261 }
4262
4263 for (int i = 0; i < nRes; i++) {
4264 Residue resI = residues[i];
4265 Rotamer[] rotI = resI.getRotamers();
4266 int nI = rotI.length;
4267 for (int ri = 0; ri < nI; ri++) {
4268 try {
4269 eE.setSelf(i, ri, 0, true);
4270 } catch (Exception e) {
4271
4272 }
4273 for (int j = i + 1; j < nRes; j++) {
4274 Residue resJ = residues[j];
4275 Rotamer[] rotJ = resJ.getRotamers();
4276 int nJ = rotJ.length;
4277 for (int rj = 0; rj < nJ; rj++) {
4278 if (i != resID && j != resID) {
4279 try {
4280 eE.set2Body(i, ri, j, rj, 0, true);
4281 } catch (Exception e) {
4282
4283 }
4284 }
4285 if (threeBodyTerm) {
4286 for (int k = j + 1; k < nRes; k++) {
4287 Residue resK = residues[k];
4288 Rotamer[] rotK = resK.getRotamers();
4289 int nK = rotK.length;
4290 for (int rk = 0; rk < nK; rk++) {
4291 try {
4292 eE.set3Body(residues, i, ri, j, rj, k, rk, 0, true);
4293 } catch (Exception e) {
4294
4295 }
4296 }
4297 }
4298 }
4299 }
4300 }
4301 }
4302 }
4303 }
4304
4305
4306
4307
4308
4309
4310
4311
4312
4313 private void testTripleEnergyElimination(Residue[] residues, int resID1, int resID2) {
4314 int nRes = residues.length;
4315
4316 if (resID1 >= nRes) {
4317 return;
4318 }
4319 if (resID2 >= nRes) {
4320 return;
4321 }
4322 if (resID1 == resID2) {
4323 return;
4324 }
4325
4326 for (int i = 0; i < nRes; i++) {
4327 Residue resI = residues[i];
4328 Rotamer[] rotI = resI.getRotamers();
4329 int nI = rotI.length;
4330 for (int ri = 0; ri < nI; ri++) {
4331 try {
4332 eE.setSelf(i, ri, 0, true);
4333 } catch (Exception e) {
4334
4335 }
4336 for (int j = i + 1; j < nRes; j++) {
4337 Residue resJ = residues[j];
4338 Rotamer[] rotJ = resJ.getRotamers();
4339 int nJ = rotJ.length;
4340 for (int rj = 0; rj < nJ; rj++) {
4341
4342
4343
4344
4345 try {
4346 eE.set2Body(i, ri, j, rj, 0, true);
4347 } catch (Exception e) {
4348
4349 }
4350 if (threeBodyTerm) {
4351 for (int k = j + 1; k < nRes; k++) {
4352 Residue resK = residues[k];
4353 Rotamer[] rotK = resK.getRotamers();
4354 int nK = rotK.length;
4355 for (int rk = 0; rk < nK; rk++) {
4356
4357 if (i != resID1 && j != resID1 && k != resID1) {
4358 try {
4359 eE.set3Body(residues, i, ri, j, rj, k, rk, 0, true);
4360 } catch (Exception e) {
4361
4362 }
4363 }
4364 if (i != resID2 && j != resID2 && k != resID2) {
4365 try {
4366 eE.set3Body(residues, i, ri, j, rj, k, rk, 0, true);
4367 } catch (Exception e) {
4368
4369 }
4370 }
4371 }
4372 }
4373 }
4374 }
4375 }
4376 }
4377 }
4378 }
4379
4380
4381
4382
4383 public enum Algorithm {
4384 INDEPENDENT,
4385 ALL,
4386 BRUTE_FORCE,
4387 WINDOW,
4388 BOX;
4389
4390 public static Algorithm getAlgorithm(int algorithm) {
4391 switch (algorithm) {
4392 case 1:
4393 return INDEPENDENT;
4394 case 2:
4395 return ALL;
4396 case 3:
4397 return BRUTE_FORCE;
4398 case 4:
4399 return WINDOW;
4400 case 5:
4401 return BOX;
4402 default:
4403 throw new IllegalArgumentException(
4404 format(" Algorithm choice was %d, not in range 1-5!", algorithm));
4405 }
4406 }
4407 }
4408
4409
4410
4411
4412 public enum DistanceMethod {
4413 ROTAMER, RESIDUE
4414 }
4415
4416 public enum Direction {
4417 FORWARD, BACKWARD
4418 }
4419
4420 private class BoxOptimization {
4421
4422
4423
4424
4425 public final int[] numXYZCells = {3, 3, 3};
4426
4427
4428
4429 public double cellBorderSize = 0;
4430
4431
4432
4433 public double approxBoxLength = 0;
4434
4435
4436
4437 public int boxInclusionCriterion = 1;
4438
4439
4440
4441 public int cellStart = 0;
4442
4443
4444
4445 public int cellEnd = -1;
4446
4447
4448
4449 public boolean manualSuperbox = false;
4450
4451
4452
4453 public double[] boxDimensions;
4454
4455
4456
4457 public double superboxBuffer = 8.0;
4458
4459
4460
4461 public int boxLoadIndex = -1;
4462
4463
4464
4465 public int[] boxLoadCellIndices;
4466
4467
4468
4469
4470
4471
4472 public double boxOptimization(List<Residue> residueList) {
4473 usingBoxOptimization = true;
4474 long beginTime = -System.nanoTime();
4475 Residue[] residues = residueList.toArray(new Residue[0]);
4476
4477
4478
4479
4480
4481
4482
4483
4484 Crystal crystal = generateSuperbox(residueList);
4485
4486
4487 int totalCells = getTotalCellCount(crystal);
4488 if (cellStart > totalCells - 1) {
4489 logIfRank0(format(" First cell out of range (%d) -- reset to first cell.", cellStart + 1));
4490 cellStart = 0;
4491 }
4492 if (cellEnd > totalCells - 1) {
4493
4494 if (cellEnd != -1 && cellEnd != Integer.MAX_VALUE) {
4495 logIfRank0(format(" Final cell out of range (%d) -- reset to last cell.", cellEnd + 1));
4496 }
4497 cellEnd = totalCells - 1;
4498 } else if (cellEnd < 0) {
4499 cellEnd = totalCells - 1;
4500 }
4501 ManyBodyCell[] cells = loadCells(crystal, residues);
4502 int numCells = cells.length;
4503 logIfRank0(format(" Optimizing cells %d to %d", (cellStart + 1), (cellEnd + 1)));
4504 for (int i = 0; i < numCells; i++) {
4505 ManyBodyCell manyBodyCell = cells[i];
4506 List<Residue> residueSubsetList = manyBodyCell.getResiduesAsList();
4507 int[] cellIndices = manyBodyCell.getABCIndices();
4508 logIfRank0(format("\n Iteration %d of cell based optimization.", (i + 1)));
4509 logIfRank0(manyBodyCell.toString());
4510 int nResidueSubset = residueSubsetList.size();
4511 if (nResidueSubset > 0) {
4512 if (rank0 && writeEnergyRestart && printFiles) {
4513 String boxHeader = format(" Box %d: %d,%d,%d", i + 1, cellIndices[0], cellIndices[1], cellIndices[2]);
4514 try {
4515 energyWriter.append(boxHeader);
4516 energyWriter.newLine();
4517 } catch (IOException ex) {
4518 logger.log(Level.SEVERE, " Exception writing box header to energy restart file.", ex);
4519 }
4520 }
4521 if (loadEnergyRestart) {
4522 boxLoadIndex = i + 1;
4523 boxLoadCellIndices = new int[3];
4524 boxLoadCellIndices[0] = cellIndices[0];
4525 boxLoadCellIndices[1] = cellIndices[1];
4526 boxLoadCellIndices[2] = cellIndices[2];
4527 }
4528
4529 long boxTime = -System.nanoTime();
4530 Residue firstResidue = residueSubsetList.getFirst();
4531 Residue lastResidue = residueSubsetList.get(nResidueSubset - 1);
4532 if (revert) {
4533 ResidueState[] coordinates = ResidueState.storeAllCoordinates(residueSubsetList);
4534 double startingEnergy = 0;
4535 double finalEnergy = 0;
4536 try {
4537 startingEnergy = currentEnergy(residueSubsetList);
4538 } catch (ArithmeticException ex) {
4539 logger.severe(format(" Exception %s in calculating starting energy of a box; FFX shutting down", ex));
4540 }
4541 globalOptimization(residueSubsetList);
4542 try {
4543 finalEnergy = currentEnergy(residueSubsetList);
4544 } catch (ArithmeticException ex) {
4545 logger.severe(format(" Exception %s in calculating starting energy of a box; FFX shutting down", ex));
4546 }
4547 if (startingEnergy <= finalEnergy) {
4548 logger.info(
4549 "Optimization did not yield a better energy. Reverting to original coordinates.");
4550 ResidueState.revertAllCoordinates(residueSubsetList, coordinates);
4551 } else {
4552
4553 int r = 0;
4554 for (Residue residue : residueSubsetList) {
4555 int index = residueList.indexOf(residue);
4556 optimum[index] = optimumSubset[r++];
4557 }
4558 }
4559 long currentTime = System.nanoTime();
4560 boxTime += currentTime;
4561 logIfRank0(format(" Time elapsed for this iteration: %11.3f sec", boxTime * 1.0E-9));
4562 logIfRank0(format(" Overall time elapsed: %11.3f sec", (currentTime + beginTime) * 1.0E-9));
4563 } else {
4564 globalOptimization(residueSubsetList);
4565
4566 int r = 0;
4567 for (Residue residue : residueSubsetList) {
4568 int index = residueList.indexOf(residue);
4569 optimum[index] = optimumSubset[r++];
4570 }
4571 long currentTime = System.nanoTime();
4572 boxTime += currentTime;
4573 logIfRank0(format(" Time elapsed for this iteration: %11.3f sec", boxTime * 1.0E-9));
4574 logIfRank0(format(" Overall time elapsed: %11.3f sec", (currentTime + beginTime) * 1.0E-9));
4575 }
4576 if (rank0 && printFiles) {
4577
4578 if (i == (numCells - 1)) {
4579 continue;
4580 }
4581 try {
4582 if (firstResidue != lastResidue) {
4583 logIfRank0(format(" File with residues %s ... %s in window written.", firstResidue, lastResidue));
4584 } else {
4585 logIfRank0(format(" File with residue %s in window written.", firstResidue));
4586 }
4587 } catch (Exception e) {
4588 logger.warning("Exception writing to file.");
4589 }
4590 }
4591 } else {
4592 logIfRank0(" Empty box: no residues found.");
4593 }
4594 }
4595 return 0.0;
4596 }
4597
4598 public void update(String boxDim) {
4599
4600 try {
4601 String[] bdTokens = boxDim.split(",+");
4602 boxDimensions = new double[6];
4603 if (bdTokens.length != 7) {
4604 logger.warning(" Improper number of arguments to boxDimensions; default settings used.");
4605 } else {
4606 for (int i = 1; i < 7; i += 2) {
4607 boxDimensions[i - 1] = parseDouble(bdTokens[i]);
4608 boxDimensions[i] = parseDouble(bdTokens[i + 1]);
4609 if (boxDimensions[i] < boxDimensions[i - 1]) {
4610 logger.info(format(" Improper dimension min %8.5f > max %8.5f; max/min reversed.",
4611 boxDimensions[i - 1], boxDimensions[i]));
4612 double temp = boxDimensions[i];
4613 boxDimensions[i] = boxDimensions[i - 1];
4614 boxDimensions[i - 1] = temp;
4615 }
4616 }
4617 superboxBuffer = parseDouble(bdTokens[0]);
4618 manualSuperbox = true;
4619 }
4620 } catch (Exception ex) {
4621 logger.warning(
4622 format(" Error in parsing box dimensions: input discarded and defaults used: %s.", ex));
4623 manualSuperbox = false;
4624 }
4625 }
4626
4627 private Crystal generateSuperbox(List<Residue> residueList) {
4628 double[] maxXYZ = new double[3];
4629 double[] minXYZ = new double[3];
4630 Crystal originalCrystal = molecularAssembly.getCrystal();
4631 if (manualSuperbox) {
4632 for (int i = 0; i < maxXYZ.length; i++) {
4633 int ii = 2 * i;
4634 minXYZ[i] = boxDimensions[ii] - superboxBuffer;
4635 maxXYZ[i] = boxDimensions[ii + 1] + superboxBuffer;
4636 }
4637 } else if (originalCrystal.aperiodic()) {
4638 if (residueList == null || residueList.isEmpty()) {
4639 throw new IllegalArgumentException(
4640 " Null or empty residue list when generating superbox.");
4641 }
4642 Atom initializerAtom = residueList.get(0).getReferenceAtom();
4643 initializerAtom.getXYZ(minXYZ);
4644 initializerAtom.getXYZ(maxXYZ);
4645 for (Residue residue : residueList) {
4646 Atom refAtom = residue.getReferenceAtom();
4647 double[] refAtomCoords = new double[3];
4648 refAtom.getXYZ(refAtomCoords);
4649 for (int i = 0; i < 3; i++) {
4650 maxXYZ[i] = Math.max(refAtomCoords[i], maxXYZ[i]);
4651 minXYZ[i] = Math.min(refAtomCoords[i], minXYZ[i]);
4652 }
4653 }
4654 for (int i = 0; i < 3; i++) {
4655 minXYZ[i] -= superboxBuffer;
4656 maxXYZ[i] += superboxBuffer;
4657 }
4658 } else {
4659 return originalCrystal;
4660 }
4661 double newA = maxXYZ[0] - minXYZ[0];
4662 double newB = maxXYZ[1] - minXYZ[1];
4663 double newC = maxXYZ[2] - minXYZ[2];
4664 if (manualSuperbox) {
4665 logger.info(format(" Manual superbox set over (minX, maxX, minY, "
4666 + "maxY, minZ, maxZ): %f, %f, %f, %f, %f, %f", minXYZ[0], maxXYZ[0], minXYZ[1],
4667 maxXYZ[1], minXYZ[2], maxXYZ[2]));
4668 } else {
4669 logger.info(
4670 " System is aperiodic: protein box generated over these coordinates (minX, maxX, minY, maxY, minZ, maxZ):");
4671 String message = " Aperiodic box dimensions: ";
4672 for (int i = 0; i < minXYZ.length; i++) {
4673 message = message.concat(format("%f,%f,", minXYZ[i], maxXYZ[i]));
4674 }
4675 message = message.substring(0, message.length() - 1);
4676 logger.info(message);
4677 }
4678 logger.info(format(" Buffer size (included in dimensions): %f\n", superboxBuffer));
4679 return new Crystal(newA, newB, newC, 90.0, 90.0, 90.0, "P1");
4680 }
4681
4682
4683
4684
4685
4686
4687
4688
4689 private int getTotalCellCount(Crystal crystal) {
4690 int numCells = 1;
4691 if (approxBoxLength > 0) {
4692 double[] boxes = new double[3];
4693 boxes[0] = crystal.a / approxBoxLength;
4694 boxes[1] = crystal.b / approxBoxLength;
4695 boxes[2] = crystal.c / approxBoxLength;
4696 for (int i = 0; i < boxes.length; i++) {
4697 if (boxes[i] < 1) {
4698 numXYZCells[i] = 1;
4699 } else {
4700 numXYZCells[i] = (int) boxes[i];
4701 }
4702 }
4703 }
4704 for (int numXYZBox : numXYZCells) {
4705 numCells *= numXYZBox;
4706 }
4707 return numCells;
4708 }
4709
4710
4711
4712
4713
4714
4715
4716
4717 @SuppressWarnings("fallthrough")
4718 private ManyBodyCell[] loadCells(Crystal crystal, Residue[] residues) {
4719 double aCellBorderFracSize = (cellBorderSize / crystal.a);
4720 double bCellBorderFracSize = (cellBorderSize / crystal.b);
4721 double cCellBorderFracSize = (cellBorderSize / crystal.c);
4722 int numCells = cellEnd - cellStart + 1;
4723 logIfRank0(format(" Number of fractional cells: %d = %d x %d x %d", numCells, numXYZCells[0], numXYZCells[1], numXYZCells[2]));
4724
4725 ManyBodyCell[] cells = new ManyBodyCell[numCells];
4726 int currentIndex = 0;
4727 int filledCells = 0;
4728 int[] xyzIndices = new int[3];
4729 boolean doBreak = false;
4730
4731
4732
4733
4734
4735
4736 for (int i = 0; i < numXYZCells[0]; i++) {
4737 if (doBreak) {
4738 break;
4739 }
4740 xyzIndices[0] = i;
4741 for (int j = 0; j < numXYZCells[1]; j++) {
4742 if (doBreak) {
4743 break;
4744 }
4745 xyzIndices[1] = j;
4746 for (int k = 0; k < numXYZCells[2]; k++) {
4747 if (currentIndex < cellStart) {
4748 ++currentIndex;
4749 continue;
4750 } else if (currentIndex > cellEnd) {
4751 doBreak = true;
4752 break;
4753 }
4754 xyzIndices[2] = k;
4755 double[] fracCoords = new double[6];
4756 fracCoords[0] = (((1.0 * i) / numXYZCells[0]) - aCellBorderFracSize);
4757 fracCoords[1] = (((1.0 * j) / numXYZCells[1]) - bCellBorderFracSize);
4758 fracCoords[2] = (((1.0 * k) / numXYZCells[2]) - cCellBorderFracSize);
4759 fracCoords[3] = (((1.0 + i) / numXYZCells[0]) + aCellBorderFracSize);
4760 fracCoords[4] = (((1.0 + j) / numXYZCells[1]) + bCellBorderFracSize);
4761 fracCoords[5] = (((1.0 + k) / numXYZCells[2]) + cCellBorderFracSize);
4762 cells[filledCells++] = new ManyBodyCell(fracCoords, xyzIndices, currentIndex);
4763 ++currentIndex;
4764 }
4765 }
4766 }
4767 assignResiduesToCells(crystal, residues, cells);
4768 for (ManyBodyCell cell : cells) {
4769 cell.sortCellResidues();
4770 }
4771 switch (direction) {
4772 case BACKWARD:
4773 ManyBodyCell[] tempCells = new ManyBodyCell[numCells];
4774 for (int i = 0; i < numCells; i++) {
4775 tempCells[i] = cells[numCells - (i + 1)];
4776 }
4777 cells = tempCells;
4778
4779 case FORWARD:
4780 default:
4781 return cells;
4782 }
4783 }
4784
4785
4786
4787
4788
4789
4790
4791
4792
4793
4794
4795 private void assignResiduesToCells(Crystal crystal, Residue[] residues, ManyBodyCell[] cells) {
4796
4797
4798 int nSymm = crystal.spaceGroup.getNumberOfSymOps();
4799
4800 for (ManyBodyCell cell : cells) {
4801 Set<Residue> toAdd = new HashSet<>();
4802 for (int iSymm = 0; iSymm < nSymm; iSymm++) {
4803 SymOp symOp = crystal.spaceGroup.getSymOp(iSymm);
4804 for (Residue residue : residues) {
4805 boolean contained;
4806 switch (boxInclusionCriterion) {
4807 default -> contained = cell.atomInsideCell(residue.getReferenceAtom(), crystal, symOp);
4808 case 2 -> contained = cell.residueInsideCell(residue, crystal, symOp, true);
4809 case 3 -> contained = cell.anyRotamerInsideCell(residue, crystal, symOp, true);
4810 }
4811 if (contained) {
4812 toAdd.add(residue);
4813 }
4814 }
4815
4816 if (toAdd.isEmpty()) {
4817 break;
4818 }
4819 }
4820 toAdd.forEach(cell::addResidue);
4821 }
4822 }
4823 }
4824
4825
4826
4827
4828
4829
4830
4831
4832
4833
4834
4835
4836
4837
4838
4839
4840
4841
4842
4843
4844
4845
4846
4847
4848
4849
4850
4851
4852
4853
4854
4855
4856
4857
4858
4859
4860
4861
4862
4863
4864
4865
4866
4867
4868
4869
4870
4871
4872
4873
4874
4875
4876
4877
4878
4879
4880
4881
4882
4883
4884
4885
4886
4887
4888
4889
4890
4891
4892
4893
4894
4895
4896
4897
4898
4899
4900
4901
4902
4903
4904
4905
4906
4907
4908
4909
4910
4911
4912
4913
4914
4915
4916
4917
4918
4919
4920
4921
4922
4923
4924
4925
4926
4927
4928
4929
4930
4931
4932
4933
4934
4935
4936
4937
4938
4939
4940
4941
4942
4943
4944
4945
4946
4947
4948
4949
4950
4951
4952
4953
4954
4955
4956
4957
4958
4959
4960
4961
4962
4963
4964
4965
4966
4967
4968
4969
4970
4971
4972
4973
4974
4975
4976
4977
4978
4979
4980
4981
4982
4983
4984
4985
4986 }