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, Dispersion, Molecule, Universe

We first create our methanol molecule. All bonded and nonbonded interactions need to be specified as normal, except the function for the Dispersion should not be passed as a parameter. This is because the ForceField provides the interaction function as well as the parameter values.

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', charge=0., atom_type=1)
C = Atom('C', position=[-0.3366, -0.1504,  0.0000], name='99', charge=0., atom_type=2)
O = Atom('O', position=[ 1.0849, -0.1713,  0.0000], name='96', charge=0., atom_type=3)
HO = Atom('H', position=[ 1.3606,  0.7699,  0.0000], name='97', charge=0., atom_type=4)

# Create the bonds with harmonic potentials
CH_bond = Bond(C, HC1)
CO_bond = Bond(C, O)
OH_bond = Bond(O, HO)

# Create the H-C-O and H-O-C bond angles
HCO_angle = BondAngle((HC1, C, O))
HOC_angle = BondAngle((HO, O, C))

# Create the H-C-O-H dihedral
HCOH_dihedral = DihedralAngle((HC1, C, O, HO))

# Duplicate the HC1 atom
HC2 = HC1.copy(position=[-0.7006,  0.3636, -0.8900])

# Create an HCH bond angle
HCH_angle = BondAngle((HC1, C, HC2))

# Duplicate the HC1 atom again
# This atom will have all existing bond (CH_bond) and bond angles (HCO_angle and HCH_angle) defined
HC3 = HC1.copy(position=[-0.7076, -1.1754,  0.0000])

# remember to add a H1-C-H3 angle! as HC3 is a copy of HC1, it only has HCH_angle
# defined between itself and HC2
H2CH3_angle = BondAngle((HC1, C, HC3))

# Create the methanol Molecule
methanol = Molecule(atoms=[HC1, HC2, HC3, C, O, HO])

# Create a universe and add the methanol
universe = Universe(10.)
universe.add_structure(methanol)

# Add dispersion interactions (no dispersion for HO)
HC_disp = Dispersion(universe, (1, 1))
C_disp = Dispersion(universe, (2, 2))
O_disp = Dispersion(universe, (3, 3))

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

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:

universe.add_force_field('OPLSAA')

This sets all of the InteractionFunction and Parameter values for each Interaction. For example, the charges of all atoms are set:

print('\n'.join(['{0} {1} charge = {2}'.format(atom.element, atom.name, atom.charge) for atom in methanol.atoms]))

H Methanol CH3- charge = 0.04 e
H Methanol CH3- charge = 0.04 e
H Methanol CH3- charge = 0.04 e
C Alcohol CH3OH & RCH2OH charge = 0.145 e
O Alcohol -OH charge = -0.683 e
H Alcohol -OH charge = 0.418 e

Similarly, the InteractionFunction and Parameter values are set for BondedInteraction. The HCO_angle now has a HarmonicPotential function, which has two Parameter values, equilibrium,_state set to 109.5 deg, and potential_strength set to 146.44 kJ / mol deg ^ 2.

HCO_angle

<BondAngle
 {function: <HarmonicPotential
 {parameters: {'equilibrium_state (#16)': <Parameter
 {ID: 16,
  type: 'equilibrium_state',
  value: 109.5 deg,
  unit: 'deg',
  fixed: False,
  constraints: None,
  interactions_name: 'BondAngle',
  functions_name: 'HarmonicPotential',
  tied: False}>, 'potential_strength (#17)': <Parameter
 {ID: 17,
  type: 'potential_strength',
  value: 146.44 kJ / mol rad ^ 2,
  unit: 'kJ / mol rad ^ 2',
  fixed: False,
  constraints: None,
  interactions_name: 'BondAngle',
  functions_name: 'HarmonicPotential',
  tied: False}>}}>,
  constrained: False}>

It is worth noting that applying a ForceField changes each Atom.name from the ForceField atom_type to the ForceField atom name. This can be seen when the charges are printed above, but to emphasize the point:

H_atom = Atom('H', position=[-0.7006,  0.3636,  0.8900], name='98', charge=0., atom_type=1)
print('Before OPLSAA applied, H_atom name is: {}'.format(H_atom.name))

uni = Universe(10.)
uni.add_structure(H_atom)
uni.add_force_field('OPLSAA')

print('After OPLSAA applied, H_atom name is: {}'.format(H_atom.name))

Before OPLSAA applied, H_atom name is: 98
Universe created with:
Dimensions [10. 10. 10.]
After OPLSAA applied, H_atom name is: Methanol CH3-

If these Atom.name were originally set to the Forcefield atom name there is no effect:

H_atom = Atom('H', position=[-0.7006,  0.3636,  0.8900], name='Methanol CH3-', charge=0., atom_type=1)
print('Before OPLSAA applied, H_atom name is: {}'.format(H_atom.name))

uni = Universe(10.)
uni.add_structure(H_atom)
uni.add_force_field('OPLSAA')

print('After OPLSAA applied, H_atom name is: {}'.format(H_atom.name))

Before OPLSAA applied, H_atom name is: Methanol CH3-
Universe created with:
Dimensions [10. 10. 10.]
After OPLSAA applied, H_atom name is: Methanol CH3-

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 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.)