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

Example 36 with Rotamer

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

the class RotamerOptimization method rotamerEnergies.

/**
 * Turn off non-bonded contributions from all residues except for one.
 * Compute the self-energy for each residue relative to the backbone
 * contribution.
 *
 * @param residues A list of residues that we undergo rotamer optimization.
 * @return template energy
 */
private double rotamerEnergies(Residue[] residues) {
    if (residues == null) {
        logger.warning(" Attempt to compute rotamer energies for an empty array of residues.");
        return 0.0;
    }
    int nResidues = residues.length;
    Atom[] atoms = molecularAssembly.getAtomArray();
    int nAtoms = atoms.length;
    allocateEliminationMemory(residues);
    if (decomposeOriginal) {
        assert library.getUsingOrigCoordsRotamer();
        for (int i = 0; i < nResidues; i++) {
            Residue resi = residues[i];
            Rotamer[] rotsi = resi.getRotamers(library);
            int lenri = rotsi.length;
            // Leave the 0'th original-coordinates rotamer alone.
            for (int ri = 1; ri < lenri; ri++) {
                eliminateRotamer(residues, i, ri, false);
            }
        }
    }
    /**
     * Initialize all atoms to be used.
     */
    for (int i = 0; i < nAtoms; i++) {
        atoms[i].setUse(true);
    }
    if (!usingBoxOptimization) {
        energiesToWrite = Collections.synchronizedList(new ArrayList<String>());
        receiveThread = new ReceiveThread(residues);
        receiveThread.start();
        if (master && writeEnergyRestart && printFiles) {
            energyWriterThread = new EnergyWriterThread(receiveThread);
            energyWriterThread.start();
        }
    }
    int loaded = 0;
    if (loadEnergyRestart) {
        if (usingBoxOptimization) {
            loaded = loadEnergyRestart(energyRestartFile, residues, boxLoadIndex, boxLoadCellIndices);
        } else {
            loaded = loadEnergyRestart(energyRestartFile, residues);
        }
    }
    long energyStartTime = System.nanoTime();
    SelfEnergyRegion singlesRegion = new SelfEnergyRegion(residues);
    TwoBodyEnergyRegion pairsRegion = new TwoBodyEnergyRegion(residues);
    ThreeBodyEnergyRegion triplesRegion = new ThreeBodyEnergyRegion(residues);
    FourBodyEnergyRegion quadsRegion = new FourBodyEnergyRegion(residues);
    try {
        if (loaded < 1) {
            selfEnergyMap.clear();
            // allocate selfEnergy
            int singleJobIndex = 0;
            selfEnergy = new double[nResidues][];
            for (int i = 0; i < nResidues; i++) {
                Residue resi = residues[i];
                Rotamer[] roti = resi.getRotamers(library);
                selfEnergy[i] = new double[roti.length];
                for (int ri = 0; ri < roti.length; ri++) {
                    if (!check(i, ri)) {
                        Integer[] selfJob = { i, ri };
                        if (decomposeOriginal && ri != 0) {
                            continue;
                        }
                        selfEnergyMap.put(singleJobIndex++, selfJob);
                    }
                }
            }
        }
        // broadcast that this proc is done with startup and allocation; ready for singles
        boolean thisProcReady = true;
        BooleanBuf thisProcReadyBuf = BooleanBuf.buffer(thisProcReady);
        multicastBuf(thisProcReadyBuf);
        // launch parallel singles calculation
        while (!readyForSingles) {
            Thread.sleep(POLLING_FREQUENCY);
        }
        logger.info(String.format(" Number of self energies to calculate: %d", selfEnergyMap.size()));
        energyWorkerTeam.execute(singlesRegion);
        long singlesTime = System.nanoTime() - energyStartTime;
        logIfMaster(format(" Time for single energies: %12.4g", (singlesTime * 1.0E-9)));
        if (loaded < 2) {
            twoBodyEnergyMap.clear();
            // allocate twoBodyEnergy and create jobs
            int pairJobIndex = 0;
            twoBodyEnergy = new double[nResidues][][][];
            for (int i = 0; i < nResidues; i++) {
                Residue resi = residues[i];
                Rotamer[] roti = resi.getRotamers(library);
                twoBodyEnergy[i] = new double[roti.length][][];
                for (int ri = 0; ri < roti.length; ri++) {
                    if (check(i, ri)) {
                        continue;
                    }
                    twoBodyEnergy[i][ri] = new double[nResidues][];
                    for (int j = i + 1; j < nResidues; j++) {
                        Residue resj = residues[j];
                        Rotamer[] rotj = resj.getRotamers(library);
                        twoBodyEnergy[i][ri][j] = new double[rotj.length];
                        for (int rj = 0; rj < rotj.length; rj++) {
                            if (checkToJ(i, ri, j, rj)) {
                                continue;
                            }
                            Integer[] pairJob = { i, ri, j, rj };
                            if (decomposeOriginal && (ri != 0 || rj != 0)) {
                                continue;
                            }
                            twoBodyEnergyMap.put(pairJobIndex++, pairJob);
                        }
                    }
                }
            }
        }
        // broadcast that this proc is done with pruning and allocation; ready for pairs
        multicastBuf(thisProcReadyBuf);
        // launch parallel pair calculation
        while (!readyFor2Body) {
            Thread.sleep(POLLING_FREQUENCY);
        }
        logger.info(String.format(" Number of 2-body energies to calculate: %d", twoBodyEnergyMap.size()));
        energyWorkerTeam.execute(pairsRegion);
        long pairsTime = System.nanoTime() - (singlesTime + energyStartTime);
        long triplesTime = 0;
        long quadsTime = 0;
        logIfMaster(format(" Time for 2-body energies:   %12.4g", (pairsTime * 1.0E-9)));
        if (threeBodyTerm) {
            if (loaded < 3) {
                threeBodyEnergyMap.clear();
                // allocate threeBodyEnergy and create jobs
                int trimerJobIndex = 0;
                threeBodyEnergy = new double[nResidues][][][][][];
                for (int i = 0; i < nResidues; i++) {
                    Residue resi = residues[i];
                    Rotamer[] roti = resi.getRotamers(library);
                    threeBodyEnergy[i] = new double[roti.length][][][][];
                    for (int ri = 0; ri < roti.length; ri++) {
                        if (check(i, ri)) {
                            continue;
                        }
                        threeBodyEnergy[i][ri] = new double[nResidues][][][];
                        for (int j = i + 1; j < nResidues; j++) {
                            Residue resj = residues[j];
                            Rotamer[] rotj = resj.getRotamers(library);
                            threeBodyEnergy[i][ri][j] = new double[rotj.length][][];
                            for (int rj = 0; rj < rotj.length; rj++) {
                                if (checkToJ(i, ri, j, rj)) {
                                    continue;
                                }
                                threeBodyEnergy[i][ri][j][rj] = new double[nResidues][];
                                for (int k = j + 1; k < nResidues; k++) {
                                    Residue resk = residues[k];
                                    Rotamer[] rotk = resk.getRotamers(library);
                                    threeBodyEnergy[i][ri][j][rj][k] = new double[rotk.length];
                                    for (int rk = 0; rk < rotk.length; rk++) {
                                        if (checkToK(i, ri, j, rj, k, rk)) {
                                            continue;
                                        }
                                        Integer[] trimerJob = { i, ri, j, rj, k, rk };
                                        if (decomposeOriginal && (ri != 0 || rj != 0 || rk != 0)) {
                                            continue;
                                        }
                                        threeBodyEnergyMap.put(trimerJobIndex++, trimerJob);
                                    }
                                }
                            }
                        }
                    }
                }
            }
            // broadcast that this proc is done with pruning and allocation; ready for trimers
            multicastBuf(thisProcReadyBuf);
            // launch parallel 3-Body calculation
            while (!readyFor3Body) {
                Thread.sleep(POLLING_FREQUENCY);
            }
            logger.info(String.format(" Number of 3-Body energies to calculate: %d", threeBodyEnergyMap.size()));
            energyWorkerTeam.execute(triplesRegion);
            triplesTime = System.nanoTime() - (pairsTime + singlesTime + energyStartTime);
            logIfMaster(format(" Time for 3-Body energies: %12.4g", (triplesTime * 1.0E-9)));
        }
        if (compute4BodyEnergy) {
            logger.info(" Creating 4-Body jobs...");
            fourBodyEnergyMap.clear();
            boolean maxedOut = false;
            // create 4-Body jobs (no memory allocation)
            int quadJobIndex = 0;
            for (int i = 0; i < nResidues; i++) {
                Residue resi = residues[i];
                Rotamer[] roti = resi.getRotamers(library);
                for (int ri = 0; ri < roti.length; ri++) {
                    if (check(i, ri)) {
                        continue;
                    }
                    for (int j = i + 1; j < nResidues; j++) {
                        Residue resj = residues[j];
                        Rotamer[] rotj = resj.getRotamers(library);
                        for (int rj = 0; rj < rotj.length; rj++) {
                            /*if (check(j, rj) || check(i, ri, j, rj)) {
                                        continue;
                                    }*/
                            if (checkToJ(i, ri, j, rj)) {
                                continue;
                            }
                            for (int k = j + 1; k < nResidues; k++) {
                                Residue resk = residues[k];
                                Rotamer[] rotk = resk.getRotamers(library);
                                for (int rk = 0; rk < rotk.length; rk++) {
                                    /*if (check(k, rk) || check(i, ri, k, rk) || check(j, rj, k, rk) || check(i, ri, j, rj, k, rk)) {
                                                continue;
                                            }*/
                                    if (checkToK(i, ri, j, rj, k, rk)) {
                                        continue;
                                    }
                                    for (int l = k + 1; l < nResidues; l++) {
                                        Residue resl = residues[l];
                                        Rotamer[] rotl = resl.getRotamers(library);
                                        for (int rl = 0; rl < rotl.length; rl++) {
                                            /*if (check(l, rl) || check(i, ri, l, rl) ||
                                                            check(j, rj, l, rl) || check(k, rk, l, rl) ||
                                                            check(i, ri, j, rj, l, rl) || check(i, ri, k, rk, l, rl) ||
                                                            check(j, rj, k, rk, l, rl)) {
                                                        continue;
                                                    }*/
                                            if (checkToL(i, ri, j, rj, k, rk, l, rl)) {
                                                continue;
                                            }
                                            Integer[] quadJob = { i, ri, j, rj, k, rk, l, rl };
                                            if (decomposeOriginal && (ri != 0 || rj != 0 || rk != 0 || rl != 0)) {
                                                continue;
                                            }
                                            fourBodyEnergyMap.put(quadJobIndex++, quadJob);
                                            if (fourBodyEnergyMap.size() > max4BodyCount) {
                                                maxedOut = true;
                                                break;
                                            }
                                        }
                                        if (maxedOut) {
                                            break;
                                        }
                                    }
                                    if (maxedOut) {
                                        break;
                                    }
                                }
                                if (maxedOut) {
                                    break;
                                }
                            }
                            if (maxedOut) {
                                break;
                            }
                        }
                        if (maxedOut) {
                            break;
                        }
                    }
                    if (maxedOut) {
                        break;
                    }
                }
                if (maxedOut) {
                    break;
                }
            }
            // broadcast that this proc is done with pruning and allocation; ready for quads
            // logger.info(format(" Proc %d broadcasting ready for quads.", world.rank()));
            multicastBuf(thisProcReadyBuf);
            // launch parallel 3-Body calculation
            int waiting = 0;
            while (!readyFor4Body) {
                Thread.sleep(POLLING_FREQUENCY);
            }
            logger.info(format(" Running quads: %d jobs.", fourBodyEnergyMap.size()));
            energyWorkerTeam.execute(quadsRegion);
            quadsTime = System.nanoTime() - (triplesTime + pairsTime + singlesTime + energyStartTime);
            logIfMaster(format(" Time for 4-Body energies:   %12.4g", quadsTime * 1.0E-9));
        }
        long allTime = singlesTime + pairsTime + triplesTime + quadsTime;
        logIfMaster(format(" Time for all energies:    %12.4g", allTime * 1.0E-9));
    } catch (Exception ex) {
        String message = " Exception computing rotamer energies in parallel.";
        logger.log(Level.SEVERE, message, ex);
    }
    // Turn on all atoms.
    for (int i = 0; i < atoms.length; i++) {
        atoms[i].setUse(true);
    }
    // Print the energy with all rotamers in their default conformation.
    if (verboseEnergies && master) {
        try {
            double defaultEnergy = currentEnergy(residues);
            logger.info(format(" Energy of the system with rotamers in their default conformation: %s", formatEnergy(defaultEnergy)));
        } catch (ArithmeticException ex) {
            logger.severe(String.format(" Exception %s in calculating default energy; FFX shutting down", ex.toString()));
        }
    }
    return backboneEnergy;
}

Example 37 with Rotamer

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

the class RotamerOptimization method reconcileNARotamersWithSubsequentResidues.

/**
 * For NA residues inside some optimization window, prune any rotamers which
 * would be incompatible with the established rotamers of downstream 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 much faster. Technically, this works by pinning it to its
 * current pucker, but if the input structure is any good, it would be the
 * current pucker anyways.
 *
 * @param residues Residues to check for incompatible rotamers.
 * @return Number of rotamers eliminated for each Residue.
 */
private void reconcileNARotamersWithSubsequentResidues(Residue[] residues, int[] eliminatedRotamers) {
    for (int i = 0; i < residues.length; i++) {
        Residue residuei = residues[i];
        switch(residuei.getResidueType()) {
            case NA:
                Rotamer[] rotamers = residues[i].getRotamers(library);
                Residue nextResidue = residuei.getNextResidue();
                if (rotamers == null || nextResidue == null) {
                    break;
                }
                boolean isInList = false;
                for (Residue residue : residues) {
                    if (residue.equals(nextResidue)) {
                        isInList = true;
                        break;
                    }
                }
                if (isInList) {
                    break;
                }
                double delta = RotamerLibrary.measureDelta(residuei);
                if (RotamerLibrary.checkPucker(delta) == 1) {
                    for (int j = 0; j < rotamers.length; j++) {
                        if (RotamerLibrary.checkPucker(rotamers[j].chi7) != 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), nextResidue.toString()));
                            }
                            eliminateRotamer(residues, i, j, print);
                            eliminatedRotamers[i]++;
                        }
                    }
                } else {
                    for (int j = 0; j < rotamers.length; j++) {
                        if (RotamerLibrary.checkPucker(rotamers[j].chi7) != 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), nextResidue.toString()));
                            }
                            eliminateRotamer(residues, i, j, print);
                            eliminatedRotamers[i]++;
                        }
                    }
                }
                break;
            case AA:
            default:
                break;
        }
    }
}

Example 38 with Rotamer

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

the class RotamerOptimization method dryRun.

/**
 * A global optimization over side-chain rotamers using a recursive
 * algorithm and information about eliminated rotamers, rotamer pairs and
 * rotamer triples.
 *
 * @param residues
 * @param i
 * @param currentRotamers
 * @return 0.
 */
private double dryRun(Residue[] residues, int i, int[] currentRotamers) {
    // This is the initialization condition.
    if (i == 0) {
        evaluatedPermutations = 0;
        evaluatedPermutationsPrint = 1000;
    }
    if (evaluatedPermutations >= evaluatedPermutationsPrint) {
        if (evaluatedPermutations % evaluatedPermutationsPrint == 0) {
            logIfMaster(format(" The permutations have reached %10.4e.", (double) evaluatedPermutationsPrint));
            evaluatedPermutationsPrint *= 10;
        }
    }
    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;
            dryRun(residues, i + 1, currentRotamers);
        }
    } 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) {
                evaluatedPermutations++;
            }
        }
    }
    return 0.0;
}

Example 39 with Rotamer

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

the class RotamerOptimization method validateDEE.

private boolean validateDEE(Residue[] residues) {
    int nres = eliminatedSingles.length;
    // Validate residues
    for (int i = 0; i < nres; i++) {
        Residue residuei = residues[i];
        int ni = eliminatedSingles[i].length;
        boolean valid = false;
        for (int ri = 0; ri < ni; ri++) {
            if (!check(i, ri)) {
                valid = true;
            }
        }
        if (!valid) {
            logger.severe(format(" Coding error: all %d rotamers for residue %s eliminated.", ni, residuei));
        }
    }
    // Validate pairs
    for (int i = 0; i < nres; i++) {
        Residue residuei = residues[i];
        Rotamer[] rotamersi = residuei.getRotamers(library);
        int ni = rotamersi.length;
        for (int j = i + 1; j < nres; j++) {
            Residue residuej = residues[j];
            Rotamer[] rotamersj = residuej.getRotamers(library);
            int nj = rotamersj.length;
            boolean valid = false;
            for (int ri = 0; ri < ni; ri++) {
                for (int rj = 0; rj < nj; rj++) {
                    if (!check(i, ri, j, rj)) {
                        valid = true;
                    }
                }
            }
            if (!valid) {
                logger.severe(format(" Coding error: all pairs for %s with residue %s eliminated.", residuei.toFormattedString(false, true), residuej));
            }
        }
    }
    return true;
}

Example 40 with Rotamer

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

the class RotamerOptimization method minMax2BodySum.

/**
 * Find the min/max of the 2-body energy.
 *
 * @param residues The residue array.
 * @param minMax The bound on the 3-body energy (minMax[0] = min, minMax[1]
 * = max.
 * @param i Residue i
 * @param ri Rotamer ri of Residue i
 * @param j Residue j
 * @param rj Rotamer rj of Residue j
 * @return true if this term is valid.
 */
private boolean minMax2BodySum(Residue[] residues, double[] minMax, int i, int ri, int j, int rj) {
    int nres = residues.length;
    double minSum = 0.0;
    double maxSum = 0.0;
    for (int k = 0; k < nres; k++) {
        if (k == i || k == j) {
            continue;
        }
        Residue residuek = residues[k];
        Rotamer[] romatersk = residuek.getRotamers(library);
        int lenrk = romatersk.length;
        double currentMin = Double.MAX_VALUE;
        double currentMax = Double.MIN_VALUE;
        for (int rk = 0; rk < lenrk; rk++) {
            if (check(k, rk)) {
                // k,rk is part of no valid phase space, so ignore it.
                continue;
            }
            if (check(i, ri, k, rk) || check(j, rj, k, rk)) {
                // Not implemented: check(i, ri, j, rj, k, rk).
                // k,rk conflicts with i,ri or j,rj, so the max is now Double.NaN. No effect on minimum.
                currentMax = Double.NaN;
            } else {
                double current = get3Body(i, ri, j, rj, k, rk);
                if (Double.isFinite(current) && current < currentMin) {
                    currentMin = current;
                }
                // Else, no new minimum found.
                if (Double.isFinite(current) && Double.isFinite(currentMax)) {
                    if (current > currentMax) {
                        currentMax = current;
                    }
                // Else, we have failed to find a new finite maximum.
                } else {
                    // The maximum is NaN.
                    currentMax = Double.NaN;
                }
            }
        }
        if (currentMin == Double.MAX_VALUE || !Double.isFinite(minSum)) {
            // We have failed to find a viable configuration for i,ri,j,rj, as it conflicts with all rk for this k.
            minMax[0] = Double.NaN;
            minMax[1] = Double.NaN;
            return false;
        } else {
            // Add finite current min to minSum.
            minSum += currentMin;
        }
        if (Double.isFinite(maxSum) && Double.isFinite(currentMax)) {
            maxSum += currentMax;
        } else {
            maxSum = Double.NaN;
        }
    }
    minMax[0] = minSum;
    minMax[1] = maxSum;
    return true;
}