Examples with Residue ffx.potential.bonded.Residue used on opensource projects

Example 71 with Residue

use of ffx.potential.bonded.Residue in project ffx by mjschnie.

the class RotamerOptimization method minMax2BodySum.

/**
 * Find the min/max of the 2-body energy.
 *
 * @param residues The residue array.
 * @param minMax The bound on the 3-body energy (minMax[0] = min, minMax[1]
 * = max.
 * @param i Residue i
 * @param ri Rotamer ri of Residue i
 * @param j Residue j
 * @param rj Rotamer rj of Residue j
 * @return true if this term is valid.
 */
private boolean minMax2BodySum(Residue[] residues, double[] minMax, int i, int ri, int j, int rj) {
    int nres = residues.length;
    double minSum = 0.0;
    double maxSum = 0.0;
    for (int k = 0; k < nres; k++) {
        if (k == i || k == j) {
            continue;
        }
        Residue residuek = residues[k];
        Rotamer[] romatersk = residuek.getRotamers(library);
        int lenrk = romatersk.length;
        double currentMin = Double.MAX_VALUE;
        double currentMax = Double.MIN_VALUE;
        for (int rk = 0; rk < lenrk; rk++) {
            if (check(k, rk)) {
                // k,rk is part of no valid phase space, so ignore it.
                continue;
            }
            if (check(i, ri, k, rk) || check(j, rj, k, rk)) {
                // Not implemented: check(i, ri, j, rj, k, rk).
                // k,rk conflicts with i,ri or j,rj, so the max is now Double.NaN. No effect on minimum.
                currentMax = Double.NaN;
            } else {
                double current = get3Body(i, ri, j, rj, k, rk);
                if (Double.isFinite(current) && current < currentMin) {
                    currentMin = current;
                }
                // Else, no new minimum found.
                if (Double.isFinite(current) && Double.isFinite(currentMax)) {
                    if (current > currentMax) {
                        currentMax = current;
                    }
                // Else, we have failed to find a new finite maximum.
                } else {
                    // The maximum is NaN.
                    currentMax = Double.NaN;
                }
            }
        }
        if (currentMin == Double.MAX_VALUE || !Double.isFinite(minSum)) {
            // We have failed to find a viable configuration for i,ri,j,rj, as it conflicts with all rk for this k.
            minMax[0] = Double.NaN;
            minMax[1] = Double.NaN;
            return false;
        } else {
            // Add finite current min to minSum.
            minSum += currentMin;
        }
        if (Double.isFinite(maxSum) && Double.isFinite(currentMax)) {
            maxSum += currentMax;
        } else {
            maxSum = Double.NaN;
        }
    }
    minMax[0] = minSum;
    minMax[1] = maxSum;
    return true;
}

Example 72 with Residue

use of ffx.potential.bonded.Residue in project ffx by mjschnie.

the class RotamerOptimization method dryRunForEnsemble.

/**
 * Finds all permutations within buffer energy of GMEC.
 *
 * @param residues
 * @param i Current depth in residue/rotamer tree.
 * @param currentRotamers Current set of rotamers at this node.
 * @param gmecEnergy Minimum energy for these residues.
 * @param permutationEnergies Energy of all permutations.
 * @param permutations Contains accepted permutations.
 * @return 0.
 */
private double dryRunForEnsemble(Residue[] residues, int i, int[] currentRotamers, double gmecEnergy, double[] permutationEnergies, int[][] permutations) {
    // This is the initialization condition.
    if (i == 0) {
        evaluatedPermutations = 0;
    }
    int nResidues = residues.length;
    Residue residuei = residues[i];
    Rotamer[] rotamersi = residuei.getRotamers(library);
    int lenri = rotamersi.length;
    if (i < nResidues - 1) {
        for (int ri = 0; ri < lenri; ri++) {
            if (check(i, ri)) {
                continue;
            }
            boolean deadEnd = false;
            for (int j = 0; j < i; j++) {
                int rj = currentRotamers[j];
                deadEnd = check(j, rj, i, ri);
                if (deadEnd) {
                    break;
                }
            }
            if (deadEnd) {
                continue;
            }
            currentRotamers[i] = ri;
            dryRunForEnsemble(residues, i + 1, currentRotamers, gmecEnergy, permutationEnergies, permutations);
        }
    } else {
        /**
         * At the end of the recursion, check each rotamer of the final
         * residue.
         */
        for (int ri = 0; ri < lenri; ri++) {
            if (check(i, ri)) {
                continue;
            }
            currentRotamers[i] = ri;
            boolean deadEnd = false;
            for (int j = 0; j < i; j++) {
                int rj = currentRotamers[j];
                deadEnd = check(j, rj, i, ri);
                if (deadEnd) {
                    break;
                }
            }
            if (deadEnd) {
                continue;
            }
            if (permutationEnergies[evaluatedPermutations] - gmecEnergy < ensembleEnergy) {
                permutations[evaluatedPermutations] = new int[nResidues];
                System.arraycopy(currentRotamers, 0, permutations[evaluatedPermutations], 0, nResidues);
            }
            evaluatedPermutations++;
        }
    }
    return 0.0;
}

Example 73 with Residue

use of ffx.potential.bonded.Residue in project ffx by mjschnie.

the class RotamerOptimization method decomposePrint.

/**
 * Prints based on a molecule-based decomposition of original energies. Would need to be re-implemented.
 * @param residues
 * @param totalEnergy
 * @param backbone
 * @param residueEnergy
 */
private void decomposePrint(Residue[] residues, double totalEnergy, double backbone, double[][] residueEnergy) {
    if (residues == null || residueEnergy == null) {
        return;
    }
    int[] molecule = molecularAssembly.getMoleculeNumbers();
    // Internal Molecule ID to Reduced Set of Molecule IDs
    HashMap<Integer, Integer> nMolecules = new HashMap<>();
    // Internal Molecule ID to Molecule Names
    HashMap<Integer, String> molNames = new HashMap<>();
    // Residue ID to Reduced Set of Molecule IDs
    HashMap<Residue, Integer> moleculeMap = new HashMap<>();
    /**
     * Sum self, 2-Body and 3-Body energies.
     */
    int nRes = residues.length;
    double sumSelf = 0;
    double sumPair = 0;
    double sumTri = 0;
    int molPointer = 0;
    int reducedIndex = 0;
    for (int i = 0; i < nRes; i++) {
        sumSelf += residueEnergy[0][i];
        sumPair += residueEnergy[1][i];
        sumTri += residueEnergy[2][i];
        // Count the number of molecules.
        Residue r = residues[i];
        Atom a0 = (Atom) r.getAtomNode(0);
        int atomIndex = a0.getIndex();
        Integer moleculeIndex = molecule[atomIndex];
        if (!nMolecules.containsKey(moleculeIndex)) {
            nMolecules.put(moleculeIndex, molPointer);
            molNames.put(moleculeIndex, r.getSegID());
            reducedIndex = molPointer;
            molPointer++;
        } else {
            reducedIndex = nMolecules.get(moleculeIndex);
        }
        moleculeMap.put(r, reducedIndex);
    }
    int nMol = nMolecules.size();
    if (nMol > 1) {
        logger.info(format("\n"));
        logger.info(format(" Molecule-Based Many-Body Energy Summation\n "));
        double[][] moleculeEnergy = new double[3][nMol];
        for (int i = 0; i < nRes; i++) {
            int iMolecule = moleculeMap.get(residues[i]);
            moleculeEnergy[0][iMolecule] += residueEnergy[0][i];
            moleculeEnergy[1][iMolecule] += residueEnergy[1][i];
            moleculeEnergy[2][iMolecule] += residueEnergy[2][i];
        }
        logger.info(format(" %9s %9s %9s %9s %9s", "Molecule", "Self", "Pair", "3-Body", "Total"));
        for (int i = 0; i < nMol; i++) {
            String molName = molNames.get(i);
            double total = moleculeEnergy[0][i] + moleculeEnergy[1][i] + moleculeEnergy[2][i];
            logger.info(format(" %9s %9.3f %9.3f %9.3f %9.3f", molName, moleculeEnergy[0][i], moleculeEnergy[1][i], moleculeEnergy[2][i], total));
        }
        logger.info(format(" %9s %9.3f %9.3f %9.3f %9.3f", "Sum", sumSelf, sumPair, sumTri, sumSelf + sumPair + sumTri));
    }
    logger.info(format("\n"));
    logger.info(format(" Residue-Based Many-Body Energy Summation\n "));
    logger.info(format(" %9s %9s %9s %9s %9s", "Residue", "Self", "Pair", "3-Body", "Total"));
    for (int i = 0; i < nRes; i++) {
        double total = residueEnergy[0][i] + residueEnergy[1][i] + residueEnergy[2][i];
        logger.info(format(" %9s %9.3f %9.3f %9.3f %9.3f", residues[i].toFormattedString(false, true), residueEnergy[0][i], residueEnergy[1][i], residueEnergy[2][i], total));
    }
    logger.info(format(" %9s %9.3f %9.3f %9.3f %9.3f", "Sum", sumSelf, sumPair, sumTri, sumSelf + sumPair + sumTri));
    logger.info(format(" Backbone:        %9.3f", backbone));
    double target = sumSelf + sumPair + sumTri + backbone;
    logger.info(format(" Expansion Total: %9.3f", target));
    logger.info(format(" True Total:      %9.3f", totalEnergy));
    logger.info(format(" Neglected:       %9.3f\n", totalEnergy - target));
}

Example 74 with Residue

use of ffx.potential.bonded.Residue in project ffx by mjschnie.

the class RotamerOptimization method deeRotamerPairElimination.

/**
 * Rotamer pair elimination driver for many-body Dead End Elimination.
 * Generally less effective than Goldstein.
 *
 * @param residues Residues under consideration.
 * @return If at least one pair eliminated.
 */
private boolean deeRotamerPairElimination(Residue[] residues) {
    int nres = residues.length;
    boolean eliminated = false;
    for (int i = 0; i < (nres - 1); i++) {
        Residue residuei = residues[i];
        Rotamer[] rotamersi = residuei.getRotamers(library);
        int lenri = rotamersi.length;
        // Minimum and maximum summation found for ri-rj pairs.
        double[][] minPairEnergies = new double[lenri][];
        double[][] maxPairEnergies = new double[lenri][];
        for (int j = i + 1; j < nres; j++) {
            Residue residuej = residues[j];
            Rotamer[] rotamersj = residuej.getRotamers(library);
            int lenrj = rotamersj.length;
            for (int ri = 0; ri < lenri; ri++) {
                if (check(i, ri)) {
                    continue;
                }
                minPairEnergies[ri] = new double[lenrj];
                maxPairEnergies[ri] = new double[lenrj];
                for (int rj = 0; rj < lenrj; rj++) {
                    if (check(j, rj) || check(i, ri, j, rj)) {
                        continue;
                    }
                    minPairEnergies[ri][rj] = getSelf(i, ri) + getSelf(j, rj) + get2Body(i, ri, j, rj);
                    maxPairEnergies[ri][rj] = minPairEnergies[ri][rj];
                    // Min and max external summations for ri-rj.
                    double[] minMax = new double[2];
                    // Add contributions from third residues k, and possibly fourth residues l.
                    if (minMaxE2(residues, minMax, i, ri, j, rj)) {
                        if (Double.isFinite(minPairEnergies[ri][rj]) && Double.isFinite(minMax[0])) {
                            minPairEnergies[ri][rj] += minMax[0];
                        } else {
                            logger.severe(String.format(" An ri-rj pair %s-%d %s-%d with NaN minimum was caught incorrectly!", residuei.toFormattedString(false, true), ri, residuej.toFormattedString(false, true), rj));
                        }
                        if (Double.isFinite(maxPairEnergies[ri][rj]) && Double.isFinite(minMax[1])) {
                            maxPairEnergies[ri][rj] += minMax[1];
                        } else {
                            // ri-rj can clash, and isn't very useful to eliminate by.
                            maxPairEnergies[ri][rj] = Double.NaN;
                        }
                    } else {
                        // A NaN minimum energy for some pair indicates it's definitely not part of the GMEC.
                        minPairEnergies[ri][rj] = Double.NaN;
                        logger.info(String.format(" Eliminating pair %s-%d %s-%d that always clashes.", residuei.toFormattedString(false, true), ri, residuej.toFormattedString(false, true), rj));
                        eliminateRotamerPair(residues, i, ri, j, rj, print);
                        eliminated = true;
                    }
                }
            }
            double pairEliminationEnergy = Double.MAX_VALUE;
            for (int ri = 0; ri < lenri; ri++) {
                if (check(i, ri)) {
                    continue;
                }
                for (int rj = 0; rj < lenrj; rj++) {
                    if (check(j, rj) || check(i, ri, j, rj)) {
                        continue;
                    }
                    if (Double.isFinite(maxPairEnergies[ri][rj]) && maxPairEnergies[ri][rj] < pairEliminationEnergy) {
                        pairEliminationEnergy = maxPairEnergies[ri][rj];
                    }
                }
            }
            if (pairEliminationEnergy == Double.MAX_VALUE) {
                logIfMaster(String.format(" All rotamer pairs for residues %s and %s have possible conflicts; cannot perform any eliminations!", residuei.toFormattedString(false, true), residuej), Level.FINE);
            } else {
                double comparisonEnergy = pairEliminationEnergy + ensembleBuffer;
                for (int ri = 0; ri < lenri; ri++) {
                    if (check(i, ri)) {
                        continue;
                    }
                    for (int rj = 0; rj < lenrj; rj++) {
                        if (check(j, rj) || check(i, ri, j, rj)) {
                            continue;
                        }
                        if (minPairEnergies[ri][rj] > comparisonEnergy) {
                            if (eliminateRotamerPair(residues, i, ri, j, rj, print)) {
                                eliminated = true;
                                logIfMaster(format(" Eliminating rotamer pair: %s %d, %s %d (%s > %s + %6.6f)", residuei.toFormattedString(false, true), ri, residuej.toFormattedString(false, true), rj, formatEnergy(minPairEnergies[ri][rj]), formatEnergy(pairEliminationEnergy), ensembleBuffer), Level.INFO);
                            } else {
                                // See above check(i, ri, j, rj) for why this should not be taken!
                                logIfMaster(format(" Already eliminated rotamer pair! %s %d, %s %d (%s > %1s + %6.6f)", residuei.toFormattedString(false, true), ri, residuej.toFormattedString(false, true), rj, formatEnergy(minPairEnergies[ri][rj]), formatEnergy(pairEliminationEnergy), ensembleBuffer), Level.WARNING);
                            }
                        }
                    }
                }
            }
            if (pairsToSingleElimination(residues, i, j)) {
                eliminated = true;
            }
        }
    }
    return eliminated;
}

Example 75 with Residue

use of ffx.potential.bonded.Residue in project ffx by mjschnie.

the class RotamerOptimization method globalOptimization.

/**
 * The main driver for optimizing a block of residues using DEE.
 *
 * @param residueList Residues to optimize.
 * @return Final energy.
 */
private double globalOptimization(List<Residue> residueList) {
    int currentEnsemble = Integer.MAX_VALUE;
    Residue[] residues = residueList.toArray(new Residue[residueList.size()]);
    int nResidues = residues.length;
    int[] currentRotamers = new int[nResidues];
    int iterations = 0;
    boolean finalTry = false;
    int bestEnsembleTargetDiffThusFar = Integer.MAX_VALUE;
    double bestBufferThusFar = ensembleBuffer;
    double startingBuffer = ensembleBuffer;
    optimum = new int[nResidues];
    if (ensembleEnergy > 0.0) {
        ensembleBuffer = ensembleEnergy;
        applyEliminationCriteria(residues, true, true);
        if (x == null) {
            Atom[] atoms = molecularAssembly.getAtomArray();
            int nAtoms = atoms.length;
            x = new double[nAtoms * 3];
        }
        /**
         * Compute the number of permutations without eliminating dead-ends
         * and compute the number of permutations using singleton
         * elimination.
         */
        double permutations = 1;
        double singletonPermutations = 1;
        for (int i = 0; i < nResidues; i++) {
            Residue residue = residues[i];
            Rotamer[] rotamers = residue.getRotamers(library);
            int nr = rotamers.length;
            if (nr > 1) {
                int nrot = 0;
                for (int ri = 0; ri < nr; ri++) {
                    if (!eliminatedSingles[i][ri]) {
                        nrot++;
                    }
                }
                permutations *= rotamers.length;
                if (nrot > 1) {
                    singletonPermutations *= nrot;
                }
            }
        }
        dryRun(residues, 0, currentRotamers);
        double pairTotalElimination = singletonPermutations - (double) evaluatedPermutations;
        double afterPairElim = singletonPermutations - pairTotalElimination;
        if (evaluatedPermutations == 0) {
            logger.severe(" No valid path through rotamer space found; try recomputing without pruning or using ensemble.");
        }
        if (master && printFiles && ensembleFile == null) {
            File file = molecularAssembly.getFile();
            String filename = FilenameUtils.removeExtension(file.getAbsolutePath());
            ensembleFile = new File(filename + ".ens");
            if (ensembleFile.exists()) {
                for (int i = 2; i < 1000; i++) {
                    ensembleFile = new File(filename + ".ens_" + i);
                    if (!ensembleFile.exists()) {
                        break;
                    }
                }
                if (ensembleFile.exists()) {
                    logger.warning(format(" Versioning failed: appending to end of file %s", ensembleFile.getName()));
                }
            }
            ensembleFilter = new PDBFilter(new File(ensembleFile.getName()), molecularAssembly, null, null);
            logger.info(format(" Ensemble file: %s", ensembleFile.getName()));
        }
        logIfMaster(format("%30s %5s %30s %5s %30s %5s", "Condition", "|", "Number of Permutations Left", "|", "Removed", "|"));
        logIfMaster(format("%s", " -------------------------------------------------------------------------------------------------------------"));
        logIfMaster(format("%30s %5s %30s %5s %30s %5s", "No Eliminations", "|", permutations, "|", "", "|"));
        logIfMaster(format("%30s %5s %30s %5s %30s %5s", "Single Eliminations", "|", singletonPermutations, "|", permutations - singletonPermutations, "|"));
        logIfMaster(format("%30s %5s %30s %5s %30s %5s", "Pair Eliminations", "|", afterPairElim, "|", pairTotalElimination, "|"));
        logIfMaster(format("%30s %5s %30s %5s %30s %5s", "Single and Pair Eliminations", "|", (double) evaluatedPermutations, "|", pairTotalElimination + (permutations - singletonPermutations), "|"));
        logIfMaster(format("%s", " -------------------------------------------------------------------------------------------------------------\n"));
        logIfMaster(format(" Energy of permutations:"));
        logIfMaster(format("%s", " ----------------------------------------------------------------------------------"));
        logIfMaster(format(" %12s %5s %25s %5s %25s %5s", "Permutation", "|", "Energy", "|", "Lowest Possible Energy", "|"));
        logIfMaster(format("%s", " ----------------------------------------------------------------------------------"));
        double e;
        if (useMonteCarlo()) {
            firstValidPerm(residues, 0, currentRotamers);
            System.arraycopy(currentRotamers, 0, optimum, 0, nResidues);
            rotamerOptimizationMC(residues, optimum, currentRotamers, nMCsteps, false, mcUseAll);
            logIfMaster(" Ensembles not currently compatible with Monte Carlo search");
        /**
         * Not currently compatible with ensembles.
         */
        } else {
            double[] permutationEnergies = new double[evaluatedPermutations];
            ensembleStates = new ArrayList<>();
            e = rotamerOptimizationDEE(molecularAssembly, residues, 0, currentRotamers, Double.MAX_VALUE, optimum, permutationEnergies);
            int[][] acceptedPermutations = new int[evaluatedPermutations][];
            for (int i = 0; i < acceptedPermutations.length; i++) {
                acceptedPermutations[i] = null;
            }
            logIfMaster(format("\n Checking permutations for distance < %5.3f kcal/mol from GMEC energy %10.8f kcal/mol", ensembleEnergy, e));
            dryRunForEnsemble(residues, 0, currentRotamers, e, permutationEnergies, acceptedPermutations);
            int numAcceptedPermutations = 0;
            for (int i = 0; i < acceptedPermutations.length; i++) {
                if (acceptedPermutations[i] != null) {
                    ++numAcceptedPermutations;
                    logIfMaster(format(" Accepting permutation %d at %8.6f < %8.6f", i, permutationEnergies[i] - e, ensembleEnergy));
                    for (int j = 0; j < nResidues; j++) {
                        Residue residuej = residues[j];
                        Rotamer[] rotamersj = residuej.getRotamers(library);
                        RotamerLibrary.applyRotamer(residuej, rotamersj[acceptedPermutations[i][j]]);
                    }
                    ResidueState[] states = ResidueState.storeAllCoordinates(residues);
                    ensembleStates.add(new ObjectPair<>(states, permutationEnergies[i]));
                    if (printFiles && master) {
                        try {
                            FileWriter fw = new FileWriter(ensembleFile, true);
                            BufferedWriter bw = new BufferedWriter(fw);
                            bw.write(format("MODEL        %d", numAcceptedPermutations));
                            for (int j = 0; j < 75; j++) {
                                bw.write(" ");
                            }
                            bw.newLine();
                            bw.flush();
                            ensembleFilter.writeFile(ensembleFile, true);
                            bw.write(format("ENDMDL"));
                            for (int j = 0; j < 64; j++) {
                                bw.write(" ");
                            }
                            bw.newLine();
                            bw.close();
                        } catch (IOException ex) {
                            logger.warning(format(" Exception writing to file: %s", ensembleFile.getName()));
                        }
                    }
                }
            }
            logIfMaster(format(" Number of permutations within %5.3f kcal/mol of GMEC energy: %6.4e", ensembleEnergy, (double) numAcceptedPermutations));
            ensembleStates.sort(null);
        }
        logIfMaster(format(" Final rotamers:"));
        logIfMaster(format("%s", " --------------------------------------------------------------------------------------------"));
        logIfMaster(format("%14s %3s %10s %3s %9s %3s %9s %3s %9s %3s", "Residue", "|", "Chi 1", "|", "Chi 2", "|", "Chi 3", "|", "Chi 4", "|"));
        logIfMaster(format("%s", " --------------------------------------------------------------------------------------------"));
        for (int i = 0; i < nResidues; i++) {
            Residue residue = residues[i];
            Rotamer[] rotamers = residue.getRotamers(library);
            int ri = optimum[i];
            Rotamer rotamer = rotamers[ri];
            logIfMaster(format(" %c (%7s,2d) | %s", residue.getChainID(), residue, ri, rotamer.toAngleString()));
            RotamerLibrary.applyRotamer(residue, rotamer);
        }
        logIfMaster(format("%s", " --------------------------------------------------------------------------------------------\n"));
        double sumSelfEnergy = 0;
        double sumPairEnergy = 0;
        double sumTrimerEnergy = 0;
        for (int i = 0; i < nResidues; i++) {
            int ri = optimum[i];
            sumSelfEnergy += getSelf(i, ri);
            logIfMaster(format(" Final self Energy (%8s,%2d): %12.4f", residues[i].toFormattedString(false, true), ri, getSelf(i, ri)));
        }
        for (int i = 0; i < nResidues - 1; i++) {
            int ri = optimum[i];
            for (int j = i + 1; j < nResidues; j++) {
                int rj = optimum[j];
                sumPairEnergy += get2Body(i, ri, j, rj);
                if (get2Body(i, ri, j, rj) > 10.0) {
                    logIfMaster(format(" Large Final Pair Energy (%8s,%2d) (%8s,%2d): %12.4f", residues[i].toFormattedString(false, true), ri, residues[j].toFormattedString(false, true), rj, get2Body(i, ri, j, rj)));
                }
            }
        }
        try {
            e = currentEnergy(residueList);
        } catch (ArithmeticException ex) {
            e = Double.NaN;
            logger.severe(String.format(" Exception %s in calculating current energy at the end of triples", ex.toString()));
        }
        logIfMaster(format(" %12s %5s %25s %5s %25s %5s", "Type", "|", "Energy", "|", "Lowest Possible Energy", "|"));
        logIfMaster(format("%s", " ----------------------------------------------------------------------------------"));
        logIfMaster(format(" %12s %5s %25f %5s %25s %5s", "Self:", "|", sumSelfEnergy, "|", "", "|"));
        logIfMaster(format(" %12s %5s %25f %5s %25s %5s", "Pair:", "|", sumPairEnergy, "|", "", "|"));
        double approximateEnergy = backboneEnergy + sumSelfEnergy + sumPairEnergy;
        if (threeBodyTerm) {
            for (int i = 0; i < nResidues - 2; i++) {
                int ri = optimum[i];
                for (int j = i + 1; j < nResidues - 1; j++) {
                    int rj = optimum[j];
                    for (int k = j + 1; k < nResidues; k++) {
                        int rk = optimum[k];
                        try {
                            sumTrimerEnergy += get3Body(i, ri, j, rj, k, rk);
                        } catch (Exception ex) {
                            logger.warning(ex.toString());
                        }
                    }
                }
            }
            approximateEnergy += sumTrimerEnergy;
            double higherOrderEnergy = e - sumSelfEnergy - sumPairEnergy - sumTrimerEnergy - backboneEnergy;
            logIfMaster(format(" %12s %5s %25f %5s %25s %5s", "Trimer:", "|", sumTrimerEnergy, "|", "", "|"));
            logIfMaster(format(" %12s %5s %25f %5s %25s %5s", "Neglected:", "|", higherOrderEnergy, "|", "", "|"));
        } else {
            double higherOrderEnergy = e - sumSelfEnergy - sumPairEnergy - backboneEnergy;
            logIfMaster(format(" %12s %5s %25f %5s %25s %5s", "Neglected:", "|", higherOrderEnergy, "|", "", "|"));
        }
        logIfMaster(format(" %12s %5s %25f %5s %25s %5s", "Approximate:", "|", approximateEnergy, "|", "", "|"));
        logIfMaster(format("%s", " ----------------------------------------------------------------------------------\n"));
        return e;
    }
    /**
     * Permutations used only to set maximum bound on ensembleNumber, thus
     * it is safe here to put that value in a 32-bit int.
     */
    int nPerms = 1;
    for (int i = 0; i < nResidues; i++) {
        Residue residue = residues[i];
        Rotamer[] rotamers = residue.getRotamers(library);
        int nr = rotamers.length;
        if (nr > 1) {
            nPerms *= rotamers.length;
        }
        if (nPerms > ensembleNumber) {
            break;
        }
    }
    if (nPerms < ensembleNumber) {
        logger.info(format(" Requested an ensemble of %d, but only %d permutations exist; returning full ensemble", ensembleNumber, nPerms));
        ensembleNumber = nPerms;
    }
    while (currentEnsemble != ensembleNumber) {
        if (monteCarlo) {
            logIfMaster(" Ensemble search not currently compatible with Monte Carlo");
            ensembleNumber = 1;
        }
        if (iterations == 0) {
            applyEliminationCriteria(residues, true, true);
        } else {
            applyEliminationCriteria(residues, false, false);
        }
        if (x == null) {
            Atom[] atoms = molecularAssembly.getAtomArray();
            int nAtoms = atoms.length;
            x = new double[nAtoms * 3];
        }
        /**
         * Compute the number of permutations without eliminating dead-ends
         * and compute the number of permutations using singleton
         * elimination.
         */
        double permutations = 1;
        double singletonPermutations = 1;
        for (int i = 0; i < nResidues; i++) {
            Residue residue = residues[i];
            Rotamer[] rotamers = residue.getRotamers(library);
            int nr = rotamers.length;
            if (nr > 1) {
                int nrot = 0;
                for (int ri = 0; ri < nr; ri++) {
                    if (!eliminatedSingles[i][ri]) {
                        nrot++;
                    }
                }
                permutations *= rotamers.length;
                if (nrot > 1) {
                    singletonPermutations *= nrot;
                }
            }
        }
        logIfMaster(format(" Collecting Permutations:"));
        logIfMaster(format("%s", " -------------------------------------------------------------------------------------------------------------"));
        dryRun(residues, 0, currentRotamers);
        double pairTotalElimination = singletonPermutations - (double) evaluatedPermutations;
        double afterPairElim = singletonPermutations - pairTotalElimination;
        currentEnsemble = (int) evaluatedPermutations;
        if (ensembleNumber == 1 && currentEnsemble == 0) {
            logger.severe(" No valid path through rotamer space found; try recomputing without pruning or using ensemble.");
        }
        if (ensembleNumber > 1) {
            if (master && printFiles && ensembleFile == null) {
                File file = molecularAssembly.getFile();
                String filename = FilenameUtils.removeExtension(file.getAbsolutePath());
                ensembleFile = new File(filename + ".ens");
                if (ensembleFile.exists()) {
                    for (int i = 2; i < 1000; i++) {
                        ensembleFile = new File(filename + ".ens_" + i);
                        if (!ensembleFile.exists()) {
                            break;
                        }
                    }
                    if (ensembleFile.exists()) {
                        logger.warning(format(" Versioning failed: appending to end of file %s", ensembleFile.getName()));
                    }
                }
                ensembleFilter = new PDBFilter(new File(ensembleFile.getName()), molecularAssembly, null, null);
                logger.info(format(" Ensemble file: %s", ensembleFile.getName()));
            }
            logIfMaster(format(" Ensemble Search Stats: (buffer: %5.3f, current: %d, target: %d)", ensembleBuffer, currentEnsemble, ensembleNumber));
        }
        if (ensembleNumber == 1 || finalTry) {
            logIfMaster(format("%30s %5s %30s %5s %30s %5s", "Condition", "|", "Number of Permutations Left", "|", "Number of Permutations Removed", "|"));
            logIfMaster(format("%s", " -------------------------------------------------------------------------------------------------------------"));
            logIfMaster(format("%30s %5s %30s %5s %30s %5s", "No Eliminations", "|", permutations, "|", "", "|"));
            logIfMaster(format("%30s %5s %30s %5s %30s %5s", "Single Eliminations", "|", singletonPermutations, "|", permutations - singletonPermutations, "|"));
            logIfMaster(format("%30s %5s %30s %5s %30s %5s", "Pair Eliminations", "|", afterPairElim, "|", pairTotalElimination, "|"));
            logIfMaster(format("%30s %5s %30s %5s %30s %5s", "Single and Pair Eliminations", "|", (double) evaluatedPermutations, "|", pairTotalElimination + (permutations - singletonPermutations), "|"));
            logIfMaster(format("%s", " -------------------------------------------------------------------------------------------------------------\n"));
            logIfMaster(format(" Energy of permutations:"));
            logIfMaster(format("%s", " ----------------------------------------------------------------------------------"));
            logIfMaster(format(" %12s %5s %25s %5s %25s %5s", "Permutation", "|", "Energy", "|", "Lowest Possible Energy", "|"));
            logIfMaster(format("%s", " ----------------------------------------------------------------------------------"));
            break;
        }
        if (Math.abs(currentEnsemble - ensembleNumber) < bestEnsembleTargetDiffThusFar) {
            bestEnsembleTargetDiffThusFar = Math.abs(currentEnsemble - ensembleNumber);
            bestBufferThusFar = ensembleBuffer;
        }
        if (currentEnsemble > ensembleNumber) {
            ensembleBuffer -= ensembleBufferStep;
            ensembleBufferStep -= (ensembleBufferStep * 0.01);
            iterations++;
        } else if (currentEnsemble < ensembleNumber) {
            ensembleBuffer += ensembleBufferStep;
            ensembleBufferStep -= (ensembleBufferStep * 0.01);
            iterations++;
        }
        if (iterations > 100) {
            if (currentEnsemble == 0) {
                // TODO: Decide whether we like these next four lines.  Has the potential to produce a crazy amount of permutations.
                logIfMaster(" Ensemble still empty; increasing buffer energy.");
                startingBuffer = 3 * startingBuffer;
                setEnsemble(10, startingBuffer);
                iterations = 0;
            } else {
                ensembleBuffer = bestBufferThusFar;
                finalTry = true;
            }
        }
    }
    if (currentEnsemble == 0) {
        logger.warning(" No valid rotamer permutations found; results will be unreliable.  Try increasing the starting ensemble buffer.");
    }
    double[] permutationEnergyStub = null;
    if (useMonteCarlo()) {
        firstValidPerm(residues, 0, currentRotamers);
        rotamerOptimizationMC(residues, optimum, currentRotamers, nMCsteps, false, mcUseAll);
    } else {
        rotamerOptimizationDEE(molecularAssembly, residues, 0, currentRotamers, Double.MAX_VALUE, optimum, permutationEnergyStub);
    }
    double[] residueEnergy = new double[nResidues];
    double sumSelfEnergy = 0;
    double sumLowSelfEnergy = 0;
    logIfMaster(format("%s", " ----------------------------------------------------------------------------------\n"));
    logIfMaster(format(" Energy contributions:"));
    logIfMaster(format("%s", " -------------------------------------------------------------------------------------"));
    logIfMaster(format(" %15s %5s %25s %5s %25s %5s", "Type", "|", "Energy", "|", "Lowest Possible Energy", "|"));
    logIfMaster(format("%s", " -------------------------------------------------------------------------------------"));
    for (int i = 0; i < nResidues; i++) {
        int ri = optimum[i];
        Residue residue = residues[i];
        Rotamer[] rotamers = residue.getRotamers(library);
        turnOnAtoms(residue);
        RotamerLibrary.applyRotamer(residue, rotamers[ri]);
        double self = getSelf(i, ri);
        residueEnergy[i] = self;
        sumSelfEnergy += self;
        double lowest = lowestSelfEnergy(residues, i);
        sumLowSelfEnergy += lowest;
        if (self - lowest > 10.0) {
            logIfMaster(format(" %15s %5s %25f %5s %25f %5s", "Self (" + residues[i] + "," + ri + "):", "|", self, "|", lowest, "|"));
        }
    }
    double sumPairEnergy = 0.0;
    double sumLowPairEnergy = 0.0;
    double[] resPairEnergy = new double[nResidues];
    double[] lowPairEnergy = new double[nResidues];
    for (int i = 0; i < nResidues - 1; i++) {
        StringBuilder sb = new StringBuilder();
        int ri = optimum[i];
        double sumPairEnergyI = 0;
        double sumLowPairEnergyI = 0;
        for (int j = i + 1; j < nResidues; j++) {
            int rj = optimum[j];
            double pair = get2Body(i, ri, j, rj);
            residueEnergy[i] += 0.5 * pair;
            residueEnergy[j] += 0.5 * pair;
            sumPairEnergy += pair;
            sumPairEnergyI += pair;
            double lowest = lowestPairEnergy(residues, i, ri, j);
            sumLowPairEnergy += lowest;
            sumLowPairEnergyI += lowest;
            resPairEnergy[i] = 0.5 * pair;
            resPairEnergy[j] = 0.5 * pair;
            lowPairEnergy[i] = 0.5 * lowest;
            lowPairEnergy[j] = 0.5 * lowest;
            if (resPairEnergy[i] - lowPairEnergy[i] > 10.0) {
                sb.append(format("  Pair Energy (%8s,%2d) (%8s,%2d): %12.4f (Lowest: %12.4f).\n", residues[i].toFormattedString(false, true), ri, residues[j].toFormattedString(false, true), rj, pair, lowest));
            }
        }
        if (sumPairEnergyI - sumLowPairEnergyI > 10.0) {
            logIfMaster(format(" %15s %5s %25f %5s %25f %5s", "Self (" + residues[i] + "," + ri + "):", "|", sumPairEnergyI, "|", sumLowPairEnergyI, "|"));
            sb.trimToSize();
            if (!sb.toString().isEmpty()) {
                logIfMaster(sb.toString());
            }
        }
    }
    double e = Double.NaN;
    try {
        e = currentEnergy(residueList);
    } catch (ArithmeticException ex) {
        logger.severe(String.format(" Exception %s in calculating current energy at the end of self and pairs", ex.toString()));
    }
    logIfMaster(format(" %15s %5s %25f %5s %25s %5s", "Backbone:", "|", backboneEnergy, "|", "", "|"));
    logIfMaster(format(" %15s %5s %25f %5s %25f %5s", "Self:", "|", sumSelfEnergy, "|", sumLowSelfEnergy, "|"));
    logIfMaster(format(" %15s %5s %25f %5s %25f %5s", "Pair:", "|", sumPairEnergy, "|", sumLowPairEnergy, "|"));
    double approximateEnergy = backboneEnergy + sumSelfEnergy + sumPairEnergy;
    double sumTrimerEnergy = 0;
    if (threeBodyTerm) {
        for (int i = 0; i < nResidues - 2; i++) {
            int ri = optimum[i];
            for (int j = i + 1; j < nResidues - 1; j++) {
                int rj = optimum[j];
                for (int k = j + 1; k < nResidues; k++) {
                    int rk = optimum[k];
                    try {
                        double triple = get3Body(i, ri, j, rj, k, rk);
                        double thirdTrip = triple / 3.0;
                        residueEnergy[i] += thirdTrip;
                        residueEnergy[j] += thirdTrip;
                        residueEnergy[k] += thirdTrip;
                        sumTrimerEnergy += triple;
                    } catch (Exception ex) {
                        logger.warning(ex.toString());
                    }
                }
            }
        }
        approximateEnergy += sumTrimerEnergy;
        double higherOrderEnergy = e - sumSelfEnergy - sumPairEnergy - sumTrimerEnergy - backboneEnergy;
        logIfMaster(format(" %15s %5s %25f %5s %25s %5s", "Trimer:", "|", sumTrimerEnergy, "|", "", "|"));
        logIfMaster(format(" %15s %5s %25f %5s %25s %5s", "Neglected:", "|", higherOrderEnergy, "|", "", "|"));
    } else {
        double higherOrderEnergy = e - sumSelfEnergy - sumPairEnergy - backboneEnergy;
        logIfMaster(format(" %15s %5s %25f %5s %25s %5s", "Neglected:", "|", higherOrderEnergy, "|", "", "|"));
    }
    logIfMaster(format(" %15s %5s %25f %5s %25s %5s", "Approximate:", "|", approximateEnergy, "|", "", "|"));
    logIfMaster(format("%s", " -------------------------------------------------------------------------------------\n"));
    logIfMaster(format(" Final rotamers:"));
    logIfMaster(format("%s", " --------------------------------------------------------------------------------------------"));
    logIfMaster(format("%17s %3s %10s %3s %9s %3s %9s %3s %9s %3s %10s %3s", "Residue", "|", "Chi 1", "|", "Chi 2", "|", "Chi 3", "|", "Chi 4", "|", "Energy", "|"));
    logIfMaster(format("%s", " --------------------------------------------------------------------------------------------"));
    for (int i = 0; i < nResidues; i++) {
        Residue residue = residues[i];
        Rotamer[] rotamers = residue.getRotamers(library);
        int ri = optimum[i];
        Rotamer rotamer = rotamers[ri];
        logIfMaster(format(" %3d %c (%7s,%2d) | %s %12.4f |", i + 1, residue.getChainID(), residue, ri, rotamer.toAngleString(), residueEnergy[i]));
        RotamerLibrary.applyRotamer(residue, rotamer);
    }
    logIfMaster(format("%s", " --------------------------------------------------------------------------------------------\n"));
    return e;
}