Search in sources :

Example 1 with EnergyException

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());
    }
}
Also used : Respa(ffx.algorithms.integrators.Respa) Stochastic(ffx.algorithms.integrators.Stochastic) EnergyException(ffx.potential.utils.EnergyException) Crystal(ffx.crystal.Crystal)

Example 2 with EnergyException

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;
}
Also used : COMRestraint(ffx.potential.nonbonded.COMRestraint) CoordRestraint(ffx.potential.nonbonded.CoordRestraint) NCSRestraint(ffx.potential.nonbonded.NCSRestraint) Atom(ffx.potential.bonded.Atom) EnergyException(ffx.potential.utils.EnergyException)

Example 3 with EnergyException

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;
}
Also used : OpenMM_Vec3(simtk.openmm.OpenMM_Vec3) OpenMM_System_addConstraint(simtk.openmm.OpenMMLibrary.OpenMM_System_addConstraint) CoordRestraint(ffx.potential.nonbonded.CoordRestraint) Atom(ffx.potential.bonded.Atom) EnergyException(ffx.potential.utils.EnergyException)

Example 4 with EnergyException

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;
}
Also used : EnergyException(ffx.potential.utils.EnergyException) EnergyException(ffx.potential.utils.EnergyException)

Example 5 with EnergyException

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;
}
Also used : EnergyException(ffx.potential.utils.EnergyException) EnergyException(ffx.potential.utils.EnergyException)

Aggregations

EnergyException (ffx.potential.utils.EnergyException)18 ForceFieldString (ffx.potential.parameters.ForceField.ForceFieldString)6 CoordRestraint (ffx.potential.nonbonded.CoordRestraint)4 PotentialsFunctions (ffx.potential.utils.PotentialsFunctions)4 PotentialsUtils (ffx.potential.utils.PotentialsUtils)4 File (java.io.File)4 Atom (ffx.potential.bonded.Atom)3 COMRestraint (ffx.potential.nonbonded.COMRestraint)3 NCSRestraint (ffx.potential.nonbonded.NCSRestraint)3 Crystal (ffx.crystal.Crystal)2 ForceFieldEnergy (ffx.potential.ForceFieldEnergy)2 Respa (ffx.algorithms.integrators.Respa)1 Stochastic (ffx.algorithms.integrators.Stochastic)1 MultiResidue (ffx.potential.bonded.MultiResidue)1 Residue (ffx.potential.bonded.Residue)1 OpenMM_System_addConstraint (simtk.openmm.OpenMMLibrary.OpenMM_System_addConstraint)1 OpenMM_Vec3 (simtk.openmm.OpenMM_Vec3)1