LAMMPS

Talk to Rishabh, Jingyang, or Ryan if you have LAMMPS questions!

More Information to come

Creating a LAMMPS data file with relevant force field parameters

1. Create the molecule in Avogadro or wherever you want to

2. Then use antechamber to get the partial charges. Antechamber is part of AmberTools. Installation instructions are available in their manual

3. The relevant command:- antechamber -fi pdb -fo prepi -i <infile.pdb> -o output.prepi -c bcc -j 4 -at gaff -nc -1 The -nc flag is the net charge in the molecule. It will error out if not invoked for ions.

   Sometimes it might be advisable to take partial charges and LJ pair coefficients directly from papers. This is especially important if you are trying to benchmark your results against previous publications. For instance, in case of TFSI-, antechamber partial charges are significantly different than partial charges available in previous publications - https://pubs.acs.org/doi/full/10.1021/jp0686893

   In short, antechamber is a good starting point, but it is always advisable to tread carefully with partial charges in MD. Also remember, that these charges are not scaled. Why is charge scaling important? Go here – ‘Kirby, B. J.; Jungwirth, P. Charge Scaling Manifesto: A Way of Reconciling the Inherently Macroscopic and Microscopic Natures of Molecular Simulations. J. Phys. Chem. Lett. 2019, 10 (23), 7531–7536. https://doi.org/10.1021/acs.jpclett.9b02652.’

4. Open the .pdb file in VMD. Use Extensions -> TK console to open the TCL scripting menu

5. Then use this sequence of commands to get the bonds, angles, dihedrals and most importantly, the box size topo retype bonds topo guessangles topo guessdihedrals set sel [atomselect top all]

   set center [measure center $sel weight none]

   $sel moveby [vecscale -1.0 $center]

   set minmax [measure minmax $sel -withradii]

   set box [vecscale 1.1 [vecsub [lindex $minmax 1] [lindex $minmax 0]]]

   pbc set $box

6. Once all these TCL commands execute without any errors, we can get the lammps data file using: - topo writelammpsdata <lammps.data> full [Note: we can change the atom type to something else like atomic depending on the system, but mostly for molecular systems, full will be the way to go]

7. This is just a skeletal data file with the relevant bonds, angles, and dihedrals. No parameters are currently in this datafile. Sadly, the grunt work starts here. On paper, we can use the moltemplate package to get all the parameters, but it has never worked reliably (at least for me) . The work around is to visit the relevant force field parameter page on the moltemplate GitHub. For instance, for GAFF go to, https://github.com/jewettaij/moltemplate/blob/master/moltemplate/force_fields/gaff.lt and extract all the parameters (LJ, bond, angle, dihedral) from there. Remember to benchmark your LJ parameters.

8. Here is what a finished file will look like

LAMMPS Description 
8 atoms
7 bonds
12 angles
9 dihedrals
0 impropers

4 atom types
3 bond types
4 angle types
1 dihedral types
0 improper types

-3.443375 3.376625  xlo xhi
-3.591275 3.749025  ylo yhi
-2.930725 3.049975  zlo zhi

Masses

1 12.010700 # C
2 18.998404 # F
3 15.999400 # O
4 32.064999 # S

Atoms # full

1 1 3 -0.63000 -1.613375 -0.108625 -1.138875 # O UNL
2 1 4 1.02000 -0.839375 0.519375 -0.059875 # S UNL
3 1 3 -0.63000 -1.476375 0.389375 1.258125 # O UNL
4 1 1 0.35000 0.748625 -0.462625 0.053125 # C UNL
5 1 2 -0.16000 1.596625 -0.000625 1.027125 # F UNL
6 1 2 -0.16000 1.473625 -0.444625 -1.110875 # F UNL
7 1 2 -0.16000 0.539625 -1.787625 0.340125 # F UNL
8 1 3 -0.63000 -0.429375 1.895375 -0.368875 # O UNL

Bonds

1 3 1 2
2 3 2 8
3 2 2 4
4 3 2 3
5 1 4 6
6 1 4 7
7 1 4 5

Angles

1 4 1 2 8
2 1 1 2 4
3 4 1 2 3
4 1 4 2 8
5 4 3 2 8
6 1 3 2 4
7 3 2 4 6
8 3 2 4 7
9 3 2 4 5
10 2 6 4 7
11 2 5 4 6
12 2 5 4 7

Dihedrals

1 1 1 2 4 6
2 1 1 2 4 7
3 1 1 2 4 5
4 1 8 2 4 6
5 1 8 2 4 7
6 1 8 2 4 5
7 1 3 2 4 6
8 1 3 2 4 7
9 1 3 2 4 5

Pair Coeffs #lj/charmm/coul/long

1 0.0660 3.50000000000 #C Kelkar Maginn 
2 0.0530 2.95000000000 #F Kelkar Maginn
3 0.2100 2.95992190114 #O Kelkar Maginn
4 0.2500 3.56359487256 #S Kelkar Maginn

Bond Coeffs #harmonic

1 356.9 1.3497 # SOURCE1_SOURCE5 2188 0.0139, this is c3-f
2 233.5 1.8075 # SOURCE1_SOURCE5 644 0.0143, this is c3-s6
3 512.7 1.4533 # SOURCE1_SOURCE5 2456 0.0105, this is o-s6

Angle Coeffs #harmonic

1 65.4 108.61   # SOURCE3_SOURCE5 1062 1.0758, this is c3-s6-o
2 70.9 107.36   # SOURCE2_SOURCE5 1178 0.5429, this is f-c3-f
3 63.1 109.68   # SOURCE4_SOURCE5 57 0.4273, this is f-c3-s6
4 73.6 120.05   # SOURCE4_SOURCE5 971 1.8153, this is o-s6-o

Dihedral Coeffs #fourier

1 1 0.14445 3 0.0 #F-C-S-O used same parameters as X-c3-s6-X

Lawrencium

module load ms/lammps/3Mar20-mpi