index | modules | gitlab | page source

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