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

%% 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           = 1E5;
OptionsStruct.DipolarPolarAngle       = 0;
OptionsStruct.DipolarAzimuthAngle     = 0;
OptionsStruct.ScatteringLength        = 86;

OptionsStruct.TrapFrequencies         = [10, 10, 72.4];
OptionsStruct.TrapDepth               = 5;
OptionsStruct.BoxSize                 = 15;
OptionsStruct.TrapPotentialType       = 'Harmonic';

OptionsStruct.NumberOfGridPoints      = [256, 512, 256];
OptionsStruct.Dimensions              = [50, 120, 150];
OptionsStruct.CutoffType              = 'Cylindrical';
OptionsStruct.SimulationMode          = 'ImaginaryTimeEvolution'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution'
OptionsStruct.TimeStepSize            = 0.005;                    % in s
OptionsStruct.NumberOfTimeSteps       = 2E6;                      % in s
OptionsStruct.EnergyTolerance         = 5E-10;
OptionsStruct.ResidualTolerance       = 1E-05;

OptionsStruct.JobNumber               = 1;
OptionsStruct.RunOnGPU                = false;
OptionsStruct.SaveData                = true;
OptionsStruct.SaveDirectory           = './Data';
options                               = Helper.convertstruct2cell(OptionsStruct);
clear OptionsStruct

sim                                   = Simulator.DipolarGas(options{:});
pot                                   = Simulator.Potentials(options{:});
sim.Potential                         = pot.trap(); % + pot.repulsive_chopstick();

%-% Run Simulation %-%
[Params, Transf, psi, V, VDk]         = sim.run();

%% - Plot numerical grid
% Plotter.visualizeSpace(Transf)
%% - Plot trap potential
Plotter.visualizeTrapPotential(sim.Potential,Params,Transf)
%% - Plot initial wavefunction
Plotter.visualizeWavefunction(psi,Params,Transf)
%% - Plot GS wavefunction
Plotter.visualizeGSWavefunction(Params.njob)

%%

% 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]

%% v_z = 500, theta = 45: a_s = 85.17

OptionsStruct = struct;

OptionsStruct.NumberOfAtoms           = 1.00e+05;     
OptionsStruct.DipolarPolarAngle       = deg2rad(45);
OptionsStruct.DipolarAzimuthAngle     = 0;
OptionsStruct.ScatteringLength        = 75.00;               

OptionsStruct.TrapFrequencies         = [0, 0, 500];
OptionsStruct.TrapPotentialType       = 'None';

OptionsStruct.NumberOfGridPoints      = [128, 128];
OptionsStruct.Dimensions              = [10, 10];      
OptionsStruct.TimeStepSize            = 0.005;           % in s
OptionsStruct.MinimumTimeStepSize     = 1E-5;            % in s
OptionsStruct.TimeCutOff              = 5E5;             % in s
OptionsStruct.EnergyTolerance         = 5E-10;
OptionsStruct.ResidualTolerance       = 1E-05;
OptionsStruct.NoiseScaleFactor        = 0.15;

OptionsStruct.MaxIterations           = 10;      
OptionsStruct.VariationalWidth        = 1;       
OptionsStruct.WidthLowerBound         = 0.01;
OptionsStruct.WidthUpperBound         = 12;
OptionsStruct.WidthCutoff             = 5e-3;

OptionsStruct.PlotLive                = true;
OptionsStruct.JobNumber               = 4;
OptionsStruct.RunOnGPU                = false;
OptionsStruct.SaveData                = true;
OptionsStruct.SaveDirectory           = './Results/Data_TiltingOfDipoles';
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();

%% v_z = 500, theta = 65: a_s = 93.41

OptionsStruct = struct;

OptionsStruct.NumberOfAtoms           = 1.00e+05;     
OptionsStruct.DipolarPolarAngle       = deg2rad(65);
OptionsStruct.DipolarAzimuthAngle     = 0;
OptionsStruct.ScatteringLength        = 81.00;        

OptionsStruct.TrapFrequencies         = [0, 0, 500];
OptionsStruct.TrapPotentialType       = 'None';

OptionsStruct.NumberOfGridPoints      = [128, 128];
OptionsStruct.Dimensions              = [10, 10];      
OptionsStruct.TimeStepSize            = 0.0025;          % in s
OptionsStruct.MinimumTimeStepSize     = 1E-5;            % in s
OptionsStruct.TimeCutOff              = 5E5;             % in s
OptionsStruct.EnergyTolerance         = 5E-10;
OptionsStruct.ResidualTolerance       = 1E-05;
OptionsStruct.NoiseScaleFactor        = 0.15;

OptionsStruct.MaxIterations           = 10;      
OptionsStruct.VariationalWidth        = 1;       
OptionsStruct.WidthLowerBound         = 0.01;
OptionsStruct.WidthUpperBound         = 12;
OptionsStruct.WidthCutoff             = 5e-3;

OptionsStruct.PlotLive                = true;
OptionsStruct.JobNumber               = 5;
OptionsStruct.RunOnGPU                = false;
OptionsStruct.SaveData                = true;
OptionsStruct.SaveDirectory           = './Results/Data_TiltingOfDipoles';
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();

%% v_z = 500, theta = 75: a_s = 98.11

OptionsStruct = struct;

OptionsStruct.NumberOfAtoms           = 1.00e+05;     
OptionsStruct.DipolarPolarAngle       = deg2rad(75);
OptionsStruct.DipolarAzimuthAngle     = 0;
OptionsStruct.ScatteringLength        = 85.00;        

OptionsStruct.TrapFrequencies         = [0, 0, 500];
OptionsStruct.TrapPotentialType       = 'None';

OptionsStruct.NumberOfGridPoints      = [128, 128];
OptionsStruct.Dimensions              = [10, 10];      
OptionsStruct.TimeStepSize            = 0.001;          % in s
OptionsStruct.MinimumTimeStepSize     = 1E-5;            % in s
OptionsStruct.TimeCutOff              = 5E5;             % in s
OptionsStruct.EnergyTolerance         = 5E-10;
OptionsStruct.ResidualTolerance       = 1E-05;
OptionsStruct.NoiseScaleFactor        = 0.05;

OptionsStruct.MaxIterations           = 10;      
OptionsStruct.VariationalWidth        = 1;       
OptionsStruct.WidthLowerBound         = 0.01;
OptionsStruct.WidthUpperBound         = 12;
OptionsStruct.WidthCutoff             = 5e-3;

OptionsStruct.PlotLive                = true;
OptionsStruct.JobNumber               = 6;
OptionsStruct.RunOnGPU                = false;
OptionsStruct.SaveData                = true;
OptionsStruct.SaveDirectory           = './Results/Data_TiltingOfDipoles';
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();

%% v_z = 500, theta = 85: a_s = 102.56

OptionsStruct = struct;

OptionsStruct.NumberOfAtoms           = 1.00e+05;     
OptionsStruct.DipolarPolarAngle       = deg2rad(85);
OptionsStruct.DipolarAzimuthAngle     = 0;
OptionsStruct.ScatteringLength        = 90.00;        

OptionsStruct.TrapFrequencies         = [0, 0, 500];
OptionsStruct.TrapPotentialType       = 'None';

OptionsStruct.NumberOfGridPoints      = [128, 128];
OptionsStruct.Dimensions              = [10, 10];      
OptionsStruct.TimeStepSize            = 0.001;          % in s
OptionsStruct.MinimumTimeStepSize     = 1E-5;            % in s
OptionsStruct.TimeCutOff              = 5E5;             % in s
OptionsStruct.EnergyTolerance         = 5E-10;
OptionsStruct.ResidualTolerance       = 1E-05;
OptionsStruct.NoiseScaleFactor        = 0.05;

OptionsStruct.MaxIterations           = 10;      
OptionsStruct.VariationalWidth        = 1;       
OptionsStruct.WidthLowerBound         = 0.01;
OptionsStruct.WidthUpperBound         = 12;
OptionsStruct.WidthCutoff             = 5e-3;

OptionsStruct.PlotLive                = true;
OptionsStruct.JobNumber               = 7;
OptionsStruct.RunOnGPU                = false;
OptionsStruct.SaveData                = true;
OptionsStruct.SaveDirectory           = './Results/Data_TiltingOfDipoles';
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();

%% - Analysis
SaveDirectory = './Results/Data_TiltingOfDipoles';
JobNumber     = 4;
Plotter.visualizeGSWavefunction2D(SaveDirectory, JobNumber)