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

Matter Modeling Asked by B. Kelly on August 19, 2021

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.

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.

Answered by Godzilla on August 19, 2021

## Related Questions

### Berry’s curvature and magnetic moment in TMDCs

2  Asked on August 19, 2021 by carmen-gonzlez

### Analog computing in matter modeling today: Any applications?

1  Asked on August 19, 2021 by ksousa

### What are the pitfalls for new users of DFT?

4  Asked on August 19, 2021

### Calculating diffusion coefficient from Mean Squared Displacement

2  Asked on August 19, 2021

### What are the types of MCSCF?

3  Asked on August 19, 2021

### What are the types of SCF?

5  Asked on August 19, 2021

### Which schedulers are compatible with a virtual machine?

1  Asked on August 19, 2021

### Is the number of possible Bravais lattices a mathematical fact?

1  Asked on August 19, 2021 by camps

### What are the great unsolved questions in Matter Modeling?

2  Asked on August 19, 2021

### Quantum Dot properties using VASP

2  Asked on August 19, 2021 by suseel-rahul

### Tools for high-throughput DFT studies?

5  Asked on August 19, 2021

### What are the types of Quantum Monte Carlo?

3  Asked on August 19, 2021 by nike-dattani

### Dealing with symmetry of ordered primitive cell during DFT structure relaxation

1  Asked on August 19, 2021 by doublekx

### Deep Neural Networks: Are they able to provide insights for the many-electron problem or DFT?

1  Asked on February 20, 2021

### Temperature effect on elastic constant using VASP

1  Asked on February 15, 2021 by niraja-moharana

### Can we consider SOC in FM materials using QE?

2  Asked on February 12, 2021 by atom

### Why linear response is absent in a non-centrosymmetric system with time reversal symmetry?

0  Asked on February 12, 2021 by walber97

### Formulate the Model Quantum Spin Hamiltonian for low dimensional (1D or 2D ) magnetic materials

0  Asked on January 26, 2021 by chumbak

### What methods are available for excited state calculations in solids?

2  Asked on January 6, 2021 by profm

### How to set vacuum space for slabs?

1  Asked on December 25, 2020 by barix