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.visualize import view
from ase.io.opls import OPLSStructure

s = OPLSStructure('172_mod.xyz')  # 172_mod.xyz if the file name for the structure above
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.