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

Example 16 with Rotamer

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

Example 17 with Rotamer

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

the class RotamerOptimization method reconcileNARotamersWithPriorResidues.

/**
 * For NA residues inside some optimization window, prune any rotamers which
 * would be incompatible with the established rotamers of upstream NA
 * residues. Could in theory be done by self energies, but every rotamer
 * which can be eliminated without calculating a self energy makes the
 * optimization significantly faster.
 *
 * @param residues Residues to check for incompatible rotamers.
 * @return Number of rotamers eliminated for each Residue.
 */
private int[] reconcileNARotamersWithPriorResidues(Residue[] residues) {
    int[] numEliminatedRotamers = new int[residues.length];
    for (int i = 0; i < residues.length; i++) {
        Residue residuei = residues[i];
        Rotamer[] rotamers = residuei.getRotamers(library);
        if (rotamers == null || residuei.getResidueType() != NA) {
            continue;
        }
        Residue prevResidue = residuei.getPreviousResidue();
        if (prevResidue == null || prevResidue.getResidueType() != NA) {
            continue;
        }
        boolean isInList = false;
        for (Residue residue : residues) {
            if (prevResidue.equals(residue)) {
                isInList = true;
                break;
            }
        }
        if (isInList) {
            continue;
        }
        double prevDelta = RotamerLibrary.measureDelta(prevResidue);
        // If between 50 and 110, assume a North pucker.
        if (RotamerLibrary.checkPucker(prevDelta) == 1) {
            for (int j = 0; j < rotamers.length; j++) {
                if (RotamerLibrary.checkPucker(rotamers[j].chi1) != 1) {
                    if (print) {
                        logIfMaster(format(" Rotamer %d of residue %s eliminated " + "for incompatibility with the sugar pucker of previous " + "residue %s outside the window.", j, residuei.toFormattedString(false, true), prevResidue.toString()));
                    }
                    eliminateRotamer(residues, i, j, print);
                    numEliminatedRotamers[i]++;
                }
            }
        } else {
            for (int j = 0; j < rotamers.length; j++) {
                if (RotamerLibrary.checkPucker(rotamers[j].chi1) != 2) {
                    if (print) {
                        logIfMaster(format(" Rotamer %d of residue %s eliminated " + "for incompatibility with the sugar pucker of previous " + "residue %s outside the window.", j, residuei.toFormattedString(false, true), prevResidue.toString()));
                    }
                    eliminateRotamer(residues, i, j, print);
                    numEliminatedRotamers[i]++;
                }
            }
        }
    // TODO: Implement support for the DNA C3'-exo pucker.
    }
    return numEliminatedRotamers;
}

Example 18 with Rotamer

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

the class RotamerOptimization method prunePairClashes.

/**
 * Prunes rotamer ri of residue i if all ri-j pair energies are worse than
 * the best i-j pair by some threshold value; additionally prunes ri-rj
 * pairs if they exceed the best i-j pair by a greater threshold value;
 * additionally performs this in reverse (searches over j-i).
 *
 * @param residues Residues whose rotamers are to be pruned.
 */
private void prunePairClashes(Residue[] residues) {
    if (!prunePairClashes) {
        return;
    }
    int nResidues = residues.length;
    for (int i = 0; i < nResidues - 1; i++) {
        Residue resi = residues[i];
        Rotamer[] rotsi = resi.getRotamers(library);
        int lenri = rotsi.length;
        for (int j = i + 1; j < nResidues; j++) {
            Residue resj = residues[j];
            Rotamer[] rotsj = resj.getRotamers(library);
            int lenrj = rotsj.length;
            double minPair = Double.MAX_VALUE;
            int minRI = -1;
            int minRJ = -1;
            for (int ri = 0; ri < lenri; ri++) {
                if (check(i, ri)) {
                    continue;
                }
                for (int rj = 0; rj < lenrj; rj++) {
                    if (check(j, rj) || check(i, ri, j, rj)) {
                        continue;
                    }
                    double pairEnergy = get2Body(i, ri, j, rj) + getSelf(i, ri) + getSelf(j, rj);
                    assert Double.isFinite(pairEnergy);
                    if (pairEnergy < minPair) {
                        minPair = pairEnergy;
                        minRI = ri;
                        minRJ = rj;
                    }
                }
            }
            // Otherwise, i and j are not on a well-packed backbone.
            assert (minRI >= 0 && minRJ >= 0);
            // Calculate all the modifiers to the pair clash elimination threshold.
            double threshold = pairClashThreshold;
            if (resi instanceof MultiResidue) {
                threshold += multiResPairClashAddn;
            }
            if (resj instanceof MultiResidue) {
                threshold += multiResPairClashAddn;
            }
            int numNARes = (resi.getResidueType() == NA ? 1 : 0) + (resj.getResidueType() == NA ? 1 : 0);
            switch(numNARes) {
                case 0:
                    break;
                case 1:
                    threshold *= pairHalfPruningFactor;
                    break;
                case 2:
                    threshold *= pruningFactor;
                    break;
                default:
                    throw new ArithmeticException(" RotamerOptimization.prunePairClashes() has somehow " + "found less than zero or more than two nucleic acid residues in a pair of" + " residues. This result should be impossible.");
            }
            double toEliminate = threshold + minPair;
            for (int ri = 0; ri < lenri; ri++) {
                if (check(i, ri)) {
                    continue;
                }
                for (int rj = 0; rj < lenrj; rj++) {
                    if (check(j, rj) || check(i, ri, j, rj)) {
                        continue;
                    }
                    double pairEnergy = get2Body(i, ri, j, rj) + getSelf(i, ri) + getSelf(j, rj);
                    assert Double.isFinite(pairEnergy);
                    if (pairEnergy > toEliminate) {
                        logIfMaster(format(" Pruning pair %s-%d %s-%d by %s-%d %s-%d; energy %s > " + "%s + %s", resi.toFormattedString(false, true), ri, resj.toFormattedString(false, true), rj, resi.toFormattedString(false, true), minRI, resj.toFormattedString(false, true), minRJ, formatEnergy(pairEnergy), formatEnergy(threshold), formatEnergy(minPair)));
                    }
                }
            }
            pairsToSingleElimination(residues, i, j);
        }
    }
}

Example 19 with Rotamer

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

the class PhMD method mcUpdate.

/**
 * The primary driver. Called by the MD engine at each dynamics step.
 */
@Override
public boolean mcUpdate(double temperature) {
    startTime = System.nanoTime();
    if (thermostat.getCurrentTemperature() > config.meltdownTemperature) {
        meltdown();
    }
    if (thermostat.getCurrentTemperature() > config.warningTemperature) {
        Atom[] atoms = mola.getAtomArray();
        logger.info(format(" System heating! Dumping atomic velocities for %d D.o.F.:", ffe.getNumberOfVariables()));
        double[] velocity = new double[3];
        for (Atom atom : atoms) {
            atom.getVelocity(velocity);
            logger.info(format(" %s: %s", atom.describe(Atom.Descriptions.Trim), Arrays.toString(velocity)));
        }
    }
    esvSystem.setTemperature(temperature);
    propagateInactiveResidues(titratingMultis);
    stepCount++;
    // Decide on the type of step to be taken.
    StepType stepType;
    if (stepCount % mcStepFrequency == 0 && stepCount % rotamerStepFrequency == 0) {
        stepType = StepType.COMBO;
    } else if (stepCount % mcStepFrequency == 0) {
        stepType = StepType.TITRATE;
    } else if (stepCount % rotamerStepFrequency == 0) {
        stepType = StepType.ROTAMER;
    } else {
        // Not yet time for an MC step, return to MD.
        if (config.logTimings) {
            long took = System.nanoTime() - startTime;
            logger.info(String.format(" CpHMD propagation time: %.6f", took * NS_TO_SEC));
        }
        return false;
    }
    logger.info(format("TitratingMultis: %d", titratingMultis.size()));
    // Randomly choose a target titratable residue to attempt protonation switch.
    int random = (config.titrateTermini) ? rng.nextInt(titratingMultis.size() + titratingTermini.size()) : rng.nextInt(titratingMultis.size());
    if (random >= titratingMultis.size()) {
        Residue target = titratingTermini.get(random - titratingMultis.size());
        boolean accepted = tryTerminusTitration((MultiTerminus) target);
        snapshotIndex++;
        if (accepted) {
            molDyn.reInit();
            previousTarget = target;
        }
        return accepted;
    }
    MultiResidue targetMulti = titratingMultis.get(random);
    // Check whether rotamer moves are possible for the selected residue.
    Residue targetMultiActive = targetMulti.getActive();
    Rotamer[] targetMultiRotamers = targetMultiActive.getRotamers(library);
    if (targetMultiRotamers == null || targetMultiRotamers.length <= 1) {
        if (stepType == StepType.ROTAMER) {
            return false;
        } else if (stepType == StepType.COMBO) {
            stepType = StepType.TITRATE;
        }
    }
    // Perform the MC move.
    boolean accepted;
    switch(stepType) {
        case TITRATE:
            accepted = tryTitrationStep(targetMulti);
            break;
        case ROTAMER:
            accepted = (config.useConformationalBias) ? tryCBMCStep(targetMulti) : tryRotamerStep(targetMulti);
            break;
        case COMBO:
            accepted = (config.useConformationalBias) ? tryCBMCStep(targetMulti) || tryTitrationStep(targetMulti) : tryComboStep(targetMulti);
            break;
        default:
            accepted = false;
            throw new IllegalStateException();
    }
    snapshotIndex++;
    if (accepted) {
        previousTarget = targetMulti;
    }
    if (config.logTimings) {
        long took = System.nanoTime() - startTime;
        logger.info(String.format(" CpHMD step time:        %.6f", took * NS_TO_SEC));
    }
    return accepted;
}

Example 20 with Rotamer

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