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