618 lines
26 KiB
Matlab
618 lines
26 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/GradientDescent/Phi050/aS_089_theta_050_phi_000_N_100000';
|
|
JobNumber = 0;
|
|
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
|
|
%%
|
|
% Parameters you can set before the loop
|
|
N = 105000;
|
|
theta = 50.0;
|
|
phi = 0;
|
|
JobNumber = 0;
|
|
|
|
% Values of aS to loop over
|
|
aS_values = [79.0 80.0 81.0 82.0 83.0 84.0 85.0 86.0 87.0 88.0 89.0 90.0];
|
|
|
|
for aS = aS_values
|
|
% Format numbers
|
|
aS_str = sprintf('%03.0f', aS);
|
|
theta_str = sprintf('%03.0f', theta);
|
|
phi_str = sprintf('%03.0f', phi);
|
|
N_str = sprintf('%d', N);
|
|
|
|
% Construct directory path
|
|
SaveDirectory = ['./Results/Data_3D/GradientDescent/Phi',theta_str, '/aS_', aS_str, ...
|
|
'_theta_', theta_str, '_phi_', phi_str, '_N_', N_str];
|
|
|
|
% Call the plotting function
|
|
try
|
|
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
|
|
catch ME
|
|
warning(ME.message)
|
|
continue;
|
|
end
|
|
|
|
% Pause to inspect plot before continuing
|
|
disp(['Plotted for aS = ', num2str(aS)])
|
|
pause(0.25)
|
|
end
|
|
|
|
|
|
%%
|
|
SaveDirectory = './Results/Data_3D/GradientDescent/';
|
|
JobNumber = 1;
|
|
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
|
|
%% Visualize phase diagram
|
|
load('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; |