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

Question

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)

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

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

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

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.

Answer 1

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<e^{-\beta \Delta H}\right>_{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.