MATLAB: Increase code efficiency in logic statements (Metropolis method)

Question

I have a code which is designed to simulate the motion of an ideal gas in a box. It is based on the metropolis method in a Monte Carlo simulation. However, I used a series of logic statements (mostly ifs) to define boundary conditions when finding adjacent particles to a randomly selected particle. The algorithm denotes a 1 for any adjacent particle inside a 1x4 matrix, and a 0 for an adjacent spot with no particle. I need the algorithm to automatically set any adjacent points outside the box to 0 for any particle on the edge of the box. Is there a way to reduce these logic statements and still get the same results? The code is below.

% define constants
beta = .01; %Inverse temperature
N=2000; %Duration of simulation
eps = -7;
mu=9;


for j=1:N
    a=randi(L);
    b=randi(L);
    c=randi(L);
    d=randi(L);
    % Calculate energy at positions
    if lattice(a,b)==1 && lattice(c,d)==0
        %If distribution (according to energy) suggests move,
              lattice(a,b)=0;
            lattice(c,d)=1;

         %Energy at random site/position
       E = n*(mu-eps)
            if (a~=1&& a~=L &&b~=1 && b~=L)
               adjacent= [lattice(a-1,b) lattice(a+1,b) lattice(a,b+1) lattice(a,b-1) 1];
            else if (a==1 && b==1)
                    adjacent= [0 lattice(a+1,b) lattice(a,b+1) 0 1];
                else if (a==L && b==L)
                        adjacent= [lattice(a-1,b) 0 0 lattice(a,b-1) 1];
                    else if(a==1&&b==L)
                            adjacent= [0 lattice(a+1,b) 0 lattice(a,b-1) 1];
                        else if (a==L &&b==1)
                                adjacent= [lattice(a-1,b) 0 lattice(a,b+1) 0 1];
                            else if (a==1 && b~=L && b~=1)
                                    adjacent= [0 lattice(a+1,b) lattice(a,b+1) lattice(a,b-1) 1];
                                 else if (a==L && b~=L && b~=1)
                                        adjacent= [lattice(a-1,b) 0 lattice(a,b+1) lattice(a,b-1) 1];
                                     else if (b==1&&a~=L&&a~=1)
                                             adjacent= [lattice(a-1,b) lattice(a+1,b) lattice(a,b+1) 0 1];
                                         else if (b==L&&a~=L&&a~=1)
                                              adjacent= [lattice(a-1,b) lattice(a+1,b) 0 lattice(a,b-1) 1];  
                                             end
                                         end
                                     end
                                end
                            end
                        end
                    end
                end
            end


            E1 = mu*sum(adjacent) + eps*sum(sum(adjacent.*adjacent)); 
            %This calculates the energy of the particle at its current
            %position

            if (c~=1&& c~=L &&d~=1 && d~=L)
               adjacent1= [lattice(c-1,d) lattice(c+1,d) lattice(c,d+1) lattice(c,d-1) 1];
            else if (c==1 && d==1)
                    adjacent1= [0 lattice(c+1,d) lattice(c,d+1) 0 1];
                else if (c==L && d==L)
                        adjacent1= [lattice(c-1,d) 0 0 lattice(c,d-1) 1];
                    else if(c==1&&d==L)
                            adjacent1= [0 lattice(c+1,d) 0 lattice(c,d-1) 1];
                        else if (c==L &&d==1)
                                adjacent1= [lattice(c-1,d) 0 lattice(c,d+1) 0 1];
                            else if (c==1 && d~=L && d~=1)
                                    adjacent1= [0 lattice(c+1,d) lattice(c,d+1) lattice(c,d-1) 1];
                                 else if (c==L && d~=L && d~=1)
                                        adjacent1= [lattice(c-1,d) 0 lattice(c,d+1) lattice(c,d-1) 1];
                                     else if (d==1&&c~=L&&c~=1)
                                             adjacent1= [lattice(c-1,d) lattice(c+1,d) lattice(c,d+1) 0 1];
                                         else if (d==L&&c~=L&&c~=1)
                                              adjacent1= [lattice(c-1,d) lattice(c+1,d) 0 lattice(c,d-1) 1];  
                                             end
                                         end
                                     end
                                end
                            end
                        end
                    end
                end
            end


            E2 = mu*sum(adjacent) + eps*sum(sum(adjacent1.*adjacent1)); 
            %Calculates the energy at randomly chosen position.

            dE = E2-E1; %Change in energy of the particle as it goes from the two locations.

            if rand<exp(-beta*dE)
                lattice(a,b)=0;
                lattice(c,d)=1;
            end



    end
 time(:,:,j)=lattice;


end

Answer 1

You can do away with all those if else statements. In effect you are only evaluating four logical expressions and selecting indices based on that. You can use the following

ab = [a-1, b; a+1, b; a, b+1; a, b-1];
abIn = [a ~= 1; a ~= L; b ~= L; b ~= 1];
adjacent = zeros(1, 5);
adjacent(abIn) = lattice(sub2ind([L, L], ab(abIn, 1), ab(abIn, 2)));
adjacent(5) = 1;

instead of the first set of if else statements, and make the name changes for the second set as well. The idea here is to only get the values of lattice where they are inside the borders. The four logical expressions are here used for logical indexing in the fourth line.