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

Example 21 with Rotamer

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

Example 22 with Rotamer

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

the class RosenbluthChiAllMove method torsionSampler.

private void torsionSampler(List<Torsion> allTors) {
    // Collects data for plot of uTors vs. theta.
    // This is for finding the appropriate offset to add to each uTors such that the minimum is zero.
    // double origChi[] = RotamerLibrary.measureRotamer(target, false);
    updateAll();
    double[] origChi = measureLysine(target, false);
    StringBuilder sb = new StringBuilder();
    sb.append(String.format(" Torsion Sampling for %s\n", target.getName()));
    sb.append(String.format("   origChi:"));
    for (int k = 0; k < origChi.length; k++) {
        sb.append(String.format(" %6.2f", origChi[k]));
    }
    sb.append("\n");
    for (int k = 0; k < origChi.length; k++) {
        Torsion tors = allTors.get(k);
        AtomType type1 = tors.getAtomArray()[0].getAtomType();
        AtomType type2 = tors.getAtomArray()[1].getAtomType();
        AtomType type3 = tors.getAtomArray()[2].getAtomType();
        AtomType type4 = tors.getAtomArray()[3].getAtomType();
        sb.append(String.format("   %d:    \"(%3d %3d %3s)  (%3d %3d %3s)  (%3d %3d %3s)  (%3d %3d %3s)\"\n", k, type1.type, type1.atomClass, type1.name, type2.type, type2.atomClass, type2.name, type3.type, type3.atomClass, type3.name, type4.type, type4.atomClass, type4.name));
    }
    logger.info(sb.toString());
    sb = new StringBuilder();
    String type = "combinatorial";
    if (type.equals("combinatorial")) {
        sb.append(String.format(" Resi chi0 chi1 chi2 chi3 List<uTors> uTorsSum uTotal\n"));
        logger.info(sb.toString());
        sb = new StringBuilder();
        boolean doChi0 = false, doChi1 = false;
        boolean doChi2 = true, doChi3 = true;
        double increment = 1.0;
        for (double chi0 = -180.0; chi0 < +180.0; chi0 += increment) {
            for (double chi1 = -180.0; chi1 <= +180.0; chi1 += increment) {
                for (double chi2 = -180.0; chi2 <= +180.0; chi2 += increment) {
                    for (double chi3 = -180.0; chi3 <= +180.0; chi3 += increment) {
                        sb.append(String.format(" %3s", target.getName()));
                        double[] newChi = new double[4];
                        if (doChi0) {
                            newChi[0] = chi0;
                        } else {
                            newChi[0] = origChi[0];
                        }
                        if (doChi1) {
                            newChi[1] = chi1;
                        } else {
                            newChi[1] = origChi[1];
                        }
                        if (doChi2) {
                            newChi[2] = chi2;
                        } else {
                            newChi[2] = origChi[2];
                        }
                        if (doChi3) {
                            newChi[3] = chi3;
                        } else {
                            newChi[3] = origChi[3];
                        }
                        Torsion chi0Tors = allTors.get(0);
                        Torsion chi1Tors = allTors.get(1);
                        Torsion chi2Tors = allTors.get(2);
                        Torsion chi3Tors = allTors.get(3);
                        Rotamer newState = createRotamer(target, newChi);
                        RotamerLibrary.applyRotamer(target, newState);
                        for (int wut = 0; wut < origChi.length; wut++) {
                            sb.append(String.format(" %3.0f", newChi[wut]));
                        }
                        writeSnapshot(String.format("%.0f", chi1), false);
                        double uTorsSum = 0.0;
                        for (int k = 0; k < allTors.size(); k++) {
                            double uTors = allTors.get(k).energy(false);
                            sb.append(String.format(" %5.2f", uTors));
                            uTorsSum += uTors;
                        }
                        sb.append(String.format(" %5.2f", uTorsSum));
                        double totalE = 1000.0;
                        try {
                            totalE = totalEnergy();
                        } catch (Exception ex) {
                        }
                        if (totalE > 1000.0) {
                            totalE = 1000.0;
                        }
                        sb.append(String.format(" %10.4f", totalE));
                        sb.append(String.format("\n"));
                        logger.info(sb.toString());
                        sb = new StringBuilder();
                        if (!doChi3) {
                            break;
                        }
                    }
                    if (!doChi2) {
                        break;
                    }
                }
                if (!doChi1) {
                    break;
                }
            }
            if (!doChi0) {
                break;
            }
        }
    } else {
        sb.append(String.format(" Resi Theta ( chi ) List<uTors,uTotal> \n"));
        logger.info(sb.toString());
        sb = new StringBuilder();
        for (double testTheta = -180.0; testTheta <= +180.0; testTheta += 0.01) {
            sb.append(String.format(" %3s %6.2f", target.getName(), testTheta));
            for (int k = 0; k < origChi.length; k++) {
                double[] newChi = new double[origChi.length];
                System.arraycopy(origChi, 0, newChi, 0, origChi.length);
                Torsion tors = allTors.get(k);
                newChi[k] = testTheta;
                Rotamer newState = createRotamer(target, newChi);
                RotamerLibrary.applyRotamer(target, newState);
                sb.append(" (");
                for (int wut = 0; wut < origChi.length; wut++) {
                    sb.append(String.format(" %6.2f", newChi[wut]));
                }
                sb.append(" )");
                writeSnapshot(String.format("%.0f", testTheta), false);
                double uTors = allTors.get(k).energy(false);
                double totalE = 0.0;
                try {
                    totalE = totalEnergy();
                } catch (Exception ex) {
                }
                sb.append(String.format(" %9.6f %9.6f", uTors, totalE));
            }
            sb.append(String.format("\n"));
            logger.info(sb.toString());
            sb = new StringBuilder();
        }
    }
    if (System.getProperty("cbmc-torsionSampler") != null) {
        try {
            File output = new File(System.getProperty("cbmc-torsionSampler"));
            BufferedWriter bw = new BufferedWriter(new FileWriter(output));
            bw.write(sb.toString());
            bw.close();
        } catch (IOException ex) {
        }
    }
    System.exit(0);
}

Example 23 with Rotamer

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

the class RosenbluthChiAllMove method cheapTorsionSet_diffs.

private TrialSet cheapTorsionSet_diffs(int setSize, String snapSuffix) {
    double origFastEnergy = fastEnergy();
    double origSlowEnergy = slowEnergy();
    double origTotalEnergy = totalEnergy();
    if (Math.abs(origTotalEnergy) > Math.abs(origFastEnergy + origSlowEnergy) + tolerance) {
        logger.warning(String.format("WAIT WHAT.  total != fast+slow : %9.4f  %9.4f  %9.4f", origTotalEnergy, origFastEnergy, origSlowEnergy));
    }
    report.append(String.format("    CheapTrialSet_Diffs (uDep  uExt  uExtBolt  running)\n"));
    TrialSet trialSet = new TrialSet(setSize);
    double[] origChi = RotamerLibrary.measureRotamer(target, false);
    double loggerW = 0.0;
    int i = 0;
    while (i < setSize) {
        double[] newChi = new double[origChi.length];
        System.arraycopy(origChi, 0, newChi, 0, origChi.length);
        Rotamer newState = null;
        for (int k = 0; k < origChi.length; k++) {
            double crit = 1.00, rng = 1.00, theta, uFastThis = 0.0, uFast;
            do {
                theta = rand.nextDouble(360.0) - 180;
                newChi[k] = theta;
                newState = createRotamer(target, newChi);
                RotamerLibrary.applyRotamer(target, newState);
                uFastThis = fastEnergy();
                uFast = uFastThis - origFastEnergy;
                if (uFast <= 0) {
                    break;
                }
                crit = FastMath.exp(-beta * uFast);
                rng = rand.nextDouble();
            } while (rng >= crit);
            report.append(String.format(" Accept-reject at movenum %d pos %d chi,fast-orig=eng,crit (rng):   %5.2f (%5.2f - %5.2f = %5.2f) %5.2f  (%3.2f)\n", moveNumber, k, theta, uFastThis, origFastEnergy, uFast, crit, rng));
        }
        // this cheap version does all thetas at once
        trialSet.theta[i] = 0.0;
        trialSet.rotamer[i] = newState;
        trialSet.uDep[i] = fastEnergy() - origFastEnergy;
        trialSet.uExt[i] = slowEnergy() - origSlowEnergy;
        loggerW += FastMath.exp(-beta * trialSet.uExt[i]);
        if (i < 4 || i > setSize - 2) {
            report.append(String.format("       %3s %d:  %9.5g  %9.5g  %9.5g  %9.5g\n", snapSuffix, i + 1, trialSet.uDep[i], trialSet.uExt[i], FastMath.exp(-beta * trialSet.uExt[i]), loggerW));
        } else if (i == 4) {
            report.append(String.format("       ...\n"));
        }
        writeSnapshot(snapSuffix, true);
        i++;
    }
    target.revertState(origState);
    updateAll();
    return trialSet;
}

Example 24 with Rotamer

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

the class RosenbluthChiAllMove method cheapTorsionSet_indiv.

/**
 * This version of the cheap method draws each INDIVIDUAL chi from its OWN
 * Boltzmann. Each member of the test set is still a full set of chis.
 */
private TrialSet cheapTorsionSet_indiv(List<Torsion> allTors, int setSize, String snapSuffix) {
    report.append(String.format("    CheapTrialSet_Indiv (uDep  uExt  uExtBolt  running)\n"));
    TrialSet trialSet = new TrialSet(setSize);
    double[] origChi = RotamerLibrary.measureRotamer(target, false);
    double loggerW = 0.0;
    int i = 0;
    while (i < setSize) {
        double[] newChi = new double[origChi.length];
        System.arraycopy(origChi, 0, newChi, 0, origChi.length);
        Rotamer newState = null;
        for (int k = 0; k < origChi.length; k++) {
            double crit, rng, theta, uTors, offset = 0.0;
            do {
                theta = rand.nextDouble(360.0) - 180;
                newChi[k] = theta;
                newState = createRotamer(target, newChi);
                RotamerLibrary.applyRotamer(target, newState);
                uTors = allTors.get(k).energy(false);
                try {
                    offset = TORSION_OFFSET_AMPRO13.valueOf(target.getName() + k).offset;
                } catch (IllegalArgumentException ex) {
                    logger.warning(ex.getMessage());
                }
                uTors += offset;
                crit = FastMath.exp(-beta * uTors);
                rng = rand.nextDouble();
            } while (rng >= crit);
        // report.append(String.format(" Accept-reject at movenum %d pos %d chi,eng,offset,crit:   %5.2f  %5.2f  %5.2f  %5.2f  %5.2f\n",
        // moveNumber, k, theta, uTors, offset, crit, rng));
        }
        double uTorsSum = 0;
        for (Torsion tors : allTors) {
            uTorsSum += tors.energy(false);
        }
        // this cheap version does all thetas at once
        trialSet.theta[i] = 0.0;
        trialSet.rotamer[i] = newState;
        trialSet.uDep[i] = uTorsSum;
        trialSet.uExt[i] = totalEnergy() - uTorsSum;
        loggerW += FastMath.exp(-beta * trialSet.uExt[i]);
        if (i < 4 || i > setSize - 2) {
            report.append(String.format("       %3s %d:  %9.5g  %9.5g  %9.5g  %9.5g\n", snapSuffix, i + 1, trialSet.uDep[i], trialSet.uExt[i], FastMath.exp(-beta * trialSet.uExt[i]), loggerW));
        } else if (i == 4) {
            report.append(String.format("       ...\n"));
        }
        writeSnapshot(snapSuffix, true);
        i++;
    }
    target.revertState(origState);
    updateAll();
    return trialSet;
}

Example 25 with Rotamer

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

the class RosenbluthChiAllMove method expensiveTorsionSet.

/**
 * Uses the accept-reject method (F/S Algorithm46) to draw new chi values
 * for the given torsion.
 */
private TrialSet expensiveTorsionSet(Torsion tors, int chiIndex, int setSize, String snapSuffix) {
    report.append(String.format("    TrialSet for Chi%d\t\t(Theta uDep uExt)\n", chiIndex));
    TrialSet trialSet = new TrialSet(setSize);
    double[] origChi = RotamerLibrary.measureRotamer(target, false);
    int i = 0;
    while (i < setSize) {
        double theta = rand.nextDouble(360.0) - 180;
        double[] newChi = new double[origChi.length];
        System.arraycopy(origChi, 0, newChi, 0, origChi.length);
        newChi[chiIndex] = theta;
        Rotamer newState = createRotamer(target, newChi);
        RotamerLibrary.applyRotamer(target, newState);
        double uTors = tors.energy(false);
        double offset = 0;
        try {
            offset = TORSION_OFFSET_AMPRO13.valueOf(target.getName() + chiIndex).offset;
        } catch (IllegalArgumentException ex) {
            logger.warning(ex.getMessage());
        }
        uTors += offset;
        double criterion = FastMath.exp(-beta * uTors);
        double rng = rand.nextDouble();
        report.append(String.format("    prop: %5.1f %.2g %.2g %.2g %.2g\n", theta, uTors, criterion, rng));
        if (rng < criterion) {
            report.append(String.format("    ^ Accepted!\n"));
            writeSnapshot(snapSuffix, false);
            trialSet.theta[i] = theta;
            trialSet.rotamer[i] = newState;
            trialSet.uDep[i] = uTors;
            // Expensive!
            trialSet.uExt[i] = totalEnergy() - uTors;
            i++;
            writeSnapshot(snapSuffix, true);
            if (i < 4 || i > setSize - 1) {
                report.append(String.format("       %3s %d:      %5.2f\t%5.2f\t%5.2f\n", snapSuffix, i, theta, trialSet.uDep[i - 1], trialSet.uExt[i - 1]));
            } else if (i == 4) {
                report.append(String.format("       ...\n"));
            }
        }
    }
    target.revertState(origState);
    updateAll();
    return trialSet;
}