diff --git a/Dipolar-Gas-Simulator/+Plotter/visualizeGSWavefunction.m b/Dipolar-Gas-Simulator/+Plotter/visualizeGSWavefunction.m
index 99dd50e..0016a12 100644
--- a/Dipolar-Gas-Simulator/+Plotter/visualizeGSWavefunction.m
+++ b/Dipolar-Gas-Simulator/+Plotter/visualizeGSWavefunction.m
@@ -12,8 +12,16 @@ function visualizeGSWavefunction(folder_path, run_index)
     Transf = Data.Transf; 
     Observ = Data.Observ; 
     
-    if isgpuarray(Data.psi), psi = gather(Data.psi); end
-    if isgpuarray(Data.Observ.residual), Observ.residual = gather(Data.Observ.residual); end
+    if isgpuarray(Data.psi)
+        psi = gather(Data.psi);
+    else
+        psi = Data.psi;
+    end
+    if isgpuarray(Data.Observ.residual)
+        Observ.residual = gather(Data.Observ.residual); 
+    else
+        Observ.residual = Data.Observ.residual; 
+    end
     
     % Axes scaling and coordinates in micrometers
     x = Transf.x * Params.l0 * 1e6; 
diff --git a/Dipolar-Gas-Simulator/+Scripts/run_locally.m b/Dipolar-Gas-Simulator/+Scripts/run_locally.m
index 8f4e809..800e105 100644
--- a/Dipolar-Gas-Simulator/+Scripts/run_locally.m
+++ b/Dipolar-Gas-Simulator/+Scripts/run_locally.m
@@ -20,7 +20,7 @@ OptionsStruct.UseApproximationForLHY  = true;
 OptionsStruct.IncludeDDICutOff        = true;
 OptionsStruct.CutoffType              = 'Cylindrical';
 OptionsStruct.SimulationMode          = 'EnergyMinimization';     % 'ImaginaryTimeEvolution' | 'RealTimeEvolution' | 'EnergyMinimization'
-OptionsStruct.GradientDescentMethod   = 'HeavyBall';           % 'HeavyBall' | 'NonLinearCGD' 
+OptionsStruct.GradientDescentMethod   = 'NonLinearCGD';           % 'HeavyBall' | 'NonLinearCGD' 
 OptionsStruct.MaxIterationsForGD      = 1000;
 OptionsStruct.TimeStepSize            = 1E-3;            % in s
 OptionsStruct.MinimumTimeStepSize     = 1E-6;            % in s
@@ -547,4 +547,8 @@ Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
 %%
 SaveDirectory = './Results/Data_3D/AnisotropicTrap/TiltedDipoles45';
 JobNumber     = 0;
+Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
+%%
+SaveDirectory = './Results/Data_3D/GradientDescent';
+JobNumber     = 0;
 Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
\ No newline at end of file
diff --git a/Dipolar-Gas-Simulator/+Scripts/run_on_cluster.m b/Dipolar-Gas-Simulator/+Scripts/run_on_cluster.m
index 34f6162..50d6588 100644
--- a/Dipolar-Gas-Simulator/+Scripts/run_on_cluster.m
+++ b/Dipolar-Gas-Simulator/+Scripts/run_on_cluster.m
@@ -1,118 +1,40 @@
-theta = 15;
+%% - Labyrinth phase
 
 OptionsStruct = struct;
 
-OptionsStruct.NumberOfAtoms           = 1E5;
-OptionsStruct.DipolarPolarAngle       = deg2rad(theta);
-OptionsStruct.DipolarAzimuthAngle     = deg2rad(0);
-OptionsStruct.ScatteringLength        = 88;
+OptionsStruct.NumberOfAtoms           = 5E5;
+OptionsStruct.DipolarPolarAngle       = deg2rad(0);
+OptionsStruct.DipolarAzimuthAngle     = 0;
+OptionsStruct.ScatteringLength        = 85;
 
-OptionsStruct.TrapFrequencies         = [50, 20, 150];
+OptionsStruct.TrapFrequencies         = [125, 125, 350];
 OptionsStruct.TrapPotentialType       = 'Harmonic';
 
-OptionsStruct.NumberOfGridPoints      = [128, 256, 64];
-OptionsStruct.Dimensions              = [50, 50, 10];
-
+OptionsStruct.NumberOfGridPoints      = [128, 128, 64];
+OptionsStruct.Dimensions              = [18, 18, 18];
 OptionsStruct.UseApproximationForLHY  = true;
 OptionsStruct.IncludeDDICutOff        = true;
-OptionsStruct.CutoffType              = 'CustomCylindrical';
-OptionsStruct.CustomCylindricalCutOffRadius = 20.0;
-OptionsStruct.CustomCylindricalCutOffHeight = 4.0;
-
-OptionsStruct.SimulationMode          = 'ImaginaryTimeEvolution'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution'
-OptionsStruct.TimeStepSize            = 1E-3;            % in s
-OptionsStruct.MinimumTimeStepSize     = 2E-6;            % in s
-OptionsStruct.TimeCutOff              = 2E6;             % in s
-OptionsStruct.EnergyTolerance         = 5E-10;
-OptionsStruct.ResidualTolerance       = 1E-05;
-OptionsStruct.NoiseScaleFactor        = 0.01;
+OptionsStruct.CutoffType              = 'Cylindrical';
+OptionsStruct.SimulationMode          = 'EnergyMinimization';  % 'ImaginaryTimeEvolution' | 'RealTimeEvolution' | 'EnergyMinimization'
+OptionsStruct.GradientDescentMethod   = 'NonLinearCGD';        % 'HeavyBall' | 'NonLinearCGD' 
+OptionsStruct.MaxIterationsForGD      = 1E5;
+OptionsStruct.NoiseScaleFactor        = 0.010;
 
 OptionsStruct.PlotLive                = false;
 OptionsStruct.JobNumber               = 0;
 OptionsStruct.RunOnGPU                = true;
 OptionsStruct.SaveData                = true;
-OptionsStruct.SaveDirectory           = sprintf('./Results/Data_3D/TiltedDipoles%s', strrep(num2str(round(rad2deg(OptionsStruct.DipolarPolarAngle),2)), '.', '_'));
+OptionsStruct.SaveDirectory           = './Results/Data_3D/GradientDescent';
 options                               = Helper.convertstruct2cell(OptionsStruct);
-clear OptionsStruct
 
 sim                                   = Simulator.DipolarGas(options{:});
 pot                                   = Simulator.Potentials(options{:});
-sim.Potential                         = pot.trap(); 
+sim.Potential                         = pot.trap();
 
 %-% Run Simulation %-%
-[Params, Transf, psi, V, VDk]         = sim.run();
+NumberOfOutputs                       = 5;
+[Params, Transf, psi, V, VDk, stats]  = Helper.runWithProfiling(@() sim.run(), NumberOfOutputs, OptionsStruct.SaveDirectory);
+fprintf('Runtime: %.3f seconds\n', stats.runtime);
+fprintf('Memory used: %.2f MB\n', stats.workspaceMemoryMB);
 
-%{
-
-%% Imaginary-Time followed by Real-time
-
-% - Imaginary-Time
-OptionsStruct = struct;
-
-OptionsStruct.NumberOfAtoms                 = 1E5;
-OptionsStruct.DipolarPolarAngle             = deg2rad(0);
-OptionsStruct.DipolarAzimuthAngle           = deg2rad(0);
-OptionsStruct.ScatteringLength              = 88;
-
-OptionsStruct.TrapFrequencies               = [50, 20, 150];
-OptionsStruct.TrapPotentialType             = 'Harmonic';
-
-OptionsStruct.NumberOfGridPoints            = [128, 256, 64];
-OptionsStruct.Dimensions                    = [50, 50, 10];
-OptionsStruct.UseApproximationForLHY        = true;
-OptionsStruct.IncludeDDICutOff              = true;
-OptionsStruct.CutoffType                    = 'CustomCylindrical';
-OptionsStruct.CustomCylindricalCutOffRadius = 20.0;
-OptionsStruct.CustomCylindricalCutOffHeight = 4.0;
-
-OptionsStruct.SimulationMode                = 'ImaginaryTimeEvolution'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution' | 'EnergyMinimization'
-OptionsStruct.TimeStepSize                  = 5E-3;           % in s
-OptionsStruct.MinimumTimeStepSize           = 1E-6;           % in s
-OptionsStruct.TimeCutOff                    = 1E2;            % in s
-OptionsStruct.EnergyTolerance               = 5E-10;
-OptionsStruct.ResidualTolerance             = 1E-08;
-OptionsStruct.NoiseScaleFactor              = 0.05;
-OptionsStruct.PlotLive                      = false;
-OptionsStruct.JobNumber                     = 0;
-OptionsStruct.RunOnGPU                      = true;
-OptionsStruct.SaveData                      = false;
-OptionsStruct.SaveDirectory                 = './Results/Data_3D/RealTimeDynamics';
-options                                     = Helper.convertstruct2cell(OptionsStruct);
-
-sim                                         = Simulator.DipolarGas(options{:});
-pot                                         = Simulator.Potentials(options{:});
-sim.Potential                               = pot.trap();
-
-%-% Run Simulation %-%
-[Params, Transf, psi, V, VDk]               = sim.run();
-
-% Save only final state as initial wavefunction for real-time propagation
-mkdir(sprintf(OptionsStruct.SaveDirectory))
-save(sprintf(strcat(OptionsStruct.SaveDirectory, '/psi_init.mat'),Params.njob),'psi','Transf','Params','VDk','V');
-
-OptionsStruct.SimulationMode                = 'RealTimeEvolution'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution' | 'EnergyMinimization'
-OptionsStruct.EquilibrationTime             = 10E-3; 
-OptionsStruct.QuenchTime                    = 30E-3;
-OptionsStruct.HoldTime                      = 50E-3;
-OptionsStruct.QuenchScatteringLength        = false;
-OptionsStruct.RotateDipoles                 = true;
-OptionsStruct.FinalScatteringLength         = 88;
-OptionsStruct.FinalDipolarPolarAngle        = deg2rad(40);
-OptionsStruct.FinalDipolarAzimuthAngle      = deg2rad(0);
-OptionsStruct.TimeStepSize                  = 1E-3;            % in s
-OptionsStruct.NoiseScaleFactor              = 0.010;
-
-OptionsStruct.PlotLive                      = false;
-OptionsStruct.JobNumber                     = 1;
-OptionsStruct.RunOnGPU                      = true;
-OptionsStruct.SaveData                      = true;
-options                                     = Helper.convertstruct2cell(OptionsStruct);
-clear OptionsStruct
-
-sim                                         = Simulator.DipolarGas(options{:});
-pot                                         = Simulator.Potentials(options{:});
-sim.Potential                               = pot.trap(); 
-
-%-% Run Simulation %-%
-[Params, Transf, psi, V, VDk]               = sim.run();
-%}
\ No newline at end of file
+clear all
\ No newline at end of file