use of ffx.potential.extended.TitrationUtils.Titration 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;
}
}
}
use of ffx.potential.extended.TitrationUtils.Titration 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;
}
use of ffx.potential.extended.TitrationUtils.Titration 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;
}
use of ffx.potential.extended.TitrationUtils.Titration 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)));
}
}
}
Aggregations