[PROGRAM: Fortran90/C/C++] METROPOLIS/MONTE CARLO simulations applied to a simple (substrate/adsorbate) model (p/s)

Considering a particle "p" chemisorbed or physisorbed at a surface "s", we can simulate the ordering pattern of a group of molecules "p" self-organizing at "s" for different coverage and temperature using the Metropolis algorithm in a Monte Carlo cycle in the canonical ensemble.

What is the canonical ensemble?

N,V and T constant: T (Temperature), N (Number of particles), V (Volume)

What is the Metropolis Algorithm?

Algorithm to generate a new configuration with probabilities.
- Eo: Energy of the system before the particle displacement
- En: New energy of the system if the particle has moved
- Kb: Boltzmann constant (eV/K)

Metropolis algorithm

Monte Carlo cycle in the canonical ensemble:

- First step is to chose the parameters of the simulation, (T, Size of supercell, coverage, number of Monte Carlo steps,interaction force between particles...).
- Second step: is the construction of the initial configuration of particles on the surface (Start random, or organized)
- Third step is the Monte Carlo simulation cycle, accepting or rejecting particle moves with respect to the Metropolis algorithm. This calculation stops at when the number of Monte Carlo steps, entered by the user in the intial conditions, is reached... (More M.C. steps = Better precision)

Diagram of a Monte Carlo cycle

How to know if we do it right?

The accuracy of the Monte Carlo program can be verified using different parameters from mathematics and themodynamics, by simulating an order-disorder phase transitions from the lowest temperature to the highest one. The parameters are the heat capacity (Cv) and the order parameter (?). Between the fully ordered state and the disorder, we find a critical phase which is determined by its single point: the critical temperature T.

Order-disorder phase transition (50% coverage)

Simulations for the simple model (p/s)

The particle "p" has only one force of interaction with its first nearest neighbor located on the next nearest hollow site. The interaction energy between 2 particles is defined as +0.1eV. The cell is a 10*10, so for 50% coverage, 50 particles are interacting with each other onto the surface "s"during the simulation.

- Conditions:

Coverage / 50%
interaction particle-particle / 0.1eV
Cell size / 10*10
Start / Random configuration
Monte Carlo Steps / 10 Million

- Order-disorder transition by Statistical Mechanics:

By increasing the temperature we can trace the Cv= f(T) curve, Starting with the perfect ordered state at low temperature, we reach the critical point (for this model) at 1350K, which is a maximum in terms of heat capacity, and then slowly get the full disorder at the highest temperatures:

Heat capacity of the (p/s) system: Cv = f(T)

- Order-disorder transition by discrete Fourier Transform:

The same evolution in terms of order parameter, gives another type of curve, if we consider that 1 is the fully ordered state and 0 the disorder, we obtain the following curve:

Order parameter of the (p/s) system: ?= f(T)

Note that the critical temperature must be the around same value for the heat capacity and the order parameter, to be sure to do it right!

Download the .pdf version of the Monte Carlo program

The link to this file is the following: MCprogram
It is written in Fortran 90 but at this stage it can be easily converted to C, C++ or any other languages...
This simple program can also be easily enhanced, I chose to write it in the simplest way for anyone to understand the different steps from the first line of code to the complete Monte Carlo program. Then the finality of this small presentation is for the user to be able to adapt this program to a more complicated case.

For more questions on Statistical mechanics and Order parameters feel free to drop me a line at: contacts