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

Example 81 with Residue

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

the class RotamerOptimization method applyEliminationCriteria.

private void applyEliminationCriteria(Residue[] residues) {
    // allocateEliminationMemory is now called for all algorithms in rotamerEnergies method.
    // allocateEliminationMemory(residues);
    /*
         * Must pin 5' ends of nucleic acids which are attached to nucleic acids
         * outside the window, to those prior residues' sugar puckers.  Then,
         * if a correction threshold is set, eliminate rotamers with excessive
         * correction vectors (up to a maximum defined by minNumberAcceptedNARotamers).
         */
    boolean containsNA = false;
    if (pruneClashes) {
        for (Residue residue : residues) {
            if (residue.getResidueType() == NA) {
                containsNA = true;
                break;
            }
        }
    }
    if (containsNA && pruneClashes) {
        logIfMaster(" Eliminating nucleic acid rotamers that conflict at their 5' end with residues outside the optimization range.");
        int[] numEliminatedRotamers = reconcileNARotamersWithPriorResidues(residues);
        // if the input structure is no good.
        if (verboseEnergies) {
            try {
                logIfMaster(format("\n Beginning Energy %s", formatEnergy(currentEnergy(residues))));
            } catch (ArithmeticException ex) {
                logger.severe(String.format(" Exception %s in calculating beginning energy; FFX shutting down.", ex.toString()));
            }
        }
        eliminateNABackboneRotamers(residues, numEliminatedRotamers);
    } else if (verboseEnergies) {
        try {
            logIfMaster(format("\n Beginning Energy %s", formatEnergy(currentEnergy(residues))));
        } catch (ArithmeticException ex) {
            logger.severe(String.format(" Exception %s in calculating beginning energy; FFX shutting down.", ex.toString()));
        }
    }
    rotamerEnergies(residues);
    if (testing) {
        int nres = residues.length;
        onlyPrunedSingles = new boolean[nres][];
        onlyPrunedPairs = new boolean[nres][][][];
        for (int i = 0; i < nres; i++) {
            Residue residuei = residues[i];
            Rotamer[] rotamersi = residuei.getRotamers(library);
            // Length rotamers i
            int lenri = rotamersi.length;
            onlyPrunedSingles[i] = new boolean[lenri];
            onlyPrunedSingles[i] = Arrays.copyOf(eliminatedSingles[i], eliminatedSingles[i].length);
            onlyPrunedPairs[i] = new boolean[lenri][][];
            // Loop over the set of rotamers for residue i.
            for (int ri = 0; ri < lenri; ri++) {
                onlyPrunedPairs[i][ri] = new boolean[nres][];
                for (int j = i + 1; j < nres; j++) {
                    Residue residuej = residues[j];
                    Rotamer[] rotamersj = residuej.getRotamers(library);
                    int lenrj = rotamersj.length;
                    onlyPrunedPairs[i][ri][j] = new boolean[lenrj];
                    onlyPrunedPairs[i][ri][j] = Arrays.copyOf(eliminatedPairs[i][ri][j], eliminatedPairs[i][ri][j].length);
                }
            }
        }
    }
    if (testSelfEnergyEliminations) {
        testSelfEnergyElimination(residues);
    } else if (testPairEnergyEliminations > -1) {
        testPairEnergyElimination(residues, testPairEnergyEliminations);
    } else if (testTripleEnergyEliminations1 > -1 && testTripleEnergyEliminations2 > -1) {
        testTripleEnergyElimination(residues, testTripleEnergyEliminations1, testTripleEnergyEliminations2);
    }
    // testSelfEnergyElimination(residues);
    // testPairEnergyElimination(residues, 19);
    // Beginning energy
    int[] currentRotamers = new int[residues.length];
    if (pruneClashes) {
        validateDEE(residues);
    }
    int i = 0;
    boolean pairEliminated;
    do {
        pairEliminated = false;
        if (useGoldstein) {
            if (selfEliminationOn) {
                i++;
                logIfMaster(format("\n Iteration %d: Applying Single Goldstein DEE conditions ", i));
                // While there are eliminated rotamers, repeatedly apply single rotamer elimination.
                while (goldsteinDriver(residues)) {
                    i++;
                    logIfMaster(this.toString());
                    logIfMaster(format("\n Iteration %d: Applying Single Rotamer Goldstein DEE conditions ", i));
                }
            }
            if (pairEliminationOn) {
                i++;
                logIfMaster(format("\n Iteration %d: Applying Rotamer Pair Goldstein DEE conditions ", i));
                // While there are eliminated rotamer pairs, repeatedly apply rotamer pair elimination.
                while (goldsteinPairDriver(residues)) {
                    pairEliminated = true;
                    i++;
                    logIfMaster(this.toString());
                    logIfMaster(format("\n Iteration %d: Applying Rotamer Pair Goldstein DEE conditions ", i));
                }
            }
        } else {
            if (selfEliminationOn) {
                i++;
                logIfMaster(format("\n Iteration %d: Applying Single DEE conditions ", i));
                // While there are eliminated rotamers, repeatedly apply single rotamer elimination.
                while (deeRotamerElimination(residues)) {
                    i++;
                    logIfMaster(toString());
                    logIfMaster(format("\n Iteration %d: Applying Single Rotamer DEE conditions ", i));
                }
            }
            if (pairEliminationOn) {
                i++;
                logIfMaster(format("\n Iteration %d: Applying Rotamer Pair DEE conditions ", i));
                // While there are eliminated rotamer pairs, repeatedly apply rotamer pair elimination.
                while (deeRotamerPairElimination(residues)) {
                    pairEliminated = true;
                    i++;
                    logIfMaster(toString());
                    logIfMaster(format("\n Iteration %d: Applying Rotamer Pair DEE conditions ", i));
                }
            }
        }
        validateDEE(residues);
        logIfMaster(toString());
    } while (pairEliminated);
    logIfMaster(" Self-consistent DEE rotamer elimination achieved.\n");
}

Example 82 with Residue

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

the class RotamerOptimization method testPairEnergyElimination.

/**
 * Test the elimination criteria by setting self and 3-body interactions to
 * zero.
 */
public void testPairEnergyElimination(Residue[] residues, int resID) {
    int nRes = residues.length;
    if (resID >= nRes) {
        return;
    }
    for (int i = 0; i < nRes; i++) {
        Residue resI = residues[i];
        Rotamer[] rotI = resI.getRotamers(library);
        int nI = rotI.length;
        for (int ri = 0; ri < nI; ri++) {
            try {
                selfEnergy[i][ri] = 0.0;
            } catch (Exception e) {
            // catch NPE.
            }
            for (int j = i + 1; j < nRes; j++) {
                Residue resJ = residues[j];
                Rotamer[] rotJ = resJ.getRotamers(library);
                int nJ = rotJ.length;
                for (int rj = 0; rj < nJ; rj++) {
                    if (i != resID && j != resID) {
                        try {
                            twoBodyEnergy[i][ri][j][rj] = 0.0;
                        } catch (Exception e) {
                        // catch NPE.
                        }
                    }
                    if (threeBodyTerm) {
                        for (int k = j + 1; k < nRes; k++) {
                            Residue resK = residues[k];
                            Rotamer[] rotK = resK.getRotamers(library);
                            int nK = rotK.length;
                            for (int rk = 0; rk < nK; rk++) {
                                try {
                                    threeBodyEnergy[i][ri][j][rj][k][rk] = 0.0;
                                } catch (Exception e) {
                                // catch NPE.
                                }
                            }
                        }
                    }
                }
            }
        }
    }
}

Example 83 with Residue

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

Example 84 with Residue

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

the class RotamerOptimization method deeRotamerElimination.

/**
 * Elimination of rotamers via the original Dead End Elimination algorithm.
 *
 * @param residues Array of residues under consideration.
 * @return True if any rotamers were eliminated.
 */
private boolean deeRotamerElimination(Residue[] residues) {
    int nres = residues.length;
    // A flag to indicate if any more rotamers or rotamer pairs were eliminated.
    boolean eliminated = false;
    // Loop over residues.
    double[] minMax = new double[2];
    double[] minEnergySingles = null;
    double[] maxEnergySingles = null;
    for (int i = 0; i < nres; i++) {
        Residue residuei = residues[i];
        Rotamer[] rotamersi = residuei.getRotamers(library);
        int lenri = rotamersi.length;
        if (minEnergySingles == null || minEnergySingles.length < lenri) {
            minEnergySingles = new double[lenri];
            maxEnergySingles = new double[lenri];
        }
        // Loop over the set of rotamers for residue i.
        for (int ri = 0; ri < lenri; ri++) {
            // Check for an eliminated single.
            if (check(i, ri)) {
                continue;
            }
            // Start the min/max summation with the self-energy.
            minEnergySingles[ri] = getSelf(i, ri);
            maxEnergySingles[ri] = minEnergySingles[ri];
            for (int j = 0; j < nres; j++) {
                if (j == i) {
                    continue;
                }
                if (minMaxPairEnergy(residues, minMax, i, ri, j)) {
                    if (Double.isFinite(minMax[0]) && Double.isFinite(minEnergySingles[ri])) {
                        minEnergySingles[ri] += minMax[0];
                    } else {
                        // There is a possibility that ri conflicts with every possible rotamer of some residue j.
                        // In that event, its minimum energy is set NaN and should be easy to eliminate.
                        minEnergySingles[ri] = Double.NaN;
                    }
                    if (Double.isFinite(minMax[0]) && Double.isFinite(maxEnergySingles[ri])) {
                        maxEnergySingles[ri] += minMax[1];
                    } else {
                        // In this branch, ri conflicts with some j,rj and cannot be practically used for elimination.
                        maxEnergySingles[ri] = Double.NaN;
                    }
                } else {
                    Residue residuej = residues[j];
                    logger.info(format(" Inconsistent Pair: %8s %2d, %8s.", residuei.toFormattedString(false, true), ri, residuej.toFormattedString(false, true)));
                // eliminateRotamer(residues, i, ri, print);
                }
            }
        }
        /**
         * Apply the singles elimination criteria to rotamers of residue i
         * by determining the most favorable maximum energy.
         */
        double eliminationEnergy = Double.MAX_VALUE;
        int eliminatingRotamer = 0;
        for (int ri = 0; ri < lenri; ri++) {
            if (check(i, ri)) {
                continue;
            }
            if (Double.isFinite(maxEnergySingles[ri]) && maxEnergySingles[ri] < eliminationEnergy) {
                eliminationEnergy = maxEnergySingles[ri];
                eliminatingRotamer = ri;
            }
        }
        if (eliminationEnergy == Double.MAX_VALUE) {
            // This branch is taken if every ri conflicts with at least one j,rj. In that case, nothing can be eliminated yet!
            logIfMaster(" Could not eliminate any i,ri because eliminationEnergy was never set!", Level.FINE);
        } else {
            /**
             * Eliminate rotamers whose minimum energy is greater than the
             * worst case for another rotamer.
             */
            for (int ri = 0; ri < lenri; ri++) {
                if (check(i, ri)) {
                    continue;
                }
                // If i,ri has a clash with all of phase space, it can be eliminated by something that doesn't clash with all phase space.
                if (!Double.isFinite(minEnergySingles[ri])) {
                    if (eliminateRotamer(residues, i, ri, print)) {
                        logIfMaster(format("  Rotamer elimination of (%8s,%2d) that always clashes.", residuei.toFormattedString(false, true), ri));
                        eliminated = true;
                    }
                }
                // Otherwise, can eliminate if its best possible energy is still worse than something else's worst possible energy.
                if (minEnergySingles[ri] > eliminationEnergy + ensembleBuffer) {
                    if (eliminateRotamer(residues, i, ri, print)) {
                        logIfMaster(format("  Rotamer elimination of (%8s,%2d) by (%8s,%2d): %12.4f > %6.4f.", residuei.toFormattedString(false, true), ri, residuei.toFormattedString(false, true), eliminatingRotamer, minEnergySingles[ri], eliminationEnergy + ensembleBuffer));
                        eliminated = true;
                    }
                }
            }
        }
    }
    return eliminated;
}

Example 85 with Residue

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

the class PhDiscount method tryContinuousTitrationMove.

/**
 * Run continuous titration MD in implicit solvent as a Monte Carlo move.
 */
private boolean tryContinuousTitrationMove(int titrationDuration, double targetTemperature) {
    long startTime = System.nanoTime();
    if (thermostat.getCurrentTemperature() > config.meltdownTemperature) {
        crashDump(new SystemTemperatureException(thermostat.getCurrentTemperature()));
    }
    propagateInactiveResidues(titratingMultiResidues, true);
    // Save the current state of the molecularAssembly. Specifically,
    // Atom coordinates and MultiResidue states : AssemblyState
    // Position, Velocity, Acceleration, etc    : DynamicsState (via MD::storeState)
    AssemblyState assemblyState = new AssemblyState(mola);
    molDyn.storeState();
    writeSnapshot(".pre-store");
    /* Assign initial titration lambdas. */
    switch(config.seedDistribution) {
        default:
        case FLAT:
            /* All lambdas have equal probability. */
            for (int i = 0; i < numESVs; i++) {
                esvSystem.setLambda(i, rng.nextDouble());
            }
            break;
        case BETA:
            /* Draw from a mathematical distribution. */
            throw new UnsupportedOperationException();
        case BOLTZMANN:
            /* Draw from a Boltzmann distribution (via Rosenbluth sets). */
            throw new UnsupportedOperationException();
        case DIRAC_CURRENT:
            /* Assign only zero or unity representing current state. */
            for (int i = 0; i < numESVs; i++) {
                Residue active = titratingMultiResidues.get(i).getActive();
                double current = (active.getAminoAcid3() == Titration.lookup(active).protForm) ? 1.0 : 0.0;
                esvSystem.setLambda(i, current);
            }
            break;
        case DIRAC_POINTFIVE:
            /* Assign half-occupancy to all titratable protons. */
            for (int i = 0; i < numESVs; i++) {
                esvSystem.setLambda(i, 0.5);
            }
            break;
    }
    /*
         * (1) Take pre-move energy.
         * (2) Hook the ExtendedSystem up to MolecularDynamics.
         * (3) Launch dynamics for nSteps = Monte-Carlo Frequency.
         * (4) Note that no callbacks to mcUpdate() occur during this period.
         * (5) Floor/ceil to discretize hydrogen occupancies.
         * (6) Take post-move energy and test on the combined Metropolis criterion.
         */
    final double Uo = currentTotalEnergy();
    // TODO: Need to ensure that we fully protonate the protein before entering continuous-protonation space.
    ffe.attachExtendedSystem(esvSystem);
    molDyn.attachExtendedSystem(esvSystem, 10);
    final double Uo_prime = currentTotalEnergy();
    logger.info(format(" %-40s %-s", "Trying continuous titration move.", format("Uo,Uo': %16.8f, %16.8f", Uo, Uo_prime)));
    molDyn.redynamic(titrationDuration, targetTemperature);
    final double Un_prime = currentTotalEnergy();
    for (int i = 0; i < esvSystem.size(); i++) {
        esvSystem.setLambda(i, Math.rint(esvSystem.getLambda(i)));
    }
    molDyn.detachExtendedSystem();
    ffe.detachExtendedSystem();
    final double Un = currentTotalEnergy();
    logger.info(format(" %-40s %-30s", "Move finished; detaching esvSystem.", format("Un',Un: %16.8f, %16.8f", Un_prime, Un)));
    /* Calculate acceptance probability from detailed balance equation. */
    final double beta = 1 / (ExtConstants.Boltzmann * thermostat.getCurrentTemperature());
    final double dgDiscrete = Un - Uo;
    final double dgContinuous = Un_prime - Uo_prime;
    final double crit = FastMath.exp(-beta * (dgDiscrete - dgContinuous));
    final double rand = rng.nextDouble();
    logger.info(format("   Un,Un',Uo',Uo:   %10.5f %10.5f %10.5f %10.5f", Un, Un_prime, Uo_prime, Uo));
    logger.info(format("   Crit,Rng:        %10.5f %10.5f", crit, rand));
    long took = System.nanoTime() - startTime;
    if (rand <= crit) {
        logger.info(format(" %-40s %-s", "Monte-Carlo accepted.", format("Wallclock: %8.3f sec", took * ns2sec)));
        writeSnapshot(".post-acpt");
        return true;
    } else {
        logger.info(format(" %-40s %-s", "Monte-Carlo denied; reverting state.", format("Wallclock: %8.3f sec", took * ns2sec)));
        writeSnapshot(".post-deny");
        assemblyState.revertState();
        ffe.reInit();
        try {
            molDyn.revertState();
        } catch (Exception ex) {
            Logger.getLogger(PhDiscount.class.getName()).log(Level.SEVERE, null, ex);
        }
        writeSnapshot(".post-rvrt");
        return false;
    }
}