Creating an atomic configuration

An MDMC simulation requires a configuration and a topology defined within a Universe object. These are orthorhombic (cubic or cuboid) boxes containing a configuration of atoms and molecules. This how-to guide will explain how to create a Universe object and an atomic configuration inside it.

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

from MDMC.MD import Universe

Creating the actual box that contains the atoms is trivial; specify a list or tuple of 3 floats, representing the x, y, and z lengths of the box respectively; these lengths are taken to be in angstroms (Å).

box_universe = Universe((10.0, 15.0, 20.0))

Universe created with:
Dimensions [10. 15. 20.]

To create a cubic Universe, only a single float needs to be specified for the dimensions; MDMC will infer that this is a cube of side length 10Å. This cube will be used as our base universe for the rest of the guide.

universe = Universe(10.0)

Universe created with:
Dimensions [10. 10. 10.]

Create an atomic configuration

Each atom is specified using an Atom object. At a minimum these can be specified just from an element symbol.

from MDMC.MD import Atom

# Create some atoms!
H1 = Atom('H')
O1 = Atom('O', position=(8.0, 8.0, 8.0), atom_type=1)
O2 = Atom(element='O', position=(9.0, 9.0, 9.0), velocity=(0., 0., 0.), charge=None)

Atoms have a few parameters here! Let’s break down some of the particularly important ones: - element: the first parameter gives the element of the atom. If atom mass is not given, it will be determined from an elemental lookup table. - position: the (x,y,z) coordinates of the atom. If not given, defaults to the origin. - velocity: the velocity of the atom along the (x,y,z) axes. Defaults to (0,0,0). - charge: the atomic charge of the atom. Defaults to None (no charge). - atom_type: an ID for the type of atom; if not given, is inferred based on the element and interactions of the atom, such that all atoms with the same element and interactions will have the same atom types. This is used to keep track of atom interactions.

It’s recommended you either define no atom types (and let MDMC do it) or define all atom types. Mixing the two may create unexpected behaviour.

Note that if set, the velocity of atoms 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 providing 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 for 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’.

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

This interaction represents a bond between two atoms. 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 represents bonds at a rotation around a central atom, and 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 interactions can have constraints imposed on them; 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; these are described at the end of this guide.

NB: “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 inRunning 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.

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

Bonded atoms can be combined into a Molecule object, which keeps track of the atoms in the molecule and the bonds between them.

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)

Add structures to a universe

Now we have created our structure objects, we need to add them to our universe box. There are two methods for adding a structure to a universe; either add_structure, which adds an individual structure to the universe, or fill, which fills the universe with copies of the structure. To fill, either specify a desired density of the structure with num_density or a desired number of structures with num_struc_units.

# 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. A list of all structures in the universe can be viewed with universe.structure_list (try it!).

Solving bond constraints

When constraining bonds, constraint algorithms such as SHAKE or RATTLE can also be passed to a universe to solve constraints imposed on bonded interactions. 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. 10. 10.]
Universe created with:
Dimensions [10. 10. 10.]

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.