How to optimize this random entropy calculation for a large number of iterations?

What my code is trying to do :

I have an initial array containing an arrangement of molecules, each in a cell and confined to moving in a 2D array ( Up Down Left Right ) ( Array size:200x200 ). At each step I take a random molecule, and move it into a random adjacent cell.

Starting from a certain number and every few iterations, I calculate the entropy of this grid. The grid is cut into small squares of 25x25. I then use the Shannon entropy to calculate the entropy of the system.

Objective:

Doing 1e8+ iterations in a decent time using an i5-6500, no GPU.

My code:

function advance_multi_lattice(grid)   #Find the next state of the system

    rnd=rand(1:count(!iszero,grid))        #Random number to be used for a random molecule.
    slots=find(!iszero,grid)    #Cells containing molecules.
    chosen_slot=find(!iszero,grid)[rnd]    #Random cell. May contain multiple molecules.
    dim=size(grid)[1]    #Need this for rnd=3,4 later.  
    grid[chosen_slot]-=1    #Remove the molecules from the cell
    rnd_arr=[1,2,3,4]     #Array to random from.


    while true


        rnd=rand(rnd_arr)   #Random number to see which side should the molecules go.

        if rnd==1    #Right for example.
            try     #In case moving right is impossible, ie: moving right gets the molecule out. Remove 1 from rnd_arr and repeat.
                grid[chosen_slot+1]+=1
                break
            catch
                filter!(e->e!=1,rnd_arr)
            end
        elseif rnd==2
            try   #Same
                grid[chosen_slot-1]+=1
                break
            catch
                filter!(e->e!=2,rnd_arr)
            end
        #Repeat for the other numbers : 3 and 4...
    return Grid
end

function S(P)   #Entropy, if no molecules then return 0.
    s=[]
    for k in P
        if k==0
            push!(s,0)
        else
            push!(s,-k*log(k))
        end

    end
    return s
end

function find_molecules(grid) #How many molecules in the array
    s=0

    for slot in grid
        s+=slot
    end
    return s
end

function entropy_scale(grid,total_molecules)    #Calculate the entropy of the grid.
    P_array=Array{Float64}([])
    for i=1:8
        for j=1:8
            push!(P_array,find_molecules(grid[(i-1)*25+1:i*25,(j-1)*25+1:j*25]))
        end
    end

    P_array=P_array./total_molecules

    return sum(S(P_array))
end

function entropy_evolution(grid,n)    #The loop function. Changes the grid and returns the entropy as a function of steps.
    t_arr=Array{Int64}([])
    S_arr=Array{Float64}([])
    p=Progress(Int(n))    #Progress bar, using ProgressMeter.
    total_molecules=find_molecules(grid)
    for k=1:1e3
        grid=advance_multi_lattice(grid)

        next!(p)
    end

    for k=1e3+1:n
        grid=advance_multi_lattice(grid)

        if  k%500==0  #Only record entropy every 500 steps
            push!(S_arr,entropy_scale(grid,totel_molecules))
        end
        next!(p)
    end

    return S_arr,grid

end        

Results for my code :

For 1e5 iterations, I get 43 seconds. Which means that if I want an interesting result ( 1e9+ ), I need a lot of time, upwards to 1hour+. Changing the entropy calculation threshold barely scratches the performance unless it's really small.