use of ffx.potential.utils.EnergyException in project ffx by mjschnie.
the class MolecularDynamics method run.
/**
* {@inheritDoc}
*/
@Override
public void run() {
done = false;
terminate = false;
/**
* Set the target temperature.
*/
thermostat.setTargetTemperature(targetTemperature);
thermostat.setQuiet(quiet);
if (integrator instanceof Stochastic) {
Stochastic stochastic = (Stochastic) integrator;
stochastic.setTemperature(targetTemperature);
}
/**
* Set the step size.
*/
integrator.setTimeStep(dt);
if (!initialized) {
/**
* Initialize from a restart file.
*/
if (loadRestart) {
Crystal crystal = molecularAssembly.getCrystal();
if (!dynFilter.readDYN(restartFile, crystal, x, v, a, aPrevious)) {
String message = " Could not load the restart file - dynamics terminated.";
logger.log(Level.WARNING, message);
done = true;
return;
} else {
molecularAssembly.getPotentialEnergy().setCrystal(crystal);
}
} else {
/**
* Initialize from using current atomic coordinates.
*/
potential.getCoordinates(x);
/**
* Initialize atomic velocities from a Maxwell-Boltzmann
* distribution or set to 0.
*/
if (initVelocities) {
thermostat.maxwell(targetTemperature);
} else {
fill(v, 0.0);
}
}
} else {
/**
* If MD has already been run (ie. Annealing or RepEx), then
* initialize velocities if requested.
*/
if (initVelocities) {
thermostat.maxwell(targetTemperature);
}
}
/**
* Compute the current potential energy.
*/
potential.setScaling(null);
try {
currentPotentialEnergy = potential.energyAndGradient(x, grad);
} catch (EnergyException ex) {
writeStoredSnapshots();
throw ex;
}
/**
* Initialize current and previous accelerations, unless they were just
* loaded from a restart file.
*/
if (!loadRestart || initialized) {
for (int i = 0; i < numberOfVariables; i++) {
a[i] = -Thermostat.convert * grad[i] / mass[i];
}
if (aPrevious != null) {
arraycopy(a, 0, aPrevious, 0, numberOfVariables);
}
}
/**
* Compute the current kinetic energy.
*/
thermostat.kineticEnergy();
currentKineticEnergy = thermostat.getKineticEnergy();
currentTemperature = thermostat.getCurrentTemperature();
currentTotalEnergy = currentKineticEnergy + currentPotentialEnergy;
initialized = true;
logger.info(format("\n %8s %12s %12s %12s %8s %8s", "Time", "Kinetic", "Potential", "Total", "Temp", "CPU"));
logger.info(format(" %8s %12s %12s %12s %8s %8s", "psec", "kcal/mol", "kcal/mol", "kcal/mol", "K", "sec"));
logger.info(format(" %8s %12.4f %12.4f %12.4f %8.2f", "", currentKineticEnergy, currentPotentialEnergy, currentTotalEnergy, currentTemperature));
/**
* Store the initialized state.
*/
storeState();
/**
* Integrate Newton's equations of motion for the requested number of
* steps, unless early termination is requested.
*/
time = System.nanoTime();
for (int step = 1; step <= nSteps; step++) {
/* Notify MonteCarlo handlers such as PhMD or rotamer drivers. */
if (monteCarloListener != null && mcNotification == MonteCarloNotification.EACH_STEP) {
monteCarloListener.mcUpdate(thermostat.getCurrentTemperature());
}
/**
* Do the half-step thermostat operation.
*/
thermostat.halfStep(dt);
/**
* Do the half-step integration operation.
*/
integrator.preForce(potential);
double priorPE = currentPotentialEnergy;
/**
* Compute the potential energy and gradients.
*/
try {
currentPotentialEnergy = potential.energyAndGradient(x, grad);
} catch (EnergyException ex) {
writeStoredSnapshots();
throw ex;
}
detectAtypicalEnergy(priorPE, defaultDeltaPEThresh);
/**
* Add the potential energy of the slow degrees of freedom.
*/
if (integrator instanceof Respa) {
Respa r = (Respa) integrator;
currentPotentialEnergy += r.getHalfStepEnergy();
}
/**
* Do the full-step integration operation.
*/
integrator.postForce(grad);
/**
* Compute the full-step kinetic energy.
*/
thermostat.kineticEnergy();
/**
* Do the full-step thermostat operation.
*/
thermostat.fullStep(dt);
/**
* Recompute the kinetic energy after the full-step thermostat
* operation.
*/
thermostat.kineticEnergy();
/**
* Remove center of mass motion ever ~100 steps.
*/
if (thermostat.getRemoveCenterOfMassMotion() && step % removeCOMMotionFrequency == 0) {
thermostat.centerOfMassMotion(true, false);
}
/**
* Collect current kinetic energy, temperature, and total energy.
*/
currentKineticEnergy = thermostat.getKineticEnergy();
currentTemperature = thermostat.getCurrentTemperature();
currentTotalEnergy = currentKineticEnergy + currentPotentialEnergy;
/**
* Update atomic velocity, acceleration and previous acceleration.
*/
potential.setVelocity(v);
potential.setAcceleration(a);
potential.setPreviousAcceleration(aPrevious);
/**
* Update extended system variables if present.
*/
if (esvSystem != null) {
esvSystem.propagateESVs(currentTemperature, dt, step * dt);
}
/**
* Log the current state every printFrequency steps.
*/
totalSimTime += dt;
if (step % printFrequency == 0) {
/**
* Original print statement
*/
time = System.nanoTime() - time;
mdTime = time;
logger.info(format(" %7.3e %12.4f %12.4f %12.4f %8.2f %8.3f", totalSimTime, currentKineticEnergy, currentPotentialEnergy, currentTotalEnergy, currentTemperature, time * NS2SEC));
time = System.nanoTime();
// Shirts et al. logging info
Crystal crystal = molecularAssembly.getCrystal();
double volume = crystal.getUnitCell().volume;
// Shirts et al. print statement
// time = System.nanoTime() - time;
// logger.info(format("Shirts %7.3e %12.4f %12.4f %12.4f %12.4f %8.2f %8.3f",
// totalSimTime, currentKineticEnergy, currentPotentialEnergy,
// currentTotalEnergy, volume, currentTemperature, time * NS2SEC));
// time = System.nanoTime();
}
if (step % printEsvFrequency == 0 && esvSystem != null) {
logger.info(format(" %7.3e %s", totalSimTime, esvSystem.getLambdaList()));
}
/**
* Write out snapshots in selected format every
* saveSnapshotFrequency steps.
*/
if (saveSnapshotFrequency > 0 && step % saveSnapshotFrequency == 0) {
for (AssemblyInfo ai : assemblies) {
if (ai.archiveFile != null && !saveSnapshotAsPDB) {
if (ai.xyzFilter.writeFile(ai.archiveFile, true)) {
logger.info(String.format(" Appended snap shot to %s", ai.archiveFile.getName()));
} else {
logger.warning(String.format(" Appending snap shot to %s failed", ai.archiveFile.getName()));
}
} else if (saveSnapshotAsPDB) {
if (ai.pdbFilter.writeFile(ai.pdbFile, false)) {
logger.info(String.format(" Wrote PDB file to %s", ai.pdbFile.getName()));
} else {
logger.warning(String.format(" Writing PDB file to %s failed.", ai.pdbFile.getName()));
}
}
}
}
/**
* Write out restart files every saveRestartFileFrequency steps.
*/
if (saveRestartFileFrequency > 0 && step % saveRestartFileFrequency == 0) {
if (dynFilter.writeDYN(restartFile, molecularAssembly.getCrystal(), x, v, a, aPrevious)) {
logger.info(String.format(" Wrote dynamics restart file to " + restartFile.getName()));
} else {
logger.info(String.format(" Writing dynamics restart file to " + restartFile.getName() + " failed"));
}
}
/**
* Notify the algorithmListener.
*/
if (algorithmListener != null && step % printFrequency == 0) {
// algorithmListener.algorithmUpdate(molecularAssembly);
for (AssemblyInfo assembly : assemblies) {
// Probably unwise to parallelize this, so that it doesn't
// hit the GUI with parallel updates.
algorithmListener.algorithmUpdate(assembly.getAssembly());
}
}
/**
* Check for a termination request.
*/
if (terminate) {
logger.info(String.format("\n Terminating after %8d time steps\n", step));
break;
}
}
/**
* Log normal completion.
*/
if (!terminate) {
logger.info(String.format(" Completed %8d time steps\n", nSteps));
}
/**
* Reset the done and terminate flags.
*/
done = true;
terminate = false;
if (monteCarloListener != null && mcNotification == MonteCarloNotification.AFTER_DYNAMICS) {
monteCarloListener.mcUpdate(thermostat.getCurrentTemperature());
}
}
use of ffx.potential.utils.EnergyException in project ffx by mjschnie.
the class ForceFieldEnergy method fillGradients.
/**
* Private method for internal use, so we don't have subclasses calling super.energy, and this class delegating to
* the subclass's getGradients method.
* @param g Gradient array to fill.
* @return Gradient array.
*/
private double[] fillGradients(double[] g) {
assert (g != null);
double[] grad = new double[3];
int n = getNumberOfVariables();
if (g.length < n) {
g = new double[n];
}
int index = 0;
for (int i = 0; i < nAtoms; i++) {
Atom a = atoms[i];
if (a.isActive()) {
a.getXYZGradient(grad);
double gx = grad[0];
double gy = grad[1];
double gz = grad[2];
if (Double.isNaN(gx) || Double.isInfinite(gx) || Double.isNaN(gy) || Double.isInfinite(gy) || Double.isNaN(gz) || Double.isInfinite(gz)) {
/*String message = format("The gradient of atom %s is (%8.3f,%8.3f,%8.3f).",
a.toString(), gx, gy, gz);*/
StringBuilder sb = new StringBuilder(format("The gradient of atom %s is (%8.3f,%8.3f,%8.3f).", a.toString(), gx, gy, gz));
double[] vals = new double[3];
a.getVelocity(vals);
sb.append(format("\n Velocities: %8.3g %8.3g %8.3g", vals[0], vals[1], vals[2]));
a.getAcceleration(vals);
sb.append(format("\n Accelerations: %8.3g %8.3g %8.3g", vals[0], vals[1], vals[2]));
a.getPreviousAcceleration(vals);
sb.append(format("\n Previous accelerations: %8.3g %8.3g %8.3g", vals[0], vals[1], vals[2]));
// logger.severe(message);
throw new EnergyException(sb.toString());
}
g[index++] = gx;
g[index++] = gy;
g[index++] = gz;
}
}
return g;
}
use of ffx.potential.utils.EnergyException in project ffx by mjschnie.
the class ForceFieldEnergyOpenMM method fillGradients.
/**
* Private method for internal use, so we don't have subclasses calling
* super.energy, and this class delegating to the subclass's getGradients
* method.
*
* @param g Gradient array to fill.
* @return Gradient array.
*/
public double[] fillGradients(double[] g) {
assert (g != null);
int n = getNumberOfVariables();
if (g.length < n) {
g = new double[n];
}
int index = 0;
Atom[] atoms = molecularAssembly.getAtomArray();
int nAtoms = atoms.length;
for (int i = 0; i < nAtoms; i++) {
Atom a = atoms[i];
if (a.isActive()) {
OpenMM_Vec3 posInNm = OpenMM_Vec3Array_get(forces, i);
/**
* Convert OpenMM Forces in KJ/Nm into an FFX gradient in
* Kcal/A.
*/
double gx = -posInNm.x * OpenMM_NmPerAngstrom * OpenMM_KcalPerKJ;
double gy = -posInNm.y * OpenMM_NmPerAngstrom * OpenMM_KcalPerKJ;
double gz = -posInNm.z * OpenMM_NmPerAngstrom * OpenMM_KcalPerKJ;
if (Double.isNaN(gx) || Double.isInfinite(gx) || Double.isNaN(gy) || Double.isInfinite(gy) || Double.isNaN(gz) || Double.isInfinite(gz)) {
/*String message = format("The gradient of atom %s is (%8.3f,%8.3f,%8.3f).",
a.toString(), gx, gy, gz);*/
StringBuilder sb = new StringBuilder(format("The gradient of atom %s is (%8.3f,%8.3f,%8.3f).", a.toString(), gx, gy, gz));
double[] vals = new double[3];
a.getVelocity(vals);
sb.append(format("\n Velocities: %8.3g %8.3g %8.3g", vals[0], vals[1], vals[2]));
a.getAcceleration(vals);
sb.append(format("\n Accelerations: %8.3g %8.3g %8.3g", vals[0], vals[1], vals[2]));
a.getPreviousAcceleration(vals);
sb.append(format("\n Previous accelerations: %8.3g %8.3g %8.3g", vals[0], vals[1], vals[2]));
// logger.severe(message);
throw new EnergyException(sb.toString());
}
a.setXYZGradient(gx, gy, gz);
g[index++] = gx;
g[index++] = gy;
g[index++] = gz;
}
}
return g;
}
use of ffx.potential.utils.EnergyException in project ffx by mjschnie.
the class QuadTopologyEnergy method energyAndGradient.
@Override
public double energyAndGradient(double[] x, double[] g, boolean verbose) {
region.setX(x);
region.setG(g);
region.setVerbose(verbose);
try {
team.execute(region);
} catch (Exception ex) {
throw new EnergyException(String.format(" Exception in calculating quad-topology energy: %s", ex.toString()), false);
}
if (verbose) {
logger.info(String.format(" Total quad-topology energy: %12.4f", totalEnergy));
}
return totalEnergy;
}
use of ffx.potential.utils.EnergyException in project ffx by mjschnie.
the class DualTopologyEnergy method energy.
@Override
public double energy(double[] x, boolean verbose) {
try {
region.setX(x);
region.setVerbose(verbose);
team.execute(region);
} catch (Exception ex) {
throw new EnergyException(String.format(" Exception in calculating dual-topology energy: %s", ex.toString()), false);
}
return totalEnergy;
}
Aggregations