% Kinetic Monte Carlo (KMC) Simulation of a Kossel Crystal with a single growth
% unit
% Jennifer Mills
% EPS 109, UC Berkeley, 2016

%Basic premise: represent crystal surface by a 2D array of integers Z, the
%value of which represents the height of the surface at that x,y point. The
%Monte Carlo algorithm consists of three steps:
% 1) Choose an (x,y) site to act upon. This is done by breaking the total
% number of sites down by the number of nearest neighbors (n - can range
% from 1 for an adatom that has no lateral neighbors, to 5 for an atom in
% the center of a terrace). The algorithm picks an n based on the
% probabilities of each of the n's and then randomly selects an (x,y) index
% from the list of sites with n nearest neighbors.
% 2) Choose whether to add a growth unit, subtract a growth unit, or have
% the growth unit diffuse on the surface (not implemented in this version
% of the code - contact the author for further details). The probability 
% of each action is proportional to the rate constant of that process.

clf ; clear ; close all

%initialize 2D matrix Z of crystal height
maxX = 40;
maxY = 40;
Z = 3*ones(maxX, maxY); %starting height of 3 for all points (perfect crystal)

% CRYSTAL SURFACE INITIAL CONDITION
% Continuous ledge
 Z(2*round(maxX/3):maxX,:) = Z(2*round(maxX/3):maxX,:) + 1;

% %crystal ledge with kink
%  Z(2*round(maxX/3):maxX,:) = Z(2*round(maxX/3):maxX,:) + 1;
%  Z(2*round(maxX/3) - 1,1:round(maxY/2)) = Z(2*round(maxX/3) - 1,1:round(maxY/2)) + 1;
%  Z(2*round(maxX/3) - 2,1:round(maxY/2)) = Z(2*round(maxX/3) - 2,1:round(maxY/2)) + 1;


% Define CRYSTAL GROWTH VARIABLES
% gamma  = average energy of a nearest-neighbor bond
eta = 4;  %gamma / kT
du = -3;   %delta mu / kT. Difference in chemical potential from equilibrium.
v = 1;    %frequency factor
%rates
kn1 = v*exp(-1*eta); %rate of removal from surface. Depends on # nearest neighbors
kn2 = v*exp(-2*eta);
kn3 = v*exp(-3*eta);
kn4 = v*exp(-4*eta);
kn5 = v*exp(-5*eta);
kp = exp(du)*kn3;    % rate of addition to the surface - independent of location

%initialize vectors for each number of nearest neighbors (Nn). Stores the
%indices of every grid point that has n nearest neighbors.
N1 = zeros(maxX*maxY, 2);
N2 = zeros(maxX*maxY, 2);
N3 = zeros(maxX*maxY, 2);
N4 = zeros(maxX*maxY, 2);
N5 = zeros(maxX*maxY, 2);

%Define simulation length
max_attempts = 20000;
attempts = 0;
f = 0; %counter for .png creation

%%%%%%%%
%%% KMC loop 
%%%%%%%%

while attempts < max_attempts

    %count nearest neighbors and assign each grid point to an Nn vector

    %clear the n counters for the next iteration.
    n1_count = 0;
    n2_count = 0;
    n3_count = 0;
    n4_count = 0;
    n5_count = 0;

    % Determine nearest neighbors
    for i = 1:maxX
        for j = 1:maxY
            %count number of nearest neighbors, taking periodic boundary
            %conditions into account.
            xR = i+1; % xRIGHT = x coordinate of my right neighbor
            xL = i-1; % xLEFT  = x coordinate of my left  neighbor
            yU = j+1; % yUP    = y coordinate of my upper neighor
            yD = j-1; % yDOWN  = y coordinate of my lower neighbor

            %Apply periodic boundary conditions to xL, xR, yU, and yD
            if xL < 1
                xL = maxX;
            end
            if xR > maxX;
                xR = 1;
            end
            if yU > maxY;
                yU = 1;
            end
            if yD < 1
                yD = maxY;
            end

            %count the number of nearest neighbors by checking to see if they 
            %have the same (or greater) Z value. n starts at one because each
            %growth unit must be placed on an existing growth unit.
            n = 1;
            if Z(xR,j) >= Z(i,j)
                n = n+1;
            end
            if Z(xL,j) >= Z(i,j)
                n = n+1;
            end
            if Z(i,yU) >= Z(i,j)
                n = n+1;
            end
            if Z(i,yD) >= Z(i,j)
                n = n+1;
            end

            %add the (i,j) point to the correct Nn vector.
            if n == 1
                n1_count = n1_count + 1;
                N1(n1_count, 1) = i;
                N1(n1_count, 2) = j;
            elseif n == 2
                n2_count = n2_count + 1;
                N2(n2_count, 1) = i;
                N2(n2_count, 2) = j;
            elseif n == 3
                n3_count = n3_count + 1;
                N3(n3_count, 1) = i;
                N3(n3_count, 2) = j;
            elseif n == 4
                n4_count = n4_count + 1;
                N4(n4_count, 1) = i;
                N4(n4_count, 2) = j;
            elseif n == 5
                n5_count = n5_count + 1;
                N5(n5_count, 1) = i;
                N5(n5_count, 2) = j;  
            end
        end
    end

    % Calculate probabilities of each n
    n_sum = (n1_count * (kn1 + kp)) + (n2_count * (kn2 + kp)) + (n3_count * (kn3 + kp)) + (n4_count * (kn4 + kp)) + (n5_count * (kn5 + kp));
    P1 = (n1_count * (kn1 + kp))/n_sum;
    P2 = (n2_count * (kn2 + kp))/n_sum;
    P3 = (n3_count * (kn3 + kp))/n_sum;
    P4 = (n4_count * (kn4 + kp))/n_sum;
    P5 = (n5_count * (kn5 + kp))/n_sum;
    Ptot = P1 + P2 + P3 + P4 + P5; 

    %check Ptot = 1 - if not, display an error message!
    if abs(Ptot - 1) >= 1e-3
        fprintf('!!Probabilities for n selection not equal to one!!');
    end

    %select an n, select the x,y point to act upon, and then select and
    %implement action (add, remove, or (optional) diffuse growth unit).
    p_check = rand();
    if p_check < P1
        select = n1_count * rand();
        x = N1(ceil(select), 1);
        y = N1(ceil(select), 2);
        %select action
        Pplus = kp / (kn1 + kp);
        Pminus = kn1 / (kn1 + kp);
        %implement!
        action_check = rand();
        if action_check < Pplus
            Z(x,y) = Z(x,y) + 1;
        else
            Z(x,y) = Z(x,y) - 1;
        end
    elseif p_check < P1 + P2
        select = n2_count * rand();
        x = N2(ceil(select), 1);
        y = N2(ceil(select), 2);
        %select action
        Pplus = kp / (kn2 + kp);
        Pminus = kn2 / (kn2 + kp);
        %implement!
        action_check = rand();
        if action_check < Pplus
            Z(x,y) = Z(x,y) + 1;
        else
            Z(x,y) = Z(x,y) - 1;
        end
    elseif p_check < P1 + P2 + P3
        select = n3_count * rand();
        x = N3(ceil(select), 1);
        y = N3(ceil(select), 2);  
        %select action
        Pplus = kp / (kn3 + kp);
        Pminus = kn3 / (kn3 + kp);
        %implement!
        action_check = rand();
        if action_check < Pplus
            Z(x,y) = Z(x,y) + 1;
        else
            Z(x,y) = Z(x,y) - 1;
        end
    elseif p_check < P1 + P2 + P3 + P4
        select = n4_count * rand();
        x = N4(ceil(select), 1);
        y = N4(ceil(select), 2);
        %select action
        Pplus = kp / (kn4 + kp);
        Pminus = kn4 / (kn4 + kp);
        %implement!
        action_check = rand();
        if action_check < Pplus
            Z(x,y) = Z(x,y) + 1;
        else
            Z(x,y) = Z(x,y) - 1;
        end
    else
        select = n5_count * rand();
        x = N5(ceil(select), 1);
        y = N5(ceil(select), 2); 
        %select action
        Pplus = kp / (kn5 + kp);
        Pminus = kn5 / (kn5 + kp);
        %implement!
        action_check = rand();
        if action_check < Pplus
            Z(x,y) = Z(x,y) + 1;
        else
            Z(x,y) = Z(x,y) - 1;
        end
    end

    attempts = attempts + 1;

    %plot the crystal surface
    %plot every 100th attempt and save the resulting image for .gif
    %creation
    if mod(attempts+99,100) == 0
        f = f+1;
        xx = 0:maxX - 1;
        yy = 0:maxY - 1;
        [X,Y] = meshgrid(xx,yy);
        surf(X,Y,Z);
        %surf(X,Y,Z, 'EdgeColor', 'none');
        caxis([0 14]);
        axis([0 maxX 0 maxY 0 14]);
        h = colorbar;
        ylabel(h, 'surface height')
        title('\Delta\mu/kT = 1', 'FontSize', 20);
        drawnow;
        %print(gcf, ['growth_u1_',num2str(f, '%03d'), '.png'], '-dpng', '-r100')
        pause(0.01);
    end

    fprintf('Finished %10.0d attempts \n', attempts);
end
        