Thermodynamic integration

Thermodynamic integration is a general scheme to compute the free energy difference F - F₀ between two systems whose potential energies are U and U₀ respectively. The free energy difference F - F₀ is the reversible work done when the potential energy function U₀ is continuously and reversibly switched to U. To do this switching, a continuously variable energy function $U_\lambda$ is defined such that for $\lambda=0, U_\lambda=U_0$ and for $\lambda=1, U_\lambda=U$. We also require $U_\lambda$ to be differentiable with respect to $\lambda$ for $0 \le \lambda \le 1$. A convenient form is:

$$U_\lambda = ( 1 - f(\lambda) ) U_0 + f(\lambda) U,$$ (1)

where $f(\lambda)$ is an arbitrary continuous and differentiable function of $\lambda$ with the property f(0)=0 and f(1)=1. The Helmholtz free energy of this hybrid system is:

$$F_\lambda = - k_{\rm B} T \ln \left\{ \frac{1}{N ! \Lambda^{3... ...ta U_\lambda ( {\bf R}_1 , \ldots {\bf R}_N ; T )} \right\},$$ (2)

Differentiating this with respect to $\lambda$ gives:

$$\frac{dF_\lambda }{d\lambda }=-k_BT\frac{\frac {1}{N!\Lambda^... ...artial U_\lambda }{\partial \lambda } \right\rangle_\lambda,$$ (3)

so

$$\Delta F=F-F_0=\int\limits_0^1d\lambda \left\langle \frac{\partial U_\lambda }{\partial \lambda }\right\rangle _\lambda .$$ (4)

A simple way of defining $U_\lambda$ is:

$$U_\lambda =\left( 1-\lambda \right) U_0 + \lambda U.$$ (5)

Differentiating U$_\lambda$ with respect to $\lambda$ and substituting into Equation 4 yields:

$$\Delta F=\int\limits_0^1d\lambda \left\langle U-U_0\right\rangle _\lambda.$$ (6)

Under the ergodicity hypothesis, thermal averages are equivalent to time averages, so we can calculate $\left\langle \cdot \right\rangle _\lambda$ using molecular dynamics (MD), taking averages over time, with the evolution of the system determined by the potential energy function $U_\lambda$. It is important to stress that the choice of the reference system does not affect the final answer for F, although it does affect the efficiency of the calculations. The latter can be understood by analysing the quantity $\langle U - U_0 \rangle_\lambda$. If this difference has large fluctuations then one would need very long simulations to calculate the average value to a sufficient statistical accuracy. Moreover, for an unwise choice of U₀ the quantity $\langle U - U_0 \rangle_\lambda$ may strongly depend on $\lambda$ so that one would need a large number of calculations at different $\lambda$'s in order to compute the integral in Eq. 6 with sufficient accuracy. It is crucial, therefore, to find a good reference system, where "good" means a system for which the fluctuations of U - U₀ are as small as possible. In fact, if the fluctuations are small enough, we can simply write $F - F_0 \simeq \langle U - U_0 \rangle_0$, with the average taken in the reference ensemble. If this is not good enough, the next approximation is readily shown to be:

$$F - F_0 \simeq \langle U - U_0 \rangle_0 - \frac{1}{2 k_{\rm... ...0 - \langle U - U_0 \rangle_0 \right]^2 \right\rangle_0.$$ (7)

This form is particularly convenient since one only needs to sample the phase space with the reference system.

In our case F is the free energy of the QMC system, and F₀ that of the DFT system, which acts as reference system. Presumably, DFT is already a "good" reference system, in the sense defined above, so that we will be able to use Eq. 7 for the difference F - F₀. The procedure will be to sample the phase space with DFT-MD, and calculate QMC energies at a number of statistically independent configurations extracted from the this simulation. This will be repeated for both solid and liquid.