View Javadoc
1   // ******************************************************************************
2   //
3   // Title:       Force Field X.
4   // Description: Force Field X - Software for Molecular Biophysics.
5   // Copyright:   Copyright (c) Michael J. Schnieders 2001-2024.
6   //
7   // This file is part of Force Field X.
8   //
9   // Force Field X is free software; you can redistribute it and/or modify it
10  // under the terms of the GNU General Public License version 3 as published by
11  // the Free Software Foundation.
12  //
13  // Force Field X is distributed in the hope that it will be useful, but WITHOUT
14  // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
15  // FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
16  // details.
17  //
18  // You should have received a copy of the GNU General Public License along with
19  // Force Field X; if not, write to the Free Software Foundation, Inc., 59 Temple
20  // Place, Suite 330, Boston, MA 02111-1307 USA
21  //
22  // Linking this library statically or dynamically with other modules is making a
23  // combined work based on this library. Thus, the terms and conditions of the
24  // GNU General Public License cover the whole combination.
25  //
26  // As a special exception, the copyright holders of this library give you
27  // permission to link this library with independent modules to produce an
28  // executable, regardless of the license terms of these independent modules, and
29  // to copy and distribute the resulting executable under terms of your choice,
30  // provided that you also meet, for each linked independent module, the terms
31  // and conditions of the license of that module. An independent module is a
32  // module which is not derived from or based on this library. If you modify this
33  // library, you may extend this exception to your version of the library, but
34  // you are not obligated to do so. If you do not wish to do so, delete this
35  // exception statement from your version.
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  * Optimize protein side-chain conformations and nucleic acid backbone conformations using rotamers.
114  *
115  * @author Michael J. Schnieders
116  * @author Jacob M. Litman
117  * @author Stephen D. LuCore
118  * @author Mallory R. Tollefson
119  * @since 1.0
120  */
121 public class RotamerOptimization implements Terminatable {
122 
123   /**
124    * Logger for this class.
125    */
126   private static final Logger logger = Logger.getLogger(RotamerOptimization.class.getName());
127   /**
128    * Fallback if there is no vdW node.
129    */
130   private static final double FALLBACK_TWO_BODY_CUTOFF = 0;
131   /**
132    * MolecularAssembly to perform rotamer optimization on.
133    */
134   protected final MolecularAssembly molecularAssembly;
135   /**
136    * The Potential to evaluate during rotamer optimization.
137    */
138   protected final Potential potential;
139   /**
140    * AlgorithmListener who should receive updates as the optimization runs.
141    */
142   protected final AlgorithmListener algorithmListener;
143   /**
144    * World Parallel Java communicator.
145    */
146   private final Comm world;
147   /**
148    * Number of Parallel Java processes.
149    */
150   private final int numProc;
151   /**
152    * Rank of this process.
153    */
154   private final int rank;
155   /**
156    * Flag to indicate if this is the master process.
157    */
158   private final boolean rank0;
159   /**
160    * Flag to control the verbosity of printing.
161    */
162   private final boolean print = false;
163   /**
164    * Flag to calculate and print additional energies (mostly for debugging).
165    */
166   private final boolean verboseEnergies = true;
167   /**
168    * If true, write out an energy restart file.
169    */
170   private final boolean writeEnergyRestart = true;
171   /**
172    * Parameters for box optimization are stored here.
173    */
174   private final BoxOptimization boxOpt;
175   /**
176    * Represents the method called to obtain the directory corresponding to the current energy; will
177    * be a simple return null for potential energy evaluations. While current energy calls will fill
178    * the rotamer list with the current rotamers of the residue, other methods may skip applying the
179    * rotamer directly.
180    */
181   private final BiFunction<List<Residue>, List<Rotamer>, File> dirSupplier;
182   /**
183    * Represents the method called to obtain energy for the current rotamer or state; defaults to the
184    * existing potential energy code. May discard the input file.
185    */
186   private final ToDoubleFunction<File> eFunction;
187   /**
188    * Flag to indicate verbose logging.
189    */
190   private final boolean verbose;
191   /**
192    * The DistanceMatrix class handles calculating distances between residues.
193    */
194   private DistanceMatrix dM;
195   /**
196    * The EnergyExpansion class compute terms in the many-body energy expansion.
197    */
198   private EnergyExpansion eE;
199   /**
200    * The EliminatedRotamers class tracks eliminated rotamers and rotamer paris.
201    */
202   private EliminatedRotamers eR;
203   /**
204    * RotamerLibrary instance.
205    */
206   protected RotamerLibrary library = RotamerLibrary.getDefaultLibrary();
207   /**
208    * Parallel evaluation of quantities used during Goldstein Pair elimination.
209    */
210   private GoldsteinPairRegion goldsteinPairRegion;
211   /**
212    * Parallel evaluation of many-body energy sums.
213    */
214   private EnergyRegion energyRegion;
215   /**
216    * Flag to indicate a request to terminate the optimization.
217    */
218   private boolean terminate = false;
219   /**
220    * Flag to indicate if the algorithm is running (done == false) or completed (done == true).
221    */
222   private boolean done = true;
223   /**
224    * Two-Body cutoff distance.
225    */
226   private double twoBodyCutoffDist;
227   /**
228    * Flag to control use of 3-body terms.
229    */
230   private boolean threeBodyTerm = false;
231   /**
232    * Three-body cutoff distance.
233    */
234   private double threeBodyCutoffDist;
235   /**
236    * Interaction partners of a Residue that come after it.
237    */
238   private int[][] resNeighbors;
239   /**
240    * All interaction partners of a Residue, including prior residues.
241    */
242   private int[][] bidiResNeighbors;
243   /**
244    * Flag to prune clashes.
245    */
246   private boolean pruneClashes = true;
247   /**
248    * Flag to prune pair clashes.
249    */
250   private boolean prunePairClashes = true;
251   /**
252    * Number of permutations whose energy is explicitly evaluated.
253    */
254   private int evaluatedPermutations = 0;
255   /**
256    * Permutations are printed when modulo this field is zero.
257    */
258   private int evaluatedPermutationsPrint = 0;
259   /**
260    * List of residues to optimize; they may not be contiguous or all members of the same chain.
261    */
262   private List<Residue> residueList;
263   /**
264    * This is the optimal rotamers corresponding to residueList.
265    */
266   private int[] optimum;
267 
268   /**
269    * Size of the sliding window.
270    */
271   private int windowSize = 7;
272   /**
273    * The distance the sliding window moves.
274    */
275   private int increment = 3;
276   /**
277    * In sliding window, whether to revert an unfavorable change.
278    */
279   private boolean revert;
280   /**
281    * The sliding window direction.
282    */
283   private Direction direction = Direction.FORWARD;
284 
285   /**
286    * The distance that the distance matrix checks for.
287    */
288   private double distance = 2.0;
289   /**
290    * Default distance method is to find the shortest distance between residues.
291    */
292   private DistanceMethod distanceMethod = DistanceMethod.RESIDUE;
293   /**
294    * The algorithm to use for rotamer optimization.
295    */
296   private Algorithm algorithm = null;
297 
298   /**
299    * Flag to indicate use of the Goldstein criteria instead of the less stringent Dead-End
300    * Elimination criteria.
301    */
302   private boolean useGoldstein = true;
303   /**
304    * The number of most-favorable structures to include as output.
305    */
306   private int ensembleNumber = 1;
307   /**
308    * The energy buffer applied to each elimination criteria to affect an ensemble.
309    */
310   private double ensembleBuffer = 0.0;
311   /**
312    * The step value of the energy buffer for use with ensemble search.
313    */
314   private double ensembleBufferStep = 0.5;
315   /**
316    * The energy boundary for structures to be included in the final ensemble.
317    */
318   private double ensembleEnergy = 0.0;
319   /**
320    * File to contain ensemble of structures.
321    */
322   private File ensembleFile;
323   /**
324    * PDBFilter to write out ensemble snapshots.
325    */
326   private PDBFilter ensembleFilter;
327   /**
328    * Flag to load the distance matrix as needed; if false, matrix is prefilled at the beginning of
329    * rotamer optimization.
330    */
331   private boolean lazyMatrix = false;
332   /**
333    * Clash energy threshold (kcal/mole).
334    */
335   private double clashThreshold = 25.0;
336   /**
337    * Clash energy threshold (kcal/mol) for MultiResidues, which can have much more variation in self
338    * and 2-Body energies.
339    */
340   private double multiResClashThreshold = 80.0;
341   /**
342    * Clash energy threshold (kcal/mole).
343    */
344   private double pairClashThreshold = 25.0;
345   /**
346    * Pair clash energy threshold (kcal/mol) for MultiResidues.
347    */
348   private double multiResPairClashAddn = 80.0;
349   /**
350    * An array of atomic coordinates (length 3 * the number of atoms).
351    */
352   private double[] x = null;
353   /**
354    * A flag to indicate use of the full N-Body AMOEBA potential energy during the rotamer
355    * optimization.
356    */
357   private boolean useForceFieldEnergy = false;
358   /**
359    * Threshold to eliminate nucleic acid Rotamers based on excessive correction distances; 0
360    * indicates the threshold is not being implemented.
361    */
362   private double nucleicCorrectionThreshold = 0;
363   /**
364    * The approximate energy from summing the backbone energy, self-energies, pair-energies, etc.
365    */
366   private double approximateEnergy = 0;
367   /**
368    * Minimum number of nucleic acid Rotamers to use for rotamer optimization should some be
369    * eliminated by the nucleic correction threshold.
370    */
371   private int minNumberAcceptedNARotamers = 10;
372   /**
373    * Factor by which to multiply the pruning constraints for nucleic acids.
374    */
375   private double nucleicPruningFactor = 10.0;
376   /**
377    * The arithmetic mean of 1.0 and the pruning factor, and is applied for AA-NA pairs.
378    */
379   private double nucleicPairsPruningFactor = ((1.0 + nucleicPruningFactor) / 2);
380   /**
381    * A list of all residues in the system, which is used to compute a distance matrix.
382    */
383   private List<Residue> allResiduesList = null;
384   /**
385    * An array of all residues in the system, which is used to compute a distance matrix.
386    */
387   private Residue[] allResiduesArray = null;
388   /**
389    * Number of residues being optimized.
390    */
391   private int nAllResidues = 0;
392   /**
393    * The array of optimum rotamers for the subset of residues being optimized during box or window
394    * optimization.
395    */
396   private int[] optimumSubset;
397   /**
398    * If true, load an energy restart file.
399    */
400   private boolean loadEnergyRestart = false;
401   /**
402    * Energy restart File instance.
403    */
404   private File energyRestartFile;
405   /**
406    * ParallelTeam instance.
407    */
408   private ParallelTeam parallelTeam;
409   /**
410    * Flag to indicate computation of 4-body energy values. This is limited to the study 4-body energy
411    * magnitudes, but is not included in the rotamer optimization.
412    */
413   private boolean compute4BodyEnergy = false;
414   /**
415    * Flag to indicate use of box optimization.
416    */
417   private boolean usingBoxOptimization = false;
418   /**
419    * If a pair of residues have two atoms closer together than the superposition threshold, the
420    * energy is set to NaN.
421    */
422   private double superpositionThreshold = 0.25;
423   /**
424    * Flag to indicate computation of a many-body expansion for original coordinates.
425    */
426   private boolean decomposeOriginal = false;
427   /**
428    * Flag to indicate use of MC optimization.
429    */
430   private boolean monteCarlo = false;
431   /**
432    * Number of MC optimization steps.
433    */
434   private int nMCSteps = 1000000;
435   /**
436    * MC temperature (K).
437    */
438   private double mcTemp = 298.15;
439   /**
440    * Check to see if proposed move has an eliminated 2-body or higher-order term.
441    */
442   private boolean mcUseAll = false;
443   /**
444    * Skips brute force enumeration in favor of pure Monte Carlo. Recommended only for testing.
445    */
446   private boolean mcNoEnum = false;
447   /**
448    * Sets whether files should be printed; true for standalone applications, false for some
449    * applications which use rotamer optimization as part of a larger process.
450    */
451   private boolean printFiles = true;
452   /**
453    * Stores states of each ensemble if printFiles is false.
454    */
455   private List<ObjectPair<ResidueState[], Double>> ensembleStates;
456   /**
457    * Maximum depth to check if a rotamer can be eliminated.
458    */
459   private int maxRotCheckDepth;
460   /**
461    * Writes energies to restart file.
462    */
463   private BufferedWriter energyWriter;
464 
465   /**
466    * False unless JUnit testing.
467    */
468   private boolean testing = false;
469   /**
470    * False unless ManyBodyTest is occurring.
471    */
472   private boolean monteCarloTesting = false;
473   /**
474    * Test Self-Energy Elimination.
475    */
476   private boolean testSelfEnergyEliminations = false;
477   /**
478    * Test Pair-Energy Elimination.
479    *
480    * <p>If greater than or equal to 0, test the specified residue.
481    */
482   private int testPairEnergyEliminations = -1;
483   /**
484    * Test Triple-Energy Elimination.
485    *
486    * <p>If greater than or equal to 0, test the specified residues.
487    */
488   private int testTripleEnergyEliminations1 = -1;
489   /**
490    * Test Triple-Energy Elimination.
491    *
492    * <p>If greater than or equal to 0, test the specified residues.
493    */
494   private int testTripleEnergyEliminations2 = -1;
495   /**
496    * Only for unit testing; turns off the singles elimination criterion.
497    */
498   private boolean selfEliminationOn = true;
499   /**
500    * Only for unit testing; turns off the pairs elimination criterion.
501    */
502   private boolean pairEliminationOn = true;
503 
504   /**
505    * RotamerOptimization constructor.
506    *
507    * @param molecularAssembly The MolecularAssembly to search rotamers for.
508    * @param potential         a {@link ffx.numerics.Potential} object.
509    * @param algorithmListener AlgorithmListener to update the GUI.
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     // Set the default 2-body Cutoff to the van der Waals cutoff.
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     // Process properties; most are seldom used and not available in ManyBodyOptions.
550     // Handle testing flags.
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     // Compute 4-body energies.
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     // Box / sliding window options.
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     // Ensemble options.
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     // Nucleic Acid Options.
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     // Superposition
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     // Multi-residue clash thresholds
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     // Monte Carlo Options.
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    * Set the "use" flag to true for all variable atoms in a residue.
688    *
689    * @param residue The residue to turn off.
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    * Set the "use" flag to true for all variable atoms in a residue.
711    *
712    * @param residue The Residue to turn on.
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    * Checks if residue i is considered to be interacting with residue j, and thus has non-null
734    * elements in the pair energies matrix.
735    *
736    * @param i A residue index.
737    * @param j A residue index j != i.
738    * @return If i and j interact.
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    * Checks if residue i is considered to be interacting with residue j, that residue k is
756    * interacting with either i or j, and thus i-j-k has non-null elements in the triple energy
757    * matrix.
758    *
759    * @param i A residue index.
760    * @param j A residue index j != i.
761    * @param k A residue index k != i, k != j.
762    * @return If i, j, and k form an interacting triple.
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    * Checks the pair elimination array to see if this permutation has been eliminated.
779    *
780    * @param i           Residue number
781    * @param ri          Rotamer number
782    * @param currentRots Array of current rotamers indeces.
783    * @return If valid
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    * Computes the environment/backbone energy, defined as energy with all sidechains under
800    * consideration turned off in their 0th rotamer.
801    *
802    * @param residues Residues under optimization.
803    * @return Backbone energy Eenv/bb.
804    * @throws ArithmeticException If an exception in calculating energy is found.
805    */
806   public double computeBackboneEnergy(Residue[] residues) throws ArithmeticException {
807     // Set all atoms to be used.
808     Atom[] atoms = molecularAssembly.getAtomArray();
809     for (Atom atom : atoms) {
810       atom.setUse(true);
811     }
812 
813     // Turn off all Residues and set them to their default conformation.
814     turnOffAllResidues(residues);
815 
816     // Compute and return the backbone energy.
817     return currentEnergy(residues);
818   }
819 
820   /**
821    * Uses existing backbone, self, 2-Body, and 3-body energies from rotamerEnergies() to calculate an
822    * approximate energy for a rotamer permutation.
823    *
824    * @param residues Current window of optimization.
825    * @param rotamers Set of rotamers to calculate an approximate energy for.
826    * @param print    Verbosity flag
827    * @return Approximate permutation energy (backbone + selfs + pairs + trimers).
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    * Calculates the energy at the current state.
867    *
868    * @param resArray Array of residues in current energy term.
869    * @return Energy of the current state.
870    */
871   public double currentEnergy(Residue[] resArray) throws ArithmeticException {
872     return currentEnergy(Arrays.asList(resArray));
873   }
874 
875   /**
876    * Wrapper intended for use with RotamerMatrixMC.
877    *
878    * @param resList Reside list.
879    * @return Returns the current energy.
880    * @throws ArithmeticException Thrown if there is an arithmetic exception computing the
881    *                             energy.
882    */
883   public double currentEnergyWrapper(List<Residue> resList) throws ArithmeticException {
884     return currentEnergy(resList);
885   }
886 
887   /**
888    * Forces the use of a ForceFieldEnergyOpenMM's underlying ForceFieldEnergy.
889    *
890    * @return Current potential energy as calculated by FFX reference platform.
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    * Utility method for formatting energies, using 16 spaces with 8 digits of precision.
903    *
904    * @param energy Energy to format.
905    * @return A string representing the energy.
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    * getEnsemble.
933    *
934    * @return a {@link java.util.List} object.
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    * setEnsemble.
948    *
949    * @param ensemble a int.
950    */
951   public void setEnsemble(int ensemble) {
952     setEnsemble(ensemble, 5.0);
953   }
954 
955   /**
956    * Return the residue list.
957    *
958    * @return a {@link java.util.List} object.
959    */
960   public List<Residue> getResidues() {
961     return residueList;
962   }
963 
964   /**
965    * Set the residue list.
966    *
967    * @param residueList a {@link java.util.List} object.
968    */
969   public void setResidues(List<Residue> residueList) {
970     this.residueList = residueList;
971   }
972 
973   /**
974    * Returns the restart file.
975    *
976    * @return energyRestartFile File with saved side-chain energies.
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       // Loop over residue k's rotamers.
994       for (int rk = 0; rk < nrk; rk++) {
995         if (eR.check(k, rk)) {
996           continue;
997         }
998         // Continue if k,rk invalid for riA/rjC.
999         if (eR.check(i, riA, k, rk) || eR.check(j, rjC, k, rk)) {
1000           // Not implemented: check(i, riA, j, rjC, k, rk).
1001           continue;
1002         }
1003         // Return false if k,rk invalid for riB/rjD.
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         // Include 3-body effects.
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           // Loop over a 4th residue l.
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             // Loop over rotamers for residue l.
1033             for (int rl = 0; rl < nrl; rl++) {
1034               // If not a part of valid phase space for riA/rjC, continue.
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                 // Not implemented: check(i, riA, j, rjC, l, rl) || check(i, riA, k, rk, l, rl) ||
1038                 // check(j, rjC, k, rk, l, rl) || check(i, riA, j, rjC, k, rk, l, rl)
1039                 continue;
1040               }
1041               if (eR.check(i, riB, l, rl) || eR.check(j, rjD, l, rl)) {
1042                 // Not implemented: check(i, riB, j, rjD, l, rl) || check(i, riB, k, rk, l, rl) ||
1043                 // check(j, rjD, k, rk, l, rl) || check(i, riB, j, rjD, k, rk, l, rl)
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    * Perform the rotamer optimization using the specified algorithm.
1090    *
1091    * @param algorithm a {@link RotamerOptimization.Algorithm} object.
1092    * @return the lowest energy found.
1093    */
1094   public double optimize(Algorithm algorithm) {
1095     this.algorithm = algorithm;
1096     return optimize();
1097   }
1098 
1099   /**
1100    * Execute the rotamer optimization.
1101    *
1102    * @return The lowest energy found.
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        * Collect all residues in the MolecularAssembly. Use all Residues with Rotamers.
1121        */
1122       allResiduesList = new ArrayList<>();
1123       // An array of polymers from the MolecularAssembly.
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       // If -DignoreNA=true, then remove nucleic acids from residue list.
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()); // for NA only
1153 
1154       nAllResidues = allResiduesList.size();
1155       allResiduesArray = allResiduesList.toArray(new Residue[nAllResidues]);
1156 
1157       /*
1158        * Distance matrix is used to add residues to the sliding window based on distance cutoff,
1159        * and for cutoff distances.
1160        *
1161        * The memory and compute overhead can be a problem for some very large structures.
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         // Allocate memory for storing optimal rotamers.
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    * A brute-force global optimization over side-chain rotamers using a recursive algorithm.
1208    *
1209    * @param molecularAssembly a {@link ffx.potential.MolecularAssembly} object.
1210    * @param residues          an array of {@link ffx.potential.bonded.Residue} objects.
1211    * @param i                 a int.
1212    * @param lowEnergy         a double.
1213    * @param optimum           an array of {@link int} objects.
1214    * @return the current energy.
1215    */
1216   public double rotamerOptimization(MolecularAssembly molecularAssembly, Residue[] residues, int i,
1217                                     double lowEnergy, int[] optimum) {
1218 
1219     // This is the initialization condition.
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       // As long as there are more residues, continue the recursion for each rotamer of the current
1232       // residue.
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        At the end of the recursion, compute the potential energy for each rotamer of the final
1252        residue. If a lower potential energy is discovered, the rotamers of each residue will be
1253        collected as the recursion returns up the chain.
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    * Sets the approximate dimensions of boxes, over-riding numXYZBoxes in determining box size.
1284    * Rounds box size up and number of boxes down to get a whole number of boxes along each axis.
1285    *
1286    * @param approxBoxLength Optional box dimensions parameter (Angstroms).
1287    */
1288   public void setApproxBoxLength(double approxBoxLength) {
1289     boxOpt.approxBoxLength = approxBoxLength;
1290   }
1291 
1292   /**
1293    * Sets the amount of overlap between adjacent boxes for box optimization.
1294    *
1295    * @param boxBorderSize Box overlap in Angstroms.
1296    */
1297   public void setBoxBorderSize(double boxBorderSize) {
1298     boxOpt.cellBorderSize = boxBorderSize;
1299   }
1300 
1301   /**
1302    * Set the ending box index.
1303    *
1304    * @param boxEnd a int.
1305    */
1306   public void setBoxEnd(int boxEnd) {
1307     // Is -1 if boxes run to completion.
1308     boxOpt.cellEnd = boxEnd;
1309   }
1310 
1311   /**
1312    * Sets behavior for how Residues are added to boxOptCells; 1 uses just reference atom (C alpha for
1313    * protein, N1/9 for nucleic acids), 2 uses any atom, 3 uses any atom in any rotamer.
1314    *
1315    * @param boxInclusionCriterion Criterion to use
1316    */
1317   public void setBoxInclusionCriterion(int boxInclusionCriterion) {
1318     boxOpt.boxInclusionCriterion = boxInclusionCriterion;
1319   }
1320 
1321   /**
1322    * Set the starting box index.
1323    *
1324    * @param boxStart a int.
1325    */
1326   public void setBoxStart(int boxStart) {
1327     boxOpt.cellStart = boxStart;
1328   }
1329 
1330   /**
1331    * setCoordinatesToEnsemble.
1332    *
1333    * @param ensnum a int.
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    * Sets the decompose-original flag.
1346    *
1347    * @param decomposeOriginal If true, decompose the energy of the structure without optimizing.
1348    */
1349   public void setDecomposeOriginal(boolean decomposeOriginal) {
1350     this.decomposeOriginal = decomposeOriginal;
1351   }
1352 
1353   /**
1354    * Set the optimization direction to forward or backward.
1355    *
1356    * @param direction a {@link RotamerOptimization.Direction} object.
1357    */
1358   public void setDirection(Direction direction) {
1359     this.direction = direction;
1360   }
1361 
1362   /**
1363    * Set the cut-off distance for inclusion of residues in sliding box and window methods.
1364    *
1365    * @param distance a double.
1366    */
1367   public void setDistanceCutoff(double distance) {
1368     this.distance = distance;
1369   }
1370 
1371   /**
1372    * Setter for the field <code>energyRestartFile</code>.
1373    *
1374    * @param file a {@link java.io.File} object.
1375    */
1376   public void setEnergyRestartFile(File file) {
1377     loadEnergyRestart = true;
1378     energyRestartFile = file;
1379   }
1380 
1381   /**
1382    * setEnsemble.
1383    *
1384    * @param ensemble       a int.
1385    * @param ensembleBuffer a double.
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    * Set the residue increment for sliding window.
1398    *
1399    * @param increment a int.
1400    */
1401   public void setIncrement(int increment) {
1402     this.increment = increment;
1403   }
1404 
1405   /**
1406    * Control the depth of self-consistency checking with a rotamer is eliminated.
1407    *
1408    * @param maxRotCheckDepth a int.
1409    */
1410   public void setMaxRotCheckDepth(int maxRotCheckDepth) {
1411     this.maxRotCheckDepth = maxRotCheckDepth;
1412   }
1413 
1414   /**
1415    * Set the minimum number of accepted nucleic acid rotamers.
1416    *
1417    * @param minNumberAcceptedNARotamers a int.
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    * Sets the option to use a number of Monte Carlo steps for final optimization.
1431    *
1432    * @param monteCarlo If Monte Carlo is to be considered
1433    * @param nMCsteps   Number of steps to be taken
1434    */
1435   public void setMonteCarlo(boolean monteCarlo, int nMCsteps) {
1436     this.monteCarlo = monteCarlo;
1437     this.nMCSteps = nMCsteps;
1438   }
1439 
1440   /**
1441    * Sets the monteCarloTesting boolean in RotamerOptimization.java to true or false. This should
1442    * only be set to true when monte carlo is being tested through the ManyBodyTest script. When true,
1443    * the method sets a seed for the pseudo-random number generator and allows the monte carlo rotamer
1444    * optimization to be deterministic.
1445    *
1446    * @param bool True or false.
1447    */
1448   public void setMonteCarloTesting(boolean bool) {
1449     this.monteCarloTesting = bool;
1450   }
1451 
1452   /**
1453    * The nucleic acid correction threshold.
1454    *
1455    * @param nucleicCorrectionThreshold a double.
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    * Also sets derivative pruning factors.
1469    *
1470    * @param nucleicPruningFactor a double.
1471    */
1472   public void setNucleicPruningFactor(double nucleicPruningFactor) {
1473     this.nucleicPruningFactor = nucleicPruningFactor;
1474     this.nucleicPairsPruningFactor = ((1.0 + nucleicPruningFactor) / 2);
1475   }
1476 
1477   /**
1478    * Sets the number of boxes in the x, y, and z axes if the box optimization is to be carried out.
1479    *
1480    * @param numXYZBoxes Int[3] of number of boxes in x, y, z.
1481    */
1482   public void setNumXYZBoxes(int[] numXYZBoxes) {
1483     arraycopy(numXYZBoxes, 0, boxOpt.numXYZCells, 0, boxOpt.numXYZCells.length);
1484   }
1485 
1486   /**
1487    * Setter for the field <code>pairClashThreshold</code>.
1488    *
1489    * @param pairClashThreshold a double.
1490    */
1491   public void setPairClashThreshold(double pairClashThreshold) {
1492     this.pairClashThreshold = pairClashThreshold;
1493   }
1494 
1495   /**
1496    * Sets whether rotamer optimization should print out any files, or act solely to optimize a
1497    * structure in memory.
1498    *
1499    * @param printFiles a boolean.
1500    */
1501   public void setPrintFiles(boolean printFiles) {
1502     this.printFiles = printFiles;
1503   }
1504 
1505   /**
1506    * Sets level of pruning: 0 for fully off, 1 for only singles, 2 for single and pair pruning.
1507    *
1508    * @param set Pruning option.
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    * Accepts a list of residues but throws out null residues. Used by the -lR flag.
1525    *
1526    * @param residues a {@link java.util.List} object.
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    * Set the algorithm to revert to starting coordinates if the energy increases.
1543    *
1544    * @param revert a boolean.
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    * setSingletonClashThreshold.
1556    *
1557    * @param singletonClashThreshold a double.
1558    */
1559   public void setSingletonClashThreshold(double singletonClashThreshold) {
1560     this.clashThreshold = singletonClashThreshold;
1561   }
1562 
1563   /**
1564    * Setter for the field <code>superpositionThreshold</code>.
1565    *
1566    * @param superpositionThreshold a double.
1567    */
1568   public void setSuperpositionThreshold(double superpositionThreshold) {
1569     this.superpositionThreshold = superpositionThreshold;
1570   }
1571 
1572   /**
1573    * Sets the threeBodyCutoffDist. All three-body energies where the rotamers have a separation
1574    * distance larger than the cutoff are set to 0.
1575    *
1576    * @param threeBodyCutoffDist Separation distance at which the interaction of three side-chains
1577    *                            is assumed to have an energy of 0.
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    * Flag to control use of 3-body energy terms.
1589    *
1590    * @param threeBodyTerm a boolean.
1591    */
1592   public void setThreeBodyEnergy(boolean threeBodyTerm) {
1593     this.threeBodyTerm = threeBodyTerm;
1594   }
1595 
1596   /**
1597    * Sets the twoBodyCutoffDist. All two-body energies where the rotamers have a separation distance
1598    * larger than the cutoff are set to 0.
1599    *
1600    * @param twoBodyCutoffDist Separation distance at which the interaction of two side-chains is
1601    *                          assumed to have an energy of 0.
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    * Specify use of Goldstein optimization.
1612    *
1613    * @param useGoldstein a boolean.
1614    */
1615   public void setUseGoldstein(boolean useGoldstein) {
1616     this.useGoldstein = useGoldstein;
1617   }
1618 
1619   /**
1620    * Setter for the field <code>windowSize</code>.
1621    *
1622    * @param windowSize a int.
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    * {@inheritDoc}
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    * {@inheritDoc}
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    * ONLY FOR UNIT TESTING. Sets a boolean to turn the pair elimination criteria off.
1692    */
1693   public void turnRotamerPairEliminationOff() {
1694     logger.info(" Turning off pair eliminations.");
1695     pairEliminationOn = false;
1696   }
1697 
1698   /**
1699    * ONLY FOR UNIT TESTING. Sets a boolean to turn the self elimination criteria off.
1700    */
1701   public void turnRotamerSingleEliminationOff() {
1702     logger.info(" Turning off single eliminations.");
1703     selfEliminationOn = false;
1704   }
1705 
1706   /**
1707    * Recursive brute-force method which uses single, pair, and potentially trimer energies to
1708    * calculate an optimum set of rotamers.
1709    *
1710    * @param residues        Optimization window
1711    * @param i               Current residue in the recursion.
1712    * @param lowEnergy       Minimum energy yet found by the recursion.
1713    * @param optimum         Optimum rotamer set yet found by the recursion.
1714    * @param currentRotamers Rotamer permutation under investigation.
1715    * @return Minimum energy found under this node in the recursion.
1716    */
1717   private double decomposedRotamerOptimization(Residue[] residues, int i, double lowEnergy,
1718                                                int[] optimum, int[] currentRotamers) {
1719 
1720     // This is the initialization condition.
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       // As long as there are more residues, continue the recursion for each rotamer of the current
1732       // residue.
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       // At the end of the recursion, compute the potential energy for each rotamer of the final
1749       // residue and update optimum[].
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            Because we print the rotamer set immediately on finding a
1763            more optimal structure, we have to reset the entire length
1764            of optimum instead of lazily doing it on the way down.
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    * Runs Monte Carlo side chain optimization using the rotamer energy matrix and potentially some
1788    * information from dead-end or Goldstein elimination. The useAllElims variable should be set false
1789    * if detailed balance is to be maintained. At present, no support for ensembles.
1790    *
1791    * @param residues      Optimization window
1792    * @param optimum       Array to store optimum rotamers
1793    * @param initialRots   Array with starting rotamers
1794    * @param maxIters      Number of MC steps to run
1795    * @param randomizeRots Scramble initialRots
1796    * @param useAllElims   Use pair/triple elimination information
1797    * @return Lowest energy found
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    * Scrambles an array of rotamers.
1854    *
1855    * @param rotamers    Array of rotamers.
1856    * @param residues    Array of residues.
1857    * @param useAllElims If trye, use all elliminations.
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    * A global optimization over side-chain rotamers using a recursive algorithm and information about
1894    * eliminated rotamers, rotamer pairs and rotamer triples
1895    *
1896    * @param molecularAssembly   MolecularAssmebly to use.
1897    * @param residues            Array of residues.
1898    * @param i                   Current permutation.
1899    * @param currentRotamers     Current array of rotamers.
1900    * @param lowEnergy           Lowest energy found.
1901    * @param optimum             Optimum set of rotamers.
1902    * @param permutationEnergies Energies of visited permutations or null.
1903    * @return current energy.
1904    */
1905   private double rotamerOptimizationDEE(MolecularAssembly molecularAssembly, Residue[] residues,
1906                                         int i, int[] currentRotamers, double lowEnergy, int[] optimum, double[] permutationEnergies) {
1907     // This is the initialization condition.
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     // As long as there are more residues, continue the recursion for each rotamer of the current
1920     // residue.
1921     if (i < nResidues - 1) {
1922       // Loop over rotamers of residue i.
1923       for (int ri = 0; ri < lenri; ri++) {
1924         // Check if rotamer ri has been eliminated by DEE.
1925         if (eR.check(i, ri)) {
1926           continue;
1927         }
1928         // Check if rotamer ri has been eliminated by an upstream rotamer (any residue's rotamer
1929         // from j = 0 .. i-1).
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        At the end of the recursion, compute the potential energy for each rotamer of the final
1960        residue. If a lower potential energy is discovered, the rotamers of each residue will be
1961        collected as the recursion returns up the chain.
1962       */
1963       for (int ri = 0; ri < lenri; ri++) {
1964         // Check if rotamer ri has been eliminated by DEE.
1965         if (eR.check(i, ri)) {
1966           continue;
1967         }
1968         currentRotamers[i] = ri;
1969         // Check if rotamer ri has been eliminated by an upstream rotamer (any residue's rotamer
1970         // from 0 .. i-1).
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         // Compute the energy based on a 3-body approximation
1984         approximateEnergy = computeEnergy(residues, currentRotamers, false);
1985         double comparisonEnergy = approximateEnergy;
1986         evaluatedPermutations++;
1987         // Compute the AMOEBA energy
1988         if (useForceFieldEnergy) {
1989           double amoebaEnergy = Double.NaN;
1990           try {
1991             // Add the rotamer pH bias to the force field energy.
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           // Log current results
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    * A global optimization over side-chain rotamers using a recursive algorithm and information about
2060    * eliminated rotamers, rotamer pairs and rotamer triples.
2061    *
2062    * @param residues        Residue array.
2063    * @param i               Current number of permutations.
2064    * @param currentRotamers Current rotamer list.
2065    * @return 0.
2066    */
2067   private double dryRun(Residue[] residues, int i, int[] currentRotamers) {
2068     // This is the initialization condition.
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       // At the end of the recursion, check each rotamer of the final residue.
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    * Finds all permutations within buffer energy of GMEC.
2131    *
2132    * @param residues            Array of residues.
2133    * @param i                   Current depth in residue/rotamer tree.
2134    * @param currentRotamers     Current set of rotamers at this node.
2135    * @param gmecEnergy          Minimum energy for these residues.
2136    * @param permutationEnergies Energy of all permutations.
2137    * @param permutations        Contains accepted permutations.
2138    * @return 0.
2139    */
2140   private double dryRunForEnsemble(Residue[] residues, int i, int[] currentRotamers,
2141                                    double gmecEnergy, double[] permutationEnergies, int[][] permutations) {
2142     // This is the initialization condition.
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       // At the end of the recursion, check each rotamer of the final residue.
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    * Return an integer array of optimized rotamers following rotamer optimization.
2200    *
2201    * @return The optimal rotamer array.
2202    */
2203   public int[] getOptimumRotamers() {
2204     return optimum;
2205   }
2206 
2207   /**
2208    * Independent optimization treats each residue sequentially.
2209    *
2210    * @param residues The list of residues to be optimized.
2211    * @return The final energy.
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           // Report energies relative to the first rotamer.
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    * Finds the first non-eliminated rotamer permutation.
2277    *
2278    * @param residues        Array of residues.
2279    * @param i               Current permutation.
2280    * @param currentRotamers Current list of rotamers.
2281    * @return If valid permutation found.
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       // At the end of the recursion, check each rotamer of the final residue.
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    * The main driver for optimizing a block of residues using DEE.
2335    *
2336    * @param residueList Residues to optimize.
2337    * @return Final energy.
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       // Compute the number of permutations without eliminating dead-ends and compute the number of
2356       // permutations using singleton elimination.
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         // Not currently compatible with ensembles.
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         // Add the force field energy to the pH bias.
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     // Permutations used only to set maximum bound on ensembleNumber, thus it is safe here to put
2568     // that value in a 32-bit int.
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       // Compute the number of permutations without eliminating dead-ends and compute the number of
2601       // permutations using singleton elimination.
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           // TODO: Decide whether we like these next four lines.  Has the potential to produce a
2689           // crazy amount of permutations.
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       // Add the force field energy to the pH bias.
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    * Use Monte Carlo if monteCarlo specified, and either skipDEE specified or nMCsteps is smaller
2846    * than the remaining permutation size.
2847    *
2848    * @return Finish DEE search with Monte Carlo.
2849    */
2850   private boolean useMonteCarlo() {
2851     return monteCarlo && (mcNoEnum || (nMCSteps < evaluatedPermutations));
2852   }
2853 
2854   /**
2855    * Performs a recursive brute-force rotamer optimization over a passed list of residues.
2856    *
2857    * @param residueList Residues to be optimized.
2858    * @return Global minimum energy conformation energy.
2859    */
2860   private double bruteForce(List<Residue> residueList) {
2861 
2862     Residue[] residues = residueList.toArray(new Residue[0]);
2863     int nResidues = residues.length;
2864 
2865     // Compute the number of permutations without eliminating dead-ends.
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       // Use a many-body energy sum to evaluate permutations.
2879       setPruning(0);
2880       rotamerEnergies(residues);
2881       int[] currentRotamers = new int[nResidues];
2882       e = decomposedRotamerOptimization(residues, 0, Double.MAX_VALUE, optimum, currentRotamers);
2883 
2884       // Set each residue to their optimal rotamer.
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       // Compute the full energy for comparison (i.e., without many-body truncation).
2893       double fullEnergy = 0;
2894       try {
2895         // Add the force field energy to the pH bias.
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         // Fall through into the FORWARD case.
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            If the window starts with a nucleic acid, and there is a 5' NA residue, ensure that 5'
3007            NA residue has been included in the window. Otherwise, that previous residue may not
3008            have had a chance to be flexible about its sugar pucker.
3009            <p>
3010            If window size is greater than increment, however, this has already been handled.
3011            Additionally, do not perform this for the first window (counter is already incremented
3012            by the time this check is performed, so first window's counter will be 2). Furthermore,
3013            do not include Residues with null Rotamer lists (this breaks things).
3014            <p>
3015            The issue: if window size = increment, the last NA residue in each window will not
3016            have flexibility about its sugar pucker, because its self-energy includes the O3' (i)
3017            to P (i+1) bond, so it must remain in the original sugar pucker to meet the i+1
3018            residue. However, this problem can be solved by ensuring that final residue is included
3019            in the next window, where it will have flexibility about its sugar pucker.
3020            <p>
3021            If you are running successive sliding window jobs on the same file, I would suggest
3022            starting the next job on the last residue of the previous job, unless you know your
3023            settings will include it.
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               // Copy sliding window optimal rotamers into the overall optimum array.
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             // Copy sliding window optimal rotamers into the overall optimum array.
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             // Don't write a file if it's the final iteration.
3095             if (windowStart + windowSize == nOptimize) {
3096               continue;
3097             }
3098             PDBFilter windowFilter = new PDBFilter(file, molecularAssembly, null, null);
3099             //   StringBuilder header = new StringBuilder(format("Iteration %d of the sliding
3100             // window\n", counter - 1));
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         // No default case.
3125         break;
3126     }
3127     return 0.0;
3128   }
3129 
3130   /**
3131    * Sorts a passed List of Residues by global index.
3132    *
3133    * @param residues List of Residues to be sorted.
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    * applyEliminationCriteria.
3143    *
3144    * @param residues    an array of {@link ffx.potential.bonded.Residue} objects.
3145    * @param getEnergies a boolean.
3146    * @param verbose     a boolean.
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             // While there are eliminated rotamers, repeatedly apply single rotamer elimination.
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             // While there are eliminated rotamer pairs, repeatedly apply rotamer pair elimination.
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             // While there are eliminated rotamers, repeatedly apply single rotamer elimination.
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             // While there are eliminated rotamer pairs, repeatedly apply rotamer pair elimination.
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; // Length rotamers i
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         // Loop over the set of rotamers for residue i.
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           // While there are eliminated rotamers, repeatedly apply single rotamer elimination.
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           // While there are eliminated rotamer pairs, repeatedly apply rotamer pair elimination.
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           // While there are eliminated rotamers, repeatedly apply single rotamer elimination.
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           // While there are eliminated rotamer pairs, repeatedly apply rotamer pair elimination.
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    * Calculates the energy at the current state.
3330    *
3331    * @param resList List of residues in current energy term.
3332    * @return Energy of the current state.
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    * Default method for obtaining energy: calculates energy of the Potential.
3343    *
3344    * @param dir Ignored, should be null
3345    * @return Current potential energy
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    * Generates list of subsequent, neighboring Residues.
3358    *
3359    * @param residues To generate neighbor lists for.
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       // Collect all neighbors.
3399       int[] nI = nearby.stream().mapToInt(Integer::intValue).toArray();
3400       bidiResNeighbors[i] = nI;
3401 
3402       // Collect only subsequent neighbors.
3403       final int fi = i; // Final copy of 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    * Turn off non-bonded contributions from all residues except for one. Compute the self-energy for
3430    * each residue relative to the backbone contribution.
3431    *
3432    * @param residues A list of residues that will undergo rotamer optimization.
3433    * @return template energy
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         // Leave the 0'th original-coordinates rotamer alone.
3458         for (int ri = 1; ri < lenri; ri++) {
3459           eR.eliminateRotamer(residues, i, ri, false);
3460         }
3461       }
3462     }
3463 
3464     // Initialize all atoms to be used.
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     // Update the EliminatedRotamers instance with the EnergyExpansion instance.
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     // Turn on all atoms.
3552     for (Atom atom : atoms) {
3553       atom.setUse(true);
3554     }
3555     // Print the energy with all rotamers in their default conformation.
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    * Applies the "default" rotamer: currently the 0'th rotamer.
3570    *
3571    * @param residue Residue to apply a default rotamer for.
3572    */
3573   private void applyDefaultRotamer(Residue residue) {
3574     RotamerLibrary.applyRotamer(residue, residue.getRotamers()[0]);
3575   }
3576 
3577   /**
3578    * Elimination of rotamers via the original Dead End Elimination algorithm.
3579    *
3580    * @param residues Array of residues under consideration.
3581    * @return True if any rotamers were eliminated.
3582    */
3583   private boolean deeRotamerElimination(Residue[] residues) {
3584     int nres = residues.length;
3585     // A flag to indicate if more rotamers or rotamer pairs were eliminated.
3586     boolean eliminated = false;
3587     // Loop over residues.
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       // Loop over the set of rotamers for residue i.
3600       for (int ri = 0; ri < lenri; ri++) {
3601         // Check for an eliminated single.
3602         if (eR.check(i, ri)) {
3603           continue;
3604         }
3605         // Start the min/max summation with the self-energy.
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               // There is a chance that ri conflicts with every possible rotamer of some
3617               // residue j.
3618               // In that event, its minimum energy is set NaN and should be easy to eliminate.
3619               minEnergySingles[ri] = Double.NaN;
3620             }
3621             if (Double.isFinite(minMax[0]) && Double.isFinite(maxEnergySingles[ri])) {
3622               maxEnergySingles[ri] += minMax[1];
3623             } else {
3624               // In this branch, ri conflicts with some j,rj and cannot be practically used for
3625               // elimination.
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             // eliminateRotamer(residues, i, ri, print);
3634           }
3635         }
3636       }
3637 
3638       /*
3639        Apply the singles elimination criteria to rotamers of residue i by determining the most
3640        favorable maximum energy.
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         // This branch is taken if every ri conflicts with at least one j,rj. In that case, nothing
3656         // can be eliminated yet!
3657         logIfRank0(" Could not eliminate any i,ri because eliminationEnergy was never set!",
3658             Level.FINE);
3659       } else {
3660         /*
3661          Eliminate rotamers whose minimum energy is greater than the worst case for another
3662          rotamer.
3663         */
3664         for (int ri = 0; ri < lenri; ri++) {
3665           if (eR.check(i, ri)) {
3666             continue;
3667           }
3668           // If i,ri has a clash with all phase space, it can be eliminated by something that
3669           // doesn't clash with all phase space.
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           // Otherwise, can eliminate if its best possible energy is still worse than something
3678           // else's worst possible energy.
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    * Rotamer pair elimination driver for many-body Dead End Elimination. Generally less effective
3696    * than Goldstein.
3697    *
3698    * @param residues Residues under consideration.
3699    * @return If at least one pair eliminated.
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       // Minimum and maximum summation found for ri-rj pairs.
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             // Min and max external summations for ri-rj.
3735             double[] minMax = new double[2];
3736 
3737             // Add contributions from third residues k, and possibly fourth residues l.
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                 // ri-rj can clash, and isn't very useful to eliminate by.
3751                 maxPairEnergies[ri][rj] = Double.NaN;
3752               }
3753             } else {
3754               // A NaN minimum energy for some pair indicates it's definitely not part of the GMEC.
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                   // See above check(i, ri, j, rj) for why this should not be taken!
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     // A flag to indicate if a rotamer is eliminated.
3829     boolean eliminated = false;
3830     // Loop over residue i.
3831     for (int i = 0; i < nres; i++) {
3832       Residue resi = residues[i];
3833       Rotamer[] roti = resi.getRotamers();
3834       int nri = roti.length;
3835       // Loop over the set of rotamers for residue i.
3836       for (int riA = 0; riA < nri; riA++) {
3837         // Continue if rotamer (i, riA) is not valid.
3838         if (eR.check(i, riA)) {
3839           continue;
3840         }
3841         for (int riB = 0; riB < nri; riB++) {
3842           // The eliminating rotamer cannot be riA and must be a valid.
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    * Attemps to eliminate rotamer riA based on riB.
3861    *
3862    * @param residues Array of Residues.
3863    * @param i        Residue i index.
3864    * @param riA      Rotamer to attempt elimination of.
3865    * @param riB      Rotamer to attempt elimination by.
3866    * @return If riA was eliminated.
3867    */
3868   private boolean goldsteinElimination(Residue[] residues, int i, int riA, int riB) {
3869     Residue resi = residues[i];
3870 
3871     // Initialize Goldstein inequality.
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     // Loop over a 2nd residue j.
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       // Loop over the rotamers for residue j.
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; // This is not a part of configuration space accessible to riA.
3896         }
3897         if (eR.check(i, riB, j, rj)) {
3898           /*
3899            This is a part of configuration space where riA is valid but not riB. Thus, if j,rj is
3900            part of the GMEC, riB is inconsistent with it. Thus, riB cannot be used to eliminate
3901            riA.
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         // Include three-body interactions.
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                If k,rk or j,rj-k,rk are not a part of valid configuration space, continue. If
3930                i,riA-k,rk or i,riA-j,rj-k,rk are not valid for riA, continue.
3931               */
3932               if (eR.check(k, rk) || eR.check(j, rj, k, rk) || eR.check(i, riA, k, rk)) {
3933                 // Not yet implemented: check(i, riA, j, rj, k, rk) because no triples get
3934                 // eliminated.
3935                 continue;
3936               }
3937               /*
3938                If i,riB-k,rk or i,riB-j,rj-k,rk are invalid for riB, there is some part of
3939                configuration space for which riA is valid but not riB.
3940               */
3941               if (eR.check(i, riB, k, rk)) {
3942                 // Not yet implemented: check(i, riB, j, rj, k, rk).
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             /* If there were no 3-body interactions with residue k, then minForResk is zero. */
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       // If there are no 2-body interactions, then minForResJ is zero.
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    * Finds and eliminates rotamer pairs according to the many-body Goldstein pairs criterion.
3995    *
3996    * @param residues Residues under consideration.
3997    * @return If any rotamer pairs were eliminated.
3998    */
3999   private boolean goldsteinPairDriver(Residue[] residues) {
4000     int nRes = residues.length;
4001     boolean eliminated = false;
4002 
4003     // First, generate pairs riA-rjC.
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         // Don't try to eliminate that which is already eliminated.
4011         if (eR.check(i, riA)) {
4012           continue;
4013         }
4014 
4015         // Residue j can be any other residue, including ones before residue i.
4016         for (int j = 0; j < nRes; j++) {
4017           // Residue j must be distinct from i.
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             // Again, no point in eliminating the already-eliminated.
4026             if (eR.check(j, rjC) || eR.check(i, riA, j, rjC)) {
4027               continue;
4028             }
4029             boolean breakOut = false;
4030 
4031             // Now, generate pairs riB-rjD. If any pair riB-rjD eliminates riA-rjC, break out of the
4032             // loop.
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                 // Do not attempt eliminating with an eliminated pair.
4045                 if (eR.check(j, rjD) || eR.check(i, riB, j, rjD)) {
4046                   continue;
4047                 }
4048                 // Do not attempt to eliminate a pair with itself.
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    * Attempt to eliminate rotamer pair (ResI-RotA, ResJ-RotC) using (ResI-RotB, ResJ-RotD).
4070    *
4071    * @param residues Array of residues.
4072    * @param i        Index of the first residue.
4073    * @param riA      Index of the first residue's rotamer to eliminate.
4074    * @param riB      Index of the first residue's rotamer to use for elimination.
4075    * @param j        Index of the 2nd residue.
4076    * @param rjC      Index of the 2nd residue's rotamer to eliminate.
4077    * @param rjD      Index of the 2nd residue's rotamer to use for elimination.
4078    * @return Return true if eliminated.
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     // Initialize the Goldstein energy.
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    * Method allows for testing of the elimination criteria by setting parameters appropriately.
4141    *
4142    * @param testing True only during RotamerOptimizationTest.
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    * Method allows for testing of the elimination criteria by setting parameters appropriately for a
4153    * specific test case.
4154    *
4155    * @param testSelfEnergyEliminations True when only self energies are calculated; pairs,
4156    *                                   triples, etc., are assumed to be 0.
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    * Method allows for testing of the elimination criteria by setting parameters appropriately for a
4171    * specific test case.
4172    *
4173    * @param testPairEnergyEliminations True when only pair energies are calculated; selves,
4174    *                                   triples, etc., are assumed to be 0.
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    * Method allows for testing of the elimination criteria by setting parameters appropriately for a
4189    * specific test case.
4190    *
4191    * @param testTripleEnergyEliminations1 True when only triple energies are calculated; selves,
4192    *                                      pairs, etc., are assumed to be 0.
4193    * @param testTripleEnergyEliminations2 True when only triple energies are calculated; selves,
4194    *                                      pairs, etc., are assumed to be 0.
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    * Test the self-energy elimination by setting 2-body and 3-body interactions to zero.
4210    *
4211    * @param residues an array of {@link ffx.potential.bonded.Residue} objects.
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               // catch NPE.
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                     // catch NPE.
4240                   }
4241                 }
4242               }
4243             }
4244           }
4245         }
4246       }
4247     }
4248   }
4249 
4250   /**
4251    * Test the elimination criteria by setting self and 3-body interactions to zero.
4252    *
4253    * @param residues an array of {@link ffx.potential.bonded.Residue} objects.
4254    * @param resID    The residue to test.
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           // catch NPE.
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                 // catch NPE.
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                     // catch NPE.
4295                   }
4296                 }
4297               }
4298             }
4299           }
4300         }
4301       }
4302     }
4303   }
4304 
4305   /**
4306    * Test the elimination criteria by setting self and 2-body interactions to zero. Two residues are
4307    * at fixed rotamers and all rotamer interactions with those two residues are calculated.
4308    *
4309    * @param residues an array of {@link ffx.potential.bonded.Residue} objects.
4310    * @param resID1   The residue number for one of two fixed residues.
4311    * @param resID2   The second residue number for one of two fixed residues.
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           // catch NPE.
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              if (i != resID1 && j != resID1) { try { twoBodyEnergy[i][ri][j][rj] = 0.0; } catch
4343              (Exception e) { // catch NPE. } }
4344             */
4345             try {
4346               eE.set2Body(i, ri, j, rj, 0, true);
4347             } catch (Exception e) {
4348               // catch NPE.
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                       // catch NPE.
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                       // catch NPE.
4369                     }
4370                   }
4371                 }
4372               }
4373             }
4374           }
4375         }
4376       }
4377     }
4378   }
4379 
4380   /**
4381    * Rotamer Optimization Methods.
4382    */
4383   public enum Algorithm {
4384     INDEPENDENT,  // 1
4385     ALL,          // 2
4386     BRUTE_FORCE,  // 3
4387     WINDOW,       // 4
4388     BOX;          // 5
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    * Allows get2BodyDistance to find the shortest distance between any two rotamers or two residues.
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      * Number of boxes for box optimization in X, Y, Z.
4424      */
4425     public final int[] numXYZCells = {3, 3, 3};
4426     /**
4427      * Box border size.
4428      */
4429     public double cellBorderSize = 0;
4430     /**
4431      * Approximate box size.
4432      */
4433     public double approxBoxLength = 0;
4434     /**
4435      * Box optimization inclusion criteria.
4436      */
4437     public int boxInclusionCriterion = 1;
4438     /**
4439      * Index of the first box to optimize.
4440      */
4441     public int cellStart = 0;
4442     /**
4443      * Index of the last box to optimize.
4444      */
4445     public int cellEnd = -1;
4446     /**
4447      * Flag to indicate manual definition of a super box.
4448      */
4449     public boolean manualSuperbox = false;
4450     /**
4451      * Dimensions of the box.
4452      */
4453     public double[] boxDimensions;
4454     /**
4455      * Buffer size for the super box.
4456      */
4457     public double superboxBuffer = 8.0;
4458     /**
4459      * Box index loaded during a restart.
4460      */
4461     public int boxLoadIndex = -1;
4462     /**
4463      * Box indeces loaded during a restart.
4464      */
4465     public int[] boxLoadCellIndices;
4466 
4467     /**
4468      * Breaks down a structure into a number of overlapping boxes for optimization.
4469      *
4470      * @return Potential energy of final structure.
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        * A new dummy Crystal will be constructed for an aperiodic system. The
4479        * purpose is to avoid using the overly large dummy Crystal used for
4480        * Ewald purposes. Atoms are not and should not be moved into the dummy
4481        * Crystal boundaries; to check if an Atom is inside a cell, use an
4482        * array of coordinates adjusted to be 0 < coordinate < 1.0.
4483        */
4484       Crystal crystal = generateSuperbox(residueList);
4485 
4486       // Cells indexed by x*(YZ divisions) + y*(Z divisions) + z.
4487       int totalCells = getTotalCellCount(crystal); // Also initializes cell count if using -bB
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         // Warn the user if the box end was explicitly set incorrectly.
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               // Copy sliding window optimal rotamers into the overall optimum array.
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             // Copy sliding window optimal rotamers into the overall optimum array.
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             // Don't write a file if it's the final iteration.
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       // String should be in format (buffer,xmin,xmax,ymin,ymax,zmin,zmax)
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 { // Crystal systems will have already returned.
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      * Returns the number of cells (boxes) for box optimization; if the -bB flag is set, sets the
4684      * final number of cells.
4685      *
4686      * @param crystal Crystal or dummy crystal being used to define boxes.
4687      * @return Total number of cells.
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      * Creates and fills cells (boxes) for box optimization.
4712      *
4713      * @param crystal  Polymer crystal or dummy crystal
4714      * @param residues All residues to be optimized
4715      * @return Filled cells.
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; // Breaks the ijk loop if the last box passed.
4730       /*
4731        * Initializes coordinates for all the boxes, indexed linearly along z,
4732        * then y, then x (so the box with xyz indices 2,3,2 in a crystal with
4733        * 4, 5, and 3 boxes along xyz would be indexed 2*5*3 + 3*3 + 2 = 41).
4734        * The int[] indices stores separate x, y, and z indices.
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           // Fall through into forward case (for now).
4779         case FORWARD:
4780         default:
4781           return cells;
4782       }
4783     }
4784 
4785     /**
4786      * Constructs the cells for box optimization and assigns them residues, presently based on C
4787      * alpha fractional coordinates; by default, cells are sorted by global index. Presently,
4788      * specifying approxBoxLength over-rides numXYZBoxes, and always rounds the number of boxes down
4789      * (to ensure boxes are always at least the specified size).
4790      *
4791      * @param crystal  Crystal group.
4792      * @param residues List of residues to be optimized.
4793      * @param cells    BoxOptCell instance.
4794      */
4795     private void assignResiduesToCells(Crystal crystal, Residue[] residues, ManyBodyCell[] cells) {
4796       // Search through residues, add them to all boxes containing their
4797       // fractional coordinates.
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           // If the identity symop produces nothing, skip checking other symops.
4816           if (toAdd.isEmpty()) {
4817             break;
4818           }
4819         }
4820         toAdd.forEach(cell::addResidue);
4821       }
4822     }
4823   }
4824 
4825   //    /**
4826   //     * Writes eliminated singles and pairs to a CSV file. Reads in .log file.
4827   //     *
4828   //     * @throws java.io.FileNotFoundException if any.
4829   //     * @throws java.io.IOException if any.
4830   //     */
4831   //    public void outputEliminated() throws FileNotFoundException, IOException {
4832   //        /**
4833   //         * eliminated.csv stores the eliminated singles and pairs.
4834   //         */
4835   //        File fileName = new File("eliminated.csv");
4836   //        String path = fileName.getCanonicalPath();
4837   //        File outputTxt = new File(path);
4838   //        Path currentRelativePath = Paths.get("eliminated.csv");
4839   //        String logPath = currentRelativePath.toAbsolutePath().toString();
4840   //        logPath = logPath.replaceAll("/eliminated.csv", "");
4841   //        File logDirectory = new File(logPath);
4842   //        /**
4843   //         * Searches for .log file within the directory.
4844   //         */
4845   //        File[] logArray = logDirectory.listFiles(new FilenameFilter() {
4846   //            public boolean accept(File dir, String filename) {
4847   //                return filename.endsWith(".log");
4848   //            }
4849   //        });
4850   //        File logFile = logArray[0];
4851   //        PrintWriter out = new PrintWriter(new FileWriter(outputTxt, true));
4852   //        BufferedReader br = new BufferedReader(new FileReader(logFile));
4853   //        String line;
4854   //        String line3;
4855   //        List<String> boxes = new ArrayList<String>();
4856   //        while ((line3 = br.readLine()) != null) {
4857   //            if (line3.contains("-a, 5")) {
4858   //                /**
4859   //                 * Stores the number of box optimization iterations and the
4860   //                 * coordinates of each box.
4861   //                 */
4862   //                while ((line = br.readLine()) != null) {
4863   //                    if (line.contains("xyz indices")) {
4864   //                        String coordinates = line.replaceAll(" Cell xyz indices:", "");
4865   //                        String nextLine = br.readLine();
4866   //                        if (!nextLine.contains("Empty box")) {
4867   //                            boxes.add(coordinates);
4868   //                            String nextLine2;
4869   //                            String nextLine4;
4870   //                            String nextLine5;
4871   //                            while (!(nextLine2 = br.readLine()).contains("Collecting
4872   // Permutations")) {
4873   //                                if (nextLine2.contains("Applying Rotamer Pair")) {
4874   //                                    int total = 0;
4875   //                                    int total2 = 0;
4876   //                                    nextLine2 = br.readLine();
4877   //                                    nextLine4 = br.readLine();
4878   //                                    nextLine5 = br.readLine();
4879   //                                    if (nextLine5.contains("Self-consistent")) {
4880   //                                        String[] split = nextLine2.split(" ");
4881   //                                        int singlesAmount = Integer.parseInt(split[1]);
4882   //                                        if (singlesIterations.size() > 0) {
4883   //                                            total =
4884   // singlesIterations.get(singlesIterations.size() - 1);
4885   //                                        }
4886   //                                        singlesIterations.add(singlesAmount + total);
4887   //                                        System.out.println(nextLine4);
4888   //                                        String[] split2 = nextLine4.split(" ");
4889   //                                        int pairsAmount = Integer.parseInt(split2[1]);
4890   //                                        if (pairsIterations.size() > 0) {
4891   //                                            total2 = pairsIterations.get(pairsIterations.size()
4892   // - 1);
4893   //                                        }
4894   //                                        pairsIterations.add(pairsAmount + total2);
4895   //                                    }
4896   //                                }
4897   //                            }
4898   //                        }
4899   //                    }
4900   //                }
4901   //                /**
4902   //                 * Writes eliminated singles and pairs to eliminated.csv, sorted
4903   //                 * by box iteration.
4904   //                 */
4905   //                out.append("Eliminated Singles\n");
4906   //                out.append("Residue,Rotamer\n");
4907   //                out.append(boxes.get(0) + "\n");
4908   //                int stopper = 0;
4909   //                for (int i = 0; i < eliminatedResidue.size(); i++) {
4910   //                    while (stopper == 0) {
4911   //                        for (int x = 0; x < (singlesIterations.size() - 1); x++) {
4912   //                            if (singlesIterations.get(x) == 0) {
4913   //                                String boxHeader = (boxes.get(x + 1));
4914   //                                out.append(boxHeader + System.lineSeparator());
4915   //                            }
4916   //                            stopper++;
4917   //                        }
4918   //                    }
4919   //                    String first = eliminatedResidue.get(i).toString();
4920   //                    String second = eliminatedRotamer.get(i).toString();
4921   //                    String resrot = first + "," + second;
4922   //                    out.append(resrot + System.lineSeparator());
4923   //                    for (int x = 0; x < (singlesIterations.size() - 1); x++) {
4924   //                        if ((i + 1) == singlesIterations.get(x)) {
4925   //                            String boxHeader = (boxes.get(x + 1));
4926   //                            out.append(boxHeader + System.lineSeparator());
4927   //                        }
4928   //                    }
4929   //                }
4930   //                out.append("Eliminated Pairs\n");
4931   //                out.append("Residue1,Rotamer1,Residue2,Rotamer2" + System.lineSeparator());
4932   //                out.append(boxes.get(0) + "\n");
4933   //                int stopper2 = 0;
4934   //                for (int i = 0; i < eliminatedResidue1.size(); i++) {
4935   //                    while (stopper2 == 0) {
4936   //                        for (int x = 0; x < (pairsIterations.size() - 1); x++) {
4937   //                            if (pairsIterations.get(x) == 0) {
4938   //                                String boxHeader = (boxes.get(x + 1));
4939   //                                out.append(boxHeader + System.lineSeparator());
4940   //                            }
4941   //                            stopper2++;
4942   //                        }
4943   //                    }
4944   //                    String first = eliminatedResidue1.get(i).toString();
4945   //                    String second = eliminatedRotamer1.get(i).toString();
4946   //                    String third = eliminatedResidue2.get(i).toString();
4947   //                    String fourth = eliminatedRotamer2.get(i).toString();
4948   //                    String resrot = first + "," + second + "," + third + "," + fourth;
4949   //                    out.append(resrot + System.lineSeparator());
4950   //                    for (int x = 0; x < (pairsIterations.size() - 1); x++) {
4951   //                        if ((i + 1) == pairsIterations.get(x)) {
4952   //                            String boxHeader = (boxes.get(x + 1));
4953   //                            out.append(boxHeader + System.lineSeparator());
4954   //                        }
4955   //                    }
4956   //                }
4957   //                /**
4958   //                 * Writes eliminated singles and pairs to eliminated.csv for
4959   //                 * global optimization.
4960   //                 */
4961   //            } else if ((line3.contains("-a, 1")) || (line3.contains("-a, 2")) ||
4962   // (line3.contains("-a, 3")) || (line3.contains("-a, 4"))) {
4963   //                out.append("Eliminated Singles\n");
4964   //                out.append("Residue,Rotamer\n");
4965   //                for (int i = 0; i < eliminatedResidue.size(); i++) {
4966   //                    String first = eliminatedResidue.get(i).toString();
4967   //                    String second = eliminatedRotamer.get(i).toString();
4968   //                    String resrot = first + "," + second;
4969   //                    out.append(resrot + System.lineSeparator());
4970   //                }
4971   //                out.append("Eliminated Pairs\n");
4972   //                out.append("Residue1,Rotamer1,Residue2,Rotamer2" + System.lineSeparator());
4973   //                for (int i = 0; i < eliminatedResidue1.size(); i++) {
4974   //                    String first = eliminatedResidue1.get(i).toString();
4975   //                    String second = eliminatedRotamer1.get(i).toString();
4976   //                    String third = eliminatedResidue2.get(i).toString();
4977   //                    String fourth = eliminatedRotamer2.get(i).toString();
4978   //                    String resrot = first + "," + second + "," + third + "," + fourth;
4979   //                    out.append(resrot + System.lineSeparator());
4980   //                }
4981   //            }
4982   //        }
4983   //        out.close();
4984   //        br.close();
4985   //    }
4986 }