Computing the heat of adsorption in a grand canonical Monte Carlo simulation

The isosteric heat of adsorption is a thermodynamic quantity that characterizes the enthalpy change associated with the adsorption of a molecule on a surface. This can be measured in the laboratory in two ways. First, one can use calorimetry to measure the heat released upon the adsorption of a given amount of gas. Second, one can look at the sensitivity of the adsorption isotherms with temperature; the Clausius-Clapeyron equation relates the isosteric heat of adsorption to where on the pressure-temperature plane an adsorbed phase in equilibrium with a gas phase lies (at a given amount of adsorption).

In molecular simulations, we typically simulate the adsorption of molecules on a surface using the Metropolis-Hastings Monte Carlo algorithm in the context of the grand canonical statistical ensemble. In this ensemble, we consider a given volume/area of the adsorbent surface in equlibrium with a bulk gas phase at constant chemical potential and temperature. This is connected to experiment by an equation of state that relates the gas phase pressure in the experiment to the chemical potential of the gas.

During the Monte Carlo simulations of the grand canonical ensemble (GCMC simulations), we compute the energy of the adsorbate molecules, $E$, at each sample. One might postulate that the isosteric heat of adsorption is related to the average value of $E/N$ during the simulation, where $N$ is the number of adsorbates in the system. This is incorrect because it is a biased estimate of the energy per adsorbed molecule; when there are fewer particles in the system, the energy of these particles are weighted more in the average.

Instead, we define the average energy of adsorption $q$ as the difference in the average energy between two systems that differ only by the average number of particles. By making the average number of particles between these two systems differ by one, this is the average energy difference upon the addition of a single particle, the average energy of adsorption:

Here, we are looking at $\langle E\rangle$ as a function of $\langle N\rangle$.

Dividing by one, which is also $N+1-N$, we see that this is approximately the derivative of $\langle E\rangle$ with respect to $\langle N \rangle$:

To get a more useful formula for estimating $q$ from GCMC simulations, we need to bring in the grand canonical partition function $\Xi(\beta, \mu, V)$, where $\beta=1/T$, $\mu$ is the chemical potential of the gas, and $V$ is the volume of adsorbent considered (corresponds to a given area of surface). Although a bit abstract, we can write the partition function as a sum over particles in the system $N$ and microstates $\nu$ with that particle number:

Here, $E_{\nu}$ is the energy of that particular microstate $\nu$.

The probability of a microstate $\nu$ with $N$ particles is then $e^{-\beta E_{\mu} + \beta \mu N} / \Xi$. From this and the expectation value formula from elementary probability theory, we can work out the relationship between the partition function and the average energy and particle number:

These two relations give us a nice way to relate $\langle E\rangle$ and $\langle N\rangle$ through the partition function:

Looking back at our equation for the average energy of adsorption $q$, we bring this quantity into $q$ using the chain rule:

We immediately recognize the derivative on the right from the fluctuation theorem in statistical mechanics. Again, this can be seen from elementary probability theory by taking the derivative of $\langle N\rangle = \frac{\partial \log \Xi}{\partial (\beta \mu)}$ once more and writing out the sum for the partition function.

As for the first derivative, we use the relationship we derived:

Here, we wrote out the expectation of $N$ from elementary probability theory using the probability of observing a given value of $N$, $p(N) = \displaystyle \sum_{\nu} \frac{e^{-\beta E_{\nu} + \beta \mu N}}{\Xi}$.

Actually taking this derivative WRT $\beta$ using the product rule ($\Xi$ is a function of $\beta$ too!),

You may be wondering why we did not consider the $e^{\beta \mu N}$ term during differentiation, but actually we are viewing $\Xi$ as $\Xi(\beta, \beta \mu, V)$ here.

We can recognize the first sum as the expected value of $EN$, $\langle E N \rangle$. The second sum involves the term $\frac{1}{\Xi} \frac{\partial \Xi}{\partial \beta}=-\langle E \rangle$. Thus,

Finally, we arrive at the average energy of adsoption as:

which are all statistical quantities we can easily measure from samples during our Monte Carlo simulations. This formula is fascinating because it relates the average energy of adsorption to fluctuations in the particle number and the covariance of the energy and particle number ($\mathrm{Cov}(E,N)=\langle EN \rangle - \langle E \rangle \langle N \rangle$); this is another fluctuation theorem in statistical mechanics.

The isosteric heat of adsorption follows from $q$ by $q_{st}= -q + RT$ since we need to account for the work to push the gas into the gas phase when it desorbs.

comments powered by Disqus