Are there any MD packages that do proper free energy sampling in the NPT ensemble?

Statistical Mechanics is the basis of molecular level calculations of properties and averages. Nowadays free energy calculations are fairly "turn-the-crank", which is not necessarily a good thing. It is easy to perform a free energy calculation and get a number.

In the molecular dynamics package GROMACS, for instance, alchemical free energy calculations are only valid if done in the NVT ensemble. This is because GROMACS samples potential energies, thus one can calculate (or press a button on pymbar or another piece of software)

$$$$\Delta A = -kT \ln \langle \exp(-\beta \Delta \mathcal{U}) \rangle$$$$

Where $$\Delta A$$ is the difference in the Helmholtz free energy between two states, $$\beta = \frac{1}{kT}$$, $$\mathcal{U}$$ is the difference in sampled potential energy between two states, and $$\langle \rangle$$ is an ensemble average.

If you run an NPT simulation you must instead sample

$$$$\Delta G = -kT \ln \frac{\langle V \exp(-\beta \Delta \mathcal{U}) \rangle}{\langle V \rangle}$$$$

While it is tempting to say that the volumes cancel so you only need the ensemble average of the NVT Boltzmann factor, this is not at all correct. GROMACS is aware that it is not correctly sampling the NPT and has tucked away a note in the user manual to the effect that for water at 298.15 K the error is small i.e., less than a kilojoule.

1. a kilojoule can make a large difference to an equilibrium constant
2. Chemistry does not happen at only 298.15! if you increase the temperature or density of the solvent (not all of us do water) the error quickly gets out of hand - this is particularly true for ionic liquids.

Question:

Is this assumption (possibly error in other codes) present in other codes? Are there any MD codes that properly sample the NPT ensemble for free energy?

Note:

While Godzilla has shown entirely correctly how one should sample the free energy in the NPT, neither his equation, or mine, which is equivalent for $$\Delta G$$, are implemented in GROMACS, and possibly others. I am looking for if any MD programs do implement the proper free energy sampling. They all properly do NPT sampling (I assume), but sampling the free energy is a step above and beyond.

• Luckily I don't yet have the critic badge, so I can't ever be blamed for it! It wasn't Godzilla either, because they got the badge in May, but still have only downvoted once. Unless maybe they un-downvoted and re-downvoted. It's one of 57 people that hold the critic badge: Not me, not you, and not Gozilla. A true mystery indeed! Jul 18 '20 at 19:59

It is straightforward to show that in a typical $$NPT$$ setting the Zwanzig equation still only depends on the energy difference and not on the volume (here I define $$H$$ to be the Hamiltonian of each system, respectively and $$x$$ to represent all variables over phase space):

$$\frac{Z_{B,NPT}}{Z_{A,NPT}} = \int \int e^{-\beta \Delta H(x)} \frac{e^{-\beta H_A(x)}}{Z_{A,NPT}} e^{-\beta PV} dxdV = \left_{A,NPT}$$

so I am not sure where you get this equation from $$-$$ it seems wrong to me. Of course, you still get $$NPT$$ contribution in this way, since the exponential averages have to be performed in the isothermal-isobaric ensemble, and these will generally be different to averages obtained in $$NVT$$ (although in many practical purposes, virtually the same, especially in dense systems that are not near any critical points).

I am not sure which part of the GROMACS manual you are referring to, but I suspect you are talking about the Berendsen barostat. The Berendsen barostat is known to be not rigorous but very good for equilibration (again, for many systems it is practically the same as a rigorous barostat and usually not the limiting error). For a production run you should use the Parrinello-Rahman barostat, which is considered to be rigorous. However, it is deterministic, which can lead to some practical problems regarding non-physical oscillations. This is the "rigorous" barostat implemented in GROMACS and to my knowledge this way of calculating free energies in $$NPT$$ is completely valid from a theoretical standpoint.

Another rigorous barostat which is also stochastic is the Monte Carlo barostat. Unfortunately, it is only implemented in AMBER and OpenMM and not in GROMACS. As far as I know, NAMD comes with a Nosé-Hoover Langevin piston barostat which is also rigorous. You can see that pretty much every MD engine comes with a Berendsen barostat and a single rigorous barostat and they all use different rigorous barostats.

As a final remark, if you are planning to run simulations at high temperatures, I'd be much more worried about the validity of your model, presumably a force field, at these extreme conditions. Also, while it is true that 1 kJ/mol can be considered significant in some cases, there are so many errors when calculating free energies, including force field accuracy, insufficient sampling, estimator convergence (e.g. Zwanzig), or even the fact that you restrict your system to be a tiny box, that I would be much less worried about the impact of the barostat, or indeed the choice of thermodynamic ensemble on my results in most conceivable applications.

• +10. Godzilla is back! Jul 17 '20 at 0:51
• I've always been here, I am just waiting to prey on molecular dynamics questions :) Jul 17 '20 at 0:51
• I saw that one, but I am afraid quantum methods are not my forte - I've never done ab initio MD in my life and I gravitate towards classical MD and free energies more than anything else Jul 17 '20 at 0:55
• Your equation is correct @Godzilla123, my question is if it is implemented in software. It is not implemented in GROMACS. Gromacs does not calculate the difference in enthalpy and put that in the boltzmann factor. If you take your equation, you get mine with manipulation. Neither is wrong. My question is if there is software doing this proper sampling - they leave out the V terms in my equation, and the PV terms in yours. The issue is not the Barostat :) Jul 18 '20 at 7:44
• Ah, I think I am getting your problem now. What I would do is that after I perform the free energy perturbation, I would go to the energy.xvg file, which contains all the modified $\Delta U$s and also the $PV$ term for each snapshot. Then I would read the file using a parser (for example alchemlyb) and modify each (reduced) $\Delta u$ to $\Delta u - \ln V + \ln\left<V\right>$ (you just need to make sure you are using the correct volumes for each $\Delta u$). And then I would pass these modified $\Delta u$s to pymbar and let it do its magic. Would that not work? Jul 18 '20 at 8:38