Get the parameters of OPLS-AA force field and set up a MD simulation for LAMMPS

Recently some people asked me how to get the parameters of OPLS-AA force field and set up a MD simulation for LAMMPS, here’s a brief tutorial mainly using an example of 10 ns NVT simulation of the system containing 1 solute (R) and 200 methanol solvent molecules. If you have other questions after reading the post, please let me know and I can add more details.

Step 1. Get the pdb files and associated parameters from LigPargen

First, making a pdb file of each molecule you want to add in the MD simulation, here in this tutorial is R.pdb and MeOH.pdb.

Go to LigPargen server, clicking choose file, and read the above pdb file, then click submit molecule, there would be many files for many different MD softwares available, download the GRO, TOP files of GROMACS and the LAMMPS file. If you get errors from LigPargen, you can try using mol files, or you can try SMILE code, for generating SMILE code, you can check this page.

Step 2. Change the formats of the parameters

To use the parameters obtained from LigPargen server with fftool which we will introduce below, we need to modify the formats of the parameters. convertLigParGen.py is a script written by Michael Robinson for this purpose. Following the information written by him, you can easily get the R.ff and MeOH.ff.

If you have ionic liquids in your system, you don’t need to do the first 2 steps, you just copy the coordinates and il.ff file file from here.

Step 3. Produce the input files using fftool

For the details about using fftool to generate LAMMPS input file and data file, you can check this fftool page of Prof. Agilio Padua, which contains very detailed information. Make sure the atomic name in the .xyz or .zmat file used in this step are consistent with that in the .ff files, for example, if H atom is H00 in the .ff file, then you need to make sure it’s also H00 in the .xyz or .zmat files. Also make sure the atomic names are different in each .ff and .xyz fiels, which means if we have a C01 atom in R.xyz and R.ff, there should be no C01 atom in the MeOH.xyz and MeOH.ff, you can change it to C001 for example.

If you want to use the polarizable force field, you could also check the tutorial of the pol_il page. Note that the packmol program will be used to get an initial system. If you want to know more details about using packmol, for example, if you want to place one type of molecules in a certain part of the simulation box, you can check this tutorial written by Dr. Tian Lu. After this step, you should be able to get the data.lmp and in.lmp.

Following information is for polarizable force field, if you are not interetsed in this, please directly move to Step 4.

For simulations with polarizable force field, there are extra files needed, here I just provide basic details, I will write more information if anyone is interested in more. Reading our supporting information might be helpful for you. Generally speaking, if you want to use the Drude polarizable force field, or more specifically, the CL&Pol force field, you need the following parameters: atomic polarizability (in alpha.ff), dipole moment and the centroid distance of dimers (in fragment.ff). Again, for the ionic liquids that are already included in the pol_il, you don’t need to do this, just copying the alpha.ff and fragment.ff into your directory and following the tutorials of the pol_il page is enough. Here we explain how to get the parameters for the species those not included there.

First for the atomic polarizability, there are at least two methods. Method 1 is that used in the 2018 PCCP paper of Schröder and co-workers. Method 2 is to separate the molecular polarizability into atomic polarizability based on the AIM volume, as used in the 2017 J. Condens. Matter Phys.paper of Holm and co-workers. For method 1, you can check the details and scripts provided in the supporting information of 2018 PCCP-1 and 2018 PCCP-2 of Schröder and co-workers. Here we present the details of method 2.

First, we need to get the moelcular polarizability, for this step you can check this c4c1im-polar.gjf input file. Here we use Sadlej’s POL (pVTZ) basis set, I put the POL.gbs here so you can copy the basis set to your local directory and use it. When the calculation is done, you can read the dipole monent $\mu$ (unit in Debye) and the molecular polarizability $\alpha$ (unit in au) from the end of the log file. The dipole moment need to be added into the fragment.ff.The atomic polarizability, $\alpha_{i}$ can be calculated as \begin{equation} \alpha_{i}=\frac{V_{i}}{V_{tot}}\alpha. \end{equation} To calculate AIM volume (unit in Bohr^3) in Multiwfn, you need to open the .wfx file generated in the calculations, then input the following 17-1-1-2(it can be any number from 1 to 4, depending what accuracy/speed you can accept)-7-2(it can be any number from 1 to 3, depending what accuracy/speed you can accept)-1 for each question, for more details, you can read section 4.17.1 of the Multiwfn manual. When this step is done, you can copy the AIM volume results into a spreadsheet like this to obtain the atomic polarizability (I will write a script for this, watch this space), then you copy the atomic polarizability into your alpha.ff file.

The other parameter which is important for the polarizable force field is the centroids distance of dimers, the dimers include all possible cation-anion, cation-neutral, anion-neutral and neutral-neutral combinations of the species in your MD simulation. So the first step is obviously get the most stable structure of each dimer. There are many different ways to do this, for example you can run a non-polarizable MD simulation and extract clusters from the MD trajectory, here I use another method, to produce dimers using the genmer module of the Molclus program and using the GFN2-xTB semi-empirical method in the xtb program to locate the most stable cluster. The most stable cluster is then further optimized by DFT method like M062X/6-31+G(d,p) level of theory, the centroids distance of each dimer can be calculated based on the optimized structure (I will write a script for this, watch this space).

Step 4. Tweak the information of the in.lmp

Techniqually after Step 3, you should be able to submit the simulation. But speficically for my research, I usually need to add electric field information in, and I also need to write the information asscoiated with polarizable force field. I write a simple script to do this automatically, see crin.py. To use it, you just need to python3 crin.py and provide your job name, the strength of the external electric field (along Z direction) and the job type. You need to put the four template files into your local directory like /scratch/x69/lx3539/file/ in my script. If you can help improve the script or/and add more functions into the script, I am very happy to hear your feedbacks.

Some tips to deal with the LAMMPS errors

Two LAMMPS errors you most frequently meet might be Errors: Out of range atoms - cannot compute PPPM and Error: Shake/bond atom missing. Generally whenever you meet whatever errors from LAMMPS, you should ask google first as LAMMPS have an email list (I recommend you to join the email list and read the discussions feed everyday to learn) containing many useful discussions so you can solve most problems youself by using the suggestions people (usually the developers of LAMMPS) provide therein.

For error 1, you can try: (1) checking topology and force field parameters (2) using smaller simulation box (3) using neigh_modify every 1 delay 0 check yes and using neighbor 3.0 bin (4) decreasing the timestep (5) re-building the system using a more reasonable volume, for example, the volume estimated using experimental density.

For error 2, you can try: (1) running an energy minimization before running MD simulation (2) sometimes reducing the number of cores is also useful (3) using a larger value for skin distance.