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

Example 76 with Residue

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

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

the class RotamerOptimization method distanceMatrix.

/**
 * Calculates a residue-residue distance matrix.
 * <p>
 * Residue-residue distance is defined as the shortest atom-atom distance in
 * any possible rotamer-rotamer pair if the residues are neighbors (central
 * atom-central atom distances are within a cutoff). Otherewise, distances
 * are set to a default of Double.MAX_VALUE.
 * </p>
 * <p>
 * The intent of using a neighbor list is to avoid tediously searching
 * rotamer- rotamer pairs when two residues are so far apart we will never
 * need the exact distance. We use the distance matrix for adding residues
 * to the sliding window and determining whether to set 3-body energy to
 * 0.0.
 * </p>
 * <p>
 * If the central atoms are too distant from each other, we can safely
 * assume no atoms will ever be close enough for addition to sliding window
 * or to cause explicit calculation of 3-body energy.
 * </p>
 */
private void distanceMatrix() {
    distanceMatrix = new double[numResidues - 1][][][];
    long numDistances = 0L;
    for (int i = 0; i < (numResidues - 1); i++) {
        Residue residuei = allResiduesArray[i];
        int lengthRi;
        try {
            if (checkIfForced(residuei)) {
                lengthRi = 1;
            } else {
                lengthRi = residuei.getRotamers(library).length;
            }
        } catch (IndexOutOfBoundsException ex) {
            if (useForcedResidues) {
                logger.warning(ex.toString());
            } else {
                logIfMaster(format(" Non-forced Residue i %s has null rotamers.", residuei.toFormattedString(false, true)), Level.WARNING);
            }
            continue;
        }
        distanceMatrix[i] = new double[lengthRi][][];
        for (int ri = 0; ri < lengthRi; ri++) {
            distanceMatrix[i][ri] = new double[numResidues][];
            for (int j = (i + 1); j < numResidues; j++) {
                Residue residuej = allResiduesArray[j];
                int lengthRj;
                try {
                    if (checkIfForced(residuej)) {
                        lengthRj = 1;
                    } else {
                        lengthRj = residuej.getRotamers(library).length;
                    }
                } catch (IndexOutOfBoundsException ex) {
                    if (useForcedResidues) {
                        logger.warning(ex.toString());
                    } else {
                        logIfMaster(format(" Residue j %s has null rotamers.", residuej.toFormattedString(false, true)));
                    }
                    continue;
                }
                distanceMatrix[i][ri][j] = new double[lengthRj];
                numDistances += lengthRj;
                if (!lazyMatrix) {
                    fill(distanceMatrix[i][ri][j], Double.MAX_VALUE);
                } else {
                    fill(distanceMatrix[i][ri][j], -1.0);
                }
            }
        }
    }
    logger.info(format(" Number of pairwise distances: %d", numDistances));
    if (!lazyMatrix) {
        ResidueState[] orig = ResidueState.storeAllCoordinates(allResiduesList);
        int nMultiRes = 0;
        /**
         * Build a list that contains one atom from each Residues: CA from
         * amino acids, C1 from nucleic acids, or the first atom otherwise.
         */
        Atom[] atoms = new Atom[numResidues];
        for (int i = 0; i < numResidues; i++) {
            Residue residuei = allResiduesArray[i];
            atoms[i] = residuei.getReferenceAtom();
            if (residuei instanceof MultiResidue) {
                ++nMultiRes;
            }
        }
        /**
         * Use of the pre-existing ParallelTeam causes a conflict when
         * MultiResidues must re-init the force field. Temporary solution
         * for sequence optimization: if > 1 residue optimized, run on only
         * one thread.
         */
        int nThreads = 1;
        if (molecularAssembly.getPotentialEnergy().getParallelTeam() != null) {
            nThreads = (nMultiRes > 1) ? 1 : molecularAssembly.getPotentialEnergy().getParallelTeam().getThreadCount();
        } else {
            // Suggested: nThreads = (nMultiRes > 1) ? 1 : ParallelTeam.getDefaultThreadCount();
            nThreads = 16;
        }
        ParallelTeam parallelTeam = new ParallelTeam(nThreads);
        Crystal crystal = molecularAssembly.getCrystal();
        int nSymm = crystal.spaceGroup.getNumberOfSymOps();
        logger.info("\n Computing Residue Distance Matrix");
        double nlistCutoff = Math.max(Math.max(distance, twoBodyCutoffDist), threeBodyCutoffDist);
        /**
         * I think this originated from the fact that side-chain
         * (and later nucleic acid) atoms could be fairly distant
         * from the reference atom.
         */
        double magicNumberBufferOfUnknownOrigin = 25.0;
        nlistCutoff += magicNumberBufferOfUnknownOrigin;
        NeighborList neighborList = new NeighborList(null, crystal, atoms, nlistCutoff, 0.0, parallelTeam);
        // Expand coordinates
        double[][] xyz = new double[nSymm][3 * numResidues];
        double[] in = new double[3];
        double[] out = new double[3];
        for (int iSymOp = 0; iSymOp < nSymm; iSymOp++) {
            SymOp symOp = crystal.spaceGroup.getSymOp(iSymOp);
            for (int i = 0; i < numResidues; i++) {
                int i3 = i * 3;
                int iX = i3 + 0;
                int iY = i3 + 1;
                int iZ = i3 + 2;
                Atom atom = atoms[i];
                in[0] = atom.getX();
                in[1] = atom.getY();
                in[2] = atom.getZ();
                crystal.applySymOp(in, out, symOp);
                xyz[iSymOp][iX] = out[0];
                xyz[iSymOp][iY] = out[1];
                xyz[iSymOp][iZ] = out[2];
            }
        }
        // Build the residue neighbor-list.
        int[][][] lists = new int[nSymm][numResidues][];
        boolean[] use = new boolean[numResidues];
        fill(use, true);
        boolean forceRebuild = true;
        boolean printLists = false;
        long neighborTime = -System.nanoTime();
        neighborList.buildList(xyz, lists, use, forceRebuild, printLists);
        neighborTime += System.nanoTime();
        logger.info(format(" Built residue neighbor list:           %8.3f sec", neighborTime * 1.0e-9));
        DistanceRegion distanceRegion = new DistanceRegion(parallelTeam.getThreadCount(), numResidues, crystal, lists, neighborList.getPairwiseSchedule());
        long parallelTime = -System.nanoTime();
        try {
            parallelTeam.execute(distanceRegion);
        } catch (Exception e) {
            String message = " Exception compting residue distance matrix.";
            logger.log(Level.SEVERE, message, e);
        }
        parallelTime += System.nanoTime();
        logger.info(format(" Pairwise distance matrix:              %8.3f sec\n", parallelTime * 1.0e-9));
        ResidueState.revertAllCoordinates(allResiduesList, orig);
        try {
            parallelTeam.shutdown();
        } catch (Exception ex) {
            logger.warning(format(" Exception shutting down parallel team for the distance matrix: %s", ex.toString()));
        }
    }
}

Example 78 with Residue

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

Example 79 with Residue

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

the class RotamerOptimization method minMaxPairEnergy.

/**
 * Computes the maximum and minimum energy i,ri might have with j, and
 * optionally (if three-body energies in use) third residues k.
 * <p>
 * The return value should be redundant with minMax[0] being NaN.
 *
 * @param residues Array of residues under consideration.
 * @param minMax Index 0 to be filled by minimum energy, index 1 filled by
 * maximum energy.
 * @param i Some residue i under consideration.
 * @param ri A rotamer for residue i.
 * @param j Some arbitrary residue i!=j.
 * @return If a valid configuration between i,ri and j could be found.
 */
private boolean minMaxPairEnergy(Residue[] residues, double[] minMax, int i, int ri, int j) {
    Residue residuej = residues[j];
    Rotamer[] rotamersj = residuej.getRotamers(library);
    int lenrj = rotamersj.length;
    boolean valid = false;
    minMax[0] = Double.MAX_VALUE;
    minMax[1] = Double.MIN_VALUE;
    // Loop over the 2nd residues' rotamers.
    for (int rj = 0; rj < lenrj; rj++) {
        // Check for an eliminated single or eliminated pair.
        if (check(i, ri) || check(j, rj) || check(i, ri, j, rj)) {
            continue;
        }
        double currMax = get2Body(i, ri, j, rj);
        // Will remain identical if truncating at 2-body.
        double currMin = currMax;
        if (threeBodyTerm) {
            double[] minMaxTriple = new double[2];
            // Loop over residue k to find the min/max 3-Body energy.
            boolean validPair = minMax2BodySum(residues, minMaxTriple, i, ri, j, rj);
            if (!validPair) {
                // Eliminate Rotamer Pair
                Residue residuei = residues[i];
                logIfMaster(format(" Inconsistent Pair: %8s %2d, %8s %2d.", residuei.toFormattedString(false, true), ri, residuej.toFormattedString(false, true), rj), Level.INFO);
                continue;
            }
            if (Double.isFinite(currMin) && Double.isFinite(minMaxTriple[0])) {
                currMin += minMaxTriple[0];
            } else {
                currMin = Double.NaN;
            }
            if (Double.isFinite(currMax) && Double.isFinite(minMaxTriple[1])) {
                currMax += minMaxTriple[1];
            } else {
                currMax = Double.NaN;
            }
        }
        valid = true;
        if (Double.isFinite(currMin) && currMin < minMax[0]) {
            minMax[0] = currMin;
        }
        if (Double.isFinite(currMax) && Double.isFinite(minMax[1])) {
            if (currMax > minMax[1]) {
                // We have a new, finite maximum.
                minMax[1] = currMax;
            }
        // Else, if currMax is finite and less than minMax[1], we do not have a new maximum.
        } else {
            // We have a non-finite maximum.
            minMax[1] = Double.NaN;
        }
    }
    // minMax[0] being set to NaN should be redundant with valid being false.
    // It would indicate i,ri clashes with something in every possible configuration.
    minMax[0] = (minMax[0] == Double.MAX_VALUE) ? Double.NaN : minMax[0];
    // minMax[1] always gets set, unless somehow everything turns up as Double.MIN_VALUE.
    return valid;
}

Example 80 with Residue

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

the class RotamerOptimization method assignResiduesToCells.

/**
 * Constructs the cells for box optimization and assigns them residues,
 * presently based on C alpha fractional coordinates; by default, cells are
 * sorted by global index. Presently, specifying approxBoxLength over-rides
 * numXYZBoxes, and always rounds the number of boxes down (to ensure boxes
 * are always at least the specified size).
 *
 * @param crystal Crystal group.
 * @param residues List of residues to be optimized.
 * @return Array of filled Cells
 */
private void assignResiduesToCells(Crystal crystal, Residue[] residues, BoxOptCell[] cells) {
    // Search through residues, add them to all boxes containing their
    // fractional coordinates.
    int numCells = cells.length;
    int nResidues = residues.length;
    for (int i = 0; i < nResidues; i++) {
        Residue residuei = residues[i];
        double[] atomFracCoords = new double[3];
        boolean[] contained;
        double[][] originalCoordinates;
        switch(boxInclusionCriterion) {
            // As case 1 is default, test other cases first.
            case 2:
                // Residue coordinates defined by any original atomic coordinate.
                originalCoordinates = residuei.storeCoordinateArray();
                contained = new boolean[numCells];
                fill(contained, false);
                // Loop over atomic coordinates in originalCoordinates.
                for (int ai = 0; ai < originalCoordinates.length; ai++) {
                    crystal.toFractionalCoordinates(originalCoordinates[ai], atomFracCoords);
                    NeighborList.moveValuesBetweenZeroAndOne(atomFracCoords);
                    for (int j = 0; j < numCells; j++) {
                        if (!contained[j] && cells[j].checkIfContained(atomFracCoords)) {
                            cells[j].addResidue(residuei);
                            contained[j] = true;
                        }
                    }
                }
                break;
            case 3:
                // Residue coordinates defined by any atomic coordinate in any rotamer.
                // originalCoordinates = storeSingleCoordinates(residuei, true);
                ResidueState origState = residuei.storeState();
                contained = new boolean[numCells];
                fill(contained, false);
                Rotamer[] rotamersi = residuei.getRotamers(library);
                for (Rotamer rotamer : rotamersi) {
                    RotamerLibrary.applyRotamer(residuei, rotamer);
                    double[][] currentCoordinates = residuei.storeCoordinateArray();
                    for (int ai = 0; ai < currentCoordinates.length; ai++) {
                        crystal.toFractionalCoordinates(currentCoordinates[ai], atomFracCoords);
                        NeighborList.moveValuesBetweenZeroAndOne(atomFracCoords);
                        for (int j = 0; j < numCells; j++) {
                            if (!contained[j] && cells[j].checkIfContained(atomFracCoords)) {
                                cells[j].addResidue(residuei);
                                contained[j] = true;
                            }
                        }
                    }
                }
                residuei.revertState(origState);
                // revertSingleResidueCoordinates(residuei, originalCoordinates, true);
                break;
            case 1:
            default:
                // Residue coordinates defined by C alpha (protein) or N1/9
                // (nucleic acids).
                double[] cAlphaCoords = new double[3];
                residuei.getReferenceAtom().getXYZ(cAlphaCoords);
                crystal.toFractionalCoordinates(cAlphaCoords, atomFracCoords);
                NeighborList.moveValuesBetweenZeroAndOne(atomFracCoords);
                for (int j = 0; j < numCells; j++) {
                    if (cells[j].checkIfContained(atomFracCoords)) {
                        cells[j].addResidue(residuei);
                    }
                }
                break;
        }
    }
}