Help coderanch get a
new server
by contributing to the fundraiser
programming forums Java Mobile Certification Databases Caching Books Engineering Micro Controllers OS Languages Paradigms IDEs Build Tools Frameworks Application Servers Open Source This Site Careers Other Pie Elite all forums
this forum made possible by our volunteer staff, including ...
Marshals:
• Campbell Ritchie
• Jeanne Boyarsky
• Ron McLeod
• Paul Clapham
• Liutauras Vilda
Sheriffs:
• paul wheaton
• Rob Spoor
• Devaka Cooray
Saloon Keepers:
• Stephan van Hulst
• Tim Holloway
• Carey Brown
• Frits Walraven
• Tim Moores
Bartenders:
• Mikalai Zaikin

What's the best way to take averages in a Monte Carlo simulation in C?

Greenhorn
Posts: 10
• Number of slices to send:
Optional 'thank-you' note:
Hi, i'm trying to perform a Monte Carlo simulation that follows these steps:

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.

This is the code I  have so far:

The problem is that for the values of yco until 0.20 i am getting peaks in the coverages and big error bars (I can see them because i am plotting coverages vs yco) for values of the probability of coadsorption different than zero. For such a probability equal to zero, i get the "base case" and everything looks as it should. I was wondering if my code follows exactly the steps i mentioned before  (which come from a paper), or if is there a better way to take averages inside my main function...
coverages-disociative.png

 Message for you sir! I think it is a tiny ad: We need your help - Coderanch server fundraiser https://coderanch.com/t/782867/Coderanch-server-fundraiser