use of ffx.potential.utils.SystemTemperatureException 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;
}
}
Aggregations