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

Example 26 with Residue

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

the class RotamerOptimization method printLargeInteractions.

/**
 * Prints a summary of pair and trimer energies above [cutoff] kcal/mol.
 *
 * @param pairCutoff
 * @param trimerCutoff
 * @param dOMode
 */
public void printLargeInteractions(double pairCutoff, double trimerCutoff, boolean dOMode) {
    Residue[] residues = residueList.toArray(new Residue[residueList.size()]);
    int nRes = residues.length;
    if (dOMode) {
        logger.info(format(" Large pair interactions (>%.2f):", pairCutoff));
        for (int i = 0; i < nRes; i++) {
            for (int j = i + 1; j < nRes; j++) {
                if (Math.abs(twoBodyEnergy[i][0][j][0]) >= pairCutoff) {
                    logger.info(format(" Large Pair %s %s:       %16.5f", residues[i].toFormattedString(false, true), residues[j].toFormattedString(false, true), twoBodyEnergy[i][0][j][0]));
                }
            }
        }
        logger.info(format("\n Large trimer interactions (>%.2f):", trimerCutoff));
        for (int i = 0; i < nRes; i++) {
            for (int j = i + 1; j < nRes; j++) {
                for (int k = j + 1; k < nRes; k++) {
                    if (Math.abs(threeBodyEnergy[i][0][j][0][k][0]) >= trimerCutoff) {
                        logger.info(format(" Large Trimer  %s %s %s:    %16.5f", residues[i].toFormattedString(false, true), residues[j].toFormattedString(false, true), residues[k].toFormattedString(false, true), threeBodyEnergy[i][0][j][0][k][0]));
                    }
                }
            }
        }
        return;
    }
    logger.info(format(" Large pair interactions (>%.2f):", pairCutoff));
    for (int i = 0; i < nRes; i++) {
        Residue resi = residues[i];
        Rotamer[] roti = resi.getRotamers(library);
        for (int ri = 0; ri < roti.length; ri++) {
            for (int j = i + 1; j < nRes; j++) {
                Residue resj = residues[j];
                Rotamer[] rotj = resj.getRotamers(library);
                for (int rj = 0; rj < rotj.length; rj++) {
                    try {
                        if (Math.abs(twoBodyEnergy[i][ri][j][rj]) >= pairCutoff) {
                            logger.info(format(" Large Pair %8s %-2d, %8s %-2d: %s", resi.toFormattedString(false, true), ri, resj.toFormattedString(false, true), rj, formatEnergy(twoBodyEnergy[i][ri][j][rj])));
                        }
                    } catch (Exception ex) {
                    }
                }
            }
        }
    }
    logger.info(format("\n Large trimer interactions (>%.2f):", trimerCutoff));
    for (int i = 0; i < nRes; i++) {
        Residue resi = residues[i];
        Rotamer[] roti = resi.getRotamers(library);
        for (int ri = 0; ri < roti.length; ri++) {
            for (int j = i + 1; j < nRes; j++) {
                Residue resj = residues[j];
                Rotamer[] rotj = resj.getRotamers(library);
                for (int rj = 0; rj < rotj.length; rj++) {
                    for (int k = j + 1; k < nRes; k++) {
                        Residue resk = residues[k];
                        Rotamer[] rotk = resk.getRotamers(library);
                        for (int rk = 0; rk < rotk.length; rk++) {
                            try {
                                if (Math.abs(threeBodyEnergy[i][ri][j][rj][k][rk]) >= trimerCutoff) {
                                    logger.info(format(" Large Trimer %8s %-2d, %8s %-2d, %8s %-2d: %s", resi.toFormattedString(false, true), ri, resj.toFormattedString(false, true), rj, resk.toFormattedString(false, true), rk, formatEnergy(threeBodyEnergy[i][ri][j][rj][k][rk])));
                                }
                            } catch (Exception ex) {
                            }
                        }
                    }
                }
            }
        }
    }
}

Example 27 with Residue

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

the class RotamerOptimization method sortResidues.

/**
 * Sorts a passed List of Residues by global index.
 *
 * @param residues List of Residues to be sorted.
 */
private void sortResidues(List<Residue> residues) {
    Comparator comparator = Comparator.comparing(Residue::getChainID).thenComparingInt(Residue::getResidueNumber);
    residues.sort(comparator);
}

Example 28 with Residue

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

the class RotamerOptimization method rotamerOptimizationDEE.

/**
 * A global optimization over side-chain rotamers using a recursive
 * algorithm and information about eliminated rotamers, rotamer pairs and
 * rotamer triples
 *
 * @param molecularAssembly
 * @param residues
 * @param i
 * @param currentRotamers
 * @param lowEnergy
 * @param optimum Optimum set of rotamers.
 * @param permutationEnergies Energies of visited permutations or null.
 * @return current energy.
 */
private double rotamerOptimizationDEE(MolecularAssembly molecularAssembly, Residue[] residues, int i, int[] currentRotamers, double lowEnergy, int[] optimum, double[] permutationEnergies) {
    // 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;
    double currentEnergy = Double.MAX_VALUE;
    List<Residue> resList = Arrays.asList(residues);
    /**
     * As long as there are more residues, continue the recursion for each
     * rotamer of the current residue.
     */
    if (i < nResidues - 1) {
        /**
         * Loop over rotamers of residue i.
         */
        for (int ri = 0; ri < lenri; ri++) {
            /**
             * Check if rotamer ri has been eliminated by DEE.
             */
            if (check(i, ri)) {
                continue;
            }
            /**
             * Check if rotamer ri has been eliminated by an upstream
             * rotamer (any residue's rotamer from j = 0 .. i-1).
             */
            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;
            }
            applyRotamer(residuei, rotamersi[ri]);
            currentRotamers[i] = ri;
            double rotEnergy = rotamerOptimizationDEE(molecularAssembly, residues, i + 1, currentRotamers, lowEnergy, optimum, permutationEnergies);
            if (rotEnergy < currentEnergy) {
                currentEnergy = rotEnergy;
            }
            if (rotEnergy < lowEnergy) {
                optimum[i] = ri;
                lowEnergy = rotEnergy;
            }
        }
    } else {
        if (ensembleStates == null) {
            ensembleStates = new ArrayList<>();
        }
        /**
         * At the end of the recursion, compute the potential energy for
         * each rotamer of the final residue. If a lower potential energy is
         * discovered, the rotamers of each residue will be collected as the
         * recursion returns up the chain.
         */
        for (int ri = 0; ri < lenri; ri++) {
            /**
             * Check if rotamer ri has been eliminated by DEE.
             */
            if (check(i, ri)) {
                continue;
            }
            currentRotamers[i] = ri;
            /**
             * Check if rotamer ri has been eliminated by an upstream
             * rotamer (any residue's rotamer from 0 .. i-1.
             */
            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;
            }
            applyRotamer(residuei, rotamersi[ri]);
            // Compute the energy based on a 3-body approximation
            double approximateEnergy = computeEnergy(residues, currentRotamers, false);
            double comparisonEnergy = approximateEnergy;
            evaluatedPermutations++;
            // Compute the AMOEBA energy
            if (useFullAMOEBAEnergy) {
                double amoebaEnergy = Double.NaN;
                try {
                    amoebaEnergy = currentEnergy(resList);
                } catch (ArithmeticException ex) {
                    logger.warning(String.format(" Exception %s in calculating full AMOEBA energy for permutation %d", ex.toString(), evaluatedPermutations));
                }
                comparisonEnergy = amoebaEnergy;
            }
            if (permutationEnergies != null) {
                permutationEnergies[evaluatedPermutations - 1] = comparisonEnergy;
            }
            if (algorithmListener != null) {
                algorithmListener.algorithmUpdate(molecularAssembly);
            }
            if (ensembleNumber > 1) {
                if (master && printFiles) {
                    try {
                        FileWriter fw = new FileWriter(ensembleFile, true);
                        BufferedWriter bw = new BufferedWriter(fw);
                        bw.write(format("MODEL        %d", evaluatedPermutations));
                        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 e) {
                        logger.warning(format("Exception writing to file: %s", ensembleFile.getName()));
                    }
                }
                ResidueState[] states = ResidueState.storeAllCoordinates(residues);
                ensembleStates.add(new ObjectPair<>(states, comparisonEnergy));
            }
            if (comparisonEnergy < currentEnergy) {
                currentEnergy = comparisonEnergy;
            }
            if (comparisonEnergy < lowEnergy) {
                lowEnergy = comparisonEnergy;
                optimum[i] = ri;
            }
            if (useFullAMOEBAEnergy) {
                // Log current results
                logIfMaster(format(" %6e AMOEBA: %12.4f 3-Body: %12.4f Neglected: %12.4f (%12.4f)", (double) evaluatedPermutations, comparisonEnergy, approximateEnergy, comparisonEnergy - approximateEnergy, lowEnergy));
            } else {
                if (threeBodyTerm) {
                    ;
                    logIfMaster(format(" %12s %5s %25f %5s %25f %5s", evaluatedPermutations, "|", approximateEnergy, "|", lowEnergy, "|"));
                // logIfMaster(format(" %6e 3-Body: %12.4f (%12.4f)",
                // (double) evaluatedPermutations, approximateEnergy, lowEnergy));
                } else {
                    logIfMaster(format(" %12s %5s %25f %5s %25f %5s", evaluatedPermutations, "|", approximateEnergy, "|", lowEnergy, "|"));
                // logIfMaster(format(" %6e 2-Body: %12.4f (%12.4f)",
                // (double) evaluatedPermutations, approximateEnergy, lowEnergy));
                }
            }
        }
        ensembleStates.sort(null);
    }
    return currentEnergy;
}

Example 29 with Residue

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

the class RotamerOptimization method reconcileNARotamersWithPriorResidues.

/**
 * For NA residues inside some optimization window, prune any rotamers which
 * would be incompatible with the established rotamers of upstream NA
 * residues. Could in theory be done by self energies, but every rotamer
 * which can be eliminated without calculating a self energy makes the
 * optimization significantly faster.
 *
 * @param residues Residues to check for incompatible rotamers.
 * @return Number of rotamers eliminated for each Residue.
 */
private int[] reconcileNARotamersWithPriorResidues(Residue[] residues) {
    int[] numEliminatedRotamers = new int[residues.length];
    for (int i = 0; i < residues.length; i++) {
        Residue residuei = residues[i];
        Rotamer[] rotamers = residuei.getRotamers(library);
        if (rotamers == null || residuei.getResidueType() != NA) {
            continue;
        }
        Residue prevResidue = residuei.getPreviousResidue();
        if (prevResidue == null || prevResidue.getResidueType() != NA) {
            continue;
        }
        boolean isInList = false;
        for (Residue residue : residues) {
            if (prevResidue.equals(residue)) {
                isInList = true;
                break;
            }
        }
        if (isInList) {
            continue;
        }
        double prevDelta = RotamerLibrary.measureDelta(prevResidue);
        // If between 50 and 110, assume a North pucker.
        if (RotamerLibrary.checkPucker(prevDelta) == 1) {
            for (int j = 0; j < rotamers.length; j++) {
                if (RotamerLibrary.checkPucker(rotamers[j].chi1) != 1) {
                    if (print) {
                        logIfMaster(format(" Rotamer %d of residue %s eliminated " + "for incompatibility with the sugar pucker of previous " + "residue %s outside the window.", j, residuei.toFormattedString(false, true), prevResidue.toString()));
                    }
                    eliminateRotamer(residues, i, j, print);
                    numEliminatedRotamers[i]++;
                }
            }
        } else {
            for (int j = 0; j < rotamers.length; j++) {
                if (RotamerLibrary.checkPucker(rotamers[j].chi1) != 2) {
                    if (print) {
                        logIfMaster(format(" Rotamer %d of residue %s eliminated " + "for incompatibility with the sugar pucker of previous " + "residue %s outside the window.", j, residuei.toFormattedString(false, true), prevResidue.toString()));
                    }
                    eliminateRotamer(residues, i, j, print);
                    numEliminatedRotamers[i]++;
                }
            }
        }
    // TODO: Implement support for the DNA C3'-exo pucker.
    }
    return numEliminatedRotamers;
}

Example 30 with Residue

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

the class RotamerOptimization method prunePairClashes.

/**
 * Prunes rotamer ri of residue i if all ri-j pair energies are worse than
 * the best i-j pair by some threshold value; additionally prunes ri-rj
 * pairs if they exceed the best i-j pair by a greater threshold value;
 * additionally performs this in reverse (searches over j-i).
 *
 * @param residues Residues whose rotamers are to be pruned.
 */
private void prunePairClashes(Residue[] residues) {
    if (!prunePairClashes) {
        return;
    }
    int nResidues = residues.length;
    for (int i = 0; i < nResidues - 1; i++) {
        Residue resi = residues[i];
        Rotamer[] rotsi = resi.getRotamers(library);
        int lenri = rotsi.length;
        for (int j = i + 1; j < nResidues; j++) {
            Residue resj = residues[j];
            Rotamer[] rotsj = resj.getRotamers(library);
            int lenrj = rotsj.length;
            double minPair = Double.MAX_VALUE;
            int minRI = -1;
            int minRJ = -1;
            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;
                    }
                    double pairEnergy = get2Body(i, ri, j, rj) + getSelf(i, ri) + getSelf(j, rj);
                    assert Double.isFinite(pairEnergy);
                    if (pairEnergy < minPair) {
                        minPair = pairEnergy;
                        minRI = ri;
                        minRJ = rj;
                    }
                }
            }
            // Otherwise, i and j are not on a well-packed backbone.
            assert (minRI >= 0 && minRJ >= 0);
            // Calculate all the modifiers to the pair clash elimination threshold.
            double threshold = pairClashThreshold;
            if (resi instanceof MultiResidue) {
                threshold += multiResPairClashAddn;
            }
            if (resj instanceof MultiResidue) {
                threshold += multiResPairClashAddn;
            }
            int numNARes = (resi.getResidueType() == NA ? 1 : 0) + (resj.getResidueType() == NA ? 1 : 0);
            switch(numNARes) {
                case 0:
                    break;
                case 1:
                    threshold *= pairHalfPruningFactor;
                    break;
                case 2:
                    threshold *= pruningFactor;
                    break;
                default:
                    throw new ArithmeticException(" RotamerOptimization.prunePairClashes() has somehow " + "found less than zero or more than two nucleic acid residues in a pair of" + " residues. This result should be impossible.");
            }
            double toEliminate = threshold + minPair;
            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;
                    }
                    double pairEnergy = get2Body(i, ri, j, rj) + getSelf(i, ri) + getSelf(j, rj);
                    assert Double.isFinite(pairEnergy);
                    if (pairEnergy > toEliminate) {
                        logIfMaster(format(" Pruning pair %s-%d %s-%d by %s-%d %s-%d; energy %s > " + "%s + %s", resi.toFormattedString(false, true), ri, resj.toFormattedString(false, true), rj, resi.toFormattedString(false, true), minRI, resj.toFormattedString(false, true), minRJ, formatEnergy(pairEnergy), formatEnergy(threshold), formatEnergy(minPair)));
                    }
                }
            }
            pairsToSingleElimination(residues, i, j);
        }
    }
}