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.

- a kilojoule can make a large difference to an equilibrium constant
- 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.

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?

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<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.

Answered by Godzilla on August 19, 2021

2 Asked on August 19, 2021 by carmen-gonzlez

2d lattice berry phase dichalcogenides transition metals valleytronics

1 Asked on August 19, 2021 by ksousa

4 Asked on August 19, 2021

2 Asked on August 19, 2021

molecular dynamics numerical convergence statistical thermodynamics

3 Asked on August 19, 2021

casscf multi reference one topic per answer reference request self consistent field

5 Asked on August 19, 2021

ab initio calculations hartree fock one topic per answer reference request self consistent field

1 Asked on August 19, 2021

1 Asked on August 19, 2021 by camps

2 Asked on August 19, 2021

2 Asked on August 19, 2021 by suseel-rahul

5 Asked on August 19, 2021

density functional theory high throughput materials design materials informatics one topic per answer

3 Asked on August 19, 2021 by nike-dattani

diffusion monte carlo monte carlo one topic per answer quantum monte carlo

1 Asked on August 19, 2021 by doublekx

density functional theory molecular dynamics software symmetry vasp

1 Asked on February 20, 2021

density functional theory kohn sham machine learning many body physics quantum monte carlo

1 Asked on February 15, 2021 by niraja-moharana

density functional theory elastic constants phonopy temperature

2 Asked on February 12, 2021 by atom

density functional theory magnetism quantum espresso spin orbit coupling

0 Asked on February 12, 2021 by walber97

0 Asked on January 26, 2021 by chumbak

2 Asked on January 6, 2021 by profm

density functional theory excited states gw bse one topic per answer td dft

1 Asked on December 25, 2020 by barix

Get help from others!

Recent Questions

Recent Answers

- Peter Machado on Why fry rice before boiling?
- Joshua Engel on Why fry rice before boiling?
- haakon.io on Why fry rice before boiling?
- Jon Church on Why fry rice before boiling?
- Lex on Does Google Analytics track 404 page responses as valid page views?

© 2022 AnswerBun.com. All rights reserved. Sites we Love: PCI Database, MenuIva, UKBizDB, Menu Kuliner, Sharing RPP, SolveDir