# Tutorial - Monte Carlo simulation additive binary gas mixture

Here we provide an example of applying the Monte Carlo method to simulate the adsorption of additive binary gas mixtures on solid surfaces [1]. The additive gas mixture of two components means that lateral interactions are considered only for the adsorbed molecules of the same type. For example, the binary system of A and B components is described by thwo types of lateral interactions: A-A and B-B, but not A-B.

## Lattice Gas Model

Please download the model `Binary_gas_mixture.zip`

.
As a model of the surface we use the square lattice with `a`

distance between nearest neighbor sites.
See the `sample_empty_unit_cell.xyz`

file. There are three state of the adsorption site in this model.
They are set up in the following xyz-files: `adsorbate_state_0.xyz`

, `adsorbate_state_1.xyz`

and `adsorbate_state_2.xyz`

for an empty site, for a site occupied with A molecule and molecule B. All states in the model
are represented by monomeric spherical particles. We assign a mark to all the possible surface states
(molecule A - atom `N`

, molecule B - `O`

, empty - `H`

). The name of the state in each xyz-file can be
found in the head of the corresponding file. The above mentioned lateral interactions between A-A, A-B and B-B molecules
are specified in the `energy`

file and equal to 4.0, 0.0 and 6.0 kJ/mol, respectively.

Please also download the following script `mc_binary_gas.py`

and read below what the script does.

```
lt = load_lattice_task('./')
lt.set_property('coverage', {'A':1.0, 'B':1.0} )
lt.set_property('A_cnt', {'A':1})
lt.set_property('B_cnt', {'B':1})
L = 36 # lattice size
T = 100. # K
muB = 0.
lt.set_ads_energy('B', muB) # kJ/mol
kB = 0.008314 # kJ/(mol*K)
```

Function `load_lattice_task`

load the lattice model into the SuSMoST.
Here we also define some properties of the site states such as coverage
that equal to 1 for both A and B adsorption complexes, and zero (by default) for
an empty site. `A_cnt`

and `B_cnt`

are for the amount of A and B adsorption complexes on the lattice.

Also we are setting up the key constants and model parameters:

`L`

is the linear size of the lattice in`a`

units;`T`

is the temperatures in the adsorption layer;`muB`

is the chemical potential of the B component that is set by`lt.set_ads_energy`

equal to the state energies of the B adsorption complexes.

as

`lt.set_ads_energy`

we set the state energies of the available adsorption complexes as`lt.set_ads_energy`

. Note that in this tutorial, the chemical potential of the components`muA`

and`muB`

means the differences between molar Gibbs free energy of the component in the gas phase and adsorption layer \(\mu_{Bgas} - \mu_{Bads} \approx-RTlnp_B\). The chemical potential of the A component will be varied in the simulation, so it will be set up further.

`kB`

is Gas constant in kJ/(mol*K). If`k_B`

is measured in kJ/K, then`k_B`

is the Boltzmann constant .

## Set the Monte Carlo Simulation

Here we perform the Monte Carlo simulation of the model in Grand Canonical Ensemble for a given amount of Monte Carlo steps at fixed values of the B component chemical potential and temperature, changing the chemical potential of the A molecules.

```
# Perform the simulation for a series of TEtB chemical potentials
for mu in np.arange(30., 0.-0.0001, -2.0 ):
lt.set_ads_energy('A', mu) # kJ/mol
m = mc.make_metropolis(lt, L, temperatures, kB)
mc.run(m, log_periods_cnt=10, log_period_steps = 1000*m.cells_count, relaxation_steps = 100000*m.cells_count, \
traj_fns = ["T={}.xyz".format(T) for T in temperatures])
```

Here we run the `m`

simulation for each value of the chemical potential of the A component in the range from 20 to -30 kJ/mol with step of -2 kJ/mol.
Each time we set new value of the A component chemical potential using `lt.set_ads_energy`

.
`make_metropolis`

function creates a task `m`

for Monte Carlo simulation. Here we perform 100 000 Monte Carlo steps (MCS) for the system relaxation
and 50 000 MCS for thermal averaging over the 50 points.
Recall that MCS is the amount of attempts to switch the state of the lattice equaled to the amount of the lattice sites.
Every 1000 MCS in the production run we save the xyz-snapshots having 50 snapshots at the end of the simulation at each value of the `muA`

.
It takes about 25 min on an average CPU.

## Calculation the thermodynamic averages

We have obtained as simple averages the total (`H`

) and potential energy (`U`

) of the adlayer per lattice site,
the total surface coverage (`coverage`

), partial coverage of the surface with A (`A_cnt`

) and B (`B_cnt`

) molecules
at the temperature 100 K changing the chemical potential of the component A from 20 to -30 with step = -2 kJ/mol.
Also we have calculated the differential heat of adsorption (`Qd`

) and heat capacity (`C`

).
All averages are calculated using 50 points. Results are saved in file.

```
E_2 = 0
covE = 0
cov_2 = 0
for T, param_log in zip(temperatures, m.full_params_log):
stat = np.mean(param_log,axis=0)[[0, -2, -3]] # coverage, H, U
for log_row in param_log:
E_2 += log_row[-3]**2.0
covE += log_row[0]*log_row[-3]
cov_2 += log_row[0]**2.0
E_2 /= len(param_log)
covE /= len(param_log)
cov_2 /= len(param_log)
Qd = -(covE - stat[0]*stat[2])/(cov_2-stat[0]**2.0)
Cp = (L^2)*(E_2 - stat[2]*stat[2])/(T*T*kB)
```

where

`E_2`

is the squared total energy per lattice site;

`covE`

is the value of surface coverage multiplied by potential energy;

`cov_2`

is the squared coverage;

`Qd`

is the differential heat of adsorption calculated with equation [3]\[Q_{d} = \frac{\langle\theta \cdot E\rangle-\langle \theta \rangle \cdot \langle E \rangle}{\langle \theta^2 \rangle-{\langle \theta \rangle}^2}\]

`C`

is the heat capacity of the molecular layer calculated as [4]\[C = L^2 \cdot \frac{\langle H^2 \rangle-\langle H \rangle^2}{T^2 R}\]

## Results

In the Fig.1 the fractional surface coverages with molecules A and B, and the total surface coverage are shown as functions of the chemical potential of the molecules A in the adlayer. It is worth to note that the equivalence of the chemical potentials in gas phase and in the adlayer is supposed. Therefore, the lower values of the chemical potential of the molecule A in the adlayer, the higher pressure of the gas A in the gas phase mixture.

Animation on the Fig.2 illustrates the loading of the adlayer with increase of the chemical potential of A component.

As it is well-known heat capacity of the system has a peak in vicinity of the phase transition. Dependence of the heat capacity of the adlayer on the chemical potential of the molecules A is shown in the Fig. 3. We can see there are few pronounced peaks on the heat capacity curve. They correspond to the continuous phase transitions those related to an appearance or disappearance of a symmetry element [5].

Detailed discussion of this and other results on adsorption of the binary gas mixture can be found in the papers [1] and [2].