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 *

universe = Universe(dimensions=21.75, constraint_algorithm=Shake(1e-4, 100), electrostatic_solver=PPPM(accuracy=1e-5))
H1 = Atom('H')
H2 = Atom('H', position=(0., 1.63298, 0.))
O = Atom('O', position=(0., 0.81649, 0.57736))
H_coulombic = Coulombic(atoms=[H1, H2], cutoff=10.)
O_coulombic = Coulombic(atoms=O, cutoff=10.)
water_mol = Molecule(position=(0, 0, 0),
                     velocity=(0, 0, 0),
                     atoms=[H1, H2, O],
                     interactions=[Bond((H1, O), (H2, O), constrained=True),
                                   BondAngle(H1, O, H2, constrained=True)],
                     name='water')
universe.fill(water_mol, num_density=0.03356718472021752)
O_dispersion = Dispersion(universe, [O.atom_type, O.atom_type], cutoff=10., vdw_tail_correction=True)
universe.add_force_field('SPCE')

Universe created with:
Dimensions [21.75 21.75 21.75]

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='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)
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
Total wall time: 0:00:00
Simulation created with lammps engine and settings:
temperature: 300.0 K
pressure: 101325.0 Pa
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.02056835377978693, -0.010000746338816117, 0.010784894567034709)

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

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 output is captured by PyLammps wrapper
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
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 K
pressure: 101325.0 Pa
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.5042565657111322, 0.0, -0.5042565657111322)

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.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.