Click here to Skip to main content
15,792,510 members
Please Sign up or sign in to vote.
0.00/5 (No votes)
See more:
I need to calculate the mean of the coverages of two trapped atoms in a square surface. Such averages are taken after certain number of steps, for example, every 10 monte carlo steps I take the average for such coverages, then I print in a file each average and its respective "y" value (while y is less or equal than 1).

What I have tried:

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>

//define the dimensions of the grid
#define MAX_X 4
#define MAX_Y 4
//define the iterations 
#define ITERATIONS (MAX_Y*MAX_X)

FILE *data;

//define the states of a cell in grid`
typedef enum  { S_EMPTY, P1_OCCUPIED, P2_OCCUPIED, S_NONE } gstate;

// help generate random coordinate of the grid
int gridrnd(int max)
{
    return (rand() % max);
}


// generates random coordinates of the grid
int generate_coords(int* j, int* i ) 
{
    if (!i || !j)
        return 1;

    *i = gridrnd(MAX_X);
    *j = gridrnd(MAX_Y);

   // printf("(%d,%d)\n\n", *j, *i);
    return 0;
}

//function to initialize the grid as empty
void grid_init(gstate grid[MAX_Y][MAX_X])
{
    for (int j = 0; j < MAX_Y; j++) {
        for (int i = 0; i < MAX_X; i++) {
            grid[j][i] = S_EMPTY;
        }
    }
}

//funstion to calculate the mean of the coverages
double calculate_mean(double sum, double count){
    double average=0.0;
    average=sum/count;
    return average;
}



int main(){
    data=fopen("data.txt","w");
    int i = 0, j = 0;
    gstate grid[MAX_Y][MAX_X];
    int particle1 = 0, particle2 = 0; // counters for the number of particle1 and particle2
    int availcells = MAX_X * MAX_Y; //first we initialize with all the cells of the matrix available
    int fullcells = 0;
    int rounds = 0;
    double N = 1.0*sizeof(grid)/sizeof(grid[0][0]); //number of the total sites in the matrix
    float r;
    double cover1 = 0.0, cover2 = 0.0, sumacov = 0.0;
    double average1 = 0.0, average2 = 0.0; //we define the average of the coverages of particle 1 and particle 2 
    double sum1 = 0.0, sum2 = 0.0; //sum of the coverages of both particle 1 and 2. useful to calculate the average of the coverages
    float Y = 0.0;
    double MCSTEPS = 10;

    srand((unsigned)time(0));
    
    // Initialize grid to be S_EMPTY
    grid_init(grid);

for(int iter = 0; iter < ITERATIONS; iter++){

        while (Y <= 1) // repeat till the grid is full
        {
            sum1=0.0;
            sum2=0.0;
                       
            for(int time = 0; time < MCSTEPS; time++ ) {
                    //LOCATE AN ENTRY OF THE MATRIX RANDOMLY
                    generate_coords(&j, &i);

                    //EVALUATE THE CHOOSEN SITE
                    switch (grid[j][i])
                    {
                        case S_EMPTY:
                            //printf("IT'S S_EMPTY, LET'S FILL IT WITH A PARTICLE. FIRST LET'S GENERATE TO DECIDE IFIT WILL BE TRAPPED\n\n");
                            
                            r = rand()/(float)RAND_MAX;
                            
                            if(r <= Y){//The particle #1 is chosen         
                              //printf("r = %lf is less than Y = %lf. We choose the particle #1\n\n", r, Y);
                                grid[j][i] = P1_OCCUPIED;
                                particle1++;
                                availcells--;
                                fullcells++;
                            }
                            else{//The particle #2 is chosen
                              //printf("r = %lf is greater than Y = %lf. We choose the particle #2\n\n", r, Y);
                                grid[j][i] = P2_OCCUPIED;
                                particle2++;
                                availcells--;
                                fullcells++;                
                            }  
                            break;
                            
                        case P1_OCCUPIED:
                            //printf("IT'S OCCUPIED WITH THE PARTICLE #1. PLEASE, GENERATE ANOTHER SITE ON THE SURFACE\n\n");
                            break;

                        case P2_OCCUPIED:
                            //printf("IT'S OCCUPIED WITH THE PARTICLE #2. PLEASE, GENERATE ANOTHER SITE ON THE SURFACE\n\n");
                            break;
                        
                        }  

                    cover1 = (double)particle1/(double)N;
                    cover2 = (double)particle2/(double)N;
                    sumacov = cover1 + cover2;// To verify if the sum of the coverages is effectively equal to 1

                    sum1+=cover1;
                    sum2+=cover2;
            }
            average1=calculate_mean(sum1,(double)MCSTEPS);
            average2=calculate_mean(sum2, (double)MCSTEPS);

            Y = Y + 1/((float)(MCSTEPS*ITERATIONS));

            fprintf(data,"%f   %f   %f\n", average1, average2, Y);
            printf("%f   %f   %f\n", average1, average2, Y);
        
        }

  }   
    printf("The process took %d rounds\n\n", rounds);
    printf("#particle1 = %d\n\n", particle1);//total of particle1 adsorbed
    printf("#particle2 = %d\n\n", particle2);//total of particle2 adsorbed
    printf("#availcells = %d\n\n",availcells);
    printf("#fullcells = %d\n\n",fullcells); 

    fclose(data);

    return 0;
}


The problem is that my code almost never selects the particle #1 and I don't know if this is related to the way I'm generating the random number "r" to select between particle #1 or #2. IN consequence the average of the coverage for the particle1 gives 0.000000.

Only to verify: Is the way I'm calculating the coverages and its averages appropiate for each particle?
Posted
Updated 22-Feb-23 3:48am
v2
Comments
CPallini 22-Feb-23 8:23am    
It depends on what 'Y' stand for.
Auyik 23-Feb-23 8:16am    
Well, I'm trying to reproduce the results of a paper from scratch and there, Y is defined as the partial pressure of particle #1, and it's in the interval [0,1]. It is used to decide which molecule will be selected to be trapped on the surface. So, we draw a random number "r" (also in [0,1]). Then, if it is less than Y, we choose particle #1, otherwise, we select particle #2.

In the mentioned paper the authors simulate the interactions between both particles. The idea is to plot Y versus both coverages and obtain the same results they obtained before. But first I'm trying to "fill" that surface with atoms, before adding interactions between them. However, I have this problem while trying to choose between particles #1 and #2.

I can't think of any other way to incorporate "Y" in the simulation to be able to choose between these two atoms. I don't have a lot of programming experience but this simulation is not supposed to be difficult, and I don't know what to do.

1 solution

You program selects, with uniform probability, the float value r between 0 and 1 (both included).
Then it compares such a number with Y value (what Y stand for?) and assigns an empty cell to particle 1 only if r <= Y. However, since Y value is slowly increasing, starting from 0, there is a very little chance for particle 1 to occupy an empty cell.
 
Share this answer
 
Comments
Auyik 22-Feb-23 9:00am    
Thanks. In my simulation "y" represents the "partial pressure" of an atom or the probability of such atom to be chosen. It's greater than zero and less than 1. Later on I will need to reproduce the dependence of the coverages with this parameter (after adding more interactions between particle #1 and particle #2).

The problem is that it's the best way I found to compare "r" with "y" so far
CPallini 22-Feb-23 9:11am    
What is, then, 'partial pressure' ?
I mean probably it should NOT be used to discriminate between particle 1 and 2.
Auyik 22-Feb-23 9:54am    
Y is defined as the partial pressure of the particle #1, and it's in the interval [0,1]. It is used to decide which molecule will be selected to adsorb. So, we draw a random number "r" (also in [0,1]) so that if it is less than Y, we choose the particle #1. Otherwise, we choose the particle #2.

This is part of the algorithm the paper I'm following. In this paper the authors simulate the interactions between both particles. The idea, is to plot Y versus both coverages, to reproduce the results of a model developed some years ago.
Auyik 23-Feb-23 19:41pm    
I think I have found a way to fix the issue. Thank you very much for your answer.
CPallini 24-Feb-23 0:25am    
Glad you found it. You are welcome.

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