Free energy perturbation for part of a molecule

GROMACS version: 2020.6
GROMACS modification: No

I am trying to perform free energy perturbation on a single functional group of a molecule. I have seen free energy can be calculated with the slow-growth method by a lambda value in the MDP options, or by manually changing the topology parameters for each step between “A” and “B” end-states.

For the former method, this would require creating a separate molecule topology description and creating a covalent bond linking this “molecule” to the remainder of the original molecule. This seems like a herculean task for such a simple requirement. Am I thinking about this the right way? Or is there a way to apply lambda to only a single residue of the molecule?

For the latter, I believe this involves changing the values of the [ atomtypes ] section epsilon and [ atoms ] section charge values for the atoms in question, such that the potential equation would look like:

U(\lambda)=4\sqrt{\epsilon_i\epsilon_j(1-\lambda)^2}[(\frac{\sigma_{ij}}{r_{ij}})^{12}-(\frac{\sigma_{ij}}{r_{ij}})^{6}]+\frac{(1-\lambda)q_iq_j}{r_{ij}}

The squared (1-\lambda) in the square-root is to escape the combination rule of \epsilon_{ij}. In short, the (1-\lambda) would be applied to \epsilon and q for each atom belonging to the functional group in question.

I would appreciate some feedback or guidance on this.

Hi, D_C!

Actually, the perturbation of a single functional group should be straightforward to handle in GROMACS. The topology of your system already contains three columns you’re not using, which define the atom type, mass and charge of this atom in the B state. If you don’t fill in any information there, GROMACS automatically guesses that state A and B are the same. Therefore, to perturb only a functional group you can use your old topology and fill in the columns for state B only for these atoms to encode the perturbation that you want to study. If the perturbation changes the number of atoms, make use of dummy atoms with no charge and vdW interactions. With this setup, you can change lambda via the MDP file and just perform the perturbation you intended to do.

Feel free to post here again if something has remained unclear. Otherwise, good luck with your project!

Best,
Sebastian

1 Like

Oh that’s amazing! I have been using GROMACS for years and there are always new features surprising me. Thanks so much, swingberm!

For anyone who may find this Q&A later, this section of the File formats documentation describes the above comment.