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

Example 41 with Rotamer

use of ffx.potential.bonded.Rotamer 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 42 with Rotamer

use of ffx.potential.bonded.Rotamer 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 43 with Rotamer

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

Example 44 with Rotamer

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

the class RotamerOptimization method firstValidPerm.

/**
 * Finds the first non-eliminated rotamer permutation.
 *
 * @param residues
 * @param i
 * @param currentRotamers
 * @return If valid permutation found.
 */
private boolean firstValidPerm(Residue[] residues, int i, int[] currentRotamers) {
    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;
            if (firstValidPerm(residues, i + 1, currentRotamers)) {
                return true;
            }
        }
    } 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;
            }
            return true;
        }
    }
    return false;
}

Example 45 with Rotamer

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

the class RotamerOptimization method allocateEliminationMemory.

protected void allocateEliminationMemory(Residue[] residues) {
    int nres = residues.length;
    eliminatedSingles = new boolean[nres][];
    eliminatedPairs = new boolean[nres][][][];
    // Loop over residues.
    for (int i = 0; i < nres; i++) {
        Residue residuei = residues[i];
        Rotamer[] rotamersi = residuei.getRotamers(library);
        // Length rotamers i
        int lenri = rotamersi.length;
        logIfMaster(format(" %3d Residue %7s with %2d rotamers.", i + 1, residuei.toFormattedString(false, true), lenri));
        eliminatedSingles[i] = new boolean[lenri];
        eliminatedPairs[i] = new boolean[lenri][][];
        // Loop over the set of rotamers for residue i.
        for (int ri = 0; ri < lenri; ri++) {
            // int npairs = nres - (i + 1);
            // TODO - reduce memory by half.
            eliminatedSingles[i][ri] = false;
            eliminatedPairs[i][ri] = new boolean[nres][];
            for (int j = i + 1; j < nres; j++) {
                Residue residuej = residues[j];
                Rotamer[] rotamersj = residuej.getRotamers(library);
                int lenrj = rotamersj.length;
                eliminatedPairs[i][ri][j] = new boolean[lenrj];
                for (int rj = 0; rj < lenrj; rj++) {
                    eliminatedPairs[i][ri][j][rj] = false;
                }
            }
        }
    }
}