Widom Particle Insertion Method¶
In a pure MD simulation, there is no direct way to access the chemical potential, since it is related to changes in the particle number that is typically fixed.
In 1963, Widom proposed a MC scheme to measure the chemical potential in a system using trial particle insertions [1].
The chemical potential for species $\alpha$ is defined as
\begin{equation}
\mu_{\alpha} = \left(\frac{\partial G}{\partial N_\alpha} \right)_{p,T,N_{\beta\neq\alpha}} = \left(\frac{\partial F}{\partial N_\alpha}\right)_{V,T,N_{\beta\neq\alpha}} \quad ,\tag{1}
\label{eq:chem_pot_def}
\end{equation}
and corresponds to the difference in the Helmholtz free energy $F$ or Gibbs free energy $G$ when adding or removing one particle in the system.
The canonical partition function is defined as
\begin{equation}
Z(N,V,T)=\frac{V^N}{\Lambda^{3N}N!} \int_{\left[0,1\right]^{3N}}\exp\left(-\beta U\left(\left\{\boldsymbol{s}^i\right\}\right)\right)\,\mathrm{d}\boldsymbol{s}^1\ldots\mathrm{d}\boldsymbol{s}^N \quad ,\tag{2}
\label{eq:partition_function}
\end{equation}
where V is the volume, $\Lambda=h/\sqrt{2\pi m k_\mathrm{B}T}$ the thermal wavelength, $\beta$ the inverse temperature, $U$ the interaction potential and $\boldsymbol{s}^i$ the rescaled (dimensionless) particle coordinates.
The Helmholtz free energy $F$ can be computed from the partition function: $F=-k_BT\ln\left(Z(N,V,T)\right)$.
For sufficiently large N, the following relation approximates the chemical potential:
\begin{align}
\mu & = -k_\mathrm{B}T \ln\left(\frac{Z(N+1,V,T)}{Z(N,V,T)}\right) \label{eq:mu_limit_largeN}\\
& = \mu_\textrm{ideal}+\mu_\textrm{ex} \nonumber\\
& = - k_\mathrm{B}T\ln\left(\frac{V}{\Lambda^{3}(N+1)}\right)+\mu_\textrm{ex}\, , \nonumber\tag{3}
\end{align}
where $\mu_\textrm{ideal}$ depends only on the density of the system, so it is a constant for systems with different interactions but the same particle density.
The interesting part is the excess chemical potential $\mu_\textrm{ex}$, which can be expressed in a form involving a canonical average:
\begin{equation}
\mu-\mu_\textrm{ideal}=\mu_\textrm{ex}=-k_\mathrm{B}T\ln\left(\int_{\left[0,1\right]^{3}}\langle\exp(-\beta\Delta U)\rangle\,\mathrm{d}\boldsymbol{s}^{N+1}\right)\quad .\tag{4}
\label{eq:mu_ex}
\end{equation}
In the above expression, $\Delta U$ is the energy difference of the system before and after the insertion of an additional particle ("test particle") into the system, and $\langle...\rangle$ is the canonical average for a system with $N$ particles.
The expression thus suggests that the excess chemical potential can be computed using the following procedure.
The canonical average $\langle\exp(-\beta\Delta U)\rangle$ can be computed by inserting a test particle into many different configurations of the system (sampled according to the canonical distribution) and then averaging the associated Boltzmann factor $\exp(-\beta\Delta U)$.
Note that this canonical average can be computed using either MC or canonical MD, in the following we will use Langevin MD for this task.
The integral over $\boldsymbol{s}^{N+1}$ corresponds to an average over different insertion positions of the test particle and can be calculated using a brute-force MC approach by simply inserting test particles at many different (random) trial positions for each sampled state of the system.
