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
.
[1]:
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)
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.
[2]:
# 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
:
[3]:
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:
[4]:
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
.
[5]:
HCO_angle
[5]:
<BondAngle
{function: <HarmonicPotential
{parameters: {'equilibrium_state (#31)': <Parameter
{ID: 31,
type: 'equilibrium_state',
value: 109.5 deg,
unit: 'deg',
fixed: False,
constraints: None,
interactions_name: 'BondAngle',
functions_name: 'HarmonicPotential',
tied: False}>, 'potential_strength (#32)': <Parameter
{ID: 32,
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:
[6]:
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:
[7]:
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 useForceField.filter_element
The latter of these methods is shown below:
[8]:
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:
[9]:
# 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.)
[ ]: