7.3 The Grand Canonical Ensemble¶
Huang, Statistical Mechanics 2ed, Section 7.3
What if you don't know how many atoms you have?¶
We've been on a journey of letting go.
In chapter 6, we held everything tight: fixed \(N\), fixed \(V\), fixed \(E\). The microcanonical ensemble. The control freak.
In section 7.1, we relaxed our grip on energy. We let the system exchange energy with a heat bath, and the Boltzmann factor fell out. Fixed \(N\), fixed \(V\), fixed \(T\). The canonical ensemble.
Now we take one more step. What if particles can also enter and leave the system? What if \(N\) itself isn't fixed?
This sounds weird at first. In your LAMMPS simulation, you created 108 atoms. They're in the box. You know exactly how many there are. Why would \(N\) fluctuate?
Because not every system works like that. And the grand canonical ensemble is how you handle the ones that don't.
Why should you care?¶
The moment you have an open system (particles can enter or leave), you need the grand canonical ensemble. In simulation, that means:
- Adsorption: Gas molecules entering a zeolite, MOF, or porous material. You don't control how many molecules go in. The chemical potential of the gas reservoir does.
- GCMC simulations: LAMMPS has
fix gcmc, which inserts and deletes particles during a simulation. It's literally sampling the grand canonical ensemble. - Interfaces: A liquid-vapor interface in equilibrium. Molecules evaporate and condense. The number in the liquid phase fluctuates.
- Chemical equilibrium: Reactions that create or destroy molecules. The number of each species isn't fixed.
If you've ever set a chemical potential instead of a particle number, you've been working in the grand canonical ensemble.
MD Connection
In LAMMPS, fix gcmc randomly inserts and deletes atoms during a simulation, accepting or rejecting each move based on the chemical potential you specify. The result is a trajectory that samples the grand canonical (\(\mu VT\)) ensemble. This is the standard method for computing adsorption isotherms: set the chemical potential (related to the gas-phase pressure), run GCMC, and count how many particles end up inside your material.
The bad intuition: "just use NVT with a big box"¶
You might think: if I want to simulate gas adsorbing into a material, I'll just put a huge gas reservoir above the surface. Make the box really big, fill it with gas, and let molecules wander into the material on their own.
That works in principle. But it's painfully slow. Most of the simulation time is wasted on gas-phase molecules bouncing around doing nothing interesting. And you need the reservoir to be large enough that removing a few molecules doesn't change the gas-phase density. That means simulating thousands of atoms when you only care about the few dozen that adsorb.
GCMC is the shortcut. Instead of simulating the reservoir explicitly, you replace it with a mathematical rule: insert particles with probability \(\propto z \cdot e^{-\beta \Delta E}\), delete them with probability \(\propto (1/z) \cdot e^{-\beta \Delta E}\). The reservoir is implicit. The chemical potential controls the flow.
But where does that insertion/deletion rule come from? From the grand canonical ensemble. Let's derive it.
The physical setup: particles wandering in and out¶
The setup is almost identical to section 7.1, with one twist.
Start with a large canonical system (fixed \(N\), \(V\), \(T\)). Now mentally draw a small subvolume \(V_1\) inside it. The rest of the system (volume \(V_2 = V - V_1\)) acts as the reservoir.
In the canonical ensemble, particles move freely throughout the whole box. At any instant, some number \(N_1\) of them happen to be inside \(V_1\), and \(N_2 = N - N_1\) are outside. No wall. No barrier. Particles wander in and out of \(V_1\) all the time.
So \(N_1\) fluctuates. Sometimes more particles drift into \(V_1\), sometimes fewer. The question is: what's the probability of finding \(N_1\) particles in \(V_1\) with specific positions and momenta \((p_1, q_1)\)?
This is the same game we played in 7.1, but now we're asking about both the energy and the number of particles in the subsystem.
Where the fugacity comes from¶
The probability of finding \(N_1\) particles in \(V_1\) with coordinates \((p_1, q_1)\) is proportional to how many ways the rest of the system can arrange itself. The rest has \(N_2 = N - N_1\) particles in volume \(V_2\), so:
The first factor is the Boltzmann weight for the subsystem (with the \(N_1!\) from correct Boltzmann counting). The second factor is the partition function of the reservoir.
Now use \(Q_{N_2} = e^{-\beta A(N_2, V_2, T)}\) and Taylor expand the free energy, since \(N_1 \ll N\) and \(V_1 \ll V\):
What's \(\frac{\partial A}{\partial N}\)? That's the chemical potential \(\mu\). And \(-\frac{\partial A}{\partial V}\)? That's the pressure \(P\). So:
Plug this back in. The ratio of partition functions becomes:
Now define the fugacity:
Don't let that word scare you. The fugacity is just a convenient way to package the chemical potential. High \(\mu\) means high \(z\) means particles are eager to enter the system. Low \(\mu\) means low \(z\) means particles would rather stay in the reservoir.
Putting it all together, the density function for the grand canonical ensemble is:
The Boltzmann factor \(e^{-\beta \mathcal{H}}\) is still here (energy weighting). But now there's a new factor: \(z^N\). More particles means a factor of \(z\) for each one. That's the chemical potential doing its job, controlling how much the system "wants" to have more particles.
Key Insight
The grand canonical ensemble has two control knobs: temperature \(T\) (how generous the heat bath is with energy) and chemical potential \(\mu\) (how generous the particle reservoir is with particles). The Boltzmann factor \(e^{-\beta E}\) handles the energy weighting. The fugacity factor \(z^N = e^{\beta \mu N}\) handles the particle number weighting. Each new ensemble relaxes one more constraint and introduces one more reservoir.
The grand partition function¶
Just like we normalized the Boltzmann factor to get \(Q_N\) in the canonical ensemble, we normalize here by summing over all possible particle numbers:
This is the grand partition function. It's a sum of canonical partition functions, each weighted by \(z^N\). If you know \(Q_N\) for every \(N\), you can build \(\mathcal{Q}\).
And just like \(A = -kT \log Q_N\) gave us all of canonical thermodynamics, the grand partition function gives us:
That's the equation of state, directly. No intermediate steps. The logarithm of the grand partition function is \(PV/kT\).
From \(\mathcal{Q}\) to everything¶
The recipe:
| Quantity | Formula |
|---|---|
| Pressure | \(PV = kT \log \mathcal{Q}\) |
| Average particle number | \(\langle N \rangle = z \dfrac{\partial}{\partial z} \log \mathcal{Q}\) |
| Internal energy | \(U = -\dfrac{\partial}{\partial \beta} \log \mathcal{Q}\) |
| Equation of state | Eliminate \(z\) between \(P\) and \(\langle N \rangle\) |
The workflow: compute \(\mathcal{Q}(z, V, T)\), take derivatives, eliminate \(z\) to get \(P(N, V, T)\). Done.
And you're probably thinking, "this looks harder than the canonical ensemble, why bother?" For most systems, you're right. The canonical ensemble is simpler. But for open systems (where \(N\) fluctuates), the grand canonical ensemble is the natural choice. Trying to force a fixed-\(N\) calculation on an open system is like trying to do NVE when you have a thermostat. You can, but you're fighting the physics.
The ensemble ladder¶
Let's step back and see the pattern. We've built three ensembles, each one relaxing one more constraint:
| Ensemble | Fixed | Fluctuates | Key function | Thermodynamic potential |
|---|---|---|---|---|
| Microcanonical | \(N, V, E\) | nothing | \(\Gamma(E)\) | \(S = k \log \Gamma\) |
| Canonical | \(N, V, T\) | energy | \(Q_N(V,T)\) | \(A = -kT \log Q_N\) |
| Grand canonical | \(\mu, V, T\) | energy + particles | \(\mathcal{Q}(z,V,T)\) | \(PV = kT \log \mathcal{Q}\) |
Each ensemble has a reservoir for the thing that fluctuates. The canonical ensemble has a heat bath (energy reservoir). The grand canonical has a heat bath and a particle reservoir.
And here's the punchline from section 7.2: for macroscopic systems, it doesn't matter which one you use. The fluctuations in energy (\(\propto 1/\sqrt{N}\)) and particle number (\(\propto 1/\sqrt{N}\)) are negligibly small. All three ensembles give the same thermodynamics. They differ only in what's convenient to calculate.
Common Mistake
The chemical potential \(\mu\) is not the same as the energy per particle. For an ideal gas, \(\mu = kT \log(n\lambda^3)\) where \(n = N/V\) is the number density and \(\lambda\) is the thermal de Broglie wavelength. At typical conditions, \(\mu\) is negative (the system prefers to have fewer particles at lower density). This confuses people. A negative \(\mu\) doesn't mean particles are "unwanted." It means the system's free energy decreases when you add a particle at the given \(T\) and \(P\). The sign depends on the balance between energy and entropy.
What this means for your simulation¶
When to use which ensemble¶
Here's a practical cheat sheet:
- NVE (microcanonical): Transport properties (diffusion, viscosity). You don't want a thermostat perturbing the dynamics.
- NVT (canonical): Most equilibrium simulations. Structure, thermodynamics, phase behavior at fixed composition.
- \(\mu\)VT (grand canonical): Adsorption, open systems, equilibrium with a reservoir at known chemical potential. Use
fix gcmcin LAMMPS.
How GCMC works (the short version)¶
In a GCMC simulation, LAMMPS alternates between regular MD steps and Monte Carlo insertion/deletion attempts:
- Insertion: Pick a random position in the box. Try to insert a particle there. Accept with probability \(\propto \frac{zV}{N+1} e^{-\beta \Delta E}\), where \(\Delta E\) is the energy change from inserting the particle.
- Deletion: Pick a random existing particle. Try to delete it. Accept with probability \(\propto \frac{N}{zV} e^{-\beta \Delta E}\).
The fugacity \(z\) (set by you via the chemical potential) controls the balance. High \(z\)? Insertions are favored, more particles accumulate. Low \(z\)? Deletions win, particles leave.
After enough cycles, the system equilibrates to the grand canonical distribution. The average number of particles \(\langle N \rangle\) emerges from the simulation. You don't set it. The chemical potential determines it.
MD Connection
In LAMMPS, a typical GCMC setup looks like fix mygcmc all gcmc 100 100 0 0.0 12345 300.0 -0.5 0.1 full_energy. The key parameter is the chemical potential (here \(-0.5\) eV). LAMMPS converts it to a fugacity internally and performs insertion/deletion attempts every 100 steps. You'll see the particle count fluctuate during the run, eventually settling around the equilibrium value for that \(\mu\) and \(T\).
Takeaway¶
The grand canonical ensemble is what you get when both energy and particles can flow between the system and a reservoir. The fugacity \(z = e^{\mu/kT}\) controls the particle flow, just like the Boltzmann factor \(e^{-E/kT}\) controls the energy flow. The grand partition function \(\mathcal{Q} = \sum z^N Q_N\) gives you \(PV/kT\) directly. Use it when you have an open system, and don't fight the physics by forcing a closed-system ensemble onto a problem that's fundamentally open.
Check Your Understanding
- Look at the GCMC insertion probability: \(\propto \frac{zV}{N+1} e^{-\beta \Delta E}\). There's a sneaky \(\frac{1}{N+1}\) sitting in there. Where does it come from? (Hint: what happens if you don't include it?)
- We've built three ensembles, each one letting go of one more thing. What if you let the volume fluctuate too? Fixed \(\mu\), \(P\), \(T\). What would bounce around? Does this ensemble actually exist, or is it nonsense?
- Someone tells you: "The chemical potential of an ideal gas at room temperature is negative. That means the system doesn't want particles." Explain why they're wrong.