Examples with RotamerLibrary.applyRotamer ffx.potential.bonded.RotamerLibrary.applyRotamer used on opensource projects

Example 1 with RotamerLibrary.applyRotamer

use of ffx.potential.bonded.RotamerLibrary.applyRotamer in project ffx by mjschnie.

the class RotamerOptimization method bruteForce.

/**
 * Performs a recursive brute-force rotamer optimization over a passed list
 * of residues.
 *
 * @param residueList Residues to be optimized.
 * @return Global minimum energy conformation energy.
 */
private double bruteForce(List<Residue> residueList) {
    Residue[] residues = residueList.toArray(new Residue[residueList.size()]);
    int nResidues = residues.length;
    int[] optimum = new int[nResidues];
    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.
     */
    double permutations = 1;
    for (int i = 0; i < nResidues; i++) {
        Residue residue = residues[i];
        Rotamer[] rotamers = residue.getRotamers(library);
        permutations *= rotamers.length;
    }
    logger.info(format(" Number of permutations: %16.8e.", permutations));
    double e;
    useFullAMOEBAEnergy = false;
    if (!useFullAMOEBAEnergy) {
        setPruning(0);
        rotamerEnergies(residues);
        int[] rotamerSet = new int[nResidues];
        fill(rotamerSet, 0);
        e = decomposedRotamerOptimization(molecularAssembly, residues, 0, Double.MAX_VALUE, optimum, rotamerSet);
        for (int i = 0; i < nResidues; i++) {
            Residue residue = residues[i];
            Rotamer[] rotamers = residue.getRotamers(library);
            RotamerLibrary.applyRotamer(residue, rotamers[optimum[i]]);
            turnOnAtoms(residue);
        }
        double fullEnergy = 0;
        try {
            fullEnergy = currentEnergy(residueList);
        } catch (Exception ex) {
            logger.severe(String.format(" Exception %s in calculating full energy; FFX shutting down", ex.toString()));
        }
        logger.info(format(" Final summation of energies:    %16.5f", e));
        logger.info(format(" Final energy of optimized structure:    %16.5f", fullEnergy));
        logger.info(format(" Neglected:    %16.5f", fullEnergy - e));
    } else {
        e = rotamerOptimization(molecularAssembly, residues, 0, Double.MAX_VALUE, optimum);
    }
    for (int i = 0; i < nResidues; i++) {
        Residue residue = residues[i];
        Rotamer[] rotamers = residue.getRotamers(library);
        int ri = optimum[i];
        Rotamer rotamer = rotamers[ri];
        logger.info(format(" %s %s (%d)", residue.getResidueNumber(), rotamer.toString(), ri));
        RotamerLibrary.applyRotamer(residue, rotamer);
        if (useFullAMOEBAEnergy) {
            try {
                e = currentEnergy(residueList);
            } catch (ArithmeticException ex) {
                logger.fine(String.format(" Exception %s in calculating full AMOEBA energy at the end of brute force", ex.toString()));
            }
        }
    }
    return e;
}

Example 2 with RotamerLibrary.applyRotamer

use of ffx.potential.bonded.RotamerLibrary.applyRotamer in project ffx by mjschnie.

the class RotamerOptimization method independent.

private double independent(List<Residue> residues) {
    if (x == null) {
        Atom[] atoms = molecularAssembly.getAtomArray();
        int nAtoms = atoms.length;
        x = new double[nAtoms * 3];
    }
    double e = Double.MAX_VALUE;
    List<Residue> rList = new ArrayList<>(Collections.nCopies(1, null));
    for (Residue residue : residues) {
        rList.set(0, residue);
        logger.info(format(" Optimizing %s side-chain.", residue));
        Rotamer[] rotamers = residue.getRotamers(library);
        potential.getCoordinates(x);
        e = Double.MAX_VALUE;
        int bestRotamer = -1;
        for (int j = 0; j < rotamers.length; j++) {
            Rotamer rotamer = rotamers[j];
            RotamerLibrary.applyRotamer(residue, rotamer);
            if (algorithmListener != null) {
                algorithmListener.algorithmUpdate(molecularAssembly);
            }
            double newE = Double.NaN;
            try {
                newE = currentEnergy(rList);
            } catch (ArithmeticException ex) {
                logger.fine(String.format(" Exception %s in energy calculations during independent for %s-%d", ex.toString(), residue, j));
            }
            if (newE < e) {
                bestRotamer = j;
            }
        }
        if (bestRotamer > -1) {
            Rotamer rotamer = rotamers[bestRotamer];
            RotamerLibrary.applyRotamer(residue, rotamer);
        }
        if (algorithmListener != null) {
            algorithmListener.algorithmUpdate(molecularAssembly);
        }
        if (terminate) {
            logger.info(format("\n Terminating after residue %s.\n", residue));
            break;
        }
    }
    return e;
}

Example 3 with RotamerLibrary.applyRotamer

use of ffx.potential.bonded.RotamerLibrary.applyRotamer in project ffx by mjschnie.

the class RotamerOptimization method evaluateDistance.

/**
 * Evaluates the pairwise distance between two residues' rotamers under any
 * symmetry operator; does "lazy loading" for the distance matrix.
 *
 * @param i Residue i
 * @param ri Rotamer for i
 * @param j Residue j
 * @param rj Rotamer for j
 * @return Shortest distance
 */
private double evaluateDistance(int i, int ri, int j, int rj) {
    Residue resi = allResiduesArray[i];
    Rotamer[] rotamersI = resi.getRotamers(library);
    Rotamer roti = rotamersI[ri];
    double[][] xi;
    if (roti.equals(resi.getRotamer())) {
        xi = resi.storeCoordinateArray();
    } else {
        ResidueState origI = resi.storeState();
        RotamerLibrary.applyRotamer(resi, roti);
        xi = resi.storeCoordinateArray();
        resi.revertState(origI);
    }
    Residue resj = allResiduesArray[j];
    Rotamer[] rotamersJ = resj.getRotamers(library);
    Rotamer rotj = rotamersJ[rj];
    double[][] xj;
    if (rotj.equals(resj.getRotamer())) {
        xj = resj.storeCoordinateArray();
    } else {
        ResidueState origJ = resj.storeState();
        RotamerLibrary.applyRotamer(resj, rotj);
        xj = resj.storeCoordinateArray();
        resj.revertState(origJ);
    }
    Crystal crystal = molecularAssembly.getCrystal();
    int nSymm = crystal.spaceGroup.getNumberOfSymOps();
    double minDist = Double.MAX_VALUE;
    for (int iSymOp = 0; iSymOp < nSymm; iSymOp++) {
        SymOp symOp = crystal.spaceGroup.getSymOp(iSymOp);
        double dist = interResidueDistance(xi, xj, symOp);
        minDist = dist < minDist ? dist : minDist;
    }
    return minDist;
}

Example 4 with RotamerLibrary.applyRotamer

use of ffx.potential.bonded.RotamerLibrary.applyRotamer in project ffx by mjschnie.

the class RotamerOptimization method setResidues.

/**
 * Set a contiguous block of residues to optimize.
 *
 * @param startResID
 * @param finalResID
 */
public void setResidues(int startResID, int finalResID) {
    Polymer polymer;
    if (chain != null) {
        polymer = molecularAssembly.getChain(chain);
    } else {
        polymers = molecularAssembly.getChains();
        polymer = polymers[0];
    }
    residueList = new ArrayList<>();
    for (int i = startResID; i <= finalResID; i++) {
        Residue residue = polymer.getResidue(i);
        if (residue != null) {
            Rotamer[] rotamers = residue.getRotamers(library);
            if (rotamers != null) {
                if (rotamers.length == 1) {
                    switch(residue.getResidueType()) {
                        case NA:
                            residue.initializeDefaultAtomicCoordinates();
                            break;
                        case AA:
                        default:
                            RotamerLibrary.applyRotamer(residue, rotamers[0]);
                            break;
                    }
                    if (addOrigRot) {
                        residueList.add(residue);
                    }
                } else {
                    residueList.add(residue);
                }
            } else if (useForcedResidues && checkIfForced(i)) {
                residueList.add(residue);
            }
        }
    }
}

Example 5 with RotamerLibrary.applyRotamer

use of ffx.potential.bonded.RotamerLibrary.applyRotamer 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;
}