Equilibrating our protein-ligand complex will be much like equilibrating any other system containing a protein in water. There are a few special considerations, in this case:

1. Applying restraints to the ligand
2. Treatment of temperature coupling groups

Restraining the Ligand

To restrain the ligand, we will need to generate a position restraint topology for it. First, create an index group for JZ4 that contains only its non-hydrogen atoms:

gmx make_ndx -f jz4.gro -o index_jz4.ndx
...
 > 0 & ! a H*
 > q

Then, execute the genrestr module and select this newly created index group (which will be group 3 in the index_jz4.ndx file):

gmx genrestr -f jz4.gro -n index_jz4.ndx -o posre_jz4.itp -fc 1000 1000 1000

Now, we need to include this information in our topology. We can do this in several ways, depending upon the conditions we wish to use. If we simply want to restrain the ligand whenever the protein is also restrained, add the following lines to your topology in the location indicated:

; Include Position restraint file
#ifdef POSRES
#include "posre.itp"
#endif

; Include ligand topology
#include "jz4.itp"

; Ligand position restraints
#ifdef POSRES
#include "posre_jz4.itp"
#endif

; Include water topology
#include "./charmm36-mar2019.ff/tip3p.itp"

Location matters! You must put the call for posre_jz4.itp in the topology as indicated. The parameters within jz4.itp define a [ moleculetype ] directive for our ligand. The moleculetype ends with the inclusion of the water topology (tip3p.itp). Placing the call to posre_jz4.itp anywhere else will attempt to apply the position restraint parameters to the wrong moleculetype.

If you want a bit more control during equilibration, i.e. restraining the protein and ligand independently, you could instead control the inclusion of the ligand position restraint file in a different #ifdef block, like so:

; Include Position restraint file
#ifdef POSRES
#include "posre.itp"
#endif

; Include ligand topology
#include "jz4.itp"

; Ligand position restraints
#ifdef POSRES_LIG
#include "posre_jz4.itp"
#endif

; Include water topology
#include "./charmm36-mar2019.ff/tip3p.itp"

In the latter case, to restrain both the protein and the ligand, we would need to specify define = -DPOSRES -DPOSRES_LIG in the .mdp file. How you want to treat your system is up to you. These examples are meant only to illustrate the flexibility GROMACS provides. For a standard equilibration procedure, restraining the protein and ligand simultaneously is probably sufficient. Your own needs may vary.

Thermostats

Proper control of temperature coupling is a sensitive issue. Coupling every moleculetype to its own thermostatting group is a bad idea. For instance, if you do the following:

tc-grps = Protein JZ4 SOL CL

Your system will probably blow up, since the temperature coupling algorithms are not stable enough to control the fluctuations in kinetic energy that groups with a few atoms (i.e., JZ4 and CL) will produce. Do not couple every single species in your system separately.

The typical approach is to set tc-grps = Protein Non-Protein and carry on. Unfortunately, the "Non-Protein" group also encompasses JZ4. Since JZ4 and the protein are physically linked very tightly, it is best to consider them as a single entity. That is, JZ4 is grouped with the protein for the purposes of temperature coupling. In the same way, the few Cl^- ions we inserted are considered part of the solvent. To do this, we need a special index group that merges the protein and JZ4. We accomplish this with the make_ndx module, supplying any coordinate file of the complete system. Here, I am using em.gro, the output (minimized) structure of our system:

gmx make_ndx -f em.gro -o index.ndx
...
  0 System              : 33506 atoms
  1 Protein             :  2614 atoms
  2 Protein-H           :  1301 atoms
  3 C-alpha             :   163 atoms
  4 Backbone            :   489 atoms
  5 MainChain           :   651 atoms
  6 MainChain+Cb        :   803 atoms
  7 MainChain+H         :   813 atoms
  8 SideChain           :  1801 atoms
  9 SideChain-H         :   650 atoms
 10 Prot-Masses         :  2614 atoms
 11 non-Protein         : 30892 atoms
 12 Other               :    22 atoms
 13 JZ4                 :    22 atoms
 14 CL                  :     6 atoms
 15 Water               : 30864 atoms
 16 SOL                 : 30864 atoms
 17 non-Water           :  2642 atoms
 18 Ion                 :     6 atoms
 19 JZ4                 :    22 atoms
 20 CL                  :     6 atoms
 21 Water_and_ions      : 30870 atoms

Merge the "Protein" and "JZ4" groups with the following, where ">" indicates the make_ndx prompt:

> 1 | 13
> q

We can now set tc-grps = Protein_JZ4 Water_and_ions to achieve our desired "Protein Non-Protein" effect.

Proceed with NVT equilibration using this .mdp file.

gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -n index.ndx -o nvt.tpr

gmx mdrun -deffnm nvt