Restraints free energy

User discussions
,
1

GROMACS version: Any version
GROMACS modification: No

Hi everyone,

thank you for your interest. I am applying very simple [ distance_restraints ] consisting in 2 distances between 2 atoms of a ligand and 2 alpha carbons of the hosting pocket. In doing so, I am annihilating the ligand, the objective being a binding free energy calculation. I do see the restraints are correctly present and act accordingly thanks to very informative content you provided in this forum: that guided me from scratch to run and successfully see the log file reporting the corresponding contribute to the energy. What I don’t see is the derivative of the potential with respect to these restraints in the dhdl.xvg output file (I am not committing the naif error of putting a series of nulls in the restraint-lambdas vector), which I need to recover the proper binding free energy. I know the manual is clear about this: that derivative will not be reported in my very case (looking at restraint-lambdas description).

The question(s):
I wonder how to deal with the above-mentioned “recovery” in this very case, or if the only alternative is the pull code. In both cases, I’d like to kindly ask for a reference/guide. In the former case: any mdp file keyword stressing that the intermolecular interactions put in the top file need to be fainted/increased with lambda? In the latter : as I’m new to the pull code, I will provide a line at the bottom exemplifying what I’d like to reproduce via the pull code itself so to ask what to do correspondingly.

““Supporting Info””:
For clarity sake, I am reporting (supposedly) useful information: disre = simple and refcoord-scaling = com is all I’ve got with respect to this restraints in my mdp file for npt conditions and, at the same time, the top file is provided with distance restraints following the scheme of this_table, which is also exemplified right below.

;  ai      aj       funct   index   funct'       low               up1              up2          factor
 1324  2136     1          0          2    0.3045544   0.4045544   0.4545544        1

Thank you in advance for your time.

Sincerely,
Jacopo

2

Got it: restrain with 1 bond, 2 angles and 3 dihedrals roots_to_Boresch, analytical formula for recovery; a very valid option

I have a technical problem, though. Putting [ angle_restraints ] below [ intermolecular_interactions ] at the end of the top file, I get:

WARNING 1 [file topol.top, line 62]:
  Some parameters for bonded interaction involving perturbed atoms are
  specified explicitly in state A, but not B - copying A to B

The line of the top file is exactly the one referring to the very first [ angle_restraints ] occurrence. If I comment that, the warning points to the second, if I comment that, too, it points to the [ dihedral_restraints ] first occurrence. The same warning does not appear if adopting solely the [ distance_restraints ]. Should I simply don’t mind, as I am annihilating and not mutating? (couple-lambda0 = vdw, couple-lambda1 = none, no hybrid or dual topology needed, and moreover all the same is working without warnings if I do not put the [ angle/dihedral_interactions ] headers at the bottom of the top file