simulation

Module for setting up and running the simulation

Classes for the simulation box, minimizer and integrator.

class MDMC.MD.simulation.ConstraintAlgorithm(accuracy: float, max_iterations: int)[source]

Class describing the algorithm and parameters which are applied to constrain BondedInteraction objects

Parameters:

  • accuracy (float) – The accuracy (tolerance) of the applied constraints

  • max_iterations (int) – The maximum number of iterations that can be used when calculating the additional force that is required to constrain the atoms to satisfy the constraints on the bonded interactions

accuracy

The accuracy (tolerance) of the applied constraints

Type:

float

property max_iterations: int

Get or set the maximum number of iterations that can be used when calculating the additional force that is required to constrain the atoms to satisfy the constraints on the bonded interactions

Returns:

The maximum number of iterations

Return type:

int

property name: str

Get the name of the class

Returns:

The name of the class

Return type:

str

class MDMC.MD.simulation.Ewald(**settings: dict)[source]

Holds the parameters that are required for the Ewald solver to be applied to both/either the electrostatic and/or dispersion interactions

Parameters:

**settings –

accuracy (float)

The relative RMS error in per-atom forces

class MDMC.MD.simulation.KSpaceSolver(**settings: dict)[source]

Class describing the k-space solver that is applied to electrostatic and/or dispersion interactions

Different MDEngine require different parameters to be specified for a k-space solver to be used. These parameters are specified in settings.

Parameters:

**settings –

accuracy (float)

The relative RMS error in per-atom forces

accuracy

The relative RMS error in per-atom forces

Type:

float

property name

Get the name of the class

Returns:

The name of the class

Return type:

str

class MDMC.MD.simulation.PPPM(**settings: dict)[source]

Holds the parameters that are required for the PPPM solver to be applied to both/either the electrostatic and/or dispersion interactions

Parameters:

**settings –

accuracy (float)

The relative RMS error in per-atom forces

class MDMC.MD.simulation.Rattle(accuracy: float, max_iterations: int)[source]

Holds the parameters which are required for the RATTLE algorithm to be applied to the constrained interactions

Parameters:

  • accuracy (float) – The accuracy (tolerance) of the applied constraints

  • max_iterations (int) – The maximum number of iterations that can be used when calculating the additional force that is required to constrain the atoms to satisfy the constraints on the bonded interactions

class MDMC.MD.simulation.Shake(accuracy: float, max_iterations: int)[source]

Holds the parameters which are required for the SHAKE algorithm to be applied to the constrained interactions

Parameters:

  • accuracy (float) – The accuracy (tolerance) of the applied constraints

  • max_iterations (int) – The maximum number of iterations that can be used when calculating the additional force that is required to constrain the atoms to satisfy the constraints on the bonded interactions

class MDMC.MD.simulation.Simulation(universe: Universe, traj_step: int, time_step: float = 1.0, engine: str = 'lammps', **settings)[source]

Sets up the molecular dynamics engine and parameters for how it should run

An ensemble is defined by whether a thermostat or barostat are present

Parameters:

  • universe (Universe) – The Universe on which the simulation is performed.

  • traj_step (int) – How many steps the simulation should take between dumping each CompactTrajectory frame. Along with time_step determines the time separation of calculated variables such as energy.

  • time_step (float, optional) – Simulation timestep in fs. Default is 1.

  • engine (str, optional) – The MDEngine used for the simulation. Default is 'lammps'.

  • **settings –

    temperature (float)

    Simulation temperature in K. Note that velocities of atoms in the MD engine are set based on the temperature. If all atoms in the universe have 0 velocity, then velocities for the MD engine will be randomly chosen from a uniform distribution and then scaled to give the correct temperature. If one or more velocities are set, then these will be scaled to give the correct temperature. In either case, only the velocities on the engine, not the universe, are affected.

    integrator (str)

    Simulation time integrator.

    lj_options (dict)

    option:value pairs for Lennard-Jones interactions.

    es_options (dict)

    option:value pairs for electrostatic interactions.

    thermostat (str)

    The name of the thermostat e.g. Nose-Hoover.

    barostat (str)

    The name of the barostat e.g. Nose-Hoover.

    pressure (float)

    Simulation pressure in Pa. This is required if a barostat is passed.

    verbose (str)

    If the output of the class instantiation should be reported, default to True.

universe

The Universe on which the simulation is performed.

Type:

Universe

traj_step

How many steps the simulation should take between dumping each CompactTrajectory frame. Along with time_step determines the time separation of calculated variables such as energy.

Type:

int

engine

A subclass of MDEngine which provides the interface to the MD library. Default is 'lammps'.

Type:

MDEngine, optional

settings

The settings passed to the Simulation. See the Parameters section for details.

Type:

dict

time_step
trajectory
auto_equilibrate(variables: list[str] = ['temp', 'pe'], eq_step: int = 10, window_size: int = 100, tolerance: float = 0.01) → Tuple[int, dict][source]

Equilibrate until the specified list of variables have stabilised. Uses the KPSS stationarity test to determine whether the variables are stationary in a given window.

Parameters:

  • variables (list[str], default ['temp', 'pe']) – The MD engine variables we use to monitor stability.

  • eq_step (int, default 10) – The number of equilibration steps between each stability check.

  • window_size (int, default 100) – The size of the rolling window of stability checks. The simulation is considered stable if it has been stationary for eq_step*window_size equilibration steps.

  • tolerance (float, 0.01) – The p-value used for the KPSS test.

Returns:

  • int – The number of equilibration steps performed.

  • vals_dict – The value of each variable per step, for graphing if desired.

minimize(n_steps: int, minimize_every: int = 10, verbose: bool = False, output_log: str = None, work_dir: str = None, **settings: dict) → None[source]

Performs an MD run intertwined with periodic structure relaxation. This way after a local minimum is found, the system is taken out of the minimum to explore a larger volume of the parameter space.

Parameters:

  • n_steps (int) – Total number of the MD run steps

  • minimize_every (int, optional) – Number of MD steps between two consecutive minimizations

  • verbose (bool, optional) – Whether to print statements when the minimization has been started and completed (including the number of minimization steps and time taken). Default is False.

  • output_log (str, optional) – Log file for the MD engine to write to. Default is None.

  • work_dir (str, optional) – Working directory for the MD engine to write to. Default is None.

  • **settings –

    etol (float)

    If the energy change between iterations is less than etol, minimization is stopped. Default depends on engine used.

    ftol (float)

    If the magnitude of the global force is less than ftol, minimization is stopped. Default depends on engine used.

    maxiter (int)

    Maximum number of iterations of a single structure relaxation procedure. Default depends on engine used.

    maxeval (int)

    Maximum number of force evaluations to perform. Default depends on engine used.

run(n_steps: int, equilibration: bool = False, verbose: bool = False, output_log: str = None, work_dir: str = None, **settings: dict) → None[source]

Runs the MD simulation for the specified number of steps. Trajectories for the simulation are only saved when equilibration is False. Additionally running equilibration for an NVE system (neither barostat nor thermostat set) will temporarily apply a Berendsen thermostat (it is removed from the simulation after the run is completed).

Parameters:

  • n_steps (int) – Number of simulation steps to run

  • equilibration (bool, optional) – If the run is for equilibration (True) or production (False). Default is False.

  • verbose (bool, optional) – Whether to print statements upon starting and completing the run. Default is False.

  • output_log (str, optional) – Log file for the MD engine to write to. Default is None.

  • work_dir (str, optional) – Working directory for the MD engine to write to. Default is None.

property time_step: float

Get or set the simulation time step in fs

Returns:

Simulation time step in fs

Return type:

float

property traj_step: int

Get or set the number of simulation steps between saving the CompactTrajectory

Returns:

Number of simulation steps that elapse between the CompactTrajectory being stored

Return type:

int

property trajectory: CompactTrajectory | None

The CompactTrajectory produced by the most recent production run of the Simulation.

Returns:

Most recent production run CompactTrajectory, or None if no production run has been performed

Return type:

CompactTrajectory

class MDMC.MD.simulation.Universe(dimensions, force_field=None, structures=None, **settings)[source]

Class where the configuration and topology are defined

Parameters:

  • dimensions (numpy.ndarray, list, float) – Dimensions of the Universe, in units of Ang. A float can be used for a cubic universe.

  • force_fields (ForceField, optional) – A force field to apply to the Universe. The force fields available are: SPC, TIP3P, OPLSAA, TIP3PFB, SPCE. Default is None.

  • structures (list, optional) – Structure objects contained in the Universe. Default is None.

  • **settings –

    kspace_solver (KSpaceSolver)

    The k-space solver to be used for both electrostatic and dispersive interactions. If this is passed then no electrostatic_solver or dispersive_solver may be passed.

    electrostatic_solver (KSpaceSolver)

    The k-space solver to be used for electrostatic interactions.

    dispersive_solver (KSpaceSolver)

    The k-space solver to be used for dispersive interactions.

    constraint_algorithm (ConstraintAlgorithm)

    The constraint algorithm which will be applied to constrained BondedInteractions.

    verbose (str)

    If the output of the class instantiation should be reported, default to True.

dimensions

Dimensions of the Universe in units of Ang.

Type:

numpy.ndarray, list, float

configuration

Stores the content, i.e. configuration of atoms etc within the universe

Type:

Configuration

force_fields

Force field applied to apply to the Universe

Type:

ForceField or None

kspace_solver

The k-space solver to be used for both electrostatic and dispersive interactions.

Type:

KSpaceSolver

electrostatic_solver

The k-space solver to be used for electrostatic interactions.

Type:

KSpaceSolver

dispersive_solver

The k-space solver to be used for dispersive interactions.

Type:

KSpaceSolver

constraint_algorithm

The constraint algorithm which will be applied to constrained BondedInteractions.

Type:

ConstraintAlgorithm

interactions
bonded_interactions
nonbonded_interactions
bonded_interaction_pairs
n_bonded
n_nonbonded
n_interactions
parameters
volume
element_list
element_dict
element_lookup
atoms
n_atoms
molecule_list
n_molecules
structure_list
top_level_structure_list
equivalent_top_level_structures_dict
force_fields
atom_types
atom_type_interactions
density
solvent_density
nbis_by_atom_type_pairs
add_bonded_interaction_pairs(*bonded_interaction_pairs: tuple[Interaction, Atom]) → None[source]

Adds one or more interaction pairs to the Universe

Parameters:

*bonded_interaction_pairs – one or more (interaction, atoms) pairs, where atoms is a tuple of all atoms for that specific bonded interaction

add_force_field(force_field: str, *interactions: Interaction, **settings: dict) → None[source]

Adds a force field to the specified interactions. If no interactions are passed, the force field is applied to all interactions in the Universe.

Parameters:

  • force_field (str) – The ForceField to parameterize interactions (if provided), or all the interactions in the Universe. The available ForceField are: SPC, TIP3P, OPLSAA, TIP3PFB, SPCE

  • *interactions – Interaction objects to parameterize with the ForceField

  • **settings –

    add_dispersions (bool or list of Atoms)

    If True, a Dispersion interaction will be added to all atoms in the Universe. If a list of Atom objects is provided, the Dispersion will be added to these instead. Any added Dispersion interactions (and any previously defined) will then be parametrized by the ForceField. The Dispersion interactions added will only be like-like. By default, no Dispersion interactions are added.

add_nonbonded_interaction(*nonbonded_interactions: tuple[Interaction, Atom]) → None[source]

Adds one or more nonbonded interactions to the Universe

Parameters:

*nonbonded_interactions –

Nonbonded interactions to be added to the Universe. Can take any number of non-bonded interactions:

Dispersion():

either Lennard-Jones or Buckingham dispersion

Coulombic():

normal or modified Coulomb interaction

with appropriate parameters for the interaction. See http://mdmcproject.org/tutorials/building-a-universe.html?highlight=interaction#Create-non-bonded-interactions for more details on non-bonded interactions.

add_structure(structure: Structure | int, force_field: str = None, center: bool = False) → None[source]

Adds a single Structure to the Universe, with optional ForceField applying only to that Structure

Parameters:

  • structure (Structure or int) – The Structure (or its ID as an int) to be added to the Universe

  • force_field (str, optional) – The force field to be applied to the structure. The available ForceField are: SPC, TIP3P, OPLSAA, TIP3PFB, SPCE

  • center (bool, optional) – Whether to center structure within the Universe as it is added

property atom_type_interactions: dict[int, list[Interaction]]

Get the atom types and the interactions for each atom type

Returns:

atom_type:interactions pairs where atom_type is a int specifying the atom type and interactions is a list of Interaction objects acting on that atom_type

Return type:

dict

property atom_types: list[Atom]

Get the atom types of atoms in the Universe

Returns:

The atom types in the Universe

Return type:

list

property atoms: list[Atom]

Get a list of the atoms in the Universe

Returns:

The atoms in the Universe

Return type:

list

property bonded_interaction_pairs: list

Get the bonded interactions and the atoms they apply to

Returns:

The (interaction, atoms) pairs in the Universe, where atoms is a tuple of all Atom objects for that specific interaction

Return type:

list

Example

For an O Atom with two bonds, one to H1 and one to H2:

>>> print(O.bonded_interaction_pairs)
[(Bond, (H1, O)), (Bond, (H2, O))]
property bonded_interactions: list

Get the bonded interactions in the Universe

Returns:

The bonded interactions in the Universe

Return type:

list

property element_dict: dict[str, Atom]

Get the elements in the Universe and example Atom objects for each element

This is required for MD engines which assign the same potential parameters for all identical element types.

Returns:

element:atom pairs, where atom is a single Atom of the specified element.

Return type:

dict

property element_list: list[str]

The elements of the atoms in the Universe

Returns:

The elements in the Universe

Return type:

list

property element_lookup: dict[str, Atom]

Get the elements by ID in the Universe

This is required for MD engines which assign the same potential parameters for all identical element names.

Returns:

element:atom pairs, where atom is a single Atom of the specified element.

Return type:

dict

property equivalent_top_level_structures_dict: dict[Structure, int]

Get a dict of equivalent top level Structure objects that exist in the Universe as keys, and the number of Structure that are equivalent to the key as a value. This does not include any Structure that are a subunit of another structure belonging to the Universe.

Returns:

The top level Structure objects in the Universe as keys, and the number of each as a value

Return type:

dict

fill(structures: Structure, force_field: str = None, num_density: float = None, num_struc_units: int = None) → None[source]

A liquid-like filling of the Universe independent of existing atoms

Adds copies of structures to existing configuration until Universe is full. As exclusion region is defined by the size of a bounding sphere, this method is most suitable for atoms or molecules with approximately equal dimensions.

Note

CURRENT APPROACH RESULTS IN NUMBER DENSITY DIFFERENT TO WHAT IS SPECIFIED DEPENDING ON HOW CLOSE CUBE ROOT OF N_MOLECULES IS TO AN int.

Note

CURRENT IMPLEMENTATION SHOULD NOT BE USED WITH NON-CUBIC UNIVERSES AS THE DENSITY MAY OR MAY NOT BE ISOTROPIC DEPENDING ON THE DIMENSIONS AND NUMBER OF UNITS.

Parameters:

  • structures (Structure or int) – The Structure with which to fill the Universe

  • force_field (str) – Applies a ForceField to the Universe. The available ForceField are: SPC, TIP3P, OPLSAA, TIP3PFB, SPCE

  • num_density (float) – Non-negative float specifying the number density of the Structure, in units of Structure / Ang ^ -3

  • num_struc_units (int) – Non-negative int specifying the number of passed Structure objects that the universe should be filled with, regardless of Universe.dimensions.

Raises:

  • ValueError – If both num_density and num_struc_units are passed

  • ValueError – If neither num_density or num_struc_units are passed

property force_fields: None

Get the ForceField acting on the Universe

The available force fields are: SPC, TIP3P, OPLSAA, TIP3PFB, SPCE

Returns:

ForceField that applies to the Universe

Return type:

list

property interactions: list

Get the interactions in the Universe

Returns:

The interactions in the Universe

Return type:

list

property molecule_list: list[Molecule]

Get a list of the Molecule objects in the Universe

Returns:

The Molecule objects in the Universe

Return type:

list

property n_atoms: int

Get the number of atoms in the Universe

Returns:

The number of atoms in the Universe

Return type:

int

property n_bonded: int

Get the number of bonded interactions in the Universe

Returns:

The number of bonded interactions in the Universe

Return type:

int

property n_interactions: int

Get the number of interactions in the Universe

Returns:

The number of interactions in the Universe

Return type:

int

property n_molecules: int

Get the number of molecules in the Universe

Returns:

The number of molecules in the Universe

Return type:

int

property n_nonbonded: int

Get the number of nonbonded interactions in the Universe

Returns:

The number of nonbonded interactions in the Universe

Return type:

int

property nbis_by_atom_type_pairs: dict[tuple[int], list[Interaction]]

Generates a dict of all nonbonded interactions possessed by all combinations of atom_type pairs in the Universe.

Returns:

A dict of {pair: interactions} pairs, where:

  • pair is a tuple for a pair of atom_types in the Universe

  • interactions is a list of nonbonded interactions that exist for this pair of atom_types

Any Dispersions in interactions are ones that exist explicity for this pair, whereas any Coulombics in interactions are ones that exist for either of the atom_types in pair.

Return type:

dict

property nonbonded_interactions: list

Get the nonbonded interactions in the Universe

Returns:

The nonbonded interactions in the Universe

Return type:

list

property parameters: Parameters

Get the parameters of the interactions that exist within the Universe

Returns:

The Parameters objects defined within Universe

Return type:

Parameters

solvate(density: float, tolerance: float = 1.0, solvent: str = 'SPCE', max_iterations: int = 100, **settings: dict) → None[source]

Fills the Universe with solvent molecules according to pre-defined coordinates.

Parameters:

  • density (float) – The desired density of the Solvent that solvates the Universe, in units of amu Ang ^ -3

  • tolerance (float, optional) – The +/- percentage tolerance of the density to be achieved. The default is 1 %. Tolerances of less than 1 % are at risk of not converging.

  • solvent (str, optional) – A str specifying an inbuilt Solvent from the following: SPC, TIP3P, TIP3PFB, SPCE. The default is ‘SPCE’.

  • max_iterations (int, optional) – The maximum number of times to try to solvate the universe to within the required density before stopping. Defaults to 100.

  • **settings –

    constraint_algorithm (ConstraintAlgorithm)

    A ConstraintAlgorithm which is applied to the Universe. If an inbuilt Solvent is selected (e.g. ‘SPCE’) and constraint_algorithm is not passed, the ConstraintAlgorithm will default to Shake(1e-4, 100).

Raises:

ValueError – If the Universe has already been solvated with a different density.

property structure_list: None

Get a list of all Structure objects that exist in the Universe. This includes all Structure that are a subunit of another structure belonging to the Universe.

Returns:

The Structure objects in the Universe

Return type:

list

property top_level_structure_list: list[Structure]

Get a list of the top level Structure objects that exist in the Universe. This does not include any Structure that are a subunit of another structure belonging to the Universe.

Returns:

The top level Structure objects in the Universe

Return type:

list