Bernard Krauth Wilson 2009

From Werner KRAUTH

(Difference between revisions)
Jump to: navigation, search
Revision as of 15:43, 29 January 2011
Werner (Talk | contribs)

← Previous diff
Revision as of 07:54, 17 October 2018
Werner (Talk | contribs)
(Paper)
Next diff →
Line 1: Line 1:
-'''E. P. Bernard, W. Krauth, D. B. Wilson ''Event-chain algorithms for hard-sphere systems'' Physical Review E ''' 80''' 056704 (2009)'''+__FORCETOC__
 +E. P. Bernard, W. Krauth, D. B. Wilson ''Event-chain algorithms for hard-sphere systems'' Physical Review E 80 056704 (2009)
-'''Abstract:''' In this paper we present the event-chain algorithms, which are fast+=Paper=
-Markov-chain Monte Carlo methods for hard spheres and related systems.+'''Description:''' In this paper, we first presented the event-chain algorithm, a fast irreversible
 +Markov-chain (initially) for hard spheres and related systems.
In a single move of these rejection-free methods, an arbitrarily long In a single move of these rejection-free methods, an arbitrarily long
chain of particles is displaced, and long-range coherent motion can chain of particles is displaced, and long-range coherent motion can
-be induced. Numerical simulations show that event-chain algorithms+be induced. Numerical simulations initially showed that event-chain algorithms
-clearly outperform the conventional Metropolis method. Irreversible+outperform the conventional Metropolis method. The method violates detailed balance but satisfied a global-balance conditions.
-versions of the algorithms, which violate detailed balance, improve the+ 
-speed of the method even further. We also compare our method with+'''Further information:''' We used this algorithm, which is about 100 times faster than the local Metropolis algorithm, for our [[Bernard Krauth 2011|work on the liquid-hexatic transition in two-dimensional hard disks.]] See [[Bernard_Krauth_b_2011|my 2011 paper, with Etienne Bernard,]] for the generalization of the event-chain algorithm to general, continuous, potentials.
-a recent implementations of the molecular dynamics algorithm.+ 
 +This algorithm can be generalized to arbitrary interaction potentials. A first version was presented in a [[Bernard Krauth 2012|paper with Etienne Bernard]], but the definitive method was worked out in collaboration with
 +[[Michel Kapfer Krauth 2013|Manon Michel and Sebastian Kapfer]], in 2014.
 + 
 +Much further research was dedicated to analyzing event-chain Monte Carlo. In [[Kapfer Krauth 2017a| later works, with Sebastian Kapfer], we showed that, in one dimensional models (which can be analyzed more fully), ECMC has much better scaling with the system size than regular Monte Carlo. This was even [[Lei Krauth 2018| proven mathematically]] in a special case, together with Ze Lei.
Line 15: Line 21:
[http://arxiv.org/pdf/0903.2954v2 Electronic version (from arXiv)] [http://arxiv.org/pdf/0903.2954v2 Electronic version (from arXiv)]
-[[Image:Event chain move2.jpg|left|150px|border|Move of the Event-chain algorithm]] +[http://pre.aps.org/abstract/PRE/v80/i1/e056704 Published version in Physical Review E (subscription required)]
 + 
 +=Illustration=
 + 
 +[[Image:Event chain move2.jpg|left|150px|border|Move of the Event-chain algorithm]] Here is the basic schema of the event-chain algorithm (the added length of the arrows equals a fixed value ''l''. The event-chain algorithm is microreversible and it has no rejections. In our applications, we use a version of the algorithm where particles move
 +either in the +''x'' or the +''y'' direction. This is also done in the below Python implementation.
<br clear="all" /> <br clear="all" />
 +
 +=Python Implementation=
 +Here is a basic python implementation of the straight event-chain algorithm for four disks in the version breaking detailed balance.
 +A randomly disk always moves in the +x direction. At each step, the configuration is translated (and flipped) so that the
 +moves starts at <math>(x,y)=(0.,0.)</math>.
 + ###========+=========+=========+=========+=========+=========+=========+=
 + ## PROGRAM : event_chain.py
 + ## TYPE : Python script (python 2.7)
 + ## PURPOSE : Event-chain simulation for 4 particles in a square box
 + ## with periodic boundary conditions.
 + ## COMMENT : L is a list of tuples, move is always in +x direction
 + ## Flip_conf exchanges x and y (effective move in +y direction)
 + ## VERSION : 10 NOV 2011
 + ##========+=========+=========+=========+=========+=========+=========+
 + from random import uniform, randint, choice, seed
 + import math, pylab, sys, cPickle
 + 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 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
 + #=========================== main program starts here =========================================
 + box = 4.0
 + data = []
 + L = [(1.,1.),(2.,2.),(3.,3.),(3.,1.)] # initial condition
 + rot = 1
 + ltilde = 0.9
 + for iter in xrange(1000000):
 + if iter%10000 == 0: print iter
 + i = randint(0,3)
 + j = (i + randint(1,3))%4
 + data.append(dist(L[i],L[j]))
 + if randint(0,1) < 1: L,rot = Flip_conf(L,rot)
 + Zero_distance_to_go = ltilde
 + next_particle = choice(L)
 + while Zero_distance_to_go > 0.0:
 + #
 + # this iteration will stop when the Zero_distance_to_go falls to zero
 + #
 + L.remove(next_particle)
 + L = [((x[0]-next_particle[0])%box,(x[1]-next_particle[1])%box) for x in L]
 + next_position,next_particle = (float("inf"),(float("inf"),0.0))
 + Current_position = 0.0
 + for x in L:
 + x_image = list(x)
 + if x[0]> box/2: x_image[0] = x[0] - box
 + if x[1]> box/2: x_image[1] = x[1] - box
 + if abs(x_image[1]) < 1.0:
 + x_dummy = x_image[0] - math.sqrt(1.0 - x_image[1]**2)
 + if x_dummy > 0.0 and x_dummy < min(Zero_distance_to_go,next_position):
 + next_position,next_particle = (x_dummy,x)
 + distance_to_next_event = next_position - Current_position
 + if Zero_distance_to_go < distance_to_next_event:
 + Current_position += Zero_distance_to_go
 + L.append((Current_position,0.0))
 + break
 + else:
 + Current_position += distance_to_next_event
 + Zero_distance_to_go -= distance_to_next_event
 + L.append((Current_position,0.0))
 + pylab.hist(data,bins=40,normed=True,alpha=0.5,cumulative=True)
 + pylab.title("Event chain for four hard disks, l = "+str(ltilde))
 + pylab.xlabel("Periodic pair distance $r_{ij}$")
 + pylab.ylabel("Integrated probabilities $\Pi(r_{ij})$")
 + pylab.savefig('Event_chain'+str(ltilde)+'.png')
 + pylab.show()
 +
 +=Comparison of output with direct-sampling algorithm=
 +[[Image:Event_direct_stepped0.9.png|left|150px|border|Comparison of Event-chain output with Direct-sampling MC]] Here is a comparison between the event-chain algorithm and a perfect-sampling approach for four disks of size 1/2 in an <math>L\times L</math> square box of length <math>L=4</math>.
 +<br clear="all" />
 +
 +
 +[[Category:Publication]] [[Category:2009]] [[Category:Algorithm]] [[Category:Hard spheres]]

Revision as of 07:54, 17 October 2018

E. P. Bernard, W. Krauth, D. B. Wilson Event-chain algorithms for hard-sphere systems Physical Review E 80 056704 (2009)

Contents

Paper

Description: In this paper, we first presented the event-chain algorithm, a fast irreversible Markov-chain (initially) for hard spheres and related systems. In a single move of these rejection-free methods, an arbitrarily long chain of particles is displaced, and long-range coherent motion can be induced. Numerical simulations initially showed that event-chain algorithms outperform the conventional Metropolis method. The method violates detailed balance but satisfied a global-balance conditions.

Further information: We used this algorithm, which is about 100 times faster than the local Metropolis algorithm, for our work on the liquid-hexatic transition in two-dimensional hard disks. See my 2011 paper, with Etienne Bernard, for the generalization of the event-chain algorithm to general, continuous, potentials.

This algorithm can be generalized to arbitrary interaction potentials. A first version was presented in a paper with Etienne Bernard, but the definitive method was worked out in collaboration with Manon Michel and Sebastian Kapfer, in 2014.

Much further research was dedicated to analyzing event-chain Monte Carlo. In [[Kapfer Krauth 2017a| later works, with Sebastian Kapfer], we showed that, in one dimensional models (which can be analyzed more fully), ECMC has much better scaling with the system size than regular Monte Carlo. This was even proven mathematically in a special case, together with Ze Lei.


Electronic version (from arXiv)

Published version in Physical Review E (subscription required)

Illustration

Move of the Event-chain algorithm
Here is the basic schema of the event-chain algorithm (the added length of the arrows equals a fixed value l. The event-chain algorithm is microreversible and it has no rejections. In our applications, we use a version of the algorithm where particles move

either in the +x or the +y direction. This is also done in the below Python implementation.

Python Implementation

Here is a basic python implementation of the straight event-chain algorithm for four disks in the version breaking detailed balance. A randomly disk always moves in the +x direction. At each step, the configuration is translated (and flipped) so that the moves starts at <math>(x,y)=(0.,0.)</math>.

###========+=========+=========+=========+=========+=========+=========+=
## PROGRAM : event_chain.py
## TYPE    : Python script (python 2.7)
## PURPOSE : Event-chain simulation for 4 particles in a square box
##           with periodic boundary conditions.
## COMMENT : L is a list of tuples, move is always in +x direction
##           Flip_conf exchanges x and y (effective move in +y direction)
## VERSION : 10 NOV 2011
##========+=========+=========+=========+=========+=========+=========+
from random import uniform, randint, choice, seed
import math, pylab, sys, cPickle
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 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
#=========================== main program starts here =========================================
box = 4.0
data = []
L = [(1.,1.),(2.,2.),(3.,3.),(3.,1.)] # initial condition
rot = 1
ltilde = 0.9
for iter in xrange(1000000):
   if iter%10000 == 0: print iter
   i = randint(0,3)
   j = (i + randint(1,3))%4
   data.append(dist(L[i],L[j]))
   if randint(0,1) < 1: L,rot = Flip_conf(L,rot)
   Zero_distance_to_go = ltilde
   next_particle = choice(L)
   while Zero_distance_to_go > 0.0:
#
#     this iteration will stop when the Zero_distance_to_go falls to zero
#
      L.remove(next_particle)
      L = [((x[0]-next_particle[0])%box,(x[1]-next_particle[1])%box) for x in L]
      next_position,next_particle = (float("inf"),(float("inf"),0.0))
      Current_position = 0.0
      for x in L:
         x_image = list(x)
         if x[0]> box/2: x_image[0] = x[0] - box
         if x[1]> box/2: x_image[1] = x[1] - box
         if abs(x_image[1]) < 1.0:
            x_dummy = x_image[0] - math.sqrt(1.0 - x_image[1]**2)
            if x_dummy > 0.0 and x_dummy < min(Zero_distance_to_go,next_position):
               next_position,next_particle = (x_dummy,x)
      distance_to_next_event = next_position - Current_position
      if Zero_distance_to_go  < distance_to_next_event:
         Current_position += Zero_distance_to_go
         L.append((Current_position,0.0))
         break
      else:
         Current_position += distance_to_next_event
         Zero_distance_to_go -= distance_to_next_event
         L.append((Current_position,0.0))
pylab.hist(data,bins=40,normed=True,alpha=0.5,cumulative=True)
pylab.title("Event chain for four hard disks, l = "+str(ltilde))
pylab.xlabel("Periodic pair distance $r_{ij}$")
pylab.ylabel("Integrated probabilities $\Pi(r_{ij})$")
pylab.savefig('Event_chain'+str(ltilde)+'.png')
pylab.show()

Comparison of output with direct-sampling algorithm

Comparison of Event-chain output with Direct-sampling MC
Here is a comparison between the event-chain algorithm and a perfect-sampling approach for four disks of size 1/2 in an <math>L\times L</math> square box of length <math>L=4</math>.


Personal tools