Click here to Skip to main content
16,017,448 members
Please Sign up or sign in to vote.
0.00/5 (No votes)
See more:
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" // Include the PCG header file
#define NMAX 64
#define MMAX 64
#define TIME 30000

enum { EMPTY, CO, O, N}; 

typedef struct {
    double CO_coverage;
} Coverage;

// Initialize the PCG RNG state
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;
}

// 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(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");
    }
}

// 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);
    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;
}

// 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(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++;
			//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(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;
		
        // 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;
    }
}

// Function to calculate the mean
double mean(double *array, int size) {
    double sum = 0.0;
    for (int i = 0; i < size; i++) {
        sum += array[i];
    }
    return sum / size;
}

// Function to calculate the standard deviation
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];  // Changed to 64x64 lattice
    int num_simulations = 50;
    int num_runs = 5; // Number of runs for each value of XCO
    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);  // 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 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;
				// Data collection phase
                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 < 50000; t++) {  // 20000 warm-up + 30000 data collection
					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) {  // Start collecting data after warm-up
						Coverage current_coverage = calculate_coverage(lattice);
						sum_coverage.CO_coverage += current_coverage.CO_coverage;
						
						data_points++;
					}
				}
                // Calculate averages
                co_coverages[run] = sum_coverage.CO_coverage / data_points;
                co2_rates[run] = sum_co2_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);
            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);
		}	
    // Write the data from the arrays to a text file
    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.
Posted
Updated 4-Jul-24 12:25pm
v2

I could only compile the provided source code with changes. Some of the errors could be the cause of problems found.
C
// Coverage coverage = { 0.0, 0.0, 0.0 };
Coverage coverage = { 0.0 }; 
    
// warning C4553: "==": Result of the expression not used. Was "=" meant?
lattice[site] == EMPTY;

// int num_simulations = 50;
#define num_simulations 50

// int num_runs = 5; // Number of runs for each value of XCO
#define num_runs 5

// missing curly bracket
  }
  return EXIT_SUCCESS;
}

The locations at which the clock marks are saved do not seem to make sense.
In particular, only the calculation time should be measured here. The initialization and writing the results to a file should not be measured.
C
int main()
{
clock_t start = clock();
// ...    
fclose(file);
printf("The results were printed in the file\n");
clock_t end = clock();

I also get additional error messages from pcg_basic.c
C
// error C4146: An unsigned type was assigned a unary minus operator. The result is still unsigned.
uint32_t rot = oldstate >> 59u;
return (xorshifted >> rot) | (xorshifted << ((-rot) & 31));

// just like here:
// uint32_t threshold = -bound % bound;
uint32_t threshold = (0x100000000ull-bound) % bound;
 
Share this answer
 
v3
Comments
yuike 4-Jul-24 18:30pm    
Thank you for pointing that out, I just updated the code. Now it should compile without problems.
merano99 4-Jul-24 20:29pm    
what about "lattice[site] == EMPTY;" - makes no sense
yuike 5-Jul-24 3:06am    
You are right. I fixed that and now the coadsorption it's affecting my results. Thank you!

If I'm having problems with other things regarding my simulation. Should I create a new thread?
merano99 5-Jul-24 3:10am    
Yes, please create a new thread with a suitable heading if you have further problems.
yuike 5-Jul-24 3:30am    
Thank you for your help!

Btw it's weird that C didn't alert me about lattice [site]==EMPTY;
Hi,

you have three statements similar to:
C
double some_variable = ldexp(pcg32_random_r(&rng), -32);
and I'm afraid they don't deliver a random double in [0,1) as you probably hope for; I expect a result in [-0.5, +0.5).

Reason: pcg32_random_r() returns an unsigned integer, whereas ldexp obviously expects a signed mantissa; and C doesn't care much about warning you for implicit signed/unsigned conversions.

One way to fix could be:
C
double some_variable = ldexp(pcg32_random_r(&rng) & 0x7FFFFFFF, -31);
making sure the mantissa is always positive.

Extra remarks:
1. the line starting with double r1 = ... seems superfluous;
2. you may avoid the ldexp() function all together by scaling up the other comparand, so it could look like:
C
unsigned int random_number = pcg32_random_r(&rng) & 0x7FFFFFFF;
if(random_number < coadsorption_probability_scaled){
assuming you precalculate once:
C
unsigned int coadsorption_probability_scaled = (unsigned int)(coadsorption_probability * pow(2, 31));
and iff coadsorption_probability is always less than 1, you could use 32 bits rather than 31.

:)

[Added]

Remark 2 is an optimisation; it basically says:
rather than comparing random / pow with testvalue
(where pow = 2^32 or 2^31)
you can equally well compare random with testvalue * pow.
That saves the calling of ldexp() and avoids the doubles.

And since you now stated coadsorption_probability must be allowed to equal 1 you can't use 2^32 (that would overflow), and that is why I showed code based on 2^31.
 
Share this answer
 
v7
Comments
yuike 4-Jul-24 11:56am    
Thank you. Yes, the probability of coadsorption should be between 0 and 1 (both inclusive). The random numbers generated should also be in this interval so that i can decide if CO or NO are adsorbed, or if the coadsortion event will occur.

Note that the probability of coadsorption is a number that I introduce by hand, In my code, I'm considering coadsorption_probability = 0.500. Do I still need to do this calculation: "unsigned int coadsorption_probability_scaled = (unsigned int)(coadsorption_probability * pow(2, 32));"?

Would you mind to elaborate a bit more about why do i need to "scale" the value of the probability of coadsorption?
Luc Pattyn 4-Jul-24 13:20pm    
The first half of my "solution 1" should fix the problems you are having. Did you try it?

Remark 2 is just an optimisation; I expanded on it inside my original message.

:)
yuike 4-Jul-24 16:54pm    
Thank you. I used "double some_variable = ldexp(pcg32_random_r(&rng) & 0x7FFFFFFF, -31);" to generate the random numbers to attempt the adsorption of CO or NO and also, for the coadsorption attemp. However, it didn't work :(

I have other functions that use the generation of random numbers as well. They are the functions find_random_site(), find_random_neighbor(int site) and shuffle(int *array, int n). But I don't know if the random numbers inside of them are being generated correctly? Is the generation of random number well implemented on them?
Luc Pattyn 4-Jul-24 17:08pm    
"it didn't work" is not informative.

did you check my statement "ldexp(pcg32_random_r(&rng), -32); will result in [-0.5, +0.5)"?
did you check my statement "ldexp(pcg32_random_r(&rng) & 0x7FFFFFFF, -31); will result in [0, 1)"?
You wrote "I started to use the PCG generator ... because the standard generator wasn't enough." So I assume you know how to check and evaluate a random generator.
yuike 4-Jul-24 18:18pm    
Oh, sorry. I just run a program to test the random values that both approaches give and they are generating random numbers between 0 and 1.

I said that the standard generator wasn't enough since I was recommended beforehand not to use it because I would need a better quality of pseudorandom numbers for a simulation like this. I implemented the generator of standard library just to check, and couldn't reproduce the "case base" of the model (when the probability of coadsorption is zero). Only when I changed to the PCG, I got the expected results.

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



CodeProject, 20 Bay Street, 11th Floor Toronto, Ontario, Canada M5J 2N8 +1 (416) 849-8900