Search in sources :

Example 1 with ResidueState

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

the class RotamerOptimization method evaluateDistance.

/**
 * Evaluates the pairwise distance between two residues' rotamers under any
 * symmetry operator; does "lazy loading" for the distance matrix.
 *
 * @param i Residue i
 * @param ri Rotamer for i
 * @param j Residue j
 * @param rj Rotamer for j
 * @return Shortest distance
 */
private double evaluateDistance(int i, int ri, int j, int rj) {
    Residue resi = allResiduesArray[i];
    Rotamer[] rotamersI = resi.getRotamers(library);
    Rotamer roti = rotamersI[ri];
    double[][] xi;
    if (roti.equals(resi.getRotamer())) {
        xi = resi.storeCoordinateArray();
    } else {
        ResidueState origI = resi.storeState();
        RotamerLibrary.applyRotamer(resi, roti);
        xi = resi.storeCoordinateArray();
        resi.revertState(origI);
    }
    Residue resj = allResiduesArray[j];
    Rotamer[] rotamersJ = resj.getRotamers(library);
    Rotamer rotj = rotamersJ[rj];
    double[][] xj;
    if (rotj.equals(resj.getRotamer())) {
        xj = resj.storeCoordinateArray();
    } else {
        ResidueState origJ = resj.storeState();
        RotamerLibrary.applyRotamer(resj, rotj);
        xj = resj.storeCoordinateArray();
        resj.revertState(origJ);
    }
    Crystal crystal = molecularAssembly.getCrystal();
    int nSymm = crystal.spaceGroup.getNumberOfSymOps();
    double minDist = Double.MAX_VALUE;
    for (int iSymOp = 0; iSymOp < nSymm; iSymOp++) {
        SymOp symOp = crystal.spaceGroup.getSymOp(iSymOp);
        double dist = interResidueDistance(xi, xj, symOp);
        minDist = dist < minDist ? dist : minDist;
    }
    return minDist;
}
Also used : SymOp(ffx.crystal.SymOp) Residue(ffx.potential.bonded.Residue) MultiResidue(ffx.potential.bonded.MultiResidue) ResidueState(ffx.potential.bonded.ResidueState) RotamerLibrary.applyRotamer(ffx.potential.bonded.RotamerLibrary.applyRotamer) Rotamer(ffx.potential.bonded.Rotamer) Crystal(ffx.crystal.Crystal)

Example 2 with ResidueState

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

the class RotamerOptimization method rotamerOptimizationDEE.

/**
 * A global optimization over side-chain rotamers using a recursive
 * algorithm and information about eliminated rotamers, rotamer pairs and
 * rotamer triples
 *
 * @param molecularAssembly
 * @param residues
 * @param i
 * @param currentRotamers
 * @param lowEnergy
 * @param optimum Optimum set of rotamers.
 * @param permutationEnergies Energies of visited permutations or null.
 * @return current energy.
 */
private double rotamerOptimizationDEE(MolecularAssembly molecularAssembly, Residue[] residues, int i, int[] currentRotamers, double lowEnergy, int[] optimum, double[] permutationEnergies) {
    // This is the initialization condition.
    if (i == 0) {
        evaluatedPermutations = 0;
    }
    int nResidues = residues.length;
    Residue residuei = residues[i];
    Rotamer[] rotamersi = residuei.getRotamers(library);
    int lenri = rotamersi.length;
    double currentEnergy = Double.MAX_VALUE;
    List<Residue> resList = Arrays.asList(residues);
    /**
     * As long as there are more residues, continue the recursion for each
     * rotamer of the current residue.
     */
    if (i < nResidues - 1) {
        /**
         * Loop over rotamers of residue i.
         */
        for (int ri = 0; ri < lenri; ri++) {
            /**
             * Check if rotamer ri has been eliminated by DEE.
             */
            if (check(i, ri)) {
                continue;
            }
            /**
             * Check if rotamer ri has been eliminated by an upstream
             * rotamer (any residue's rotamer from j = 0 .. i-1).
             */
            boolean deadEnd = false;
            for (int j = 0; j < i; j++) {
                int rj = currentRotamers[j];
                deadEnd = check(j, rj, i, ri);
                if (deadEnd) {
                    break;
                }
            }
            if (deadEnd) {
                continue;
            }
            applyRotamer(residuei, rotamersi[ri]);
            currentRotamers[i] = ri;
            double rotEnergy = rotamerOptimizationDEE(molecularAssembly, residues, i + 1, currentRotamers, lowEnergy, optimum, permutationEnergies);
            if (rotEnergy < currentEnergy) {
                currentEnergy = rotEnergy;
            }
            if (rotEnergy < lowEnergy) {
                optimum[i] = ri;
                lowEnergy = rotEnergy;
            }
        }
    } else {
        if (ensembleStates == null) {
            ensembleStates = new ArrayList<>();
        }
        /**
         * At the end of the recursion, compute the potential energy for
         * each rotamer of the final residue. If a lower potential energy is
         * discovered, the rotamers of each residue will be collected as the
         * recursion returns up the chain.
         */
        for (int ri = 0; ri < lenri; ri++) {
            /**
             * Check if rotamer ri has been eliminated by DEE.
             */
            if (check(i, ri)) {
                continue;
            }
            currentRotamers[i] = ri;
            /**
             * Check if rotamer ri has been eliminated by an upstream
             * rotamer (any residue's rotamer from 0 .. i-1.
             */
            boolean deadEnd = false;
            for (int j = 0; j < i; j++) {
                int rj = currentRotamers[j];
                deadEnd = check(j, rj, i, ri);
                if (deadEnd) {
                    break;
                }
            }
            if (deadEnd) {
                continue;
            }
            applyRotamer(residuei, rotamersi[ri]);
            // Compute the energy based on a 3-body approximation
            double approximateEnergy = computeEnergy(residues, currentRotamers, false);
            double comparisonEnergy = approximateEnergy;
            evaluatedPermutations++;
            // Compute the AMOEBA energy
            if (useFullAMOEBAEnergy) {
                double amoebaEnergy = Double.NaN;
                try {
                    amoebaEnergy = currentEnergy(resList);
                } catch (ArithmeticException ex) {
                    logger.warning(String.format(" Exception %s in calculating full AMOEBA energy for permutation %d", ex.toString(), evaluatedPermutations));
                }
                comparisonEnergy = amoebaEnergy;
            }
            if (permutationEnergies != null) {
                permutationEnergies[evaluatedPermutations - 1] = comparisonEnergy;
            }
            if (algorithmListener != null) {
                algorithmListener.algorithmUpdate(molecularAssembly);
            }
            if (ensembleNumber > 1) {
                if (master && printFiles) {
                    try {
                        FileWriter fw = new FileWriter(ensembleFile, true);
                        BufferedWriter bw = new BufferedWriter(fw);
                        bw.write(format("MODEL        %d", evaluatedPermutations));
                        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 e) {
                        logger.warning(format("Exception writing to file: %s", ensembleFile.getName()));
                    }
                }
                ResidueState[] states = ResidueState.storeAllCoordinates(residues);
                ensembleStates.add(new ObjectPair<>(states, comparisonEnergy));
            }
            if (comparisonEnergy < currentEnergy) {
                currentEnergy = comparisonEnergy;
            }
            if (comparisonEnergy < lowEnergy) {
                lowEnergy = comparisonEnergy;
                optimum[i] = ri;
            }
            if (useFullAMOEBAEnergy) {
                // Log current results
                logIfMaster(format(" %6e AMOEBA: %12.4f 3-Body: %12.4f Neglected: %12.4f (%12.4f)", (double) evaluatedPermutations, comparisonEnergy, approximateEnergy, comparisonEnergy - approximateEnergy, lowEnergy));
            } else {
                if (threeBodyTerm) {
                    ;
                    logIfMaster(format(" %12s %5s %25f %5s %25f %5s", evaluatedPermutations, "|", approximateEnergy, "|", lowEnergy, "|"));
                // logIfMaster(format(" %6e 3-Body: %12.4f (%12.4f)",
                // (double) evaluatedPermutations, approximateEnergy, lowEnergy));
                } else {
                    logIfMaster(format(" %12s %5s %25f %5s %25f %5s", evaluatedPermutations, "|", approximateEnergy, "|", lowEnergy, "|"));
                // logIfMaster(format(" %6e 2-Body: %12.4f (%12.4f)",
                // (double) evaluatedPermutations, approximateEnergy, lowEnergy));
                }
            }
        }
        ensembleStates.sort(null);
    }
    return currentEnergy;
}
Also used : ResidueState(ffx.potential.bonded.ResidueState) FileWriter(java.io.FileWriter) RotamerLibrary.applyRotamer(ffx.potential.bonded.RotamerLibrary.applyRotamer) Rotamer(ffx.potential.bonded.Rotamer) IOException(java.io.IOException) BufferedWriter(java.io.BufferedWriter) Residue(ffx.potential.bonded.Residue) MultiResidue(ffx.potential.bonded.MultiResidue)

Example 3 with ResidueState

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

the class PhMD method tryComboStep.

/**
 * Attempt a combination titration/rotamer MC move.
 *
 * @param targetMulti
 * @return accept/reject
 */
private boolean tryComboStep(MultiResidue targetMulti) {
    if (CAUTIOUS) {
        throw new UnsupportedOperationException();
    }
    // Record the pre-change total energy.
    double previousTotalEnergy = currentTotalEnergy();
    double previousElectrostaticEnergy = currentElectrostaticEnergy();
    // Write the pre-combo snapshot.
    writeSnapshot(true, StepType.COMBO, config.snapshots);
    String startString = targetMulti.toString();
    String startName = targetMulti.getActive().getName();
    // Choose from the list of available titrations for the active residue.
    List<Titration> avail = titrationMap.get(targetMulti.getActive());
    Titration titration = avail.get(rng.nextInt(avail.size()));
    // Perform the chosen titration.
    TitrationType titrationType = performTitration(targetMulti, titration, config.inactivateBackground);
    // Change rotamer state, but first save coordinates so we can return to them if rejected.
    Residue residue = targetMulti.getActive();
    ArrayList<Atom> atoms = residue.getAtomList();
    ResidueState origState = residue.storeState();
    double[] chi = new double[4];
    RotamerLibrary.measureAARotamer(residue, chi, false);
    AminoAcid3 aa = AminoAcid3.valueOf(residue.getName());
    Rotamer origCoordsRotamer = new Rotamer(aa, origState, chi[0], 0, chi[1], 0, chi[2], 0, chi[3], 0);
    // Swap to the new rotamer.
    Rotamer[] rotamers = residue.getRotamers(library);
    int rotaRand = rng.nextInt(rotamers.length);
    RotamerLibrary.applyRotamer(residue, rotamers[rotaRand]);
    // Write the post-combo snapshot.
    writeSnapshot(false, StepType.COMBO, config.snapshots);
    // Evaluate both MC criteria.
    String endName = targetMulti.getActive().getName();
    // Evaluate the titration probability of the step.
    double pKaref = titration.pKa;
    double dG_ref = titration.refEnergy;
    double temperature = thermostat.getCurrentTemperature();
    double kT = BOLTZMANN * temperature;
    double dG_elec = currentElectrostaticEnergy() - previousElectrostaticEnergy;
    if (config.zeroReferenceEnergies) {
        dG_ref = 0.0;
    }
    double prefix = Math.log(10) * kT * (pH - pKaref);
    if (titrationType == TitrationType.DEPROT) {
        prefix = -prefix;
    }
    double postfix = dG_elec - dG_ref;
    double dG_titr = prefix + postfix;
    double titrCriterion = exp(-dG_titr / kT);
    // Evaluate the rotamer probability of the step.
    double dG_rota = currentTotalEnergy() - previousTotalEnergy;
    double rotaCriterion = exp(-dG_rota / kT);
    StringBuilder sb = new StringBuilder();
    sb.append(String.format(" Assessing possible MC combo step:\n"));
    sb.append(String.format("     dG_elec: %16.8f\n", dG_elec));
    sb.append(String.format("     dG_titr: %16.8f\n", dG_titr));
    sb.append(String.format("     dG_rota: %16.8f\n", dG_rota));
    sb.append(String.format("     -----\n"));
    // Automatic acceptance if both energy changes are favorable.
    if (dG_titr < 0 && dG_rota < 0 && config.mcOverride != MCOverride.REJECT) {
        sb.append(String.format("     Accepted!"));
        logger.info(sb.toString());
        numMovesAccepted++;
        propagateInactiveResidues(titratingMultis, false);
        return true;
    } else {
        // Conditionally accept based on combined probabilities.
        if (dG_titr < 0 || config.mcOverride == MCOverride.ACCEPT) {
            titrCriterion = 1.0;
        }
        if (dG_rota < 0) {
            rotaCriterion = 1.0;
        }
        if (config.mcOverride == MCOverride.REJECT) {
            titrCriterion = 0.0;
        }
        double metropolis = random();
        double comboCriterion = titrCriterion * rotaCriterion;
        sb.append(String.format("     titrCrit:   %9.4f\n", titrCriterion));
        sb.append(String.format("     rotaCrit:   %9.4f\n", rotaCriterion));
        sb.append(String.format("     criterion:  %9.4f\n", comboCriterion));
        sb.append(String.format("     rng:        %9.4f\n", metropolis));
        if (metropolis < comboCriterion) {
            sb.append(String.format("     Accepted!"));
            logger.info(sb.toString());
            numMovesAccepted++;
            propagateInactiveResidues(titratingMultis, false);
            return true;
        } else {
            // Move was denied.
            sb.append(String.format("     Denied."));
            logger.info(sb.toString());
            // Undo both pieces of the rejected move IN THE RIGHT ORDER.
            RotamerLibrary.applyRotamer(residue, origCoordsRotamer);
            performTitration(targetMulti, titration, config.inactivateBackground);
            ffe.reInit();
            molDyn.reInit();
            return false;
        }
    }
}
Also used : ResidueState(ffx.potential.bonded.ResidueState) AminoAcid3(ffx.potential.bonded.ResidueEnumerations.AminoAcid3) TitrationType(ffx.potential.extended.TitrationUtils.TitrationType) Rotamer(ffx.potential.bonded.Rotamer) Atom(ffx.potential.bonded.Atom) TitrationUtils.inactivateResidue(ffx.potential.extended.TitrationUtils.inactivateResidue) MultiResidue(ffx.potential.bonded.MultiResidue) Residue(ffx.potential.bonded.Residue) Titration(ffx.potential.extended.TitrationUtils.Titration) TitrationUtils.performTitration(ffx.potential.extended.TitrationUtils.performTitration)

Example 4 with ResidueState

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

the class PhMD method tryRotamerStep.

/**
 * Attempt a rotamer MC move.
 *
 * @param targetMulti
 * @return accept/reject
 */
private boolean tryRotamerStep(MultiResidue targetMulti) {
    if (CAUTIOUS) {
        throw new UnsupportedOperationException();
    }
    // Record the pre-change total energy.
    double previousTotalEnergy = currentTotalEnergy();
    // Write the before-step snapshot.
    writeSnapshot(true, StepType.ROTAMER, config.snapshots);
    // Save coordinates so we can return to them if move is rejected.
    Residue residue = targetMulti.getActive();
    ArrayList<Atom> atoms = residue.getAtomList();
    ResidueState origState = residue.storeState();
    double[] chi = new double[4];
    RotamerLibrary.measureAARotamer(residue, chi, false);
    AminoAcid3 aa = AminoAcid3.valueOf(residue.getName());
    Rotamer origCoordsRotamer = new Rotamer(aa, origState, chi[0], 0, chi[1], 0, chi[2], 0, chi[3], 0);
    // Select a new rotamer and swap to it.
    Rotamer[] rotamers = residue.getRotamers(library);
    int rotaRand = rng.nextInt(rotamers.length);
    RotamerLibrary.applyRotamer(residue, rotamers[rotaRand]);
    // Write the post-rotamer change snapshot.
    writeSnapshot(false, StepType.ROTAMER, config.snapshots);
    // Check the MC criterion.
    double temperature = thermostat.getCurrentTemperature();
    double kT = BOLTZMANN * temperature;
    double postTotalEnergy = currentTotalEnergy();
    double dG_tot = postTotalEnergy - previousTotalEnergy;
    double criterion = exp(-dG_tot / kT);
    StringBuilder sb = new StringBuilder();
    sb.append(String.format(" Assessing possible MC rotamer step:\n"));
    sb.append(String.format("     prev:   %16.8f\n", previousTotalEnergy));
    sb.append(String.format("     post:   %16.8f\n", postTotalEnergy));
    sb.append(String.format("     dG_tot: %16.8f\n", dG_tot));
    sb.append(String.format("     -----\n"));
    // Automatic acceptance if energy change is favorable.
    if (dG_tot < 0) {
        sb.append(String.format("     Accepted!"));
        logger.info(sb.toString());
        numMovesAccepted++;
        propagateInactiveResidues(titratingMultis, true);
        return true;
    } else {
        // Conditional acceptance if energy change is positive.
        double metropolis = random();
        sb.append(String.format("     criterion:  %9.4f\n", criterion));
        sb.append(String.format("     rng:        %9.4f\n", metropolis));
        if (metropolis < criterion) {
            sb.append(String.format("     Accepted!"));
            logger.info(sb.toString());
            numMovesAccepted++;
            propagateInactiveResidues(titratingMultis, true);
            return true;
        } else {
            // Move was denied.
            sb.append(String.format("     Denied."));
            logger.info(sb.toString());
            // Undo the rejected move.
            RotamerLibrary.applyRotamer(residue, origCoordsRotamer);
            return false;
        }
    }
}
Also used : TitrationUtils.inactivateResidue(ffx.potential.extended.TitrationUtils.inactivateResidue) MultiResidue(ffx.potential.bonded.MultiResidue) Residue(ffx.potential.bonded.Residue) ResidueState(ffx.potential.bonded.ResidueState) AminoAcid3(ffx.potential.bonded.ResidueEnumerations.AminoAcid3) Rotamer(ffx.potential.bonded.Rotamer) Atom(ffx.potential.bonded.Atom)

Example 5 with ResidueState

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

Aggregations

Residue (ffx.potential.bonded.Residue)11 ResidueState (ffx.potential.bonded.ResidueState)11 MultiResidue (ffx.potential.bonded.MultiResidue)10 Rotamer (ffx.potential.bonded.Rotamer)8 Atom (ffx.potential.bonded.Atom)6 RotamerLibrary.applyRotamer (ffx.potential.bonded.RotamerLibrary.applyRotamer)6 NACorrectionException (ffx.potential.bonded.NACorrectionException)5 IOException (java.io.IOException)5 Crystal (ffx.crystal.Crystal)3 PDBFilter (ffx.potential.parsers.PDBFilter)3 File (java.io.File)3 ArrayList (java.util.ArrayList)3 SymOp (ffx.crystal.SymOp)2 AminoAcid3 (ffx.potential.bonded.ResidueEnumerations.AminoAcid3)2 TitrationUtils.inactivateResidue (ffx.potential.extended.TitrationUtils.inactivateResidue)2 BufferedWriter (java.io.BufferedWriter)2 FileWriter (java.io.FileWriter)2 ParallelTeam (edu.rit.pj.ParallelTeam)1 ROLS (ffx.potential.bonded.ROLS)1 Torsion (ffx.potential.bonded.Torsion)1