6 Free Energy Methods

In this lecture, we will consider aspects of computing free energy differences from molecular simulations. These methods include the Widom test particle insertion method, the method of overlapping distributions, umbrella sampling, and the non-equilibrium Jarszinsky relation method. In the interest of time, we will consider only two of these techniques in detail: The Widom method and thermodynamic integration.

6.1 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\) in F&S, is defined by

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: \begin{equation} \mu = \frac {\partial F}{\partial N}_{VT} \end{equation}

For large \(N\),

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: \begin{equation} \mu _{\rm ex} = -k_BT \ln \int d{\bf s}_{N+1}\left \langle \exp \left (-\beta \Delta \mathscr {U}\right )\right \rangle _N \end{equation} 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 (Eq. 83):

  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 labeled \(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 )\). Figure 7.1 in F&S reports results of \(\mu _{\rm ex}\left (\rho ,2.0\right )\) using the Widom method for \(\rho \in \left \{0.4,0.5,0.6,0.7\right \}\), and compares to results from Grand canonical simulations. Just for fun, I repeat this exercise for \(T\) = 3.0. The \(\mu VT\) simulations were carried out using the code of Case Study 9 of F&S, first discussed in Sec. 6.1. Here, simulations of 40,000 MC cycles were performed at each state point for which the ideal gas pressure, \(P_{\rm id}\) of the bath is chosen from 0.016, 0.032, 0.064, 0.128, 0.15, 0.20, 0.25, 0.50, 0.75, 1.0, 1.5, 2.0, 3.0, 4.0, 6.0, and the maximum displacement was 0.33 for all runs. Each run has a 100-cycle equilibration, after which it is sampled every 2 steps. The system contains \(N\) = 108 particles. For the Widom method, I considered densities \(\rho \in \left \{0.4, 0.5, 0.6, 0.7, 0.8\right \}\). The system is equilibrated for 100 MC cycles, and is then sampled every two cycles for another 40,000 cycles, using a constant maximum displacement of 0.4. Below is a plot of \(\mu _{\rm ex}\) vs. \(\rho \) at \(T\) = 3.0:

Clearly, I have not taken pains to make sure the two methods agree exactly. It is indeed unclear from the text what these pains might be. But the results agree relatively well, and we might think that, given more cycles in both cases, or by using the same maximum displacement, we would converge to the same answer using both techniques. (Is this so, or have I left a bug in one of my codes...? Perhaps one of you will find out...)

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) using the code for Case Study 15 from book’s website.

6.2 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 Reference [15].

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 \begin{equation} \mathscr {U}\left (\lambda \right ) = \left (1-\lambda \right )\mathscr {U}_{I} + \lambda \mathscr {U}_{II} \end{equation}

Let us consider the canonical partition function for a system obeying a general potential \(\mathscr {U}\left (\lambda \right )\): \begin{equation} 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 ] \end{equation} 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 \):

The free energy difference between I and II is given by: \begin{equation} F_{II} - F{I} = \int _{\lambda =0}^{\lambda =1}\left \langle \frac {\partial \mathscr {U}}{\partial \lambda }\right \rangle _\lambda d\lambda \end{equation} where, \(\left \langle \frac {\partial \mathscr {U}}{\partial \lambda }\right \rangle _\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:

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: \begin{equation} \label {eq:lambda_lj} U_{\rm LJ,\lambda }\left (r;\lambda \right ) = 4\left (\lambda ^5r^{-12} - \lambda ^3r^{-6}\right ) \end{equation} So, the total potential is given by \begin{equation} \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 ) \end{equation} 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 \langle \partial \mathscr {U}/\partial \lambda \right \rangle _\lambda \) vs. \(\lambda \) which can be integrated to yield a single value for \(\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. 6.1), for at least low to moderate densities.

I have done a rough comparison of the thermodynamic integration method described above to the grand canonical MC simulation technique described in Sec. 5.1. Below is a plot of \(\left \langle \partial \mathscr {U}/\partial \lambda \right \rangle _\lambda \) vs. \(\lambda \) for three densities \(\rho \in \left \{0.5, 0.7, 0.9\right \}\). Each point is computed from a single MC simulation using the code mclj_ti.c. The temperature was \(T\) = 2.0, and run for 10\(^6\) cycles for \(N\) = 216. We see that the data is not terribly smooth; it is not clear how many more cycles would result in smoother data.

However, integration is not too sensitive to these fluctuations. As we see below, 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 simulations.

The grand canonical simulations are of course much cheaper, as it requires only a single MC simulation to give a value of \(\mu _{\rm ex}\). However, thermodynamic integration is a generally applicable technique for computing free energy differences.

Questions:

1. Can we make the agreement better by running more MC cycles? How about by using more values of \(\lambda \)?
2. How does this compare to the Widom method?