What is the proper way of doing Free energy calculations in gromacs?

User discussions
1

GROMACS version: 2021.5
GROMACS modification: No
Here post your question

Hi All!

I aim to calculate the change in free energy when charging a molecule. Initially, my system has zero net charge, and I want to assign it a +1 charge in the final state. Based on my review of related literature, I have a few questions:

  1. Some studies report using only a few hundred picoseconds of MD runs with free energy tags. How should I determine the appropriate simulation time for my calculations?
  2. Is it valid to fully equilibrate the system over a few hundred nanoseconds and then apply free energy tags for a shorter production run (e.g., a few hundred picoseconds)?
  3. Should I equilibrate the system with lambda tags before the production run with free energy tags?
  4. If the system ends up with a net charge of +1 in the final state, would that be considered unphysical? Should I include counterions to neutralize the system? If so, how should I handle the charging process of the counterion (e.g., transitioning from a dummy state)?
  5. How can I determine the optimal number of lambda points required for accurate free energy calculations?

Thank you in advance for your guidance and suggestions! I look forward to hearing your insights.

2

This is not a simple question. My reply should be taken as my personal advice, but definitely not as “the truth”.

  1. What I would recommend is to see if there is a trend when you extend the simulations. If you want to check if 500 ps (per lambda point) is enough, try simulating for 1000, or 1500 ps, and see if the results are different and/or if there is a trend in the results when you analyze the first 500 ps, 500-1000 ps and 1000-1500 ps. However, even this way, you might be simulating too short to discover trends etc. In general, I’m very suspicious of running only for a few hundred ps per lambda point, especially if that also includes equilibration time.
  2. I would recommend equilibrating fully in one of the end states. Then you can set-up a sequence of simulations, each starting from the end of the previous one. That means that you can reduce the equilibration time in each lambda state, but you will not avoid it. You will have to experiment with the -b option to gmx bar to see when you discard enough. Using the methods in point 1, above, might help. The problem with this solution is that you can’t run all lambda state simulations in parallel. The other alternative is to run all in parallel, but with significantly longer equilibration times in each lambda state, especially for the states very different from the fully equilibrated state. You can also combine these, i.e., run a sequence of fairly short equilibration simulations, followed by proper full-length production simulations (in parallel). Regarding equilibration, my impression is that the required equilibration time per lambda point is often underestimated. Depending on the transformation you are doing you might need quite a few ns of equilibration at each lambda state, especially if you do not make a sequence of equilibration simulations.
  3. See answer 2.
  4. It is usually a good idea to alchemically turn on a counter ion. I would recommend doing that with a dual-topology approach.
  5. You can use the stdev column in the gmx bar output as a guide. This is an estimate of the quality of the sampling. You want this value to be close to 1.00 for all lambda state transitions. If the value is higher than 1.00 you want a tighter spacing, in general, in that lambda state transition. A value of 2.00 means that a good guess is to put a new lambda point in the middle of that transition. This also means that you can get an estimate of the number of “optimal” (roughly estimated) from the sum of the stdev column. I.e, = sum(stdev) + 1. However, if you add lambda points you will need to rerun the neighboring lambda states as well. So, you don’t want to modify the lambda point distribution unless you really see that it’s needed. I wouldn’t bother unless there are at least a few transitions with stdev higher than 1.3, and perhaps not even then. A rough guide is that you can usually manage with reasonably few (6-11) linearly distributed Coulomb lambda points, but for Lennard-Jones decoupling you might need more, and with a tighter distribution close to, at least one of the two, end states.