# Algorithm for Grand Canonical Monte Carlo (GCMC) Molecular Simulation - Initialization

I am trying to make sense of the GCMC algorithm in Frenkel and Smit. I do not understand the conditional statement if a particle is to be displaced or exchanged(removed or added) in a Monte Carlo cycle.

Algorithm 12 in Frenkel and Smit, PER CYCLE:

ran = int(ranf()*(npav + nexc) + 1)
IF (ran < number_of_particles) THEN
'displace a particle'
ELSE
END


It is mentioned in the book that:

Per cycle we perform on average npav attempts to displace particles and nexc attempts to exchange particles with the reservoir.

How large should the numbers npav and nexc be? From my understanding, since for GCMC we only impose volume, temperature, chemical potential, but not number of particles, how can I determine npav and nexc?

I do not have here my copy of the Frenkel and Smit book. However, the usual way of implementing Metropolis Monte Carlo simulations is to have cycles of proposed new configurations equal to the number of particles. Notice that this is just a conventional choice. Nothing prevents having a cycle with a number of elementary steps different from the number of atoms. Then, in GCMC one has to decide how many of these steps have to be displacements and how many add or delete particles. Provided that npav and nexc are not zero, their choice has no direct relationship with the physical conditions of the simulation. Their value impacts the efficiency of the algorithm, exactly like the maximum displacement in the case of displacements moves. Precisely like the maximum displacement in canonical MC, in principle, one could optimize it by looking for the value minimizing the statistical uncertainty on the averages over a given total number of steps. However, again like in the case of the maximum displacement, the exact values are not critical, and the much simpler criterion of obtaining a number of moves of each kind large enough to have reasonable statistics works pretty well.

• Hi GiorgioP, thanks for explaining this! I have attached a screenshot of the algorithm in Frenkel and Smit in my post. I have a couple of questions after reading your explanation and I will post it as an answer below to avoid things getting messy here. (Sorry I need help - I am kind of left alone to do training on molecular dynamics/simulations) Apr 10, 2021 at 16:59
• Hi GiorgioP, my questions are posted. Apr 10, 2021 at 18:15

Previously when I was working on Metropolis Monte Carlo NVT, my algorithm per cycle, for a system of 200 particles, is such that:

for i = 1:200
1) suggest a trial move for particle i
2) apply boundary conditions
3) calculate the change in energy due to this trial move
4) Acceptance criteria
IF change_in_energy < 0
accept trial move
ELSEIF random_number < exp(-change_in_energy/kT)
accept trial move


I think this matches what you mentioned about "implementing Metropolis Monte Carlo simulations is to have cycles of proposed new configurations equal to the number of particles." - I just loop through every particle in my system.

As for GCMC, the values for npav and nexc are determined from trial and error?

My task is to write a GCMC program and perform a simulation for the insertion of methane molecules in an empty cubic box with length equal to $$24\;Å$$, $$T=300\;K$$, pressure range $$10^0 - 10^4\;kPa$$.

I suppose that one way to determine the convergence of GCMC is by looking at the number of particles in the system - it will increase and eventually fluctuate around an average value?