Ligand position restraints inclusion statement place

GROMACS version:2020.3
GROMACS modification: Yes/No
I used to run my protein- ligand simulations using the following modifications to the topol.top file:
Before solvation step:

1)Cut [ atomtypes ] section from ligand.itp file, save the file and then paste it under
; Include forcefield parameters
2) Include ligand topology in topol.top file after the position restraint as follows:
; Include ligand topology
#include “ligand.itp”
3) Add in the [ molecules ] section in topol.top :
ligand 1

Then after creating index for ligand and generating position restraint file for ligand after EM:

4)Include posre_ligand.itp in topol.top file after ligand topology and before water topology as follows:
; Ligand position restraints
#ifdef POSRES
#include “posre_ligand.itp”
#endif

I then found others adding an include statement of the ligand topology (ligand.itp) after the forcefield parameters in the system topology file (topol.top), then add the include statement of the ligand position restraints at the end of the ligand.itp file without explicitly writing it in the topol.top file.

I tried that method without changing anything in my mdp files and then I found in the final simulation run for 500 ns that the ligand has escaped from its binding site and is moving freely away from the protein.

So I just wanna know if that has something to do with the place of including the ligand position restraint file, or should I had added some other parameters in my mdp files, if so what are they?

I’m attaching 2 npt.mdp files the one I used in both methods and the other one which was used by others in the second method which I didn’t use, maybe anyone can help me understand what went wrong.
npt_others.mdp (4.1 KB)

npt_mine.mdp (2.6 KB)

Hi,
In your mdp file define = -DPOSRES will trigger the inclusion of posre_ligand.itp into your topology, used for implementing position restraints for the ligand, but not for the protein.
I wonder if it is possible that what you observe is the protein moving out of the ligand (while the ligand is position restraints)?
Best regards
\Alessandra

Hi,
thanks for your reply, but I re-checked and it’s obvious that the ligand suddenly bounce out of its binding site and move freely away from the protein, while the protein remains restrained allover the simulation.

Hi,
I see. Thus in the topology file you have more than one #ifdef POSRE statements.
Probably it helps if you upload a link to your the topology files.
\Alessandra

Hi,
Here is the problematic topology file, where the #ifdef POSRE statement of the ligand is put at the end of the ligand.itp file.
I attached also another topology file, the way I used to do it before.
topol_problematic.top (2.5 MB)
topol_correct.top (700.6 KB)

Please provide ligand.itp as well. Hard to tell anything about the problematic topology without it.

Here it is the ligand.itp file, I just changed the extension now to .top to be able to attach it here.
ligand.top (24.0 KB)

I can’t see any reason why this isn’t working. An #include statement simply means “paste the content of this file here” so the two approaches should be equivalent. Could the displacement be simply a periodicity issue?

I’ve already removed periodicity with gmx trjconv!
The ligand bounces out of the binding site after about 200 ns from 500 ns.
I thought that it maybe due to not adding the #ifdef POSRE statement explicitly at the lower part of the topol.top, not inside the ligand.itp file, but I’m not sure why is that happening, and actually I repeated the simulation 2 more times manually, as I used to do and the same ligand stayed at its binding site over the whole 500 ns!

The #ifdef statement needs to immediately follow the [moleculetype] definition for the ligand, whether that is inside the ligand’s .itp file or immediately following it, makes no difference. But if you did the calculation three times and only had a problem in one, that suggests incorrect use of a define statement in one input file rather than a problem in the topology. Note that #ifdef POSRE would not match define = -DPOSRES but I’m assuming that you’ve just got a typo here in your post.

Sorry, I didn’t get what does that mean, or why could it happen?

Yes, my bad, it’s a typo.

If you have three systems that all use the same topology, you need to check to make sure you are actually invoking the position restraints correctly. If the restraints are active, there will be an energy contribution that will be printed to the .log and .edr file. The other possibility is that this is normal behavior - you mention very long simulations that are performed with the ligand restrained; perhaps the ligand is somewhat weakly bound and the protein to which it is bound simply drifts over time, leaving the ligand behind. But nothing in the files you’ve provided suggests a problem.