Kapfer Krauth 2016

From Werner KRAUTH

(Difference between revisions)
Jump to: navigation, search
Revision as of 22:02, 25 September 2016
Werner (Talk | contribs)
(→Paper)
← Previous diff		 Revision as of 22:08, 25 September 2016
Werner (Talk | contribs)
(→Python implementation)
Next diff →
Line 17: Line 17:
The paper contains a Python demo implementation of the cell-veto algorithm (that is also available on [https://github.com/Cell-veto github]). The paper contains a Python demo implementation of the cell-veto algorithm (that is also available on [https://github.com/Cell-veto github]).
- 
- import math, random, sys, pylab 
- from copy import deepcopy 
- def CellRate(TargetCellLL, ActiveCellLL, L, CellBoundary): # always >= 0+ import math, random, sys
- Rate = 0.0+ import numpy as np
- for a in CellBoundary:+
- for b in CellBoundary:+
- del_x = TargetCellLL[0] - ActiveCellLL[0] + (a[0] - b[0]) / L+
- del_y = TargetCellLL[1] - ActiveCellLL[1] + (a[1] - b[1]) / L+
- Rate = max(Rate, PairRate(del_x, del_y, k_max))+
- return Rate+
- def PairRate(del_x, del_y, k_max): # Here for 1/r potential in 2d+ def norm (x, y):
- q = 0.0+ """norm of a two-dimensional vector"""
- for ky in range(-k_max, k_max + 1):+ return (x*x + y*y) ** 0.5
- for kx in range(-k_max, k_max + 1):+
- q += (del_x + kx ) / (+
- (del_x + kx) ** 2 + (del_y + ky) ** 2) ** (3.0 / 2.0)+
- q += 1.0 / ((del_x + kx + 0.5) ** 2 + (del_y + ky) ** 2) ** (1.0 / 2.0)+
- q -= 1.0 / ((del_x - kx - 0.5) ** 2 + (del_y + ky) ** 2) ** (1.0 / 2.0)+
- return max(0.0, q)+
- def CellTranslate(TargetCell, ActiveCell): # TargetCell -> 0 as CellTranslate -> ActiveCell+ def dist (a, b):
- kt_y = (TargetCell // L) % L+ """periodic distance between two two-dimensional points a and b"""
- kt_x = (TargetCell - L * kt_y) % L+ delta_x = (a[0] - b[0] + 2.5) % 1.0 - 0.5
- ka_y = (ActiveCell // L) % L+ delta_y = (a[1] - b[1] + 2.5) % 1.0 - 0.5
- ka_x = (ActiveCell - L * ka_y) % L+ return norm (delta_x, delta_y)
- delx = (kt_x + ka_x) % L+
- dely = (kt_y + ka_y) % L+
- return delx + dely * L+
- def CellIt(a): # Cell index for x,y positions+ def random_exponential (rate):
- return int((a[0] % L) * L) + L * int((a[1] % L) * L)+ """sample an exponential random number with given rate parameter"""
 + return -math.log (random.uniform (0.0, 1.0)) / rate
- def CellLimit(CellNumber): # Returns Rightmost x position. Can be used for x and y+ def pair_event_rate (delta_x, delta_y):
- return (CellNumber % L + 1) / float(L)+ """compute the particle event rate for the 1/r potential in 2D (lattice-screened version)"""
 + q = 0.0
 + for ky in range (-k_max, k_max + 1):
 + for kx in range (-k_max, k_max + 1):
 + q += (delta_x + kx) / norm (delta_x + kx, delta_y + ky) ** 3
 + q += 1.0 / norm (delta_x + k_max + 0.5, delta_y + ky)
 + q -= 1.0 / norm (delta_x - k_max - 0.5, delta_y + ky)
 + return max (0.0, q)
 +
 + def translated_cell (target_cell, active_cell):
 + """translate target_cell with respect to active_cell"""
 + kt_y = target_cell // L
 + kt_x = target_cell % L
 + ka_y = active_cell // L
 + ka_x = active_cell % L
 + del_x = (kt_x + ka_x) % L
 + del_y = (kt_y + ka_y) % L
 + return del_x + L*del_y
- def dist(x,y):+ def cell_containing (a):
- """periodic distance between two two-dimensional points x and y"""+ """return the index of the cell which contains the point a"""
- d_x= abs(x[0] - y[0]) % 1.0+ k_x = int (a[0] * L)
- d_x = min(d_x, 1.0 - d_x)+ k_y = int (a[1] * L)
- d_y= abs(x[1] - y[1]) % 1.0+ return k_x + L*k_y
- d_y = min(d_y, 1.0 - d_y)+
- return math.sqrt(d_x**2 + d_y**2)+
- def WalkerSet(pi_in):+ def walker_setup (pi):
- pi = deepcopy(pi_in)+ """compute the lookup table for Walker's algorithm"""
N_walker = len(pi) N_walker = len(pi)
walker_mean = sum(a[0] for a in pi) / float(N_walker) walker_mean = sum(a[0] for a in pi) / float(N_walker)
Line 71: Line 69:
for p in pi: for p in pi:
if p[0] > walker_mean: if p[0] > walker_mean:
- long_s.append(p)+ long_s.append (p[:])
else: else:
- short_s.append(p)+ short_s.append (p[:])
walker_table = [] walker_table = []
for k in range(N_walker - 1): for k in range(N_walker - 1):
Line 90: Line 88:
return N_walker, walker_mean, walker_table return N_walker, walker_mean, walker_table
- def WalkerSample(walker_table, walker_mean, N_walker):+ def sample_cell_veto (active_cell):
- Upsilon = random.uniform(0.0, walker_mean)+ """determine the cell which raised the cell veto"""
- i = random.randint(0, N_walker - 1)+ # first sample the distance vector using Walker's algorithm
 + i = random.randint (0, N_walker - 1)
 + Upsilon = random.uniform (0.0, walker_mean)
if Upsilon < walker_table[i][0]: if Upsilon < walker_table[i][0]:
- return walker_table[i][1]+ veto_offset = walker_table[i][1]
- else: return walker_table[i][2]+ else:
 + veto_offset = walker_table[i][2]
 + # translate with respect to active cell
 + veto_rate = Q_cell[veto_offset][0]
 + vetoing_cell = translated_cell (veto_offset, active_cell)
 + return vetoing_cell, veto_rate
 +
 +
 + N = 40
 + k_max = 3 # extension of periodic images.
 + chain_ell = 0.18 # displacement during one chain
 + L = 10 # number of cells along each dimension
 + density = N / 1.
 + cell_side = 1.0 / L
 + # precompute the cell-veto rates
 + cell_boundary = []
 + cb_discret = 10 # going around the boundary of a cell (naive)
 + for i in range (cb_discret):
 + x = i / float (cb_discret)
 + cell_boundary += [(x*cell_side, 0.0), (cell_side, x*cell_side),
 + (cell_side - x*cell_side, cell_side),
 + (0.0, cell_side - x*cell_side)]
 +
 + excluded_cells = [ del_x + L*del_y for del_x in (0, 1, L-1) \
 + for del_y in (0, 1, L-1) ]
 + Q_cell = []
 +
 + for del_y in xrange (L):
 + for del_x in xrange (L):
 + k = del_x + L*del_y
 + Q = 0.0
 + # "nearby" cells have no cell vetos
 + if k not in excluded_cells:
 + # scan the cell boundaries of both active and target cells
 + # to find the maximum of event rate
 + for delta_a in cell_boundary:
 + for delta_t in cell_boundary:
 + delta_x = del_x*cell_side + delta_t[0] - delta_a[0]
 + delta_y = del_y*cell_side + delta_t[1] - delta_a[1]
 + Q = max (Q, pair_event_rate (delta_x, delta_y))
 + Q_cell.append ([Q, k])
 +
 + Q_tot = sum (a[0] for a in Q_cell)
 + N_walker, walker_mean, walker_table = walker_setup (Q_cell)
 +
 + # histogram for computing g(r)
 + hbins = 50
 + histo = np.zeros (hbins)
 + histo_binwid = .5 / hbins
 + hsamples = 0
 +
 + # random initial configuration
 + particles = [ (random.uniform (0.0, 1.0), random.uniform (0.0, 1.0))
 + for _ in xrange (N) ]
- CellBoundary = []+ for iter in xrange (250000):
- NStep = 10 # going around the boundary of a cell (naive)+ if iter % 100 == 0:
- for i in range(NStep):+ print iter
- x = i / float(NStep)+
- CellBoundary += [(x, 0.0), (1.0, x), (1.0 - x, 1.0), (0.0, 1.0 - x)]+
- histo = []+ # possibly exchange x and y coordinates for ergodicity
- L = 10+
- NCell = L ** 2+
- NPart = 40+
- beta = 1.0+
- k_max = 2 # extension of periodic images.+
- CellLL = [(x / float(L), y / float(L)) for y in range(L) for x in range(L)]+
- CellExclude = [0, 1, L - 1, L, L + 1, 2 * L - 1, NCell - L, NCell - L + 1, NCell - 1]+
- QCell = []+
- for k in range(NCell):+
- if k in CellExclude: QCell.append([0, k])+
- else: Dummy = QCell.append([CellRate(CellLL[k], (0.0, 0.0), L,+
- CellBoundary), k])+
- QPrime = sum(a[0] for a in QCell) # cell rate+
- N_walker, walker_mean, walker_table = WalkerSet(QCell)+
- Particles = []+
- for k in range(NPart):+
- a = (random.uniform(0.0, 1.0), random.uniform(0.0, 1.0))+
- Particles.append(a)+
- for iter in range(10000):+
if random.randint(0,1) == 1: if random.randint(0,1) == 1:
- for k in range(NPart):+ particles = [ (y,x) for (x,y) in particles ]
- Particles[k] = (Particles[k][1], Particles[k][0])+ # pick active particle for first move
- Surplus = []+ active_particle = random.choice (particles)
- CellOcc = {}+ particles.remove (active_particle)
- ltilde = 0.18+ active_cell = cell_containing (active_particle)
- for k in range(NPart):+ # put particles into cells
- a = Particles[k]+ surplus = []
- n_cell = CellIt(a)+ cell_occupant = [ None ] * L * L
- if not CellOcc.has_key(n_cell):+ for part in particles:
- CellOcc[n_cell] = a+ k = cell_containing (part)
 + if cell_occupant[k] is None:
 + cell_occupant[k] = part
else: else:
- Surplus.append(a)+ surplus.append (part)
- if random.uniform(0.0, 1.0) < len(Surplus) / float(NPart):+
- ActiveParticle = random.choice(Surplus)+ # run one event chain
- Surplus.remove(ActiveParticle) # Active particle into Cell, not into Surplus+ distance_to_go = chain_ell
- ActiveCell = CellIt(ActiveParticle)+
- if CellOcc.has_key(ActiveCell):+
- Surplus.append(CellOcc.pop(ActiveCell))+
- CellOcc[ActiveCell] = ActiveParticle[:]+
- else:+
- while True:+
- ActiveCell = random.randint(0, NCell - 1)+
- if CellOcc.has_key(ActiveCell):+
- ActiveParticle = CellOcc[ActiveCell]+
- break+
- ActiveCellLimit = CellLimit(ActiveCell)+
- distance_to_go = ltilde+
while distance_to_go > 0.0: while distance_to_go > 0.0:
- PossibleActiveParticle= ActiveParticle[:]+ planned_event_type = 'end-of-chain'
- Possible_distance_to_go = distance_to_go+ planned_displacement = distance_to_go
- ActiveCellChange = False+ target_particle = None
- Lifting = False+ target_cell = None
- while True:+
- DistanceLimit = PossibleActiveParticle[0] + Possible_distance_to_go+ active_cell_limit = cell_side * (active_cell % L + 1)
- DelS =-math.log(random.uniform(0.0, 1.0)) / QPrime+ if active_cell_limit - active_particle[0] <= planned_displacement:
- if DistanceLimit < ActiveCellLimit and PossibleActiveParticle[0] + DelS > DistanceLimit:+ planned_event_type = 'active-cell-change'
- PossibleActiveParticle = (DistanceLimit, PossibleActiveParticle[1])+ planned_displacement = active_cell_limit - active_particle[0]
- Possible_distance_to_go = 0.0+
- break # Distance-to-go break+ delta_s = random_exponential (Q_tot)
- elif PossibleActiveParticle[0] + DelS > ActiveCellLimit:+ while delta_s < planned_displacement:
- Possible_distance_to_go -= (ActiveCellLimit - PossibleActiveParticle[0])+ vetoing_cell, veto_rate = sample_cell_veto (active_cell)
- PossibleActiveParticle = (ActiveCellLimit % 1.0, PossibleActiveParticle[1])+ part = cell_occupant[vetoing_cell]
- ActiveCellChange = True+ if part is not None:
- break #AC break+ Ratio = pair_event_rate (part[0] - active_particle[0] - delta_s, \
- else:+ part[1] - active_particle[1]) \
- PossibleActiveParticle = (PossibleActiveParticle[0] + DelS, PossibleActiveParticle[1])+ / veto_rate
- Possible_distance_to_go -= DelS+ if random.uniform (0.0, 1.0) < Ratio:
- TargetCell = WalkerSample(walker_table, walker_mean, N_walker)+ planned_event_type = 'particle'
- cell_rate_active_target = QCell[TargetCell][0]+ planned_displacement = delta_s
- TargetCell = CellTranslate(TargetCell, ActiveCell)+ target_particle = part
- if CellOcc.has_key(TargetCell):+ target_cell = vetoing_cell
- TargetParticle = CellOcc[TargetCell]+
- Ratio = PairRate(TargetParticle[0] - PossibleActiveParticle[0], TargetParticle[1] -+
- PossibleActiveParticle[1], k_max) / cell_rate_active_target+
- delx = (TargetParticle[0] - PossibleActiveParticle[0]) % 1.0+
- dely = (TargetParticle[1] - PossibleActiveParticle[1]) % 1.0+
- if random.uniform(0.0, 1.0) < Ratio:+
- Lifting = True+
- break # Lifting break+
- #+
- # here the sr (naive and a bit approximate, as we suppose a constant rate)+
- #+
- ToBeChecked = Surplus[:]+
- for ECell in CellExclude:+
- DummyCell = CellTranslate(ECell, ActiveCell)+
- if CellOcc.has_key(DummyCell):+
- ToBeChecked.append(CellOcc[DummyCell])+
- ToBeChecked.remove(ActiveParticle)+
- DelSMax = PossibleActiveParticle[0] - ActiveParticle[0]+
- for PossibleTargetParticle in ToBeChecked:+
- QRateLoc = PairRate(PossibleTargetParticle[0] - ActiveParticle[0],+
- PossibleTargetParticle[1] - ActiveParticle[1], k_max)+
- if QRateLoc > 0.0:+
- DelS =-math.log(random.uniform(0.0, 1.0)) / QRateLoc+
- if DelS < DelSMax: # Displacement cannot+
- # interfere with cell boundaries or distance_to_go+
- Lifting = True+
- ActiveCellChange = False+
- Possible_distance_to_go = distance_to_go - DelS+
- DelSMax = DelS+
- TargetParticle = PossibleTargetParticle[:] # have to take into+
- # account that the TargetParticle may be a Surplus one...+
- PossibleActiveParticle = (ActiveParticle[0] + DelS,+
- PossibleActiveParticle[1])+
- ActiveParticle = PossibleActiveParticle[:] #First move, then lift+
- distance_to_go = Possible_distance_to_go+
- if ActiveCellChange:+
- NewActiveCell = CellIt((ActiveParticle[0] + 0.001, ActiveParticle[1])) # Naive+
- if CellOcc.has_key(NewActiveCell):+
- Surplus.append(CellOcc.pop(NewActiveCell))+
- CellOcc[NewActiveCell] = ActiveParticle[:]+
- CellOcc.pop(ActiveCell) # Active cell occupied by active particle+
- for a in Surplus:+
- if CellIt(a) == ActiveCell:+
- CellOcc[ActiveCell] = a+
- Surplus.remove(a)+
break break
- ActiveCell = NewActiveCell+ delta_s += random_exponential (Q_tot)
- ActiveCellLimit = CellLimit(ActiveCell)+
- else:+ # compile the list of particles that need separate treatment
- CellOcc[ActiveCell] = ActiveParticle[:]+ extra_particles = surplus[:]
- if Lifting:+ for k in excluded_cells:
- if TargetParticle in Surplus:+ part = cell_occupant[translated_cell (k, active_cell)]
- TargetCell = CellIt(TargetParticle)+ if part is not None:
- Surplus.remove(TargetParticle)+ extra_particles.append (part)
- if CellOcc.has_key(TargetCell):+
- Surplus.append(CellOcc.pop(TargetCell))+ # naive version of the short-range code by discretization
- CellOcc[TargetCell] = TargetParticle[:]+ delta_s = 0.0
- ActiveParticle = TargetParticle[:]+ short_range_step = 1e-3
- ActiveCell = CellIt(ActiveParticle) # Naive , zu verbessern+ while delta_s < planned_displacement:
- ActiveCellLimit = CellLimit(ActiveCell)+ for possible_target_particle in extra_particles:
- # Naive, Particles vector for x <-> y transfer+ # this supposes a constant event rate over the time interval
- Particles = []+ # [delta_s:delta_s+short_range_step]
- for k in range(NCell):+ q = pair_event_rate (possible_target_particle[0] - active_particle[0] - delta_s,
- if CellOcc.has_key(k):+ possible_target_particle[1] - active_particle[1])
- Particles.append(CellOcc[k])+ if q > 0.0:
- Particles += Surplus+ event_time = random_exponential (q)
- for k in range(NPart):+ if event_time < short_range_step and delta_s + event_time < planned_displacement:
- for l in range(k):+ planned_event_type = 'particle'
- histo.append(dist(Particles[k], Particles[l]))+ planned_displacement = delta_s + event_time
- pylab.hist(histo, bins=100, range=(0.0, 1.0), normed=True)+ target_particle = possible_target_particle
- pylab.title('Demo_cell, ECMC, $k_{\max}$ = ' + str(k_max) + ' $NPart$ = ' ++ target_cell = cell_containing (target_particle)
- str(NPart))+ break
- pylab.savefig('eventchain.png')+ delta_s += short_range_step
- pylab.show()+
 + # advance active particle
 + distance_to_go -= planned_displacement
 + new_x = active_particle[0] + planned_displacement
 + active_particle = (new_x % 1.0, active_particle[1])
 +
 + if planned_event_type == 'active-cell-change':
 + ac_x = (active_cell_limit + 0.5*cell_side) % 1.0
 + active_cell = cell_containing ([ac_x, active_particle[1]])
 + active_particle = (active_cell % L * cell_side, active_particle[1])
 +
 + elif planned_event_type == 'particle':
 + # remove newly active particle from store
 + if target_particle in surplus:
 + surplus.remove (target_particle)
 + else:
 + cell_occupant[target_cell] = None
 + # put the previously active particle in the store
 + if cell_occupant[active_cell] is not None:
 + surplus.append (active_particle)
 + else:
 + cell_occupant[active_cell] = active_particle
 + active_particle = target_particle
 + active_cell = cell_containing (active_particle)
 +
 + # restore particles vector for x <-> y transfer
 + particles = [ active_particle ]
 + particles += [ part for part in cell_occupant if part is not None ]
 + particles += surplus
 + 
 + # form histogram for computing radial distribution function g(r)
 + for k in range (len (particles)):
 + for l in range (k):
 + ibin = int (dist (particles[k], particles[l]) / histo_binwid)
 + if ibin < len (histo):
 + histo[ibin] += 1
 + hsamples += 1
 +
 + # compute g(r) from histogram
 + half_bin = .5 * histo_binwid
 + r = np.arange (0., hbins) * histo_binwid + half_bin
 + g_of_r = histo / density / hsamples * 2
 + g_of_r /= math.pi * ((r+half_bin)**2 - (r-half_bin)**2)
 + # save g(r)
 + np.savetxt ('cvmc-radial-distr-func.dat', zip (r, g_of_r))

Revision as of 22:08, 25 September 2016

S. C. Kapfer, W. Krauth Cell-veto Monte Carlo algorithm for long-range systems Physical Review E 94, 031302(R) (2016)

Paper

Abstract We present a rigorous efficient event-chain Monte Carlo algorithm for long-range interacting particle systems. Using a cell-veto scheme within the factorized Metropolis algorithm, we compute each single-particle move with a fixed number of operations. For slowly decaying potentials such as Coulomb interactions, screening line charges allow us to take into account periodic boundary conditions. We discuss the performance of the cell-veto Monte Carlo algorithm for general inverse-power-law potentials, and illustrate how it provides a new outlook on one of the prominent bottlenecks in large-scale atomistic Monte Carlo simulations.

Electronic version (from arXiv)

Journal version (subscription required)

Illustration

Here are pictures of three bosons on an permutation cycle.
The cell-veto algorithm is an application of the event-chain paradigm that was developed in 2009, with E. P. Bernard and D. B. Wilson, and much extended together with M. Michel and S. C. Kapfer, where we introduced the factorized Metropolis algorithm. In this algorithm, unlike in 99.9% of all simulation algorithms (Monte Carlo or Molecular dynamics), one does not compute the system energy in order to decide on a change of the physical system, but rather looks at all the interactions separately. So, if a particle a (the active particle) wants to move, it has to ask all its partners t_1, t_2, .... (the target particles). If there is only a single veto, the move is rejected. In the cell-veto algorithm (see the right side of the figure), the identification of the rejecting particle is preceeded by that of a veto cell. The advantage of this is that cell vetos can be identified immediately (in a constant number of operations, that is, in O(1)), and then instantly confirmed or infirmed on the particle level.

Python implementation

The paper contains a Python demo implementation of the cell-veto algorithm (that is also available on github). 

 import math, random, sys
import numpy as np

def norm (x, y):
    """norm of a two-dimensional vector"""
    return (x*x + y*y) ** 0.5

def dist (a, b):
   """periodic distance between two two-dimensional points a and b"""
   delta_x = (a[0] - b[0] + 2.5) % 1.0 - 0.5
   delta_y = (a[1] - b[1] + 2.5) % 1.0 - 0.5
   return norm (delta_x, delta_y)

def random_exponential (rate):
    """sample an exponential random number with given rate parameter"""
    return -math.log (random.uniform (0.0, 1.0)) / rate

def pair_event_rate (delta_x, delta_y):
    """compute the particle event rate for the 1/r potential in 2D (lattice-screened version)"""
    q = 0.0
    for ky in range (-k_max, k_max + 1):
        for kx in range (-k_max, k_max + 1):
            q += (delta_x + kx) / norm (delta_x + kx, delta_y + ky) ** 3
        q += 1.0 / norm (delta_x + k_max + 0.5, delta_y + ky)
        q -= 1.0 / norm (delta_x - k_max - 0.5, delta_y + ky)
    return max (0.0, q)

def translated_cell (target_cell, active_cell):
    """translate target_cell with respect to active_cell"""
    kt_y = target_cell // L
    kt_x = target_cell  % L
    ka_y = active_cell // L
    ka_x = active_cell  % L
    del_x = (kt_x + ka_x) % L
    del_y = (kt_y + ka_y) % L
    return del_x + L*del_y

def cell_containing (a):
    """return the index of the cell which contains the point a"""
    k_x = int (a[0] * L)
    k_y = int (a[1] * L)
    return k_x + L*k_y

def walker_setup (pi):
    """compute the lookup table for Walker's algorithm"""
    N_walker = len(pi)
    walker_mean = sum(a[0] for a in pi) / float(N_walker)
    long_s = []
    short_s = []
    for p in pi:
        if p[0] > walker_mean:
            long_s.append (p[:])
        else:
            short_s.append (p[:])
    walker_table = []
    for k in range(N_walker - 1):
        e_plus = long_s.pop()
        e_minus = short_s.pop()
        walker_table.append((e_minus[0], e_minus[1], e_plus[1]))
        e_plus[0] = e_plus[0] - (walker_mean - e_minus[0])
        if e_plus[0] < walker_mean:
            short_s.append(e_plus)
        else:
            long_s.append(e_plus)
    if long_s != []:
        walker_table.append((long_s[0][0], long_s[0][1], long_s[0][1]))
    else:
        walker_table.append((short_s[0][0], short_s[0][1], short_s[0][1]))
    return N_walker, walker_mean, walker_table

def sample_cell_veto (active_cell):
    """determine the cell which raised the cell veto"""
    # first sample the distance vector using Walker's algorithm
    i = random.randint (0, N_walker - 1)
    Upsilon = random.uniform (0.0, walker_mean)
    if Upsilon < walker_table[i][0]:
        veto_offset = walker_table[i][1]
    else:
        veto_offset = walker_table[i][2]
    # translate with respect to active cell
    veto_rate = Q_cell[veto_offset][0]
    vetoing_cell = translated_cell (veto_offset, active_cell)
    return vetoing_cell, veto_rate


N = 40
k_max = 3 # extension of periodic images.
chain_ell = 0.18  # displacement during one chain
L = 10  # number of cells along each dimension
density = N / 1.
cell_side = 1.0 / L
# precompute the cell-veto rates
cell_boundary = []
cb_discret = 10 # going around the boundary of a cell (naive)
for i in range (cb_discret):
    x = i / float (cb_discret)
    cell_boundary += [(x*cell_side, 0.0), (cell_side, x*cell_side),
                      (cell_side - x*cell_side, cell_side),
                      (0.0, cell_side - x*cell_side)]

excluded_cells = [ del_x + L*del_y for del_x in (0, 1, L-1) \
                                   for del_y in (0, 1, L-1) ]
Q_cell = []

for del_y in xrange (L):
    for del_x in xrange (L):
        k = del_x + L*del_y
        Q = 0.0
        # "nearby" cells have no cell vetos
        if k not in excluded_cells:
            # scan the cell boundaries of both active and target cells
            # to find the maximum of event rate
            for delta_a in cell_boundary:
                for delta_t in cell_boundary:
                    delta_x = del_x*cell_side + delta_t[0] - delta_a[0]
                    delta_y = del_y*cell_side + delta_t[1] - delta_a[1]
                    Q = max (Q, pair_event_rate (delta_x, delta_y))
        Q_cell.append ([Q, k])

Q_tot = sum (a[0] for a in Q_cell)
N_walker, walker_mean, walker_table = walker_setup (Q_cell)

# histogram for computing g(r)
hbins = 50
histo = np.zeros (hbins)
histo_binwid = .5 / hbins
hsamples = 0

# random initial configuration
particles = [ (random.uniform (0.0, 1.0), random.uniform (0.0, 1.0))
    for _ in xrange (N) ]

for iter in xrange (250000):
    if iter % 100 == 0:
        print iter

    # possibly exchange x and y coordinates for ergodicity
    if random.randint(0,1) == 1:
        particles = [ (y,x) for (x,y) in particles ]
    # pick active particle for first move
    active_particle = random.choice (particles)
    particles.remove (active_particle)
    active_cell = cell_containing (active_particle)
    # put particles into cells
    surplus = []
    cell_occupant = [ None ] * L * L
    for part in particles:
        k = cell_containing (part)
        if cell_occupant[k] is None:
            cell_occupant[k] = part
        else:
            surplus.append (part)

    # run one event chain
    distance_to_go = chain_ell
    while distance_to_go > 0.0:
        planned_event_type = 'end-of-chain'
        planned_displacement = distance_to_go
        target_particle = None
        target_cell = None

        active_cell_limit = cell_side * (active_cell % L + 1)
        if active_cell_limit - active_particle[0] <= planned_displacement:
            planned_event_type = 'active-cell-change'
            planned_displacement = active_cell_limit - active_particle[0]

        delta_s = random_exponential (Q_tot)
        while delta_s < planned_displacement:
            vetoing_cell, veto_rate = sample_cell_veto (active_cell)
            part = cell_occupant[vetoing_cell]
            if part is not None:
                Ratio = pair_event_rate (part[0] - active_particle[0] - delta_s, \
                                         part[1] - active_particle[1])           \
                        / veto_rate
                if random.uniform (0.0, 1.0) < Ratio:
                    planned_event_type = 'particle'
                    planned_displacement = delta_s
                    target_particle = part
                    target_cell = vetoing_cell
                    break
            delta_s += random_exponential (Q_tot)

        # compile the list of particles that need separate treatment
        extra_particles = surplus[:]
        for k in excluded_cells:
            part = cell_occupant[translated_cell (k, active_cell)]
            if part is not None:
                extra_particles.append (part)

        # naive version of the short-range code by discretization
        delta_s = 0.0
        short_range_step = 1e-3
        while delta_s < planned_displacement:
            for possible_target_particle in extra_particles:
                # this supposes a constant event rate over the time interval
                # [delta_s:delta_s+short_range_step]
                q = pair_event_rate (possible_target_particle[0] - active_particle[0] - delta_s,
                                     possible_target_particle[1] - active_particle[1])
                if q > 0.0:
                    event_time = random_exponential (q)
                    if event_time < short_range_step and delta_s + event_time < planned_displacement:
                        planned_event_type = 'particle'
                        planned_displacement = delta_s + event_time
                        target_particle = possible_target_particle
                        target_cell = cell_containing (target_particle)
                        break
            delta_s += short_range_step

        # advance active particle
        distance_to_go -= planned_displacement
        new_x = active_particle[0] + planned_displacement
        active_particle = (new_x % 1.0, active_particle[1])

        if planned_event_type == 'active-cell-change':
            ac_x = (active_cell_limit + 0.5*cell_side) % 1.0
            active_cell = cell_containing ([ac_x, active_particle[1]])
            active_particle = (active_cell % L * cell_side, active_particle[1])

        elif planned_event_type == 'particle':
            # remove newly active particle from store
            if target_particle in surplus:
                surplus.remove (target_particle)
            else:
                cell_occupant[target_cell] = None
            # put the previously active particle in the store
            if cell_occupant[active_cell] is not None:
                surplus.append (active_particle)
            else:
                cell_occupant[active_cell] = active_particle
            active_particle = target_particle
            active_cell = cell_containing (active_particle)

    # restore particles vector for x <-> y transfer
    particles = [ active_particle ]
    particles += [ part for part in cell_occupant if part is not None ]
    particles += surplus

    # form histogram for computing radial distribution function g(r)
    for k in range (len (particles)):
        for l in range (k):
            ibin = int (dist (particles[k], particles[l]) / histo_binwid)
            if ibin < len (histo):
                histo[ibin] += 1
        hsamples += 1

# compute g(r) from histogram
half_bin = .5 * histo_binwid
r = np.arange (0., hbins) * histo_binwid + half_bin
g_of_r = histo / density / hsamples * 2
g_of_r /= math.pi * ((r+half_bin)**2 - (r-half_bin)**2)
# save g(r)
np.savetxt ('cvmc-radial-distr-func.dat', zip (r, g_of_r))
Retrieved from "http://www.lps.ens.fr/~krauth/index.php/Kapfer_Krauth_2016"
Personal tools