An Assessment of Approximate Methods for Anharmonic Free Energies

TL;DR

This paper evaluates the accuracy and efficiency of approximate methods for calculating anharmonic free energies in materials, comparing them to reference path integral calculations across various systems.

Contribution

It provides a comprehensive benchmark of vibrational self-consistent field and self-consistent phonon methods against path integral results for diverse materials.

Findings

Methods work well for simple solids like carbon allotropes.

Methods are unreliable for complex molecular systems like ice and paracetamol.

No systematic improvement over harmonic approximation for relative free energies.

Abstract

Quantitative evaluations of the free energy of materials must take into account thermal and zero-point energy fluctuations. While these effects can easily be estimated within a harmonic approximation, corrections arising from the anharmonic nature of the interatomic potential are often crucial and require computationally costly path integral simulations. Consequently, different approximate frameworks for computing affordable estimates of the anharmonic free energies have been developed over the years. Understanding which of the approximations involved are justified for a given system is complicated by the lack of comparative benchmarks. To facilitate this choice we assess the accuracy and efficiency of some of the commonly used approximate methods -- the vibrational self-consistent field and self-consistent phonons -- by comparing the anharmonic correction to the Helmholtz free energy against reference path integral calculations. These benchmarks are performed for a diverse set of systems, ranging from simple quasi-harmonic solids to flexible molecular crystals with freely-rotating units. Our results suggest that for simple solids such as allotropes of carbon these methods yield results that are in excellent agreement with the reference calculations, at a considerably lower computational cost. For more complex molecular systems such as polymorphs of ice and paracetamol the methods do not consistently provide a reliable approximation of the anharmonic correction. Despite substantial cancellation of errors when comparing the stability of different phases, we do not observe a systematic improvement over the harmonic approximation even for relative free-energies. Our results suggest that efforts towards obtaining computationally-feasible anharmonic free-energies for flexible molecular solids should therefore be directed towards reducing the expense of path integral methods

Figures

Equations

$$\hat{H}=-\sum _{p,i}\frac{{\hslash }^2}{2m_i}{\nabla }_{{\mathbf{r}}_{pi}}^2+V\left(\left\{{\mathbf{r}}_{pi}\right\}\right)$$

$$\hat{H}=-\frac{{\hslash }^2}{2}{\nabla }_{\mathbf{r}}{\mathsf{M}}^{-1}{\nabla }_{\mathbf{r}}^T+V\left(\mathbf{r},\mathsf{h}\right)$$

$$Z(N,\mathcal{V},\beta )=\text{Tr}\left[\text{exp}\left(-\beta \hat{H}\right)\right],$$

$$A(N,\mathcal{V},\beta )=-{\beta }^{-1}\ln Z(N,\mathcal{V},\beta ).$$

$$V^{\text{har}}\left(\mathbf{r},\mathsf{h}\right)=V^{(0)}+\frac{1}{2}(\mathbf{r}-{\mathbf{r}}^0)\mathsf{K}(\mathbf{r}-{\mathbf{r}}^0{)}^T$$

$$\mathsf{K}={\mathsf{M}}^{\frac{1}{2}}\mathsf{U}{\mathbf{\Omega }}^2{\mathsf{U}}^T{\mathsf{M}}^{\frac{1}{2}}=\tilde{\mathsf{U}}{\mathbf{\Omega }}^2{\tilde{\mathsf{U}}}^T$$

$$\begin{array}{rl}{\hat{H}}^{\text{har}}& =-\frac{{\hslash }^2}{2}{\nabla }_{\mathbf{q}}^2+\frac{1}{2}\mathbf{q}{\mathbf{\Omega }}^2{\mathbf{q}}^T+V({\mathbf{r}}^0)\\ & =V^{(0)}+{\hat{H}}^{\text{cm}}+\sum _{i=1}^d\left[-\frac{{\hslash }^2}{2}{\nabla }_{q_i}^2+\frac{1}{2}{\omega }_i^2{q}_i^2\right]\end{array}$$

$$|{\Psi }_{\left(\mathbf{s}\right)}^{\text{har}}⟩=\prod _{i=1}^d\otimes |{\varphi }_{i,(s_i)}⟩$$

$$A^{\text{har}}(N,\mathcal{V},\beta )=V^{(0)}+A^{\text{cm}}+\sum _{i=1}^d\left[\frac{\hslash {\omega }_i}{2}+{\beta }^{-1}\ln \left(1-e^{-\beta \hslash {\omega }_i}\right)\right].$$

$$V(\mathbf{q})=V^0+\sum _i^d{V}^{(1)}(q_i)+\frac{1}{2}\sum _i^d\sum _{j\ne i}^d{V}^{(2)}(q_i,q_j)+\cdots ,$$

$$V^{(1)}(q_i)=V(0,\dots ,q_i,\dots ,0)-V^{(0)},$$

$$\begin{array}{rl}{V}^{(2)}(q_i,q_j)=& V(0,\dots ,q_i,\dots ,q_j,\dots ,0)\\ & -V^{(1)}(q_i)-V^{(1)}(q_j)-V^{(0)}.\end{array}$$

$${\hat{H}}^{\text{imf}}=V^{(0)}+\sum _i^d\left[-\frac{{\hslash }^2}{2}{\nabla }_{q_i}^2+V^{(1)}(q_i)\right].$$

$$A^{\text{imf}}(N,\mathcal{V},\beta )=V^{(0)}+A^{\text{cm}}-\sum _{i=1}^d\left[{\beta }^{-1}\ln \sum _{s_i}\text{exp}\left(-\beta E_{i,(s_i)}^{\text{imf}}\right)\right].$$

$$\left[-\frac{{\hslash }^2}{2}{\nabla }_{\mathbf{q}}^2+V(\mathbf{q})\right]|{\Psi }_{\left(\mathbf{s}\right)}^{\text{vscf}}⟩=E_{\left(\mathbf{s}\right)}^{\text{vscf}}|{\Psi }_{\left(\mathbf{s}\right)}^{\text{vscf}}⟩$$

$$\left[-\frac{{\hslash }^2}{2}{\nabla }_{q_i}^2+{\bar{V}}_i(q_i)\right]|{\psi }_{i,(s_i)}^{\text{vscf}}⟩=E_i^{\text{vscf}}|{\psi }_{i,(s_i)}^{\text{vscf}}⟩$$

$${\bar{V}}_i(q_i)=\sum _{j\ne i}\rho (q_j)V(\mathbf{q})$$

$$\rho (q_j)=\frac{{\sum }_{s_j}\text{exp}\left(-\beta E_{j,(s_j)}^{\text{vscf}}\right)|{\psi }_{j,(s_j)}^{\text{vscf}}⟩⟨{\psi }_{j,(s_j)}^{\text{vscf}}|}{{\sum }_{s_j}\text{exp}\left(-\beta E_{j,(s_j)}^{\text{vscf}}\right)}.$$

$$A^{\text{vscf}}(N,\mathcal{V},\beta )=V^{(0)}+A^{\text{cm}}-{\beta }^{-1}\ln \sum _{\mathbf{s}}\text{exp}\left(-\beta \sum _i{E}_{i,(s_i)}^{\text{vscf}}\right).$$

$$E_{\left(\mathbf{s}\right)}^{\text{vscf},(2)}=\sum _{{\mathbf{s}}^{\prime }\ne \mathbf{s}}\frac{\text{matrixelement}{\Psi }_{\left({\mathbf{s}}^{\prime }\right)}^{\text{vscf}}V(\mathbf{q})-{\sum }_i{\bar{V}}_i(q_i){{\Psi }_{\left(\mathbf{s}\right)}^{\text{vscf}}}^2}{E_{\mathbf{s}}^{\text{vscf},(1)}-E_{{\mathbf{s}}^{\prime }}^{\text{vscf},(1)}}$$

$$\begin{array}{r}{A}^{\text{mp2}}(N,\mathcal{V},\beta )=A^{\text{vscf}}(N,\mathcal{V},\beta )-{\beta }^{-1}\ln \sum _{\mathbf{s}}\text{exp}\left(-\beta E_{\mathbf{s}}^{\text{vscf},(2)}\right)\end{array}$$

$$A<A^{\text{scp}}=\text{expectationvalue}{\hat{H}+{\beta }^{-1}\ln {\hat{\rho }}^{\text{scp}}}_{{\hat{H}}^{\text{scp}}};{\hat{\rho }}^{\text{scp}}=\frac{\text{exp}\left(-\beta {\hat{H}}^{\text{scp}}\right)}{\text{Tr}\left[\text{exp}\left(-\beta {\hat{H}}^{\text{scp}}\right)\right]}$$

$${\hat{\rho }}^{\text{scp}}(\mathbf{r})=(2\pi \mathsf{D}{)}^{-\frac{1}{2}}\text{exp}\left(-\frac{1}{2}(\mathbf{r}-{\mathbf{r}}^{\text{scp}}){\mathsf{D}}^{-1}(\mathbf{r}-{\mathbf{r}}^{\text{scp}}{)}^T\right)$$

$${\widetilde{q_i}}^{\text{scp}}(T)=\sqrt{\frac{\hslash }{2{\omega }_i^{\text{scp}}}\coth \frac{\hslash {\omega }_i^{\text{scp}}}{2k_{B}T}}$$

$$\begin{array}{rl}\text{expectationvalue}{{\mathbf{f}}_{\mathbf{r}}(\mathbf{r})}_{H^{\text{scp}}}& =0\\ \text{expectationvalue}{\mathsf{K}(\mathbf{r})}_{H^{\text{scp}}}& ={\mathsf{K}}^{\text{scp}}\end{array}$$

$$\begin{array}{rl}{\mathsf{K}}_{{\hat{H}}_{l+1}^{\text{scp}}}^{\text{scp}}& =\text{expectationvalue}{\mathsf{K}(\mathbf{r})}_{{\hat{H}}_l^{\text{scp}}}\\ {\mathbf{r}}_{{\hat{H}}_{l+1}^{\text{scp}}}^{\text{scp}}& ={\mathbf{r}}_{{\hat{H}}_l^{\text{scp}}}^{\text{scp}}+{{\mathsf{K}}_{{\hat{H}}_{l+1}^{\text{scp}}}^{\text{scp}}}^{-1}\text{expectationvalue}{{\mathbf{f}}_{\mathbf{r}}}_{{\hat{H}}_l^{\text{scp}}}\end{array}$$

$$\begin{array}{rl}{A}^{\text{scp}}& =A^{\text{cm}}+\left[\frac{\hslash {\omega }_i^{\text{scp},l}}{2}+{\beta }^{-1}\ln \left(1-e^{-\beta \hslash {\omega }_i^{\text{scp},l}}\right)\right]\\ & +\text{expectationvalue}{V(\mathbf{r})-\frac{1}{2}(\mathbf{r}-{\mathbf{r}}_l^{\text{scp}}){\mathsf{K}}_l^{\text{scp}}(\mathbf{r}-{\mathbf{r}}_l^{\text{scp}}{)}^T}_{{\hat{H}}_l^{\text{scp}}}.\end{array}$$

$$\Delta A_{\text{cl}}=A_{\text{cl}}-A_{\text{cl}}^{\text{har}}={\int }_0^1\text{d}\lambda \left(\frac{\partial A}{\partial \lambda }\right)={\int }_0^1\text{d}\lambda \text{expectationvalue}{V-V^{\text{har}}}_{H^{\lambda }},$$

$$\Delta A_{\text{cl}}=-T{\int }_0^T\text{d}\tilde{T}\frac{\text{expectationvalue}{V-V^{(0)}-\frac{3N}{2}{k}_B\tilde{T}}_{\tilde{T}}}{{\tilde{T}}^2},$$

$$\Delta A_{\text{qn}}=A_{\text{qn}}-A_{\text{cl}}={\int }_0^1\text{d}g{g}^{-1}\text{expectationvalue}{\hat{T}-T_{\text{cl}}}_{{\hat{H}}^g}.$$

Full text

1 Introduction

The free energy is a key thermodynamic quantity which provides a measure of phase stability. Knowledge of the free energy and its derivatives with respect to temperature and applied fields can, in principle, be used to calculate every other thermodynamic observable. Reliable predictions of free energies from atomistic simulations remain a challenge because they require an accurate description of inter-atomic interactions, as well as proper treatment of the statistical mechanics of the nuclear degrees of freedom. The availability of computational resources combined with developments in electronic structure theory 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 have made it possible to calculate the Born-Oppenheimer (BO) surfaces that govern nuclear motion routinely and accurately. Consequently, the accuracy of free energy calculations is often limited by the statistical sampling of the nuclear degrees of freedom 11. This is most commonly performed within a harmonic approximation, which is reasonable for quasi-harmonic systems, such as metals at low temperatures 12, but fails close to the melting temperature 13, 12 and in the presence of defects 14. The problem is exacerbated for the case of organic solids, which require a proper description of anharmonicity arising from quantum nuclear motion even at room temperature 15, 16, 17, 18. This is similarly true for (i) systems containing light elements, such as hydrogen 19, 20, 21, helium 22 and water 23, ice 24, 25, 26, metal organic frameworks  27, 28, the record high-$T_{c}$ conventional superconductor SH3 29, (ii) systems of reduced effective dimensionality such as graphene 30, 31 and (iii) molecular systems such as paracetamol 18.

Within the BO approximation32 and given the BO potential, exact anharmonic free energies can be calculated using approaches based on imaginary time path integral (PI) simulations 33, 34. However, the required number of force evaluations, when combined with accurate electronic structure calculations, render these approaches impractical for any but the smallest systems. Consequently, a small zoo of frameworks has been developed which approximately account for quantum anharmonic motion at a much lower computational cost. These invoke different approximations and exhibit different scaling behaviour with system size.

In this work, we present an extensive benchmark of the accuracy of some of the most common approximate techniques — such as the harmonic approximation (HAR) 35, self-consistent phonons (SCP) 36, 37, the independent mode framework (IMF) 38, and the vibrational self-consistent field (VSCF) 38 — against reference results obtained using PI thermodynamic integration (QTI) 39, 18. Computationally efficient algorithms for these methods have been developed and implemented in the universal force engine i-PI 40. The accuracy and the computational efficiency of the methods is tested on a set of solids ranging from from simple allotropes of carbon, anharmonic but relatively rigid polymorphs of ice to polymorphs of paracetamol that contain (nearly) freely rotating internal degrees of freedom.

2 Theory

To briefly outline the different free energy methods we discuss, we consider a three-dimensional periodic system, whose minimum potential energy, “equilibrium” atomic positions form a Bravais lattice, noting that finite and aperiodic systems simply represent the limit of infinite period. The full ionic Hamiltonian of such a system is

[TABLE]

where $m_{i}$ is the mass of nucleus $i$, $V$ is the BO potential governing nuclear motion, and $p$ and $i$ run over the Bravais points and the nuclei within a unit cell, respectively. In practice we perform supercell simulations using periodic Born-von Karman simulation cells $\boldsymbol{\mathsf{h}}$ consisting of $N_{a}\times N_{b}\times N_{c}$ replicas of the unit cell and with cell vectors $N_{a}{\bf R}_{a}$, $N_{b}{\bf R}_{b}$, and $N_{c}{\bf R}_{c}$.

In the following we only consider $\Gamma$-point vibrational motion within the simulation cell. We thereby sample those ${\bf K}$-points within the first vibrational Brillouin Zone (BZ) of the underlying unit cell, for which $\textrm{exp}\left(-i{\bf K}\cdot{\bf R}\right)=1\quad\forall\quad{\bf R}=N_{a}{\bf R}_{a}+N_{b}{\bf R}_{b}+N_{c}{\bf R}_{c}$. The Hamiltonian of the system is then uniquely defined given the positions of the $N$ particles within the simulation cell $\boldsymbol{\mathsf{h}}$:

[TABLE]

where $(\nabla_{{\bf r}},{\bf r})\equiv(\{\nabla_{r_{1}},\dots,\nabla_{r_{3N}}\},\{r_{1},\dots,r_{3N}\})$ denotes the momenta and positions of the $3N$ degrees of freedom associated with the $N$ particles and $\boldsymbol{\mathsf{M}}=\text{Diag}\left[m_{1},\dots,m_{3N}\right]$. The canonical partition function of the system at inverse temperature $\beta=(k_{B}T)^{-1}$ and volume $\mathcal{V}=\textrm{Det}\left[\boldsymbol{\mathsf{h}}\right]$ is defined as

[TABLE]

where the trace can be performed over any complete basis set. In the thermodynamic limit, the Helmholtz free energy of the system is

[TABLE]

Direct computation of $A$ is hindered by the computational complexity of solving the Schrödinger equation associated with the Hamiltonian $\hat{H}$, motivating approximate but computationally more affordable approaches.

2.1 Harmonic Approximation

For small displacements, ${\bf r}-{\bf r}^{0}$, of the particles from their equilibrium positions, ${\bf r}^{0}\equiv\textrm{argmin}_{{\bf r}}~{}V({\bf r},\boldsymbol{\mathsf{h}})$, the potential can be Taylor expanded. Truncation after the quadratic term amounts to the harmonic approximation

[TABLE]

with $V^{(0)}\equiv V({\bf r}^{0})$ and $\boldsymbol{\mathsf{K}}=\left.\nabla^{2}V({\bf r},\boldsymbol{\mathsf{h}})\right|_{{\bf r}={\bf r}^{0}}$. The spectral decomposition of the Hessian,

[TABLE]

provides the unitary matrix $\boldsymbol{\mathsf{U}}$, the mass-scaled transformation matrix $\tilde{\boldsymbol{\mathsf{U}}}$, and the diagonal matrix containing the normal mode frequencies $\boldsymbol{\Omega}=\text{Diag}\left[\omega_{1},\dots,\omega_{3N}\right]$. After transformation to the normal mode coordinates $\nabla_{{\bf q}}\equiv\tilde{\boldsymbol{\mathsf{U}}}^{T}\nabla_{{\bf r}}$ and ${\bf q}\equiv\tilde{\boldsymbol{\mathsf{U}}}({\bf r}-{\bf r}^{0})$, the Hamiltonian

[TABLE]

separates into $V^{(0)}$, the centre of mass term $\hat{H}^{\textrm{cm}}$, and a term describing a system of $d=3N-3$ independent simple harmonic oscillators (SHO) whose energies and wave functions for a given excitation state $s_{i}$, $E_{i,(s_{i})}^{\textrm{har}}=(s_{i}+1/2)\hbar\omega_{i}$ and $|\phi^{(s_{i})}_{i}\rangle$, are known analytically. In finite systems, global rotations decouple analogously. The centre of mass contribution to the free energy $A^{\textrm{cm}}$ is that of a free particle in a three dimensional box with a volume and shape equal to that of the Wigner-Seitz cell of the system, while the contribution from the free rotations of finite systems can be computed within the rigid rotor approximation41.

For the remaining system of harmonic oscillators, the $d$-body wave function of the global state described by the $d$-tuple ${\bf s}=\left(s_{1},\dots,s_{d}\right)$ is a Hartree product of the independent normal mode wave functions:

[TABLE]

and the free energy is:

[TABLE]

2.2 Independent Mode Framework

A first approximation to anharmonic quantum nuclear motion is detailed in the work of Monserrat et al. 38. The potential is expanded in terms of the normal mode coordinates

[TABLE]

where

[TABLE]

is the (anharmonic) independent mode term and

[TABLE]

describes pairwise coupling between normal modes. This expansion can be continued for more general $n$-body terms $V^{(n)}$. Since one starts with the harmonic approximation, in which the normal modes are non-interacting, the hope is that higher-order terms decrease in size with increasing $n$. The validity of this assumption is discussed in section 4. Truncation after $V^{(1)}$ amounts to the independent mode approximation with the Hamiltonian

[TABLE]

Despite the presence of anharmonicity the normal modes remain independent and a Hartree product analogous to Eq. (8) of anharmonic normal mode wave functions solves the Schrödinger equation yielding the eigenvalues $E^{\textrm{imf}}_{i,(s_{i})}$. The Helmholtz free energy is

[TABLE]

2.3 Vibrational Self-consistent Field

Retaining terms involving $V^{(2)}$ (and / or higher order terms) leads to coupling of the previously independent normal modes and complicates the solution of the Schrödinger equation. Monserrat et al. 38 solve the equation

[TABLE]

within the iterative Vibrational Self-Consistent Field (VSCF) approach, where $V({\bf q})$ represents the truncated form of Eq. (10). Using a Hartree product trial wavefunction amounts to a mean-field (MF) treatment and leads to the VSCF equations

[TABLE]

where ${\bar{V}}_{i}(q_{i})$ is the mean-field potential experienced by normal mode $i$,

[TABLE]

with

[TABLE]

To lowest order the VSCF free energy becomes

[TABLE]

A perturbation theory can be constructed in terms of the (assumed to be small) difference between the mapped out and the MF potential, $V({\bf q})-\sum_{i}{\bar{V}}_{i}(q_{i})$, leading to a second-order MP2 correction to the energy of state ${\bf s}$ given by

[TABLE]

and the approximate free energy

[TABLE]

2.4 Self-consistent Phonons

Another way of calculating an anharmonic correction on top of a harmonic approximation exploits the Gibbs-Bogoliubov inequality42, which states that the true free energy of a system is always bounded from above by the free energy $A^{\textrm{scp}}$ computed using a trial density matrix, $\hat{\rho}^{\textrm{scp}}$:

[TABLE]

where $\expectationvalue{\square}_{\hat{H}^{\textrm{scp}}}=\textrm{Tr}\left[\rho^{\textrm{scp}}~{}\square\right]$ is an ensemble average defined by the the trial density matrix $\hat{\rho}^{\textrm{scp}}$. Within the self-consistent phonons method43, 44, $\hat{\rho}^{\textrm{scp}}$ is the density matrix of a harmonic Hamiltonian with Hessian $\boldsymbol{\mathsf{K}}^{\textrm{scp}}$ and equilibrium positions ${\bf r}^{\textrm{scp}}$:

[TABLE]

where, $\boldsymbol{\mathsf{D}}=\boldsymbol{\mathsf{M}}^{-\frac{1}{2}}\boldsymbol{\mathsf{U}}\boldsymbol{\mathsf{Q}}^{2}\boldsymbol{\mathsf{U}}^{T}\boldsymbol{\mathsf{M}}^{-\frac{1}{2}}$ and $\boldsymbol{\mathsf{Q}}$ is a diagonal matrix containing the root-mean-square (RMS) displacements 36:

[TABLE]

of the normal modes. The lowest upper bound to the true free energy is obtained by minimizing the free energy with respect to ${\bf r}^{\textrm{scp}}$ and $\boldsymbol{\mathsf{K}}^{\textrm{scp}}$. This leads to the steady state conditions36

[TABLE]

where ${\bf f}_{{\bf r}}({\bf r})$ correspond to the forces of the potential $V({{\bf r}})$. The solution is obtained in a self consistent manner by starting with educated guesses of $\left({\bf r}^{\textrm{scp}},\boldsymbol{\mathsf{K}}^{\textrm{scp}}\right)=\left({\bf r}^{\textrm{scp}}_{0},\boldsymbol{\mathsf{K}}^{\textrm{scp}}_{0}\right)$ – which are in practice chosen to be those obtained within the harmonic approximation – and updating

[TABLE]

until convergence is achieved. Here $\hat{H}^{\textrm{scp}}_{l}$ denotes the the trial Hamiltonian of the $l$-th SCP iteration. The resultant free energy at the $l$-th iteration is calculated as:

[TABLE]

2.5 Thermodynamic Integration

Within the thermodynamic integration scheme the free energy differences between two states is calculated as the work to reversibly transform one state into the other 45, 46, 47. For solids this method can be used to calculate the classical anharmonic correction to the harmonic Helmholtz free energy as the reversible work done while “switching on" the anharmonic part of the potential 48, 49. In the Hamiltonian $H^{\lambda}=(1-\lambda)H^{\textrm{har}}+\lambda~{}H$ the Kirkwoord coupling parameter $\lambda$ smoothly switches the potential from harmonic $(\lambda=0)$ to fully anharmonic $(\lambda=1)$. The free energy difference is obtained by computing the integral of the thermodynamic force along the switching path:

[TABLE]

where, $\expectationvalue{\square}_{H^{\lambda}}$ represents an average over the classical canonical ensemble sampled by the intermediate Hamiltonian and $A_{\textrm{cl}}^{\textrm{har}}$ is the classical harmonic free energy. Setting aside statistical errors, $\Delta A_{\textrm{cl}}$ can be computed exactly by sampling the thermodynamic forces at multiple values of $\lambda\in[0,1]$ using molecular dynamics simulations. Alternatively, $\Delta A_{\textrm{cl}}$ can also be calculated by performing a thermodynamic integration from a low temperature harmonic solid 50, 51 as

[TABLE]

where $\expectationvalue{\square}_{\tilde{T}}$ is an average over the $N\mathcal{V}\tilde{T}$ ensemble.

To include quantum anharmonic corrections due to zero- point energy, tunnelling, etc., a second thermodynamic integration must be performed to calculate the work required to reversibly transform the particles from classical to quantum 52, 53. This can be achieved by defining the Hamiltonian $\hat{H}^{g}=-g\frac{\hbar^{2}}{2}\nabla_{{\bf r}}\boldsymbol{\mathsf{M}}^{-1}\nabla_{{\bf r}}^{T}+V\left({\bf r},\boldsymbol{\mathsf{h}}\right)$, where $g$ scales the mass of the particles39, 16, 54. As $g$ is varied from 1 to 0, the de Broglie wavelength of the particles smoothly drops from its physical value to zero, yielding the desired transformation from quantum to classical particles. The corresponding free energy difference is:

[TABLE]

where $\expectationvalue{\hat{T}}_{\hat{H}^{g}}$ represents the average kinetic energy for the intermediate Hamiltonian and $T_{\textrm{cl}}$ is the classical kinetic energy which is independent of the mass of the system. Eq. (30) can be computed exactly (modulo statistical error) by sampling the quantum canonical ensembles for $g\in[0,1]$ using PI molecular dynamics (PIMD). The difference between the classical and quantum kinetic energy can be computed directly using a centroid-virial kinetic energy estimator 55, 22. The total anharmonic free energy is computed as:

[TABLE]

3 Implementation

In order to perform a direct comparison of the different approximate methods, IMF, VSCF, and SCP were implemented within i-PI 40, an open-source python package for atomistic simulations, which couples to a variety of density-functional-theory (DFT) and empirical and machine-learning potential codes. The IMF, VSCF, and SCP implementations are schematically shown in Figs. 1 to 3. The reference free energies can simply be evaluated on the basis of (PI)MD simulations and require no dedicated implementation.

3.1 Independent Mode Framework

Given the transformation matrix, ${\tilde{\bf U}}$, and normal-mode frequencies from the harmonic approximation, we perform single-point energy and force evaluations for equally spaced configurations $q_{i}^{j}=jf\tilde{q_{i}}(T)$ along each normal mode $i$ where $\tilde{q_{i}}$ is the RMS displacement of the normal mode at a target temperature $T$. We increase $j$ by one at a time until the sampled energy $V^{(1)}(q_{i}^{j})$ exceeds a user-defined multiple $n_{E}$ of the thermal harmonic energy $V^{(1)}(q_{i}^{j})>n_{E}E^{\textrm{har}}_{i}(T)$. This ensures that the potential is always mapped out far enough into the classically forbidden region (but only as far as necessary) to localize the nuclear density, at temperatures lower or equal than the chosen target. The independent mode potential $\sum_{i}V^{(1)}(q_{i})$ is then reconstructed by fitting cubic splines to $\{(q_{i}^{j},V^{(1)}(q_{i}^{j}))\}$. The corresponding independent mode Hamiltonian is expanded in a basis of SHO eigenstates and diagonalized to evaluate the independent mode anharmonic Helmholtz free energy. The Helmholtz free energy is converged with respect to the density of the frozen-phonon samples $q_{i}^{j}$ by repeatedly halving $f$ and supplementing the already collected $\{(q_{i}^{j},V^{(1)}(q_{i}^{j}))\}$ with corresponding samples, until the required convergence threshold is met. For each $f$ the Helmholtz free energy is converged with respect to the size of the SHO basis.

3.2 Vibrational Self-Consistent Field

The implementation of the VSCF framework is split into two modules: one for mapping the potential energy surface (PES), and one for solving the VSCF problem. The mapping strategy mirrors that employed in the IMF module. In a first loop over normal modes, we collect $\{(jf\tilde{q_{i}}(T),V^{(1)}(jf\tilde{q_{i}}(T))\}$ until the sampled potential exceeds a user-defined multiple of the harmonic energy $V^{(1)}(jf\tilde{q_{i}}(T))>n_{E}E^{\textrm{har}}_{i}(T)$, thereby also determining the sampling range for the coupling corrections. In a second loop over $n$-tuples of normal modes, we then sample $\{((q_{i_{1}}^{j_{1}},\ldots,q_{i_{n}}^{j_{n}}),V^{(n)}(q_{i_{1}}^{j_{1}},\ldots,q_{i_{n}}^{j_{n}}))\}$ in a similar fashion and extract the coupling corrections $V^{(n)}(q_{i_{1}},\ldots,q_{i_{n}})$ using cubic spline fits. Currently, sampling and fitting of $n=2,3$ are implemented. The extracted coupling corrections are stored for use within the VSCF solver module.

The module for solving the VSCF problem consists of two submodules, the first of which performs the VSCF calculation itself. The thermal density determining the mean-field potentials $\{{\bar{V}}_{i}(q_{i})\}$ is initialised as the IMF thermal density. Within a VSCF step the MF independent mode Hamiltonians for the given MF potentials are constructed, expanded in an SHO basis, and diagonalized to determine the updated MF thermal vibrational density and the free energy $A^{\textrm{vscf}}$. To stabilize the VSCF convergence, 50% of the thermal density resulting from the previous VSCF iteration are mixed in before the mean-field potentials $\{{\bar{V}}_{i}(q_{i})\}$ are updated and the next VSCF step is initiated. This is repeated until self-consistency has been reached as indicated by convergence of the associated free energy $A^{\textrm{vscf}}$ to within the required threshold.

The second submodule allows the calculation of an MP2 correction on top of the MF eigenstates and -energies by looping over pairs of eigenstates $({\bf s},{\bf s}^{\prime})$ to evaluate the MP2 corrections in Eq. (20) on a real-space grid of predefined density. Care is taken to only consider eigenstates of the self-consistent MF description with eigenenergies $E^{\textrm{vscf}}_{{\bf s}}$ within a set multiple of $k_{B}T$.

3.3 Self-Consistent Phonons

Our implementation of the SCP method is schematically shown in Fig. 3. In the first loop over SCP steps, we construct the trial density matrix $\hat{\rho}^{\textrm{scp}}({\bf r})$ using the mean position ${\bf r}^{\textrm{scp}}$ and Hessian $\boldsymbol{\mathsf{K}}^{\textrm{scp}}$ obtained from the previous step. In the first step $\hat{\rho}^{\textrm{scp}}({\bf r})$ are those obtained within the harmonic approximation. In the second loop we calculate the ensemble averages of the forces and the Hessian, necessary to perform the optimization steps described in equation 26. These are realized as Gaussian integrals and computed using MC importance sampling as

[TABLE]

where $N_{s}$ is the number of samples. Samples are generated by translating $3(N-1)$-tuples of (quasi-)random numbers on the interval $[0,1]$ into atomic displacements from the mean position ${\bf r}^{\textrm{scp}}$ using the inverse cumulative distribution function of $\rho^{\textrm{scp}}({\bf r})$ with a Beasley-Springer-Moro algorithm 56, 57. To speed up the convergence of the averages with respect to the number of samples, we employ the following tricks:

For small system sizes, instead of drawing pseudo random numbers, we use low-discrepancy quasi-random numbers – specifically Sobol sequences 58 – as was done in the implementation of Brown and coworkers 36. This leads to a more uniform sampling, so that error in equation 32 decays as $\mathcal{O}\left(\ln(N_{s})^{d}/N_{s}\right)$ which becomes $\sim\mathcal{O}\left(1/N\right)$ for low dimensional integrals59. For large system sizes, we resort to pseudo random numbers as the performance of Sobol sequences degrades60. We use the FORTRAN implementation of Burkardt 61 to generate Sobol sequences. 2. 2.

As was done in the implementation of Errea and co-workers62, we re-use samples from previous SCP iterations via a reweighting scheme. Given the updated trial density $\hat{\rho}^{\textrm{scp}}_{l}$ at the $l$-th SCF iteration, the reweighted average using the $N_{s}$ samples $\{{\bf r}\}^{k}$ drawn from the trial density $\hat{\rho}^{\textrm{scp}}_{k}$ at the $k$-th SCP iteration is

[TABLE]

We minimize the error in the global estimates $\expectationvalue{\square}_{l}$ at the $l$-th SCP iteration by weighting samples drawn in the $k$-th SCP iteration according to a “batch weight”

[TABLE]

where $\square^{k}\equiv\square(\{{\bf r}\}^{k})$ and the variance of a generic observable over samples from $k$-th SCP iteration is 63

[TABLE]

where $w_{l}^{k}\equiv w_{l}^{k}\left(\{{\bf r}\}^{k}\right)$, provided both $\square$ and $-\ln{w_{l}^{k}}$ are normally distributed. Neglecting $\expectationvalue{-\square\ln{w_{l}^{k}}}$ renders the batch weights independent of the observable being considered

[TABLE]

and thereby suitable for both Hessians and forces. 3. 3.

Taking inspiration from stochastic over-relaxation algorithms64, 65, we always draw pairs of configurations $({\bf q}^{i},{\bf q}^{i+1})$, where ${\bf q}^{i+1}=-{\bf q}^{i}$, ensuring that forces from the symmetric part of $V$ cancel out exactly. 4. 4.

To compute the average Hessian we use an integral by parts, as suggested in reference 36, and to further reduce the variance, we express it in terms of the difference between the harmonic and the anharmonic forces:

[TABLE]

Samples are drawn in sets of $N_{s}$ until at least one component of the average forces (in terms of normal mode coordinates) is statistically significant, as assessed by whether the average over samples is larger than the standard deviation.

We only update ${\bf r}^{\textrm{scp}}$ along those normal modes which exhibit a statistically significant net force, $\expectationvalue{{f_{{\bf q}}}_{i}}$, where $\expectationvalue{missing}={\tilde{\bf U}}^{T}\expectationvalue{{\bf f}_{{\bf r}}({\bf r})}$. Direct application of Eq. 26 (b) in Cartesian space may lead to instability due to the residual statistical errors. The optimization continues until no statistically significant force component remain or the “batch weights” become smaller than a preset threshold, at which point a new SCP iteration begins.

While all modes must be real upon convergence, insufficient statistics may lead to spurious imaginary modes with $\omega_{i}^{2}<0$ before convergence is achieved. Such imaginary modes are treated by setting $\omega_{i}^{2}=-\omega_{i}^{2}$ in the effective harmonic description.

4 Results and discussion

We first describe the systems that have been studied and the potentials that have been used to describe inter-particle interactions. We then investigate the scaling of the computational cost of the methods with respect to the system size, before assessing their accuracy by systematically comparing the approximate free energies to reference quantum thermodynamic integrations. We neglect the centre of mass contribution to the free energy throughout as it cancels out when comparing methods. We converge all results with respect to the sampling of the vibrational BZ by increasing the simulation cell size, allowing us to compare the different methods in equivalent and physically meaningful conditions.

4.1 Systems and Computational Details

Three sets of materials are studied in this work in order of increasing complexity. As a first example we consider the diamond 66 and lonsdaleite 67 allotropes of carbon. These differ only in the stacking of hexagonal bilayers of tetrahedrally-coordinated carbon atoms. Their room temperature densities are identical to within experimental error 66, 67, and equal to $3.51\textrm{g}/\textrm{cm}^{3}$. We consider simulation cells containing up to 64 atoms, starting from the two- and four-atom primitive cells for diamond and lonsdaleite respectively. All the cells were designed to be as close to cubic as possible to render the effective sampling of the vibrational BZ as uniform as possible. Inter-atomic interactions are modelled using the Gaussian Approximation Potential (GAP) of Deringer and Csányi 68, which is based on LDA DFT calculations for configurations from MD simulations of liquid and amorphous carbon. For crystalline carbon (including diamond and graphite) it has been shown to reproduce DFT energies and forces to within RMS errors of 2 meV/atom and 0.1 eV/Angstrom, respectively.

As a second benchmark, we discuss two proton-ordered polymorphs of ice, hexagonal (XIh) 69, 70 and cubic (XIc) 71 ice. These become thermodynamically (meta-)stable below the experimental transition temperature for proton-disordering of 72 K 69, 70. For XIc we assume $I4_{1}/amd$ symmetry, noting that the true experimental structure of XIc is still under debate 71. In direct analogy to the above carbon allotropes the oxygen sublattices of XIh and XIc only differ in the stacking of bilayers of tetrahedrally-coordinated oxygen atoms. In view of the absence of experimental data for sufficiently pure XIc we take its density to be identical to that of XIh. We use the experimental density at ambient pressure and 10 K of $0.93\textrm{g}/\textrm{cm}^{3}$, noting that the thermal expansion of ice XIh between 10 K and 70 K is less than 0.5 %. We use simulation cells containing up to 16 molecules to allow for the possibility of coupling between pseudo-translations, which are not present at the $\Gamma$-point of the unit cell, and librational, bending and O–H bond stretching modes. The interatomic interactions are described using a Behler-Parinello type neural network (NN) 72, based on B3LYP+D3 DFT reference calculations for around 20,000 liquid water and hexagonal ice configurations from MD and PIMD trajectories. This potential successfully reproduces the density of states, pair correlation functions and energy fluctuations of B3LYP+D3 liquid water 73 and has been used to study the quantum kinetic energy, proton momentum distribution, and vibrational density of states of solid and liquid water 74, 75.

Finally we analyze two polymorphs of paracetamol (N-acetyl-p-aminophenol), the monoclinic form I 76 and the orthorhombic form II 77. The two forms differ in the packing of hydrogen bonded sheets of molecules – zig-zag for form I and almost planar for form II. We consider the conventional unit cells containing four and eight formula units for forms I and II respectively, at room temperature experimental densities. Inter-atomic interactions are described on the basis of the Merk Molecular Force field also used in Ref. 18. The accuracy parameter of the PPPM method78 used for calculating electrostatic interactions was set to $10^{-6}$ fractional error in the individual force components, which is smaller than the value usually required, to ensure a smooth PES. While this simple potential contains harmonic terms for bonds and angles, it remains highly anharmonic as the dihedral interaction term describes a (almost) free rotation of the methyl groups at room temperature.

4.2 Computational cost

We define the computational cost of an approximate method ($\square$) as the minimum number of energy or force evaluations required to converge the per atom anharmomic free energy ($A^{\square}-A^{\textrm{har}}$) to within 10% of the reference value for the largest system size. For the case of diamond, this tolerance is equal to a stringent 0.2 meV/atom. Fig. 4 shows how the cost of these methods and the reference QTI scales with the number of atoms in the simulation cell.

The cost of IMF calculations depends linearly on the number of normal modes – which scales linearly with system size – and the number of points sampled along each mode. In our implementation, the later remains weakly dependent on the potential due to the variable, dynamically optimized sampling point density.

Analogously the VSCF approach exhibits a rough scaling behaviour of $N^{m}$ where $m$ is the dimensionality of the potential surfaces that are being sampled. For the case of diamond we use $m=2$ and therefore observe a $N^{2}$ dependence for large $N$. As the anharmonicity of the potential for diamond is very much dependent on BZ sampling (see SI S.1.1), the cost for the primitive cell is an outlier.

The cost for the SCP scheme, using pseudo random numbers, arises from the use of Monte Carlo importance sampling of the optimal effective harmonic description, which scales independently of system size. The statistical reweighting scheme reduces the cost for small sizes but becomes increasingly less beneficial as $N$ is increased63. The net result is a near linear scaling behaviour for the system sizes that we have considered. We note that the use of Sobol sequences improves the convergence of the MC integrals for small system sizes, and thereby reduces the cost of the SCP, but leads to an unfavourable exponential scaling (see SI S1.3) for large $N$.

The reference calculations (QTI) were performed using a combination of a TI from the harmonic reference to the anharmonic potential using classical MD and a quantum TI over mass using PIMD. One should note that the variance of the integrands and the cost of performing one molecular dynamics step is different for classical and quantum MD. Thus the minimum number of force evaluations required to reduce the error to within the tolerance is an optimization problem detailed in the supporting information (see SI S4). Since both the integrals are effective energies and the fluctuations of the potential energy in the canonical ensemble display a $~{}\sim\frac{1}{\sqrt{N}}$ behaviour with respect to its mean, QTI also displays a $~{}\sim\frac{1}{\sqrt{N}}$ scaling behaviour.

In summary, in the limit of small system size, SCP and IMF display the most favourable scaling. The reference technique QTI displays a $\mathcal{O}(1/\sqrt{N})$ behaviour, making it the least expensive method in the limit of large system size. For carbon, however, as well as for all the system discussed in this work, we do not reach this limit, and QTI requires a substantially larger number of force evaluations than either SCP or IMF. It is worth noting that at fixed cell size (unlike all other methods) IMF and VSCF also provide the temperature dependence of the free energy without any additional force evaluations.

In the current implementation, none of the above free energy methods exploit crystal symmetries. Exploiting crystal symmetries in HAR, IMF and VSCF is straightforward and the reduction in computational cost is simply related to the reduction in the number of independent normal modes. Crystal symmetries can similarly, albeit not quite as trivially, be exploited in SCP 62 and the other methods. However, crystal symmetries do not affect the overall scaling behaviour with respect to the number of degrees of freedom considered in a given calculation and have therefore not been considered in the benchmarks for the computational cost.

4.3 Accuracy

We gauge the accuracy of the approximate methods by studying the error incurred in the absolute anharmonic free energy, and in the free energy differences between two phases of the same material.

4.3.1 Allotropes of Carbon

Diamond and lonsdaleite are mildly anharmonic systems which serve as excellent starting points for our study. We find that 32-atom simulation cells suffice to converge the free energy difference between diamond and lonsdaleite, with respect to BZ sampling. A detailed description of the workflow and system size convergence can be found in the SI (see S1.1).

As shown in panels (a) and (b) of Fig. 5, the quantum anharmonic contribution to the free energy of both diamond and lonsdaleite is approximately 2 meV/atom. IMF, which considers anharmonicity only along normal modes, underestimates the anharmonic free energy by around 1 meV/atom. Including pairwise mean-field coupling using VSCF leads to a large over correction that increases the error to over 2 meV/atom, while SCP (which also includes a mean-field anharmonic corrections within Gaussian statistics) gives excellent results in comparison to the reference. This indicates that the error in VSCF arises from the truncation of the potential.

We also study the accuracy of the methods in reproducing the free energy difference between diamond and lonsdaleite, as shown in panel (c) of Fig. 5. Notably, the free energy contribution from quantum anharmonic motion for the two allotropes is almost identical so that the anharmonic free energy difference is almost identical to the harmonic free energy difference. Fortuitously IMF and VSCF benefit from large amounts of error cancellation and reproduce the exact result within the errors in the anharmonic free energies. Overall, all approximate methods perform reasonably well at reproducing both the (very small) anharmonic corrections and the free energy difference.

4.3.2 Polymorphs of Ice

Ice XIh and ice XIc are a more challenging test-case because of the simultaneous large anharmonic free energy due to the pronounced anharmonicity of the O–H bond, the coupling between the high and low frequency modes79, and the small free energy difference between the hexagonal and cubic polymorphs 25, 26. Supercells containing 16 molecules of water suffice to converge the free energy difference for all methods. Details of the calculations can be found in the SI (see S2.1)

As shown in panels (a) and (b) of Fig. 6, the overall contribution from quantum anharmonicity to the free energy is around 25 meV/molecule for both systems. Contrary to the case of carbon, we find the approximate methods do not accurately reproduce the reference anharmonic free energy. For instance, the IMF technique produces qualitatively incorrect anharmonic corrections. The VSCF approach with pairwise couplings of normal modes provides the best approximation, but remains off by over 10 mev/molecule. The SCP scheme incurs errors of around 20 meV/molecule.

In line with previous path integral calculations on hexagonal and cubic ice26, we find the free energy difference between the polymorphs of ice XI to be almost zero, as shown in panel (c) of Fig. 6. IMF predicts the hexagonal form to be more stable by around 7 meV/molecule. After adding mean field coupling corrections within VSCF the margin of stability reduces to around 5 meV/molecule. The SCP scheme benefits from cancellation of errors and fortuitously gives the correct result within 1 mev/molecule.

4.3.3 Polymorphs of Paracetamol

As a final test, we consider forms 1 and 2 of crystalline paracetamol. These are more complex molecular crystals, for which free energy calculations are complicated by the presence of quasi-free rotations of the methyl groups. Reference free energies are obtained by classical TI with respect to temperature 51 from 10 K to 300 K followed by quantum TI over masses. To ensure that all of the three degenerate rotational conformers of the paracetamol molecules, are explored evenly at low temperatures, the classical TI is performed using replica exchange molecular dynamics 80. At 10 K classical anharmonicity is almost completely suppressed and the classical anharmonic free energy can be accurately estimated using a harmonic approximation (accounting for the three degenerate conformers). For reference, the free energies were recalculated using the TI route employed in Ref. 18 which proves to reproduce the same result.

As shown in Figs. 7 (a) and (b), the overall quantum anharmonic corrections for forms I and II are around -58 and -46 meV/molecule. The degeneracy of the conformers contributes the dominant part, $-k_{B}T\ln(3)\approx-28$ meV / molecule. All approximate anharmonic methods produce qualitatively incorrect free energy corrections due to anharmonicity. Furthermore, just as the harmonic approximation, they are unable to account for the degeneracy of the three conformers, that would have to be identified and corrected for manually.

As shown in Fig. 7 (c), the free energy difference between the two forms is around 12 meV/molecule. The difference with respect to Ref.18 arises due to the use of slightly different lattice constants (see SI S3.1.4), a more accurate path integral sampling technique73 and a finer PPPM mesh for the Ewald summation of electrostatics. As in the cases of carbon and ice, the IMF and SCP benefit from significant error cancellation. The former correctly predicts form II to be more stable but performs worse than a harmonic approximation in getting the correct magnitude. The latter also predicts the correct sign but fortuitously estimates the magnitude to within 5 mev/molecule of the exact result. VSCF doesn’t benefit from error cancellation to the same extent and overestimates the stability of form I by over 80 meV/molecule.

The failure of normal mode based approaches for paracetamol is unsurprising, as the description of quasi-free rotations requires curvilinear coordinates. In paracetamol the potential energy barrier for rotational motion correspond to approximately 200 K$k_{B}$ implies (even classically) quasi-free rotation of the methyl groups at room temperature. For the force field used to describe paracetamol, the potential governing rotation and breathing of methyl groups can be extracted explicitly (neglecting coupling to the remainder of the molecule) and takes the simple form

[TABLE]

where $k=53.114$ eV/Å2, $r_{0}=1.09$ Å, and $V_{\theta}=8.54$ meV. This allows us to study the failure of the above approximate free energy methods with the activation of angular motion on the basis of a simple toy model,

[TABLE]

which can easily be studied over a range of temperatures. The exact solution for this simple model is obtained by exact diagonalization (ED) of the Hamiltonian matrix on a regular, square two dimensional real space grid of $256\times 256$ points spanning from $(x,y)=(-1.5r_{0},-1.5r_{0})$ to $(x,y)=(1.5r_{0},1.5r_{0})$. We find that the reference free energy is converged to within 0.2 meV.

The temperature dependence of the free energy of the model system is shown in Fig. 8. Since none of the approximate methods are able to account for the three fold degeneracy of the system, we remove $-k_{B}T\ln{3}$ from the exact results and study the anharmonic free energy of one of the minima. The results show that above around 75 K the VSCF approximation becomes increasingly inaccurate as the amplitude of angular motion of methyl groups increases and the vibrational density delocalizes over the three equivalent potential energy minima (see Fig. 8). Furthermore, the harmonic, IMF, and SCP approximations severely overestimate of the free energy throughout. For the harmonic and IMF approximation this can simply be explained by the fact that linear coordinates mix angular and much higher frequency radial motion, so that the effective “angular mode” is stiffened substantially, while the radial mode retains the true harmonic frequency.

5 Conclusions

Diamond and lonsdaleite, as examples of simple quasi-harmonic solids, highlight the utility of approximate free energy methods. While the accuracy of the approximate Helmholtz free energies varies, all approaches achieve sub-2 meV/atom accuracy and, more importantly, consistently yield a systematic improvement over the harmonic approximation at a substantially lower computational cost than the reference QTI. On the other hand, ice and paracetamol, as examples of more complex, molecular crystals, highlight the limitations of approximate techniques. The free energies of the molecular crystals are substantially overestimated due to the inherent limitations of normal modes based descriptions in the presence of large-amplitude curvilinear librational or quasi-free rotational motion. The simple model description of the rotation of the methyl group in paracetamol demonstrates that SCP, IMF and VSCF artificially stiffen the rotational modes, leading to an overestimation of the total free energy. The failure of these methods for ice can also be understood along the same lines: at larger displacements along the normal modes initially corresponding to librational motion, O–H bonds are stretched and bent, leading to an overestimation of the effective frequency of librational motion and consequently the free energy. This is confirmed by the blue shifts of the librational modes observed in the case of IMF and SCP with respect to the harmonic approximation (see SI S2.3). Consequently, these methods do not consistently yield systematic improvements over the harmonic approximation for free energies of solids that possess high amplitude librational or quasi-rotational modes. However, we expect these techniques to perform well for atomic and ionic solids, where the point-like nature of the atomic/ionic building blocks suppresses large-amplitude curvilinear motion. It is worth mentioning that (in suitable applications) normal modes based approximate methods lend themselves to identifying the atomistic/structural origins of anharmonicity and facilitate analyses of, for example, spectral properties of strongly anharmonic phonons, as probed by inelastic scattering processes, the formation charge-density-waves, and ferroelectric instabilities 81.

We also find that the approximate methods methods benefit from error cancellation, leading to errors in free energy differences that are consistently smaller than the errors in the absolute Helmholtz free energies. However, it is worth noting that we have compared systems with very similar local environments. In general such beneficial cancellation of errors is not guaranteed. We demonstrate this in section S5 of the SI by studying the free energy of a few high density phases of ice (II, IX and XV) relative to that of XIh. As shown in in Fig. S5, SCP – that benefits from error cancellation when comparing the cubic and hexagonal forms of ice – does not yield a qualitatively correct order of relative free energies.

As the approximate results can vary from almost quantitatively accurate to qualitatively incorrect results, QTI is the only free energy methods among those considered in this work that provides reliable anharmonic free energies for large and complex organic solids. Given that it displays a $\mathcal{O}(1/\sqrt{N})$ computational cost, it may furthermore require comparable or fewer force evaluations than a SCP or VSCF calculation for systems of interest, in particular when considering biological or pharmaceutical compounds that involve large unit cells with flexible molecular units. It is further worth noting that QTI (and other statistical sampling methods) are substantially less susceptible to noise in the underlying PES than the harmonic approximation, IMF and VSCF. While random noise largely cancels out in the ensemble averages calculated in statistical sampling methods, especially in combination with stochastic thermostats 82, 83, the harmonic approximation relies on the ability to determine a meaningful dynamical matrix and thus a differentiable PES, and the IMF and VSCF require an interpolatable PES. This is demonstrated by performing the Ewald summation in the description of paracetamol using a coarser PPPM mesh, which leads to discontinuities in the PES of 0.50 meV (see SI S6). While the resultant QTI and SCP free energies remain largely unaffected, the free energy estimates obtained from the analytic methods – at least in the implementation we discuss here – depend heavily on the size of the finite displacements underlying the mapping of the PES and the VSCF in particular eventually fails to converge altogether.

Efforts towards obtaining a computationally feasible anharmonic free energy should therefore be channelled towards reducing the cost of performing a QTI or, at least, a classical-nuclei TI with nuclear quantum contributions evaluated at a more approximate level 18, 28. This includes streamlining hierarchical frameworks84 that perform the full free energy calculations using inexpensive bespoke potentials85 or cheaper basis sets12, reducing the cost to that of reversibly switching an ab initio potential. Machine learning potentials offer exceptional promise to provide ab-initio-quality potential energy surfaces to evaluate the anharmonic free energy 26, and approximate methods could also constitute an effective sampling approach to generate data to train and validate such ML potentials.

6 Acknowledgements

VK and EAE contributed equally to the work. We would like to thank Francesco Mauri, Marco Cherubini for helping us validate our SCP implementation. Bartomeu Monserrat, is also acknowledged for providing a critical reading of an early version of the manuscript. MC and EAE, acknowledge funding by the European Research Council under the European Union’s Horizon 2020 research and innovation program (grant no. 677013-HBMAP). MR acknowledges funding from BiGmax, the Max Planck research network on big-data-driven materials science.

7 Supporting Information

Details of calculations on allotropes of carbon, ice XI and paracetamol; Analysis of efficiency of the thermodynamic integration route; Relative stability of some of the high density proton ordered phases of ice; Influence of noise in the potential energy landscapes on IMF