Defining molecule interactions

In this guide we will define non-bonded interactions between molecules in an MDMC universe, as well as add solvers for long-range interactions. This tutorial focuses on adding interactions defined directly by parameters; for applying force fields, please see this how-to guide.

Firstly, we need to create a universe. For this, we use the universe.fill method as described in creating atomic configurations. This universe is a cube of side length 10Å, filled with Hydrogen atoms to a density of 0.01 atoms per cubic angstrom.

from MDMC.MD import Universe, Atom

universe = Universe(10.0)
universe.fill(Atom('H', atom_type=1), num_density = 0.01)

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

We now apply a non-bonded interaction to the universe. Non-bonded interactions are applied on a per-species (atom_type) basis, rather than to individual atoms. They must also have a Universe specified, so that they know which atoms they apply to. For example, the following code creates a Lennard-Jones dispersive interaction between hydrogen atoms.

from MDMC.MD import Dispersion, LennardJones

interaction = Dispersion(universe=universe,
                        atom_types=(1, 1),
                        function=LennardJones(epsilon=0.65, sigma=3.),
                        cutoff=10.)

The parameters here are: - universe: The universe that the interaction applies to. - atom_types: Which atom types the interaction applies between. If we wanted to have interaction between, say, our hydrogen atoms and some oxygen atoms with atom_type=2, we would use atom_types=[1, 2] here. - function: The function used to calculate the dispersive potential energy between atoms. For Lennard-Jones, they are epsilon and sigma. - cutoff: The cutoff distance beyond which the atoms are assumed to not interact.

Note that the parameters for the interactions (here, epsilon and sigma) are used for the refinement of the simulation.

Most other interactions featured in MDMC (for a full list, see the API reference for interaction functions) work this way, with the exception of the Coulombic (Coulomb’s Law) interaction. This can be defined in two ways; either by atom types, or directly between atoms.

from MDMC.MD import Coulombic

coulombic_universe = Coulombic(universe, atom_types=[1], charge=0.42)  # apply to a universe, on atom type 1

H1 = Atom('H')
H2 = Atom('H', position=[0., 1.63298, 0.])

coulombic_atoms = Coulombic(atoms=[H1, H2], charge=0.42)  # no universe required!

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; charge=0.42 is automatically interpreted as charge=Coulomb(0.42).

Solving long-range interactions

More difficult than determining close-range potential energy between atoms is the contribution of non-bonded interactions at longer distances. Calculating this contribution can be very slow - the time taken to do so quadruples when the number of atoms in the system doubles. To do this more efficiently, MD engines use algorithms known as K-space solvers.

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 you would like to calculate the long range contribution to the non-bonded energy 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 electrostatic 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. 10. 10.]
Universe created with:
Dimensions [10. 10. 10.]
Universe created with:
Dimensions [10. 10. 10.]

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.