Monte Carlo Simulation of a liquid - Need help computing g(r)

[JRPB]

EDIT: I'm not in my sharpest moment. I just found a bunch of posts that discuss this. I'll read them and update this post if I find an answer.

I'm working on a hard-sphere MC simulation (for a class). I want to compute the radial distribution function $g(r)$. To put you in context, as to my non-expertise in statistical mechanics, I'm a physics senior undergrad taking the last couple of courses.

My program generates valid states for a non-periodic square box. States are computed with a restricted random walker: particles stay inside the box and do not overlap each other. I want to compute thermodynamic quantities. My ignorance led me to believe (taking equiprobability of micro-states for granted) that sampling such states is equivalent to working with an NVE ensemble. I asked one of my professors about this connection and realized that it's not that simple.

I learned that you can use the RDF to compute some of these quantities. This leads me to my current predicament. How do I compute this function? I tried the literature and couldn't make heads or tails of the definition. As a stab-in-the-dark I tried my own method to see if I could reproduce some of the good looking RDFs around. I came up with an RDFesque function. Far from smooth.

My algorithm:

1) Surround the box with 8 replicas (same particle state, just displaced)
2) Pick a random particle from the original set
3) Compute the maximum circumference (centered at the particle) that fits inside the 9 boxes.
4) Create a circular shell of radii: r, r+dr and count all the particles' centers that lie within this shell.
5) Divide this number by the shell's area
6) Repeat from step 2 and average over N calculations

Attached is an image of my result. The labels are in spanish (it's still pretty clear what is what). The top right chart is a time series of the accept-rate of my random walker. Ignore the x-scale on the g(r). I'm plotting distance from particle v.s. g(r). I tested different shell widths $dr$. $\frac{\sigma}{10}$ seems ok. Anything smaller gets noisy very fast.

Any suggestions or pointers to references on how to correctly implement this will be greatly appreciated.

 

[klawlor419]

I don't know that this will help, but the radial distribution function is just simply a tally of the number of particles surrounding each of the particles in a system, from right next to particle to a large distance off it.

For a liquid, like a Lennard-Jones model fluid or something like it, what is typically done is the create the tally for each of the particles and average the tallies. This gives more smoother RDF.

Your statement will be something like this, if a particle center is within r and r+dr to the particle you are computing the radial distribution function for, tally = tally + 1 . There are some averaging schemes were you can perform the averaging as you go, otherwise you can you just add all the tallies together and divide by N-number of particles.

This is roughly what the RDF is.

 

[JRPB]

I don't know that this will help, but the radial distribution ...
This is roughly what the RDF is.

Thanks for the explanation. I'll post the solution for anyone interested (as soon as I finish it, I have some school work I have to catch up to first)

 

[JRPB]

Update

This is an update on my progress. I found an undergrad thesis work on hard disks, I read it and got a much better looking g(r). It still needs improvement. I'm calculating shells+neighbors for the entire set of disks. Which obviously leaves me with some very unbalanced g(r)'s to average over (those disks near the edges of the box will have half-empty shells -or half-full, aha!-).

I thought about adding clones of the simulation box and surround the orginal one. And THEN survey the one box. However, I'm afraid it'll induce bogus results because of the noticeable structure that border effects create. Any stats. mechs. with useful suggestions?

Attachment: an example run of my simulator.

 

[JRPB]

Alright. What I was trying to do won't work well on a non-periodic box. I got a suggestion (from an experienced Stat. Mech.) to calculate a density profile instead of the g(r). Something along the lines of: count how many particles are in your grid's cells, add up over entire simulation, normalize and fit a polynomial whose "contact value" (at x = 0) is a decent estimate of the pressure in the bulk.

Repeating over different volume fractions one could, in principle, sample points from the equation of state of the "hard disk fluid". That's it for me, I won't be posting more updates (I'm taking a different path).