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

Example 61 with Residue

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

the class RotamerOptimization method loadEnergyRestart.

private int loadEnergyRestart(File restartFile, Residue[] residues, int boxIteration, int[] cellIndices) {
    try {
        int nResidues = residues.length;
        Path path = Paths.get(restartFile.getCanonicalPath());
        List<String> lines = Files.readAllLines(path, StandardCharsets.UTF_8);
        List<String> linesThisBox = new ArrayList<>();
        try {
            backboneEnergy = computeBackboneEnergy(residues);
        } catch (ArithmeticException ex) {
            logger.severe(String.format(" Exception %s in calculating backbone energy; FFX shutting down.", ex.toString()));
        }
        if (usingBoxOptimization && boxIteration >= 0) {
            boolean foundBox = false;
            for (int i = 0; i < lines.size(); i++) {
                String line = lines.get(i);
                if (line.startsWith("Box")) {
                    String[] tok = line.replaceAll("Box", "").replaceAll(":", ",").replaceAll(" ", "").split(",");
                    int readIteration = Integer.parseInt(tok[0]);
                    int readCellIndexX = Integer.parseInt(tok[1]);
                    int readCellIndexY = Integer.parseInt(tok[2]);
                    int readCellIndexZ = Integer.parseInt(tok[3]);
                    if (readIteration == boxIteration && readCellIndexX == cellIndices[0] && readCellIndexY == cellIndices[1] && readCellIndexZ == cellIndices[2]) {
                        foundBox = true;
                        for (int j = i + 1; j < lines.size(); j++) {
                            String l = lines.get(j);
                            if (l.startsWith("Box")) {
                                break;
                            }
                            linesThisBox.add(l);
                        }
                        break;
                    }
                }
            }
            if (!foundBox) {
                logIfMaster(format(" Didn't find restart energies for Box %d: %d,%d,%d", boxIteration, cellIndices[0], cellIndices[1], cellIndices[2]));
                return 0;
            } else if (linesThisBox.size() == 0) {
                return 0;
            } else {
                lines = linesThisBox;
            }
        }
        List<String> singleLines = new ArrayList<>();
        List<String> pairLines = new ArrayList<>();
        List<String> tripleLines = new ArrayList<>();
        for (String line : lines) {
            String[] tok = line.split("\\s");
            if (tok[0].startsWith("Self")) {
                singleLines.add(line);
            } else if (tok[0].startsWith("Pair")) {
                pairLines.add(line);
            } else if (tok[0].startsWith("Triple")) {
                tripleLines.add(line);
            }
        }
        int loaded = 0;
        if (tripleLines.size() > 0) {
            loaded = 3;
        } else if (pairLines.size() > 0) {
            loaded = 2;
        } else if (singleLines.size() > 0) {
            loaded = 1;
        } else {
            logger.warning(format(" Empty or unreadable energy restart file: %s.", restartFile.getCanonicalPath()));
        }
        if (loaded >= 1) {
            selfEnergyMap.clear();
            // allocate selfEnergy array and create self jobs
            HashMap<String, Integer> reverseJobMapSingles = new HashMap<>();
            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++) {
                    Integer[] selfJob = { i, ri };
                    if (decomposeOriginal && ri != 0) {
                        continue;
                    }
                    selfEnergyMap.put(singleJobIndex, selfJob);
                    String revKey = format("%d %d", i, ri);
                    reverseJobMapSingles.put(revKey, singleJobIndex);
                    singleJobIndex++;
                }
            }
            // fill in self-energies from file while removing the corresponding jobs from selfEnergyMap
            for (String line : singleLines) {
                try {
                    String[] tok = line.replace(",", "").replace(":", "").split("\\s+");
                    int i;
                    if (tok[1].contains("-")) {
                        i = nameToNumber(tok[1], residues);
                    } else {
                        i = Integer.parseInt(tok[1]);
                    }
                    int ri = Integer.parseInt(tok[2]);
                    double energy = Double.parseDouble(tok[3]);
                    try {
                        selfEnergy[i][ri] = energy;
                        if (verbose) {
                            logIfMaster(format(" From restart file: Self energy %3d (%8s,%2d): %s", i, residues[i].toFormattedString(false, true), ri, formatEnergy(energy)));
                        }
                    } catch (Exception e) {
                        if (verbose) {
                            logIfMaster(format(" Restart file out-of-bounds index: %s", line));
                        }
                    }
                    // remove that job from the pool
                    String revKey = format("%d %d", i, ri);
                    Integer[] ret = selfEnergyMap.remove(reverseJobMapSingles.get(revKey));
                    if (ret == null) {
                    // logIfMaster(format("(sdl %d) Restart file contained unnecessary value for %s", BOXNUM, revKey));
                    }
                } catch (NumberFormatException ex) {
                    logger.log(Level.WARNING, format(" Unparsable line in energy restart file: \n%s", line), ex);
                }
            }
            logIfMaster(" Loaded self energies from restart file.");
            // Pre-Prune if self-energy is Double.NaN.
            for (int i = 0; i < residues.length; i++) {
                Residue residue = residues[i];
                Rotamer[] rotamers = residue.getRotamers(library);
                int nrot = rotamers.length;
                for (int ri = 0; ri < nrot; ri++) {
                    if (!check(i, ri) && Double.isNaN(getSelf(i, ri))) {
                        logIfMaster(format(" Rotamer (%7s,%2d) self-energy %12.4f pre-pruned since energy is NaN.", residue, ri, getSelf(i, ri)));
                        eliminateRotamer(residues, i, ri, false);
                    }
                }
            }
            // prune singles
            if (pruneClashes) {
                pruneSingleClashes(residues);
            }
        }
        /**
         * Remap to sequential integer keys.
         */
        condenseEnergyMap(selfEnergyMap);
        if (loaded >= 2) {
            if (selfEnergyMap.size() > 0) {
                logIfMaster(" Double-check that parameters match original run due to missing self-energies.");
            }
            twoBodyEnergyMap.clear();
            // allocated twoBodyEnergy array and create pair jobs
            HashMap<String, Integer> reverseJobMapPairs = new HashMap<>();
            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 (check(j, rj) || check(i, ri, j, rj)) {
                                continue;
                            }
                            Integer[] pairJob = { i, ri, j, rj };
                            if (decomposeOriginal && (ri != 0 || rj != 0)) {
                                continue;
                            }
                            twoBodyEnergyMap.put(pairJobIndex, pairJob);
                            String revKey = format("%d %d %d %d", i, ri, j, rj);
                            reverseJobMapPairs.put(revKey, pairJobIndex);
                            pairJobIndex++;
                        }
                    }
                }
            }
            // fill in pair-energies from file while removing the corresponding jobs from twoBodyEnergyMap
            for (String line : pairLines) {
                try {
                    String[] tok = line.replace(",", "").replace(":", "").split("\\s+");
                    int i;
                    if (tok[1].contains("-")) {
                        i = nameToNumber(tok[1], residues);
                    } else {
                        i = Integer.parseInt(tok[1]);
                    }
                    int ri = Integer.parseInt(tok[2]);
                    int j;
                    if (tok[3].contains("-")) {
                        j = nameToNumber(tok[3], residues);
                    } else {
                        j = Integer.parseInt(tok[3]);
                    }
                    int rj = Integer.parseInt(tok[4]);
                    double energy = Double.parseDouble(tok[5]);
                    try {
                        twoBodyEnergy[i][ri][j][rj] = energy;
                        if (verbose) {
                            logIfMaster(format(" From restart file: Pair energy [(%8s,%2d),(%8s,%2d)]: %12.4f", residues[i].toFormattedString(false, true), ri, residues[j].toFormattedString(false, true), rj, energy));
                        }
                    } catch (Exception e) {
                        if (verbose) {
                            logIfMaster(format(" Restart file out-of-bounds index: %s", line));
                        }
                    }
                    // remove that job from the pool
                    String revKey = format("%d %d %d %d", i, ri, j, rj);
                    Integer[] ret = twoBodyEnergyMap.remove(reverseJobMapPairs.get(revKey));
                } catch (NumberFormatException ex) {
                    logger.log(Level.WARNING, format("Unparsable line in energy restart file: \n%s", line), ex);
                }
            }
            logIfMaster(" Loaded 2-body energies from restart file.");
            // Loop over first residue.
            for (int i = 0; i < nResidues - 1; i++) {
                Residue resi = residues[i];
                Rotamer[] roti = resi.getRotamers(library);
                int ni = roti.length;
                // Loop over second residue.
                for (int j = i + 1; j < nResidues; j++) {
                    Residue resj = residues[j];
                    Rotamer[] rotj = resj.getRotamers(library);
                    int nj = rotj.length;
                    // Loop over the rotamers for residue i.
                    for (int ri = 0; ri < ni; ri++) {
                        if (!validRotamer(residues, i, ri)) {
                            continue;
                        }
                        // Loop over rotamers for residue j.
                        for (int rj = 0; rj < nj; rj++) {
                            if (!validRotamer(residues, j, rj) || check(i, ri, j, rj)) {
                                continue;
                            }
                            if (!check(i, ri, j, rj) && Double.isNaN(get2Body(i, ri, j, rj))) {
                                logIfMaster(format(" Rotamer Pair (%7s,%2d) (%7s,%2d) 2-body energy %12.4f pre-pruned since energy is NaN.", i, ri, j, rj, get2Body(i, ri, j, rj)));
                                eliminateRotamerPair(residues, i, ri, j, rj, print);
                            }
                        }
                    }
                }
            }
            // prune pairs
            if (prunePairClashes) {
                prunePairClashes(residues);
            }
        }
        /**
         * Remap to sequential integer keys.
         */
        condenseEnergyMap(twoBodyEnergyMap);
        if (loaded >= 3) {
            if (twoBodyEnergyMap.size() > 0) {
                if (master) {
                    logger.warning("Double-check that parameters match original run!  Found trimers in restart file, but pairs job queue is non-empty.");
                }
            }
            HashMap<String, Integer> reverseJobMapTrimers = new HashMap<>();
            threeBodyEnergyMap.clear();
            // allocate 3-Body energy array, fill in 3-Body energies from the restart file.
            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 (check(j, rj) || check(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 (check(k, rk) || check(i, ri, k, rk) || check(j, rj, k, rk)) {
                                        // Not implemented: check(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);
                                    String revKey = format("%d %d %d %d %d %d", i, ri, j, rj, k, rk);
                                    reverseJobMapTrimers.put(revKey, trimerJobIndex);
                                    trimerJobIndex++;
                                }
                            }
                        }
                    }
                }
            }
            // fill in 3-Body energies from file while removing the corresponding jobs from threeBodyEnergyMap
            for (String line : tripleLines) {
                try {
                    String[] tok = line.replace(",", "").replace(":", "").split("\\s+");
                    int i;
                    if (tok[1].contains("-")) {
                        i = nameToNumber(tok[1], residues);
                    } else {
                        i = Integer.parseInt(tok[1]);
                    }
                    int ri = Integer.parseInt(tok[2]);
                    int j;
                    if (tok[3].contains("-")) {
                        j = nameToNumber(tok[3], residues);
                    } else {
                        j = Integer.parseInt(tok[3]);
                    }
                    int rj = Integer.parseInt(tok[4]);
                    int k;
                    if (tok[5].contains("-")) {
                        k = nameToNumber(tok[5], residues);
                    } else {
                        k = Integer.parseInt(tok[5]);
                    }
                    int rk = Integer.parseInt(tok[6]);
                    double energy = Double.parseDouble(tok[7]);
                    try {
                        threeBodyEnergy[i][ri][j][rj][k][rk] = energy;
                    } catch (ArrayIndexOutOfBoundsException ex) {
                        if (verbose) {
                            logIfMaster(format(" Restart file out-of-bounds index: %s", line));
                        }
                    } catch (NullPointerException npe) {
                        if (verbose) {
                            logIfMaster(format(" NPE in loading 3-body energies: pruning " + "likely changed! 3-body %s-%d %s-%d %s-%d", residues[i].toFormattedString(false, true), ri, residues[j], rj, residues[k], rk));
                        }
                    }
                    if (verbose) {
                        logIfMaster(format(" From restart file: Trimer energy %3d %-2d, %3d %-2d, %3d %-2d: %s", i, ri, j, rj, k, rk, formatEnergy(energy)));
                    }
                    // remove that job from the pool
                    String revKey = format("%d %d %d %d %d %d", i, ri, j, rj, k, rk);
                    Integer[] ret = threeBodyEnergyMap.remove(reverseJobMapTrimers.get(revKey));
                    if (ret == null) {
                    // logIfMaster(format("(sdl %d) Restart file contained unnecessary value for %s", BOXNUM, revKey));
                    }
                } catch (NumberFormatException ex) {
                    logger.log(Level.WARNING, format("Unparsable line in energy restart file: \n%s", line), ex);
                }
            }
            logIfMaster(" Loaded trimer energies from restart file.");
        }
        /**
         * Remap to sequential integer keys.
         */
        condenseEnergyMap(threeBodyEnergyMap);
        return loaded;
    } catch (IOException ex) {
        logger.log(Level.WARNING, "Exception while loading energy restart file.", ex);
    }
    return 0;
}

Example 62 with Residue

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

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

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

the class RotamerOptimization method boxOptimization.

/**
 * Breaks down a structure into a number of overlapping boxes for
 * optimization.
 *
 * @return Potential energy of final structure.
 */
private double boxOptimization(ArrayList<Residue> residueList) {
    this.usingBoxOptimization = true;
    long beginTime = -System.nanoTime();
    Residue[] residues = residueList.toArray(new Residue[residueList.size()]);
    boolean firstCellSaved = false;
    /*
         * A new dummy Crystal will be constructed for an aperiodic system. The
         * purpose is to avoid using the overly large dummy Crystal used for
         * Ewald purposes. Atoms are not and should not be moved into the dummy
         * Crystal boundaries; to check if an Atom is inside a cell, use an
         * array of coordinates adjusted to be 0 < coordinate < 1.0.
         */
    Crystal crystal = generateSuperbox(residueList);
    // Cells indexed by x*(YZ divisions) + y*(Z divisions) + z.
    // Also initializes cell count if using -bB
    int totalCells = getTotalCellCount(crystal);
    if (boxStart > totalCells - 1) {
        logger.severe(format(" FFX shutting down: Box optimization start is out of range of total boxes: %d > %d", (boxStart + 1), totalCells));
    }
    if (boxEnd > totalCells - 1) {
        boxEnd = totalCells - 1;
        logIfMaster(" Final box out of range: reset to last possible box.");
    } else if (boxEnd < 0) {
        boxEnd = totalCells - 1;
    }
    BoxOptCell[] cells = loadCells(crystal, residues);
    int numCells = cells.length;
    logIfMaster(format(" Optimizing boxes  %d  to  %d", (boxStart + 1), (boxEnd + 1)));
    for (int i = 0; i < numCells; i++) {
        BoxOptCell celli = cells[i];
        ArrayList<Residue> residuesList = celli.getResiduesAsList();
        int[] cellIndices = celli.getXYZIndex();
        logIfMaster(format("\n Iteration %d of the box optimization.", (i + 1)));
        logIfMaster(format(" Cell index (linear): %d", (celli.getLinearIndex() + 1)));
        logIfMaster(format(" Cell xyz indices: x = %d, y = %d, z = %d", cellIndices[0] + 1, cellIndices[1] + 1, cellIndices[2] + 1));
        int nResidues = residuesList.size();
        if (nResidues > 0) {
            readyForSingles = false;
            finishedSelf = false;
            readyFor2Body = false;
            finished2Body = false;
            readyFor3Body = false;
            finished3Body = false;
            energiesToWrite = Collections.synchronizedList(new ArrayList<String>());
            receiveThread = new ReceiveThread(residuesList.toArray(new Residue[1]));
            receiveThread.start();
            if (master && writeEnergyRestart && printFiles) {
                if (energyWriterThread != null) {
                    int waiting = 0;
                    while (energyWriterThread.getState() != java.lang.Thread.State.TERMINATED) {
                        try {
                            if (waiting++ > 20) {
                                logger.log(Level.SEVERE, " ReceiveThread/EnergyWriterThread from previous box locked up.");
                            }
                            logIfMaster(" Waiting for previous iteration's communication threads to shut down... ");
                            Thread.sleep(10000);
                        } catch (InterruptedException ex) {
                        }
                    }
                }
                energyWriterThread = new EnergyWriterThread(receiveThread, i + 1, cellIndices);
                energyWriterThread.start();
            }
            if (loadEnergyRestart) {
                boxLoadIndex = i + 1;
                boxLoadCellIndices = new int[3];
                boxLoadCellIndices[0] = cellIndices[0];
                boxLoadCellIndices[1] = cellIndices[1];
                boxLoadCellIndices[2] = cellIndices[2];
            }
            long boxTime = -System.nanoTime();
            Residue firstResidue = residuesList.get(0);
            Residue lastResidue = residuesList.get(nResidues - 1);
            if (firstResidue != lastResidue) {
                logIfMaster(format(" Residues %s ... %s", firstResidue.toString(), lastResidue.toString()));
            } else {
                logIfMaster(format(" Residue %s", firstResidue.toString()));
            }
            if (revert) {
                ResidueState[] coordinates = ResidueState.storeAllCoordinates(residuesList);
                // If x has not yet been constructed, construct it.
                if (x == null) {
                    Atom[] atoms = molecularAssembly.getAtomArray();
                    int nAtoms = atoms.length;
                    x = new double[nAtoms * 3];
                }
                double startingEnergy = 0;
                double finalEnergy = 0;
                try {
                    startingEnergy = currentEnergy(residuesList);
                } catch (ArithmeticException ex) {
                    logger.severe(String.format(" Exception %s in calculating starting energy of a box; FFX shutting down", ex.toString()));
                }
                globalOptimization(residuesList);
                try {
                    finalEnergy = currentEnergy(residuesList);
                } catch (ArithmeticException ex) {
                    logger.severe(String.format(" Exception %s in calculating starting energy of a box; FFX shutting down", ex.toString()));
                }
                if (startingEnergy <= finalEnergy) {
                    logger.warning("Optimization did not yield a better energy. Reverting to orginal coordinates.");
                    ResidueState.revertAllCoordinates(residuesList, coordinates);
                }
                long currentTime = System.nanoTime();
                boxTime += currentTime;
                logIfMaster(format(" Time elapsed for this iteration: %11.3f sec", boxTime * 1.0E-9));
                logIfMaster(format(" Overall time elapsed: %11.3f sec", (currentTime + beginTime) * 1.0E-9));
            } else {
                globalOptimization(residuesList);
                long currentTime = System.nanoTime();
                boxTime += currentTime;
                logIfMaster(format(" Time elapsed for this iteration: %11.3f sec", boxTime * 1.0E-9));
                logIfMaster(format(" Overall time elapsed: %11.3f sec", (currentTime + beginTime) * 1.0E-9));
            }
            if (master && printFiles) {
                String filename = FilenameUtils.removeExtension(molecularAssembly.getFile().getAbsolutePath()) + ".partial";
                File file = new File(filename);
                if (firstCellSaved) {
                    file.delete();
                }
                // Don't write a file if it's the final iteration.
                if (i == (numCells - 1)) {
                    continue;
                }
                PDBFilter windowFilter = new PDBFilter(file, molecularAssembly, null, null);
                try {
                    windowFilter.writeFile(file, false);
                    if (firstResidue != lastResidue) {
                        logIfMaster(format(" File with residues %s ... %s in window written.", firstResidue.toString(), lastResidue.toString()));
                    } else {
                        logIfMaster(format(" File with residue %s in window written.", firstResidue.toString()));
                    }
                    firstCellSaved = true;
                } catch (Exception e) {
                    logger.warning(format("Exception writing to file: %s", file.getName()));
                }
            }
        /*for (Residue residue : residueList) {
                    if (residue instanceof MultiResidue) {
                        ((MultiResidue) residue).setDefaultResidue();
                        residue.reInitOriginalAtomList();
                    }
                }*/
        } else {
            logIfMaster(format(" Empty box: no residues found."));
        }
    }
    return 0.0;
}

Example 65 with Residue

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