I am trying to reproduce the following in C language:

1. The simulation starts by assuming an infinite reservoir compound of CO and NO molecules with partial pressures 𝑦CO and 𝑦NO(= 1 − 𝑦CO), respectively.

2. The reservoir is in contact with a surface which is simulated by means of a square lattice of 64×64. Periodic boundary conditions are applied to avoid edge effects.

3. The equilibrium coverages are measured as a function of 𝑦CO. To find the critical points, a run each up to 50000 Monte Carlo cycles (MCS) needs to be carried out.If the run completes 50000 MCS without the lattice getting poisoned (surface saturation), the particular point is considered to be within SRS.

4. In order to obtain the coverages of CO, N and O, the initial 20000 MCS are disregarded for equilibrium and averages are taken over the subsequent 30000 MCS. The values of coverages/production rate are obtained after 10 MCS, so that the final values are an average of 3000 configurations.

5. The CO or the NO molecules are selected randomly with relative probabilities 𝑦CO and 𝑦NO, respectively.

6. A surface site is picked up at random and the steps involved in the simulation are as follows:

(a) If the selected molecule is NO and the picked

site is empty, then the four nearest neighbours of the selected site are scanned at random. If all the sites areoccupied, the trial ends. In the case that any adjacent site is vacant, NO dissociates into N and O which occupy the vacant sites with equal probability. The nearest neighbouring sites of adsorbed N and O atoms in step are then searched for the presence of CO and N respectively. If there is a CO (N), then CO2 (N2) is formed, leaving two vacant sites.

(b) In the case that CO is selected, one of the following events occurs: (i) if the picked site is vacant, then CO is adsorbed on the chosen vacant site

as in step. The nearest neighbours of this adsorbed CO molecule are scanned for the presence of an adsorbed O atom in order to complete reaction step. In the case that the O atom is found, CO2 is formed which immediately desorbs leaving two empty sites from the surface and the trial ends. (ii) If the picked site is occupied by the O atom, then CO is co-adsorbed on that O atom with probability 𝑝. This coadsorbed CO then reacts with O at that site with unit probability to form CO2.

**What I have tried:**
This is the code I've written so far:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <assert.h>
#include <stdbool.h>
#include <omp.h>
#include "pcg_basic.h"
#define NMAX 64
#define MMAX 64
#define TIME 30000
enum { EMPTY, CO, O, N};
typedef struct {
double CO_coverage;
} Coverage;
pcg32_random_t rng;
Coverage calculate_coverage(int lattice[NMAX * MMAX]) {
Coverage coverage = {0.0};
for (int i = 0; i < NMAX * MMAX; i++) {
if (lattice[i] == CO) coverage.CO_coverage += 1.0;
}
coverage.CO_coverage /= (NMAX * MMAX);
return coverage;
}
int find_random_site() {
return pcg32_boundedrand_r(&rng, NMAX * MMAX);
}
void print_lattice(int lattice[NMAX * MMAX]){
for (int i = 0; i < NMAX; i++) {
for (int j = 0; j < MMAX; j++) {
printf("%d ", lattice[i * MMAX + j]);
}
printf("\n");
}
}
int find_random_neighbor(int site) {
int x = site / MMAX;
int y = site % MMAX;
int direction = pcg32_boundedrand_r(&rng, 4);
if (direction == 0) {
x = (x - 1 + NMAX) % NMAX;
} else if (direction == 1) {
x = (x + 1) % NMAX;
} else if (direction == 2) {
y = (y - 1 + MMAX) % MMAX;
} else {
y = (y + 1) % MMAX;
}
return x * MMAX + y;
}
void find_nearest_neighbors(int site, int exclude, int neighbors[4]) {
int x = site / MMAX;
int y = site % MMAX;
int x_exclude = exclude / MMAX;
int y_exclude = exclude % MMAX;
neighbors[0] = ((x - 1 + NMAX) % NMAX) * MMAX + y;
neighbors[1] = ((x + 1) % NMAX) * MMAX + y;
neighbors[2] = x * MMAX + (y - 1 + MMAX) % MMAX;
neighbors[3] = x * MMAX + (y + 1) % MMAX;
}
void shuffle(int *array, int n) {
for (int i = n - 1; i > 0; i--) {
int j = pcg32_boundedrand_r(&rng, i + 1);
int temp = array[i];
array[i] = array[j];
array[j] = temp;
}
}
int coadsorption_events = 0;
void add_CO(int lattice[NMAX * MMAX], int *co2_count, double coadsorption_probability) {
int site = find_random_site();
if (lattice[site] == EMPTY) {
lattice[site] = CO;
int neighbors[4];
find_nearest_neighbors(site, -1, neighbors);
int order[4] = {0, 1, 2, 3};
shuffle(order, 4);
for (int i = 0; i < 4; i++) {
int neighbor_site = neighbors[order[i]];
if (lattice[neighbor_site] == O) {
lattice[site] = EMPTY;
lattice[neighbor_site] = EMPTY;
(*co2_count)++;
return;
}
}
}
else if (lattice[site] == O){
double random_number = ldexp(pcg32_random_r(&rng), -32);
if(random_number < coadsorption_probability){
lattice[site] == EMPTY;
(*co2_count)++;
coadsorption_events++;
return;
}
}
}
void add_NO_dissociate(int lattice[NMAX * MMAX], int *co2_count, int *n2_count) {
int site = find_random_site();
int adj_site = find_random_neighbor(site);
if (lattice[site] == EMPTY && lattice[adj_site] == EMPTY) {
lattice[site] = N;
lattice[adj_site] = O;
bool reacted_with_co = false;
int neighbors_3[4];
find_nearest_neighbors(adj_site, -1, neighbors_3);
int order_3[4] = {0, 1, 2, 3};
shuffle(order_3, 4);
for (int j = 0; j < 4; j++) {
int neighbor_site_3 = neighbors_3[order_3[j]];
if (lattice[neighbor_site_3] == CO) {
lattice[adj_site] = EMPTY;
lattice[neighbor_site_3] = EMPTY;
(*co2_count)++;
reacted_with_co = true;
break;
}
}
int neighbors_1[4];
find_nearest_neighbors(site, adj_site, neighbors_1);
int order_1[4] = {0, 1, 2, 3};
shuffle(order_1, 4);
for (int i = 0; i < 4; i++) {
int neighbor_site = neighbors_1[order_1[i]];
if (lattice[neighbor_site] == N) {
lattice[site] = EMPTY;
lattice[neighbor_site] = EMPTY;
(*n2_count)++;
break;
}
}
return;
}
}
double mean(double *array, int size) {
double sum = 0.0;
for (int i = 0; i < size; i++) {
sum += array[i];
}
return sum / size;
}
double std_dev(double *array, int size) {
double avg = mean(array, size);
double sum_sq_diff = 0.0;
for (int i = 0; i < size; i++) {
double diff = array[i] - avg;
sum_sq_diff += diff * diff;
}
return sqrt(sum_sq_diff / size);
}
int main() {
clock_t start = clock();
pcg32_srandom_r(&rng, time(NULL) ^ (intptr_t)&printf, (intptr_t)&rng);
int lattice[64 * 64];
int num_simulations = 50;
int num_runs = 5;
double XCO_values[num_simulations];
double co_coverages_mean[num_simulations];
double co2_rates_mean[num_simulations];
double co_coverages_std_dev[num_simulations];
double co2_rates_std_dev[num_simulations];
double coadsorption_probability = 0.50;
for (int i = 0; i < num_simulations; i++) {
double XCO = i / (double)(num_simulations - 1);
XCO_values[i] = XCO;
double co_coverages[num_runs];
double co2_rates[num_runs];
for (int run = 0; run < num_runs; run++) {
for (int j = 0; j < 64 * 64; j++) {
lattice[j] = EMPTY;
}
int co2_count = 0, n2_count = 0;
Coverage sum_coverage = {0};
double sum_co2_rate = 0, sum_n2_rate = 0;
int data_points = 0;
for (int t = 0; t < 50000; t++) {
for (int k = 0; k < 64 * 64; k++) {
double r = ldexp(pcg32_random_r(&rng), -32);
if (r < XCO) {
double r1 = ldexp(pcg32_random_r(&rng), -32);
add_CO(lattice, &co2_count, coadsorption_probability);
} else {
add_NO_dissociate(lattice, &co2_count, &n2_count);
}
}
if (t >= 20000 && (t - 20000 + 1) % 10 == 0) {
Coverage current_coverage = calculate_coverage(lattice);
sum_coverage.CO_coverage += current_coverage.CO_coverage;
data_points++;
}
}
co_coverages[run] = sum_coverage.CO_coverage / data_points;
co2_rates[run] = sum_co2_rate / data_points;
}
co_coverages_mean[i] = mean(co_coverages, num_runs);
co2_rates_mean[i] = mean(co2_rates, num_runs);
co_coverages_std_dev[i] = std_dev(co_coverages, num_runs);
co2_rates_std_dev[i] = std_dev(co2_rates, num_runs);
}
FILE *file = fopen("results.txt", "w");
if (file == NULL) {
printf("Error opening file!\n");
return EXIT_FAILURE;
}
for (int i = 0; i < num_simulations; i++) {
fprintf(file, "%f\t%f\t%f\t%f\t%f\n", XCO_values[i], co_coverages_mean[i], co2_rates_mean[i], co_coverages_std_dev[i], co2_rates_std_dev[i]);
}
fclose(file);
printf("The results were printed in the file\n");
clock_t end = clock();
double cpu_time_used = ((double) (end - start)) / CLOCKS_PER_SEC;
printf("Total runtime: %f seconds\n", cpu_time_used);
return EXIT_SUCCESS;
}

I know it's huge, but I have included absolutely all the parts of it for more context. In this version, I am considering only the coverage of CO and production of CO2 to simplify everything a bit.

My problem is that no matter what the value of the probability of coadsorption I use, it doesn't affect my results and I get the same every time I run it for every value of this probability. I don't know if it has to do with the way i am collecting the data. I am really confused and don't know what to do next.

Note: I started to use the PCG generator, the basic version from

https://www.pcg-random.org/download.html#id3, because the standard generator wasn't enough.