diff --git a/Dipolar-Gas-Simulator/+Scripts/determineDensityModulation.m b/Dipolar-Gas-Simulator/+Scripts/determineDensityModulation.m new file mode 100644 index 0000000..7baa605 --- /dev/null +++ b/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 + + + + + + \ No newline at end of file diff --git a/Dipolar-Gas-Simulator/+Scripts/run_hybrid_worker_for_phase_boundary.m b/Dipolar-Gas-Simulator/+Scripts/run_hybrid_worker_for_phase_boundary.m new file mode 100644 index 0000000..15b7c2f --- /dev/null +++ b/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 \ No newline at end of file diff --git a/Dipolar-Gas-Simulator/+Scripts/run_hybrid_worker_for_phase_boundary_wrapper.m b/Dipolar-Gas-Simulator/+Scripts/run_hybrid_worker_for_phase_boundary_wrapper.m new file mode 100644 index 0000000..cb64dfe --- /dev/null +++ b/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 diff --git a/Dipolar-Gas-Simulator/+Scripts/run_hybrid_worker_wrapper.m b/Dipolar-Gas-Simulator/+Scripts/run_hybrid_worker_wrapper.m index fa41e35..d920a38 100644 --- a/Dipolar-Gas-Simulator/+Scripts/run_hybrid_worker_wrapper.m +++ b/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) diff --git a/Dipolar-Gas-Simulator/+Scripts/run_locally.m b/Dipolar-Gas-Simulator/+Scripts/run_locally.m index c26fdc8..debd816 100644 --- a/Dipolar-Gas-Simulator/+Scripts/run_locally.m +++ b/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; -theta = 30.0; +N = 100000; +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; @@ -615,4 +615,78 @@ yticklabels(string(ScatteringLengths)); 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; \ No newline at end of file +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 \ 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 4927b21..dc4a3de 100644 --- a/Dipolar-Gas-Simulator/+Scripts/run_on_cluster.m +++ b/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; diff --git a/Dipolar-Gas-Simulator/+Scripts/saveUpdatedScatteringLength.m b/Dipolar-Gas-Simulator/+Scripts/saveUpdatedScatteringLength.m new file mode 100644 index 0000000..67a3b76 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/Dipolar-Gas-Simulator/submit_jobs.sh b/Dipolar-Gas-Simulator/submit_jobs.sh index 7cbe973..431df34 100644 --- a/Dipolar-Gas-Simulator/submit_jobs.sh +++ b/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 -----------