Vfrdtyky

I'm getting traceback issues with my kmc.time_step(), self.evolve(events_index, dt) and my ran_int = randint(len(index)) saying that zip has no length but have no clue how to go about solving this.

Im trying to model a kinetic monte carlo simulation by checking nearest neighbours of an array to see if they have any occupied states and then either "diffusing" through if they aren't occupied (moving to the unoccupied point). If the point is occupied then diffusion cant happen so aggregation occurs. I have a bunch of conditional if statements to check these but then why I try to run I get traceback errors and have no clue why.

from numpy import *

# from scipy.signal import correlate2d

from numpy.random import randint, choice, uniform

from matplotlib.pyplot import *





class KineticMonteCarlo(object):



    def __init__(self, Model):



        """Initialise a Kinetic Monte Carlo evolution method for

        the given model. The model object requires the following

        vectorised methods:



        Model.evolve : <method>

            evolves model according to given events

        Model.get_rates : <method>

            returns event rates

        """



        # reference state

        self.Model = Model

        self.state = Model.state



        # event rates and evolution rule

        self.evolve = Model.evolve

        self.get_rates = Model.get_rates



    def time_step(self):

        """Perform parallel Kinetic Monte-Carlo update on model"""



        # calculate transition rate fractions

        rates = self.get_rates()  # a k-dimensional array

        probabilities = cumsum(rates, axis=0)  # a k-dim array



        total_rate = probabilities[-1]

        probabilities /= total_rate



        # choose events according to rates

        u = uniform(0, 1)

        events_index = sum(probabilities < u)



        # generate waiting time

        v = uniform(0, 1)

        dt = -log(v)/total_rate



        # carry out events

        self.evolve(events_index, dt)





class LatticeGas(object):



    def __init__(self, n_sites, density):



        # initialisation of lattice gas with given density

        self.size = 2*int(ceil(sqrt(n_sites)/2))

        self.state = choice([0, 1], size=(self.size, self.size), p=[1-density, density])

        self.time = 0.

        self.aggcnt = 0

        self.diffcnt = 0

        self.n = 0



        ''' this return a a matrix of 0 and 1; 1 for occupied sites'''



    def get_rates(self):

        temp = 100

        kb = 1.38

        beta = 1/(kb * temp)

        ads_rate = 0.001 * exp(- beta * 0.01)

        des_rate = 0.001 * exp(- beta * 0.01)

        diff_rate = 0.1 * exp(- beta * 0.35)



        return array([ads_rate, des_rate, diff_rate])



    def evolve(self, events_index, dt):

        non_zero_m = self.state != 0

        non_zeroi, non_zeroj = where(non_zero_m)

        index = zip(non_zeroi, non_zeroj)



        ran_ind = randint(len(index))

        randi, randj = index[ran_ind]



        n1 = (randi - 1) % self.size

        n2 = (randj - 2) % self.size

        e1 = (randj + 1) % self.size

        e2 = (randj + 2) % self.size

        w1 = (randj - 1) % self.size

        w2 = (randj - 2) % self.size

        s1 = (randi + 1) % self.size

        s2 = (randi + 2) % self.size



        if events_index == 0:

            self.state[index[randint(len(index))]]



        elif events_index == 1:

            index = zip(*where(self.state == 0))

            self.state[index[randint(len(index))]]



        elif events_index == 2:

            while True:

                rand_direc = randint(4)



                if rand_direc == 3 and self.state[n1, randj] == 0:

                    self.state[n1, randj] = 1

                    self.state[randi, randj] = 0

                    if self.state[n2, randj % self.size == 0] and self.state[n1, w1] == 0 and self.state[n1, e1] == 0:

                        self.diffcnt += 1

                    else:

                        self.aggcnt == 1

                    break



                if rand_direc == 2 and self.state[randi, e1] == 0:

                    self.state[randi, e1] = 1

                    self.state[randi, randj] = 0

                    if self.state[n1, e1] and self.state[s1, e1] == 0 and self.state[randi, e2] == 0:

                        self.diffcnt += 1

                    else:

                        self.aggcnt += 1

                    break



                elif rand_direc == 1 and self.state[s1, randj] == 0:

                    self.state[s1, randj] = 1

                    self.state[randi, randj] = 0

                    if self.state[s1, e1] == 0 and self.state[s1, w1] == 0 and self.state[s2, randj] == 0:

                        self.diffcnt += 1

                    else:

                        self.aggcount += 1

                    break



                elif rand_direc == 0 and self.state[randi, w1] == 0:

                    self.state[randi, w1] = 1

                    self.state[randi, randj] = 0

                    if self.state[n1, w1] == 0 and self.state[s1, w1] == 0 and self.state[randi, w2] == 0:

                        self.diffcnt += 1

                    else:

                        self.aggcnt += 1

                    break



    def show(self):

        # Display current lattice configuration, with black

        # and white representing whether there is a particle or not



        figure(figsize=(5, 5))

        imshow(self.state, interpolation="none", cmap='bone_r', vmin=0, vmax=1)

        axis('off')

        show()





# example of how code should work

model = LatticeGas(50, 0.1)

kmc = KineticMonteCarlo(model)

steps = 100



for i in range(steps):

    kmc.time_step()

    model.show()