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

Example 16 with MultiResidue

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

Example 17 with MultiResidue

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

the class PhMD method tryCBMCStep.

private boolean tryCBMCStep(MultiResidue targetMulti) {
    if (CAUTIOUS) {
        throw new UnsupportedOperationException();
    }
    List<Residue> targets = new ArrayList<>();
    targets.add(targetMulti.getActive());
    // cost still scales with this, unfortunately
    int trialSetSize = 5;
    // irrelevant for manual step call
    int mcFrequency = 1;
    boolean writeSnapshots = false;
    System.setProperty("cbmc-type", "CHEAP");
    RosenbluthCBMC cbmc = new RosenbluthCBMC(mola, mola.getPotentialEnergy(), null, targets, mcFrequency, trialSetSize, writeSnapshots);
    boolean accepted = cbmc.cbmcStep();
    if (config.logTimings) {
        long took = System.nanoTime() - startTime;
        logger.info(String.format(" CBMC time: %1.3f", took * NS_TO_SEC));
    }
    return accepted;
}

Example 18 with MultiResidue

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

the class PhMD method tryTitrationStep.

/**
 * Perform a titration MC move.
 *
 * @param targetMulti
 * @return accept/reject
 */
private boolean tryTitrationStep(Residue target) {
    boolean terminus = false;
    MultiResidue targetMulti = null;
    MultiTerminus targetTerm = null;
    if (target instanceof MultiResidue) {
        targetMulti = (MultiResidue) target;
        terminus = false;
    } else if (target instanceof MultiTerminus) {
        targetTerm = (MultiTerminus) target;
        terminus = true;
    } else {
        logger.warning("Improper method call.");
    }
    // Record the pre-change electrostatic energy.
    ffe.energy(false, false);
    final double previousElectrostaticEnergy = ffe.getElectrostaticEnergy();
    // Write the pre-titration change snapshot.
    writeSnapshot(true, StepType.TITRATE, config.snapshots);
    String startString = target.toString();
    String startName = target.getName();
    double pKaref = 0;
    double dG_ref = 0;
    Titration titration = null;
    final TitrationType titrationType;
    if (terminus) {
        if (targetTerm.end == MultiTerminus.END.NTERM) {
            pKaref = 10.0;
            dG_ref = 0.0;
        } else {
            pKaref = 3.0;
            dG_ref = 0.0;
        }
        titrationType = targetTerm.titrateTerminus_v1(thermostat.getCurrentTemperature());
    } else {
        logger.info(format("targetMulti:  %s", targetMulti.toString()));
        logger.info(format("getActive:    %s", targetMulti.getActive().toString()));
        logger.info(format("titrationMap: %s", Arrays.toString(titrationMap.get(targetMulti.getActive()).toArray())));
        // Choose from the list of available titrations for the active residue.
        List<Titration> avail = titrationMap.get(targetMulti.getActive());
        titration = avail.get(rng.nextInt(avail.size()));
        // Perform the chosen titration.
        titrationType = performTitration(targetMulti, titration, config.inactivateBackground);
        reInitialize(true, true);
        // Test the MC criterion for a titration step.
        pKaref = titration.pKa;
        dG_ref = titration.refEnergy;
    }
    // Write the post-titration change snapshot.
    writeSnapshot(true, StepType.TITRATE, config.snapshots);
    double temperature = thermostat.getCurrentTemperature();
    double kT = BOLTZMANN * temperature;
    ffe.energy(false, false);
    final double currentElectrostaticEnergy = ffe.getElectrostaticEnergy();
    final double dG_elec = currentElectrostaticEnergy - previousElectrostaticEnergy;
    if (config.zeroReferenceEnergies) {
        dG_ref = 0.0;
    }
    if (config.refOverride.isPresent()) {
        dG_ref = config.refOverride.getAsDouble();
    }
    /**
     * dG_elec = electrostatic energy component of the titratable residue
     * dG_ref = electrostatic component of the transition energy for the
     * reference compound
     */
    double prefix = Math.log(10) * kT * (pH - pKaref);
    if (titrationType == TitrationType.DEPROT) {
        prefix = -prefix;
    }
    // Either positive ref == deprotonation or == standard -> nonstandard transition.
    if (titrationType == TitrationType.PROT) {
        dG_ref = -dG_ref;
    }
    double postfix = dG_elec - dG_ref;
    double dG_MC = prefix + postfix;
    StringBuilder sb = new StringBuilder();
    sb.append(String.format(" Assessing possible MC protonation step:\n"));
    sb.append(String.format("     %s --> %s\n", startString, target.toString()));
    sb.append(String.format("     dG_ref:  %7.2f                pKaref:  %7.2f\n", dG_ref, pKaref));
    sb.append(String.format("     pH_term: %9.4f              elec_term: %10.4f\n", prefix, postfix));
    sb.append(String.format("     dG_elec: %9.4f              dG_MC:     %10.4f\n", dG_elec, dG_MC));
    sb.append(String.format("     -----\n"));
    // Test Monte-Carlo criterion.
    if (dG_MC < 0 && config.mcOverride != MCOverride.REJECT) {
        sb.append(String.format("     Accepted!"));
        logger.info(sb.toString());
        numMovesAccepted++;
        return true;
    }
    double criterion = exp(-dG_MC / kT);
    double metropolis = random();
    sb.append(String.format("     crit:    %9.4f              rng:       %10.4f\n", criterion, metropolis));
    if ((metropolis < criterion && config.mcOverride != MCOverride.REJECT) || config.mcOverride == MCOverride.ACCEPT) {
        numMovesAccepted++;
        molDyn.reInit();
        long took = System.nanoTime() - startTime;
        sb.append(String.format("     Accepted!                                                %1.3f", took * NS_TO_SEC));
        logger.info(sb.toString());
        return true;
    }
    // Move was rejected, undo the titration state change.
    performTitration(targetMulti, titration, config.inactivateBackground);
    reInitialize(true, true);
    long took = System.nanoTime() - startTime;
    sb.append(String.format("     Denied.                                                  %1.3f", took * NS_TO_SEC));
    logger.info(sb.toString());
    return false;
}

Example 19 with MultiResidue

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

the class PhMD method recursiveMap.

/**
 * Finds Titration definitions for the given Residue and adds them to the given MultiResidue.
 * For three-state transitions, simply populate the enumeration with multiple titrations
 * from a shared state and this will include them in MultiResidue construction.
 */
private void recursiveMap(Residue member, MultiResidue multiRes) {
    // Map titrations for this member.
    Titration[] titrations = Titration.multiLookup(member);
    titrationMap.put(member, Arrays.asList(titrations));
    // For each titration, check whether it needs added as a MultiResidue option.
    for (Titration titration : titrations) {
        // Allow manual override of Histidine treatment.
        if ((titration.deprotForm == AminoAcid3.HID && config.histidineMode == HistidineMode.HIE_ONLY) || (titration.deprotForm == AminoAcid3.HIE && config.histidineMode == HistidineMode.HID_ONLY)) {
            continue;
        }
        // Find all the choices currently available to this MultiResidue.
        List<AminoAcid3> choices = new ArrayList<>();
        for (Residue choice : multiRes.getConsideredResidues()) {
            choices.add(choice.getAminoAcid3());
        }
        // If this Titration target is not a choice for the MultiResidue, then add it.
        if (!choices.contains(titration.protForm) || !(choices.contains(titration.deprotForm))) {
            String targetName = (member.getAminoAcid3() == titration.protForm) ? titration.deprotForm.toString() : titration.protForm.toString();
            int resNumber = member.getResidueNumber();
            ResidueType resType = member.getResidueType();
            Residue newChoice = new Residue(targetName, resNumber, resType);
            multiRes.addResidue(newChoice);
            titrationMap.put(newChoice, Arrays.asList(Titration.multiLookup(newChoice)));
        }
    }
}