diff --git a/Dipolar-Gas-Simulator/+Scripts/run_locally.m b/Dipolar-Gas-Simulator/+Scripts/run_locally.m index 32a7e59..1c9acfc 100644 --- a/Dipolar-Gas-Simulator/+Scripts/run_locally.m +++ b/Dipolar-Gas-Simulator/+Scripts/run_locally.m @@ -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 @@ -345,4 +354,5 @@ SaveDirectory = './Results/Data_TiltingOfDipoles/HarmonicTrap/AspectRatio/AR3_7' JobNumber = 2; % 82 % JobNumber = 3; % 83 Plotter.visualizeGSWavefunction2D(SaveDirectory, JobNumber) -[contrast, period_X, period_Y] = Scripts.analyzeGSWavefunction_in_plane_trap(SaveDirectory, JobNumber); \ No newline at end of file +[contrast, period_X, period_Y] = Scripts.analyzeGSWavefunction_in_plane_trap(SaveDirectory, JobNumber); + diff --git a/Dipolar-Gas-Simulator/+Scripts/run_on_cluster.m b/Dipolar-Gas-Simulator/+Scripts/run_on_cluster.m index 558dc31..0477329 100644 --- a/Dipolar-Gas-Simulator/+Scripts/run_on_cluster.m +++ b/Dipolar-Gas-Simulator/+Scripts/run_on_cluster.m @@ -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) + + a = LatticeSpacing(k); + lx = a; + ly = sqrt(3)*a; + + % Initialize OptionsStruct + OptionsStruct = struct; -%% AR = 4.0 + % Assign values to OptionsStruct + OptionsStruct.NumberOfAtoms = ppum * (lx*ly); + OptionsStruct.DipolarPolarAngle = deg2rad(0); + OptionsStruct.DipolarAzimuthAngle = 0; + OptionsStruct.ScatteringLength = 75.00; -OptionsStruct = struct; + OptionsStruct.TrapFrequencies = [0, 0, 500]; + OptionsStruct.TrapPotentialType = 'None'; -OptionsStruct.NumberOfAtoms = 5E5; -OptionsStruct.DipolarPolarAngle = 0; -OptionsStruct.DipolarAzimuthAngle = 0; -OptionsStruct.ScatteringLength = 78.0; + 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 = 1.10; + OptionsStruct.WidthLowerBound = 0.01; + OptionsStruct.WidthUpperBound = 12; + OptionsStruct.WidthCutoff = 1e-2; -AspectRatio = 4.0; -HorizontalTrapFrequency = 125; -VerticalTrapFrequency = HorizontalTrapFrequency * AspectRatio; -OptionsStruct.TrapFrequencies = [HorizontalTrapFrequency, HorizontalTrapFrequency, VerticalTrapFrequency]; -OptionsStruct.TrapPotentialType = 'Harmonic'; + 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))); -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; + options = Helper.convertstruct2cell(OptionsStruct); + + solver = VariationalSolver2D.DipolarGas(options{:}); + pot = VariationalSolver2D.Potentials(options{:}); + solver.Potential = pot.trap(); -OptionsStruct.MaxIterations = 10; -OptionsStruct.VariationalWidth = 2.0; -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 - -solver = VariationalSolver2D.DipolarGas(options{:}); -pot = VariationalSolver2D.Potentials(options{:}); -solver.Potential = pot.trap(); - -%-% Run Solver %-% -[Params, Transf, psi, V, VDk] = solver.run(); \ No newline at end of file + % Run Solver + [Params, Transf, psi, V, VDk] = solver.run(); + +end \ No newline at end of file diff --git a/Dipolar-Gas-Simulator/+Scripts/run_on_cluster_optimal_system_size.m b/Dipolar-Gas-Simulator/+Scripts/run_on_cluster_optimal_system_size.m index a0ea3e2..e68b4f5 100644 --- a/Dipolar-Gas-Simulator/+Scripts/run_on_cluster_optimal_system_size.m +++ b/Dipolar-Gas-Simulator/+Scripts/run_on_cluster_optimal_system_size.m @@ -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(); -end + [Params, Transf, psi, V, VDk] = solver.run(); + +end \ No newline at end of file diff --git a/Dipolar-Gas-Simulator/+VariationalSolver2D/@Calculator/calculateVDkWithCutoff.m b/Dipolar-Gas-Simulator/+VariationalSolver2D/@Calculator/calculateVDk.m similarity index 60% rename from Dipolar-Gas-Simulator/+VariationalSolver2D/@Calculator/calculateVDkWithCutoff.m rename to Dipolar-Gas-Simulator/+VariationalSolver2D/@Calculator/calculateVDk.m index 7c49361..386fc76 100644 --- a/Dipolar-Gas-Simulator/+VariationalSolver2D/@Calculator/calculateVDkWithCutoff.m +++ b/Dipolar-Gas-Simulator/+VariationalSolver2D/@Calculator/calculateVDk.m @@ -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 \ No newline at end of file diff --git a/Dipolar-Gas-Simulator/+VariationalSolver2D/@Calculator/calculateVDkWithoutCutoff.m b/Dipolar-Gas-Simulator/+VariationalSolver2D/@Calculator/calculateVDkWithoutCutoff.m deleted file mode 100644 index 4f15d39..0000000 --- a/Dipolar-Gas-Simulator/+VariationalSolver2D/@Calculator/calculateVDkWithoutCutoff.m +++ /dev/null @@ -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 \ No newline at end of file diff --git a/Dipolar-Gas-Simulator/+VariationalSolver2D/@Calculator/calculateVariationalEnergy.m b/Dipolar-Gas-Simulator/+VariationalSolver2D/@Calculator/calculateVariationalEnergy.m index 4a780f3..66e3bec 100644 --- a/Dipolar-Gas-Simulator/+VariationalSolver2D/@Calculator/calculateVariationalEnergy.m +++ b/Dipolar-Gas-Simulator/+VariationalSolver2D/@Calculator/calculateVariationalEnergy.m @@ -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); diff --git a/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/initialize.m b/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/initialize.m index 6e63bee..fc75799 100644 --- a/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/initialize.m +++ b/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/initialize.m @@ -4,11 +4,7 @@ function [psi,V,VDk] = initialize(this,Params,VParams,Transf) assert(~anynan(V), 'Potential not defined! Specify as .Potential = .trap() + .'); % == 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); diff --git a/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/run.m b/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/run.m index f8e7c57..a26f588 100644 --- a/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/run.m +++ b/Dipolar-Gas-Simulator/+VariationalSolver2D/@DipolarGas/run.m @@ -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 ---