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

Latest working version - modified calculation of VDk to have the flag to include cutoff within it so that there is a single universal function that calculated the DDI potential for the solver.

Browse Source
This commit is contained in:
Karthik 2025-03-19 14:56:26 +01:00
parent 4af0ce726a
commit f6394cd6a0
8 changed files with 124 additions and 137 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -326,11 +326,19 @@ JobNumber = 1;
[contrast, period_X, period_Y] = Scripts.analyzeGSWavefunction_in_plane_trap(SaveDirectory, JobNumber);
%%
SaveDirectory = './Results/Data_TiltingOfDipoles/TransitionAngle/OptimalSystemSize/Hz500/Degree15';
SaveDirectory = './Results/Data_TiltingOfDipoles/TransitionAngle/OptimalSystemSize/Hz500/Degree0';
% Define the desired range for Lx and Ly
Lx_Range = 4:0.4:11; % Extend Lx from 5 to 10
Ly_Range = 4:0.4:11; % Extend Ly from 5 to 10
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';
% Define the desired range
LatticeSpacing = 1.0:0.05:4.0;
MinEnergyDataArray = Scripts.analyzeGSWavefunction_constrained_optimal_system_size(SaveDirectory, LatticeSpacing);
% Plotter.visualizeGSWavefunction2D(SaveDirectory, 1)
%% - Analysis
SaveDirectory = './Results/Data_TiltingOfDipoles/HarmonicTrap/AspectRatio/AR2_8';
JobNumber = 0; % 79
@ -338,6 +346,7 @@ JobNumber = 0; % 79
% 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
@ -346,3 +355,4 @@ JobNumber = 2; % 82
% JobNumber = 3; % 83
Plotter.visualizeGSWavefunction2D(SaveDirectory, JobNumber)
[contrast, period_X, period_Y] = Scripts.analyzeGSWavefunction_in_plane_trap(SaveDirectory, JobNumber);

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

@ -1,89 +1,74 @@
%% AR = 2.8
%% Tilting of the dipoles
OptionsStruct = struct;
ppum = 1250; % Atom Number Density in per micrometers
OptionsStruct.NumberOfAtoms = 5E5;
OptionsStruct.DipolarPolarAngle = 0;
OptionsStruct.DipolarAzimuthAngle = 0;
OptionsStruct.ScatteringLength = 78.0;
%% theta = 0
AspectRatio = 2.8;
HorizontalTrapFrequency = 125;
VerticalTrapFrequency = HorizontalTrapFrequency * AspectRatio;
OptionsStruct.TrapFrequencies = [HorizontalTrapFrequency, HorizontalTrapFrequency, VerticalTrapFrequency];
OptionsStruct.TrapPotentialType = 'Harmonic';
LatticeSpacing = 1.0:0.05:4.0;
totalIterations = numel(LatticeSpacing);
OptionsStruct.NumberOfGridPoints = [256, 256];
OptionsStruct.Dimensions = [18, 18];
OptionsStruct.TimeStepSize = 1E-3; % in s
OptionsStruct.MinimumTimeStepSize = 1E-5; % in s
OptionsStruct.TimeCutOff = 2E6; % in s
OptionsStruct.EnergyTolerance = 5E-10;
OptionsStruct.ResidualTolerance = 1E-04;
OptionsStruct.NoiseScaleFactor = 0.05;
% create a local cluster object
cluster = parcluster('local');
OptionsStruct.MaxIterations = 10;
OptionsStruct.VariationalWidth = 2.0;
OptionsStruct.WidthLowerBound = 0.01;
OptionsStruct.WidthUpperBound = 12;
OptionsStruct.WidthCutoff = 1e-2;
% get the number of dedicated cores from environment
nprocs = str2num(getenv('SLURM_NPROCS'));
OptionsStruct.PlotLive = false;
OptionsStruct.JobNumber = 0;
OptionsStruct.RunOnGPU = true;
OptionsStruct.SaveData = true;
OptionsStruct.SaveDirectory = './Results/Data_TiltingOfDipoles/HarmonicTrap/AspectRatio/AR2_8';
options = Helper.convertstruct2cell(OptionsStruct);
clear OptionsStruct
% you may explicitly set the JobStorageLocation to the tmp directory that is unique to each cluster job (and is on local, fast scratch)
parpool_tmpdir = [getenv('TMP'),'/.matlab/local_cluster_jobs/slurm_jobID_',getenv('SLURM_JOB_ID')];
mkdir(parpool_tmpdir);
cluster.JobStorageLocation = parpool_tmpdir;
solver = VariationalSolver2D.DipolarGas(options{:});
pot = VariationalSolver2D.Potentials(options{:});
solver.Potential = pot.trap();
% start the parallel pool
parpool(cluster, nprocs)
%-% Run Solver %-%
[Params, Transf, psi, V, VDk] = solver.run();
% Parallel loop over all combinations of i, j
parfor (k = 1:totalIterations, cluster)
%% AR = 4.0
a = LatticeSpacing(k);
lx = a;
ly = sqrt(3)*a;
OptionsStruct = struct;
% Initialize OptionsStruct
OptionsStruct = struct;
OptionsStruct.NumberOfAtoms = 5E5;
OptionsStruct.DipolarPolarAngle = 0;
OptionsStruct.DipolarAzimuthAngle = 0;
OptionsStruct.ScatteringLength = 78.0;
% Assign values to OptionsStruct
OptionsStruct.NumberOfAtoms = ppum * (lx*ly);
OptionsStruct.DipolarPolarAngle = deg2rad(0);
OptionsStruct.DipolarAzimuthAngle = 0;
OptionsStruct.ScatteringLength = 75.00;
AspectRatio = 4.0;
HorizontalTrapFrequency = 125;
VerticalTrapFrequency = HorizontalTrapFrequency * AspectRatio;
OptionsStruct.TrapFrequencies = [HorizontalTrapFrequency, HorizontalTrapFrequency, VerticalTrapFrequency];
OptionsStruct.TrapPotentialType = 'Harmonic';
OptionsStruct.TrapFrequencies = [0, 0, 500];
OptionsStruct.TrapPotentialType = 'None';
OptionsStruct.NumberOfGridPoints = [256, 256];
OptionsStruct.Dimensions = [18, 18];
OptionsStruct.TimeStepSize = 1E-3; % in s
OptionsStruct.MinimumTimeStepSize = 1E-5; % in s
OptionsStruct.TimeCutOff = 2E6; % in s
OptionsStruct.EnergyTolerance = 5E-10;
OptionsStruct.ResidualTolerance = 1E-04;
OptionsStruct.NoiseScaleFactor = 0.05;
OptionsStruct.NumberOfGridPoints = [128, 128];
OptionsStruct.Dimensions = [lx, ly];
OptionsStruct.TimeStepSize = 5E-4; % 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.IncludeDDICutOff = false;
OptionsStruct.MaxIterations = 10;
OptionsStruct.VariationalWidth = 2.0;
OptionsStruct.WidthLowerBound = 0.01;
OptionsStruct.WidthUpperBound = 12;
OptionsStruct.WidthCutoff = 1e-2;
OptionsStruct.MaxIterations = 10;
OptionsStruct.VariationalWidth = 1.10;
OptionsStruct.WidthLowerBound = 0.01;
OptionsStruct.WidthUpperBound = 12;
OptionsStruct.WidthCutoff = 1e-2;
OptionsStruct.PlotLive = false;
OptionsStruct.JobNumber = 0;
OptionsStruct.RunOnGPU = true;
OptionsStruct.SaveData = true;
OptionsStruct.SaveDirectory = './Results/Data_TiltingOfDipoles/HarmonicTrap/AspectRatio/AR4_0';
options = Helper.convertstruct2cell(OptionsStruct);
clear OptionsStruct
OptionsStruct.PlotLive = false;
OptionsStruct.JobNumber = k;
OptionsStruct.RunOnGPU = true;
OptionsStruct.SaveData = true;
OptionsStruct.SaveDirectory = sprintf('./Results/Data_TiltingOfDipoles/TransitionAngle/OptimalSystemSize/Hz500/Degree%i', round(rad2deg(OptionsStruct.DipolarPolarAngle)));
solver = VariationalSolver2D.DipolarGas(options{:});
pot = VariationalSolver2D.Potentials(options{:});
solver.Potential = pot.trap();
options = Helper.convertstruct2cell(OptionsStruct);
%-% Run Solver %-%
[Params, Transf, psi, V, VDk] = solver.run();
solver = VariationalSolver2D.DipolarGas(options{:});
pot = VariationalSolver2D.Potentials(options{:});
solver.Potential = pot.trap();
% Run Solver
[Params, Transf, psi, V, VDk] = solver.run();
end

53
Dipolar-Gas-Simulator/+Scripts/run_on_cluster_optimal_system_size.m
View File

@ -1,31 +1,43 @@
%% Tilting of the dipoles
% Atom Number Density = 1250 ppum
ppum = 1250;
ppum = 1250; % Atom Number Density in per micrometers
%% theta = 0
Lx = 8:0.4:11;
Ly = 4:0.4:11;
cluster = parcluster;
Lx = 4:0.4:11;
Ly = 4:0.4:11;
% Create a combined index for all i, j pairs
[I, J] = ndgrid(1:numel(Lx), 1:numel(Ly));
totalIterations = numel(I);
[I, J] = ndgrid(1:numel(Lx), 1:numel(Ly));
totalIterations = numel(I);
% create a local cluster object
cluster = parcluster('local');
% get the number of dedicated cores from environment
nprocs = str2num(getenv('SLURM_NPROCS'));
% you may explicitly set the JobStorageLocation to the tmp directory that is unique to each cluster job (and is on local, fast scratch)
parpool_tmpdir = [getenv('TMP'),'/.matlab/local_cluster_jobs/slurm_jobID_',getenv('SLURM_JOB_ID')];
mkdir(parpool_tmpdir);
cluster.JobStorageLocation = parpool_tmpdir;
% start the parallel pool
parpool(cluster, nprocs)
% Parallel loop over all combinations of i, j
parfor (k = 1:totalIterations, cluster)
i = I(k);
j = J(k);
lx = Lx(i);
ly = Ly(j);
i = I(k);
j = J(k);
lx = Lx(i);
ly = Ly(j);
% Initialize OptionsStruct
OptionsStruct = struct;
OptionsStruct = struct;
% Assign values to OptionsStruct
OptionsStruct.NumberOfAtoms = ppum * (lx*ly);
OptionsStruct.DipolarPolarAngle = deg2rad(15);
OptionsStruct.DipolarPolarAngle = deg2rad(0);
OptionsStruct.DipolarAzimuthAngle = 0;
OptionsStruct.ScatteringLength = 75.00;
@ -49,17 +61,18 @@ parfor (k = 1:totalIterations, cluster)
OptionsStruct.WidthCutoff = 1e-2;
OptionsStruct.PlotLive = false;
OptionsStruct.JobNumber = (i - 1) * num_inner + j;
OptionsStruct.JobNumber = k;
OptionsStruct.RunOnGPU = true;
OptionsStruct.SaveData = true;
OptionsStruct.SaveDirectory = sprintf('./Results/Data_TiltingOfDipoles/TransitionAngle/OptimalSystemSize/Hz500/Degree%i', round(rad2deg(OptionsStruct.DipolarPolarAngle)));
options = Helper.convertstruct2cell(OptionsStruct);
options = Helper.convertstruct2cell(OptionsStruct);
solver = VariationalSolver2D.DipolarGas(options{:});
pot = VariationalSolver2D.Potentials(options{:});
solver.Potential = pot.trap();
solver = VariationalSolver2D.DipolarGas(options{:});
pot = VariationalSolver2D.Potentials(options{:});
solver.Potential = pot.trap();
% Run Solver
[Params, Transf, psi, V, VDk] = solver.run();
[Params, Transf, psi, V, VDk] = solver.run();
end

20
Dipolar-Gas-Simulator/+VariationalSolver2D/@Calculator/calculateVDkWithCutoff.m → Dipolar-Gas-Simulator/+VariationalSolver2D/@Calculator/calculateVDk.m
View File

@ -1,4 +1,4 @@
function VDk = calculateVDkWithCutoff(~, Transf, Params, ell)
function VDk = calculateVDk(~, Transf, Params, ell, IncludeDDICutOff)
% == Calculating the DDI potential in Fourier space with appropriate cutoff == %
% Interaction in K space
@ -17,13 +17,15 @@ function VDk = calculateVDkWithCutoff(~, Transf, Params, ell)
% Full DDI
VDk = Fpar*sin(Params.theta)^2 + Fperp*cos(Params.theta)^2;
% Implementing a cutoff:
VDr = ifftn(VDk);
VDr = fftshift(VDr);
rcut = min(Transf.Xmax,Transf.Ymax);
R = sqrt(Transf.X.^2+Transf.Y.^2) < 0.9*rcut;
VDr = VDr.*double(R);
VDr = ifftshift(VDr);
VDk = fftn(VDr);
if IncludeDDICutOff
% Implementing a cutoff:
VDr = ifftn(VDk);
VDr = fftshift(VDr);
rcut = min(Transf.Xmax,Transf.Ymax);
R = sqrt(Transf.X.^2+Transf.Y.^2) < 0.9*rcut;
VDr = VDr.*double(R);
VDr = ifftshift(VDr);
VDk = fftn(VDr);
end
end

19
Dipolar-Gas-Simulator/+VariationalSolver2D/@Calculator/calculateVDkWithoutCutoff.m
View File

@ -1,19 +0,0 @@
function VDk = calculateVDkWithoutCutoff(~, Transf, Params, ell)
% == Calculating the DDI potential in Fourier space with appropriate cutoff == %
% Interaction in K space
QX = Transf.KX*ell/sqrt(2);
QY = Transf.KY*ell/sqrt(2);
Qsq = QX.^2 + QY.^2;
absQ = sqrt(Qsq);
QDsq = QX.^2*cos(Params.eta)^2 + QY.^2*sin(Params.eta)^2;
% Bare interaction
Fpar = -1 + 3*sqrt(pi)*QDsq.*erfcx(absQ)./absQ; % Scaled complementary error function erfcx(x) = e^(x^2) * erfc(x)
Fperp = 2 - 3*sqrt(pi).*absQ.*erfcx(absQ);
Fpar(absQ == 0) = -1;
% Full DDI
VDk = Fpar*sin(Params.theta)^2 + Fperp*cos(Params.theta)^2;
end

4
Dipolar-Gas-Simulator/+VariationalSolver2D/@Calculator/calculateVariationalEnergy.m
View File

@ -1,4 +1,4 @@
function E = calculateVariationalEnergy(this, psi, Params, ell, Transf, V)
function E = calculateVariationalEnergy(this, psi, Params, ell, Transf, V, IncludeDDICutOff)
g_eff = Params.gs * (1/(sqrt(2*pi)*ell));
gamma_eff = Params.gammaQF * (sqrt(2/5)/(pi^(3/4)*ell^(3/2)));
@ -9,7 +9,7 @@ function E = calculateVariationalEnergy(this, psi, Params, ell, Transf, V)
normfac = Params.Lx*Params.Ly/numel(psi);
% DDIs
VDk = this.calculateVDkWithCutoff(Transf, Params, ell); % VDk depends on the variational parameters. THIS HAS TO BE RECALCULATED HERE FOR THE CONSTRAINED OPTIMIZATION TO WORK!
VDk = this.calculateVDk(Transf, Params, ell, IncludeDDICutOff); % VDk depends on the variational parameters. THIS HAS TO BE RECALCULATED HERE FOR THE CONSTRAINED OPTIMIZATION TO WORK!
frho = fftn(abs(psi).^2);
Phi = real(ifftn(frho.*VDk));
Eddi = 0.5*Params.gdd*Phi.*abs(psi).^2/(sqrt(2*pi)*ell);

6
Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/initialize.m
View File

@ -4,11 +4,7 @@ function [psi,V,VDk] = initialize(this,Params,VParams,Transf)
assert(~anynan(V), 'Potential not defined! Specify as <SimulatorObject>.Potential = <PotentialsObject>.trap() + <AdditionalTerms>.');
% == Calculating the DDIs == %
if this.IncludeDDICutOff
VDk = this.Calculator.calculateVDkWithCutoff(Transf, Params, VParams.ell);
else
VDk = this.Calculator.calculateVDkWithoutCutoff(Transf, Params, VParams.ell);
end
VDk = this.Calculator.calculateVDk(Transf, Params, VParams.ell, this.IncludeDDICutOff);
% == Setting up the initial wavefunction == %
psi = this.setupWavefunction(Params,Transf);

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

@ -28,7 +28,7 @@ function [Params, Transf, psi, V, VDk] = run(this)
[psi,V,VDk] = this.initialize(Params,VParams,Transf);
ells(1) = VParams.ell;
E_Var = @(x) this.Calculator.calculateVariationalEnergy(psi, Params, x, Transf, V)/Params.N;
E_Var = @(x) this.Calculator.calculateVariationalEnergy(psi, Params, x, Transf, V, this.IncludeDDICutOff)/Params.N;
E_vs_iter(1) = E_Var(ells(1));
t_idx = 1; firstpoint = 1;
@ -61,7 +61,7 @@ function [Params, Transf, psi, V, VDk] = run(this)
psi = gather(psi);
% --- Constrained minimization ---
E_Var = @(x) this.Calculator.calculateVariationalEnergy(psi, Params, x, Transf, V)/Params.N;
E_Var = @(x) this.Calculator.calculateVariationalEnergy(psi, Params, x, Transf, V, this.IncludeDDICutOff)/Params.N;
VParams.ell = fmincon(E_Var,VParams.ell,[],[],[],[],Params.ell_lower,Params.ell_upper,[],fminconoptions);
% --- Convergence check ---