745 lines
31 KiB
Matlab
745 lines
31 KiB
Matlab
%% This script is testing the functionalities of the Dipolar Gas Simulator
|
|
%
|
|
% Important: Run only sectionwise!!
|
|
|
|
%% - Create Simulator, Potential and Calculator object with specified options
|
|
|
|
OptionsStruct = struct;
|
|
|
|
OptionsStruct.NumberOfAtoms = 40000;
|
|
OptionsStruct.DipolarPolarAngle = deg2rad(0);
|
|
OptionsStruct.DipolarAzimuthAngle = 0;
|
|
OptionsStruct.ScatteringLength = 95;
|
|
|
|
OptionsStruct.TrapFrequencies = [30, 60, 90];
|
|
OptionsStruct.TrapPotentialType = 'Harmonic';
|
|
|
|
OptionsStruct.NumberOfGridPoints = [128, 64, 64];
|
|
OptionsStruct.Dimensions = [30, 20, 20];
|
|
OptionsStruct.UseApproximationForLHY = true;
|
|
OptionsStruct.IncludeDDICutOff = true;
|
|
OptionsStruct.CutoffType = 'Cylindrical';
|
|
OptionsStruct.SimulationMode = 'EnergyMinimization'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution' | 'EnergyMinimization'
|
|
OptionsStruct.GradientDescentMethod = 'NonLinearCGD'; % 'HeavyBall' | 'NonLinearCGD'
|
|
OptionsStruct.MaxIterationsForGD = 1000;
|
|
OptionsStruct.TimeStepSize = 1E-3; % in s
|
|
OptionsStruct.MinimumTimeStepSize = 1E-6; % in s
|
|
OptionsStruct.TimeCutOff = 2E6; % in s
|
|
OptionsStruct.EnergyTolerance = 5E-10;
|
|
OptionsStruct.ResidualTolerance = 1E-08;
|
|
OptionsStruct.NoiseScaleFactor = 0.010;
|
|
|
|
OptionsStruct.PlotLive = true;
|
|
OptionsStruct.JobNumber = 0;
|
|
OptionsStruct.RunOnGPU = false;
|
|
OptionsStruct.SaveData = true;
|
|
OptionsStruct.SaveDirectory = './Results/Data_3D/GradientDescent'; % './Results/Data_3D/AnisotropicTrap/Tilted0'
|
|
options = Helper.convertstruct2cell(OptionsStruct);
|
|
|
|
sim = Simulator.DipolarGas(options{:});
|
|
pot = Simulator.Potentials(options{:});
|
|
sim.Potential = pot.trap(); % + pot.repulsive_chopstick();
|
|
|
|
%-% Run Simulation %-%
|
|
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);
|
|
|
|
%% - Plot numerical grid
|
|
% Plotter.visualizeSpace(Transf)
|
|
%% - Plot trap potential
|
|
Plotter.visualizeTrapPotential(sim.Potential,Params,Transf)
|
|
%% - Plot initial wavefunction
|
|
Plotter.visualizeWavefunction(psi,Params,Transf)
|
|
|
|
%% Imaginary-Time followed by Real-time
|
|
|
|
% - Imaginary-Time
|
|
OptionsStruct = struct;
|
|
|
|
OptionsStruct.NumberOfAtoms = 1E5;
|
|
OptionsStruct.DipolarPolarAngle = deg2rad(0);
|
|
OptionsStruct.DipolarAzimuthAngle = deg2rad(0);
|
|
OptionsStruct.ScatteringLength = 110;
|
|
|
|
OptionsStruct.TrapFrequencies = [50, 20, 150];
|
|
OptionsStruct.TrapPotentialType = 'Harmonic';
|
|
|
|
OptionsStruct.NumberOfGridPoints = [64, 128, 32];
|
|
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 = true;
|
|
OptionsStruct.JobNumber = 0;
|
|
OptionsStruct.RunOnGPU = false;
|
|
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 = 10E-3;
|
|
OptionsStruct.HoldTime = 10E-3;
|
|
OptionsStruct.QuenchScatteringLength = false;
|
|
OptionsStruct.RotateDipoles = true;
|
|
OptionsStruct.FinalScatteringLength = 88;
|
|
OptionsStruct.FinalDipolarPolarAngle = deg2rad(35);
|
|
OptionsStruct.FinalDipolarAzimuthAngle = deg2rad(0);
|
|
OptionsStruct.TimeStepSize = 1E-3; % in s
|
|
OptionsStruct.NoiseScaleFactor = 0.010;
|
|
|
|
OptionsStruct.PlotLive = true;
|
|
OptionsStruct.JobNumber = 0;
|
|
OptionsStruct.RunOnGPU = false;
|
|
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();
|
|
|
|
%%
|
|
% To reproduce results from the Blair Blakie paper:
|
|
% (n*add^2, as/add)
|
|
% Critical point: (0.0978, 0.784); Triangular phase: (0.0959, 0.750); Stripe phase: (0.144, 0.765); Honeycomb phase: (0.192, 0.780)
|
|
|
|
% N = ((n*add^2)/Params.add^2) * (Params.Lx *1E-6)^2
|
|
% Critical point: N = 2.0427e+07; Triangular phase: N = 2.0030e+07; Stripe phase: N = 3.0077e+07; Honeycomb phase: N = 4.0102e+07 for dimensions fixed to 100
|
|
|
|
% as = ((as/add)*Params.add)/Params.a0
|
|
% Critical point: 102.5133; Triangular phase: 98.0676; Stripe phase: 100.0289; Honeycomb phase: 101.9903
|
|
|
|
%% - Create Variational2D and Calculator object with specified options
|
|
|
|
OptionsStruct = struct;
|
|
|
|
OptionsStruct.NumberOfAtoms = 2.0030e+07;
|
|
OptionsStruct.DipolarPolarAngle = 0;
|
|
OptionsStruct.DipolarAzimuthAngle = 0;
|
|
OptionsStruct.ScatteringLength = 98.0676;
|
|
|
|
OptionsStruct.TrapFrequencies = [0, 0, 72.4];
|
|
OptionsStruct.TrapPotentialType = 'None';
|
|
|
|
OptionsStruct.NumberOfGridPoints = [128, 128];
|
|
OptionsStruct.Dimensions = [100, 100];
|
|
OptionsStruct.TimeStepSize = 0.005; % in s
|
|
OptionsStruct.MinimumTimeStepSize = 1E-5; % in s
|
|
OptionsStruct.TimeCutOff = 1E6; % in s
|
|
OptionsStruct.EnergyTolerance = 5E-10;
|
|
OptionsStruct.ResidualTolerance = 1E-05;
|
|
OptionsStruct.NoiseScaleFactor = 0.05;
|
|
|
|
OptionsStruct.MaxIterations = 10;
|
|
OptionsStruct.VariationalWidth = 5.0;
|
|
OptionsStruct.WidthLowerBound = 1;
|
|
OptionsStruct.WidthUpperBound = 12;
|
|
OptionsStruct.WidthCutoff = 5e-3;
|
|
|
|
OptionsStruct.PlotLive = true;
|
|
OptionsStruct.JobNumber = 1;
|
|
OptionsStruct.RunOnGPU = false;
|
|
OptionsStruct.SaveData = true;
|
|
OptionsStruct.SaveDirectory = './Results/Data_TriangularPhase';
|
|
options = Helper.convertstruct2cell(OptionsStruct);
|
|
clear OptionsStruct
|
|
|
|
solver = VariationalSolver2D.DipolarGas(options{:});
|
|
pot = VariationalSolver2D.Potentials(options{:});
|
|
solver.Potential = pot.trap();
|
|
|
|
%-% Run Solver %-%
|
|
[Params, Transf, psi, V, VDk] = solver.run();
|
|
%%
|
|
% Solve BdG equations
|
|
% Load data
|
|
Data = load(sprintf(horzcat(solver.SaveDirectory, '/Run_%03i/psi_gs.mat'),solver.JobNumber),'psi','Transf','Params','VParams','V');
|
|
|
|
Transf = Data.Transf;
|
|
Params = Data.Params;
|
|
VParams = Data.VParams;
|
|
V = Data.V;
|
|
|
|
if isgpuarray(Data.psi)
|
|
psi = gather(Data.psi);
|
|
else
|
|
psi = Data.psi;
|
|
end
|
|
|
|
% == DDI potential == %
|
|
VDk = solver.Calculator.calculateVDkWithCutoff(Transf, Params, VParams.ell);
|
|
|
|
% == Chemical potential == %
|
|
muchem = solver.Calculator.calculateChemicalPotential(psi,Params,VParams,Transf,VDk,V);
|
|
|
|
[evals, modes] = BdGSolver2D.solveBogoliubovdeGennesIn2D(psi, Params, VDk, VParams, Transf, muchem);
|
|
|
|
% Save the eigenvalues and eigenvectors to a .mat file
|
|
save(sprintf(strcat(solver.SaveDirectory, '/Run_%03i/bdg_eigen_data.mat'),solver.JobNumber), 'evals', 'modes', '-v7.3');
|
|
|
|
%% - Create Variational2D and Calculator object with specified options
|
|
|
|
OptionsStruct = struct;
|
|
|
|
OptionsStruct.NumberOfAtoms = 3.0077e+07;
|
|
OptionsStruct.DipolarPolarAngle = 0;
|
|
OptionsStruct.DipolarAzimuthAngle = 0;
|
|
OptionsStruct.ScatteringLength = 100.0289;
|
|
|
|
OptionsStruct.TrapFrequencies = [0, 0, 72.4];
|
|
OptionsStruct.TrapPotentialType = 'None';
|
|
|
|
OptionsStruct.NumberOfGridPoints = [128, 128];
|
|
OptionsStruct.Dimensions = [100, 100];
|
|
OptionsStruct.TimeStepSize = 0.005; % in s
|
|
OptionsStruct.MinimumTimeStepSize = 1E-5; % in s
|
|
OptionsStruct.TimeCutOff = 1E6; % in s
|
|
OptionsStruct.EnergyTolerance = 5E-10;
|
|
OptionsStruct.ResidualTolerance = 1E-05;
|
|
OptionsStruct.NoiseScaleFactor = 0.05;
|
|
|
|
OptionsStruct.MaxIterations = 10;
|
|
OptionsStruct.VariationalWidth = 5;
|
|
OptionsStruct.WidthLowerBound = 1;
|
|
OptionsStruct.WidthUpperBound = 12;
|
|
OptionsStruct.WidthCutoff = 5e-3;
|
|
|
|
OptionsStruct.PlotLive = true;
|
|
OptionsStruct.JobNumber = 1;
|
|
OptionsStruct.RunOnGPU = false;
|
|
OptionsStruct.SaveData = true;
|
|
OptionsStruct.SaveDirectory = './Results/Data_StripePhase';
|
|
options = Helper.convertstruct2cell(OptionsStruct);
|
|
clear OptionsStruct
|
|
|
|
solver = VariationalSolver2D.DipolarGas(options{:});
|
|
pot = VariationalSolver2D.Potentials(options{:});
|
|
solver.Potential = pot.trap();
|
|
|
|
%-% Run Solver %-%
|
|
[Params, Transf, psi, V, VDk] = solver.run();
|
|
|
|
%% - Create Variational2D and Calculator object with specified options
|
|
|
|
OptionsStruct = struct;
|
|
|
|
OptionsStruct.NumberOfAtoms = 4.2e+07;
|
|
OptionsStruct.DipolarPolarAngle = 0;
|
|
OptionsStruct.DipolarAzimuthAngle = 0;
|
|
OptionsStruct.ScatteringLength = 101.35;
|
|
|
|
OptionsStruct.TrapFrequencies = [0, 0, 72.4];
|
|
OptionsStruct.TrapPotentialType = 'None';
|
|
|
|
OptionsStruct.NumberOfGridPoints = [128, 128];
|
|
OptionsStruct.Dimensions = [100, 100];
|
|
OptionsStruct.TimeStepSize = 0.005; % in s
|
|
OptionsStruct.MinimumTimeStepSize = 1E-5; % in s
|
|
OptionsStruct.TimeCutOff = 2E6; % in s
|
|
OptionsStruct.EnergyTolerance = 5E-10;
|
|
OptionsStruct.ResidualTolerance = 1E-05;
|
|
OptionsStruct.NoiseScaleFactor = 0.05;
|
|
|
|
OptionsStruct.MaxIterations = 10;
|
|
OptionsStruct.VariationalWidth = 6;
|
|
OptionsStruct.WidthLowerBound = 1;
|
|
OptionsStruct.WidthUpperBound = 12;
|
|
OptionsStruct.WidthCutoff = 5e-3;
|
|
|
|
OptionsStruct.PlotLive = true;
|
|
OptionsStruct.JobNumber = 1;
|
|
OptionsStruct.RunOnGPU = false;
|
|
OptionsStruct.SaveData = true;
|
|
OptionsStruct.SaveDirectory = './Results/Data_HoneycombPhase';
|
|
options = Helper.convertstruct2cell(OptionsStruct);
|
|
clear OptionsStruct
|
|
|
|
solver = VariationalSolver2D.DipolarGas(options{:});
|
|
pot = VariationalSolver2D.Potentials(options{:});
|
|
solver.Potential = pot.trap();
|
|
|
|
%-% Run Solver %-%
|
|
[Params, Transf, psi, V, VDk] = solver.run();
|
|
|
|
%% - Plot numerical grid
|
|
% Plotter.visualizeSpace2D(Transf)
|
|
%% - Plot trap potential
|
|
% Plotter.visualizeTrapPotential2D(solver.Potential,Params,Transf)
|
|
%% - Plot initial wavefunction
|
|
Plotter.visualizeWavefunction2D(psi,Params,Transf)
|
|
%% - Plot GS wavefunction
|
|
Plotter.visualizeGSWavefunction2D(solver.SaveDirectory, solver.JobNumber)
|
|
|
|
%% - Analysis
|
|
SaveDirectory = './Results/Data_TriangularPhase';
|
|
JobNumber = 1;
|
|
Plotter.visualizeGSWavefunction2D(SaveDirectory, JobNumber)
|
|
|
|
%% - Analysis
|
|
SaveDirectory = './Results/Data_StripePhase';
|
|
JobNumber = 2;
|
|
Plotter.visualizeGSWavefunction2D(SaveDirectory, JobNumber)
|
|
|
|
%% - Analysis
|
|
SaveDirectory = './Results/Data_HoneycombPhase';
|
|
JobNumber = 3;
|
|
Plotter.visualizeGSWavefunction2D(SaveDirectory, JobNumber)
|
|
|
|
%% Tilting of the dipoles
|
|
% Atom Number = 1.00e+05
|
|
% System size = [10, 10]
|
|
|
|
% - Analysis
|
|
SaveDirectory = './Results/Data_TiltingOfDipoles';
|
|
JobNumber = 5;
|
|
Plotter.visualizeGSWavefunction2D(SaveDirectory, JobNumber)
|
|
|
|
%% Tilting of the dipoles
|
|
% Atom Number = 1.00e+05
|
|
% System size = [5 * l_rot, 5 * l_rot]
|
|
|
|
%% - Analysis
|
|
SaveDirectory = './Results/Data_TiltingOfDipoles/AdjustedSystemSize/Hz500';
|
|
JobNumber = 2;
|
|
Plotter.visualizeGSWavefunction2D(SaveDirectory, JobNumber)
|
|
[contrast, period_X, period_Y] = Scripts.analyzeGSWavefunction(SaveDirectory, JobNumber);
|
|
|
|
%% - Analysis
|
|
SaveDirectory = './Results/Data_TiltingOfDipoles/AdjustedSystemSize/Hz750';
|
|
JobNumber = 1;
|
|
% Plotter.visualizeGSWavefunction2D(SaveDirectory, JobNumber)
|
|
[contrast, period_X, period_Y] = Scripts.analyzeGSWavefunction(SaveDirectory, JobNumber);
|
|
|
|
%% - Analysis
|
|
SaveDirectory = './Results/Data_TiltingOfDipoles/AdjustedSystemSize/Hz1000';
|
|
JobNumber = 1;
|
|
% Plotter.visualizeGSWavefunction2D(SaveDirectory, JobNumber)
|
|
[contrast, period_X, period_Y] = Scripts.analyzeGSWavefunction(SaveDirectory, JobNumber);
|
|
|
|
%% - Analysis
|
|
SaveDirectory = './Results/Data_TiltingOfDipoles/AdjustedSystemSize/Hz2000';
|
|
JobNumber = 1;
|
|
% Plotter.visualizeGSWavefunction2D(SaveDirectory, JobNumber)
|
|
[contrast, period_X, period_Y] = Scripts.analyzeGSWavefunction(SaveDirectory, JobNumber);
|
|
|
|
%% - Analysis
|
|
SaveDirectory = './Results/Data_TiltingOfDipoles/TransitionAngle/Hz500';
|
|
NumberOfJobs = 16;
|
|
SaveOption = 'video';
|
|
[contrast_array, periodX_array, periodY_array] = Scripts.analyzeGSWavefunction_multipleruns(SaveDirectory, NumberOfJobs, SaveOption);
|
|
|
|
%% - Analysis
|
|
SaveDirectory = './Results/Data_TiltingOfDipoles/TransitionAngle/Hz750';
|
|
NumberOfJobs = 16;
|
|
SaveOption = 'video';
|
|
[contrast_array, periodX_array, periodY_array] = Scripts.analyzeGSWavefunction_multipleruns(SaveDirectory, NumberOfJobs, SaveOption);
|
|
|
|
%% - Analysis
|
|
SaveDirectory = './Results/Data_TiltingOfDipoles/TransitionAngle/Hz1000';
|
|
NumberOfJobs = 16;
|
|
SaveOption = 'video';
|
|
[contrast_array, periodX_array, periodY_array] = Scripts.analyzeGSWavefunction_multipleruns(SaveDirectory, NumberOfJobs, SaveOption);
|
|
|
|
%% - Analysis
|
|
SaveDirectory = './Results/Data_TiltingOfDipoles/TransitionAngle/Hz2000';
|
|
NumberOfJobs = 16;
|
|
SaveOption = 'video';
|
|
[contrast_array, periodX_array, periodY_array] = Scripts.analyzeGSWavefunction_multipleruns(SaveDirectory, NumberOfJobs, SaveOption);
|
|
%% - Analysis
|
|
SaveDirectory = './Results/Data_TiltingOfDipoles/TransitionAngle_newruns/Hz1000';
|
|
JobNumber = 0;
|
|
% Plotter.visualizeGSWavefunction2D(SaveDirectory, JobNumber)
|
|
[contrast, period_X, period_Y] = Scripts.analyzeGSWavefunction(SaveDirectory, JobNumber);
|
|
%% - Analysis
|
|
SaveDirectory = './Results/Data_TiltingOfDipoles/HarmonicTrap/Hz500';
|
|
JobNumber = 0;
|
|
Plotter.visualizeGSWavefunction2D(SaveDirectory, JobNumber)
|
|
[contrast, period_X, period_Y] = Scripts.analyzeGSWavefunction_in_plane_trap(SaveDirectory, JobNumber);
|
|
|
|
%% - Analysis
|
|
SaveDirectory = './Results/Data_TiltingOfDipoles/HarmonicTrap/Hz750';
|
|
JobNumber = 1;
|
|
% Plotter.visualizeGSWavefunction2D(SaveDirectory, JobNumber)
|
|
[contrast, period_X, period_Y] = Scripts.analyzeGSWavefunction_in_plane_trap(SaveDirectory, JobNumber);
|
|
|
|
%% - Analysis
|
|
SaveDirectory = './Results/Data_TiltingOfDipoles/HarmonicTrap/Hz1000';
|
|
JobNumber = 1;
|
|
% Plotter.visualizeGSWavefunction2D(SaveDirectory, JobNumber)
|
|
[contrast, period_X, period_Y] = Scripts.analyzeGSWavefunction_in_plane_trap(SaveDirectory, JobNumber);
|
|
|
|
%% - Analysis
|
|
SaveDirectory = './Results/Data_TiltingOfDipoles/HarmonicTrap/Hz2000';
|
|
JobNumber = 1;
|
|
% Plotter.visualizeGSWavefunction2D(SaveDirectory, JobNumber)
|
|
[contrast, period_X, period_Y] = Scripts.analyzeGSWavefunction_in_plane_trap(SaveDirectory, JobNumber);
|
|
|
|
%%
|
|
SaveDirectory = './Results/Data_TiltingOfDipoles/TransitionAngle/OptimalSystemSize/Hz500/Degree0';
|
|
% Define the desired range for Lx and Ly
|
|
Lx_Range = 4:0.2:6; % Extend Lx from 5 to 10
|
|
Ly_Range = 4:0.2:6; % Extend Ly from 5 to 10
|
|
MinEnergyDataArray = Scripts.analyzeGSWavefunction_optimal_system_size(SaveDirectory, Lx_Range, Ly_Range);
|
|
|
|
%%
|
|
% SaveDirectory = './Results/Data_TiltingOfDipoles/TransitionAngle/OptimalSystemSize/Hz500/Degree0';
|
|
SaveDirectory = './Results/Data_TiltingOfDipoles/TransitionAngle/OptimalSystemSizeWithBiasAnsatz/Hz500/Degree0';
|
|
% Define the desired range
|
|
% LatticeSpacing = 1.0:0.05:4.0;
|
|
LatticeSpacing = linspace(1.6,2.25,61);
|
|
MinEnergyDataArray = Scripts.analyzeGSWavefunction_constrained_optimal_system_size(SaveDirectory, LatticeSpacing);
|
|
[val, idx] = min(MinEnergyDataArray(:,3));
|
|
OptimalSize = MinEnergyDataArray(idx,1:2);
|
|
Plotter.visualizeGSWavefunction2D(SaveDirectory, idx)
|
|
[contrast, periodX, periodY] = Scripts.analyzeGSWavefunction(SaveDirectory, idx, false);
|
|
%%
|
|
% SaveDirectory = './Results/Data_TiltingOfDipoles/TransitionAngle/OptimalSystemSize/Hz500/Degree5';
|
|
SaveDirectory = './Results/Data_TiltingOfDipoles/TransitionAngle/OptimalSystemSizeWithBiasAnsatz/Hz500/Degree5';
|
|
% Define the desired range
|
|
% LatticeSpacing = 1.0:0.05:4.0;
|
|
LatticeSpacing = linspace(1.6,2.25,61);
|
|
MinEnergyDataArray = Scripts.analyzeGSWavefunction_constrained_optimal_system_size(SaveDirectory, LatticeSpacing);
|
|
[val, idx] = min(MinEnergyDataArray(:,3));
|
|
OptimalSize = MinEnergyDataArray(idx,1:2);
|
|
Plotter.visualizeGSWavefunction2D(SaveDirectory, idx)
|
|
[contrast, periodX, periodY] = Scripts.analyzeGSWavefunction(SaveDirectory, idx, false);
|
|
|
|
%%
|
|
% SaveDirectory = './Results/Data_TiltingOfDipoles/TransitionAngle/OptimalSystemSize/Hz500/Degree7_5';
|
|
SaveDirectory = './Results/Data_TiltingOfDipoles/TransitionAngle/OptimalSystemSizeWithBiasAnsatz/Hz500/Degree7_5';
|
|
% Define the desired range
|
|
% LatticeSpacing = 1.0:0.05:4.0;
|
|
LatticeSpacing = linspace(1.6,2.25,61);
|
|
MinEnergyDataArray = Scripts.analyzeGSWavefunction_constrained_optimal_system_size(SaveDirectory, LatticeSpacing);
|
|
[val, idx] = min(MinEnergyDataArray(:,3));
|
|
OptimalSize = MinEnergyDataArray(idx,1:2);
|
|
Plotter.visualizeGSWavefunction2D(SaveDirectory, idx)
|
|
[contrast, periodX, periodY] = Scripts.analyzeGSWavefunction(SaveDirectory, idx, false);
|
|
|
|
%%
|
|
% SaveDirectory = './Results/Data_TiltingOfDipoles/TransitionAngle/OptimalSystemSize/Hz500/Degree10';
|
|
SaveDirectory = './Results/Data_TiltingOfDipoles/TransitionAngle/OptimalSystemSizeWithBiasAnsatz/Hz500/Degree10';
|
|
% Define the desired range
|
|
% LatticeSpacing = 1.0:0.05:4.0;
|
|
LatticeSpacing = linspace(1.6,2.25,61);
|
|
MinEnergyDataArray = Scripts.analyzeGSWavefunction_constrained_optimal_system_size(SaveDirectory, LatticeSpacing);
|
|
[val, idx] = min(MinEnergyDataArray(:,3));
|
|
OptimalSize = MinEnergyDataArray(idx,1:2);
|
|
Plotter.visualizeGSWavefunction2D(SaveDirectory, idx)
|
|
[contrast, periodX, periodY] = Scripts.analyzeGSWavefunction(SaveDirectory, idx, false);
|
|
|
|
%%
|
|
% SaveDirectory = './Results/Data_TiltingOfDipoles/TransitionAngle/OptimalSystemSize/Hz500/Degree15';
|
|
SaveDirectory = './Results/Data_TiltingOfDipoles/TransitionAngle/OptimalSystemSizeWithBiasAnsatz/Hz500/Degree15';
|
|
% Define the desired range
|
|
% LatticeSpacing = 1.0:0.05:4.0;
|
|
LatticeSpacing = linspace(1.6,2.25,61);
|
|
MinEnergyDataArray = Scripts.analyzeGSWavefunction_constrained_optimal_system_size(SaveDirectory, LatticeSpacing);
|
|
[val, idx] = min(MinEnergyDataArray(:,3));
|
|
OptimalSize = MinEnergyDataArray(idx,1:2);
|
|
Plotter.visualizeGSWavefunction2D(SaveDirectory, idx)
|
|
[contrast, periodX, periodY] = Scripts.analyzeGSWavefunction(SaveDirectory, idx, false);
|
|
|
|
%% - Analysis
|
|
SaveDirectory = './Results/Data_TiltingOfDipoles/HarmonicTrap/AspectRatio/AR2_8';
|
|
JobNumber = 0; % 79
|
|
% JobNumber = 1; % 80
|
|
% JobNumber = 2; % 81
|
|
Plotter.visualizeGSWavefunction2D(SaveDirectory, JobNumber)
|
|
[contrast, period_X, period_Y] = Scripts.analyzeGSWavefunction_in_plane_trap(SaveDirectory, JobNumber);
|
|
|
|
%% - Analysis
|
|
SaveDirectory = './Results/Data_TiltingOfDipoles/HarmonicTrap/AspectRatio/AR3_7';
|
|
% JobNumber = 0; % 80
|
|
% JobNumber = 1; % 81
|
|
JobNumber = 2; % 82
|
|
% JobNumber = 3; % 83
|
|
Plotter.visualizeGSWavefunction2D(SaveDirectory, JobNumber)
|
|
[contrast, period_X, period_Y] = Scripts.analyzeGSWavefunction_in_plane_trap(SaveDirectory, JobNumber);
|
|
%% - Plot GS wavefunction
|
|
% 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/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';
|
|
% SaveDirectory = './Results/Data_3D/ApproximateLHY/AspectRatio3_9';
|
|
JobNumber = 3;
|
|
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)
|
|
%%
|
|
SaveDirectory = './Results/Data_3D/TiltedDipoles0';
|
|
JobNumber = 0;
|
|
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
|
|
%%
|
|
SaveDirectory = './Results/Data_3D/TiltedDipoles15';
|
|
JobNumber = 2;
|
|
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
|
|
%%
|
|
SaveDirectory = './Results/Data_3D/TiltedDipoles30';
|
|
JobNumber = 2;
|
|
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
|
|
%%
|
|
SaveDirectory = './Results/Data_3D/TiltedDipoles45';
|
|
JobNumber = 2;
|
|
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
|
|
%%
|
|
SaveDirectory = './Results/Data_3D/RealTimeDynamics';
|
|
JobNumber = 0;
|
|
Plotter.makeMovie(SaveDirectory, JobNumber)
|
|
%%
|
|
SaveDirectory = './Results/Data_3D/AnisotropicTrap/TiltedDipoles15';
|
|
JobNumber = 0;
|
|
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
|
|
%%
|
|
SaveDirectory = './Results/Data_3D/AnisotropicTrap/TiltedDipoles30';
|
|
JobNumber = 0;
|
|
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
|
|
%%
|
|
SaveDirectory = './Results/Data_3D/AnisotropicTrap/TiltedDipoles35';
|
|
JobNumber = 0;
|
|
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
|
|
%%
|
|
SaveDirectory = './Results/Data_3D/AnisotropicTrap/TiltedDipoles40';
|
|
JobNumber = 0;
|
|
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
|
|
%%
|
|
SaveDirectory = './Results/Data_3D/AnisotropicTrap/TiltedDipoles45';
|
|
JobNumber = 0;
|
|
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
|
|
%%
|
|
SaveDirectory = './Results/Data_3D/PhaseDiagram/GradientDescent/';
|
|
JobNumber = 0;
|
|
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
|
|
%%
|
|
SaveDirectory = './Results/Data_3D/PhaseDiagram/ImagTimePropagation/';
|
|
JobNumber = 0;
|
|
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
|
|
%% Generate lists
|
|
|
|
% Set display format
|
|
format longG
|
|
|
|
% Generate the lists
|
|
list1 = linspace(1E5, 5E6, 25);
|
|
list2 = linspace(80, 105, 25);
|
|
|
|
% Convert to strings with no scientific notation
|
|
str_list1 = compose('%.0f', list1);
|
|
str_list2 = compose('%.2f', list2);
|
|
|
|
% Join as space-separated strings
|
|
row1 = strjoin(str_list1', ' ');
|
|
row2 = strjoin(str_list2', ' ');
|
|
|
|
% Display results
|
|
disp(row1)
|
|
disp(row2)
|
|
|
|
%% Use space-separated floating-point/integer values
|
|
SCATTERING_LENGTH_RANGE = "[80.00 81.04 82.08 83.12 84.17 85.21 86.25 87.29 88.33 89.38 90.42 91.46 92.50 93.54 94.58 95.62 96.67 97.71 98.75 99.79 100.83 101.88 102.92 103.96 105.00]";
|
|
NUM_ATOMS_LIST = "[100000 304167 508333 712500 916667 1120833 1325000 1529167 1733333 1937500 2141667 2345833 2550000 2754167 2958333 3162500 3366667 3570833 3775000 3979167 4183333 4387500 4591667 4795833 5000000]";
|
|
|
|
%% Explore the phase diagram
|
|
SCATTERING_LENGTH_RANGE = [91.46 92.50 93.54 94.58 95.62 96.67 97.71];
|
|
NUM_ATOMS_LIST = [100000 304167 508333 712500 916667 1120833 1325000 1529167 1733333 1937500 2141667 2345833 2550000 2754167 2958333 3162500 3366667 3570833 3775000 3979167 4183333 4387500 4591667 4795833 5000000];
|
|
|
|
JobNumber = 0;
|
|
|
|
for i = 1:length(SCATTERING_LENGTH_RANGE)
|
|
aS = SCATTERING_LENGTH_RANGE(i);
|
|
for j = 1:length(NUM_ATOMS_LIST)
|
|
N = NUM_ATOMS_LIST(j);
|
|
|
|
SaveDirectory = sprintf('D:/Results-Numerics/PhaseDiagram/ImagTimePropagation/aS_%.6e_theta_000_phi_000_N_%d', aS, N);
|
|
fprintf('Processing JobNumber %d: %s\n', JobNumber, SaveDirectory);
|
|
|
|
% Call the plotting function
|
|
try
|
|
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
|
|
catch ME
|
|
warning(ME.message)
|
|
continue;
|
|
end
|
|
|
|
% Pause to inspect plot before continuing
|
|
pause(1.5)
|
|
end
|
|
end
|
|
|
|
%% Parameters that need to be run again
|
|
% Cell array of the input strings
|
|
strs = {
|
|
'aS_9.875000e+01_theta_000_phi_000_N_2958333'
|
|
'aS_9.562000e+01_theta_000_phi_000_N_4183333'
|
|
'aS_9.562000e+01_theta_000_phi_000_N_1325000'
|
|
'aS_9.458000e+01_theta_000_phi_000_N_508333'
|
|
'aS_9.458000e+01_theta_000_phi_000_N_5000000'
|
|
'aS_9.458000e+01_theta_000_phi_000_N_1937500'
|
|
'aS_9.354000e+01_theta_000_phi_000_N_4795833'
|
|
'aS_9.250000e+01_theta_000_phi_000_N_5000000'
|
|
'aS_9.250000e+01_theta_000_phi_000_N_4183333'
|
|
'aS_9.146000e+01_theta_000_phi_000_N_1937500'
|
|
'aS_8.833000e+01_theta_000_phi_000_N_2754167'
|
|
'aS_8.833000e+01_theta_000_phi_000_N_1325000'
|
|
'aS_8.729000e+01_theta_000_phi_000_N_3775000'
|
|
'aS_8.625000e+01_theta_000_phi_000_N_712500'
|
|
'aS_8.625000e+01_theta_000_phi_000_N_4795833'
|
|
'aS_8.625000e+01_theta_000_phi_000_N_3979167'
|
|
'aS_8.521000e+01_theta_000_phi_000_N_1937500'
|
|
'aS_8.208000e+01_theta_000_phi_000_N_3979167'
|
|
'aS_105_theta_000_phi_000_N_916667'
|
|
'aS_1.039600e+02_theta_000_phi_000_N_3366667'
|
|
'aS_1.008300e+02_theta_000_phi_000_N_3979167'
|
|
'aS_080_theta_000_phi_000_N_5000000'
|
|
};
|
|
|
|
% Extract values using regex
|
|
tokens = regexp(strs, 'aS_([0-9.e+-]+)_theta_000_phi_000_N_([0-9]+)', 'tokens');
|
|
|
|
% Convert to numeric array
|
|
pairs = cellfun(@(x) [str2double(x{1}), str2double(x{2})], [tokens{:}], 'UniformOutput', false);
|
|
pairs = vertcat(pairs{:});
|
|
|
|
% Sort by a_s (first column)
|
|
pairs = sortrows(pairs, 1);
|
|
|
|
% Display as a table
|
|
disp(array2table(pairs, 'VariableNames', {'a_s', 'N'}));
|
|
|
|
%% Phase diagram for untilted case
|
|
N = [1E5, 3.04E5, 5.08E5, 7.125E5, 9.16E5, 1.12E6, 1.325E6, 1.529E6, ...
|
|
1.733E6, 1.9375E6, 2.141E6, 2.345E6, 2.55E6, 2.75E6, 2.95E6, ...
|
|
3.1625E6, 3.367E6, 3.57E6, 3.775E6, 3.979E6, 4.183E6, 4.3875E6, ...
|
|
4.591E6, 4.795E6, 5E6];
|
|
|
|
as_UB = [88.33, 94.58, 95.62, 96.67, 97.71, 97.71, 97.71, 97.71, 97.71, 97.71, 97.71, 97.71, 97.71, ...
|
|
97.71, 97.71, 97.71, 97.71, 97.71, 97.71, 97.71, 97.71, 97.71, 97.71, 97.71, 97.71];
|
|
|
|
as_LB = [87.29, 93.54, 94.58, 95.62, 96.67, 96.67, 96.67, 96.67, 96.67, 96.67, 96.67, 96.67, 96.67, ...
|
|
96.67, 96.67, 96.67, 96.67, 96.67, 96.67, 96.67, 96.67, 96.67, 96.67, 96.67, 96.67];
|
|
|
|
% Filter only rows with non-NaN data
|
|
valid_idx = ~isnan(as_LB) & ~isnan(as_UB);
|
|
N_valid = N(valid_idx);
|
|
LB_valid = as_LB(valid_idx);
|
|
UB_valid = as_UB(valid_idx);
|
|
|
|
% Create shaded area between LB and UB
|
|
x_fill = [N_valid, fliplr(N_valid)];
|
|
y_fill = [LB_valid, fliplr(UB_valid)];
|
|
|
|
% Plot settings
|
|
figure(1);
|
|
set(gcf,'Position',[100 100 950 750])
|
|
fill(x_fill, y_fill, [0.8 0.8 1], 'EdgeColor', 'none'); % Light blue shade
|
|
% Axes settings
|
|
set(gca, 'XScale', 'log');
|
|
xlim([1E4, 1E7]);
|
|
ylim([79, 106]);
|
|
xlabel('Atom number', 'Interpreter', 'latex', 'FontSize', 16);
|
|
ylabel('Scattering length $a_s$ ($a_0$)', 'Interpreter', 'latex', 'FontSize', 16);
|
|
grid on;
|
|
set(gca,'FontSize',16,'Box','On','Linewidth',2);
|
|
|
|
%% Visualize phase diagram
|
|
load('./Results/Data_3D/GradientDescent/phase_diagram_matrix_theta_30.mat')
|
|
PhaseDiagramMatrix = M;
|
|
ScatteringLengths = SCATTERING_LENGTH_RANGE;
|
|
NumberOfAtoms = NUM_ATOMS_LIST * 1E-5;
|
|
xlen = length(ScatteringLengths);
|
|
ylen = length(NumberOfAtoms);
|
|
|
|
figure(1);
|
|
set(gcf,'Position',[100 100 950 750])
|
|
set(gca,'FontSize',16,'Box','On','Linewidth',2);
|
|
hImg = imagesc(M);
|
|
set(gca, 'YDir', 'normal');
|
|
colormap(parula);
|
|
ax = gca; % Get current axes
|
|
ax.FontSize = 16; % Set tick label font size (adjust as needed)
|
|
axis equal tight;
|
|
xticks(1:xlen);
|
|
yticks(1:ylen);
|
|
xticklabels(string(NumberOfAtoms));
|
|
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;
|
|
|
|
%% Density modulation determination
|
|
|
|
% SaveDirectory = './Results/Data_3D/GradientDescent/Phi020/aS_090_theta_020_phi_000_N_100000';
|
|
% JobNumber = 0;
|
|
|
|
SaveDirectory = './Results/Data_3D/PhaseDiagram';
|
|
JobNumber = 1;
|
|
|
|
% 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
|