Applying OPLSAA in lammps

Moltemplate can generate LAMMPS data files with the correct topology
and force-field parameters for OPLSAA. There are several examples you
can download at www.moltemplate.org for building molecules which show
how to use the "oplsaa_moltemplate.py" tool to build molecules using
OPLSAA.

ethylene+benzene_PACKMOL/
waterSPCE+methane/

These are located in the
"all_atom_examples/OPLSAA_force_field_examples/" subdirectory.

There are README files in these directories which explain how to do this.

As for the force-field parameters, they come from TINKER's
implementation of OPLSAA. You must download those parameters from
this web-site, (and edit them, and use them with the
"oplsaa_moltemplate.py" program to generate a moltemplate input file.
See the README file for details)
http://dasher.wustl.edu/tinker/distribution/params/oplsaa.prm

Unfortunately moltemplate currently has two bugs:
(which will be fixed soon)

1) Improper interactions have atoms which may be listed in the wrong
order. (After carefully looking at what went wrong, I believe that
this has the minor effect of effectively changing the strength of the
improper forces by approximately factor of 2.) If you are simulating
simple molecules like alkane chains which do not contain impropers,
then this bug does not effect you.

2) Currently there is currently a bug in moltemplate preventing it
from reading PDB files (or XYZ files). If you don't need to read the
coordinates from PDB or XYZ files, then this bug does not effect you
either.

I have spent a few weeks fixing moltemplate bugs (including these two
bugs) and I will post a major update in the next few days with one or
two new OPLSAA examples. Meanwhile, feel free to post questions to
the LAMMPS mailing list.

Hope this helps

Andrew

Bug fix report for MOLTEMPLATE
  (OPLSAA, AMBER/GAFF, PDB/XYZ files)

  I just wanted to wrap up this bug report.

The two issues I mentioned in my previous email on this thread have
been fixed (regarding OPLSAA impropers and PDB/XYZ file input).
  I also fixed many bugs and added a few new examples, features, and
speed improvements. (I am actually pretty confident now that the
OPLSAA and AMBER force-fields are behaving correctly.)
  Please download the newest version of moltemplate from www.moltemplate.org.

Cheers
Andrew

Hi Andrew,

When following the OPLS-AA procedure of moltemplate to prepare the the data file and input files for using in lammps, what unit should be used? In other words, what is the best unit to be used in the input file, when we have prepared the data file based on OPLS-AA by moltemplate (as the oplsaa.prm file was downloaded from the websites suggested)?
I used real, but I think it is not correct, as the residues don’t have enough mobility.

Thanks
-X

Dear Xiaolin

The oplsaa_moltemplate.py script copies the force field parameters directly out of the “oplsaa_subset.prm” file.
(the “oplsaa.prm” file)
It does not modify them. Incidentally this is the file distributed on the TINKER website that you download at the beginning of the procedure.

What units does “oplsaa.prm” use? I double checked the units tonight by looking at the “vdw” (LJ) parameters for TIP3P water (3.15061 & 0.1521. This is for atom type#63) Those parameters are definitely in Angstroms and kcal/mole. (What LAMMPS calls “real” units.)

I don’t know why your residues are not behaving correctly. Be sure to make sure that the LAMMPS input script you are using contains these lines:

include system.in.init
read_data system.data
include system.in.settings
include system.in.charges

The last file contains atomic charges. Make sure you include it.

units, force-field styles, mixing rules, special_bonds settings, cutuff distances, kspace settings, etc are specified in the system.in.init file. They have all been chosen carefully and (by this point), are hopefully correct. But I would not be shocked if they are not. (It would not be the first time.)

I don’t run all-atom simulations. If you find out what the problem is, please post another message.

Cheers

Andrew

Incidentally, I am toying with the idea of eventually creating an “opls.lt” file that would eliminate the need to use oplsaa_moltemplate.py script. Its low on the priority list right now.

Thanks for your response.
I think there shouldn’t be a problem in my input file, and I guess the difference in the behavior of the residues are related to the nature of the OPLS-AA force field when compared to others like PCFF.

This is my input file:

boundary p p p
units real
atom_style full
bond_style harmonic #I didn’t use in hybrid form
angle_style harmonic #I didn’t use in hybrid form
dihedral_style opls #I didn’t use in hybrid form
improper_style harmonic #I didn’t use in hybrid form
pair_style soft 12
pair_style lj/cut/coul/long 10.0 10.0 #I didn’t use in hybrid form

read_data system.data

group water type 1 2

set type 1 charge -0.847600
set type 2 charge 0.423800
.
.
.

include “system.in.settings”

pair_modify mix geometric tail yes table 12

kspace_style ewald 1e-4 #used ewald instead of pppm
special_bonds lj/coul 0.0 0.0 0.5
neighbor 1 bin
neigh_modify delay 0 every 1 check yes page 1000000
comm_modify mode single cutoff 5 vel no

Thanks, by the way, for your earlier help finding bugs in this script, by the way.

I’m not sure if my reply was very helpful.

Actually, the most likely culprit is probably the dihedral parameters generated by oplsaa_moltemplate.py. Its hard to tell, since Jason and I are getting these numbers by reverse-engineering the format of the oplsaa.prm file which we don’t know. Some force-field file formats assume that the user knows they are supposed to divide the coefficients by 2 (or n). But I did not do that with oplsaa.prm.

Please double check the dihedral parameters and coefficients generated by oplsaa_moltemplat.py are correctly. Let me know if you find anything.

Thanks again

Andrew

Incidentally, I did a quick comparison of one of the dihedral parameters generated by oplsaa_moltemplate.py, and the coreesponding one in “gaff.lt”, for the same 4 atoms. Alas, I could not see the discrepancy. (I’m looking at " @dihedral:hc-c3-c3-hc on line 11621 of “gaff.lt”, and the torsion parameters on line 3669 of “oplsaa.prm”. “gaff.lt” uses 0.15, and “oplsaa.prm” uses 0.3. The factor of two differences is accounted for by the difference between “dihedral_style fourier”, and “dihedral_style opls”. (See the LAMMPS documentation. For the record, “gaff.lt” currently uses “fourier”, and oplsaa_moltemplate.py uses “opls”)

Too lazy to check the other dihedrals.