SwissParam

How to set up a GROMACS simulation with a molecule parametrized in SwissParam.

Created and maintained by the Molecular modeling group, SIB.

This tutorial shows how to use SwissParam to setup a molecular dynamics 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, pleace cite [2].

As a particular example, we look at the protein trypsin with its aminomethylcyclohexane (AMC) ligand, as taken from the PDB entry 1TNG.

Separate the original pdb file into two pdb files, one for the protein and one for the small molecule.
```
  protein.pdb
  ligand_raw.pdb
```
Open the ligand_raw.pdb file in UCSF Chimera. Choose Tools > Structure Editing > addH. Save as .mol2 file
```
  ligand.mol2
```
For alternative methods to get a .mol2 file for your ligand, see the SwissParam webpage

Go to SwissParam and submit the ligand.mol2 file. You will get the output in your email. Uncompress the archive :
```
  unzip ligand.zip
```
For this tutorial using GROMACS, you will only need the pdb file (with hydrogens) and the .itp file with the GROMACS topology of the ligand :
```
  ligand.pdb
  ligand.itp 
```
Generate a GROMACS topology for the protein without the ligand. This particular system includes a calcium ion. It has to be renamed in the pdb file to match the CHARMM terminology. In the protein.pdb file, change the calcium atom name and the 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:
```
  pdb2gmx -f protein.pdb -ff charmm27 -water tip3p -ignh -o conf.pdb -nochargegrp
```

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 "ligand.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 :
```
  ligand        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.
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 ion CAL. Atoms will be automatically renumbered in the next step.

Now we are back to the usual GROMACS setup procedure.

  editconf -f conf.pdb -o boxed.pdb -c -d 1.2 -bt octahedron

  genbox -cs -cp boxed.pdb -o solvated.pdb -p topol.top

If you have an appropriate input file em.mdp for an energy minimization in GROMACS, you can then run

  grompp -f em.mdp -c solvated.pdb -p topol.top

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
```

Notes

To verify that your solvated molecule looks good in its octahedric box, you have to remember that gromacs always works internally with an equivalent triclinic unit cell. To see the octahedron, you will have to transform your coordinates (after having run grompp), using :
```
  trjconv -s topol.tpr -f solvated.pdb -o compact.pdb -ur compact -pbc mol
```
It can happen that there is more that one small molecule in your system. Due to the fact that 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 topol.top file

  	; Include forcefield parameters
	#include "charmm27.ff/forcefield.itp"

	; Include chain topologies
	#include "ligand.itp"
	#include "topol_Protein_chain_A.itp"
	#include "topol_Ion_chain_A2.itp"

	; Include water topology
	#include "charmm27.ff/tip3p.itp"

	#ifdef POSRES_WATER
	; Position restraint for each water oxygen
	[ position_restraints ]
	;  i funct       fcx        fcy        fcz
	   1    1       1000       1000       1000
	#endif

	; Include topology for ions
	#include "charmm27.ff/ions.itp"

	[ system ]
	; Name
	TRYPSIN in water

	[ molecules ]
	; Compound        #mols
	Protein_chain_A     1
	Ion_chain_A2        1
	ligand		    1
	SOL               163
	SOL             11430

References

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), pp 459–466.
http://old.swissparam.ch / V. Zoete, M. A. Cuendet, A. Grosdidier, O. Michielin, SwissParam, a Fast Force Field Generation Tool For Small Organic Molecules, to be submitted