New scripts to determine phase transition boundary between modulated and unmodulated state

Dipolar-Gas-Simulator/+Scripts/determineDensityModulation.m

@@ -0,0 +1,47 @@
+function [ModulationFlag] = determineDensityModulation(psi, Params, Transf)
+    % Axes scaling and coordinates in micrometers
+    x               = Transf.x * Params.l0 * 1e6; 
+    dx              = x(2)-x(1); 
+    % Compute probability density |psi|^2
+    n               = abs(psi).^2;
+    nyz             = squeeze(trapz(n*dx,1));
+    
+    densityProfile  = sum(nyz,2);
+    
+    % We need to first smoothen the density profile and subtract the smoothened profile from the
+    % original density profile to - 
+    % 1. Remove the dominant low-frequency content.
+    % 2. Isolate the high-frequency components (like a sinusoidal modulation).
+    % FFT of the residual then highlights the periodic part, making the
+    % modulation easy to detect.
+
+    % Step 1: Smooth the density profile (Gaussian smoothing)
+    smoothedProfile = smooth(densityProfile, 10);
+    
+    % Step 2: Compute the residual (original - smoothed)
+    residual        = densityProfile - smoothedProfile; % We do this 
+    
+    % Step 3: Compute the Fourier Transform of the residual
+    N               = length(residual);
+    Y               = fft(residual);
+    P2              = abs(Y/N);       % Two-sided spectrum
+    P1              = P2(1:N/2+1);    % Single-sided spectrum
+    P1(2:end-1)     = 2*P1(2:end-1);  % Correct for the energy in the negative frequencies
+    
+    % Step 4: Check for significant peaks in the Fourier spectrum
+    % We check if the peak frequency is above a certain threshold
+    threshold       = 1E-3;           % This can be adjusted based on the expected modulation strength
+    peakValue       = max(P1);
+    
+    if peakValue > threshold
+        ModulationFlag = true;   % Indicates sinusoidal modulation
+    else
+        ModulationFlag = false;  % Indicates otherwise
+    end
+end
+
+
+
+    
+    
+            

Dipolar-Gas-Simulator/+Scripts/run_hybrid_worker_for_phase_boundary.m

@@ -0,0 +1,115 @@
+function run_hybrid_worker_for_phase_boundary(batchParams, batchIdx)
+    % Set up local cluster for parallel pool
+    cluster                    = parcluster('local');
+    nprocs                     = str2double(getenv('SLURM_CPUS_PER_TASK'));
+    if isnan(nprocs), nprocs   = feature('numcores'); end
+
+    tmpdir                     = fullfile(getenv('TMPDIR'), sprintf('matlab_job_%d', batchIdx));
+    if ~exist(tmpdir, 'dir'); mkdir(tmpdir); end
+    cluster.JobStorageLocation = tmpdir;
+
+    pool = parpool(cluster, nprocs);
+
+    nJobs = size(batchParams, 1);
+    parfor k = 1:nJobs
+        % Unpack parameter tuple
+        initial_a_s  = batchParams(k, 1);
+        theta_deg    = batchParams(k, 2);
+        phi_deg      = batchParams(k, 3);
+        N_atoms      = batchParams(k, 4);
+    
+        theta_rad    = deg2rad(theta_deg);
+        phi_rad      = deg2rad(phi_deg);
+    
+        adjusted_a_s = initial_a_s;
+    
+        % Create unique save directory
+        jobName = sprintf('aS_%03d_theta_%03d_phi_%03d_N_%d', initial_a_s, theta_deg, phi_deg, N_atoms);
+        saveDir = fullfile('./Results/Data_3D/GradientDescent', jobName);
+        if ~exist(saveDir, 'dir')
+            mkdir(saveDir);
+        end
+        
+        srcFile = './Results/Data_3D/GradientDescent/psi_init.mat';
+        if exist(srcFile, 'file')
+            copyfile(srcFile, fullfile(saveDir, 'psi_init.mat'));
+        end
+    
+        MaxAttempts  = 10;
+        SuccessFlag  = false;
+        AttemptCount = 0;
+    
+        while ~SuccessFlag && AttemptCount < MaxAttempts
+            AttemptCount = AttemptCount + 1;
+    
+            % Options for this run
+            OptionsStruct = struct;
+            OptionsStruct.NumberOfAtoms           = N_atoms;
+            OptionsStruct.DipolarPolarAngle       = theta_rad;
+            OptionsStruct.DipolarAzimuthAngle     = phi_rad;
+            OptionsStruct.ScatteringLength        = adjusted_a_s;
+    
+            OptionsStruct.TrapFrequencies         = [50, 20, 150];
+            OptionsStruct.TrapPotentialType       = 'Harmonic';
+    
+            OptionsStruct.NumberOfGridPoints      = [128, 256, 128];
+            OptionsStruct.Dimensions              = [30, 50, 30];
+            OptionsStruct.UseApproximationForLHY  = true;
+            OptionsStruct.IncludeDDICutOff        = true;
+            OptionsStruct.CutoffType              = 'CustomCylindrical';
+            OptionsStruct.CustomCylindricalCutOffRadius = 12;
+            OptionsStruct.CustomCylindricalCutOffHeight = 10;
+            OptionsStruct.SimulationMode          = 'ImaginaryTimeEvolution';
+            OptionsStruct.TimeStepSize            = 5E-4;
+            OptionsStruct.MinimumTimeStepSize     = 2E-6;
+            OptionsStruct.TimeCutOff              = 1E5;
+            OptionsStruct.EnergyTolerance         = 5E-10;
+            OptionsStruct.ResidualTolerance       = 1E-05;
+            OptionsStruct.NoiseScaleFactor        = 0.010;
+    
+            OptionsStruct.PlotLive                = false;
+            OptionsStruct.JobNumber               = 0;
+            OptionsStruct.RunOnGPU                = true;
+            OptionsStruct.SaveData                = true;
+            OptionsStruct.SaveDirectory           = saveDir;
+    
+            options = Helper.convertstruct2cell(OptionsStruct);
+    
+            sim = Simulator.DipolarGas(options{:});
+            pot = Simulator.Potentials(options{:});
+            sim.Potential = pot.trap();
+    
+            NumberOfOutputs = 5;
+            try
+                [Params, Transf, psi, ~, ~, stats] = Helper.runWithProfiling(@() sim.run(), NumberOfOutputs, saveDir);
+    
+                if Scripts.determineDensityModulation(psi, Params, Transf)
+                    SuccessFlag = true;
+                    runDir  = fullfile(saveDir, sprintf('Run_%03d', OptionsStruct.JobNumber));
+                    psiFile = fullfile(runDir, 'psi_gs.mat');
+    
+                    if exist(psiFile, 'file')
+                        Scripts.saveUpdatedScatteringLength(psiFile, adjusted_a_s, AttemptCount);
+                    else
+                        warning('Expected file %s not found. Cannot save final a_s.', psiFile);
+                    end
+                else
+                    adjusted_a_s = adjusted_a_s - 0.2;  % Tweak as per your needs
+                end
+    
+            catch ME
+                fprintf('ERROR in job %d (attempt %d):\n%s\n', k, AttemptCount, getReport(ME, 'extended'));
+                adjusted_a_s = adjusted_a_s - 0.2;
+            end
+        end
+    
+        if SuccessFlag
+            fprintf('Batch %d | Job %d: a_s = %.2f, theta = %d°, phi = %d°, N = %d | Time = %.2f s\n', ...
+                batchIdx, k, adjusted_a_s, theta_deg, phi_deg, N_atoms, stats.runtime);
+        else
+            fprintf('Batch %d | Job %d FAILED after %d tries.\n', batchIdx, k, MaxAttempts);
+        end
+    end
+
+    delete(pool);
+end

Dipolar-Gas-Simulator/+Scripts/run_hybrid_worker_for_phase_boundary_wrapper.m

@@ -0,0 +1,32 @@
+function run_hybrid_worker_for_phase_boundary_wrapper(batchIdx)
+    a_s_list     = parse_environmental_variable('SCATTERING_LENGTH_RANGE', 85);  % Scattering length list
+    theta        = parse_environmental_variable('POLAR_ANGLE_RANGE', 0);         % Single polar angle
+    phi          = parse_environmental_variable('AZIMUTHAL_ANGLE_RANGE', 0);     % Single azimuthal angle
+    N_atoms_list = parse_environmental_variable('NUM_ATOMS_LIST', 90000);        % Atom number list
+    chunkSize    = str2double(getenv('CHUNK_SIZE'));
+
+    % Validate matching list lengths
+    if numel(a_s_list) ~= numel(N_atoms_list)
+        error('Length of a_s_list and N_atoms_list must match.');
+    end
+
+    totalJobs = numel(a_s_list);
+    totalBatches = ceil(totalJobs / chunkSize);
+
+    if batchIdx > totalBatches
+        error('Batch index %d exceeds total batches (%d)', batchIdx, totalBatches);
+    end
+
+    firstIdx = (batchIdx - 1) * chunkSize + 1;
+    lastIdx  = min(batchIdx * chunkSize, totalJobs);
+
+    % Construct the batch parameters
+    batchParams = zeros(lastIdx - firstIdx + 1, 4); % Columns: a_s, theta, phi, N_atoms
+    for i = 1:size(batchParams, 1)
+        idx = firstIdx + i - 1;
+        batchParams(i, :) = [a_s_list(idx), theta, phi, N_atoms_list(idx)];
+    end
+
+    % Call worker with this batch of parameters
+    Scripts.run_hybrid_worker_for_phase_boundary(batchParams, batchIdx);
+end

Dipolar-Gas-Simulator/+Scripts/run_hybrid_worker_wrapper.m

@@ -1,5 +1,5 @@
 function run_hybrid_worker_wrapper(batchIdx)
-    a_s_list     = parse_environmental_variable('SCATTERING_LENGTH_RANGE', 85);         % Scattering length
+    a_s_list     = parse_environmental_variable('SCATTERING_LENGTH_RANGE', 85);  % Scattering length
     theta_list   = parse_environmental_variable('POLAR_ANGLE_RANGE', 0);         % Polar angle
     phi_list     = parse_environmental_variable('AZIMUTHAL_ANGLE_RANGE', 0);     % Azimuthal angle
     N_atoms_list = parse_environmental_variable('NUM_ATOMS_LIST', 90000);        % Atom number
@@ -20,7 +20,7 @@ function run_hybrid_worker_wrapper(batchIdx)
 
     % Call worker with this batch of parameters
     batchParams = paramGrid(firstIdx:lastIdx, :);
-    Scripts.run_hybrid_worker(batchParams, batchIdx);
+    Scripts.run_hybrid_worker_for_phase_boundary(batchParams, batchIdx);
 end
 
 function vals = parse_environmental_variable(varName, default)

Dipolar-Gas-Simulator/+Scripts/run_locally.m

@@ -549,13 +549,13 @@ SaveDirectory = './Results/Data_3D/AnisotropicTrap/TiltedDipoles45';
 JobNumber     = 0;
 Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
 %%
-SaveDirectory = './Results/Data_3D/GradientDescent/Phi050/aS_089_theta_050_phi_000_N_100000';
+SaveDirectory = './Results/Data_3D/GradientDescent/Phi020/aS_090_theta_020_phi_000_N_100000';
 JobNumber     = 0;
 Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
 %%
 % Parameters you can set before the loop
-N         = 165000;
+N         = 100000;
-theta     = 30.0;
+theta     = 20.0;
 phi       = 0;
 JobNumber = 0;
 
@@ -592,7 +592,7 @@ JobNumber     = 1;
 Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
 
 %% Visualize phase diagram
-load('phase_diagram_matrix_theta_30.mat')
+load('./Results/Data_3D/GradientDescent/phase_diagram_matrix_theta_30.mat')
 PhaseDiagramMatrix   = M;
 ScatteringLengths    = SCATTERING_LENGTH_RANGE;
 NumberOfAtoms        = NUM_ATOMS_LIST * 1E-5;
@@ -616,3 +616,77 @@ xlabel('Number of Atoms (x 10^5)', 'Interpreter', 'tex', FontSize=16);
 ylabel('Scattering Length (a0)', FontSize=16);
 % title('Zero-temperature Phase Diagram for \theta = 0', 'Interpreter', 'tex', FontSize=16);
 grid on;
+
+%% Density modulation determination
+
+SaveDirectory = './Results/Data_3D/GradientDescent/Phi030/aS_079_theta_030_phi_000_N_100000';
+JobNumber     = 0;
+
+% Load data
+Data          = load(sprintf(horzcat(SaveDirectory, '/Run_%03i/psi_gs.mat'),JobNumber),'psi','Params','Transf','Observ');
+
+Params        = Data.Params;
+Transf        = Data.Transf; 
+Observ        = Data.Observ; 
+
+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
+
+ModulationFlag = determineDensityModulation(psi, Params, Transf);
+
+if ModulationFlag
+    disp('The state is modulated');
+else
+    disp('The state is not modulated');
+end
+
+function ModulationFlag = determineDensityModulation(psi, Params, Transf)
+    
+    % Axes scaling and coordinates in micrometers
+    x               = Transf.x * Params.l0 * 1e6; 
+    dx              = x(2)-x(1); 
+    % Compute probability density |psi|^2
+    n               = abs(psi).^2;
+    nyz             = squeeze(trapz(n*dx,1));
+    
+    densityProfile  = sum(nyz,2);
+    
+    % We need to first smoothen the density profile and subtract the smoothened profile from the
+    % original density profile to - 
+    % 1. Remove the dominant low-frequency content.
+    % 2. Isolate the high-frequency components (like a sinusoidal modulation).
+    % FFT of the residual then highlights the periodic part, making the
+    % modulation easy to detect.
+
+    % Step 1: Smooth the density profile (Gaussian smoothing)
+    smoothedProfile = smooth(densityProfile, 10);
+    
+    % Step 2: Compute the residual (original - smoothed)
+    residual        = densityProfile - smoothedProfile; % We do this 
+    
+    % Step 3: Compute the Fourier Transform of the residual
+    N               = length(residual);
+    Y               = fft(residual);
+    P2              = abs(Y/N);       % Two-sided spectrum
+    P1              = P2(1:N/2+1);    % Single-sided spectrum
+    P1(2:end-1)     = 2*P1(2:end-1);  % Correct for the energy in the negative frequencies
+    
+    % Step 4: Check for significant peaks in the Fourier spectrum
+    % We check if the peak frequency is above a certain threshold
+    threshold       = 1E-3;           % This can be adjusted based on the expected modulation strength
+    peakValue       = max(P1);
+    
+    if peakValue > threshold
+        ModulationFlag = true;   % Indicates sinusoidal modulation
+    else
+        ModulationFlag = false;  % Indicates otherwise
+    end
+end

Dipolar-Gas-Simulator/+Scripts/run_on_cluster.m

@@ -13,21 +13,26 @@ OptionsStruct.Dimensions              = [30, 50, 30];
 OptionsStruct.UseApproximationForLHY  = true;
 OptionsStruct.IncludeDDICutOff        = true;
 OptionsStruct.CutoffType              = 'Cylindrical';
-OptionsStruct.SimulationMode          = 'EnergyMinimization';     % 'ImaginaryTimeEvolution' | 'RealTimeEvolution' | 'EnergyMinimization'
+OptionsStruct.SimulationMode          = 'ImaginaryTimeEvolution'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution' | 'EnergyMinimization'
 OptionsStruct.GradientDescentMethod   = 'NonLinearCGD';           % 'HeavyBall' | 'NonLinearCGD' 
 OptionsStruct.MaxIterationsForGD      = 15000;
+OptionsStruct.TimeStepSize            = 5E-4;            % in s
+OptionsStruct.MinimumTimeStepSize     = 2E-6;            % in s
+OptionsStruct.TimeCutOff              = 1E6;             % in s
+OptionsStruct.EnergyTolerance         = 5E-10;
+OptionsStruct.ResidualTolerance       = 1E-05;
 OptionsStruct.NoiseScaleFactor        = 0.010;
 
 OptionsStruct.PlotLive                = false;
 OptionsStruct.JobNumber               = 0;
 OptionsStruct.RunOnGPU                = true;
 OptionsStruct.SaveData                = true;
-OptionsStruct.SaveDirectory           = './Results/Data_3D/GradientDescent'; % './Results/Data_3D/AnisotropicTrap/Tilted0'
+OptionsStruct.SaveDirectory           = './Results/Data_3D/GradientDescent';
 options                               = Helper.convertstruct2cell(OptionsStruct);
 
 sim                                   = Simulator.DipolarGas(options{:});
 pot                                   = Simulator.Potentials(options{:});
-sim.Potential                         = pot.trap(); % + pot.repulsive_chopstick();
+sim.Potential                         = pot.trap();
 
 %-% Run Simulation %-%
 NumberOfOutputs                       = 5;

Dipolar-Gas-Simulator/+Scripts/saveUpdatedScatteringLength.m

@@ -0,0 +1,3 @@
+function saveUpdatedScatteringLength(filename, final_a_s, num_attempts)
+    save(filename, 'final_a_s', 'num_attempts');
+end

Dipolar-Gas-Simulator/submit_jobs.sh

@@ -2,9 +2,9 @@
 
 # Use space-separated floating-point/integer values
 SCATTERING_LENGTH_RANGE="[79.0 80.0 81.0 82.0 83.0 84.0 85.0 86.0 87.0 88.0 89.0 90.0]"
-POLAR_ANGLE_RANGE="[20.0 30.0 40.0 50.0]"
+POLAR_ANGLE_RANGE="[0.0 20.0 30.0]"
 AZIMUTHAL_ANGLE_RANGE="[0.0]"
-NUM_ATOMS_LIST="[50000 55000 60000 65000 70000 75000 80000 85000 90000 95000 100000 105000]"
+NUM_ATOMS_LIST="[50000 55000 60000 65000 70000 75000 80000 85000 90000 95000 100000 105000 110000 115000 120000 125000 130000 135000 140000 145000 150000 155000 160000 165000]"
 CHUNK_SIZE=4
 
 # ----------- Count total combinations for SLURM array -----------