Building an MDMC Universe

An MDMC simulation requires a configuration and a topology defined within a Universe object. Although MDMC is somewhat flexible with regards to the order in which these are created, a suggested approach is as follows:

Restriction on creating a Universe object: MDMC currently only supports orthorhombic boxed Universes; the dimensions must be specified when creating (initialising) the Universe object:

# Import the Universe class
from MDMC.MD import Universe
# Initialise a Universe with dimensions in Ang
universe = Universe([10.0, 15.0, 20.0])

Supported DL_POLY version 5.0
Universe created with:
  Dimensions       [10.0, 15.0, 20.0]
  Force field                    None
  Number of atoms                   0

To create a cubic Universe, only a single float needs to be specified for the dimensions:

universe = Universe(10.0)

Universe created with:
  Dimensions       [10.0, 10.0, 10.0]
  Force field                    None
  Number of atoms                   0

Create an atomic configuration

Configurations can either be specified by the user or read from a CIF file.

Create an atom

Each atom is specified using an Atom object, which possesses (amongst other things) a position, velocity, element, mass, charge, and atom_type. At a minimum the element must be specified when creating an Atom object:

# Import the Atom class
from MDMC.MD import Atom
# Create a hydrogen atom
H1 = Atom('H')

This will create a Hydrogen atom with mass determined from an elemental lookup table, position set to the origin, velocity set to (0., 0., 0.), no charge and no atom_type. Note that if set, the velocity of atoms in the MD engine will be scaled when creating a Simulation in order to ensure the temperature is accurate. Otherwise, if the velocities of all Atom objects in a Simulation are 0, then the velocities of the atoms in the MD engine will be randomly chosen in order to provide an accurate temperature. For more details see Running a Simulation.

Atoms can also be created by copying another atom and passing the position of the new atom:

H2 = H1.copy(position=[1., 1., 1.])

The copied atom will have identical properties to the original attribute, except with a different ID (which is unique for all structural units) and a different position. This includes interactions, which will apply to the copied atom in the same way as the original atom e.g. if H1 is bonded to an atom O, H2 will also be bonded to atom O. An example of this is shown in the section ‘Create bonded interactions’.

Atom type

Each atom has an atom_type, which is used for applying non-bonded interactions between atoms; all atoms of the same atom_type will have the same non-bonded interactions applying to them. The atom_type can be specified when the atom is created:

O1 = Atom('O', position=[8.0, 8.0, 8.0], atom_type=1)
# and to add the created atom to universe do
universe.add_structure(O1)

Alternatively, the atom_type can be inferred by MDMC. This occurs when an atom with no defined atom_type is added to a Universe:

O2 = Atom('O', position=[9.0, 9.0, 9.0])
universe.add_structure(O2)

Assuming no other atoms have been added to the universe, this sets O2.atom_type equal to 1. MDMC infers the atom_type of each atom based on its element and its interactions: all atoms with the same element and the same interactions when they are added to the universe will have the same atom_type. Once an atom’s atom_type is set, it is immutable i.e. it cannot be changed. It is recommended that either all atom types are specified when atoms are created or none are (i.e. either MDMC is left to infer and assign all atom types or no atom types. Only specifying some atom types could result in unexpected behaviour; if it is absolutely necessary then ensure that all atoms with specified atom types are added to the universe first).

To see what atoms have been added to the universe:

universe.structure_list

[<Atom
  {name: 'O',
   ID: 3,
   element: 'O',
   position: UnitNDArray([8., 8., 8.]) Ang,
   velocity: UnitNDArray([0., 0., 0.]) Ang / fs}>,
 <Atom
  {name: 'O',
   ID: 4,
   element: 'O',
   position: UnitNDArray([9., 9., 9.]) Ang,
   velocity: UnitNDArray([0., 0., 0.]) Ang / fs}>]

You can see other attributes of the universe (or any other MDMC containers) with:

help(universe)

Help on Universe in module MDMC.MD.simulation object:

class Universe(MDMC.MD.container.AtomContainer)
 |  Universe(dimensions, force_field=None, structures=None, **settings)
 |
 |  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:
 |      TIP3P, SPC, TIP3PFB, OPLSAA, 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.
 |
 |  Attributes
 |  ----------
 |  dimensions : numpy.ndarray, list, float
 |      Dimensions of the ``Universe`` in units of ``Ang``.
 |  configuration : Configuration
 |      Stores the content, i.e. configuration of atoms etc within the universe
 |  force_fields : ForceField or None
 |      Force field applied to apply to the Universe
 |  kspace_solver : KSpaceSolver
 |      The k-space solver to be used for both electrostatic and dispersive
 |      interactions.
 |  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``.
 |  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
 |
 |  Method resolution order:
 |      Universe
 |      MDMC.MD.container.AtomContainer
 |      abc.ABC
 |      builtins.object
 |
 |  Methods defined here:
 |
 |  __eq__(self, other) -> bool
 |      Return self==value.
 |
 |  __init__(self, dimensions, force_field=None, structures=None, **settings)
 |      Initialize self.  See help(type(self)) for accurate signature.
 |
 |  __repr__(self)
 |
 |  __str__(self) -> str
 |      Return str(self).
 |
 |  add_bonded_interaction_pairs(self, *bonded_interaction_pairs: 'tuple[Interaction, Atom]') -> None
 |      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(self, force_field: str, *interactions: 'Interaction', **settings: dict) -> None
 |      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:
 |          TIP3P, SPC, TIP3PFB, OPLSAA, 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(self, *nonbonded_interactions: 'tuple[Interaction, Atom]') -> None
 |      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(self, structure: Union[MDMC.MD.structures.Structure, int], force_field: str = None, center: bool = False) -> None
 |      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:
 |          TIP3P, SPC, TIP3PFB, OPLSAA, SPCE
 |      center : bool, optional
 |          Whether to center `structure` within the Universe as it is
 |          added
 |
 |  fill(self, structures: MDMC.MD.structures.Structure, force_field: str = None, num_density: float = None, num_struc_units: int = None) -> None
 |      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:
 |          TIP3P, SPC, TIP3PFB, OPLSAA, 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
 |
 |  solvate(self, density: float, tolerance: float = 1.0, solvent: str = 'SPCE', max_iterations: int = 100, **settings: dict) -> None
 |      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:
 |          TIP3P, TIP3PFB, SPCE, SPC.
 |          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.
 |
 |  ----------------------------------------------------------------------
 |  Readonly properties defined here:
 |
 |  atom_type_interactions
 |      Get the atom types and the interactions for each atom type
 |
 |      Returns
 |      -------
 |      dict
 |          ``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``
 |
 |  atom_types
 |      Get the atom types of atoms in the ``Universe``
 |
 |      Returns
 |      -------
 |      list
 |          The atom types in the ``Universe``
 |
 |  atoms
 |      Get a list of the atoms in the ``Universe``
 |
 |      Returns
 |      -------
 |      list
 |          The atoms in the ``Universe``
 |
 |  bonded_interaction_pairs
 |      Get the bonded interactions and the atoms they apply to
 |
 |      Returns
 |      -------
 |      list
 |          The (``interaction``, ``atoms``) pairs in the ``Universe``, where
 |          ``atoms`` is a `tuple` of all ``Atom`` objects for that specific
 |          ``interaction``
 |
 |      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))]
 |
 |  bonded_interactions
 |      Get the bonded interactions in the ``Universe``
 |
 |      Returns
 |      -------
 |      list
 |          The bonded interactions in the ``Universe``
 |
 |  density
 |
 |  element_dict
 |      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
 |      -------
 |      dict
 |          ``element:atom pairs``, where ``atom`` is a single ``Atom`` of the
 |          specified ``element``.
 |
 |  element_list
 |      The elements of the atoms in the ``Universe``
 |
 |      Returns
 |      -------
 |      list
 |          The elements in the ``Universe``
 |
 |  element_lookup
 |      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
 |      -------
 |      dict
 |          ``element:atom pairs``, where ``atom`` is a single ``Atom`` of the
 |          specified ``element``.
 |
 |  equivalent_top_level_structures_dict
 |      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
 |      -------
 |      dict
 |          The top level ``Structure`` objects in the ``Universe`` as
 |          keys, and the number of each as a value
 |
 |  force_fields
 |      Get the ``ForceField`` acting on the ``Universe``
 |
 |      The available force fields are:
 |      TIP3P, SPC, TIP3PFB, OPLSAA, SPCE
 |
 |      Returns
 |      -------
 |      list
 |          ``ForceField`` that applies to the ``Universe``
 |
 |  interactions
 |      Get the interactions in the ``Universe``
 |
 |      Returns
 |      -------
 |      list
 |          The interactions in the ``Universe``
 |
 |  molecule_list
 |      Get a list of the ``Molecule`` objects in the ``Universe``
 |
 |      Returns
 |      -------
 |      list
 |          The ``Molecule`` objects in the ``Universe``
 |
 |  n_atoms
 |      Get the number of atoms in the ``Universe``
 |
 |      Returns
 |      -------
 |      int
 |          The number of atoms in the ``Universe``
 |
 |  n_bonded
 |      Get the number of bonded interactions in the ``Universe``
 |
 |      Returns
 |      -------
 |      int
 |          The number of bonded interactions in the ``Universe``
 |
 |  n_interactions
 |      Get the number of interactions in the ``Universe``
 |
 |      Returns
 |      -------
 |      int
 |          The number of interactions in the ``Universe``
 |
 |  n_molecules
 |      Get the number of molecules in the ``Universe``
 |
 |      Returns
 |      -------
 |      int
 |          The number of molecules in the ``Universe``
 |
 |  n_nonbonded
 |      Get the number of nonbonded interactions in the ``Universe``
 |
 |      Returns
 |      -------
 |      int
 |          The number of nonbonded interactions in the ``Universe``
 |
 |  nbis_by_atom_type_pairs
 |      Generates a `dict` of all nonbonded interactions possessed by all
 |      combinations of ``atom_type`` pairs in the ``Universe``.
 |
 |      Returns
 |      -------
 |      dict
 |          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``.
 |
 |  nonbonded_interactions
 |      Get the nonbonded interactions in the ``Universe``
 |
 |      Returns
 |      -------
 |      list
 |          The nonbonded interactions in the ``Universe``
 |
 |  parameters
 |      Get the parameters of the interactions that exist within the
 |      ``Universe``
 |
 |      Returns
 |      -------
 |      Parameters
 |          The ``Parameters`` objects defined within ``Universe``
 |
 |  solvent_density
 |
 |  structure_list
 |      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
 |      -------
 |      list
 |          The ``Structure`` objects in the ``Universe``
 |
 |  top_level_structure_list
 |      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
 |      -------
 |      list
 |          The top level ``Structure`` objects in the ``Universe``
 |
 |  volume
 |
 |  ----------------------------------------------------------------------
 |  Data descriptors defined here:
 |
 |  dimensions
 |
 |  ----------------------------------------------------------------------
 |  Data and other attributes defined here:
 |
 |  __abstractmethods__ = frozenset()
 |
 |  __hash__ = None
 |
 |  ----------------------------------------------------------------------
 |  Methods inherited from MDMC.MD.container.AtomContainer:
 |
 |  __getitem__(self, index: Union[int, slice]) -> Union[ForwardRef('Atom'), ForwardRef('list[Atom]')]
 |      Returns
 |      -------
 |      Atom, list
 |          The atom (or atoms) for the specified index (or slice)
 |
 |  ----------------------------------------------------------------------
 |  Data descriptors inherited from MDMC.MD.container.AtomContainer:
 |
 |  __dict__
 |      dictionary for instance variables (if defined)
 |
 |  __weakref__
 |      list of weak references to the object (if defined)

The same applies for methods, for example:

help(universe.add_structure)

Help on method add_structure in module MDMC.MD.simulation:

add_structure(structure: Union[MDMC.MD.structures.Structure, int], force_field: str = None, center: bool = False) -> None method of MDMC.MD.simulation.Universe instance
    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:
        TIP3P, SPC, TIP3PFB, OPLSAA, SPCE
    center : bool, optional
        Whether to center `structure` within the Universe as it is
        added

Create bonded interactions

There are three bonded interaction types within MDMC: Bond, BondAngle and Dihedral. Each Interaction must have an InteractionFunction which describes the Interaction.

Bond

Below this is demonstrated where a Bond with a HarmonicPotential function is created.

# Import Bond and HarmonicPotential
from MDMC.MD import Bond
from MDMC.MD import HarmonicPotential
# Create a Bond with a HarmonicPotential
# The first argument in the HarmonicPotential is the equilibrium state and the second is the potential strength
HH_bond = Bond(H1, H2, function=HarmonicPotential(1., 100., interaction_type='bond'))

To see which units are supported:

from MDMC.common.units import SYSTEM
SYSTEM

{'LENGTH': 'Ang',
 'TIME': 'fs',
 'MASS': 'amu',
 'CHARGE': 'e',
 'ANGLE': 'deg',
 'TEMPERATURE': 'K',
 'ENERGY': 'kJ / mol',
 'FORCE': 'kJ / Ang mol',
 'PRESSURE': 'Pa',
 'ENERGY_TRANSFER': 'meV',
 'ARBITRARY': 'arb'}

BondAngle

A BondAngle is created in the same manner except it requires a minimum of three atoms. The central atom should be the 2nd atom. For example, a water molecule would have:

from MDMC.MD import BondAngle

H1 = Atom('H')
H2 = Atom('H', position=[0., 1.63298, 0.])
O = Atom('O', position=[0., 0.81649, 0.57736])
HOH_angle = BondAngle(H1, O, H2)

# The following is equivalent
HOH_angle = BondAngle(H2, O, H1)

Currently HarmonicPotential is the only InteractionFunction that can be applied to either Bond or BondAngle.

All Bond and BondAngle can be constrained; this can either be set when creating the Bond (or BondAngle) or afterwards:

HH_bond.constrained = True
HOH_angle = BondAngle(H1, O, H2, constrained=True)

For a constraint to be applied during a simulation, the Universe must have a ConstraintAlgorithm.

DihedralAngle

A DihedralAngle is also created in the same manner, except it requires four atoms. A DihedralAngle can be proper or improper, as specified by DihedralAngle.improper. By default a DihedralAngle is proper.

The atoms in a proper DihedralAngle must be specified so that the 2nd and 3rd atoms are the central two atoms:

from MDMC.MD import DihedralAngle

C1 = Atom('C', position=[5.12033922, 4.63847287, 4.94075943])
N = Atom('N', position=[5.12991894, 3.78609704, 3.79996577])
C2 = Atom('C', position=[4.91462725, 4.08992816, 2.5091264])
O = Atom('O', position=[4.67405373, 5.24130678, 2.1180462])
proper = DihedralAngle(atoms=[C1, N, C2, O])

# For a proper DihedralAngle, the equivalent atom order is:
proper = DihedralAngle(atoms=[O, C2, N, C1])

The atoms in an improper DihedralAngle must be specific so that the 1st atom is the central atom.

N = Atom('N', position=[4.97080909, 2.91075722, 1.57280005])
C1 = Atom('C', position=[5.12033922, 4.63847287, 4.94075943])
C2 = Atom('C', position=[5.12991894, 3.78609704, 3.79996577])
C3 = Atom('C', position=[4.91462725, 4.08992816, 2.5091264 ])
improper = DihedralAngle(atoms=[N, C1, C2, C3], improper=True)

# The following are some of the equivalent permutations
# The only unique atom location is the first location (central atom)
improper = DihedralAngle(atoms=[N, C2, C1, C3], improper=True)
improper = DihedralAngle(atoms=[N, C3, C1, C2], improper=True)
improper = DihedralAngle(atoms=[N, C1, C3, C2], improper=True)

Currently Periodic is the only InteractionFunction that can be applied to DihedralAngle interactions.

Copying bonded atoms

As mentioned above in the ‘Create an atom’ section, if an Atom with a BondedInteraction is copied, the new atom will also have the same BondedInteraction (and be bonded to the same atom or atoms as the original). For example:

wH1 = Atom('H', position=(5., 5., 5.))
wO = Atom('O', position=(5., 6.63298, 5.))
wBond = Bond((wH1, wO), function=HarmonicPotential(1., 100., interaction_type='bond'))
wH2 = wH1.copy(position=(0., 0.81649, 0.57736))

Both atoms wH1 and wH2 are bonded to wO:

wBond.atoms

[(<Atom
   {name: 'H',
    ID: 16,
    element: 'H',
    position: UnitNDArray([5., 5., 5.]) Ang,
    velocity: UnitNDArray([0., 0., 0.]) Ang / fs}>,
  <Atom
   {name: 'O',
    ID: 17,
    element: 'O',
    position: UnitNDArray([5.     , 6.63298, 5.     ]) Ang,
    velocity: UnitNDArray([0., 0., 0.]) Ang / fs}>),
 (<Atom
   {name: 'H',
    ID: 18,
    element: 'H',
    position: UnitNDArray([0.     , 0.81649, 0.57736]) Ang,
    velocity: UnitNDArray([0., 0., 0.]) Ang / fs}>,
  <Atom
   {name: 'O',
    ID: 17,
    element: 'O',
    position: UnitNDArray([5.     , 6.63298, 5.     ]) Ang,
    velocity: UnitNDArray([0., 0., 0.]) Ang / fs}>)]

Create a molecule

A molecule consists of two or more atoms and at least one bonded interaction:

# Import the Atom, Molecule, Bond, BondAngle and HarmonicPotential classes
from MDMC.MD import Molecule, BondAngle
# Create a H2 molecule
H1.position = [0., 0., 0.]
H2.position = [1., 1., 1.]
H_mol = Molecule(position=[2., 1.5, 1.], atoms=[H1, H2], interactions=[HH_bond])

When a Molecule is created, the position of the atoms relative to one another is fixed. The atoms are then moved so that the position of the molecular centre of mass is what was passed when creating the molecule. In the example above, the atoms were at [0., 0., 0.] and [1., 1., 1.] before H_mol was created; therefore they will always be separated by 1.0 Ang in each dimension, no matter where the molecule is moved to. The molecular centre of mass is set to [2., 1.5, 1.], so the atom positions are changed to [2., 1.5, 1.] and [3., 2.5, 2.] respectively. It is also possible to copy molecules:

H_mol2 = H_mol.copy(position=[5., 5., 5.])

When a Molecule is copied, each Atom is copied, as are all of the bonded interactions between these atoms (and all of the non-bonded interactions).

One method for building molecules is to copy atoms, as the interactions are also copied. For example, to build methane:

m_C = Atom('C')
# Define the H atom positions relative to a C at [0,0,0]
d = 0.629118
H_pos = [[d, d, d], [-d, -d, d], [d, -d, -d], [-d, d, -d]]
m_H = Atom('H', position=H_pos[0])
# Add a bond between
CH_bonds = Bond(m_C, m_H, function=HarmonicPotential(1.09, 100., interaction_type='bond'))
# Make three H atom copies (which are therefore all bonded to m_C)
H_atoms = [m_H]
for pos in H_pos[1:]:
    H_atoms.append(m_H.copy(pos))
# The next two lines simply create a list of unique HCH triplets
# e.g. [(m_H1, m_C, m_H2), (m_H1, m_C, m_H3), ..., (m_H3, m_C, m_H4)].
from itertools import combinations
HCH_triplets = [(i[0], m_C, i[1]) for i in combinations(H_atoms, 2)]
# Unpack the list of triplets with * notation i.e. [(...), (...), (...)] becomes (...), (...), (...)
HCH_angles = BondAngle(*HCH_triplets, function=HarmonicPotential(109.5, 10., interaction_type='angle'))
# Create a methane molecule by adding C atom to list of H atoms
methane = Molecule(atoms=[m_C]+H_atoms)

Creating molecules using configuration files

In addition to creating molecules manually, it is also possible to create them from atomic configuration files (e.g CIF files). Please see the tutorial Reading atoms from configuration files for a detailed description on how to do this.

Add structural units to a universe

There are two methods for adding a structural unit to a universe:

# Add an individual structural unit to the universe
universe.add_structure(H_mol)
# Fill the universe with the structural unit repeated on a cubic lattice with a specific number density
universe.fill(H_mol, num_density = 0.01)

Currently the fill command cannot be used in conjunction with add_structure. When using it with a cubic universe the density will be isotropic, however the exact number of structural units added by fill may be lower than expected as it will always add a cube number. Using fill with a non-cubic universe is not recommended as the density may or may not be isotropic depending on the dimensions of the universe and the number of units.

Create non-bonded interactions

Non-bonded interactions are applied to atoms based on their atom_type, rather than to individual atoms. They must also have a Universe specified, so that they know to which atoms they apply (atoms must have an atom_type once they have been added to a universe). For example, to create a Dispersion interaction with a LennardJones interaction function between two atom types:

# Import the Dispersion interaction and Lennard-Jones function
from MDMC.MD import Dispersion
from MDMC.MD import LennardJones
# Create a dispersion interaction with a Lennard-Jones function between atoms with atom_type 1 (O) and atom_type 2 (H)
# The first LJ parameter is epsilon and the second is sigma
# The cutoff for the dispersion interaction can also be set (in Ang)
LJ_HO = Dispersion(universe, [1, 2], function=LennardJones(0.65, 3.), cutoff=10.)

The exception to this is Coulombic interactions, which can be applied either to a list of atoms or to a list of atom types. If the Coulombic interaction is applied to a list of atoms, the universe does not need to be specified:

# Import the Coulombic interaction
from MDMC.MD import Coulombic
# Create a Coulombic interaction with a Coulomb potential and a charge of 0.42
c_H = Coulombic(atoms=[H1, H2], charge=0.42)

As Coulombic interactions typically have the same interaction function (i.e. a Coulomb function, where the force results in Coulomb’s law), Coulombic interactions do not need to be specified with an interaction function (although other functions can be provided); to set the function as Coulomb, a value for the charge can be passed. The warning highlights the Coulomb interaction function has been automatically set. This is equivalent to:

from MDMC.MD import Coulomb
c_H = Coulombic(atoms=[H1, H2], function=Coulomb(0.42))

As with all non-bonded interactions, a Coulombic interaction can also be created by specifying the atom types:

c_O = Coulombic(universe, atom_types=[1], charge=0.42)

In this case a universe must be provided.

Applying a ForceField

As well as manually setting the Parameter values for each interaction (e.g. Coulombic(atoms=[H1, H2], function=Coulomb(0.42)) to set the charge), it is also possible to apply a ForceField to a Universe, in order to set the Parameter values. Please see the tutorial Applying a ForceField for a detailed description on how to do this.

Adding kspace solvers to a universe

There are several solvers for determining the long range energy contribution for non-bonded interactions, including Ewald, particle-particle particle-mesh (PPPM), and particle-mesh Ewald (PME). If the long range contribution to the non-bonded energy (i.e. > cutoff) is to be calculated during a simulation, a kspace solver has to be added to the universe. A solver can be specified for either electrostatic or dispersive interactions, or both. This can either be during universe initialisation or afterwards:

# Import Ewald and PPPM kspace solvers
from MDMC.MD import Ewald, PPPM
# Create an ewald solver
ewald = Ewald(accuracy=1e-5)
# Initialise a universe with an Ewald solver for both electrostatics and dispersive interactions
uni1 = Universe(10., kspace_solver=ewald)
# Initialise a universe and then add a PPPM solver for electrostatic interactions
uni2 = Universe(10.)
pppm = PPPM(accuracy=1e-4)
uni2.electrostatic_solver = pppm
# Initialise a universe with a PPPM solver for dispersive interactions
uni3 = Universe(10., dispersive_solver=PPPM)

Universe created with:
  Dimensions       [10.0, 10.0, 10.0]
  Force field                    None
  Number of atoms                   0

Universe created with:
  Dimensions       [10.0, 10.0, 10.0]
  Force field                    None
  Number of atoms                   0

Universe created with:
  Dimensions       [10.0, 10.0, 10.0]
  Force field                    None
  Number of atoms                   0

Not all kspace solvers are implemented for all MD engines, and they may also require different parameters to be specified - see the MD engine documentation for more information.

Adding constraint algorithms to a universe

In a similar manner to kspace solvers, constraint algorithms (e.g. SHAKE, RATTLE) can also be passed to a universe. A constraint algorithm is required if any of the bonded interactions are constrained.

# Import Shake and Rattle
from MDMC.MD import Shake, Rattle
# Initialise a universe with the Shake algorithm
# The first Shake parameter is the accuracy and the second is the maximum number of iterations used for any constraint calculation
shake = Shake(1e-4, 100)
uni4 = Universe(10., constraint_algorithm=shake)
# Add Rattle after universe initialisation
rattle = Rattle(1e-5, 1000)
uni5 = Universe(10.)
uni5.constraint_algorithm = rattle

Universe created with:
  Dimensions       [10.0, 10.0, 10.0]
  Force field                    None
  Number of atoms                   0

Universe created with:
  Dimensions       [10.0, 10.0, 10.0]
  Force field                    None
  Number of atoms                   0

Not all constraint algorithms are implemented for all MD engines, and they may also required different parameters to be specified - see the MD engine documentation for more information.

Example Universe filled with water

For water, SPCE, SPC, TIP3P and TIP3PFB forcefields are predefined, so parameters do not need to be set when interactions are defined; instead parameters are set by adding a forcefield to the universe.

Note that all water models have (harmonic) potential strengths defined for their Bonds and BondAngles. In order to create a rigid water molecule in accordance with the models, a constraint_algorithm should be passed to the Universe and constrained=True passed to the relevant interactions, as shown below.

“Constraining” the ability of the bonds to oscillate during MD in this way should not be confused with constraining the value of an MDMC Parameter to a certain numerical range during refinement, as described in Running a Refinement. It is entirely possible to have a rigid bond but allow the length of that bond to change between refinement steps, or conversely have a bond that is free to oscillate during MD but the equilibrium length is not altered as part of the refinement.

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]
  Force field                       None
  Number of atoms                      0

Accessing information

After creating a Universe, the various properties and attributes of it can be accessed as follows:

Geometry

print('The dimension(s) of the Universe are {0}, giving a volume of {1}'
      ''.format(universe.dimensions, universe.volume))

The dimension(s) of the Universe are [21.75 21.75 21.75] Ang, giving a volume of 10289.109375 Ang ^ 3

Forcefield

print('The current force field applied to the Universe is:\n{}'
      ''.format(universe.force_fields))

The current force field applied to the Universe is:
<SPCE
 {interaction_dictionary: {(<class 'MDMC.MD.interactions.Coulombic'>, ('O',)): <Coulomb
 {parameters: {'charge (#54)': <Parameter
 {ID: 54,
  type: 'charge',
  value: -0.8476 e,
  unit: 'e',
  fixed: False,
  constraints: None,
  interactions_name: None,
  functions_name: None,
  tied: False}>}}>, (<class 'MDMC.MD.interactions.Coulombic'>, ('H',)): <Coulomb
 {parameters: {'charge (#55)': <Parameter
 {ID: 55,
  type: 'charge',
  value: 0.4238 e,
  unit: 'e',
  fixed: False,
  constraints: None,
  interactions_name: None,
  functions_name: None,
  tied: False}>}}>, (<class 'MDMC.MD.interactions.Dispersion'>, ('O', 'O')): <LennardJones
 {parameters: {'epsilon (#56)': <Parameter
 {ID: 56,
  type: 'epsilon',
  value: 0.6502 kJ / mol,
  unit: 'kJ / mol',
  fixed: False,
  constraints: None,
  interactions_name: None,
  functions_name: None,
  tied: False}>, 'sigma (#57)': <Parameter
 {ID: 57,
  type: 'sigma',
  value: 3.166 Ang,
  unit: 'Ang',
  fixed: False,
  constraints: None,
  interactions_name: None,
  functions_name: None,
  tied: False}>}}>, (<class 'MDMC.MD.interactions.Bond'>, ('H', 'O')): <HarmonicPotential
 {parameters: {'equilibrium_state (#58)': <Parameter
 {ID: 58,
  type: 'equilibrium_state',
  value: 1.0 Ang,
  unit: 'Ang',
  fixed: False,
  constraints: None,
  interactions_name: None,
  functions_name: None,
  tied: False}>, 'potential_strength (#59)': <Parameter
 {ID: 59,
  type: 'potential_strength',
  value: 4637.0 kJ / mol Ang ^ 2,
  unit: 'kJ / mol Ang ^ 2',
  fixed: False,
  constraints: None,
  interactions_name: None,
  functions_name: None,
  tied: False}>}}>, (<class 'MDMC.MD.interactions.BondAngle'>, ('H', 'O', 'H')): <HarmonicPotential
 {parameters: {'equilibrium_state (#60)': <Parameter
 {ID: 60,
  type: 'equilibrium_state',
  value: 109.47 deg,
  unit: 'deg',
  fixed: False,
  constraints: None,
  interactions_name: None,
  functions_name: None,
  tied: False}>, 'potential_strength (#61)': <Parameter
 {ID: 61,
  type: 'potential_strength',
  value: 383.0 kJ / mol rad ^ 2,
  unit: 'kJ / mol rad ^ 2',
  fixed: False,
  constraints: None,
  interactions_name: None,
  functions_name: None,
  tied: False}>}}>},
  n_body: 3}>

Configuration of atoms and molecules

print('There are {0} atoms in the Universe\n'
      ''.format(universe.n_atoms))

# For a complete list run `print(universe.atoms)`

There are 1029 atoms in the Universe

print('There are {0} molecules in the Universe\n'
      ''.format(universe.n_molecules))

# For a complete list run `print(universe.molecules)`

There are 343 molecules in the Universe

Configuration of interactions

print('There are {0} total interactions in the Universe\n'
      ''.format(universe.n_interactions))

# For a complete list run `print(universe.interactions)`

There are 1032 total interactions in the Universe

print('There are {0} bonded interactions in the Universe\n'
      ''.format(universe.n_bonded))

# For a complete list run `print(universe.bonded_interactions)`

There are 1029 bonded interactions in the Universe

print('There are {0} nonbonded interactions in the Universe\n'
      ''.format(universe.n_nonbonded))

# For a complete list run `print(universe.nonbonded_interactions)`

There are 3 nonbonded interactions in the Universe

Constraints

print('The current constraint algorithm applied to the Universe is {}'
      ''.format(universe.constraint_algorithm.name))

The current constraint algorithm applied to the Universe is Shake

Solvers

Note that a kspace_solver is mutually excusive with the other solver types.

print('The kspace_solver is {}'.format(universe.kspace_solver and universe.kspace_solver.name))
print('The electrostatic_solver is {}'.format(universe.electrostatic_solver and universe.electrostatic_solver.name))
print('The dispersive_solver is {}'.format(universe.dispersive_solver and universe.dispersive_solver.name))

The kspace_solver is None
The electrostatic_solver is PPPM
The dispersive_solver is None

Parameters

print('The MDMC parameters of the Universe are:\n{}'
      ''.format(universe.parameters))

The MDMC parameters of the Universe are:
{'equilibrium_state (#26)': <Parameter
 {ID: 26,
  type: 'equilibrium_state',
  value: 1.0 Ang,
  unit: 'Ang',
  fixed: False,
  constraints: None,
  interactions_name: 'Bond',
  functions_name: 'HarmonicPotential',
  tied: False}>, 'potential_strength (#27)': <Parameter
 {ID: 27,
  type: 'potential_strength',
  value: 4637.0 kJ / mol Ang ^ 2,
  unit: 'kJ / mol Ang ^ 2',
  fixed: False,
  constraints: None,
  interactions_name: 'Bond',
  functions_name: 'HarmonicPotential',
  tied: False}>, 'equilibrium_state (#36)': <Parameter
 {ID: 36,
  type: 'equilibrium_state',
  value: 109.47 deg,
  unit: 'deg',
  fixed: False,
  constraints: None,
  interactions_name: 'BondAngle',
  functions_name: 'HarmonicPotential',
  tied: False}>, 'potential_strength (#37)': <Parameter
 {ID: 37,
  type: 'potential_strength',
  value: 383.0 kJ / mol rad ^ 2,
  unit: 'kJ / mol rad ^ 2',
  fixed: False,
  constraints: None,
  interactions_name: 'BondAngle',
  functions_name: 'HarmonicPotential',
  tied: False}>, 'epsilon (#48)': <Parameter
 {ID: 48,
  type: 'epsilon',
  value: 0.6502 kJ / mol,
  unit: 'kJ / mol',
  fixed: False,
  constraints: None,
  interactions_name: 'Dispersion',
  functions_name: 'LennardJones',
  tied: False}>, 'sigma (#49)': <Parameter
 {ID: 49,
  type: 'sigma',
  value: 3.166 Ang,
  unit: 'Ang',
  fixed: False,
  constraints: None,
  interactions_name: 'Dispersion',
  functions_name: 'LennardJones',
  tied: False}>, 'charge (#15)': <Parameter
 {ID: 15,
  type: 'charge',
  value: 0.4238 e,
  unit: 'e',
  fixed: False,
  constraints: None,
  interactions_name: 'Coulombic',
  functions_name: 'Coulomb',
  tied: False}>, 'charge (#38)': <Parameter
 {ID: 38,
  type: 'charge',
  value: -0.8476 e,
  unit: 'e',
  fixed: False,
  constraints: None,
  interactions_name: 'Coulombic',
  functions_name: 'Coulomb',
  tied: False}>}