Search in sources :

Example 11 with MultiResidue

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

the class TitrationUtils method titratingMultiresidueFactory.

/**
 * Create a MultiResidue from the given Residue by adding its alternated
 * protonation state(s) as alternate possibilities.
 */
public static MultiResidue titratingMultiresidueFactory(MolecularAssembly mola, Residue res) {
    ForceField ff = mola.getForceField();
    Potential potential = mola.getPotentialEnergy();
    if (!(potential instanceof ForceFieldEnergy)) {
        logger.warning(String.format("TitrationFactory only supported by ForceFieldEnergy potentials."));
        throw new IllegalStateException();
    }
    ForceFieldEnergy ffe = (ForceFieldEnergy) potential;
    /* Create new titration state. */
    Titration titration = Titration.lookup(res);
    String targetName = (titration.protForm != res.getAminoAcid3()) ? titration.protForm.toString() : titration.deprotForm.toString();
    int resNumber = res.getResidueNumber();
    Residue.ResidueType resType = res.getResidueType();
    Residue newRes = new Residue(targetName, resNumber, resType);
    /* Wrap both states in a MultiResidue. */
    MultiResidue multiRes = new MultiResidue(res, ff, ffe);
    Polymer polymer = findResiduePolymer(res, mola);
    polymer.addMultiResidue(multiRes);
    multiRes.addResidue(newRes);
    /* Begin in protonated state by default. */
    multiRes.setActiveResidue(titration.protForm);
    propagateInactiveResidues(multiRes, false);
    ffe.reInit();
    return multiRes;
}
Also used : MultiResidue(ffx.potential.bonded.MultiResidue) Residue(ffx.potential.bonded.Residue) Polymer(ffx.potential.bonded.Polymer) ForceField(ffx.potential.parameters.ForceField) ForceFieldEnergy(ffx.potential.ForceFieldEnergy) Potential(ffx.numerics.Potential) MultiResidue(ffx.potential.bonded.MultiResidue)

Example 12 with MultiResidue

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

the class RotamerOptimization method slidingWindowOptimization.

private double slidingWindowOptimization(ArrayList<Residue> residueList, int windowSize, int increment, boolean revert, double distance, Direction direction) {
    long beginTime = -System.nanoTime();
    boolean incrementTruncated = false;
    boolean firstWindowSaved = false;
    int counter = 1;
    int windowEnd;
    int nOptimize = residueList.size();
    if (nOptimize < windowSize) {
        windowSize = nOptimize;
        logger.info(" Window size too small for given residue range; truncating window size.");
    }
    switch(direction) {
        case BACKWARD:
            ArrayList<Residue> temp = new ArrayList<>();
            for (int i = nOptimize - 1; i >= 0; i--) {
                temp.add(residueList.get(i));
            }
            residueList = temp;
        // Fall through into the FORWARD case.
        case FORWARD:
            for (int windowStart = 0; windowStart + (windowSize - 1) < nOptimize; windowStart += increment) {
                long windowTime = -System.nanoTime();
                windowEnd = windowStart + (windowSize - 1);
                logIfMaster(format("\n Iteration %d of the sliding window.\n", counter++));
                Residue firstResidue = residueList.get(windowStart);
                Residue lastResidue = residueList.get(windowEnd);
                if (firstResidue != lastResidue) {
                    logIfMaster(format(" Residues %s ... %s", firstResidue.toString(), lastResidue.toString()));
                } else {
                    logIfMaster(format(" Residue %s", firstResidue.toString()));
                }
                ArrayList<Residue> currentWindow = new ArrayList<>();
                // Not filled if useForcedResidues == false.
                ArrayList<Residue> onlyRotameric = new ArrayList<>();
                for (int i = windowStart; i <= windowEnd; i++) {
                    Residue residue = residueList.get(i);
                    if (useForcedResidues && residue.getRotamers(library) != null) {
                        onlyRotameric.add(residue);
                    }
                    currentWindow.add(residueList.get(i));
                }
                if (distance > 0) {
                    for (int i = windowStart; i <= windowEnd; i++) {
                        Residue residuei = residueList.get(i);
                        int indexI = allResiduesList.indexOf(residuei);
                        int lengthRi;
                        if (checkIfForced(residuei)) {
                            lengthRi = 1;
                        } else {
                            lengthRi = residuei.getRotamers(library).length;
                        }
                        for (int ri = 0; ri < lengthRi; ri++) {
                            for (int j = 0; j < numResidues; j++) {
                                Residue residuej = allResiduesArray[j];
                                Rotamer[] rotamersj = residuej.getRotamers(library);
                                if (currentWindow.contains(residuej) || rotamersj == null) {
                                    continue;
                                }
                                int lengthRj = rotamersj.length;
                                for (int rj = 0; rj < lengthRj; rj++) {
                                    double rotamerSeparation = get2BodyDistance(indexI, ri, j, rj);
                                    // if (distanceMatrix[indexI][ri][j][rj] <= distance) {
                                    if (rotamerSeparation <= distance) {
                                        if (!currentWindow.contains(residuej)) {
                                            logIfMaster(format(" Adding residue %s at distance %16.8f Ang from %s %d.", residuej.toFormattedString(false, true), rotamerSeparation, residuei.toFormattedString(false, true), ri));
                                            currentWindow.add(residuej);
                                            if (useForcedResidues) {
                                                onlyRotameric.add(residuej);
                                            }
                                        }
                                        break;
                                    }
                                }
                            }
                        }
                    }
                }
                /**
                 * If the window starts with a nucleic acid, and there is a
                 * 5' NA residue, ensure that 5' NA residue has been
                 * included in the window. Otherwise, that previous residue
                 * may not have had a chance to be flexible about its sugar
                 * pucker.
                 *
                 * If window size is greater than increment, however, this
                 * has already been handled. Additionally, do not perform
                 * this for the first window (counter is already incremented
                 * by the time this check is performed, so first window's
                 * counter will be 2). Furthermore, do not include Residues
                 * with null Rotamer lists (this breaks things).
                 *
                 * The issue: if window size = increment, the last NA
                 * residue in each window will not have flexibility about
                 * its sugar pucker, because its self-energy includes the
                 * O3' (i) to P (i+1) bond, so it must remain in the
                 * original sugar pucker to meet the i+1 residue. However,
                 * this problem can be solved by ensuring that final residue
                 * is included in the next window, where it will have
                 * flexibility about its sugar pucker.
                 *
                 * If you are running successive sliding window jobs on the
                 * same file, I would suggest starting the next job on the
                 * last residue of the previous job, unless you know your
                 * settings will include it.
                 */
                if (counter > 2 && windowSize <= increment && firstResidue.getResidueType() == NA) {
                    Residue prevResidue = firstResidue.getPreviousResidue();
                    if (prevResidue != null && prevResidue.getResidueType() == NA && !currentWindow.contains(prevResidue) && prevResidue.getRotamers(library) != null) {
                        logIfMaster(format(" Adding nucleic acid residue 5' of window start %s to give it flexibility about its sugar pucker.", prevResidue.toString()));
                        currentWindow.add(prevResidue);
                        if (useForcedResidues) {
                            onlyRotameric.add(prevResidue);
                        }
                    }
                }
                if (useForcedResidues) {
                    sortResidues(onlyRotameric);
                } else {
                    sortResidues(currentWindow);
                }
                if (revert) {
                    ResidueState[] coordinates = ResidueState.storeAllCoordinates(currentWindow);
                    // 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 = Double.NaN;
                    try {
                        startingEnergy = currentEnergy(currentWindow);
                    } catch (ArithmeticException ex) {
                        logger.severe(String.format(" Exception %s in calculating starting energy of a window; FFX shutting down", ex.toString()));
                    }
                    if (useForcedResidues) {
                        if (onlyRotameric.size() < 1) {
                            logger.info(" Window has no rotameric residues.");
                            ResidueState.revertAllCoordinates(currentWindow, coordinates);
                        } else {
                            globalOptimization(onlyRotameric);
                            double finalEnergy = Double.NaN;
                            try {
                                finalEnergy = currentEnergy(currentWindow);
                            } catch (ArithmeticException ex) {
                                logger.severe(String.format(" Exception %s in calculating final energy of a window; FFX shutting down", ex.toString()));
                            }
                            if (startingEnergy <= finalEnergy) {
                                logger.warning("Optimization did not yield a better energy. Reverting to orginal coordinates.");
                                ResidueState.revertAllCoordinates(currentWindow, coordinates);
                            }
                        }
                    } else {
                        globalOptimization(currentWindow);
                        double finalEnergy = Double.NaN;
                        try {
                            finalEnergy = currentEnergy(currentWindow);
                        } catch (ArithmeticException ex) {
                            logger.severe(String.format(" Exception %s in calculating final energy of a window; FFX shutting down", ex.toString()));
                        }
                        if (startingEnergy <= finalEnergy) {
                            logger.warning("Optimization did not yield a better energy. Reverting to orginal coordinates.");
                            ResidueState.revertAllCoordinates(currentWindow, coordinates);
                        }
                    }
                } else if (useForcedResidues) {
                    if (onlyRotameric.size() < 1) {
                        logger.info(" Window has no rotameric residues.");
                    } else {
                        globalOptimization(onlyRotameric);
                    }
                } else {
                    globalOptimization(currentWindow);
                }
                if (!incrementTruncated) {
                    if (windowStart + (windowSize - 1) + increment > nOptimize - 1) {
                        increment = nOptimize - windowStart - windowSize;
                        if (increment == 0) {
                            break;
                        }
                        logger.warning(" Increment truncated in order to optimize entire residue range.");
                        incrementTruncated = true;
                    }
                }
                if (master && printFiles) {
                    File file = molecularAssembly.getFile();
                    if (firstWindowSaved) {
                        file.delete();
                    }
                    // Don't write a file if its the final iteration.
                    if (windowStart + windowSize == nOptimize) {
                        continue;
                    }
                    // String filename = FilenameUtils.removeExtension(file.getAbsolutePath());
                    // File tempFile = new File(filename + ".win");
                    // PDBFilter windowFilter = new PDBFilter(new File(filename + ".win"), molecularAssembly, null, null);
                    PDBFilter windowFilter = new PDBFilter(file, molecularAssembly, null, null);
                    // StringBuilder header = new StringBuilder(format("Iteration %d of the sliding window\n", counter - 1));
                    try {
                        windowFilter.writeFile(file, false);
                        if (firstResidue != lastResidue) {
                            logger.info(format(" File with residues %s ... %s in window written to.", firstResidue.toString(), lastResidue.toString()));
                        } else {
                            logger.info(format(" File with residue %s in window written to.", firstResidue.toString()));
                        }
                    } catch (Exception e) {
                        logger.warning(format("Exception writing to file: %s", file.getName()));
                    }
                    firstWindowSaved = true;
                }
                long currentTime = System.nanoTime();
                windowTime += currentTime;
                logIfMaster(format(" Time elapsed for this iteration: %11.3f sec", windowTime * 1.0E-9));
                logIfMaster(format(" Overall time elapsed: %11.3f sec", (currentTime + beginTime) * 1.0E-9));
            /*for (Residue residue : residueList) {
                        if (residue instanceof MultiResidue) {
                            ((MultiResidue) residue).setDefaultResidue();
                            residue.reInitOriginalAtomList();
                        }
                    }*/
            }
            break;
        default:
            // No default case.
            break;
    }
    return 0.0;
}
Also used : ResidueState(ffx.potential.bonded.ResidueState) ArrayList(java.util.ArrayList) RotamerLibrary.applyRotamer(ffx.potential.bonded.RotamerLibrary.applyRotamer) Rotamer(ffx.potential.bonded.Rotamer) Atom(ffx.potential.bonded.Atom) IOException(java.io.IOException) NACorrectionException(ffx.potential.bonded.NACorrectionException) Residue(ffx.potential.bonded.Residue) MultiResidue(ffx.potential.bonded.MultiResidue) File(java.io.File) PDBFilter(ffx.potential.parsers.PDBFilter)

Example 13 with MultiResidue

use of ffx.potential.bonded.MultiResidue 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;
}
Also used : ResidueState(ffx.potential.bonded.ResidueState) ArrayList(java.util.ArrayList) Atom(ffx.potential.bonded.Atom) IOException(java.io.IOException) NACorrectionException(ffx.potential.bonded.NACorrectionException) Residue(ffx.potential.bonded.Residue) MultiResidue(ffx.potential.bonded.MultiResidue) File(java.io.File) PDBFilter(ffx.potential.parsers.PDBFilter) Crystal(ffx.crystal.Crystal)

Example 14 with MultiResidue

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

the class RotamerOptimization method pruneSingleClashes.

/**
 * Uses calculated energies to prune rotamers based on a threshold distance
 * from that residue's minimum energy rotamer (by default 20 kcal/mol). The
 * threshold can be modulated by presence of nucleic acids or MultiResidues,
 * which require more generous pruning criteria.
 *
 * @param residues Residues to prune rotamers over.
 */
private void pruneSingleClashes(Residue[] residues) {
    if (!pruneClashes) {
        return;
    }
    for (int i = 0; i < residues.length; i++) {
        Residue residue = residues[i];
        Rotamer[] rotamers = residue.getRotamers(library);
        int nrot = rotamers.length;
        double minEnergy = Double.MAX_VALUE;
        int minRot = -1;
        for (int ri = 0; ri < nrot; ri++) {
            if (!check(i, ri) && getSelf(i, ri) < minEnergy) {
                minEnergy = getSelf(i, ri);
                minRot = ri;
            }
        }
        /**
         * Regular: ep = minEnergy + clashThreshold Nucleic acids: ep =
         * minEnergy + (clashThreshold * factor * factor) MultiResidues: ep
         * = minEnergy + multiResClashThreshold
         *
         * Nucleic acids are bigger than amino acids, and MultiResidues can
         * have wild swings in energy on account of chemical perturbation.
         */
        double energyToPrune = (residue instanceof MultiResidue) ? multiResClashThreshold : clashThreshold;
        energyToPrune = (residue.getResidueType() == NA) ? energyToPrune * singletonNAPruningFactor * pruningFactor : energyToPrune;
        energyToPrune += minEnergy;
        for (int ri = 0; ri < nrot; ri++) {
            if (!check(i, ri) && (getSelf(i, ri) > energyToPrune)) {
                if (eliminateRotamer(residues, i, ri, print)) {
                    logIfMaster(format("  Rotamer (%7s,%2d) self-energy %s pruned by (%7s,%2d) %s.", residue, ri, formatEnergy(getSelf(i, ri)), residue, minRot, formatEnergy(minEnergy)));
                }
            }
        }
    }
}
Also used : Residue(ffx.potential.bonded.Residue) MultiResidue(ffx.potential.bonded.MultiResidue) RotamerLibrary.applyRotamer(ffx.potential.bonded.RotamerLibrary.applyRotamer) Rotamer(ffx.potential.bonded.Rotamer) MultiResidue(ffx.potential.bonded.MultiResidue)

Example 15 with MultiResidue

use of ffx.potential.bonded.MultiResidue 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()));
        }
    }
}
Also used : SymOp(ffx.crystal.SymOp) ParallelTeam(edu.rit.pj.ParallelTeam) NeighborList(ffx.potential.nonbonded.NeighborList) ResidueState(ffx.potential.bonded.ResidueState) Atom(ffx.potential.bonded.Atom) IOException(java.io.IOException) NACorrectionException(ffx.potential.bonded.NACorrectionException) Residue(ffx.potential.bonded.Residue) MultiResidue(ffx.potential.bonded.MultiResidue) MultiResidue(ffx.potential.bonded.MultiResidue) Crystal(ffx.crystal.Crystal)

Aggregations

MultiResidue (ffx.potential.bonded.MultiResidue)19 Residue (ffx.potential.bonded.Residue)18 Atom (ffx.potential.bonded.Atom)10 Rotamer (ffx.potential.bonded.Rotamer)6 TitrationUtils.inactivateResidue (ffx.potential.extended.TitrationUtils.inactivateResidue)6 ArrayList (java.util.ArrayList)6 ResidueState (ffx.potential.bonded.ResidueState)5 AminoAcid3 (ffx.potential.bonded.ResidueEnumerations.AminoAcid3)4 NACorrectionException (ffx.potential.bonded.NACorrectionException)3 Polymer (ffx.potential.bonded.Polymer)3 RotamerLibrary.applyRotamer (ffx.potential.bonded.RotamerLibrary.applyRotamer)3 Titration (ffx.potential.extended.TitrationUtils.Titration)3 TitrationUtils.performTitration (ffx.potential.extended.TitrationUtils.performTitration)3 File (java.io.File)3 IOException (java.io.IOException)3 Crystal (ffx.crystal.Crystal)2 MultiTerminus (ffx.potential.bonded.MultiTerminus)2 TitrationType (ffx.potential.extended.TitrationUtils.TitrationType)2 PDBFilter (ffx.potential.parsers.PDBFilter)2 ParallelTeam (edu.rit.pj.ParallelTeam)1