Free energy calculation with fixed relative position

User discussions
,
1

GROMACS version: 2024/1
GROMACS modification: No

Hello everyone!

I need to do an alchemical solvation calculation of a molecule. However, I also need to keep the relative position of the coupled molecule and the center of mass of another group of molecules constant. Could I do that by combining the pull code with restraint_lambdas and the regular vdw lambdas? If so, how can I make sure that the totally decoupled state, with no pull constant k, doesn’t crash due to pull distance higher than 0.49 box size?

I’m doing this so I can evaluate the free energy contribution of the position restraint.

I think this is similar to whats done in this link at the very end, but I couldn’t find a mdp file for this calculation or the analytical expression mentioned to post-process.

https://tutorials.gromacs.org/docs/free-energy-of-solvation.html

Thanks in advance.

2

I think what you describe sounds like Boresch-like restraints (one bond, two angles and three dihedrals). Although in that case, the restraints apply to certain atoms between your coupled molecule and the other group; and you can evaluate the contribution to the free energy analytically. To apply these intermolecular interaction, you can play with the topology. For example, at the end of your topol.top file you can define intermolecular interactions between atoms using absolute numbering (in this case, you must define the type of potential you want to use, as well as the force constants and equilibrium values for when your system is in the coupled or uncoupled state). Naive example:

[ system ]
Receptor and molecule

[ molecules ]
;molecule name nr.
REC 1
MOL 1

[ intermolecular_interactions]
[ bonds ]
; ai aj type bA kA bB kB
1 2 6 0.30 0.0 0.30 4000.0

[ angles ]
; ai aj ak type thA fcA thB fcB
3 1 2 1 120.2 0.0 120.2 41.84
1 2 4 1 101.8 0.0 101.8 41.84

[ dihedrals ]
; ai aj ak al type thA fcA thB fcB
5 3 1 2 2 -110.5 0.0 -110.5 41.84
3 1 2 4 2 50.9 0.0 50.9 41.84
1 2 4 6 2 3.9 0.0 3.9 41.84

3

Thank you!

I’ll look into it.

I’m curious how this could also be done by using the pull feature as mentioned in the above tutorial. Could anyone give me an example of how to do this and then evaluate the free energy of the pull part?