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

Example 6 with ResidueState

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

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

the class RotamerOptimization method globalOptimization.

/**
 * The main driver for optimizing a block of residues using DEE.
 *
 * @param residueList Residues to optimize.
 * @return Final energy.
 */
private double globalOptimization(List<Residue> residueList) {
    int currentEnsemble = Integer.MAX_VALUE;
    Residue[] residues = residueList.toArray(new Residue[residueList.size()]);
    int nResidues = residues.length;
    int[] currentRotamers = new int[nResidues];
    int iterations = 0;
    boolean finalTry = false;
    int bestEnsembleTargetDiffThusFar = Integer.MAX_VALUE;
    double bestBufferThusFar = ensembleBuffer;
    double startingBuffer = ensembleBuffer;
    optimum = new int[nResidues];
    if (ensembleEnergy > 0.0) {
        ensembleBuffer = ensembleEnergy;
        applyEliminationCriteria(residues, true, true);
        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
         * and compute the number of permutations using singleton
         * elimination.
         */
        double permutations = 1;
        double singletonPermutations = 1;
        for (int i = 0; i < nResidues; i++) {
            Residue residue = residues[i];
            Rotamer[] rotamers = residue.getRotamers(library);
            int nr = rotamers.length;
            if (nr > 1) {
                int nrot = 0;
                for (int ri = 0; ri < nr; ri++) {
                    if (!eliminatedSingles[i][ri]) {
                        nrot++;
                    }
                }
                permutations *= rotamers.length;
                if (nrot > 1) {
                    singletonPermutations *= nrot;
                }
            }
        }
        dryRun(residues, 0, currentRotamers);
        double pairTotalElimination = singletonPermutations - (double) evaluatedPermutations;
        double afterPairElim = singletonPermutations - pairTotalElimination;
        if (evaluatedPermutations == 0) {
            logger.severe(" No valid path through rotamer space found; try recomputing without pruning or using ensemble.");
        }
        if (master && printFiles && ensembleFile == null) {
            File file = molecularAssembly.getFile();
            String filename = FilenameUtils.removeExtension(file.getAbsolutePath());
            ensembleFile = new File(filename + ".ens");
            if (ensembleFile.exists()) {
                for (int i = 2; i < 1000; i++) {
                    ensembleFile = new File(filename + ".ens_" + i);
                    if (!ensembleFile.exists()) {
                        break;
                    }
                }
                if (ensembleFile.exists()) {
                    logger.warning(format(" Versioning failed: appending to end of file %s", ensembleFile.getName()));
                }
            }
            ensembleFilter = new PDBFilter(new File(ensembleFile.getName()), molecularAssembly, null, null);
            logger.info(format(" Ensemble file: %s", ensembleFile.getName()));
        }
        logIfMaster(format("%30s %5s %30s %5s %30s %5s", "Condition", "|", "Number of Permutations Left", "|", "Removed", "|"));
        logIfMaster(format("%s", " -------------------------------------------------------------------------------------------------------------"));
        logIfMaster(format("%30s %5s %30s %5s %30s %5s", "No Eliminations", "|", permutations, "|", "", "|"));
        logIfMaster(format("%30s %5s %30s %5s %30s %5s", "Single Eliminations", "|", singletonPermutations, "|", permutations - singletonPermutations, "|"));
        logIfMaster(format("%30s %5s %30s %5s %30s %5s", "Pair Eliminations", "|", afterPairElim, "|", pairTotalElimination, "|"));
        logIfMaster(format("%30s %5s %30s %5s %30s %5s", "Single and Pair Eliminations", "|", (double) evaluatedPermutations, "|", pairTotalElimination + (permutations - singletonPermutations), "|"));
        logIfMaster(format("%s", " -------------------------------------------------------------------------------------------------------------\n"));
        logIfMaster(format(" Energy of permutations:"));
        logIfMaster(format("%s", " ----------------------------------------------------------------------------------"));
        logIfMaster(format(" %12s %5s %25s %5s %25s %5s", "Permutation", "|", "Energy", "|", "Lowest Possible Energy", "|"));
        logIfMaster(format("%s", " ----------------------------------------------------------------------------------"));
        double e;
        if (useMonteCarlo()) {
            firstValidPerm(residues, 0, currentRotamers);
            System.arraycopy(currentRotamers, 0, optimum, 0, nResidues);
            rotamerOptimizationMC(residues, optimum, currentRotamers, nMCsteps, false, mcUseAll);
            logIfMaster(" Ensembles not currently compatible with Monte Carlo search");
        /**
         * Not currently compatible with ensembles.
         */
        } else {
            double[] permutationEnergies = new double[evaluatedPermutations];
            ensembleStates = new ArrayList<>();
            e = rotamerOptimizationDEE(molecularAssembly, residues, 0, currentRotamers, Double.MAX_VALUE, optimum, permutationEnergies);
            int[][] acceptedPermutations = new int[evaluatedPermutations][];
            for (int i = 0; i < acceptedPermutations.length; i++) {
                acceptedPermutations[i] = null;
            }
            logIfMaster(format("\n Checking permutations for distance < %5.3f kcal/mol from GMEC energy %10.8f kcal/mol", ensembleEnergy, e));
            dryRunForEnsemble(residues, 0, currentRotamers, e, permutationEnergies, acceptedPermutations);
            int numAcceptedPermutations = 0;
            for (int i = 0; i < acceptedPermutations.length; i++) {
                if (acceptedPermutations[i] != null) {
                    ++numAcceptedPermutations;
                    logIfMaster(format(" Accepting permutation %d at %8.6f < %8.6f", i, permutationEnergies[i] - e, ensembleEnergy));
                    for (int j = 0; j < nResidues; j++) {
                        Residue residuej = residues[j];
                        Rotamer[] rotamersj = residuej.getRotamers(library);
                        RotamerLibrary.applyRotamer(residuej, rotamersj[acceptedPermutations[i][j]]);
                    }
                    ResidueState[] states = ResidueState.storeAllCoordinates(residues);
                    ensembleStates.add(new ObjectPair<>(states, permutationEnergies[i]));
                    if (printFiles && master) {
                        try {
                            FileWriter fw = new FileWriter(ensembleFile, true);
                            BufferedWriter bw = new BufferedWriter(fw);
                            bw.write(format("MODEL        %d", numAcceptedPermutations));
                            for (int j = 0; j < 75; j++) {
                                bw.write(" ");
                            }
                            bw.newLine();
                            bw.flush();
                            ensembleFilter.writeFile(ensembleFile, true);
                            bw.write(format("ENDMDL"));
                            for (int j = 0; j < 64; j++) {
                                bw.write(" ");
                            }
                            bw.newLine();
                            bw.close();
                        } catch (IOException ex) {
                            logger.warning(format(" Exception writing to file: %s", ensembleFile.getName()));
                        }
                    }
                }
            }
            logIfMaster(format(" Number of permutations within %5.3f kcal/mol of GMEC energy: %6.4e", ensembleEnergy, (double) numAcceptedPermutations));
            ensembleStates.sort(null);
        }
        logIfMaster(format(" Final rotamers:"));
        logIfMaster(format("%s", " --------------------------------------------------------------------------------------------"));
        logIfMaster(format("%14s %3s %10s %3s %9s %3s %9s %3s %9s %3s", "Residue", "|", "Chi 1", "|", "Chi 2", "|", "Chi 3", "|", "Chi 4", "|"));
        logIfMaster(format("%s", " --------------------------------------------------------------------------------------------"));
        for (int i = 0; i < nResidues; i++) {
            Residue residue = residues[i];
            Rotamer[] rotamers = residue.getRotamers(library);
            int ri = optimum[i];
            Rotamer rotamer = rotamers[ri];
            logIfMaster(format(" %c (%7s,2d) | %s", residue.getChainID(), residue, ri, rotamer.toAngleString()));
            RotamerLibrary.applyRotamer(residue, rotamer);
        }
        logIfMaster(format("%s", " --------------------------------------------------------------------------------------------\n"));
        double sumSelfEnergy = 0;
        double sumPairEnergy = 0;
        double sumTrimerEnergy = 0;
        for (int i = 0; i < nResidues; i++) {
            int ri = optimum[i];
            sumSelfEnergy += getSelf(i, ri);
            logIfMaster(format(" Final self Energy (%8s,%2d): %12.4f", residues[i].toFormattedString(false, true), ri, getSelf(i, ri)));
        }
        for (int i = 0; i < nResidues - 1; i++) {
            int ri = optimum[i];
            for (int j = i + 1; j < nResidues; j++) {
                int rj = optimum[j];
                sumPairEnergy += get2Body(i, ri, j, rj);
                if (get2Body(i, ri, j, rj) > 10.0) {
                    logIfMaster(format(" Large Final Pair Energy (%8s,%2d) (%8s,%2d): %12.4f", residues[i].toFormattedString(false, true), ri, residues[j].toFormattedString(false, true), rj, get2Body(i, ri, j, rj)));
                }
            }
        }
        try {
            e = currentEnergy(residueList);
        } catch (ArithmeticException ex) {
            e = Double.NaN;
            logger.severe(String.format(" Exception %s in calculating current energy at the end of triples", ex.toString()));
        }
        logIfMaster(format(" %12s %5s %25s %5s %25s %5s", "Type", "|", "Energy", "|", "Lowest Possible Energy", "|"));
        logIfMaster(format("%s", " ----------------------------------------------------------------------------------"));
        logIfMaster(format(" %12s %5s %25f %5s %25s %5s", "Self:", "|", sumSelfEnergy, "|", "", "|"));
        logIfMaster(format(" %12s %5s %25f %5s %25s %5s", "Pair:", "|", sumPairEnergy, "|", "", "|"));
        double approximateEnergy = backboneEnergy + sumSelfEnergy + sumPairEnergy;
        if (threeBodyTerm) {
            for (int i = 0; i < nResidues - 2; i++) {
                int ri = optimum[i];
                for (int j = i + 1; j < nResidues - 1; j++) {
                    int rj = optimum[j];
                    for (int k = j + 1; k < nResidues; k++) {
                        int rk = optimum[k];
                        try {
                            sumTrimerEnergy += get3Body(i, ri, j, rj, k, rk);
                        } catch (Exception ex) {
                            logger.warning(ex.toString());
                        }
                    }
                }
            }
            approximateEnergy += sumTrimerEnergy;
            double higherOrderEnergy = e - sumSelfEnergy - sumPairEnergy - sumTrimerEnergy - backboneEnergy;
            logIfMaster(format(" %12s %5s %25f %5s %25s %5s", "Trimer:", "|", sumTrimerEnergy, "|", "", "|"));
            logIfMaster(format(" %12s %5s %25f %5s %25s %5s", "Neglected:", "|", higherOrderEnergy, "|", "", "|"));
        } else {
            double higherOrderEnergy = e - sumSelfEnergy - sumPairEnergy - backboneEnergy;
            logIfMaster(format(" %12s %5s %25f %5s %25s %5s", "Neglected:", "|", higherOrderEnergy, "|", "", "|"));
        }
        logIfMaster(format(" %12s %5s %25f %5s %25s %5s", "Approximate:", "|", approximateEnergy, "|", "", "|"));
        logIfMaster(format("%s", " ----------------------------------------------------------------------------------\n"));
        return e;
    }
    /**
     * Permutations used only to set maximum bound on ensembleNumber, thus
     * it is safe here to put that value in a 32-bit int.
     */
    int nPerms = 1;
    for (int i = 0; i < nResidues; i++) {
        Residue residue = residues[i];
        Rotamer[] rotamers = residue.getRotamers(library);
        int nr = rotamers.length;
        if (nr > 1) {
            nPerms *= rotamers.length;
        }
        if (nPerms > ensembleNumber) {
            break;
        }
    }
    if (nPerms < ensembleNumber) {
        logger.info(format(" Requested an ensemble of %d, but only %d permutations exist; returning full ensemble", ensembleNumber, nPerms));
        ensembleNumber = nPerms;
    }
    while (currentEnsemble != ensembleNumber) {
        if (monteCarlo) {
            logIfMaster(" Ensemble search not currently compatible with Monte Carlo");
            ensembleNumber = 1;
        }
        if (iterations == 0) {
            applyEliminationCriteria(residues, true, true);
        } else {
            applyEliminationCriteria(residues, false, false);
        }
        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
         * and compute the number of permutations using singleton
         * elimination.
         */
        double permutations = 1;
        double singletonPermutations = 1;
        for (int i = 0; i < nResidues; i++) {
            Residue residue = residues[i];
            Rotamer[] rotamers = residue.getRotamers(library);
            int nr = rotamers.length;
            if (nr > 1) {
                int nrot = 0;
                for (int ri = 0; ri < nr; ri++) {
                    if (!eliminatedSingles[i][ri]) {
                        nrot++;
                    }
                }
                permutations *= rotamers.length;
                if (nrot > 1) {
                    singletonPermutations *= nrot;
                }
            }
        }
        logIfMaster(format(" Collecting Permutations:"));
        logIfMaster(format("%s", " -------------------------------------------------------------------------------------------------------------"));
        dryRun(residues, 0, currentRotamers);
        double pairTotalElimination = singletonPermutations - (double) evaluatedPermutations;
        double afterPairElim = singletonPermutations - pairTotalElimination;
        currentEnsemble = (int) evaluatedPermutations;
        if (ensembleNumber == 1 && currentEnsemble == 0) {
            logger.severe(" No valid path through rotamer space found; try recomputing without pruning or using ensemble.");
        }
        if (ensembleNumber > 1) {
            if (master && printFiles && ensembleFile == null) {
                File file = molecularAssembly.getFile();
                String filename = FilenameUtils.removeExtension(file.getAbsolutePath());
                ensembleFile = new File(filename + ".ens");
                if (ensembleFile.exists()) {
                    for (int i = 2; i < 1000; i++) {
                        ensembleFile = new File(filename + ".ens_" + i);
                        if (!ensembleFile.exists()) {
                            break;
                        }
                    }
                    if (ensembleFile.exists()) {
                        logger.warning(format(" Versioning failed: appending to end of file %s", ensembleFile.getName()));
                    }
                }
                ensembleFilter = new PDBFilter(new File(ensembleFile.getName()), molecularAssembly, null, null);
                logger.info(format(" Ensemble file: %s", ensembleFile.getName()));
            }
            logIfMaster(format(" Ensemble Search Stats: (buffer: %5.3f, current: %d, target: %d)", ensembleBuffer, currentEnsemble, ensembleNumber));
        }
        if (ensembleNumber == 1 || finalTry) {
            logIfMaster(format("%30s %5s %30s %5s %30s %5s", "Condition", "|", "Number of Permutations Left", "|", "Number of Permutations Removed", "|"));
            logIfMaster(format("%s", " -------------------------------------------------------------------------------------------------------------"));
            logIfMaster(format("%30s %5s %30s %5s %30s %5s", "No Eliminations", "|", permutations, "|", "", "|"));
            logIfMaster(format("%30s %5s %30s %5s %30s %5s", "Single Eliminations", "|", singletonPermutations, "|", permutations - singletonPermutations, "|"));
            logIfMaster(format("%30s %5s %30s %5s %30s %5s", "Pair Eliminations", "|", afterPairElim, "|", pairTotalElimination, "|"));
            logIfMaster(format("%30s %5s %30s %5s %30s %5s", "Single and Pair Eliminations", "|", (double) evaluatedPermutations, "|", pairTotalElimination + (permutations - singletonPermutations), "|"));
            logIfMaster(format("%s", " -------------------------------------------------------------------------------------------------------------\n"));
            logIfMaster(format(" Energy of permutations:"));
            logIfMaster(format("%s", " ----------------------------------------------------------------------------------"));
            logIfMaster(format(" %12s %5s %25s %5s %25s %5s", "Permutation", "|", "Energy", "|", "Lowest Possible Energy", "|"));
            logIfMaster(format("%s", " ----------------------------------------------------------------------------------"));
            break;
        }
        if (Math.abs(currentEnsemble - ensembleNumber) < bestEnsembleTargetDiffThusFar) {
            bestEnsembleTargetDiffThusFar = Math.abs(currentEnsemble - ensembleNumber);
            bestBufferThusFar = ensembleBuffer;
        }
        if (currentEnsemble > ensembleNumber) {
            ensembleBuffer -= ensembleBufferStep;
            ensembleBufferStep -= (ensembleBufferStep * 0.01);
            iterations++;
        } else if (currentEnsemble < ensembleNumber) {
            ensembleBuffer += ensembleBufferStep;
            ensembleBufferStep -= (ensembleBufferStep * 0.01);
            iterations++;
        }
        if (iterations > 100) {
            if (currentEnsemble == 0) {
                // TODO: Decide whether we like these next four lines.  Has the potential to produce a crazy amount of permutations.
                logIfMaster(" Ensemble still empty; increasing buffer energy.");
                startingBuffer = 3 * startingBuffer;
                setEnsemble(10, startingBuffer);
                iterations = 0;
            } else {
                ensembleBuffer = bestBufferThusFar;
                finalTry = true;
            }
        }
    }
    if (currentEnsemble == 0) {
        logger.warning(" No valid rotamer permutations found; results will be unreliable.  Try increasing the starting ensemble buffer.");
    }
    double[] permutationEnergyStub = null;
    if (useMonteCarlo()) {
        firstValidPerm(residues, 0, currentRotamers);
        rotamerOptimizationMC(residues, optimum, currentRotamers, nMCsteps, false, mcUseAll);
    } else {
        rotamerOptimizationDEE(molecularAssembly, residues, 0, currentRotamers, Double.MAX_VALUE, optimum, permutationEnergyStub);
    }
    double[] residueEnergy = new double[nResidues];
    double sumSelfEnergy = 0;
    double sumLowSelfEnergy = 0;
    logIfMaster(format("%s", " ----------------------------------------------------------------------------------\n"));
    logIfMaster(format(" Energy contributions:"));
    logIfMaster(format("%s", " -------------------------------------------------------------------------------------"));
    logIfMaster(format(" %15s %5s %25s %5s %25s %5s", "Type", "|", "Energy", "|", "Lowest Possible Energy", "|"));
    logIfMaster(format("%s", " -------------------------------------------------------------------------------------"));
    for (int i = 0; i < nResidues; i++) {
        int ri = optimum[i];
        Residue residue = residues[i];
        Rotamer[] rotamers = residue.getRotamers(library);
        turnOnAtoms(residue);
        RotamerLibrary.applyRotamer(residue, rotamers[ri]);
        double self = getSelf(i, ri);
        residueEnergy[i] = self;
        sumSelfEnergy += self;
        double lowest = lowestSelfEnergy(residues, i);
        sumLowSelfEnergy += lowest;
        if (self - lowest > 10.0) {
            logIfMaster(format(" %15s %5s %25f %5s %25f %5s", "Self (" + residues[i] + "," + ri + "):", "|", self, "|", lowest, "|"));
        }
    }
    double sumPairEnergy = 0.0;
    double sumLowPairEnergy = 0.0;
    double[] resPairEnergy = new double[nResidues];
    double[] lowPairEnergy = new double[nResidues];
    for (int i = 0; i < nResidues - 1; i++) {
        StringBuilder sb = new StringBuilder();
        int ri = optimum[i];
        double sumPairEnergyI = 0;
        double sumLowPairEnergyI = 0;
        for (int j = i + 1; j < nResidues; j++) {
            int rj = optimum[j];
            double pair = get2Body(i, ri, j, rj);
            residueEnergy[i] += 0.5 * pair;
            residueEnergy[j] += 0.5 * pair;
            sumPairEnergy += pair;
            sumPairEnergyI += pair;
            double lowest = lowestPairEnergy(residues, i, ri, j);
            sumLowPairEnergy += lowest;
            sumLowPairEnergyI += lowest;
            resPairEnergy[i] = 0.5 * pair;
            resPairEnergy[j] = 0.5 * pair;
            lowPairEnergy[i] = 0.5 * lowest;
            lowPairEnergy[j] = 0.5 * lowest;
            if (resPairEnergy[i] - lowPairEnergy[i] > 10.0) {
                sb.append(format("  Pair Energy (%8s,%2d) (%8s,%2d): %12.4f (Lowest: %12.4f).\n", residues[i].toFormattedString(false, true), ri, residues[j].toFormattedString(false, true), rj, pair, lowest));
            }
        }
        if (sumPairEnergyI - sumLowPairEnergyI > 10.0) {
            logIfMaster(format(" %15s %5s %25f %5s %25f %5s", "Self (" + residues[i] + "," + ri + "):", "|", sumPairEnergyI, "|", sumLowPairEnergyI, "|"));
            sb.trimToSize();
            if (!sb.toString().isEmpty()) {
                logIfMaster(sb.toString());
            }
        }
    }
    double e = Double.NaN;
    try {
        e = currentEnergy(residueList);
    } catch (ArithmeticException ex) {
        logger.severe(String.format(" Exception %s in calculating current energy at the end of self and pairs", ex.toString()));
    }
    logIfMaster(format(" %15s %5s %25f %5s %25s %5s", "Backbone:", "|", backboneEnergy, "|", "", "|"));
    logIfMaster(format(" %15s %5s %25f %5s %25f %5s", "Self:", "|", sumSelfEnergy, "|", sumLowSelfEnergy, "|"));
    logIfMaster(format(" %15s %5s %25f %5s %25f %5s", "Pair:", "|", sumPairEnergy, "|", sumLowPairEnergy, "|"));
    double approximateEnergy = backboneEnergy + sumSelfEnergy + sumPairEnergy;
    double sumTrimerEnergy = 0;
    if (threeBodyTerm) {
        for (int i = 0; i < nResidues - 2; i++) {
            int ri = optimum[i];
            for (int j = i + 1; j < nResidues - 1; j++) {
                int rj = optimum[j];
                for (int k = j + 1; k < nResidues; k++) {
                    int rk = optimum[k];
                    try {
                        double triple = get3Body(i, ri, j, rj, k, rk);
                        double thirdTrip = triple / 3.0;
                        residueEnergy[i] += thirdTrip;
                        residueEnergy[j] += thirdTrip;
                        residueEnergy[k] += thirdTrip;
                        sumTrimerEnergy += triple;
                    } catch (Exception ex) {
                        logger.warning(ex.toString());
                    }
                }
            }
        }
        approximateEnergy += sumTrimerEnergy;
        double higherOrderEnergy = e - sumSelfEnergy - sumPairEnergy - sumTrimerEnergy - backboneEnergy;
        logIfMaster(format(" %15s %5s %25f %5s %25s %5s", "Trimer:", "|", sumTrimerEnergy, "|", "", "|"));
        logIfMaster(format(" %15s %5s %25f %5s %25s %5s", "Neglected:", "|", higherOrderEnergy, "|", "", "|"));
    } else {
        double higherOrderEnergy = e - sumSelfEnergy - sumPairEnergy - backboneEnergy;
        logIfMaster(format(" %15s %5s %25f %5s %25s %5s", "Neglected:", "|", higherOrderEnergy, "|", "", "|"));
    }
    logIfMaster(format(" %15s %5s %25f %5s %25s %5s", "Approximate:", "|", approximateEnergy, "|", "", "|"));
    logIfMaster(format("%s", " -------------------------------------------------------------------------------------\n"));
    logIfMaster(format(" Final rotamers:"));
    logIfMaster(format("%s", " --------------------------------------------------------------------------------------------"));
    logIfMaster(format("%17s %3s %10s %3s %9s %3s %9s %3s %9s %3s %10s %3s", "Residue", "|", "Chi 1", "|", "Chi 2", "|", "Chi 3", "|", "Chi 4", "|", "Energy", "|"));
    logIfMaster(format("%s", " --------------------------------------------------------------------------------------------"));
    for (int i = 0; i < nResidues; i++) {
        Residue residue = residues[i];
        Rotamer[] rotamers = residue.getRotamers(library);
        int ri = optimum[i];
        Rotamer rotamer = rotamers[ri];
        logIfMaster(format(" %3d %c (%7s,%2d) | %s %12.4f |", i + 1, residue.getChainID(), residue, ri, rotamer.toAngleString(), residueEnergy[i]));
        RotamerLibrary.applyRotamer(residue, rotamer);
    }
    logIfMaster(format("%s", " --------------------------------------------------------------------------------------------\n"));
    return e;
}

Example 8 with ResidueState

use of ffx.potential.bonded.ResidueState 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 9 with ResidueState

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

Example 10 with ResidueState

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

the class RotamerOptimization method eliminateNABackboneRotamers.

/**
 * Eliminates NA backbone rotamers with corrections greater than threshold.
 * The int[] parameter allows the method to know how many Rotamers for each
 * residue have previously been pruned; currently, this means any Rotamer
 * pruned by reconcileNARotamersWithPriorResidues.
 * <p>
 * A nucleic correction threshold of 0 skips the entire method; this check
 * is presently being performed inside the method in case it is called again
 * at some point.
 *
 * @param residues Residues to eliminate bad backbone rotamers over.
 * @param numEliminatedRotamers Number of previously eliminated rotamers per
 * residue.
 */
private void eliminateNABackboneRotamers(Residue[] residues, int[] numEliminatedRotamers) {
    /* Atom atoms[] = molecularAssembly.getAtomArray();
         int nAtoms = atoms.length;
         String begin[] = new String[nAtoms];
         for (int i = 0; i < nAtoms; i++) {
         begin[i] = atoms[i].toString();
         } */
    if (nucleicCorrectionThreshold != 0) {
        logIfMaster(format(" Eliminating nucleic acid rotamers with correction vectors larger than %5.3f A", nucleicCorrectionThreshold));
        logIfMaster(format(" A minimum of %d rotamers per NA residue will carry through to energy calculations.", minNumberAcceptedNARotamers));
        ArrayList<Residue> resList = new ArrayList<>();
        resList.addAll(Arrays.asList(residues));
        ResidueState[] origCoordinates = ResidueState.storeAllCoordinates(resList);
        for (int j = 0; j < residues.length; j++) {
            Residue nucleicResidue = residues[j];
            Rotamer[] rotamers = nucleicResidue.getRotamers(library);
            if (nucleicResidue.getResidueType() == NA && rotamers != null) {
                int nrotamers = rotamers.length;
                // Default to all rotamers that have not previously been
                // eliminated; subtract as rotamers are rejected.
                int numAcceptedRotamers = nrotamers - numEliminatedRotamers[j];
                if (minNumberAcceptedNARotamers >= numAcceptedRotamers) {
                    continue;
                }
                ArrayList<DoubleIndexPair> rejectedRotamers = new ArrayList<>();
                for (int i = 0; i < nrotamers; i++) {
                    if (!check(j, i)) {
                        try {
                            RotamerLibrary.applyRotamer(nucleicResidue, rotamers[i], nucleicCorrectionThreshold);
                        } catch (NACorrectionException error) {
                            double rejectedCorrection = error.getCorrection();
                            numAcceptedRotamers--;
                            DoubleIndexPair rejected = new DoubleIndexPair(i, rejectedCorrection);
                            rejectedRotamers.add(rejected);
                        }
                    }
                }
                int numAdditionalRotamersToAccept = minNumberAcceptedNARotamers - numAcceptedRotamers;
                if (numAdditionalRotamersToAccept > 0) {
                    DoubleIndexPair[] rejectedArray = new DoubleIndexPair[rejectedRotamers.size()];
                    for (int i = 0; i < rejectedArray.length; i++) {
                        rejectedArray[i] = rejectedRotamers.get(i);
                    }
                    Arrays.sort(rejectedArray);
                    rejectedRotamers = new ArrayList<>();
                    rejectedRotamers.addAll(Arrays.asList(rejectedArray));
                    for (int i = 0; i < numAdditionalRotamersToAccept; i++) {
                        rejectedRotamers.remove(0);
                    }
                }
                for (DoubleIndexPair rotToReject : rejectedRotamers) {
                    eliminateRotamer(residues, j, rotToReject.getIndex(), print);
                    logIfMaster(format(" Correction magnitude was %6.4f A > %5.3f A", rotToReject.getDoubleValue(), nucleicCorrectionThreshold));
                }
            }
            nucleicResidue.revertState(origCoordinates[j]);
        // revertSingleResidueCoordinates(nucleicResidue, originalCoordinates[j]);
        }
    }
}