How to simulate a grand-canonical ensemble?¶
Metropolis-Hastings Algorithm¶
The grand-canonical distribution (see eq. (1)) can be sampled using the Metropolis Monte Carlo method, similiar to the canonical ensemble. Before we derive the acceptance criterion for the GCMC method, let us shortly recall the generic Metropolis-Hastings algorithm.
We want to calculate averages of observables $A$ which are distributed according to some distribution $p^\mathrm{eq}_s$ (e.g. the canonical or grand-canonical distribution)
$$
\langle A\rangle = \frac{\sum_{s\in\Gamma} A_s p^\mathrm{eq}_s}{\sum_{s\in\Gamma} p^\mathrm{eq}_s}. \quad \tag{2}
\label{eq:can_average_mc}
$$
Here, $s$ labels the different states of the systems, $A_s$ is the value of the observable $A$ in the state $s$ and the sum runs over the set of all states $\Gamma$.
The sum can also be an integral for systems with continuous degrees of freedom, for this derivation it does not matter.
For almost all systems, the ensemble average (2) can not be evaluated exactly and only numerical approximations can be obtained.
The basic idea of all MC methods is to randomly sample only a subset of states.
In the most naive approach (simple sampling), states are simply selected with a uniform probability from the space $\Gamma$ of states $s$, resulting in a sample $\gamma\subset\Gamma$.
The canonical average can then be approximated by a sum over $\gamma$:
$$
\langle A\rangle \approx \frac{\sum_{s\in\gamma} A_s p^\mathrm{eq}_s}{\sum_{s\in\gamma} p^\mathrm{eq}_s} \quad \tag{3}
$$
This method is highly inefficient, as most states will typically not contribute significantly to the average.
To get a more efficient sampling method, let us select the states $s$ according to a non-uniform probability distribution $p_s$.
Then, the average can be written as
$$
\langle A\rangle^{\mathrm{can}} \approx \frac{\sum_{s\in\gamma} A_s p^\mathrm{eq}_s/p_s}{\sum_{s\in\gamma} p^\mathrm{eq}_s/p_s}, \quad \tag{4}
\label{eq:can_average_mc_prob}
$$
where we have to divide in the numerator and in the denominator by $p_s$ to compensate for the fact that the samples are distributed non-uniformly.
It is obvious that we recover (3) from (4) as the special case of a uniform distribution.
To sample the averages more efficiently, we simply pick
$$
p_s \propto p^\mathrm{eq}_s \quad \tag{5}
\label{eq:mc_prob}
$$
such that the selected samples are representative of the desired equilibrium distribution.
In this case, the estimators of the averages are thus simply given by
$$
\langle A\rangle \approx \frac{1}{N}\sum_{s\in\gamma}A_s.
$$
The remaining problem is now to find a method which allows us to sample states according to (5).
This can be done in terms of so-called Markov chains.
A Markov chain is a stochastic process without memory.
In our case, the Markov chain is simply a sequence $\{s_i\}_{i=1,...N}$ of $N$ states $s_i$ of the system and the fact that the Markov chain has no memory means that the transition probability from a state $s_i$ to a state $s_j$ can be given in terms of a transition matrix $W_{i\rightarrow j}$ which is constant.
The equilibrium probability distribution has to be invariant under $W$. Mathematically this condition can be expressed in the form
$$
\sum_{i}W_{i\rightarrow j}p_i^\mathrm{eq} = p_j^\mathrm{eq}
$$
A sufficient constraint on $W_{i\rightarrow j}$ which ensures the validity of this condition
is the so-called criterion of detailed balance:
$$
\frac{W_{i\rightarrow j}}{W_{j\rightarrow i}} = \frac{p_j^\mathrm{eq}} {p_i^\mathrm{eq}}.
$$
Typically, the transition matrix elements are factorized into a proposal probability $g_{i\rightarrow j}$ and acceptance probability $P_{i\rightarrow j}^\mathrm{acc}$:
$$
W_{i\rightarrow j} = g_{i\rightarrow j} P_{i\rightarrow j}^\mathrm{acc}.
$$
For the symmetric choice $g_{i\rightarrow j}=g_{j\rightarrow i}$, the condition of detailed balance reduces to
$$
\frac{P_{i\rightarrow j}^\mathrm{acc}}{P_{j\rightarrow i}^\mathrm{acc}} = \frac{p_j^\mathrm{eq}}{p_i^\mathrm{eq}}.
\quad \tag{6}
$$
By checking the cases $p_j^\mathrm{eq}>p_i^\mathrm{eq}$ and $p_j^\mathrm{eq} \leq p_i^\mathrm{eq}$ one can see that the acceptance probability
$$
P_{i\rightarrow j}^\mathrm{acc,\text{ }Metropolis} = \min\left(1,\frac{p_j^\mathrm{eq}}{p_i^\mathrm{eq}}\right) \quad \tag{7}
$$
fulfills equation (6). With this, we can now formulate the Metropolis-Hastings algorithm [2].
We start with an arbitrary initial state $s_1$ and generate a sequence $\{s_i\}_{i=1,...N}$ of $N$ states.
Each new state (n) is generated from the previous state (o) in the following way.
First, we propose a new state according to the symmetric proposal probability $g_{\mathrm{o}\rightarrow\mathrm{n}}$.
Then, we accept the new configuration with a probability given by equation (7).
If the proposed new state is rejected, the old state is retained.
Applying this algorithm many times results in the desired Markov chain of representative configurations.
GCMC method¶
To simulate a system in a grand-canonical ensemble the described algorithm can be applied the grand-canonical distribution [1]. When considering the possible states of the system, it becomes clear that every state can be reached by a combination of two kinds of moves:
Displacement moves that leave the particle numbers fixed and only change positions. This is equivalent to sampling in the canonical ensemble with the current particle number. Due to the ergodic hypothesis this operation can be performed either using canonical Monte Carlo methods or canonical molecular dynamic simulations (e. g. Langevin dynamics). In this tutorial we choose the latter.
Insertion and Deletion moves that add or delete a single particle while keeping all other particle positions fixed.
In the GCMC method the insertion and deletion moves are implemented in the following fashion. First of all it is determined randomly with equal probability if an insertion or deletion move is performed. In the case of an insertion the additional particle is placed at a uniformly chosen random position in the simulation box. Then the proposed new state is accepted according to the criterion
$$
\begin{split}
P_{N_i\rightarrow N_i+1}^\mathrm{acc,\text{ }GCMC} &= \min\left(1,\frac{1}{N_i+1}\left(\frac{V}{\Lambda_i^3}\right)\times \exp\left(-\beta \left(U^{N_i+1}(\mathbf{q})-U^{N_i}(\mathbf{q})\right)+\beta \mu_i\right)\right).
\end{split}
$$
In the case of a deletion a randomly selected particle is removed from the system.
We have the analogous result
$$
\begin{split}
P_{N_i\rightarrow N_i-1}^\mathrm{acc,\text{ }GCMC} &= \min\left(1,N_i\left(\frac{\Lambda_i^3}{V}\right)\times \exp\left(-\beta \left(U^{N_i-1}(\mathbf{q})-U^{N_i}(\mathbf{q})\right)-\beta \mu_i\right)\right).
\end{split}
$$
