Source code for MMTK.Dynamics

# This module implements MD integrators
#
# Written by Konrad Hinsen
#

"""
Molecular Dynamics integrators
"""

__docformat__ = 'restructuredtext'

from MMTK import Environment, Features, Trajectory, Units
import MMTK_dynamics
from Scientific import N

#
# Integrator base class
#
class Integrator(Trajectory.TrajectoryGenerator):

    def __init__(self, universe, options):
        Trajectory.TrajectoryGenerator.__init__(self, universe, options)

    default_options = {'first_step': 0, 'steps': 100, 'delta_t': 1.*Units.fs,
                       'background': False, 'threads': None,
                       'mpi_communicator': None, 'actions': []}

    available_data = ['configuration', 'velocities', 'gradients',
                      'energy', 'thermodynamic', 'time', 'auxiliary']

    restart_data = ['configuration', 'velocities', 'energy',
                    'thermodynamic', 'auxiliary']

    def __call__(self, options):
        raise AttributeError

#
# Velocity-Verlet integrator
#
[docs]class VelocityVerletIntegrator(Integrator): """ Velocity-Verlet molecular dynamics integrator The integrator can handle fixed atoms, distance constraints, a thermostat, and a barostat, as well as any combination. It is fully thread-safe. The integration is started by calling the integrator object. All the keyword options (see documnentation of __init__) can be specified either when creating the integrator or when calling it. The following data categories and variables are available for output: - category "time": time - category "configuration": configuration and box size (for periodic universes) - category "velocities": atomic velocities - category "gradients": energy gradients for each atom - category "energy": potential and kinetic energy, plus extended-system energy terms if a thermostat and/or barostat are used - category "thermodynamic": temperature, volume (if a barostat is used) and pressure - category "auxiliary": extended-system coordinates if a thermostat and/or barostat are used """ def __init__(self, universe, **options): """ :param universe: the universe on which the integrator acts :type universe: :class:`~MMTK.Universe.Universe` :keyword steps: the number of integration steps (default is 100) :type steps: int :keyword delta_t: the time step (default is 1 fs) :type delta_t: float :keyword actions: a list of actions to be executed periodically (default is none) :type actions: list :keyword threads: the number of threads to use in energy evaluation (default set by MMTK_ENERGY_THREADS) :type threads: int :keyword background: if True, the integration is executed as a separate thread (default: False) :type background: bool :keyword mpi_communicator: an MPI communicator object, or None, meaning no parallelization (default: None) :type mpi_communicator: Scientific.MPI.MPICommunicator """ Integrator.__init__(self, universe, options) self.features = [Features.FixedParticleFeature, Features.NoseThermostatFeature, Features.AndersenBarostatFeature, Features.DistanceConstraintsFeature] def __call__(self, **options): """ Run the integrator. The keyword options are the same as described under __init__. """ self.setCallOptions(options) used_features = Features.checkFeatures(self, self.universe) configuration = self.universe.configuration() velocities = self.universe.velocities() if velocities is None: raise ValueError("no velocities") masses = self.universe.masses() fixed = self.universe.getAtomBooleanArray('fixed') nt = self.getOption('threads') comm = self.getOption('mpi_communicator') evaluator = self.universe.energyEvaluator(threads=nt, mpi_communicator=comm) evaluator = evaluator.CEvaluator() constraints, const_distances_sq, c_blocks = \ _constraintArrays(self.universe) type = 'NVE' if Features.NoseThermostatFeature in used_features: type = 'NVT' thermostat = self.universe.environmentObjectList( Environment.NoseThermostat)[0] t_parameters = thermostat.parameters t_coordinates = thermostat.coordinates else: t_parameters = N.zeros((0,), N.Float) t_coordinates = N.zeros((2,), N.Float) if Features.AndersenBarostatFeature in used_features: if self.universe.cellVolume() is None: raise ValueError("Barostat requires finite volume universe") if type == 'NVE': type = 'NPH' else: type = 'NPT' barostat = self.universe.environmentObjectList( Environment.AndersenBarostat)[0] b_parameters = barostat.parameters b_coordinates = barostat.coordinates else: b_parameters = N.zeros((0,), N.Float) b_coordinates = N.zeros((1,), N.Float) args = (self.universe, configuration.array, velocities.array, masses.array, fixed.array, evaluator, constraints, const_distances_sq, c_blocks, t_parameters, t_coordinates, b_parameters, b_coordinates, self.getOption('delta_t'), self.getOption('first_step'), self.getOption('steps'), self.getActions(), type + ' dynamics trajectory with ' + self.optionString(['delta_t', 'steps'])) return self.run(MMTK_dynamics.integrateVV, args) # # Velocity scaling, removal of global translation/rotation #
[docs]class VelocityScaler(Trajectory.TrajectoryAction): """ Periodic velocity scaling action A VelocityScaler object is used in the action list of a VelocityVerletIntegrator. It rescales all atomic velocities by a common factor to make the temperature of the system equal to a predefined value. """ def __init__(self, temperature, window=0., first=0, last=None, skip=1): """ :param temperature: the temperature value to which the velocities should be scaled :type temperature: float :param window: the deviation from the ideal temperature that is tolerated in either direction before rescaling takes place :type window: float :param first: the number of the first step at which the action is executed :type first: int :param last: the number of the first step at which the action is no longer executed. A value of None indicates that the action should be executed indefinitely. :type last: int :param skip: the number of steps to skip between two applications of the action :type skip: int """ Trajectory.TrajectoryAction.__init__(self, first, last, skip) self.parameters = N.array([temperature, window], N.Float) self.Cfunction = MMTK_dynamics.scaleVelocities
[docs]class Heater(Trajectory.TrajectoryAction): """ Periodic heating action A Heater object us used in the action list of a VelocityVerletIntegrator. It scales the velocities to a temperature that increases over time. """ def __init__(self, temp1, temp2, gradient, first=0, last=None, skip=1): """ :param temp1: the temperature value to which the velocities should be scaled initially :type temp1: float :param temp2: the final temperature value to which the velocities should be scaled :type temp2: float :param gradient: the temperature gradient (in K/ps) :type gradient: float :param first: the number of the first step at which the action is executed :type first: int :param last: the number of the first step at which the action is no longer executed. A value of None indicates that the action should be executed indefinitely. :type last: int :param skip: the number of steps to skip between two applications of the action :type skip: int """ Trajectory.TrajectoryAction.__init__(self, first, last, skip) self.parameters = N.array([temp1, temp2, gradient], N.Float) self.Cfunction = MMTK_dynamics.heat
[docs]class BarostatReset(Trajectory.TrajectoryAction): """ Barostat reset action A BarostatReset object is used in the action list of a VelocityVerletIntegrator. It resets the barostat coordinate to zero. """ def __init__(self, first=0, last=None, skip=1): """ :param first: the number of the first step at which the action is executed :type first: int :param last: the number of the first step at which the action is no longer executed. A value of None indicates that the action should be executed indefinitely. :type last: int :param skip: the number of steps to skip between two applications of the action :type skip: int """ Trajectory.TrajectoryAction.__init__(self, first, last, skip) self.parameters = N.zeros((0,), N.Float) self.Cfunction = MMTK_dynamics.resetBarostat
[docs]class TranslationRemover(Trajectory.TrajectoryAction): """ Action that eliminates global translation A TranslationRemover object is used in the action list of a VelocityVerletIntegrator. It subtracts the total velocity of the system from each atomic velocity. """ def __init__(self, first=0, last=None, skip=1): """ :param first: the number of the first step at which the action is executed :type first: int :param last: the number of the first step at which the action is no longer executed. A value of None indicates that the action should be executed indefinitely. :type last: int :param skip: the number of steps to skip between two applications of the action :type skip: int """ Trajectory.TrajectoryAction.__init__(self, first, last, skip) self.parameters = None self.Cfunction = MMTK_dynamics.removeTranslation
[docs]class RotationRemover(Trajectory.TrajectoryAction): """ Action that eliminates global rotation A RotationRemover object is used in the action list of a VelocityVerletIntegrator. It adjusts the atomic velocities such that the total angular momentum is zero. """ def __init__(self, first=0, last=None, skip=1): """ :param first: the number of the first step at which the action is executed :type first: int :param last: the number of the first step at which the action is no longer executed. A value of None indicates that the action should be executed indefinitely. :type last: int :param skip: the number of steps to skip between two applications of the action :type skip: int """ Trajectory.TrajectoryAction.__init__(self, first, last, skip) self.parameters = None self.Cfunction = MMTK_dynamics.removeRotation # # Construct constraint arrays #
def _constraintArrays(universe): nc = universe.numberOfDistanceConstraints() constraints = N.zeros((nc, 2), N.Int) const_distances_sq = N.zeros((nc, ), N.Float) if nc > 0: i = 0 c_blocks = [0] for o in universe.objectList(): for c in o.distanceConstraintList(): constraints[i, 0] = c[0].index constraints[i, 1] = c[1].index const_distances_sq[i] = c[2]*c[2] i = i + 1 c_blocks.append(i) c_blocks = N.array(c_blocks) else: c_blocks = N.zeros((1,), N.Int) return constraints, const_distances_sq, c_blocks # # Enforce distance constraints for current configuration # def enforceConstraints(universe, configuration=None): constraints, const_distances_sq, c_blocks = _constraintArrays(universe) if len(constraints) == 0: return if configuration is None: configuration = universe.configuration() MMTK_dynamics.enforceConstraints(universe._spec, configuration.array, universe.masses().array, constraints, const_distances_sq, c_blocks) # # Project velocity vector onto the constraint surface # def projectVelocities(universe, velocities): constraints, const_distances_sq, c_blocks = _constraintArrays(universe) if len(constraints) == 0: return MMTK_dynamics.projectVelocities(universe._spec, universe.configuration().array, velocities.array, universe.masses().array, constraints, const_distances_sq, c_blocks)