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

Example 1 with MultiTerminus

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

the class PhMD method readyup.

private void readyup() {
    // Create MultiTerminus objects to wrap termini.
    if (config.titrateTermini) {
        for (Residue res : mola.getResidueList()) {
            if (res.getPreviousResidue() == null || res.getNextResidue() == null) {
                MultiTerminus multiTerminus = new MultiTerminus(res, ff, ffe, mola);
                Polymer polymer = findResiduePolymer(res, mola);
                polymer.addMultiTerminus(res, multiTerminus);
                reInitialize(true, false);
                titratingTermini.add(multiTerminus);
                logger.info(String.format(" Titrating: %s", multiTerminus));
            }
        }
    }
    /* Create containers for titratables: MultiResidues for discrete, ExtendedVariables for continuous. */
    if (distribution == Distribution.CONTINUOUS) {
        esvSystem = new ExtendedSystem(mola);
        esvSystem.setConstantPh(pH);
        for (Residue res : chosenResidues) {
            MultiResidue multi = TitrationUtils.titratingMultiresidueFactory(mola, res);
            TitrationESV esv = new TitrationESV(esvSystem, multi);
            titratingESVs.add(esv);
            for (Residue background : multi.getInactive()) {
                inactivateResidue(background);
            }
            esvSystem.addVariable(esv);
        }
        ffe.attachExtendedSystem(esvSystem);
        logger.info(format(" Continuous pHMD readied with %d residues.", titratingESVs.size()));
    } else {
        for (Residue res : chosenResidues) {
            // Create MultiResidue objects to wrap titratables.
            MultiResidue multiRes = new MultiResidue(res, ff, ffe);
            Polymer polymer = findResiduePolymer(res, mola);
            polymer.addMultiResidue(multiRes);
            recursiveMap(res, multiRes);
            // Switch back to the original form and ready the ForceFieldEnergy.
            multiRes.setActiveResidue(res);
            reInitialize(true, false);
            titratingMultis.add(multiRes);
            logger.info(String.format(" Titrating: %s", multiRes));
        }
        logger.info(format(" Discrete MCMD readied with %d residues.", titratingMultis.size()));
    }
    switch(distribution) {
        default:
        case DISCRETE:
            molDyn.setMonteCarloListener(this, MonteCarloNotification.EACH_STEP);
            break;
        case CONTINUOUS:
            ffe.attachExtendedSystem(esvSystem);
            molDyn.attachExtendedSystem(esvSystem, 100);
            break;
    }
}

Example 2 with MultiTerminus

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