Monte Carlo simulation for two dimensional particles in the NVT ensemble-Lennard Jones interaction

A Monte Carlo simulation for two dimensional particles in the NVT ensemble-Lennard Jones interaction

Monte-Carlo in NVT ensemble for Lennard-Jonse potantioal in 2D
calculate the energy after 'Nstep' Monte-Carlo runs (Metropolis algorithm)
for a system of 'N' particles in 2D interacting via Lennard-Jonse pair
potential. periodic boundary conditions (PBC) are apllied.

inputs:
N - number of particles
T - reduced temperature
Nsteps - number of steps
maxdr - maximum particle displacement
initialConfig - initial configuration of particles (2 by N matrix)
initialU - initial energy of the configuration
rCutoff - the cutoff distance for the energy
optional: the 'm' in the interaction (U = 4((1/r)^12 - (1/r)^-m)),
default is 6.

outputs:
finalU - the energy in each step (1 by Nstep matrix)
finalConfigurations - the coordinates of all particles in each step (2 by N
by Nstep matrix)
finalDistances - the pair distances of all particles in each step (N by N
by Nstep matrix)
moveCount - counts accepted moves

the potential of Lennard-Jonse in reduced units:
U = 4*[(1/r)^12 - (1/r)^6]
(you can change the 6 with the optional input 'm')

the potential of Lennard-Jonse in non-reduced units:
U = 4*epsilon*[(sigma/r)^12 - (sigma/r)^6]

reduced units:
T(reduced) = kT/epsilon | r(reduced) = r/sigma | U(reduced) = U/epsilon

usage examples:

create an initial configuration of 100 particles in a 100*100 box:
N = 100;
L = 100;
initialConfig = L*rand(2,N);
rho = N/L^2;

choose simulation parameters:
T = 1; Nsteps = 1000; maxdr = 1; rCutoff = 2.5;

calculate the pair distances:
initialDistances = sqrt(bsxfun(@(x1,x2) (x1-x2).^2 ,...
initialConfig(1,:),initialConfig(1,:)')...
+bsxfun(@(x3,x4) (x4-x3).^2 ,...
initialConfig(2,:),initialConfig(2,:)'));

calculate the initial energy (up to the cutoff):
d = initialDistances(and(initialDistances < rCutoff,initialDistances > 0));
initialU = 4*sum(sum(d.^(-12)-d.^(-6)));

calculate energy
[finalU,~,~,finalConfiguration,finalDistances,moveCount] = ...
MonteCarlo2DLJHeart(N,T,rho,Nsteps,maxdr,initialConfig,rCutoff...
,initialDistances,initialU)

also calculate the pressure.
first calculate the initial virial:
virial = -(2*rho/N)*(sum(sum(6*d.^(-6) - 12*d.^(-12))));

now calculate the energy and pressure using monte carlo
[finalU,finalVirial,finalPressure,finalConfiguration,finalDistances,moveCount] = ...
MonteCarlo2DLJHeart(N,T,rho,Nsteps,maxdr,initialConfig,rCutoff...
,initialDistances,initialU,'virial',virial);