diff --git a/Dipolar-Gas-Simulator/+Plotter/visualizeGSWavefunction.m b/Dipolar-Gas-Simulator/+Plotter/visualizeGSWavefunction.m index 52b053b..1bb9ef4 100644 --- a/Dipolar-Gas-Simulator/+Plotter/visualizeGSWavefunction.m +++ b/Dipolar-Gas-Simulator/+Plotter/visualizeGSWavefunction.m @@ -5,7 +5,6 @@ function visualizeGSWavefunction(folder_path, run_index) format long - % Load data Data = load(sprintf(horzcat(folder_path, '/Run_%03i/psi_gs.mat'),run_index),'psi','Params','Transf','Observ'); diff --git a/Dipolar-Gas-Simulator/+Scripts/run_locally.m b/Dipolar-Gas-Simulator/+Scripts/run_locally.m index 11754b4..b499fcd 100644 --- a/Dipolar-Gas-Simulator/+Scripts/run_locally.m +++ b/Dipolar-Gas-Simulator/+Scripts/run_locally.m @@ -17,17 +17,17 @@ VerticalTrapFrequency = AspectRatio * HorizontalTrapFrequency; OptionsStruct.TrapFrequencies = [HorizontalTrapFrequency, HorizontalTrapFrequency, VerticalTrapFrequency]; OptionsStruct.TrapPotentialType = 'Harmonic'; -OptionsStruct.NumberOfGridPoints = [64, 64, 32]; +OptionsStruct.NumberOfGridPoints = [128, 128, 64]; OptionsStruct.Dimensions = [18, 18, 18]; OptionsStruct.UseApproximationForLHY = true; OptionsStruct.IncludeDDICutOff = true; OptionsStruct.CutoffType = 'Cylindrical'; -OptionsStruct.SimulationMode = 'ImaginaryTimeEvolution'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution' -OptionsStruct.TimeStepSize = 0.002; % in s -OptionsStruct.MinimumTimeStepSize = 1E-6; % in s +OptionsStruct.SimulationMode = 'EnergyMinimization'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution' | 'EnergyMinimization' +OptionsStruct.TimeStepSize = 1E-4; % in s +OptionsStruct.MinimumTimeStepSize = 2E-10; % in s OptionsStruct.TimeCutOff = 2E6; % in s OptionsStruct.EnergyTolerance = 5E-10; -OptionsStruct.ResidualTolerance = 1E-05; +OptionsStruct.ResidualTolerance = 1E-08; OptionsStruct.NoiseScaleFactor = 0.01; OptionsStruct.PlotLive = true; @@ -51,19 +51,34 @@ sim.Potential = pot.trap(); % + pot.repulsive_chopstick( Plotter.visualizeTrapPotential(sim.Potential,Params,Transf) %% - Plot initial wavefunction Plotter.visualizeWavefunction(psi,Params,Transf) +%% +SaveDirectory = './Results/Data_3D/AspectRatio2'; +JobNumber = 0; +Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber) %% - Plot GS wavefunction -% SaveDirectory = './Results/Data_3D/AspectRatio2_8'; -SaveDirectory = './Results/Data_3D/AspectRatio3_7'; -% SaveDirectory = './Results/Data_3D/AspectRatio3_8'; +% SaveDirectory = './Results/Data_3D/CompleteLHY/AspectRatio2_8'; +% SaveDirectory = './Results/Data_3D/CompleteLHY/AspectRatio3_7'; +SaveDirectory = './Results/Data_3D/CompleteLHY/AspectRatio3_8'; JobNumber = 0; Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber) %% -% SaveDirectory = './Results/Data_3D/BeyondSSD_SSD'; -SaveDirectory = './Results/Data_3D/BeyondSSD_Labyrinth'; -% SaveDirectory = './Results/Data_3D/BeyondSSD_Honeycomb'; +% SaveDirectory = './Results/Data_3D/CompleteLHY/BeyondSSD_SSD'; +% SaveDirectory = './Results/Data_3D/CompleteLHY/BeyondSSD_Labyrinth'; +SaveDirectory = './Results/Data_3D/CompleteLHY/BeyondSSD_Honeycomb'; +JobNumber = 0; +Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber) +%% - Plot GS wavefunction +% SaveDirectory = './Results/Data_3D/ApproximateLHY/AspectRatio2_8'; +% SaveDirectory = './Results/Data_3D/ApproximateLHY/AspectRatio3_7'; +SaveDirectory = './Results/Data_3D/ApproximateLHY/AspectRatio3_8'; +JobNumber = 0; +Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber) +%% +% SaveDirectory = './Results/Data_3D/ApproximateLHY/BeyondSSD_SSD'; +SaveDirectory = './Results/Data_3D/ApproximateLHY/BeyondSSD_Labyrinth'; +% SaveDirectory = './Results/Data_3D/ApproximateLHY/BeyondSSD_Honeycomb'; JobNumber = 0; Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber) - %% % To reproduce results from the Blair Blakie paper: diff --git a/Dipolar-Gas-Simulator/+Scripts/run_on_cluster.m b/Dipolar-Gas-Simulator/+Scripts/run_on_cluster.m index 79f3aa0..8487045 100644 --- a/Dipolar-Gas-Simulator/+Scripts/run_on_cluster.m +++ b/Dipolar-Gas-Simulator/+Scripts/run_on_cluster.m @@ -5,7 +5,7 @@ OptionsStruct = struct; OptionsStruct.NumberOfAtoms = 4E5; OptionsStruct.DipolarPolarAngle = deg2rad(0); OptionsStruct.DipolarAzimuthAngle = 0; -OptionsStruct.ScatteringLength = 87; +OptionsStruct.ScatteringLength = 85; AspectRatio = 2.0; HorizontalTrapFrequency = 125; @@ -13,21 +13,22 @@ VerticalTrapFrequency = AspectRatio * HorizontalTrapFrequency; OptionsStruct.TrapFrequencies = [HorizontalTrapFrequency, HorizontalTrapFrequency, VerticalTrapFrequency]; OptionsStruct.TrapPotentialType = 'Harmonic'; -OptionsStruct.NumberOfGridPoints = [256, 256, 128]; +OptionsStruct.NumberOfGridPoints = [128, 128, 64]; OptionsStruct.Dimensions = [18, 18, 18]; +OptionsStruct.UseApproximationForLHY = true; OptionsStruct.IncludeDDICutOff = true; OptionsStruct.CutoffType = 'Cylindrical'; OptionsStruct.SimulationMode = 'ImaginaryTimeEvolution'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution' -OptionsStruct.TimeStepSize = 2E-3; % in s +OptionsStruct.TimeStepSize = 1E-3; % in s OptionsStruct.MinimumTimeStepSize = 1E-6; % in s OptionsStruct.TimeCutOff = 2E6; % in s OptionsStruct.EnergyTolerance = 5E-10; OptionsStruct.ResidualTolerance = 1E-05; OptionsStruct.NoiseScaleFactor = 0.01; -OptionsStruct.PlotLive = false; +OptionsStruct.PlotLive = true; OptionsStruct.JobNumber = 0; -OptionsStruct.RunOnGPU = true; +OptionsStruct.RunOnGPU = false; OptionsStruct.SaveData = true; OptionsStruct.SaveDirectory = './Results/Data_3D/BeyondSSD_Labyrinth'; options = Helper.convertstruct2cell(OptionsStruct); diff --git a/Dipolar-Gas-Simulator/+Scripts/run_on_cluster_beyondSSD.m b/Dipolar-Gas-Simulator/+Scripts/run_on_cluster_beyondSSD.m index 3a31e3a..c5438c7 100644 --- a/Dipolar-Gas-Simulator/+Scripts/run_on_cluster_beyondSSD.m +++ b/Dipolar-Gas-Simulator/+Scripts/run_on_cluster_beyondSSD.m @@ -13,16 +13,17 @@ VerticalTrapFrequency = AspectRatio * HorizontalTrapFrequency; OptionsStruct.TrapFrequencies = [HorizontalTrapFrequency, HorizontalTrapFrequency, VerticalTrapFrequency]; OptionsStruct.TrapPotentialType = 'Harmonic'; -OptionsStruct.NumberOfGridPoints = [256, 256, 128]; +OptionsStruct.NumberOfGridPoints = [128, 128, 64]; OptionsStruct.Dimensions = [18, 18, 18]; +OptionsStruct.UseApproximationForLHY = true; OptionsStruct.IncludeDDICutOff = true; OptionsStruct.CutoffType = 'Cylindrical'; OptionsStruct.SimulationMode = 'ImaginaryTimeEvolution'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution' -OptionsStruct.TimeStepSize = 2E-3; % in s -OptionsStruct.MinimumTimeStepSize = 1E-6; % in s +OptionsStruct.TimeStepSize = 1E-4; % in s +OptionsStruct.MinimumTimeStepSize = 2E-10; % in s OptionsStruct.TimeCutOff = 2E6; % in s OptionsStruct.EnergyTolerance = 5E-10; -OptionsStruct.ResidualTolerance = 1E-05; +OptionsStruct.ResidualTolerance = 1E-08; OptionsStruct.NoiseScaleFactor = 0.01; OptionsStruct.PlotLive = false; @@ -55,16 +56,17 @@ VerticalTrapFrequency = AspectRatio * HorizontalTrapFrequency; OptionsStruct.TrapFrequencies = [HorizontalTrapFrequency, HorizontalTrapFrequency, VerticalTrapFrequency]; OptionsStruct.TrapPotentialType = 'Harmonic'; -OptionsStruct.NumberOfGridPoints = [256, 256, 128]; +OptionsStruct.NumberOfGridPoints = [128, 128, 64]; OptionsStruct.Dimensions = [18, 18, 18]; +OptionsStruct.UseApproximationForLHY = true; OptionsStruct.IncludeDDICutOff = true; OptionsStruct.CutoffType = 'Cylindrical'; OptionsStruct.SimulationMode = 'ImaginaryTimeEvolution'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution' -OptionsStruct.TimeStepSize = 2E-3; % in s -OptionsStruct.MinimumTimeStepSize = 1E-6; % in s +OptionsStruct.TimeStepSize = 1E-4; % in s +OptionsStruct.MinimumTimeStepSize = 2E-10; % in s OptionsStruct.TimeCutOff = 2E6; % in s OptionsStruct.EnergyTolerance = 5E-10; -OptionsStruct.ResidualTolerance = 1E-05; +OptionsStruct.ResidualTolerance = 1E-08; OptionsStruct.NoiseScaleFactor = 0.01; OptionsStruct.PlotLive = false; @@ -97,16 +99,17 @@ VerticalTrapFrequency = AspectRatio * HorizontalTrapFrequency; OptionsStruct.TrapFrequencies = [HorizontalTrapFrequency, HorizontalTrapFrequency, VerticalTrapFrequency]; OptionsStruct.TrapPotentialType = 'Harmonic'; -OptionsStruct.NumberOfGridPoints = [256, 256, 128]; +OptionsStruct.NumberOfGridPoints = [128, 128, 64]; OptionsStruct.Dimensions = [18, 18, 18]; +OptionsStruct.UseApproximationForLHY = true; OptionsStruct.IncludeDDICutOff = true; OptionsStruct.CutoffType = 'Cylindrical'; OptionsStruct.SimulationMode = 'ImaginaryTimeEvolution'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution' -OptionsStruct.TimeStepSize = 2E-3; % in s -OptionsStruct.MinimumTimeStepSize = 1E-6; % in s +OptionsStruct.TimeStepSize = 1E-4; % in s +OptionsStruct.MinimumTimeStepSize = 2E-10; % in s OptionsStruct.TimeCutOff = 2E6; % in s OptionsStruct.EnergyTolerance = 5E-10; -OptionsStruct.ResidualTolerance = 1E-05; +OptionsStruct.ResidualTolerance = 1E-08; OptionsStruct.NoiseScaleFactor = 0.01; OptionsStruct.PlotLive = false; diff --git a/Dipolar-Gas-Simulator/+Scripts/run_on_cluster_vary_aspect_ratio.m b/Dipolar-Gas-Simulator/+Scripts/run_on_cluster_vary_aspect_ratio.m index 56dc6b8..1005272 100644 --- a/Dipolar-Gas-Simulator/+Scripts/run_on_cluster_vary_aspect_ratio.m +++ b/Dipolar-Gas-Simulator/+Scripts/run_on_cluster_vary_aspect_ratio.m @@ -13,16 +13,17 @@ VerticalTrapFrequency = AspectRatio * HorizontalTrapFrequency; OptionsStruct.TrapFrequencies = [HorizontalTrapFrequency, HorizontalTrapFrequency, VerticalTrapFrequency]; OptionsStruct.TrapPotentialType = 'Harmonic'; -OptionsStruct.NumberOfGridPoints = [256, 256, 128]; +OptionsStruct.NumberOfGridPoints = [128, 128, 64]; OptionsStruct.Dimensions = [18, 18, 18]; +OptionsStruct.UseApproximationForLHY = true; OptionsStruct.IncludeDDICutOff = true; OptionsStruct.CutoffType = 'Cylindrical'; OptionsStruct.SimulationMode = 'ImaginaryTimeEvolution'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution' -OptionsStruct.TimeStepSize = 2E-3; % in s -OptionsStruct.MinimumTimeStepSize = 1E-6; % in s +OptionsStruct.TimeStepSize = 1E-4; % in s +OptionsStruct.MinimumTimeStepSize = 2E-10; % in s OptionsStruct.TimeCutOff = 2E6; % in s OptionsStruct.EnergyTolerance = 5E-10; -OptionsStruct.ResidualTolerance = 1E-05; +OptionsStruct.ResidualTolerance = 1E-08; OptionsStruct.NoiseScaleFactor = 0.01; OptionsStruct.PlotLive = false; @@ -55,16 +56,60 @@ VerticalTrapFrequency = AspectRatio * HorizontalTrapFrequency; OptionsStruct.TrapFrequencies = [HorizontalTrapFrequency, HorizontalTrapFrequency, VerticalTrapFrequency]; OptionsStruct.TrapPotentialType = 'Harmonic'; -OptionsStruct.NumberOfGridPoints = [256, 256, 128]; +OptionsStruct.NumberOfGridPoints = [128, 128, 64]; OptionsStruct.Dimensions = [18, 18, 18]; +OptionsStruct.UseApproximationForLHY = true; OptionsStruct.IncludeDDICutOff = true; OptionsStruct.CutoffType = 'Cylindrical'; OptionsStruct.SimulationMode = 'ImaginaryTimeEvolution'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution' -OptionsStruct.TimeStepSize = 2E-3; % in s -OptionsStruct.MinimumTimeStepSize = 1E-6; % in s +OptionsStruct.TimeStepSize = 1E-4; % in s +OptionsStruct.MinimumTimeStepSize = 2E-10; % in s OptionsStruct.TimeCutOff = 2E6; % in s OptionsStruct.EnergyTolerance = 5E-10; -OptionsStruct.ResidualTolerance = 1E-05; +OptionsStruct.ResidualTolerance = 1E-08; +OptionsStruct.NoiseScaleFactor = 0.01; + +OptionsStruct.PlotLive = false; +OptionsStruct.JobNumber = 0; +OptionsStruct.RunOnGPU = true; +OptionsStruct.SaveData = true; +OptionsStruct.SaveDirectory = sprintf('./Results/Data_3D/AspectRatio%s', strrep(num2str(AspectRatio), '.', '_')); +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(); + +%% - Aspect Ratio: 3.8 + +OptionsStruct = struct; + +OptionsStruct.NumberOfAtoms = 5E5; +OptionsStruct.DipolarPolarAngle = deg2rad(0); +OptionsStruct.DipolarAzimuthAngle = 0; +OptionsStruct.ScatteringLength = 85; + +AspectRatio = 3.8; +HorizontalTrapFrequency = 125; +VerticalTrapFrequency = AspectRatio * HorizontalTrapFrequency; +OptionsStruct.TrapFrequencies = [HorizontalTrapFrequency, HorizontalTrapFrequency, VerticalTrapFrequency]; +OptionsStruct.TrapPotentialType = 'Harmonic'; + +OptionsStruct.NumberOfGridPoints = [128, 128, 64]; +OptionsStruct.Dimensions = [18, 18, 18]; +OptionsStruct.UseApproximationForLHY = true; +OptionsStruct.IncludeDDICutOff = true; +OptionsStruct.CutoffType = 'Cylindrical'; +OptionsStruct.SimulationMode = 'ImaginaryTimeEvolution'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution' +OptionsStruct.TimeStepSize = 1E-4; % in s +OptionsStruct.MinimumTimeStepSize = 2E-10; % in s +OptionsStruct.TimeCutOff = 2E6; % in s +OptionsStruct.EnergyTolerance = 5E-10; +OptionsStruct.ResidualTolerance = 1E-08; OptionsStruct.NoiseScaleFactor = 0.01; OptionsStruct.PlotLive = false; diff --git a/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/DipolarGas.m b/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/DipolarGas.m index 49a785f..c30b382 100644 --- a/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/DipolarGas.m +++ b/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/DipolarGas.m @@ -54,7 +54,7 @@ classdef DipolarGas < handle & matlab.mixin.Copyable addParameter(p, 'Dimensions', 10 * ones(1,3),... @(x) assert(isnumeric(x) && isvector(x) && all(x > 0))); addParameter(p, 'SimulationMode', 'ImaginaryTimeEvolution',... - @(x) assert(any(strcmpi(x,{'ImaginaryTimeEvolution','RealTimeEvolution'})))); + @(x) assert(any(strcmpi(x,{'ImaginaryTimeEvolution','RealTimeEvolution', 'EnergyMinimization'})))); addParameter(p, 'CutOffType', 'Cylindrical',... @(x) assert(any(strcmpi(x,{'None', 'Cylindrical','CylindricalInfiniteZ', 'Spherical', 'CustomCylindrical'})))); addParameter(p, 'TimeStepSize', 5E-4,... diff --git a/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/minimizeEnergyFunctional.m b/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/minimizeEnergyFunctional.m new file mode 100644 index 0000000..d4cfd63 --- /dev/null +++ b/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/minimizeEnergyFunctional.m @@ -0,0 +1,105 @@ +function [psi] = minimizeEnergyFunctional(this,psi,Params,Transf,VDk,V) + + format long; + + % Define the function handle + f = @(X) this.Calculator.calculateTotalEnergy(X, Params, Transf, VDk, V)/Params.N; + + % Convergence Criteria: + Epsilon = 1E-5; + + % Gradient Calculation + J = compute_gradient(psi, Params, Transf, VDk, V); % Gradient at initial point + S = -J; % Initial search direction + + % Iteration Counter: + i = 1; + + % Minimization + while norm(S(:)) > Epsilon % Halting Criterion + % Step Size Optimization (Line Search) + alpha = optimize_step_size(f, psi, S, Params, Transf, VDk, V); + + % Update the solution + x_o_new = psi + alpha * S; + + % Update the gradient and search direction + J_old = J; + J = compute_gradient(x_o_new, Params, Transf, VDk, V); + S = update_search_direction(S, J, J_old); + + % Update for next iteration + psi = x_o_new; + Norm = sum(abs(psi(:)).^2)*Transf.dx*Transf.dy*Transf.dz; + psi = sqrt(Params.N)*psi/sqrt(Norm); + i = i + 1; + end + + disp('Minimum found!'); + fprintf('Number of Iterations for Convergence: %d\n\n', i); + + muchem = this.Calculator.calculateChemicalPotential(psi,Params,Transf,VDk,V); + + disp('Saving data...'); + save(sprintf(strcat(this.SaveDirectory, '/Run_%03i/psi_gs.mat'),Params.njob),'psi','muchem','Transf','Params','VDk','V'); + disp('Save complete!'); + +end + +%% Helper functions +% Numerical Gradient Calculation using the finite differences method +function J = compute_gradient(psi, Params, Transf, VDk, V) + + % Kinetic energy + KEop = 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2); + H_o = KEop + V; + + %Contact interactions + Hint = Params.gs*abs(psi).^2; + + % DDIs + frho = fftn(abs(psi).^2); + Phi = real(ifftn(frho.*VDk)); + Hddi = Params.gdd*Phi; + + %Quantum fluctuations + Hqf = Params.gammaQF*abs(psi).^3; + + H = H_o + Hint + Hddi + Hqf; + + J = H.*psi; + +end + +% Line Search (Step Size Optimization) +function alpha = optimize_step_size(f, X, S, Params, Transf, VDk, V) + alpha = 1; % Initial step size + rho = 0.5; % Step size reduction factor + c = 1E-4; % Armijo condition constant + max_iter = 100; % Max iterations for backtracking + tol = 1E-8; % Tolerance for stopping + + grad = compute_gradient(X, Params, Transf, VDk, V); % Compute gradient once + f_X = f(X); % Evaluate f(X) once + + for k = 1:max_iter + % Evaluate Armijo condition with precomputed f(X) and grad + if f(X + alpha * S) <= f_X + c * alpha * (S(:)' * grad(:)) + break; + else + alpha = rho * alpha; % Reduce the step size + end + + % Early stopping if step size becomes too small + if alpha < tol + break; + end + end + +end + +% Update Search Direction (Polak-Ribiere method) +function S_new = update_search_direction(S, J_new, J_old) + beta = (norm(J_new(:))^2) / (norm(J_old(:))^2); + S_new = -J_new + beta * S; +end \ No newline at end of file diff --git a/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/run.m b/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/run.m index 3f4dcf9..b4d0643 100644 --- a/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/run.m +++ b/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/run.m @@ -19,5 +19,9 @@ function [Params, Transf, psi,V,VDk] = run(this) % --- Run Simulation --- mkdir(sprintf(this.SaveDirectory)) mkdir(sprintf(strcat(this.SaveDirectory, '/Run_%03i'),Params.njob)) - [psi] = this.propagateWavefunction(psi,Params,Transf,VDk,V,t_idx,Observ); + if strcmp(this.SimulationMode, 'EnergyMinimization') + [psi] = this.minimizeEnergyFunctional(psi,Params,Transf,VDk,V); + else + [psi] = this.propagateWavefunction(psi,Params,Transf,VDk,V,t_idx,Observ); + end end