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

Example 31 with Rotamer

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

the class RotamerOptimization method currentEnergy.

/**
 * Calculates the energy at the current state.
 *
 * @param resList List of residues in current energy term.
 * @return Energy of the current state.
 */
private double currentEnergy(List<Residue> resList) throws ArithmeticException {
    List<Rotamer> rots = resList.stream().filter(res -> res != null).map(Residue::getRotamer).collect(Collectors.toList());
    File energyDir = dirSupplier.apply(resList, rots);
    return eFunction.applyAsDouble(energyDir);
}

Example 32 with Rotamer

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

the class RotamerOptimization method goldsteinElimination.

/**
 * Attemps to eliminate rotamer riA based on riB.
 *
 * @param residues
 * @param i
 * @param riA Rotamer to attempt elimination of.
 * @param riB Rotamer to attempt elimination by.
 * @return If riA was eliminated.
 */
private boolean goldsteinElimination(Residue[] residues, int i, int riA, int riB) {
    int nres = residues.length;
    Residue resi = residues[i];
    // Initialize Goldstein inequality.
    double selfDiff = getSelf(i, riA) - getSelf(i, riB);
    double goldsteinEnergy = selfDiff;
    double sumPairDiff = 0.0;
    double sumTripleDiff = 0.0;
    // Loop over a 2nd residue j.
    for (int j = 0; j < nres; j++) {
        if (j == i) {
            continue;
        }
        Residue resj = residues[j];
        Rotamer[] rotj = resj.getRotamers(library);
        int nrj = rotj.length;
        double minForResJ = Double.MAX_VALUE;
        double minPairDiff = 0.0;
        double minTripleDiff = 0.0;
        int rjEvals = 0;
        // Loop over the rotamers for residue j.
        for (int rj = 0; rj < nrj; rj++) {
            if (check(j, rj)) {
                continue;
            }
            if (check(i, riA, j, rj)) {
                // This is not a part of configuration space accessible to riA.
                continue;
            }
            if (check(i, riB, j, rj)) {
                /**
                 * This is a part of configuration space where riA is valid
                 * but not riB. Thus, if j,rj is part of the GMEC, riB is
                 * inconsistent with it. Thus, riB cannot be used to
                 * eliminate riA.
                 */
                return false;
            }
            double pairI = get2Body(i, riA, j, rj);
            double pairJ = get2Body(i, riB, j, rj);
            double pairDiff = pairI - pairJ;
            rjEvals++;
            // Include three-body interactions.
            double tripleDiff = 0.0;
            if (threeBodyTerm) {
                for (int k = 0; k < nres; k++) {
                    if (k == i || k == j) {
                        continue;
                    }
                    Residue resk = residues[k];
                    Rotamer[] rotk = resk.getRotamers(library);
                    int nrk = rotk.length;
                    int rkEvals = 0;
                    double minForResK = Double.MAX_VALUE;
                    for (int rk = 0; rk < nrk; rk++) {
                        /**
                         * If k,rk or j,rj-k,rk are not a part of valid
                         * configuration space, continue. If i,riA-k,rk or
                         * i,riA-j,rj-k,rk are not valid for riA, continue.
                         */
                        if (check(k, rk) || check(j, rj, k, rk) || check(i, riA, k, rk)) {
                            // Not yet implemented: check(i, riA, j, rj, k, rk) because no triples get eliminated.
                            continue;
                        }
                        /**
                         * If i,riB-k,rk or i,riB-j,rj-k,rk are invalid for
                         * riB, there is some part of configuration space
                         * for which riA is valid but not riB.
                         */
                        if (check(i, riB, k, rk)) {
                            // Not yet implemented: check(i, riB, j, rj, k, rk).
                            return false;
                        }
                        rkEvals++;
                        double e = get3Body(i, riA, j, rj, k, rk) - get3Body(i, riB, j, rj, k, rk);
                        if (e < minForResK) {
                            minForResK = e;
                        }
                    }
                    /**
                     * If there were no 3-body interactions with residue k,
                     * then minForResk is zero.
                     */
                    if (rkEvals == 0) {
                        minForResK = 0.0;
                    }
                    tripleDiff += minForResK;
                }
            }
            double sumOverK = pairDiff + tripleDiff;
            if (sumOverK < minForResJ) {
                minForResJ = sumOverK;
                minPairDiff = pairDiff;
                minTripleDiff = tripleDiff;
            }
        }
        // If there are no 2-body interactions, then minForResJ is zero.
        if (rjEvals == 0) {
            minForResJ = 0.0;
            minPairDiff = 0.0;
            minTripleDiff = 0.0;
        }
        sumPairDiff += minPairDiff;
        sumTripleDiff += minTripleDiff;
        goldsteinEnergy += minForResJ;
    }
    if (goldsteinEnergy > ensembleBuffer) {
        if (eliminateRotamer(residues, i, riA, print)) {
            logIfMaster(format("  Rotamer elimination of (%8s,%2d) by (%8s,%2d): %12.4f > %6.4f.", resi.toFormattedString(false, true), riA, resi.toFormattedString(false, true), riB, goldsteinEnergy, ensembleBuffer));
            logIfMaster(format("   Self: %12.4f, Pairs: %12.4f, Triples: %12.4f.", selfDiff, sumPairDiff, sumTripleDiff));
            return true;
        }
    }
    return false;
}

Example 33 with Rotamer

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

the class RotamerOptimization method minMaxE3.

/**
 * Calculates the minimum and maximum summations over additional residues
 * for some rotamer triples ri-rj-rk.
 *
 * @param residues Residues under consideration.
 * @param minMax Result array: 0 is min summation, 1 max summation.
 * @param i Residue i.
 * @param ri Rotamer for residue i.
 * @param j Residue j!=i.
 * @param rj Rotamer for residue j.
 * @param k Residue k!=j and k!=i.
 * @param rk Rotamer for residue k.
 * @return False if ri-rj-rk always clashes with other residues.
 * @throws IllegalArgumentException if there are pre-existing eliminations
 * in ri-rj-rk.
 */
private boolean minMaxE3(Residue[] residues, double[] minMax, int i, int ri, int j, int rj, int k, int rk) throws IllegalArgumentException {
    Residue resi = residues[i];
    Residue resj = residues[j];
    Residue resk = residues[k];
    if (check(i, ri) || check(j, rj) || check(k, rk) || check(i, ri, j, rj) || check(i, ri, k, rk) || check(j, rj, k, rk)) {
        // Not implemented: check(i, ri, j, rj, k, rk).
        throw new IllegalArgumentException(String.format(" Called for minMaxE2 on an eliminated triple %s-%d %s-%d %s-%d", resi.toFormattedString(false, true), ri, resj.toFormattedString(false, true), rj, resk.toFormattedString(false, true), rk));
    }
    // These two are a summation of mins/maxes over all fourth residues l.
    minMax[0] = 0;
    minMax[1] = 0;
    int nRes = residues.length;
    for (int l = 0; l < nRes; l++) {
        if (l == i || l == j || l == k) {
            continue;
        }
        Residue resl = residues[l];
        Rotamer[] rotsl = resl.getRotamers(library);
        int lenrl = rotsl.length;
        // Find min/max rl for residue l.
        double currentMax = Double.MIN_VALUE;
        double currentMin = Double.MAX_VALUE;
        for (int rl = 0; rl < lenrl; rl++) {
            if (check(l, rl) || check(k, rk, l, rl)) {
                // Not valid phase space for anything.
                continue;
            }
            double current;
            if (check(i, ri, l, rl) || check(j, rj, l, rl)) {
                // Not implemented: checking ri-rj-rl, ri-rk-rl, rj-rk-rl, or ri-rj-rk-rl.
                current = Double.NaN;
            } else {
                // ri-rj-rl is accounted for at a different part of the summation as ri-rj-rk.
                current = get3Body(i, ri, k, rk, l, rl) + get3Body(j, rj, k, rk, l, rl);
            }
            // minMaxE4(args)
            if (Double.isFinite(current) && current < currentMin) {
                // rl forms a more favorable 3-body than any prior rl for this residue l.
                currentMin = current;
            }
            if (Double.isFinite(current) && Double.isFinite(currentMax)) {
                if (current > currentMax) {
                    currentMax = current;
                }
            // Else, no new finite max found.
            } else {
                currentMax = Double.NaN;
            }
        }
        if (Double.isFinite(currentMin)) {
            minMax[0] += currentMin;
        } else {
            // Else, ri-rj-rk inevitably conflicts with l.
            minMax[0] = Double.NaN;
            minMax[1] = Double.NaN;
            return false;
        }
        if (Double.isFinite(currentMax) && Double.isFinite(minMax[1])) {
            minMax[1] += currentMax;
        } else {
            minMax[1] = Double.NaN;
        }
    // Finished with this residue l.
    }
    return Double.isFinite(minMax[0]);
}

Example 34 with Rotamer

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

the class RotamerOptimization method decomposedRotamerOptimization.

/**
 * Recursive brute-force method which uses single, pair, and potentially
 * trimer energies to calculate an optimum set of rotamers.
 *
 * @param molecularAssembly
 * @param residues Optimization window
 * @param i Current residue in the recursion.
 * @param lowEnergy Minimum energy yet found by the recursion.
 * @param optimum Optimum rotamer set yet found by the recursion.
 * @param currentRotamers Rotamer permutation under investigation.
 * @return Minimum energy found under this node in the recursion.
 */
private double decomposedRotamerOptimization(MolecularAssembly molecularAssembly, Residue[] residues, int i, double lowEnergy, int[] optimum, int[] currentRotamers) {
    // This is the initialization condition.
    if (i == 0) {
        evaluatedPermutations = 0;
    }
    int nResidues = residues.length;
    Residue current = residues[i];
    Rotamer[] rotamers = current.getRotamers(library);
    int lenri = rotamers.length;
    double currentEnergy = Double.MAX_VALUE;
    if (i < nResidues - 1) {
        /**
         * As long as there are more residues, continue the recursion for
         * each rotamer of the current residue.
         */
        for (int ri = 0; ri < lenri; ri++) {
            if (check(i, ri)) {
                continue;
            }
            currentRotamers[i] = ri;
            double rotEnergy = decomposedRotamerOptimization(molecularAssembly, residues, i + 1, lowEnergy, optimum, currentRotamers);
            if (rotEnergy < lowEnergy) {
                lowEnergy = rotEnergy;
            }
            if (rotEnergy < currentEnergy) {
                currentEnergy = rotEnergy;
            }
        }
    } else {
        /**
         * At the end of the recursion, compute the potential energy for
         * each rotamer of the final residue and update optimum[].
         */
        for (int ri = 0; ri < lenri; ri++) {
            if (check(i, ri)) {
                continue;
            }
            currentRotamers[i] = ri;
            double rotEnergy = computeEnergy(residues, currentRotamers, false);
            ++evaluatedPermutations;
            if (rotEnergy < currentEnergy) {
                currentEnergy = rotEnergy;
            }
            if (rotEnergy < lowEnergy) {
                /* Because we print the rotamer set immediately on finding a
                     * more optimal structure, we have to reset the entire length
                     * of optimum instead of lazily doing it on the way down.
                     */
                System.arraycopy(currentRotamers, 0, optimum, 0, optimum.length);
                if (evaluatedPermutations > 1) {
                    logger.info(format(" Minimum energy update: %f < %f, permutation %d", rotEnergy, lowEnergy, evaluatedPermutations));
                    String permutation = " Rotamer permutation: " + optimum[0];
                    for (int j = 1; j < nResidues; j++) {
                        permutation = permutation.concat(", " + optimum[j]);
                    }
                    logger.info(permutation);
                } else {
                    logger.info(format(" First minimum energy (permutation 1): %f", rotEnergy));
                }
                lowEnergy = rotEnergy;
            }
        }
    }
    return currentEnergy;
}

Example 35 with Rotamer

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

the class RotamerOptimization method pruneSingleClashes.

/**
 * Uses calculated energies to prune rotamers based on a threshold distance
 * from that residue's minimum energy rotamer (by default 20 kcal/mol). The
 * threshold can be modulated by presence of nucleic acids or MultiResidues,
 * which require more generous pruning criteria.
 *
 * @param residues Residues to prune rotamers over.
 */
private void pruneSingleClashes(Residue[] residues) {
    if (!pruneClashes) {
        return;
    }
    for (int i = 0; i < residues.length; i++) {
        Residue residue = residues[i];
        Rotamer[] rotamers = residue.getRotamers(library);
        int nrot = rotamers.length;
        double minEnergy = Double.MAX_VALUE;
        int minRot = -1;
        for (int ri = 0; ri < nrot; ri++) {
            if (!check(i, ri) && getSelf(i, ri) < minEnergy) {
                minEnergy = getSelf(i, ri);
                minRot = ri;
            }
        }
        /**
         * Regular: ep = minEnergy + clashThreshold Nucleic acids: ep =
         * minEnergy + (clashThreshold * factor * factor) MultiResidues: ep
         * = minEnergy + multiResClashThreshold
         *
         * Nucleic acids are bigger than amino acids, and MultiResidues can
         * have wild swings in energy on account of chemical perturbation.
         */
        double energyToPrune = (residue instanceof MultiResidue) ? multiResClashThreshold : clashThreshold;
        energyToPrune = (residue.getResidueType() == NA) ? energyToPrune * singletonNAPruningFactor * pruningFactor : energyToPrune;
        energyToPrune += minEnergy;
        for (int ri = 0; ri < nrot; ri++) {
            if (!check(i, ri) && (getSelf(i, ri) > energyToPrune)) {
                if (eliminateRotamer(residues, i, ri, print)) {
                    logIfMaster(format("  Rotamer (%7s,%2d) self-energy %s pruned by (%7s,%2d) %s.", residue, ri, formatEnergy(getSelf(i, ri)), residue, minRot, formatEnergy(minEnergy)));
                }
            }
        }
    }
}