The Method of Overlapping Distributions

One interesting feature of the Widom method is that the only trial move is insertion; however, the free-energy difference between an $N$-particle system and an $N+1$-particle system should not depend on which direction the trial moves take. If we imagine a “Widom real-particle removal” method, we'd write the chemical potential as

$$\mu = +k_BT\ln\left(Q_N/Q_{N+1}\right) = \mu_{\rm id} + k_BT\ln\langle\exp\left(+\beta\Delta\mathscr{U}\right)\rangle$$ (327)

Sampling $\langle\exp\left(+\beta\Delta\mathscr{U}\right)\rangle$ in a straightforward NVT MC simulation won't work, however, because...

There is, however, a right way to use bidirectional energy changes to compute free-energy differences, termed the “overlapping distribution method” and attributed to Bennett [44]. Consider two systems 0 and 1, obeying potentials $\mathscr{U}_0$ and $\mathscr {U}_1$, respectively. Let $q_i$ be the scaled configurational integral of the Boltzmann factor:

$$q_i = \int d{{\boldsymbol s}^N} e^{-\beta\mathscr{U}_i}$$ (328)

We can then express the free energy difference between these systems as (assuming for simplicity they have the same volumes):

$$\beta\Delta F = -\ln\frac{q_1}{q_0}$$ (329)

Consider next we run an NVT MC simulation on $\mathscr {U}_1$ and sample $\Delta\mathscr{U}\equiv\mathscr{U}_1-\mathscr{U}_0$. Formally, the probability density of $\Delta\mathscr{U}$ from this simulation is

$$p_1(\Delta\mathscr{U}) = \frac{\displaystyle \int d{{\boldsymbol s}^N} e^{-\beta\mathscr{U}_1}\delta(\mathscr{U}_1-\mathscr{U}_0-\Delta\mathscr{U})}{q_1}$$ (330)

$$= \frac{\displaystyle \int d{{\boldsymbol s}^N} e^{-\beta(\mathsc... ...+\Delta\mathscr{U})}\delta(\mathscr{U}_1-\mathscr{U}_0-\Delta\mathscr{U})}{q_1}$$ (331)

$$= \frac{q_0}{q_1}e^{-\beta\Delta\mathscr{U}}\frac{1}{q_0}\int d{\... ...}^Ne^{-\beta\mathscr{U}_0}\delta(\mathscr{U}_1-\mathscr{U}_0-\Delta\mathscr{U})$$ (332)

$$= \frac{q_0}{q_1}e^{-\beta\Delta\mathscr{U}}p_0(\Delta\mathscr{U})$$ (333)

In going from Eq. 329 to 330 we have used the fact that the Dirac delta function will only permit one value of $\Delta\mathscr{U}$ to survive integration. Taking the log of both sides gives

$$\ln p_1(\Delta\mathscr{U}) = -\ln\frac{q_1}{q_0} - \beta\Delta\mathscr{U} + \ln p_0(\Delta\mathscr{U})$$ (334)

$$= \beta(\Delta F - \Delta\mathscr{U}) + \ln p_0(\Delta\mathscr{U})$$ (335)

This equation provides a way to estimate $\Delta F$ at any one value of $\Delta\mathscr{U}$ so long as good estimates of both $p_1(\Delta\mathscr{U})$ and $p_0(\Delta\mathscr{U})$ are available. That means we must be able to sample a sufficiently large domain of $\Delta\mathscr{U}$ from both a simulation run on $\mathscr {U}_1$ and another run on $\mathscr{U}_0$. That is, there must be a domain of $\Delta\mathscr{U}$ where $p_1$ and $p_0$ overlap.

Bennett[44] suggests the following transformation of $p_1$ and $p_0$ to permit easy calculation of $\Delta F$. Letting

$$f_0(\Delta\mathscr{U}) \equiv \ln p_0(\Delta\mathscr{U}) - \frac{\beta\Delta\mathscr{U}}{2},$$ (336)

$$f_1(\Delta\mathscr{U}) \equiv \ln p_1(\Delta\mathscr{U}) + \frac{\beta\Delta\mathscr{U}}{2}$$ (337)

gives

$$\beta\Delta F = f_1(\Delta\mathscr{U})-f_0(\Delta\mathscr{U}).$$ (338)

This means we can measure $f_1$ and $f_0$ in separate simulations, and then observe a constant offset between them to be $\beta\Delta F$.

Suppose we now take the example of system 1 with $N$ real particles and system 0 with $N-1$ real particles and one ideal-gas particle. The free-energy change from 0 to 1 is the excess chemical potential (yet again!). Fig. 41 illustrates using Bennett's method to compute $\mu _{\rm ex}$ of the Lennard-Jones fluid at $T$ = 1.2 for a few different densities. For each density, two simulations were run: simulation-0 computes the distribution of $\Delta\mathscr{U}$, the energy associated with converting the ideal-gas particle to a real particle, while simulation-1 computes the same distribution for converting a randomly chosen particle from being an ideal-gas particle to being a real particle. This latter $\Delta\mathscr{U}$ is easily computed using the single-particle energy function e_i. It is important to note that the direction of the $\Delta$ is from ideal-gas to real for both simulations. Note too that since we sample $\Delta{\mathscr{U}}$ for particle insertion in simulation-0, we can just as easily compute the expectation $\langle\exp(-\beta\Delta\mathscr{U})\rangle$ and thereby get a direct estimate of $\beta \mu _{\rm ex}$.

At the moderately low density of $\rho$ = 0.7, we see a clear constant offset $\beta \mu _{\rm ex}$ between $f_0$ and $f_1$. Note clear agreement between the offset over a finite-size domain of $\mathscr\Delta{U}$ and the single-point Widom estimate. For the somewhat higher density of 0.9, the offset is a bit noisier, reflecting somewhat poorer sampling. For the highest density, the sampling in simulation-0 is so poor that it is nearly impossible to detect an overlap domain.

Figure 41: $p_0$ and $p_1$ vs $\Delta\mathscr{U}$ (left), and $f_0$, $f_1$, $f_1$ - $f_0$ vs $\Delta\mathscr{U}$ (right) computed using NVT MC simulations at $T$ = 1.2 of system-0 and system-1, for densities of 0.7 (top), 0.9 (middle), and 1.0 (bottom). Estimates of $\beta \mu _{\rm ex}$ from Widom test-particle insertion are shown as red horizontal lines in the plots on the right.