Examples with AssemblyState ffx.potential.AssemblyState used on opensource projects

Example 1 with AssemblyState

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