Excess Chemical Potential via the Widom Method

We recall that the free energy of the canonical ensemble, termed the Helmholtz free energy and denoted $F$, is defined by

$$F = -k_BT \ln Q\left(N,V,T\right)$$ (308)

$$= -k_BT \ln\left(\left\{\frac{V^N}{\Lambda^{3N}N!}\right\} \left\{V... ...bf r}^N \exp\left[-\beta\mathscr{U}\left({\bf r}^N\right)\right]\right\}\right)$$ (309)

$$= -k_BT\ln\left(\frac{V^N}{\Lambda^{3N}N!}\right) - k_BT\ln\left(\int d{\bf s}^N \exp\left[-\beta\mathscr{U}\left({\bf s}^N;L\right)\right]\right)$$ (310)

$$\equiv F_{\rm id}\left(N,V,T\right) + F_{\rm ex}\left(N,V,T\right)$$ (311)

Here, $F_{\rm id}$ is the “ideal gas” free energy, and $F_{\rm ex}$ is the “excess” free energy. The chemical potential is defined as the change in free energy upon addition of a particle:

$$\mu = \left(\frac{\partial F}{\partial N}\right)_{VT}$$ (312)

For large $N$,

$$\mu = -k_BT \ln\left(Q_{N+1}/Q_N\right)$$ (313)

$$= -k_BT\ln\left(\frac{V/\Lambda^d}{N+1}\right) - k_BT\ln\left( \fra... ...cr{U}\left({\bf s}^{N};L\right)\right] \end{displaymath}\end{minipage}}}\right)$$ (314)

$$= \mu_{\rm id}\left(\rho\right) + \mu_{\rm ex}$$ (315)

which defines the excess chemical potential, $\mu _{\rm ex}$. We see that we can express the quotient of configurational integrals in $\mu _{\rm ex}$ as an integration of the ensemble average of $\Delta\mathscr{U} \equiv \mathscr{U}\left({\bf s}^{N+1}\right) - \mathscr{U}\left({\bf s}^N\right)$ over ${\bf s}_{N+1}$, the scaled coordinates of the $(N+1)$'th particle, or “test” particle:

$$\mu_{\rm ex} = -k_BT \ln \int d{\bf s}_{N+1}\left<\exp\left(-\beta\Delta\mathscr{U}\right)\right>_N$$ (316)

This equation implies that we can measure $\mu _{\rm ex}$ by performing a brute force sampling of $\exp\left(-\beta\Delta\mathscr{U}\right)$ in an otherwise normal NVT MC simulation. That is, we can at frequent intervals in a normal MC program “create” a test particle with a position sampled uniformly in the box, compute $\mathscr{U}\left({\bf s}^{N+1}\right) - \mathscr{U}\left({\bf s}^N\right)$, and accumulate $\exp\left(-\beta\Delta\mathscr{U}\right)$. This is the Widom method.

The code mclj_widom.c implements the Widom method for the Lennard-Jones fluid in an NVT simulation. Below is a code fragment for sampling $\Delta\mathscr{U}$ using the Lennard-Jones pair potential 89:

  rx[N]=(gsl_rng_uniform(r)-0.5)*L;
  ry[N]=(gsl_rng_uniform(r)-0.5)*L;
  rz[N]=(gsl_rng_uniform(r)-0.5)*L;

  for (j=0;j<N;j++) {
    dx  = rx[N]-rx[j];
    dy  = ry[N]-ry[j];
    dz  = rz[N]-rz[j];
    r2  = dx*dx + dy*dy + dz*dz;
    r6i = 1.0/(r2*r2*r2);
    du += 4*(r6i*r6i - r6i);
  }

The particle with index $N$ is assumed to be the “test particle”; the other particles are indexed 0 to $N-1$. In the first bit, the position of the test particle is generated as a uniformly random location inside a cubic box of side length $L$. Then we loop over the particles 0 to $N-1$ and accumulate $\Delta\mathscr{U}$.

Using the code mclj_widom.c, we can measure $\mu_{\rm ex}\left(\rho,T\right)$ in NVT MC simulations. This represents an alternate way of computing $\mu _{\rm ex}$ to that of $\mu VT$ MC, in which $\rho$ is an observable. In Fig. 39, I show $\mu _{\rm ex}$ vs. $\rho$ at $T$ = 3.0 computed using both $\mu VT$ MC and the NVT Widom method. The $\mu VT$ simulations were initialized with 512 particles initially, while the Widom simulations were run with $N$ = 216 particles. Each point is an average of three indenpemdent calculations.

Figure 39: $\mu ^{\rm ex}$ vs. $\rho$ for the Lennard-Jones fluid at $T$ = 3.0 computed using a grand canonical $\mu$VT Monte Carlo simulation and an NVT simulation with the Widom sampling method.

It would be useful to know how to determine which of these apparently competing methods is best for computing $\mu _{\rm ex}$. They are both similar in computational requirements (this is not further qualified here; if someone wants to make this comparison, he or she is welcome to do this as a project). On the one hand, we have an inherent limitation of the grand canonical simulation: one cannot specify the system density exactly; rather it is an observable with some mean and fluctuations. The Widom method does allow one to specify the density precisely, and in this regard, it is probably more trustworthy in computing $\mu _{\rm ex}$. On the other hand, the Widom method suffers the limitation that it is not generally applicable to systems with any potential energy function. For example, for hard-sphere systems, the Widom method would always predict that $\mu _{\rm ex}$ is 0, a clearly nonsensical answer. The “overlapping distribution method” of Bennett, discussed in Section 7.2.3 of F&S, offers a means to overcome this particular limitation. We do not cover this method in lecture, but you are encouraged to explore the overlapping distribution method on your own (maybe as a project).