Thermodynamic Integration

Thermodynamic integration is a conceptually simple, albeit expensive, way to calculate free energy differences from MC or MD simulations. In this example, we will consider the calculation (again) of chemical potential in a Lennard-Jones fluid at a given temperature and density, a task performed very well already by the Widom method (so long as the densities are not too high.) More details of the method can be found in the work of Tironi and van Gunsteren [43].

We begin with the relation derived in the book for a free energy difference, $F_{II} - F_{I}$, between two systems which are identical (same number of particles, density, temperature, etc.) except that they obey two different potentials. System I obeys $\mathscr{U}_I$ and System II $\mathscr{U}_{II}$. To measure this free energy difference, we must integrate along a reversible path from I to II. So let us suppose that we can write a “metapotential” that uses a switching parameter, $\lambda$, to measure distance along this path. So, when $\lambda = 0$, we are in System I, and when $\lambda = 1$ we are in System II. One way we might encode this (though this is not necessarily a general splitting, as we shall see below) is

$$\mathscr{U}\left(\lambda\right) = \left(1-\lambda\right)\mathscr{U}_{I} + \lambda\mathscr{U}_{II}$$ (317)

Let us consider the canonical partition function for a system obeying a general potential $\mathscr{U}\left(\lambda\right)$:

$$Q\left(N,V,T,\lambda\right) = \frac{1}{\Lambda^{3N}N!}\int d{\bf r}^N\exp\left[-\beta\mathscr{U}\left(\lambda\right)\right]$$ (318)

Recalling that the free energy is given by $F = -k_BT\ln Q$, we can express the derivative of the Helmholtz free energy with respect to $\lambda$:

$$\left(\frac{\partial F\left(\lambda\right)}{\partial\lambda}\right)_{N,V,T} = -\frac{1}{\beta}\frac{\partial}{\partial\lambda}\ln Q\left(N,V,T,\lambda\right)$$ (319)

$$= -\frac{1}{\beta Q\left(N,V,T,\lambda\right)}\frac{\partial Q\left(N,V,T,\lambda\right)}{\partial\lambda}$$ (320)

$$= \frac{ \begin{minipage}{6cm} \begin{displaymath} \int d{\bf r}^N ... ...t[-\beta\mathscr{U}\left(\lambda\right)\right] \end{displaymath}\end{minipage}}$$ (321)

The free energy difference between I and II is given by:

$$F_{II} - F_{I} = \int_{\lambda=0}^{\lambda=1}\left<\frac{\partial\mathscr{U}}{\partial\lambda}\right>_\lambda d\lambda$$ (322)

where, $\left<\frac{\partial\mathscr{U}}{\partial\lambda}\right>_\lambda$ is the ensemble average of the derivative of $\mathscr{U}$ with respect to $\lambda$.

To compute $\mu _{\rm ex}$, we imagine two systems: System I has $N-1$ “real” particles, and 1 ideal gas particle, and system II has $N$ real particles. The two free energies can be written:

$$F_{\rm I} = F_{\rm id}\left(N,V,T\right) + F_{\rm ex}\left(N-1,V,T\right)$$ (323)

$$F_{\rm II} = F_{\rm id}\left(N,V,T\right) + F_{\rm ex}\left(N,V,T\right)$$ (324)

For large values of $N$, we see that $\mu_{\rm ex} = F_{\rm II} - F_{\rm I}$. So, we have another route to compute $\mu _{\rm ex}$. First, we tag a particle $i_\lambda$, call it the “$\lambda$-particle”, and apply the following modified potential to its pairwise interactions:

$$U_{\rm LJ,\lambda}\left(r;\lambda\right) = 4\left(\lambda^5r^{-12} - \lambda^3r^{-6}\right)$$ (325)

So, the total potential is given by

$$\mathscr{U} = {\sum_{i<j}} ^\prime U_{\rm LJ}\left(r_{ij}\right) + \sum_i U_{\rm LJ,\lambda}\left(r_{i,i_\lambda};\lambda\right)$$ (326)

where the prime denotes that we ignore particle $i_\lambda$ in the sum. When $\lambda = 1$, all particles interact fully, and we have System II.

Next, we conduct many independent MC simulations at various values of $\lambda$ and a given value of $\rho$ and $T$, generating for each $\left(T,\rho\right)$ a table of $\left<\partial\mathscr{U}/\partial\lambda\right>_\lambda$ vs. $\lambda$ which can be integrated to yield a single value for $\mu _{\rm ex}$. The code mclj_ti.c implements this sampling when a value for $\lambda$ is specified. To demonstrate its use to compute $\mu _{\rm ex}$, the following protocol was used:

1. Run 600,000-cycle NVT $\lambda$-MC simulations at given density $\rho$ and temperature $T$ (with $N$ = 216 particles) for values of $\lambda\in\left\{0,0.1,0.2,\dots,1\right\}$ to compute $\langle d\mathscr{U}/d\lambda\rangle$. Three independent simulations are run for each $N$-$\rho$-$T$-$\lambda$ point, and $\langle d\mathscr{U}/d\lambda\rangle$ is an average over these three.
2. Use Simpson's rule in scipy.integrate.simpson to numerically integrate $\langle d\mathscr{U}/d\lambda\rangle$ vs. $\lambda$ to obtain $\mu _{\rm ex}$.

This turns out to be an expensive way to compute the chemical potential for a Lennard-Jones fluid, compared to the Widom method (Sec. 9.1) or grand canonical MC (Sec. 5.1), for at least low to moderate densities. At very high densities, however, particle insertion moves in grand canonical and Widom-method simulations become difficult. Fig. 40 shows a plot of $\langle d\mathscr {U}/d\lambda$ vs. $\lambda$ for various densities, all at $T$ = 3.0 (left), with values of $\mu _{\rm ex}$ found from integrating those curves shown together with the data from Fig. 39 showing $\mu _{\rm ex}$ from both grand canonical MC and the Widom method.

Figure 40: (left) $\langle d\mathscr {U}/d\lambda$ vs. $\lambda$ for the Lennard-Jones fluid for three values of $\rho$ at $T$ = 3.0, computed from MC simulations of 216 particles for 6$\times$10$^6$ cycles. Five independent simulations were run for each $(\rho ,\lambda )$, and values of $\langle d\mathscr {U}/d\lambda$ are averaged over these five. Standard deviations a smaller than the symbol size. (right) Overlay of $\mu _{\rm ex}$ vs $\rho$ at $T$ = 3.0 using grand canonical MC, the Widom method, and thermodynamic integration.

The curves of $\langle d\mathscr{U}/d\lambda\rangle$ vs $\lambda$ are not completely noise-free, but integrating each of these curves to produce a single value of $\mu _{\rm ex}$ produces values that are not too off from the grand canonical and Widom-method simulations.