How to use SwissParam files for molecular dynamics simulations with GROMACS

This tutorial shows how to use SwissParam to setup a molecular dynamics (MD) simulation of a protein with a small-molecule ligand in GROMACS, using the CHARMM forcefield. When using the CHARMM forcefield in GROMACS, please cite [1]. When using SwissParam, please cite [2] and [3].

Generally, SwissParam uses the CHARMM36 force field for ligands, so the same should be used for the system in MD simulations. The necessary CHARMM36 force field files for GROMACS are available here and should be copied to $GMXPREFIX/share/gromacs/top/charmm36-*.ff. With the command-line tool of SwissParam, it is also possible to set up ligands with the CHARMM27 force field for use in MD simulations with this older force field, which is included in current GROMACS distributions.

As an example, we look at the protein trypsin with the aminomethylcyclohexane (AMC) ligand, as taken from the PDB entry 1tng.

1. Separate the original pdb file into two pdb files, one for the protein and one for the small molecule.
protein.pdb
ligand_raw.pdb


2. Open the ligand_raw.pdb file in UCSF Chimera or UCSF ChimeraX. Add hydrogens (Tools > Structure Editing > AddHydrogens), check if the protonation state is the desired one, and save in mol2 format.
lig.mol2
For alternative methods to get a Mol2 file for your ligand, see How to obtain a correct Mol2 file.


3. Go to SwissParam and submit the lig.mol2 file. You will get the output in the results page. Uncompress the archive:
> tar -xzf results.tar.gz
For using GROMACS, you will only need the pdb file (with hydrogens) and the itp file with the GROMACS topology of the ligand:
lig.pdb
lig.itp


4. Generate a GROMACS topology for the protein without the ligand. This particular system includes a calcium ion, which has to be renamed to match the CHARMM terminology. In the protein.pdb file, change the calcium atom name and residue name to CAL:
HETATM 2013 CAL CAL A 901 37.738 24.784 36.891 1.00 18.31 CA
Then run the pdb2gmx program from the GROMACS suite:
> gmx pdb2gmx -f protein.pdb -water tip3p -ignh -o conf.pdb


5. Modify the topol.top file generated by pdb2gmx. In the "Include chain topologies" section, add the following statement, before all other lines referring to the protein or ions :
> #include "lig.itp"
This statement has to be at this precise position in the file because it contains [atomtypes] and [pairtypes] sections which have to appear in the topology before any [moleculetype] section. Then, in the [molecules] section, right after lines indicating the number of protein and ion molecules, add the following line :
> LIG 1
If there are crystal water molecules in the structure, there will be a line indicating the number of solvent molecules (SOL) in the topology. In this case, the line above should be inserted before any SOL lines. See below for an example topol.top file.


6. Merge the protein and ligand coordinates. Insert the ATOM lines from the ligand.pdb into the conf.pdb file, right after the line for the calcium ion. Atoms will be automatically renumbered in the next step. Assure that there are no big clashes between the ligand and the protein.
   
   > cat conf.pdb lig.pdb | grep -v END > system.pdb
   
   
7. Now we can use the usual GROMACS setup procedure:
> gmx editconf -f system.pdb -c -d 1.2 -o boxed.pdb
> gmx solvate -cs -cp boxed.pdb -p topol.top -o solvated.pdb


8. You can then run an energy minimization in GROMACS, for example using the em.mdp input file given below:
   > gmx grompp -f em.mdp -c solvated.pdb -p topol.top -o em.tpr
   > gmx mdrun -s em.tpr
   
   
9. It is possible to include position restraints on the ligand by using the POSRES_LIGAND flag. In the case where both protein and ligand atoms should be restrained, the corresponding line in the em.mdp file would read :
> define = -DPOSRES -DPOSRES_LIGAND


10. It is possible to add ions at a physiological concentration and to obtain a net zero charge of the system.
    > gmx genion -s em.tpr -p topol.top -neutral -conc 0.15 -o system-with-ions.pdb

Note

• There are cases where you want to include more that one small molecule in your system. Since all [atomtypes] and [pairtypes] sections have to come before any [moleculetype] section, it is not possible to sequentially include two or more itp files created by SwissParam. Instead, you will have to manually merge the contents of those itp files, such that the right order of sections is maintained.

Example em.mdp file

Example topol.top file

References

1. P. Bjelkmar, P. Larsson, M. A. Cuendet, B. Hess, E. Lindahl, Implementation of the CHARMM Force Field in GROMACS: Analysis of Protein Stability Effects from Correction Maps, Virtual Interaction Sites, and Water Models, J. Chem. Theory Comput., 2010, 6(2), 459-466.
2. Bugnon M, Goullieux M, Röhrig UF, Perez MAS, Daina A, Michielin O, Zoete V. SwissParam 2023: a modern web-based tool for efficient small molecule parameterization. J. Chem. Inf. Model., 2023
3. Zoete V, Cuendet MA, Grosdidier A, Michielin O. SwissParam: a fast force field generation tool for small organic molecules. J. Comput. Chem., 2011, 32(11), 2359-68.