Running a Simulation in MDMC

We first create an example universe, filled with water with an SPCE forcefield. If you’d like to learn more about any of this, please read the following:

• Creating Atomic Configurations

• Defining Molecule Interactions

• Applying a Forcefield

from MDMC.MD import *
from MDMC.MD.force_fields.three_site_water import ThreeSiteWater, add_three_site_water_ff

universe = Universe(dimensions=24.83602653)
universe.fill(ThreeSiteWater(), num_density=0.03356718472021752)
add_three_site_water_ff(universe, cutoff=10.0, ewald=1e-4)

Universe created with:
Dimensions [24.84 24.84 24.84]

Creating a Simulation object

Simulations in MDMC are run using external MD engines (e.g. LAMMPS). First create a universe, as above. This universe object must then be passed when creating a Simulation object, along with simulation properties:

# Import the Simulation class
from MDMC.MD import Simulation

# Create an NPT simulation
simulation = Simulation(universe, engine='openmm', time_step=1., temperature=300., traj_step=10)

Simulation created with openmm engine and settings:
temperature: 300.0 K

MDMC allows detailed control of the atomic velocities when creating the Universe. In the case where no velocities were provided (i.e. all Atom objects have the default velocity of 0) then the starting velocities are determined by the MD engine (randomly chosen from a uniform distribution, and then scaled so that the velocities are consistent with the temperature provided to the Simulation). If some or all of the atoms have been set with MDMC then these velocities will be scaled to the correct temperature by the MD engine. In both cases, only the velocities of atoms within the MD engine are affected, and the state of the original Universe is unchanged:

# Check the simulation object, which was created from an input Universe where all atom velocities were equal to zero
print(f'Velocity of first atom in MDMC universe is {universe.atoms[0].velocity}')
state = simulation.engine.openmm_simulation.context.getState(velocities=True)
print(f'Velocity of first atom in MD engine is {state.getVelocities()[0]}')

# In comparision, create a new simulation object where (artifcially and for demonstration purposes) one atom has non-zero velocity
velocity = (1, 0, -1)
print(f'Changing the velocity of the first atom to {velocity}')
universe.atoms[0].velocity = velocity
simulation_2 = Simulation(universe, engine='openmm', time_step=1., temperature=300.,
                          pressure=101325., traj_step=10, thermostat='nose',
                          barostat='nose', t_damp=100, p_damp=1000, openmm_platform="CPU")
print(f'Velocity of first atom in MDMC universe is {universe.atoms[0].velocity}')
state = simulation_2.engine.openmm_simulation.context.getState(velocities=True)
print(f'Velocity of first atom in MD engine is {state.getVelocities()[0]}')

# Reset atom velocity back to zero
universe.atoms[0].velocity = (0, 0, 0)

Velocity of first atom in MDMC universe is [0. 0. 0.] Ang / fs
Velocity of first atom in MD engine is Vec3(x=-0.846017370436946, y=1.1498262718456465, z=-1.7074638472513446) nm/ps
Changing the velocity of the first atom to (1, 0, -1)
Simulation created with openmm engine and settings:
temperature: 300.0 K
pressure: 101325.0 Pa
thermostat: nose
barostat: nose
t_damp: 100
p_damp: 1000
openmm_platform: CPU


Velocity of first atom in MDMC universe is [ 1  0 -1] Ang / fs
Velocity of first atom in MD engine is Vec3(x=0.1, y=0.0, z=-0.1) nm/ps

Energy minimisation and running a simulation

The universe energy can be minimised by:

# Minimise the system during a 100-step MD run, minimising every 10 steps
simulation.minimize(100, minimize_every=10)

The simulation can be equilibrated by (this will take ~30s):

simulation.run(1000, equilibration=True)

During the equilibration phase, statistics and trajectories are not captured. If you have not explicitly specified a thermostat or barostat in your calculation, a Berendsen thermostat will be used for the duration of the equilibration. This Berendsen thermostat will not be used during the production run.

The simulation can be run like so. (this will take ~60s):

# Run the simulation for 2000 steps
simulation.run(2000)

Trajectory

An MDMC CompactTrajectory object can be created following a simulation run using:

trajectory = simulation.trajectory

The times of all of the steps of trajectory can be accessed with the trajectory.times attribute:

trajectory.times[1] - trajectory.times[0]

np.float64(10.000000000000002)

Variation between MD engines

In theory, all simulations on a MDMC Universe should be able to use any MD engine, although in practice this is limited by whether a particular MD engine supports a specific feature and if it has been implemented in the MD engine interface. If a feature is not supported by or implemented for a specific MD engine, MDMC will raise a NotImplementedError.