Running a Simulation in MDMC

This tutorial continues from the Building a Universe tutorial, using the example universe filled with SPCE water.

The %%capture and %run commands below simply executes the Building a Universe notebook and captures the variables into this notebook. They are only valid if they are executed in the same folder as the Building a Universe notebook. Otherwise, please copy the last section of the Building a Universe to set the same state.

%%capture
# Run Building a universe notebook and hide output
%run "building-a-universe.ipynb"

Creating a Simulation object

Simulations in MDMC are run using external MD engines (e.g. LAMMPS). First create a universe (see Building a Universe - an example universe can be found at the bottom of the page). 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='lammps', time_step=1., temperature=300.,
                        pressure=101325., traj_step=10, thermostat='nose',
                        barostat='nose', t_damp=100, p_damp=1000)

LAMMPS (29 Sep 2021 - Update 3)
LAMMPS output is captured by PyLammps wrapper
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:98)
  using 1 OpenMP thread(s) per MPI task
LAMMPS (29 Sep 2021 - Update 3)
LAMMPS output is captured by PyLammps wrapper
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:98)
  using 1 OpenMP thread(s) per MPI task
Total wall time: 0:00:00
Simulation created with lammps engine and settings:
  temperature     300.0
  pressure     101325.0
  thermostat       nose
  barostat         nose
  t_damp            100
  p_damp           1000

The NPT simulation is created with a Nose-Hoover thermostat and a Berendsen barostat, with a temperature of 300 K, a pressure of 101325 Pa, and a time_step of 1.0 fs. The trajectory is recorded every 10 simulation steps. The number of steps over which the temperature and pressure are damped for the thermostat and barostats is 100 and 1000 respectively.

NVE and NVT simulations can be created in a similar manner by omitting one or both of the thermostat and barostat. If both are omitted the temperature must still be specified as this initialises the atomic velocities.

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('Velocity of first atom in MDMC universe is {}'.format(universe.atoms[0].velocity))
print('Velocity of first atom in MD engine is {}'.format(simulation.engine.lmp.atoms[0].velocity))

# In comparision, create a new simulation object where (artifcially and for demonstration purposes) one atom has non-zero velocity
velocity = (1, 0, -1)
print('\nChanging the velocity of the first atom to {}\n'.format(velocity))
universe.atoms[0].velocity = velocity
simulation_2 = Simulation(universe, engine='lammps', time_step=1., temperature=300.,
                          pressure=101325., traj_step=10, thermostat='nose',
                          barostat='nose', t_damp=100, p_damp=1000)
print('Velocity of first atom in MDMC universe is {}'.format(universe.atoms[0].velocity))
print('Velocity of first atom in MD engine is {}'.format(simulation_2.engine.lmp.atoms[0].velocity))

# 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 (-0.02117368674793425, 0.0006781080459690333, -0.02126741525190225)

Changing the velocity of the first atom to (1, 0, -1)

LAMMPS (29 Sep 2021 - Update 3)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:98)
  using 1 OpenMP thread(s) per MPI task
LAMMPS output is captured by PyLammps wrapper
LAMMPS (29 Sep 2021 - Update 3)
LAMMPS output is captured by PyLammps wrapper
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:98)
  using 1 OpenMP thread(s) per MPI task
Total wall time: 0:00:00
Some but not all atom velocities set. Atoms with non-zero velocity will be re-scaled to match target temperature, atoms with zero velocity will remain stationary.
Simulation created with lammps engine and settings:
  temperature     300.0
  pressure     101325.0
  thermostat       nose
  barostat         nose
  t_damp            100
  p_damp           1000

Velocity of first atom in MDMC universe is [ 1  0 -1] Ang / fs
Velocity of first atom in MD engine is (0.5042465603227623, 0.0, -0.5042465603227623)

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)

Trajectories are not stored during equilibration.

The simulation can be run by (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]

10.0

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.