Setting up an OPLS force field calculation#
In order to facilitate the definition of structures for the use of OPLS force fields, there are some helper classes.
Modified xyz#
Suppose, we define the ethanal molecule as extended xyz file
(172_ext.xyz
):
7
Lattice="6.0 0.0 0.0 0.0 6.0 0.0 0.0 0.0 6.0" Properties=species:S:1:pos:R:3:molid:I:1:type:S:1 pbc="T T T"
O 1.613900000000000 -0.762100000000000 -0.000000000000000 1 O
C -0.327900000000000 0.522700000000000 0.000000000000000 1 CT
C 0.392400000000000 -0.722900000000000 -0.000000000000000 1 C
H -0.960000000000000 0.580900000000000 0.887500000000000 1 HC
H -0.960000000000000 0.580900000000000 -0.887500000000000 1 HC
H 0.346400000000000 1.382000000000000 0.000000000000000 1 HC
H -0.104900000000000 -1.581400000000000 -0.000000000000000 1 H1
Then we can read and view the structure using:
from ase.io.opls import OPLSStructure
from ase.visualize import view
# 172_mod.xyz if the file name for the structure above
s = OPLSStructure('172_mod.xyz')
view(s) # view with real elements
elements = {'CT': 'Si', 'HC': 'H', 'H1': 'He'}
view(s.colored(elements)) # view with fake elements
Defining the force field#
The definitions of the force field can be stored in an Amber like style
(172_defs.par
):
# the blocks are separated by empty lines
# comments are allowed
#
# one body - LJ-parameters and charge
CT 0.0028619844 3.50 0.000
C 0.0028619844 3.50 0.000
O 0.0073717780 3.12 -0.683
HC 0.0013009018 2.50 0.000
H1 0 0 0
# bonds
C -CT 13.75 1.522 JCC,7,(1986),230; AA
C -O 24.72 1.229 JCC,7,(1986),230; AA,CYT,GUA,THY,URA
CT-HC 14.74 1.090 changed from 331 bsd on NMA nmodes; AA, SUGARS
C -H1 14.74 1.090
# angles
CT-C -O 3.47 120.40
HC-CT-HC 1.52 109.50
CT-C -H1 3.04 117.00
# dihedrals
O -C -CT-HC 0.0 0.0 0.0 0.0
# cutoffs
C -CT 2.0
C -O 1.8
CT-HC 1.4
C -H1 1.4
C3-O1 1.8 # extra stuff, should not bother
We can write LAMMPS input using the information above:
from ase.io.opls import OPLSff, OPLSStructure
s = OPLSStructure('172_ext.xyz')
with open('172_defs.par') as fd:
opls = OPLSff(fd)
opls.write_lammps(s, prefix='lmp')
which writes the LAMMPS input files lmp_atoms
defining atoms, bonds,
etc., and lmp_opls
defining the corresponding OPLS force field. A
rudimentary lmp_in
is also written.