Restraints with FEP

User discussions
1

GROMACS version: All
GROMACS modification: No

Hi,

I am starting FEP calculations and I am a little bit confused with the restraints. I want to compute the absolute binding free energy of a very small ligand to a target protein, and I need to add restraints otherwise when both coulombic and vdw interactions are set to zero the ligand wanders. I will probably use 2 or 3 harmonic potentials. What is the recommanded strategy with Gromacs: a section “[ intermolecular_interactions ]” in the topology file, or using the pull code in the .mdp file?

Then, I was thinking using restraint-lambdas in the .mdp file, but I see in this tutorial Absolute Binding Free Energy - Gromacs 2016 - AlchemistryWiki that they are using bonded-lambdas which confuses me. Maybe they are using bonded-lambdas because they are using intermolecular_interactions?

Finally, I don’t understand how Gromacs can be aware of which restraints to add/remove: if I use the pull code to keep my ligand in place, but also to add another restraint far from the binding site, when the restraints are added with restraint-lambdas how Gromacs knows which one to add? Is it because in
he .mdp I say the I want do apply FEP to the molecule LIG and then Gromacs checks which coordinates from the pull code involve atoms from LIG? (and the same applies when intermolecular_interactions are used).

Thank you for your help

Nicolas

2

Hi,

A little follow-up following some tests.

I first tried to use “restraint_lambdas” in the mdp file, and nothing has worked. I tried:

  • to use the pull code (with pull-coord-type being either umbrella or constraint), since the user guide says that for restraint_lambdas “Only the restraint interactions: dihedral restraints, and the pull code restraints are controlled with this component of the lambda vector.”
  • to use the intermolecular_interactions section in the topology with a section distance_restraints where I gave r0/r1/r2/fac for the state A and the state B. Gromacs said “Incorrect number of parameters - found 10, expected 6 or 6 for Dis. Rest. (after the function type)”
  • to use the intermolecular_interactions section in the topology with the bonds section (giving r0 and k for state A and state B)

What worked was using the intermolecular_interactions section in the topology with the bonds section with bonded_lambdas in the mdp (and not with restraint_lambdas), as done in the tutorial linked above.

Thus, I really don’t understand how one can use restraint_lambdas in mdp.

Nicolas

3

What I did not understand is the values r0=0.654 , thA= 88.8, thB= 32.9 they have used. And how to get these values and from where did they arrive on the r0 value of a bond which does not exist.