Selecting the algorithms¶

Integrator: running the simulation¶

An integrator is an algorithm responsible for updating the positions and the velocities of the current frame of the simulation.

Verlet integrators¶

Verlet integrators are based on Taylor expensions of Newton’s second law. They provide a simple way to integrate the movement, and conserve the energy if a sufficently small timestep is used. Assuming the absence of barostat and thermostat, they provide a NVE integration.

type Verlet¶: The Verlet algorithm is described here for example. The main constructor for this integrator is Verlet(timestep), where timestep is the timestep in femtosecond.

type VelocityVerlet¶: The velocity-Verlet algorithm is descibed extensively in the literature, for example in this webpages. The main constructor for this integrator is VelocityVerlet(timestep), where timestep is the integration timestep in femtosecond. This is the default integration algorithm in Jumos.

Checking the simulation consistency¶

Molecular dynamic is usually a garbage in, garbage out set of algorithms. The numeric and physical issues are not caught by the algorithm themselves, and the physical (and chemical) consistency of the simulation should be checked often.

In Jumos, this is achieved by the Check algorithms, which are presented in this section. Checking algorithms can be added to a simulation by using the add_check() function.

Existing checks¶

type GlobalVelocityIsNull¶: This algorithm checks if the global velocity (the total moment of inertia) is null for the current simulation. The absolute tolerance is \(10^{-5}\ A/fs\).

type TotalForceIsNull¶: This algorithm checks if the sum of the forces is null for the current simulation. The absolute tolerance is \(10^{-5}\ uma \cdot A/fs^2\).

type AllPositionsAreDefined¶

This algorithm checks is all the positions and all the velocities are defined numbers, i.e. if all the values are not infinity or the NaN (not a number) values.

This algorithm is used by default by all the molecular dynamic simulation.

Controlling the simulation¶

While running a simulation, we often want to have control over some simulation parameters: the temperature, the pressure, … This is the goal of the Control algorithms.

Such algorithms are subtypes of BaseControl, and can be added to a simulation using the add_control() function:

Controlling the temperature: Thermostats¶

Various algorithms are available to control the temperature of a simulation and perform pseudo NVT simulations. The following thermostating algorithms are currently implemented:

type VelocityRescaleThermostat¶

The velocity rescale algorithm controls the temperature by rescaling all the velocities when the temperature differs exceedingly from the desired temperature.

The constructor takes two parameters: the desired temperature and a tolerance interval. If the absolute difference between the current temperature and the desired temperature is larger than the tolerance, this algorithm rescales the velocities.

sim = MolecularDynamic(2.0)

# This sets the temperature to 300K, with a tolerance of 50K
thermostat = VelocityRescaleThermostat(300, 50)

add_control(sim, thermostat)

type BerendsenThermostat¶

The berendsen thermostat sets the simulation temperature by exponentially relaxing to a desired temperature. A more complete description of this algorithm can be found in the original article [1].

The constructor takes as parameters the desired temperature, and the coupling parameter, expressed in simulation timestep units. A coupling parameter of 100, will give a coupling time of \(150\ fs\) if the simulation timestep is \(1.5\ fs\), and a coupling time of \(200\ fs\) if the timestep is \(2.0\ fs\).

BerendsenThermostat(T[, coupling])

Creates a Berendsen thermostat at the temperature T with a coupling parameter of coupling. The default values for coupling is \(100\).

sim = MolecularDynamic(2.0)

# This sets the temperature to 300K
thermostat = BerendsenThermostat(300)

add_control(sim, thermostat)

[1]	H.J.C. Berendsen, et al. J. Chem Phys 81, 3684 (1984); doi: 10.1063/1.448118

Controlling the pressure: Barostats¶

type BerendsenBarostat¶: TODO

Other controls¶

type WrapParticles¶

This control wraps the positions of all the particles inside the unit cell.

This control is present by default in the molecular dynamic simulations.

Functions for algorithms selection¶

The six following functions are used to to select specific algorithms for the simulation. They allow to add and change all the algorithms, even in the middle of a run.

set_integrator(sim, integrator)¶

Sets the simulation integrator to integrator.

Usage example:

# Creates the integrator directly in the function
set_integrator(sim, Verlet(2.5))

# Binds the integrator to a variable if you want to change a parameter
integrator = Verlet(0.5)
set_integrator(sim, integrator)
run!(sim, 300)   # Run with a 0.5 fs timestep
integrator.timestep = 1.5
run!(sim, 3000)  # Run with a 1.5 fs timestep

set_forces_computation(sim, forces_computer)¶: Sets the simulation algorithm for forces computation to forces_computer.

add_check(sim, check)¶

Adds a check to the simulation check list and issues a warning if the check is already present.

Usage example:

# Note the parentheses, needed to instanciate the new check.
add_check(sim, AllPositionsAreDefined())

add_control(sim, control)¶

Adds a control algorithm to the simulation list. If the control algorithm is already present, a warning is issued.

Usage example:

add_control(sim, RescaleVelocities(300, 100))

add_compute(sim, compute)¶

Adds a compute algorithm to the simulation list. If the algorithm is already present, a warning is issued.

Usage example:

# Note the parentheses, needed to instanciate the new compute algorithm.
add_compute(sim, EnergyCompute())

add_output(sim, output)¶

Adds an output algorithm to the simulation list. If the algorithm is already present, a warning is issued.

Usage example:

add_output(sim, TrajectoryOutput("mytraj.xyz"))