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

Example 11 with ResidueState

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

the class RosenbluthOBMC method mcStep.

/**
 * Does all the work for a move. Moveset is a continuous 360 degree spin of
 * the chi[0] torsion. U_or in Frenkel's notation (uDep here) is the
 * associated torsion energy. Evaluation criterion: P_accept = Min( 1,
 * (Wn/Wo)*exp(-beta(U[n]-U[o]) )
 */
private boolean mcStep() {
    numMovesProposed++;
    boolean accepted;
    // Select a target residue.
    int which = ThreadLocalRandom.current().nextInt(targets.size());
    Residue target = targets.get(which);
    ResidueState origState = target.storeState();
    List<ROLS> chiList = target.getTorsionList();
    Torsion chi0 = getChiZeroTorsion(target);
    writeSnapshot("orig");
    /* Create old and new trial sets, calculate Wn and Wo, and choose a move bn.
            When doing strictly chi[0] moves, Frenkel/Smit's 'old' and 'new' configurations
            are the same state. The distinction is made here only to aid in future generalization.
         */
    List<MCMove> oldTrialSet = createTrialSet(target, origState, trialSetSize - 1);
    List<MCMove> newTrialSet = createTrialSet(target, origState, trialSetSize);
    report = new StringBuilder();
    report.append(String.format(" Rosenbluth Rotamer MC Move: %4d\n", numMovesProposed));
    report.append(String.format("    residue:   %s\n", target.toString()));
    report.append(String.format("    chi0:      %s\n", chi0.toString()));
    MCMove proposal = calculateRosenbluthFactors(target, chi0, origState, oldTrialSet, origState, newTrialSet);
    /* Calculate the independent portion of the total old-conf energy.
            Then apply the move and calculate the independent total new-conf energy.
         */
    setState(target, origState);
    writeSnapshot("uIndO");
    double uIndO = getTotalEnergy() - getTorsionEnergy(chi0);
    proposal.move();
    writeSnapshot("uIndN");
    double uIndN = getTotalEnergy() - getTorsionEnergy(chi0);
    // Apply acceptance criterion.
    double temperature = thermostat.getCurrentTemperature();
    double beta = 1.0 / (BOLTZMANN * temperature);
    double dInd = uIndN - uIndO;
    double dIndE = FastMath.exp(-beta * dInd);
    double criterion = (Wn / Wo) * FastMath.exp(-beta * (uIndN - uIndO));
    double metropolis = Math.min(1, criterion);
    double rng = ThreadLocalRandom.current().nextDouble();
    report.append(String.format("    theta:     %3.2f\n", ((RosenbluthChi0Move) proposal).theta));
    report.append(String.format("    criterion: %1.4f\n", criterion));
    report.append(String.format("       Wn/Wo:     %.2f\n", Wn / Wo));
    report.append(String.format("       uIndN,O:  %7.2f\t%7.2f\n", uIndN, uIndO));
    report.append(String.format("       dInd(E):  %7.2f\t%7.2f\n", dInd, dIndE));
    report.append(String.format("    rng:       %1.4f\n", rng));
    if (rng < metropolis) {
        numMovesAccepted++;
        report.append(String.format(" Accepted.\n"));
        accepted = true;
    } else {
        proposal.revertMove();
        report.append(String.format(" Denied.\n"));
        accepted = false;
    }
    logger.info(report.toString());
    // Cleanup.
    Wn = 0.0;
    Wo = 0.0;
    return accepted;
}