Examples with TitrationType ffx.potential.extended.TitrationUtils.TitrationType used on opensource projects

Example 1 with TitrationType

use of ffx.potential.extended.TitrationUtils.TitrationType in project ffx by mjschnie.

the class PhMD method tryComboStep.

/**
 * Attempt a combination titration/rotamer MC move.
 *
 * @param targetMulti
 * @return accept/reject
 */
private boolean tryComboStep(MultiResidue targetMulti) {
    if (CAUTIOUS) {
        throw new UnsupportedOperationException();
    }
    // Record the pre-change total energy.
    double previousTotalEnergy = currentTotalEnergy();
    double previousElectrostaticEnergy = currentElectrostaticEnergy();
    // Write the pre-combo snapshot.
    writeSnapshot(true, StepType.COMBO, config.snapshots);
    String startString = targetMulti.toString();
    String startName = targetMulti.getActive().getName();
    // Choose from the list of available titrations for the active residue.
    List<Titration> avail = titrationMap.get(targetMulti.getActive());
    Titration titration = avail.get(rng.nextInt(avail.size()));
    // Perform the chosen titration.
    TitrationType titrationType = performTitration(targetMulti, titration, config.inactivateBackground);
    // Change rotamer state, but first save coordinates so we can return to them if 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);
    // Swap to the new rotamer.
    Rotamer[] rotamers = residue.getRotamers(library);
    int rotaRand = rng.nextInt(rotamers.length);
    RotamerLibrary.applyRotamer(residue, rotamers[rotaRand]);
    // Write the post-combo snapshot.
    writeSnapshot(false, StepType.COMBO, config.snapshots);
    // Evaluate both MC criteria.
    String endName = targetMulti.getActive().getName();
    // Evaluate the titration probability of the step.
    double pKaref = titration.pKa;
    double dG_ref = titration.refEnergy;
    double temperature = thermostat.getCurrentTemperature();
    double kT = BOLTZMANN * temperature;
    double dG_elec = currentElectrostaticEnergy() - previousElectrostaticEnergy;
    if (config.zeroReferenceEnergies) {
        dG_ref = 0.0;
    }
    double prefix = Math.log(10) * kT * (pH - pKaref);
    if (titrationType == TitrationType.DEPROT) {
        prefix = -prefix;
    }
    double postfix = dG_elec - dG_ref;
    double dG_titr = prefix + postfix;
    double titrCriterion = exp(-dG_titr / kT);
    // Evaluate the rotamer probability of the step.
    double dG_rota = currentTotalEnergy() - previousTotalEnergy;
    double rotaCriterion = exp(-dG_rota / kT);
    StringBuilder sb = new StringBuilder();
    sb.append(String.format(" Assessing possible MC combo step:\n"));
    sb.append(String.format("     dG_elec: %16.8f\n", dG_elec));
    sb.append(String.format("     dG_titr: %16.8f\n", dG_titr));
    sb.append(String.format("     dG_rota: %16.8f\n", dG_rota));
    sb.append(String.format("     -----\n"));
    // Automatic acceptance if both energy changes are favorable.
    if (dG_titr < 0 && dG_rota < 0 && config.mcOverride != MCOverride.REJECT) {
        sb.append(String.format("     Accepted!"));
        logger.info(sb.toString());
        numMovesAccepted++;
        propagateInactiveResidues(titratingMultis, false);
        return true;
    } else {
        // Conditionally accept based on combined probabilities.
        if (dG_titr < 0 || config.mcOverride == MCOverride.ACCEPT) {
            titrCriterion = 1.0;
        }
        if (dG_rota < 0) {
            rotaCriterion = 1.0;
        }
        if (config.mcOverride == MCOverride.REJECT) {
            titrCriterion = 0.0;
        }
        double metropolis = random();
        double comboCriterion = titrCriterion * rotaCriterion;
        sb.append(String.format("     titrCrit:   %9.4f\n", titrCriterion));
        sb.append(String.format("     rotaCrit:   %9.4f\n", rotaCriterion));
        sb.append(String.format("     criterion:  %9.4f\n", comboCriterion));
        sb.append(String.format("     rng:        %9.4f\n", metropolis));
        if (metropolis < comboCriterion) {
            sb.append(String.format("     Accepted!"));
            logger.info(sb.toString());
            numMovesAccepted++;
            propagateInactiveResidues(titratingMultis, false);
            return true;
        } else {
            // Move was denied.
            sb.append(String.format("     Denied."));
            logger.info(sb.toString());
            // Undo both pieces of the rejected move IN THE RIGHT ORDER.
            RotamerLibrary.applyRotamer(residue, origCoordsRotamer);
            performTitration(targetMulti, titration, config.inactivateBackground);
            ffe.reInit();
            molDyn.reInit();
            return false;
        }
    }
}

Example 2 with TitrationType

use of ffx.potential.extended.TitrationUtils.TitrationType in project ffx by mjschnie.

the class MultiTerminus method titrateTerminus_v1.

/**
 * Changes the charge state of this MultiTerminus.
 * Keep existing Atom objects but updates types, bonded terms, and builds new proton if necessary.
 */
public TitrationType titrateTerminus_v1(double temperature) {
    logger.info(String.format(" Titrating residue %s (currently %d).", this.toString(), (isCharged ? 1 : 0)));
    // StringBuilder sb = new StringBuilder();
    // sb.append(" Contents of children: ");
    // for (Atom atom : this.getAtomList()) {
    // sb.append(String.format("%s, ", atom.getName()));
    // }
    // logger.info(sb.toString());
    /**
     * Get references to the backbone atoms.
     */
    TitrationType titrationType = (isCharged) ? TitrationType.DEPROT : TitrationType.PROT;
    Atom N = getBBAtom("N");
    Atom CA = getBBAtom("CA");
    Atom C = getBBAtom("C");
    Atom O = getBBAtom("O");
    Atom OXT = getBBAtom("OXT");
    Atom H1 = getBBAtom("H1");
    Atom H2 = getBBAtom("H2");
    Atom H3 = getBBAtom("H3");
    Atom HA = getBBAtom("HA");
    Atom OH = getBBAtom("OH");
    Atom HO = getBBAtom("HO");
    String resName = C.getResidueName();
    int resSeq = C.getResidueNumber();
    Character chainID = C.getChainID();
    List<Atom> typeChanged = new ArrayList<>();
    if (DEBUG) {
        printBonds();
    }
    if (end == END.NTERM) {
        if (isCharged) {
            if (rolsH3.isEmpty()) {
                for (Angle a : H3.getAngles()) {
                    rolsH3.add(a);
                }
                for (Torsion t : H3.getTorsions()) {
                    rolsH3.add(t);
                }
            }
            N.setAtomType(forceField.getAtomType(Integer.toString(BB_TYPE.N.neutType)));
            H1.setAtomType(forceField.getAtomType(Integer.toString(BB_TYPE.H1.neutType)));
            H2.setAtomType(forceField.getAtomType(Integer.toString(BB_TYPE.H2.neutType)));
            // intxyz(H1, N, 1.02, CA, 109.5, C, 180.0, 0);
            // intxyz(H2, N, 1.02, CA, 109.5, C, 0.0, 0);
            bondH3 = N.getBond(H3);
            bondH3.removeFromParent();
            for (Angle a : H3.getAngles()) {
                a.removeFromParent();
            }
            for (Torsion t : H3.getTorsions()) {
                t.removeFromParent();
            }
            H3.removeFromParent();
            H3.setParent(null);
            H3.setUse(false);
            uberH3 = H3;
            typeChanged.add(N);
            typeChanged.add(H1);
            typeChanged.add(H2);
        // logger.info(String.format(" Finished titration. H3 status: %b %b %b",
        // N.getBond(H3) == null, this.getAtomNode().contains(H3) == null, H3.getParent() == null));
        } else {
            if (H3 != null) {
                logger.severe("N-terminal found in incorrect charge state.");
            }
            if (uberH3 == null || bondH3 == null) {
                logger.severe("Please start with (N-)termini in the charged state.");
            }
            H3 = uberH3;
            N.setAtomType(forceField.getAtomType(Integer.toString(BB_TYPE.N.chrgType)));
            H1.setAtomType(forceField.getAtomType(Integer.toString(BB_TYPE.H1.chrgType)));
            H2.setAtomType(forceField.getAtomType(Integer.toString(BB_TYPE.H2.chrgType)));
            intxyz(H1, N, 1.02, CA, 109.5, C, 180.0, 0);
            intxyz(H2, N, 1.02, CA, 109.5, C, 60.0, 0);
            intxyz(H3, N, 1.02, CA, 109.5, C, -60.0, 0);
            maxwellMe(H3, temperature);
            this.getAtomNode().add(H3);
            this.add(bondH3);
            // logger.info("Number of rolsH3 terms: " + rolsH3.size());
            for (ROLS rols : rolsH3) {
                if (rols instanceof Angle) {
                    this.add((Angle) rols);
                } else if (rols instanceof Torsion) {
                    this.add((Torsion) rols);
                }
            }
            H3.setParent(this.getAtomNode());
            H3.setUse(true);
            typeChanged.add(N);
            typeChanged.add(H1);
            typeChanged.add(H2);
            typeChanged.add(H3);
        // logger.info(String.format(" Finished titration. H3 statuses: "
        // + "(They have each other: %b %b) (I have them: %b %b) (I have bond: %b)",
        // N.getBond(H3) != null, H3.getBond(N) != null,
        // this.getAtomNode().contains(H3) != null, H3.getParent() == this.getAtomNode(),
        // this.getBondList().contains(bondH3)));
        // logger.info(String.format(" Bonds from H3: %s %s",
        // H3.getBonds().get(0).get1_2(H3).describe(Atom.Descriptions.INDEX_NAME),
        // H3.getBonds().get(0).get1_2(H3).getBonds().get(0).get1_2(H3.getBonds().get(0).get1_2(H3)).describe(Atom.Descriptions.INDEX_NAME)));
        }
    } else if (end == END.CTERM) {
        if (isCharged) {
            OXT.setName("OH");
            OH = OXT;
            if (HO != null) {
                logger.warning("C-terminal in unusual charge state.");
            }
            C.setAtomType(forceField.getAtomType(Integer.toString(BB_TYPE.C.neutType)));
            O.setAtomType(forceField.getAtomType(Integer.toString(BB_TYPE.O.neutType)));
            OH.setAtomType(forceField.getAtomType(Integer.toString(BB_TYPE.OH.neutType)));
            if (uberHO == null) {
                // Gotta build the HO and all its bonded terms.
                uberHO = new Atom(mola.getAtomArray().length, "HO", OH.getAltLoc(), new double[3], resName, resSeq, chainID, OH.getOccupancy(), OH.getTempFactor(), OH.getSegID(), true);
                uberHO.setAtomType(forceField.getAtomType(Integer.toString(BB_TYPE.HO.neutType)));
                bondHO = new Bond(OH, uberHO);
                bondHO.setBondType(forceField.getBondType(String.format("%d %d", OH.getAtomType().atomClass, uberHO.getAtomType().atomClass)));
                Angle a1 = Angle.angleFactory(bondHO, OH.getBond(C), forceField);
                Torsion t1 = Torsion.torsionFactory(O.getBond(C), C.getBond(OH), bondHO, forceField);
                Torsion t2 = Torsion.torsionFactory(CA.getBond(C), C.getBond(OH), bondHO, forceField);
                this.add(a1);
                this.add(t1);
                this.add(t2);
            }
            HO = uberHO;
            intxyz(HO, OXT, 1.02, C, 109.5, CA, -1.7, 0);
            maxwellMe(HO, temperature);
            this.getAtomNode().add(HO);
            this.add(bondHO);
            HO.setParent(this.getAtomNode());
            HO.setUse(true);
            // logger.info("Number of rolsHO terms: " + rolsHO.size());
            for (ROLS rols : rolsHO) {
                if (rols instanceof Angle) {
                    this.add((Angle) rols);
                } else if (rols instanceof Torsion) {
                    this.add((Torsion) rols);
                }
            }
            typeChanged.add(C);
            typeChanged.add(O);
            typeChanged.add(OH);
            typeChanged.add(HO);
        } else {
            if (rolsHO.isEmpty()) {
                rolsHO = new ArrayList<>();
                for (Angle a : HO.getAngles()) {
                    rolsHO.add(a);
                }
                for (Torsion t : HO.getTorsions()) {
                    rolsHO.add(t);
                }
            }
            C.setAtomType(forceField.getAtomType(Integer.toString(BB_TYPE.C.chrgType)));
            O.setAtomType(forceField.getAtomType(Integer.toString(BB_TYPE.O.chrgType)));
            OH.setAtomType(forceField.getAtomType(Integer.toString(BB_TYPE.OH.chrgType)));
            OH.setName("OXT");
            OXT = OH;
            bondHO = OH.getBond(HO);
            bondHO.removeFromParent();
            for (Angle a : HO.getAngles()) {
                a.removeFromParent();
            }
            for (Torsion t : HO.getTorsions()) {
                t.removeFromParent();
            }
            HO.removeFromParent();
            HO.setParent(null);
            HO.setUse(false);
            uberHO = HO;
            typeChanged.add(C);
            typeChanged.add(O);
            typeChanged.add(OH);
        }
    }
    updateGeometry();
    updateBondedTerms();
    isCharged = !isCharged;
    forceFieldEnergy.reInit();
    if (DEBUG) {
        printBonds();
    }
    return titrationType;
}

Example 3 with TitrationType

use of ffx.potential.extended.TitrationUtils.TitrationType 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 4 with TitrationType

use of ffx.potential.extended.TitrationUtils.TitrationType in project ffx by mjschnie.

the class PhMD method tryTerminusTitration.

/**
 * Perform a titration MC move.
 *
 * @param targetMulti
 * @return accept/reject
 */
private boolean tryTerminusTitration(MultiTerminus target) {
    if (CAUTIOUS) {
        throw new UnsupportedOperationException();
    }
    // Record the pre-change electrostatic energy.
    double previousElectrostaticEnergy = currentElectrostaticEnergy();
    // 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;
    TitrationType type = null;
    if (target.end == MultiTerminus.END.NTERM) {
        pKaref = 8.23;
        dG_ref = 85.4929;
        type = target.isCharged ? TitrationType.DEPROT : TitrationType.PROT;
    } else if (target.end == MultiTerminus.END.CTERM) {
        pKaref = 3.55;
        dG_ref = -61.3825;
        type = target.isCharged ? TitrationType.PROT : TitrationType.DEPROT;
    }
    boolean beganCharged = target.isCharged;
    target.titrateTerminus_v1(thermostat.getCurrentTemperature());
    // Write the post-titration change snapshot.
    writeSnapshot(true, StepType.TITRATE, config.snapshots);
    double temperature = thermostat.getCurrentTemperature();
    double kT = BOLTZMANN * temperature;
    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 pHterm = Math.log(10) * kT * (pH - pKaref);
    if (type == TitrationType.DEPROT) {
        pHterm = -pHterm;
    }
    // Either positive ref == deprotonation or == standard -> nonstandard transition.
    if (type == TitrationType.PROT) {
        dG_ref = -dG_ref;
    }
    double ddGterm = dG_elec - dG_ref;
    double dG_MC = pHterm + ddGterm;
    // StringBuilder sb = new StringBuilder();
    // sb.append(String.format(" Assessing possible MC protonation step:\n"));
    // sb.append(String.format("     %s --> %s\n", startString, targetMulti.toString()));
    // sb.append(String.format("     pKaref:  %7.2f\n", pKaref));
    // sb.append(String.format("     dG_ref:  %7.2f\n", dG_ref));
    // sb.append(String.format("     dG_elec: %16.8f\n", dG_elec));
    // sb.append(String.format("     dG_MC:   %16.8f\n", dG_MC));
    // sb.append(String.format("     -----\n"));
    StringBuilder sb = new StringBuilder();
    sb.append(String.format(" Assessing possible MC protonation step:\n"));
    if (beganCharged) {
        sb.append(String.format("     %sc --> %sn\n", startString, target.toString()));
    } else {
        sb.append(String.format("     %sn --> %sc\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", pHterm, ddGterm));
    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 || config.mcOverride == MCOverride.ONCE) {
        numMovesAccepted++;
        long took = System.nanoTime() - startTime;
        sb.append(String.format("     Accepted!                                                %1.3f", took * NS_TO_SEC));
        logger.info(sb.toString());
        if (config.mcOverride == MCOverride.ONCE) {
            config.mcOverride = MCOverride.REJECT;
        }
        return true;
    }
    // Move was rejected, undo the titration state change.
    target.titrateTerminus_v1(thermostat.getCurrentTemperature());
    long took = System.nanoTime() - startTime;
    sb.append(String.format("     Denied.                                                  %1.3f", took * NS_TO_SEC));
    logger.info(sb.toString());
    return false;
}