Examples with Stochastic ffx.algorithms.integrators.Stochastic used on opensource projects

Example 1 with Stochastic

use of ffx.algorithms.integrators.Stochastic 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());
    }
}