karthik/Calculations
1
0
Fork 0
Code Issues Pull Requests repo.project_board Releases Wiki Activity

Benchmarking complete - Code reliably reproduces expected ground states, script added to begin investigation in to tilt of the dipoles.

Browse Source
This commit is contained in:
Karthik 2025-01-23 20:50:06 +01:00
parent 5e7940a1b1
commit 59fa679549
7 changed files with 238 additions and 144 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

62
Dipolar-Gas-Simulator/+Scripts/run_locally.m
View File

@ -76,7 +76,7 @@ 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.TimeCutOff = 1E6; % in s
OptionsStruct.EnergyTolerance = 5E-10;
OptionsStruct.ResidualTolerance = 1E-05;
OptionsStruct.NoiseScaleFactor = 0.05;
@ -91,7 +91,7 @@ OptionsStruct.PlotLive = true;
OptionsStruct.JobNumber = 1;
OptionsStruct.RunOnGPU = false;
OptionsStruct.SaveData = true;
OptionsStruct.SaveDirectory = './Data_TriangularPhase';
OptionsStruct.SaveDirectory = './Results/Data_TriangularPhase';
options = Helper.convertstruct2cell(OptionsStruct);
clear OptionsStruct
@ -144,7 +144,7 @@ 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.TimeCutOff = 1E6; % in s
OptionsStruct.EnergyTolerance = 5E-10;
OptionsStruct.ResidualTolerance = 1E-05;
OptionsStruct.NoiseScaleFactor = 0.05;
@ -159,7 +159,7 @@ OptionsStruct.PlotLive = true;
OptionsStruct.JobNumber = 1;
OptionsStruct.RunOnGPU = false;
OptionsStruct.SaveData = true;
OptionsStruct.SaveDirectory = './Data_StripePhase';
OptionsStruct.SaveDirectory = './Results/Data_StripePhase';
options = Helper.convertstruct2cell(OptionsStruct);
clear OptionsStruct
@ -174,7 +174,7 @@ solver.Potential = pot.trap();
OptionsStruct = struct;
OptionsStruct.NumberOfAtoms = 4.148e+07;
OptionsStruct.NumberOfAtoms = 4.2e+07;
OptionsStruct.DipolarPolarAngle = 0;
OptionsStruct.DipolarAzimuthAngle = 0;
OptionsStruct.ScatteringLength = 101.35;
@ -201,7 +201,7 @@ OptionsStruct.PlotLive = true;
OptionsStruct.JobNumber = 1;
OptionsStruct.RunOnGPU = false;
OptionsStruct.SaveData = true;
OptionsStruct.SaveDirectory = './Data_HoneycombPhase';
OptionsStruct.SaveDirectory = './Results/Data_HoneycombPhase';
options = Helper.convertstruct2cell(OptionsStruct);
clear OptionsStruct
@ -222,17 +222,61 @@ Plotter.visualizeWavefunction2D(psi,Params,Transf)
Plotter.visualizeGSWavefunction2D(solver.SaveDirectory, solver.JobNumber)
%% - Analysis
SaveDirectory = './Data_TriangularPhase';
SaveDirectory = './Results/Data_TriangularPhase';
JobNumber = 1;
Plotter.visualizeGSWavefunction2D(SaveDirectory, JobNumber)
%% - Analysis
SaveDirectory = './Data_StripePhase';
SaveDirectory = './Results/Data_StripePhase';
JobNumber = 2;
Plotter.visualizeGSWavefunction2D(SaveDirectory, JobNumber)
%% - Analysis
SaveDirectory = './Data_HoneycombPhase';
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 = 0: a_s = 76.41
OptionsStruct = struct;
OptionsStruct.NumberOfAtoms = 1.00e+05;
OptionsStruct.DipolarPolarAngle = 0;
OptionsStruct.DipolarAzimuthAngle = 0;
OptionsStruct.ScatteringLength = 76.41;
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 = 2E6; % in s
OptionsStruct.EnergyTolerance = 5E-10;
OptionsStruct.ResidualTolerance = 1E-05;
OptionsStruct.NoiseScaleFactor = 0.05;
OptionsStruct.MaxIterations = 10;
OptionsStruct.VariationalWidth = 2;
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_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();

143
Dipolar-Gas-Simulator/+Scripts/run_on_cluster.m
View File

@ -1,29 +1,20 @@
%%
% 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
%% Tilting of the dipoles
% Atom Number = 1.00e+05
% System size = [10, 10]
%% v_z = 500, theta = 0: a_s = 76.41
OptionsStruct = struct;
OptionsStruct.NumberOfAtoms = 2.0030e+07;
OptionsStruct.NumberOfAtoms = 1.00e+05;
OptionsStruct.DipolarPolarAngle = 0;
OptionsStruct.DipolarAzimuthAngle = 0;
OptionsStruct.ScatteringLength = 98.0676;
OptionsStruct.ScatteringLength = 76.41;
OptionsStruct.TrapFrequencies = [0, 0, 72.4];
OptionsStruct.TrapFrequencies = [0, 0, 500];
OptionsStruct.TrapPotentialType = 'None';
OptionsStruct.NumberOfGridPoints = [128, 128];
OptionsStruct.Dimensions = [100, 100];
OptionsStruct.Dimensions = [10, 10];
OptionsStruct.TimeStepSize = 0.005; % in s
OptionsStruct.MinimumTimeStepSize = 1E-5; % in s
OptionsStruct.TimeCutOff = 2E6; % in s
@ -32,7 +23,7 @@ OptionsStruct.ResidualTolerance = 1E-05;
OptionsStruct.NoiseScaleFactor = 0.05;
OptionsStruct.MaxIterations = 10;
OptionsStruct.VariationalWidth = 5;
OptionsStruct.VariationalWidth = 2;
OptionsStruct.WidthLowerBound = 1;
OptionsStruct.WidthUpperBound = 12;
OptionsStruct.WidthCutoff = 5e-3;
@ -41,7 +32,7 @@ OptionsStruct.PlotLive = false;
OptionsStruct.JobNumber = 1;
OptionsStruct.RunOnGPU = true;
OptionsStruct.SaveData = true;
OptionsStruct.SaveDirectory = './Data_TriangularPhase';
OptionsStruct.SaveDirectory = './Results/Data_TiltingOfDipoles';
options = Helper.convertstruct2cell(OptionsStruct);
clear OptionsStruct
@ -51,117 +42,3 @@ 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 = 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 = 2E6; % 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 = false;
OptionsStruct.JobNumber = 2;
OptionsStruct.RunOnGPU = true;
OptionsStruct.SaveData = true;
OptionsStruct.SaveDirectory = './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.0102e+07;
OptionsStruct.DipolarPolarAngle = 0;
OptionsStruct.DipolarAzimuthAngle = 0;
OptionsStruct.ScatteringLength = 101.9903;
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 = false;
OptionsStruct.JobNumber = 3;
OptionsStruct.RunOnGPU = true;
OptionsStruct.SaveData = true;
OptionsStruct.SaveDirectory = './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();
%%
%{
% 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');
%}

166
Dipolar-Gas-Simulator/+Scripts/run_on_cluster_benchmarking.m Normal file
View File

@ -0,0 +1,166 @@
%%
% 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 = 2E6; % 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 = false;
OptionsStruct.JobNumber = 1;
OptionsStruct.RunOnGPU = true;
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();
%% - 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 = 2E6; % 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 = false;
OptionsStruct.JobNumber = 2;
OptionsStruct.RunOnGPU = true;
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 = false;
OptionsStruct.JobNumber = 3;
OptionsStruct.RunOnGPU = true;
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();
%%
%{
% 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');
%}

3
Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/propagateWavefunction.m
View File

@ -1,4 +1,7 @@
function [psi, Observ] = propagateWavefunction(this, psi, Params, VParams, Transf, VDk, V, t_idx, Observ, vrun)
format long
set(0,'defaulttextInterpreter','latex')
set(groot, 'defaultAxesTickLabelInterpreter','latex'); set(groot, 'defaultLegendInterpreter','latex');

1
Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/run.m
View File

@ -1,4 +1,5 @@
function [Params, Transf, psi, V, VDk] = run(this)
format long;
% --- Obtain simulation parameters ---

5
Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/setupParameters.m
View File

@ -1,4 +1,7 @@
function [Params] = setupParameters(this)
format long
CONSTANTS = Helper.PhysicsConstants;
hbar = CONSTANTS.PlanckConstantReduced; % [J.s]
kbol = CONSTANTS.BoltzmannConstant; % [J/K]
@ -6,7 +9,7 @@ function [Params] = setupParameters(this)
muB = CONSTANTS.BohrMagneton; % [J/T]
a0 = CONSTANTS.BohrRadius; % [m]
m0 = CONSTANTS.AtomicMassUnit; % [kg]
w0 = 2*pi*61.6316; % Angular frequency unit [s^-1]
w0 = 2*pi*61.658214297935530; % Angular frequency unit [s^-1]
mu0factor = 0.3049584233607396; % =(m0/me)*pi*alpha^2 -- me=mass of electron, alpha=fine struct. const.
% mu0=mu0factor *hbar^2*a0/(m0*muB^2)
% Number of points in each direction

2
Dipolar-Gas-Simulator/bwhpc_matlab_gpe_sim_gpu.slurm
View File

@ -6,7 +6,7 @@
#SBATCH --nodes=1
#SBATCH --ntasks-per-node=1
#SBATCH --gres=gpu:A40:1
#SBATCH --mem=8G
#SBATCH --mem=2G
# Estimated wallclock time for job
#SBATCH --time=03:00:00
#SBATCH --job-name=simulation