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

Example 11 with Residue

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

the class RotamerOptimization method goldsteinPairDriver.

/**
 * Finds and eliminates rotamer pairs according to the many-body Goldstein
 * pairs criterion.
 *
 * @param residues Residues under consideration.
 * @return If any rotamer pairs were eliminated.
 */
private boolean goldsteinPairDriver(Residue[] residues) {
    int nRes = residues.length;
    boolean eliminated = false;
    // First, generate pairs riA-rjC.
    for (int i = 0; i < nRes; i++) {
        Residue resi = residues[i];
        Rotamer[] rotsi = resi.getRotamers(library);
        int lenri = rotsi.length;
        for (int riA = 0; riA < lenri; riA++) {
            // Don't try to eliminate that which is already eliminated.
            if (check(i, riA)) {
                continue;
            }
            // Residue j can be any other residue, including ones before residue i.
            for (int j = 0; j < nRes; j++) {
                // Residue j must be distinct from i.
                if (i == j) {
                    continue;
                }
                Residue resj = residues[j];
                Rotamer[] rotsj = resj.getRotamers(library);
                int lenrj = rotsj.length;
                for (int rjC = 0; rjC < lenrj; rjC++) {
                    // Again, no point in eliminating the already-eliminated.
                    if (check(j, rjC) || check(i, riA, j, rjC)) {
                        continue;
                    }
                    boolean breakOut = false;
                    // Now, generate pairs riB-rjD. If any pair riB-rjD eliminates riA-rjC, break out of the loop.
                    for (int riB = 0; riB < lenri; riB++) {
                        if (breakOut) {
                            break;
                        }
                        if (check(i, riB)) {
                            continue;
                        }
                        for (int rjD = 0; rjD < lenrj; rjD++) {
                            if (breakOut) {
                                break;
                            }
                            // Do not attempt eliminating with an eliminated pair.
                            if (check(j, rjD) || check(i, riB, j, rjD)) {
                                continue;
                            }
                            // Do not attempt to eliminate a pair with itself.
                            if (riA == riB && rjC == rjD) {
                                continue;
                            }
                            if (goldsteinPairElimination(residues, i, riA, riB, j, rjC, rjD)) {
                                breakOut = true;
                                eliminated = true;
                            }
                        }
                    }
                }
                if (pairsToSingleElimination(residues, i, j)) {
                    eliminated = true;
                }
            }
        }
    }
    return eliminated;
}

Example 12 with Residue

use of ffx.potential.bonded.Residue 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 13 with Residue

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

the class RotamerOptimization method minMaxE2.

/**
 * Calculates the minimum and maximum summations over additional residues
 * for some pair ri-rj.
 *
 * @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.
 * @return False if ri-rj always clashes with other residues.
 * @throws IllegalArgumentException If ri, rj, or ri-rj eliminated.
 */
private boolean minMaxE2(Residue[] residues, double[] minMax, int i, int ri, int j, int rj) throws IllegalArgumentException {
    Residue resi = residues[i];
    Residue resj = residues[j];
    if (check(i, ri) || check(j, rj) || check(i, ri, j, rj)) {
        throw new IllegalArgumentException(String.format(" Called for minMaxE2 on an eliminated pair %s-%d %s-%d", resi.toFormattedString(false, true), ri, resj.toFormattedString(false, true), rj));
    }
    // Minimum summation over third residues k.
    minMax[0] = 0;
    // Maximum summation over third residues k.
    minMax[1] = 0;
    int nRes = residues.length;
    for (int k = 0; k < nRes; k++) {
        if (k == i || k == j) {
            continue;
        }
        Residue resk = residues[k];
        Rotamer[] rotsk = resk.getRotamers(library);
        int lenrk = rotsk.length;
        double[] minMaxK = new double[2];
        minMaxK[0] = Double.MAX_VALUE;
        minMaxK[1] = Double.MIN_VALUE;
        for (int rk = 0; rk < lenrk; rk++) {
            if (check(k, rk)) {
                // Not a valid part of phase space.
                continue;
            }
            if (check(i, ri, k, rk) || check(j, rj, k, rk)) {
                // Not implemented: check(i, ri, j, rj, k, rk).
                // i,ri or j,rj clashes with this rotamer, max will be NaN.
                // Minimum for this rk will be a clash, which is never a minimum.
                minMaxK[1] = Double.NaN;
            } else {
                // Min and max summations over 4th residues l, plus the ri-rk and rj-rk interactions.
                // If no 3-body term, just the ri-rk and rj-rk interactions.
                double currentMin = get2Body(i, ri, k, rk) + get2Body(j, rj, k, rk);
                double currentMax = currentMin;
                if (threeBodyTerm) {
                    // If the 3-Body eliminated, would fill max to Double.NaN.
                    currentMin += get3Body(i, ri, j, rj, k, rk);
                    currentMax = currentMin;
                    // Obtain min and max summations over l.
                    double[] minMaxTriple = new double[2];
                    if (minMaxE3(residues, minMaxTriple, i, ri, j, rj, k, rk)) {
                        // A non-finite triples minimum should have the code taking the else branch.
                        assert (Double.isFinite(minMaxTriple[0]) && minMaxTriple[0] != Double.MAX_VALUE);
                        // Add the min and max summations over all 4th residues l.
                        currentMin += minMaxTriple[0];
                        if (Double.isFinite(currentMax) && Double.isFinite(minMaxTriple[1])) {
                            currentMax += minMaxTriple[1];
                        } else {
                            currentMax = Double.NaN;
                        }
                    } else {
                        // i, ri, j, rj, k, rk creates an inevitable clash with some residue l.
                        currentMin = Double.NaN;
                        currentMax = Double.NaN;
                    }
                }
                assert (threeBodyTerm || currentMax == currentMin);
                // Now check if rk displaces previously searched rk for min/max over this k.
                if (Double.isFinite(currentMin) && currentMin < minMaxK[0]) {
                    // rk has a more favorable minimum than previously searched rk.
                    minMaxK[0] = currentMin;
                }
                if (Double.isFinite(currentMax) && Double.isFinite(minMaxK[1])) {
                    // rk has a less favorable maximum than previously searched rk.
                    minMaxK[1] = (currentMax > minMaxK[1]) ? currentMax : minMaxK[1];
                } else {
                    // Our maximum is a NaN.
                    minMaxK[1] = Double.NaN;
                }
            }
        }
        if (Double.isFinite(minMaxK[0])) {
            // Add the minimum contribution from this k to the summation.
            minMax[0] += minMaxK[0];
        } else {
            // Else, ri-rj conflicts with all rk for this k, and can be swiftly eliminated.
            minMax[0] = Double.NaN;
            minMax[1] = Double.NaN;
            return false;
        }
        if (Double.isFinite(minMaxK[1]) && Double.isFinite(minMax[1])) {
            // Add the max contribution from this k to the summation.
            minMax[1] += minMaxK[1];
        } else {
            // Otherwise, the max for ri-rj is a clash.
            minMax[1] = Double.NaN;
        }
    }
    return Double.isFinite(minMax[0]);
}

Example 14 with Residue

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

the class RotamerOptimization method goldsteinPairSumOverK.

private double goldsteinPairSumOverK(Residue[] residues, int lb, int ub, int i, int riA, int riB, int j, int rjC, int rjD, ArrayList<Residue> blockedResidues) {
    double sumOverK = 0.0;
    int nres = residues.length;
    for (int k = lb; k <= ub; k++) {
        if (k == j || k == i) {
            continue;
        }
        double minForResK = Double.MAX_VALUE;
        Residue resk = residues[k];
        Rotamer[] rotk = resk.getRotamers(library);
        int nrk = rotk.length;
        int rkEvals = 0;
        // Loop over residue k's rotamers.
        for (int rk = 0; rk < nrk; rk++) {
            if (check(k, rk)) {
                continue;
            }
            // Continue if k,rk invalid for riA/rjC.
            if (check(i, riA, k, rk) || check(j, rjC, k, rk)) {
                // Not implemented: check(i, riA, j, rjC, k, rk).
                continue;
            }
            // Return false if k,rk invalid for riB/rjD.
            if (check(i, riB, k, rk) || check(j, rjD, k, rk)) {
                blockedResidues.add(resk);
                return Double.NaN;
            }
            rkEvals++;
            double currentResK = get2Body(i, riA, k, rk) - get2Body(i, riB, k, rk) + get2Body(j, rjC, k, rk) - get2Body(j, rjD, k, rk);
            // Include 3-body effects.
            if (threeBodyTerm) {
                double sumOverL = (get3Body(i, riA, j, rjC, k, rk) - get3Body(i, riB, j, rjD, k, rk));
                // Loop over a 4th residue l.
                for (int l = 0; l < nres; l++) {
                    if (l == k || l == i || l == j) {
                        continue;
                    }
                    Residue residuel = residues[l];
                    Rotamer[] rotamersl = residuel.getRotamers(library);
                    int nrl = rotamersl.length;
                    int rlEvaluations = 0;
                    double minForResL = Double.MAX_VALUE;
                    // Loop over rotamers for residue l.
                    for (int rl = 0; rl < nrl; rl++) {
                        // If not a part of valid phase space for riA/rjC, continue.
                        if (check(l, rl) || check(k, rk, l, rl) || check(i, riA, l, rl) || check(j, rjC, l, rl)) {
                            // Not implemented: check(i, riA, j, rjC, l, rl) || check(i, riA, k, rk, l, rl) || check(j, rjC, k, rk, l, rl) || check(i, riA, j, rjC, k, rk, l, rl)
                            continue;
                        }
                        if (check(i, riB, l, rl) || check(j, rjD, l, rl)) {
                            // Not implemented: check(i, riB, j, rjD, l, rl) || check(i, riB, k, rk, l, rl) || check(j, rjD, k, rk, l, rl) || check(i, riB, j, rjD, k, rk, l, rl)
                            blockedResidues.add(residuel);
                            return Double.NaN;
                        }
                        rlEvaluations++;
                        double e = get3Body(i, riA, k, rk, l, rl) - get3Body(i, riB, k, rk, l, rl) + get3Body(j, rjC, k, rk, l, rl) - get3Body(j, rjD, k, rk, l, rl);
                        if (e < minForResL) {
                            minForResL = e;
                        }
                    }
                    if (rlEvaluations == 0) {
                        minForResL = 0.0;
                    }
                    sumOverL += minForResL;
                }
                currentResK += sumOverL;
            }
            if (currentResK < minForResK) {
                minForResK = currentResK;
            }
        }
        if (rkEvals == 0) {
            minForResK = 0.0;
            blockedResidues.add(resk);
        }
        sumOverK += minForResK;
    }
    return sumOverK;
}

Example 15 with Residue

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

the class RotamerOptimization method generateSuperbox.

/**
 * Returns the superbox used to generate the boxes for sliding box. If
 * superbox coordinates manually set, uses them plus the defined buffer.
 * Else, if an aperiodic system, uses maximum and minimum C alpha (or N1/9)
 * coordinates plus superboxBuffer (by default, 8A, longer than a lysine
 * side chain or N1/N9 distance to any other atom). Else, it just uses the
 * ordinary crystal.
 *
 * @param residueList List of residues to incorporate.
 * @return Superbox crystal.
 */
private Crystal generateSuperbox(List<Residue> residueList) {
    double[] maxXYZ = new double[3];
    double[] minXYZ = new double[3];
    Crystal originalCrystal = molecularAssembly.getCrystal();
    if (manualSuperbox) {
        for (int i = 0; i < maxXYZ.length; i++) {
            int ii = 2 * i;
            minXYZ[i] = boxDimensions[ii] - superboxBuffer;
            maxXYZ[i] = boxDimensions[ii + 1] + superboxBuffer;
        }
    } else if (originalCrystal.aperiodic()) {
        if (residueList == null || residueList.isEmpty()) {
            throw new IllegalArgumentException(" Null or empty residue list when generating superbox.");
        }
        Atom initializerAtom = residueList.get(0).getReferenceAtom();
        initializerAtom.getXYZ(minXYZ);
        initializerAtom.getXYZ(maxXYZ);
        for (Residue residue : residueList) {
            Atom refAtom = residue.getReferenceAtom();
            double[] refAtomCoords = new double[3];
            refAtom.getXYZ(refAtomCoords);
            for (int i = 0; i < 3; i++) {
                maxXYZ[i] = (refAtomCoords[i] > maxXYZ[i] ? refAtomCoords[i] : maxXYZ[i]);
                minXYZ[i] = (refAtomCoords[i] < minXYZ[i] ? refAtomCoords[i] : minXYZ[i]);
            }
        }
        for (int i = 0; i < 3; i++) {
            minXYZ[i] -= superboxBuffer;
            maxXYZ[i] += superboxBuffer;
        }
    } else {
        return originalCrystal;
    }
    double newA = maxXYZ[0] - minXYZ[0];
    double newB = maxXYZ[1] - minXYZ[1];
    double newC = maxXYZ[2] - minXYZ[2];
    if (manualSuperbox) {
        logger.info(format(" Manual superbox set over (minX, maxX, minY, " + "maxY, minZ, maxZ): %f, %f, %f, %f, %f, %f", minXYZ[0], maxXYZ[0], minXYZ[1], maxXYZ[1], minXYZ[2], maxXYZ[2]));
        logger.info(format(" Buffer size (included in dimensions): %f\n", superboxBuffer));
    } else {
        // Crystal systems will have already returned.
        logger.info(" System is aperiodic: protein box generated over these coordinates (minX, maxX, minY, maxY, minZ, maxZ):");
        String message = " Aperiodic box dimensions: ";
        for (int i = 0; i < minXYZ.length; i++) {
            message = message.concat(format("%f,%f,", minXYZ[i], maxXYZ[i]));
        }
        message = message.substring(0, message.length() - 1);
        logger.info(message);
        logger.info(format(" Buffer size (included in dimensions): %f\n", superboxBuffer));
    }
    return new Crystal(newA, newB, newC, 90.0, 90.0, 90.0, "P1");
}