Can I improve the results obtained from a monte carlo simulation in C?

Question

0.00/5 (No votes)

See more:

These are the steps of the simulation I am performing:

1. We assume 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 have so far:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <assert.h>
#include <stdbool.h>
#include "pcg_basic.h" // Include the PCG header file

#define NMAX 64
#define MMAX 64
#define TIME 50000
#define WARMUP 20000
#define num_simulations 50
#define num_runs 5 // Number of runs for each value of XCO

typedef enum { EMPTY, CO, O, N } LAT_ELEM;

typedef struct {
    double CO_coverage;
    double O_coverage;
    double N_coverage;
} Coverage;

// Initialize the PCG RNG state
pcg32_random_t rng;


///////////////////Function Prototypes ////////////////////////////
Coverage calculate_coverage(LAT_ELEM lattice[NMAX * MMAX]);
int find_random_site();
void print_lattice(LAT_ELEM lattice[NMAX * MMAX]);
int find_random_neighbor(int site);
void find_nearest_neighbors(int site, int exclude, int neighbors[4]);
void shuffle(int *array, int n);
void add_CO(LAT_ELEM lattice[NMAX * MMAX], int *co2_count, double coadsorption_probability);
void add_NO_dissociate(LAT_ELEM lattice[NMAX * MMAX], int *co2_count, int *n2_count);
void write_lattice_to_csv(int N, int M, LAT_ELEM lattice[N * M], double XCO);
double mean(double *array, int size);
double std_dev(double *array, int size);
//////////////////////////////////////////////////////////////////

//Calculate coverages function
Coverage calculate_coverage(LAT_ELEM lattice[NMAX * MMAX]) {
    Coverage coverage = {0.0, 0.0, 0.0};
    for (int i = 0; i < NMAX * MMAX; i++) {
        if (lattice[i] == CO) coverage.CO_coverage += 1.0;
        else if (lattice[i] == O) coverage.O_coverage += 1.0;
        else if (lattice[i] == N) coverage.N_coverage += 1.0;
    }
    coverage.CO_coverage /= (NMAX * MMAX);
    coverage.O_coverage /= (NMAX * MMAX);
    coverage.N_coverage /= (NMAX * MMAX);
    return coverage;
}

// Function to find a random site on a grid of size NMAX x MMAX
int find_random_site() {
    return pcg32_boundedrand_r(&rng, NMAX * MMAX);
}

// Print out the lattice
void print_lattice(LAT_ELEM 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");
    }
}

// Function to find a random neighbor of a given site on a grid of size NMAX x MMAX
int find_random_neighbor(int site) {
    int x = site / MMAX;
    int y = site % MMAX;
    int direction = pcg32_boundedrand_r(&rng, 4);
	
	switch (direction) {
	case 0:
		x = (x - 1 + NMAX) % NMAX;
		break;
	case 1:
		x = (x + 1) % NMAX;
		break;
	case 2:
		y = (y - 1 + MMAX) % MMAX;
		break;
	default:
		y = (y + 1) % MMAX;
	}
	return x * MMAX + y;
}

// Function to find the nearest neighbors of a given site, excluding a specific site, on a grid of size NMAX x MMAX
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;
}

// Function to shuffle an array of integers of size n
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;

// Function to add a CO molecule to the lattice
void add_CO(LAT_ELEM 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) & 0x7FFFFFFF, -31);
		if(random_number < random_number < coadsorption_probability){
			lattice[site] = EMPTY;
			(*co2_count)++;
			coadsorption_events++;
			//printf("%d\n",coadsorption_events);
             return;
		}
	}
}

// Function to add an NO molecule to the lattice when it dissociates into N and O
void add_NO_dissociate(LAT_ELEM 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;
			
		
        // Check nearest neighbors of O
        bool reacted_with_co = false; // Flag to check if O reacted with CO
        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; // Set flag to true
                break;
            }
        }
		
        // Check nearest neighbors of N
        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;
    }
}

void write_lattice_to_csv(int N, int M, LAT_ELEM lattice[N * M], double XCO) {
    char filename[50];
    sprintf(filename, "lattice_XCO_%0.2f.csv", XCO);

    FILE *file = fopen(filename, "w");
    if (file == NULL) {
        printf("Error opening file!\n");
        return;
    }

    for (int i = 0; i < N; i++) {
        for (int j = 0; j < M; j++) {
            fprintf(file, "%d", lattice[i * M + j]);
            if (j < M - 1) {
                fprintf(file, ",");
            }
        }
        fprintf(file, "\n");
    }

    fclose(file);
}

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();
	
	// Seed the PCG RNG
	pcg32_srandom_r(&rng, time(NULL) ^ (intptr_t)&printf, (intptr_t)&rng);

    LAT_ELEM lattice[NMAX * MMAX]; 

    double XCO_values[num_simulations];
    double co_coverages_mean[num_simulations];
    double o_coverages_mean[num_simulations];
    double n_coverages_mean[num_simulations];
    double co2_rates_mean[num_simulations];
    double n2_rates_mean[num_simulations];
    double co_coverages_std_dev[num_simulations];
    double o_coverages_std_dev[num_simulations];
    double n_coverages_std_dev[num_simulations];
    double co2_rates_std_dev[num_simulations];
    double n2_rates_std_dev[num_simulations];
	
    double coadsorption_probability = 0.005000;
        
    for (int i = 0; i < num_simulations; i++) {
		double XCO = i / (double)(num_simulations - 1);  // This ensures XCO ranges from 0 to 1
        XCO_values[i] = XCO;

        // Arrays to store the results of multiple runs for the same value of XCO
        double co_coverages[num_runs];
        double o_coverages[num_runs];
        double n_coverages[num_runs];
        double co2_rates[num_runs];
        double n2_rates[num_runs];

        for (int run = 0; run < num_runs; run++) {
			//Initialize the lattice as empty 
            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;
				
			// Combined warm-up and data collection phase
			for (int t = 0; t < TIME; t++) {  // 20000 warm-up + 30000 data collection
				for (int k = 0; k < NMAX * MMAX; k++) {
					double r = ldexp(pcg32_random_r(&rng) & 0x7FFFFFFF, -31);
					if (r < XCO) {
						add_CO(lattice, &co2_count, coadsorption_probability);
					} else {
						add_NO_dissociate(lattice, &co2_count, &n2_count);                
					}		
				}

				if (t >= WARMUP && (t - WARMUP + 1) % 10 == 0) {  // Start collecting data after warm-up
					Coverage current_coverage = calculate_coverage(lattice);
					sum_coverage.CO_coverage += current_coverage.CO_coverage;
					sum_coverage.O_coverage += current_coverage.O_coverage;
					sum_coverage.N_coverage += current_coverage.N_coverage;
					sum_co2_rate += co2_count;
					sum_n2_rate += n2_count;
					data_points++;
				}
			}
			// Calculate averages
            co_coverages[run] = sum_coverage.CO_coverage / data_points;
            o_coverages[run] = sum_coverage.O_coverage / data_points;
            n_coverages[run] = sum_coverage.N_coverage / data_points;
            co2_rates[run] = sum_co2_rate / data_points;
            n2_rates[run] = sum_n2_rate / data_points;
        }

        // Calculate the mean and standard deviation of the coverages and rates of production for this value of XCO
        co_coverages_mean[i] = mean(co_coverages, num_runs);
        o_coverages_mean[i] = mean(o_coverages, num_runs);
        n_coverages_mean[i] = mean(n_coverages, num_runs);
        co2_rates_mean[i] = mean(co2_rates, num_runs);
        n2_rates_mean[i] = mean(n2_rates, num_runs);

        co_coverages_std_dev[i] = std_dev(co_coverages, num_runs);
        o_coverages_std_dev[i] = std_dev(o_coverages, num_runs);
        n_coverages_std_dev[i] = std_dev(n_coverages, num_runs);
        co2_rates_std_dev[i] = std_dev(co2_rates, num_runs);
        n2_rates_std_dev[i] = std_dev(n2_rates, num_runs);
    }
  
    // Write the data from the arrays to a text file
    FILE *file = fopen("3-coadsorption-disasociative-PCO-0-500.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\t%f\t%f\t%f\t%f\t%f\t%f\n", XCO_values[i], co_coverages_mean[i], o_coverages_mean[i], n_coverages_mean[i], co2_rates_mean[i], n2_rates_mean[i], co_coverages_std_dev[i], o_coverages_std_dev[i], n_coverages_std_dev[i], co2_rates_std_dev[i], n2_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;
}

The standard deviations of the coverage of CO, for low yco values (between 0 and 0.22) are very large and I shouldn't have that. In fact, for a yco value around 0.25, the CO coverage should "jump" from about 0.1 to 0.9 and I don't have that. I don't know if I'm implementing the data collection as indicated in the simulation steps I shared above ("the values of coverages/production rate are obtained after 10 MCS, so that the final values are an average of 3000 configurations").

I there a better way to collect my data and get better results?

Posted 5-Jul-24 8:48am

yuike

Updated 6-Jul-24 11:21am

v3

Add a Solution

2 solutions

Add a Solution

This content, along with any associated source code and files, is licensed under The Code Project Open License (CPOL)

merano99 · Answer 1 · 2024-07-06T09:09:00

I was not interested in whether the simulation calculates useful things or whether algorithms can be improved. Here I am going to address source code problems that are immediately obvious and sometimes prevent the source code from being compiled at all on some compilers.

1. constants in C should all be at the top, where you can easily find them.
2. the enum is provided, but is not used and does not even have a name with which it could be used.
3. there should be prototypes for all functions to prevent misunderstandings and assumptions by the compiler. In particular, the data type for lattice should be used consistently.

C

#define NMAX 64
#define MMAX 64
#define TIME 30000

#define  num_simulations  50
#define  num_runs  5          /* Number of runs for each value of XCO */

typedef struct {
    double CO_coverage;
} Coverage;

// enum { EMPTY, CO, O, N };
typedef enum { EMPTY, CO, O, N } LAT_ELEM;

// Prototypes of the functions
Coverage calculate_coverage(LAT_ELEM lattice[NMAX * MMAX]);
void print_lattice(LAT_ELEM lattice[NMAX * MMAX]);
void add_CO(LAT_ELEM lattice[NMAX * MMAX], int* co2_count, double coadsorption_probability);
void add_NO_dissociate(LAT_ELEM lattice[NMAX * MMAX], int* co2_count, int* n2_count);
int find_random_site();
void find_nearest_neighbors(int site, int exclude, int neighbors[4]);
void shuffle(int* array, int n);
double mean(double* array, int size);
double std_dev(double* array, int size);

Instead of multiple jumps with an if statement, a switch-case would be more effective.

C

// if (direction == 0) { ...
switch (direction) {
case 0:
    x = (x - 1 + NMAX) % NMAX;
    break;
case 1:
    x = (x + 1) % NMAX;
    break;
case 2:
    y = (y - 1 + MMAX) % MMAX;
    break;
default:
    y = (y + 1) % MMAX;

The time measurement should only measure the pure computing time and, if possible, no outputs.

The header for omp is available, but omp is not used. The complete simulation is currently only calculated on one core, which of course does not make optimum use of the computing time. In addition to the header, the appropriate flag must also be set depending on the compiler:

C

int max_threads = omp_get_max_threads();
int* thread_printed = calloc(max_threads, sizeof(int));

printf("using %d threads\n", max_threads);

clock_t start = clock();

int i;
#pragma omp parallel for
for ( i = 0; i < num_simulations; i++) {
    int thread_id = omp_get_thread_num();
    if (!thread_printed[thread_id]) {
        printf("Thread %d is processing\n", thread_id);
        thread_printed[thread_id] = 1;
    }
    // ...

Rick York · Answer 2 · 2024-07-06T10:19:00

I agree with the previous poster's recommendations although I prefer to order the functions so that no prototypes are necessary. I will also add a few others. There are still too many constants for my tastes. For example - the number 4 seems to be a magic value. I think that should be a constant or a definition. Also, the expression NMAX * MMAX appears many times but I have a better solution to that issue which will follow. One other thing - it makes little sense to have a structure with a single data member. You can make it a type definition if you think it needs to be but it shouldn't be a structure.

You go through a lot of "gyrations" to make the one-dimensional lattice array act like a two-dimensional matrix. By gyrations I mean look at what happens in find_nearest_neighbors and find_random_site. Instead you could make it an actual matrix. Here is how I would do that :

C

#define LAT_ROWS 64
#define LAT_COLS 64

typedef LAT_ELEM  LatRow[ LAT_COLS ];   // a single row of LAT_COLS columns
typedef LatRow    Lattice[ LAT_ROWS ];  // the lattice, an array of rows.

Using this structure, you can access members of it like this :

C

LAT_ELEM element = Lattice[ row ][ col ];
Lattic[ row ][ col ] = new_element;

Passing it to functions is easy. The key thing to remember is the number of rows does not need to defined to the function being called but it does need to know the size of each row. This means function prototypes and calling the function look like this :

C

void print_lattice( LatRow lattice[] )
{
    for( int row = 0; row < LAT_ROWS; ++row )
    {
        for( int col = 0; col < LAT_COLS; ++col )
        {
            printf( "%d ", lattice[ row ][ col ] );
        }
        printf("\n");
    }
}

// a declration and a function call

Lattice lattice;
print_lattice( lattice );

You could declare the function as taking a Lattice but I wanted to make the point of how one uses functions to accept arrays as arguments.

Lastly, I would consider using C++. I find it to be much easier to deal with because it has fewer constraints. You could still write your code exactly as you are but you have more capabilities and options. For example, instead of defining macros for values like NMAX and MMAX, you could make them constants :

C++

const int LAT_ROWS = 64;
const int LAT_COLS = 64;

There are many advantages to using constants over macros and I will not list them here.

Can I improve the results obtained from a monte carlo simulation in C?

2 solutions

Solution 1

Solution 2

Add your solution here

Preview 0

Can I improve the results obtained from a monte carlo simulation in C?

2 solutions

Solution 1

Solution 2

Add your solution here

Preview 0

Existing Members

...or Join us