Applying a ForceField

As well as defining Parameter values manually, it is also possible to apply a ForceField (e.g. OPLSAA).

To apply a ForceField, it is necessary to specify the ForceField atom_type. It is important to note that this is different from the ``Atom.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.

from MDMC.MD import Atom, Bond, BondAngle, DihedralAngle, Dispersion, Molecule, Universe

# Alternatively, could simply import all (as shown below).
#from MDMC.MD import *

# Define the unique atoms using the ForceField atom_type
# These can be seen in the oplsaa.dat file (MDMC/MD/force_fields/data/oplsaa.dat)
# The H1 atom will be copied after the bond and bond angles have been defined
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)

Supported DL_POLY version 5.0

All bonded and nonbonded interactions need to be specified as normal, except the function should not be passed as a parameter. This is because the ForceField provides the InteractionFunction as well as the Parameter values.

# 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 bond (CH_bond) and bond angles (HCO_angle and HCH_angle) defined
HC3 = HC1.copy(position=[-0.7076, -1.1754,  0.0000])

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

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 (#5)': <Parameter
 {ID: 5,
  type: 'equilibrium_state',
  value: 109.5 deg,
  unit: 'deg',
  fixed: False,
  constraints: None,
  interactions_name: 'BondAngle',
  functions_name: 'HarmonicPotential',
  tied: False}>, 'potential_strength (#6)': <Parameter
 {ID: 6,
  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.0, 10.0, 10.0]
  Force field                    None
  Number of atoms                   0

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

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