Bernard Krauth 2012
From Werner KRAUTH
(Difference between revisions)
| Revision as of 22:07, 29 November 2011 Werner (Talk | contribs) (→Illustration) ← Previous diff |
Revision as of 22:08, 29 November 2011 Werner (Talk | contribs) (→Illustration) Next diff → |
||
| Line 16: | Line 16: | ||
| =Illustration= | =Illustration= | ||
| - | [[Image:Event stepped move.jpg|200px|]] Here's an illustration of the algorithm for the [http://en.wikipedia.org/wiki/Tower_of_Hanoi Tower-of-Hanoi potential] | + | [[Image:Event stepped move.jpg|200px|]] Here's an illustration of the algorithm for the [http://en.wikipedia.org/wiki/Tower_of_Hanoi Tower-of-Hanoi potential]. |
| <br clear="all" /> | <br clear="all" /> | ||
Revision as of 22:08, 29 November 2011
E. P. Bernard, W. Krauth Event-driven Monte Carlo algorithm for general potentials (Preprint ArXiv (2011))
Contents |
Paper
Abstract: We extend the event-chain Monte Carlo algorithm from hard-sphere interactions to the micro- canonical ensemble (constant potential energy) for general potentials. This event-drivenMonte Carlo algorithm is non-local, rejection-free, and allows for the breaking of detailed balance. The algorithm uses a discretized potential, but its running speed is asymptotically independent of the discretization. We implement the algorithm for the cut-off linear potential, and discuss its possible implementation directly in the continuum limit.
Electronic version (from arXiv, original version)
Illustration
Here's an illustration of the algorithm for the Tower-of-Hanoi potential.
Algorithm
Here is a naive Python implementation of the algorithm
#!/usr/bin/python
###========+=========+=========+=========+=========+=========+=========+=
## PROGRAM : event_stepped.py
## TYPE : Python script (python 2.7)
## PURPOSE : Event-chain simulation for N particles with a stepped
## potential in a square with periodic boundary conditions.
## COMMENT : L is a list of tuples
## VERSION : 29 Nov 2011
##========+=========+=========+=========+=========+=========+=========+
from random import uniform, randint, choice, seed
import math, pylab, sys, cPickle
Energy_check_iteration = 0
def translate(x,del_x):
x_trans=((x[0]-del_x[0])%box,((x[1]-del_x[1])%box))
return x_trans
def dist(x,y):
"""periodic distance between two two-dimensional points x and y"""
d_x= abs(x[0]-y[0])%box
d_x = min(d_x,box-d_x)
d_y= abs(x[1]-y[1])%box
d_y = min(d_y,box-d_y)
return math.sqrt(d_x**2 + d_y**2)
def potential(r,k):
"""particles have a radius of interaction sigma = 1. The interaction passes
in k steps from 1 (at r=1) to zero (at r=0). Conditions on interaction range
and box size
2 sigma < box
assure that two particles interact only through one path."""
if r > sigma:
pot = 0
else:
pot = k - int(r*k)
return pot
def Energy_calc(L):
"""elementary calculation of the energy for the configuration in list L """
Energy=0
for k in range(1,N):
for l in range(k):
Energy += potential(dist(L[k],L[l]),Nstep)
return Energy
def x_image_calc(x_vec):
"""x_vec describes the position of a particle inside the box. This
function computes the closest "image" of x_vec to the point (0,0). The
conditions on the box size box, the interaction range sigma and the
length of the interval l_max
l_max < box/2.
l_max + 2 sigma < box
assure that the image position is closest to the entire interval [(0,0)
(l_max,0)."""
x = x_vec[0]
y = x_vec[1]
if x < 0 or x > box or y < 0 or y > box: print 'error',sys.exit()
if x < box/2: x_image = x
else : x_image = x - box
if y < box/2: y_image = y
else: y_image = y-box
return (x_image, y_image)
def x_intersection(x_vec,d):
"""compute intersections between the circle of radius d, with center x_vec, and the x-axis.
We suppose that abs(x_vec[1]) < d, so that two such intersections exist."""
del_x= math.sqrt(d**2 - x_vec[1]**2)
return (x_vec[0]-del_x, x_vec[0]+del_x)
def Flip_conf(L,rot):
L_flip=[]
for (a,b) in L:
if rot == 1: L_flip.append((box - b,a))
else: L_flip.append((b,box - a))
return L_flip,-rot
def Energy_check(L,message_string):
global Energy_check_iteration,Energy_level
Energy_check_iteration += 1
a= Energy_calc(L)
if a == Energy_level:
return ('OK',a,Energy_check_iteration)
else:
return ('Error',a,Energy_level,Energy_check_iteration)
#==============================================================================================
#=========================== main program starts here =========================================
#==============================================================================================
#
seed('Test')
Event_eps=1.e-12 # disallow events within too close bounds.
Zero = 0.0
One = 1.0
box = 4.0
sigma = 1.0
N = 4
Nstep = 100
data = []
while True:
L = [(uniform(0.,box),uniform(0.,box)) for k in range(N)]
Energy_level_initial = Energy_calc(L)
if Energy_level_initial == 90: break
Energy_level = Energy_level_initial
rot = 1
ltilde = 0.03
for iter in xrange(100000):
L_save = L[:]
if iter % 10000 ==0: print iter
i=randint(0,N-1)
j=(i+randint(1,N-1))%N
data.append(dist(L[i],L[j]))
if randint(0,1) < 1: L,rot = Flip_conf(L,rot)
#
# this iteration will be used up when the Zero_distance_to_go falls to zero
#
Zero_distance_to_go = ltilde
next_particle = choice(L)
while Zero_distance_to_go > Zero:
Total_distance_to_go = min(box/2,box/2-sigma)*.3971
L.remove(next_particle)
L = [translate(x,next_particle) for x in L]
next_events = []
Current_position = Zero
for x in L:
x_image = x_image_calc(x)
for k in range(1,Nstep+1):
distance = k/float(Nstep)
if abs(x_image[1]) < distance:
x_dummy = x_intersection(x_image,distance)
if x_dummy[0] > Event_eps and x_dummy[0] < Total_distance_to_go:
next_events.append((x_dummy[0],1,x))
if x_dummy[1] > Event_eps and x_dummy[1] < Total_distance_to_go:
next_events.append((x_dummy[1],-1,x))
next_events.append((float("inf"),0,(float("inf"),Zero))) # Final event, at infinity
next_events.sort(reverse=True)
while min(Total_distance_to_go,Zero_distance_to_go) > Zero:
#
# this single-particle loop leads to either:
# 1. A particle collision
# 2. Zero_distance_to_go = 0
# 3. Total_distance_to_go = 0
#
next_position,next_energy_level,next_particle = next_events.pop()
distance_to_next_event = next_position - Current_position
if Energy_level < Energy_level_initial:
if Total_distance_to_go < distance_to_next_event:
Current_position += Total_distance_to_go
next_particle = (Current_position,Zero)
L.append(next_particle)
break
else: # go to event, but energy will not be too high: remain in loop of one particle...
Current_position += distance_to_next_event
Energy_level += next_energy_level
Total_distance_to_go -= distance_to_next_event
else: # Energy_level == Energy_level_intial
min_dist= min(Total_distance_to_go, Zero_distance_to_go)
if min_dist < distance_to_next_event:
Current_position += min_dist
Total_distance_to_go -= min_dist
Zero_distance_to_go -= min_dist
next_particle = (Current_position,Zero)
L.append(next_particle)
break # we can break here, because either Tdtg or Zdtg are zero, but next particle is old one
else:
#
# We go to a new event.
#
Current_position += distance_to_next_event
Total_distance_to_go -= distance_to_next_event
Zero_distance_to_go -= distance_to_next_event
if next_energy_level < 0: # we went to event, but energy is not too high, remain in loop of one particle
Energy_level += next_energy_level
else: # exit one-particle loop: use the next particle as given in the list...
L.append((Current_position,Zero))
break
if Energy_calc(L) != 90:
L = L_save[:]
print 'error at ', iter
f=open("event_stepped_13.data","w")
cPickle.dump(data,f)
f.close()
pylab.hist(data,bins=40,normed=True)
pylab.title("Event_stepped ltilde= "+str(ltilde)+" N= "+str(N)+" Energy = "+str(Energy_level))
pylab.savefig('Event_stepped.png')
pylab.show()
