Workaround for initial high memory usage?

I'm currently working on a master thesis in statistical mechanics where I'm looking at the 3D ising model. This uses the Metropolis Monte Carlo methods, and I generate quite alot of data. My problem is that when computing on my university's cluster, the initial memory usage, which can be of order 2GB (6%), really dampens how many cores I can use without blocking others out from the cluster. But after some time, the memory usage drops down to 1%. Is there a way to not spike the memory usage in the beginning of the script? Maybe something like running a small system in order to have the routines set up before embarking on a large system?

Cheers!

EDIT: Been asked to submit code, so it is found below. I have not shown every function, but below the script are the functions that handles most of the data.

clear variables
close all
%Constants (boltzmann constant set to k = 1)
%Following sets are as follows: [Jz,Jz',J1,J2], negative signs = AFM
JKobayashi = [0.6015,0,-0.0508,-0.038];
JCabrera1 = [2.19,-0.29,-0.031,0.77*0.031];
JCabrera2 = [2.4,-0.36,-0.036,0.75*0.036];
Jsets = [JKobayashi/JKobayashi(1);JCabrera1/JCabrera1(1);JCabrera2/JCabrera2(1)];
Temp = [0.1,0.15,linspace(0.2,0.8,100),linspace(0.85,2.0,10)]; %Temperature measured in units of J [K];
Sweeps = 400000; %Number of sweeps done for measurements
NBoot = 100; %Number of means calculated in error analysis
EquilibriumTime = 25000; %Twice the number of sweeps before equilibrium
UseSet = 2;

Constants = Jsets(UseSet,:);
LatticeSize = 16; %Specifying the lattice as a LatticeSize X LatticeSize X    LatticeSize

%colors = ["red";"blue";"green";"black";"cyan";"magenta"];

States = zeros(LatticeSize,LatticeSize,LatticeSize,length(Temp));
[E,E2,M,M2,Mabs,M4,HeatCap,MagSus,...
  SEHeatCap,SEMagSus,SEMag,SEMagAbs,Corrtimes]=...
  deal(zeros(length(Temp),1));

parfor i=1:length(Temp)
  %Setting up a random state
  rng(123)
  IsingState = (rand(LatticeSize,LatticeSize,LatticeSize) > 0.5)*2-1;

    %Equilibration of the random state
    IsingState = ...
        Equilibration(IsingState,Constants,Temp(i),EquilibriumTime*LatticeSize^3);

    %Metropolis
    [MagnetizationArray,EnergyArray,IsingState]=...
        Metropolis(IsingState,Constants,Temp(i),Sweeps);
    States(:,:,:,i) = IsingState;

    %Calculating the correlation time and gathering uncorrelated
    %measurements
     [EnergyArray,MagnetizationArray,Corrtimes(i)] = ...
         Measurements(EnergyArray,MagnetizationArray,EquilibriumTime);

    %Calculating means and standard error in the means
    [ E(i),E2(i),M(i),M2(i),Mabs(i),M4(i),HeatCap(i),...
    MagSus(i),SEHeatCap(i),SEMagSus(i),SEMag(i),SEMagAbs(i) ] =...
    Observables(EnergyArray,MagnetizationArray,Temp(i),NBoot,LatticeSize);

end

titles = ...
  {sprintf('Kobayashi %gx%gx%g %g',LatticeSize,LatticeSize,LatticeSize,Sweeps)...
  ,sprintf('Cabrera et al 1 %gx%gx%g %g',LatticeSize,LatticeSize,LatticeSize,Sweeps)...
  ,sprintf('Cabrera et al 2 %gx%gx%g %g',LatticeSize,LatticeSize,LatticeSize,Sweeps)};
figure(1)
subplot(2,2,1)
errorbar(Temp,HeatCap,SEHeatCap,'O')
xlabel('kT/J')
ylabel('C_v(T)')
title('Heat Capacity')
subplot(2,2,2)
errorbar(Temp,MagSus,SEMagSus,'O')
title('Magnetic Susceptibility')
xlabel('kT/J')
ylabel('chi(T)')
subplot(2,2,3)
errorbar(Temp,M,SEMag,'O')
title('Magnetization')
xlabel('kT/J')
ylabel('Magnetization')
subplot(2,2,4)
errorbar(Temp,Mabs,SEMagAbs,'O')
title('|Magnetization|')
xlabel('kT/J')
ylabel('|Magnetization|')
suptitle(char(titles(UseSet)))
saveas(figure(1),char(titles(UseSet)))
save(char(titles(UseSet)))
% figure()
% LatticePlot(States(:,:,:,15))
% view(0,87)

Metropolis, Observables, Measurements:

function [ MagnetizationArray,EnergyArray,IsingState ]...
  = Metropolis( IsingState,J,Temp,Sweeps )
%UNTITLED2 Summary of this function goes here
%   Detailed explanation goes here
Emptymatrix = zeros(10,1);
LatticeSize = length(IsingState);
Magnetization = sum(sum(sum(IsingState)));
Energy = -sum(sum(sum(IsingState.*...
  (J(4)*(circshift(IsingState,[1,0,0])...
  +circshift(IsingState,[-1,0,0])...
  +circshift(IsingState,[-1,-1,0])...
  +circshift(IsingState,[1,1,0]))...
  +J(3)*(circshift(IsingState,[0,1,0])...
  +circshift(IsingState,[0,-1,0]))...
  +J(1)*(circshift(IsingState,[0,0,-1])...
  +circshift(IsingState,[0,0,1]))...
  +J(2)*(circshift(IsingState,[0,0,-2])...
  +circshift(IsingState,[0,0,2]))...
  ))))/2;
State = IsingState;
EnergyArray = zeros(Sweeps,1);
MagnetizationArray = zeros(Sweeps,1);
for l = 1:Sweeps
  spinsites = randi(LatticeSize^3,LatticeSize^3,1);
  randomprob = rand(LatticeSize^3,1);
  for j =1:LatticeSize^3
      DeltaE = EnergyDifference( IsingState,J,spinsites(j),Emptymatrix,LatticeSize);
      if DeltaE <= 0 || exp(-DeltaE/Temp) >= randomprob(j)
          IsingState(spinsites(j)) = -IsingState(spinsites(j));
          Energy = Energy + DeltaE;
          Magnetization = Magnetization...
                  +2*IsingState(spinsites(j));
      end
      State = State + IsingState;
  end
  EnergyArray(l) = Energy;
  MagnetizationArray(l)= Magnetization;
end
IsingState = (State/(Sweeps*LatticeSize^3+1)>zeros(LatticeSize,LatticeSize,LatticeSize))...
  + (State/(Sweeps*LatticeSize^3+1)<zeros(LatticeSize,LatticeSize,LatticeSize))*(-1);
%MagnetizationArray = abs(MagnetizationArray)/LatticeSize^3;
end

function [ E,E2,M,M2,Mabs,M4,HeatCap,...
  MagSus,SEHeatCap,SEMagSus,SEMag,SEMagAbs ] =...
  Observables( Energy,Magnetization,Temp,NBoot,LatticeSize)
%UNTITLED4 Summary of this function goes here
%   Detailed explanation goes here
E = mean(Energy); 
E2 = mean(Energy.^2);
M = mean(Magnetization);
M2 = mean(Magnetization.^2);
Mabs = mean(abs(Magnetization)); 
M4 = mean(Magnetization.^4);

HeatCap = 1/(LatticeSize^3*Temp^2)*(E2-E^2);
MagSus = 1/(LatticeSize^3*Temp)*(M2-M^2);
SEHeatCap = ErrorHeatCap(Energy,NBoot,LatticeSize,Temp);
SEMagSus = ErrorMagSus(Magnetization,NBoot,LatticeSize,Temp);

M = M/LatticeSize^3;
Mabs = Mabs/LatticeSize^3;
SEMag = ErrorMag(Magnetization,NBoot,LatticeSize);
SEMagAbs = ErrorMag(abs(Magnetization),NBoot,LatticeSize);

end

function [ SEHeatCap ] = ErrorHeatCap(Energy,Nboot,LatticeSize,Temp)
%UNTITLED4 Summary of this function goes here
%   Detailed explanation goes here
n = length(Energy);
HeatCapTemp = zeros(Nboot,1);
for i = 1:Nboot
  j = randi(n,1,n);
  Etemp = Energy(j);
  E = mean(Etemp);
  E2 = mean(Etemp.^2);
  HeatCapTemp(i) = 1/(LatticeSize^3*Temp^2)*(E2-E^2);
end
SEHeatCap = sqrt(mean(HeatCapTemp.^2)-mean(HeatCapTemp)^2);
end

function [SEMagsus] = ErrorMagSus(Magnetization,Nboot,LatticeSize,Temp)
%UNTITLED4 Summary of this function goes here
%   Detailed explanation goes here
n = length(Magnetization);
MagSusTemp = zeros(Nboot,1);
for i = 1:Nboot
  j = randi(n,1,n);
  Mtemp = Magnetization(j);
  M = mean(Mtemp);
  M2 = mean(Mtemp.^2);
  MagSusTemp(i) = 1/(LatticeSize^3*Temp)*(M2-M^2);
end
SEMagsus = sqrt(mean(MagSusTemp.^2)-mean(MagSusTemp)^2);
end

function [ SEMag ] = ErrorMag( Magnetization,Nboot,LatticeSize )
%UNTITLED4 Summary of this function goes here
%   Detailed explanation goes here
n = length(Magnetization);
Magtemp = zeros(Nboot,1);
for i = 1:Nboot
  j = randi(n,1,n);
  Magtemp(i) = mean(Magnetization(j))/LatticeSize^3;
end
SEMag = sqrt(mean(Magtemp.^2)-mean(Magtemp)^2);
end

function [ EnergyArray,MagnetizationArray,CorrTime ] = Measurements( Energy,Magnetization,EquilibriumTime )
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here
Sweeps = length(Magnetization);
CorrFunc = autocorr(Magnetization,Sweeps-1);
zerocorr = find(CorrFunc<0,1);

if isempty(zerocorr)==1
  CorrTime = 1;
elseif zerocorr >= EquilibriumTime/2
  CorrTime = 2*ceil(trapz(...
  CorrFunc(1:EquilibriumTime/2)));
else
  CorrTime = 2*ceil(trapz(...
  CorrFunc(1:zerocorr)));
end

if isnan(CorrTime)==1
  CorrTime = 1;
end
EnergyArray = Energy(1:CorrTime:Sweeps);
MagnetizationArray = Magnetization(1:CorrTime:Sweeps);

 end

[Jan]

Without seeing the code and knowing any details of what you are doing, suggesting a modification will be based on pure guessing. Please edit the question and provide some useful details.

[Carl]

Hi Aleksander, thanks for providing the code. It would be helpful if you could pinpoint exactly which lines are causing the high initial memory usage. I would recommend tracking your system's memory as you step through your code line-by-line (unfortunately, there is no supported way of doing memory profiling in MATLAB). You can step through your code by setting breakpoints and using the debugger.

You can use something like the 'memory' command in MATLAB, an analogous system command, your OS task manager, or some other tool to track memory usage. After doing this, can you see which lines are causing the high memory use?

[Aleksander Seland]

Thanks! I found out there was an array being created that was unnecessarily long.