Applying a force field

As well as defining interactions by parameters (e.g. the Lennard-Jones potential), it is also possible to define atomic interaction via a force field. These define the interactions between atoms based on existing experimental or calculated data rather than calculating them during the simulation, and are an attractive choice for accurate simulation.

To apply a ForceField, it is necessary to specify the ForceField’s atom_type. It is important to note that this is different from the Atom’s atom_type which every atom possesses, regardless of whether or not a ForceField is used.

This can be done by either passing the ForceField atom_type or the ForceField atom name as the Atom.name.

It is also important to note that if you require the charge to be set from the ForceField, you should pass this when creating each Atom. Any float can be passed as this will be changed when the ForceField is applied.

Below is an example using methanol and the OPLSAA ForceField. Any other force field has an analogous procedure.

from MDMC.MD import Atom, Bond, BondAngle, DihedralAngle, Molecule, Universe
from MDMC.MD.force_fields.OPLSAA import add_opls_force_field

We first create our methanol molecule.

Note that we are setting the name parameter of the atoms to match the ‘atom type’ IDs used in OPLSAA for the same atoms. This is important, as MDMC will use these names to define the force field interactions.

# Create the atoms
HC1 = Atom("H", position=[-0.7006,  0.3636,  0.8900], name="98")
HC2 = Atom("H", position=[-0.7006,  0.3636, -0.8900], name="98")
HC3 = Atom("H", position=[-0.7076, -1.1754,  0.0000], name="98")
C = Atom("C", position=[-0.3366, -0.1504,  0.0000], name="99")
O = Atom("O", position=[ 1.0849, -0.1713,  0.0000], name="96")
HO = Atom("H", position=[ 1.3606,  0.7699,  0.0000], name="97")

# Create the methanol Molecule
methanol = Molecule(
    atoms=[HC1, HC2, HC3, C, O, HO],
    interactions=[
        # Create the bonds with harmonic potentials
        Bond((C, HC1), (C, HC2), (C, HC3)),
        Bond((O, HO)),
        Bond((C, O), constrained=True),
        # Create the H-C-O bond angles
        BondAngle((HC1, C, O), (HC2, C, O), (HC3, C, O)),
        # Create an HCH bond angle
        BondAngle((HC1, C, HC2), (HC2, C, HC3), (HC3, C, HC1)),
        # Create the H-O-C bond angle
        BondAngle((HO, O, C)),
        # Create the H-C-O-H dihedral
        DihedralAngle((HC1, C, O, HO), (HC2, C, O, HO), (HC3, C, O, HO))
    ]
)

# Create a universe and add the methanol
universe = Universe(dimensions=15.0)
universe.fill(methanol, num_density=0.01)

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

Now that all of the interactions have been defined and the methanol has been added to a Universe, a ForceField can be applied to the Universe:

add_opls_force_field(universe, cutoff=6.0, ewald=1e-4)

This sets all of the InteractionFunction and Parameter values for each Interaction.

Determining the ForceField atom_type or atom name

To determine the correct ForceField atom_type or name for each Atom, there are two methods:

• Search through the .dat file for the ForceField (MDMC/MD/force_fields/data/oplsaa.dat)

• Import the ForceField and use ForceField.filter_element

The latter of these methods is shown below:

from MDMC.MD.force_fields.OPLSAA import OPLSAA
oplsaa = OPLSAA()

# If we wanted to determine the correct type of a Chlorine atom
chlorines = oplsaa.filter_element('Cl')
print(chlorines.to_string())

     atom_type  atom_group                     name  charge element    mass  nbonds
44          45          21  Methylene Chloride (UA)  -0.250      Cl  35.453       1
46          47          21    Chloroform CHCl3 (UA)  -0.140      Cl  35.453       1
48          49          21     Carbon Tetrachloride  -0.062      Cl  35.453       1
167        168          21      Chloroalkene Cl-CH=  -0.120      Cl  35.453       1
205        206          21       Chlorobenzene C-Cl  -0.180      Cl  35.453       1
340        341          21      Chloroalkene Cl2-C=  -0.060      Cl  35.453       1
343        344          21         Chloride Ion Cl-  -1.000      Cl  35.453       0
649        650          21      Cl..CH3..Cl- Sn2 TS  -0.628      Cl  35.453       1
799        800          21      Alkyl Chloride C-Cl  -0.200      Cl  35.453       1
810        811          21     Acyl Chloride Cl-C=O  -0.080      Cl  35.453       1
875        876          21      Chloride Ion (GBSA)  -1.000      Cl  35.453       0

So if the Cl atom is part of a chloroalkene, the correct atom could be created with either of the following:

# Specifying the ForceField atom_type:
Cl = Atom('Cl', name='168', charge=0.)

# Is equivalent to specifying the ForceField atom name
Cl = Atom('Cl', name='Chloroalkene Cl-CH=', charge=0.)