The toy model of statistical entropy that I talked about the other day is the sort of thing that, were I a good computational physicist, I would’ve banged out very quickly. I’m not a good computational physicist, but by cargo-culting my way through some of the VPython examples, I managed to get something that mostly works:
The graph at the bottom of that window is the entropy versus “time” for a lattice of 20 sites with a 25% hopping probability (either left or right). The window with the colored balls at the upper left is a graphical representation– red dots are “occupied” sites, white “unoccupied.” The VPython code to do this is here: entropytestN2.py, from which you can see that I’m not much of a computational physicist, let alone a programmer (there are some really kludgey initialization steps that exist because it kept throwing strange errors when I did simpler things that should’ve worked, and I didn’t have the patience to figure out the underlying problem).
It is, as I expected, kind of fun to play around with. I think there might be some bias built in that ends up pushing things to the right, but I’m not sure whether that’s real or just my inability to recognize randomness. As expected, for small numbers of sites, you see frequent fluctuations down to zero, and as you increase the number of sites, significant decreases in entropy become much rarer.
I haven’t attempted to do anything quantitative with this, because I don’t know how to write data out from VPython into a form that any of my other analysis programs might use. If I think of a good way to do it, I might try to quantify this a little more, but for now, I’m happy just to have a working toy.
And, thinking about it on the way home, I realized that there is, in fact, a subtle rightward bias to the motion of the “atoms” in this toy model, and I know what I need to do to fix it. But like a good faculty member, I’ll leave it as an exercise to the reader to identify what the problem is.
I assume you figured out that this happens because Site[-1] calls the last item in the array (i.e. the rightmost site)so you have atoms on the left jumping to the right.
I glanced at the code and it looks to me that the reason for the biass is that right jumps are resolved before left ones.
So if there are 2 atoms with one vacant space between them the chance that it will be the right atom jumping left is lower since by the time this atom is considered the left atom might have already occupied the vacant space.
It’s not just the order of resolving of right-left jumps which is the problem, but also the fact that you always sweep through the array from left (1) to right (NSites). I kludge-fixed your python by mirroring the bit that does the right-left moves, choosing randomly between right-left or left-right, and then replicating the whole lot again with a reverse sweep through the array and choosing randomly whether to do forward or backward sweeps. It’s ugly, but it works.
A better way would be to respect detailed balance properly by choosing occupied sites at random and moving left or right with equal probability (throwing away those moves which are disallowed). You may also want to use a better quality random number generator than random().
The direction of the sweep seems to be the main issue– when I changed it to go from right to left rather than left to right, it mostly took care of the problem. (It tends to have slightly more atoms in the left half than the right at equilibrium, but nothing like the imbalance I got going left to right.)
Weirdly, alternating from-the-left and from-the-right updates didn’t make that much of a difference. It slowed the process down, but if it ran long enough (~10,000 time steps), it still ended up with nearly everything in the right half. From-the-right only is more equal, though that doesn’t make any sense to me.
Quick and dirty 2D version, with no graph of entropy, but a counter:
http://www.openprocessing.org/visuals/?visualID=50154
Thanks for the fun idea.
The problem i mentioned was related to sweep direction. If we have a single vacant space and 2 atoms on the sides that both want to jump there then the left atom will do it because of the sweep direction and the way jumps are resolved.
But now after checking the code more througoutly I see that there is even bigger bias with sweep direction – an atom can make multiple jumps in one sweep as long as they are in the direction of the sweep.
However since the atoms can pass through the edges of the array (ie atoms can exit on the right side and enter on the left side and vice versa) the bias in direction of the sweep in itself should not lead to atoms concentrating in any particular place. Rather the bias makes atoms more likely to travel together in the direction of the sweep at some speed determined by parameters of the simulation and it’s this speed combined with simulation time which decides where they will end up.
Anyway as far as I can see making the sweeps alternately from left and from right should alleviate the problem by making the biasses cancel each other.
If you replace
for i in range(NSites):
by
if n//2==n/2:
R=range(NSites)
else:
R=range(NSites-1,-1,-1)
for i in R:
that should do the trick.
The best way however would be to update atoms in random order.
I’m not sure that sweep direction is the full explanation of this problem; if it were, then alternating the sweep direction should have fixed the problem (and according to Chad, it didn’t). In any case, the bias is less than posters 3 and 4 indicate. Paul is correct that if two atoms have an empty space between them, then the one that wants to move in the direction of the sweep can move there and prevent the other from moving in there. There is an opposite effect when two atoms occupy neighboring spaces: both of them can move against the sweep (the first one shifts, freeing a space for the second one to shift into), but only the second one can move in the direction of the sweep (it prevents the first one from moving that direction). My suggestion would be to do it as a two-stage process (similar to the algorithm in the game Diplomacy): first come up with a list of proposed moves, then disallow moves that would conflict with other moves (or with atoms not moving). I understand that people who do N-body simulations for a living do something similar: first calculate the forces on each particle (or grid cell) from the current distribution of particles, then move the particles according to those forces.
In order to generate proper canonical averages, it’s not sufficient just to alternate the direction of the sweeps; you need to choose one direction or other at random at each step to ensure true microscopic reversibility. With this measure in place, it doesn’t matter if sites are occasionally blocked in one direction or the other (although it will reduce efficiency of the algorithm).
Chad – when you switched the direction of the sweep to right to left, did you still have probleft=0.25 as in your original code posted above? This would explain why it only (sort of) works that way, because you have two biases in opposite directions that almost cancel each other out.
Regarding getting data out of the program, if your VPython syntax is relatively close to recent versions of regular Python, try something like
######### #Before the loop outfile = open("/home/outfilename.tab", 'w') #Within the while loop outfile.write( 't'.join( [ str(i) for i in statearray ] ) + 'n') #After the loop outfile.close() ######Should get you a tab-separated data-table of the statearray, with one line for each of outfile.write() statements.
The join() method of a string is the Python idiom for concatenating a list of strings, and the “[ … for … in …]” syntax is a list comprehension, which make dealing with lists *much* easier.
As I notice from the screenshot you’re using Windows, I should also say that I believe you can also use Windows style path names to specify the file:
outfile = open(r”C:Documentsoutfile.tab”,’w’)
The r is there before the quotes to make it a “raw” string, so the backslashes aren’t interpreted as escape characters. (e.g. there’s a newline in “C:Documentsnewfile.tab”)