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

Latest updated version of the full 3D simulator.

Browse Source
This commit is contained in:
Karthik 2025-03-26 17:20:28 +01:00
parent 95781545d4
commit 488c6c9f0b
19 changed files with 462 additions and 552 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

74
Dipolar-Gas-Simulator/+Plotter/plotLive.m
View File

@ -3,42 +3,74 @@ function plotLive(psi,Params,Transf,Observ)
set(groot, 'defaultAxesTickLabelInterpreter','latex'); set(groot, 'defaultLegendInterpreter','latex');
format long
x = Transf.x*Params.l0*1e6;
y = Transf.y*Params.l0*1e6;
z = Transf.z*Params.l0*1e6;
% Axes scaling and coordinates in micrometers
x = Transf.x * Params.l0 * 1e6;
y = Transf.y * Params.l0 * 1e6;
z = Transf.z * Params.l0 * 1e6;
dx = x(2)-x(1); dy = y(2)-y(1); dz = z(2)-z(1);
%Plotting
subplot(2,3,1)
% Compute probability density |psi|^2
n = abs(psi).^2;
%Plotting
figure(1);
clf
set(gcf,'Position', [100, 100, 1600, 900])
t = tiledlayout(2, 3, 'TileSpacing', 'compact', 'Padding', 'compact'); % 2x3 grid
nexttile;
nxz = squeeze(trapz(n*dy,2));
nyz = squeeze(trapz(n*dx,1));
nxy = squeeze(trapz(n*dz,3));
plotxz = pcolor(x,z,nxz');
set(plotxz, 'EdgeColor', 'none');
xlabel('$x$ [$\mu$m]'); ylabel('$z$ [$\mu$m]');
subplot(2,3,2)
cbar1 = colorbar;
cbar1.Label.Interpreter = 'latex';
% cbar1.Ticks = []; % Disable the ticks
colormap(gca, Helper.Colormaps.plasma())
xlabel('$x$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14)
ylabel('$z$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14)
title('$|\Psi(x,y)|^2$', 'Interpreter', 'latex', 'FontSize', 14)
nexttile;
plotyz = pcolor(y,z,nyz');
set(plotyz, 'EdgeColor', 'none');
xlabel('$y$ [$\mu$m]'); ylabel('$z$ [$\mu$m]');
cbar1 = colorbar;
cbar1.Label.Interpreter = 'latex';
% cbar1.Ticks = []; % Disable the ticks
colormap(gca, Helper.Colormaps.plasma())
xlabel('$y$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14)
ylabel('$z$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14)
title('$|\Psi(x,y)|^2$', 'Interpreter', 'latex', 'FontSize', 14)
subplot(2,3,3)
nexttile;
plotxy = pcolor(x,y,nxy');
set(plotxy, 'EdgeColor', 'none');
xlabel('$x$ [$\mu$m]'); ylabel('$y$ [$\mu$m]');
cbar1 = colorbar;
cbar1.Label.Interpreter = 'latex';
% cbar1.Ticks = []; % Disable the ticks
colormap(gca, Helper.Colormaps.plasma())
xlabel('$x$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14)
ylabel('$y$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14)
title('$|\Psi(x,y)|^2$', 'Interpreter', 'latex', 'FontSize', 14)
subplot(2,3,4)
nexttile;
plot(-log10(Observ.residual),'-b')
ylabel('$-\mathrm{log}_{10}(r)$'); xlabel('steps');
subplot(2,3,5)
ylabel('$-\mathrm{log}_{10}(r)$', 'FontSize', 14); xlabel('Time steps', 'FontSize', 14);
title('Residual', 'FontSize', 14);
grid on
nexttile;
plot(Observ.EVec,'-b')
ylabel('$E$'); xlabel('steps');
ylabel('$E_{tot}$', 'FontSize', 14); xlabel('Time steps', 'FontSize', 14);
title('Total Energy', 'FontSize', 14);
grid on
subplot(2,3,6)
nexttile;
plot(Observ.mucVec,'-b')
ylabel('$\mu$'); xlabel('steps');
ylabel('$\mu$', 'FontSize', 14); xlabel('Time steps', 'FontSize', 14);
title('Chemical Potential', 'FontSize', 14);
grid on
end

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

@ -6,29 +6,30 @@
OptionsStruct = struct;
OptionsStruct.NumberOfAtoms = 1E5;
OptionsStruct.DipolarPolarAngle = 0;
OptionsStruct.NumberOfAtoms = 5E5;
OptionsStruct.DipolarPolarAngle = deg2rad(0);
OptionsStruct.DipolarAzimuthAngle = 0;
OptionsStruct.ScatteringLength = 86;
OptionsStruct.ScatteringLength = 85;
OptionsStruct.TrapFrequencies = [10, 10, 72.4];
OptionsStruct.TrapDepth = 5;
OptionsStruct.BoxSize = 15;
OptionsStruct.TrapFrequencies = [125, 125, 350];
OptionsStruct.TrapPotentialType = 'Harmonic';
OptionsStruct.NumberOfGridPoints = [256, 512, 256];
OptionsStruct.Dimensions = [50, 120, 150];
OptionsStruct.NumberOfGridPoints = [128, 128, 64];
OptionsStruct.Dimensions = [30, 30, 15];
OptionsStruct.CutoffType = 'Cylindrical';
OptionsStruct.SimulationMode = 'ImaginaryTimeEvolution'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution'
OptionsStruct.TimeStepSize = 0.005; % in s
OptionsStruct.NumberOfTimeSteps = 2E6; % in s
OptionsStruct.TimeStepSize = 0.002; % in s
OptionsStruct.MinimumTimeStepSize = 1E-6; % in s
OptionsStruct.TimeCutOff = 1E6; % in s
OptionsStruct.EnergyTolerance = 5E-10;
OptionsStruct.ResidualTolerance = 1E-05;
OptionsStruct.NoiseScaleFactor = 0.05;
OptionsStruct.PlotLive = true;
OptionsStruct.JobNumber = 1;
OptionsStruct.RunOnGPU = false;
OptionsStruct.SaveData = true;
OptionsStruct.SaveDirectory = './Data';
OptionsStruct.SaveDirectory = './Results/Data_3D';
options = Helper.convertstruct2cell(OptionsStruct);
clear OptionsStruct

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

@ -1,254 +1,82 @@
%% Tilting of the dipoles
% Atom Number = 1250 ppum
% System size = [sf * unitcell_x, sf * unitcell_x]
%% - Aspect Ratio: 2.8
ppum = 1250; % Atom Number Density in per micrometers
OptionsStruct = struct;
%% v_z = 500, theta = 0: a_s = 75.00
a = 1.8058;
scalingfactor = 5;
lx = scalingfactor*a;
ly = scalingfactor*sqrt(3)*a;
% Initialize OptionsStruct
OptionsStruct = struct;
% Assign values to OptionsStruct
OptionsStruct.NumberOfAtoms = ppum * (lx*ly);
OptionsStruct.NumberOfAtoms = 5E5;
OptionsStruct.DipolarPolarAngle = deg2rad(0);
OptionsStruct.DipolarAzimuthAngle = 0;
OptionsStruct.ScatteringLength = 75.00;
OptionsStruct.ScatteringLength = 85;
OptionsStruct.TrapFrequencies = [0, 0, 500];
OptionsStruct.TrapPotentialType = 'None';
AspectRatio = 2.8;
HorizontalTrapFrequency = 125;
VerticalTrapFrequency = AspectRatio * HorizontalTrapFrequency;
OptionsStruct.TrapFrequencies = [HorizontalTrapFrequency, HorizontalTrapFrequency, VerticalTrapFrequency];
OptionsStruct.TrapPotentialType = 'Harmonic';
OptionsStruct.NumberOfGridPoints = [128, 128];
OptionsStruct.Dimensions = [lx, ly];
OptionsStruct.TimeStepSize = 1E-3; % in s
OptionsStruct.MinimumTimeStepSize = 1E-5; % in s
OptionsStruct.TimeCutOff = 2E6; % in s
OptionsStruct.NumberOfGridPoints = [256, 256, 128];
OptionsStruct.Dimensions = [30, 30, 15];
OptionsStruct.CutoffType = 'Cylindrical';
OptionsStruct.SimulationMode = 'ImaginaryTimeEvolution'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution'
OptionsStruct.TimeStepSize = 0.002; % in s
OptionsStruct.MinimumTimeStepSize = 1E-6; % in s
OptionsStruct.TimeCutOff = 1E6; % in s
OptionsStruct.EnergyTolerance = 5E-10;
OptionsStruct.ResidualTolerance = 1E-05;
OptionsStruct.NoiseScaleFactor = 0.05;
OptionsStruct.IncludeDDICutOff = false;
OptionsStruct.MaxIterations = 10;
OptionsStruct.VariationalWidth = 1.15;
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/AdjustedSystemSize/Hz500';
OptionsStruct.SaveDirectory = sprintf('./Results/Data_3D/AspectRatio%s', strrep(num2str(AspectRatio), '.', '_'));
options = Helper.convertstruct2cell(OptionsStruct);
clear OptionsStruct
solver = VariationalSolver2D.DipolarGas(options{:});
pot = VariationalSolver2D.Potentials(options{:});
solver.Potential = pot.trap();
sim = Simulator.DipolarGas(options{:});
pot = Simulator.Potentials(options{:});
sim.Potential = pot.trap(); % + pot.repulsive_chopstick();
%-% Run Solver %-%
[Params, Transf, psi, V, VDk] = solver.run();
%-% Run Simulation %-%
[Params, Transf, psi, V, VDk] = sim.run();
%% v_z = 500, theta = 5: a_s = 75.00
%% - Aspect Ratio: 3.7
a = 1.795;
scalingfactor = 5;
lx = scalingfactor*a;
ly = scalingfactor*sqrt(3)*a;
OptionsStruct = struct;
% Initialize OptionsStruct
OptionsStruct = struct;
% Assign values to OptionsStruct
OptionsStruct.NumberOfAtoms = ppum * (lx*ly);
OptionsStruct.DipolarPolarAngle = deg2rad(5);
OptionsStruct.NumberOfAtoms = 5E5;
OptionsStruct.DipolarPolarAngle = deg2rad(0);
OptionsStruct.DipolarAzimuthAngle = 0;
OptionsStruct.ScatteringLength = 75.00;
OptionsStruct.ScatteringLength = 85;
OptionsStruct.TrapFrequencies = [0, 0, 500];
OptionsStruct.TrapPotentialType = 'None';
AspectRatio = 3.7;
HorizontalTrapFrequency = 125;
VerticalTrapFrequency = 125;
OptionsStruct.TrapFrequencies = [HorizontalTrapFrequency, HorizontalTrapFrequency, VerticalTrapFrequency];
OptionsStruct.TrapPotentialType = 'Harmonic';
OptionsStruct.NumberOfGridPoints = [128, 128];
OptionsStruct.Dimensions = [lx, ly];
OptionsStruct.TimeStepSize = 1E-3; % in s
OptionsStruct.MinimumTimeStepSize = 1E-5; % in s
OptionsStruct.TimeCutOff = 2E6; % in s
OptionsStruct.NumberOfGridPoints = [256, 256, 128];
OptionsStruct.Dimensions = [30, 30, 15];
OptionsStruct.CutoffType = 'Cylindrical';
OptionsStruct.SimulationMode = 'ImaginaryTimeEvolution'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution'
OptionsStruct.TimeStepSize = 0.002; % in s
OptionsStruct.MinimumTimeStepSize = 1E-6; % in s
OptionsStruct.TimeCutOff = 1E6; % in s
OptionsStruct.EnergyTolerance = 5E-10;
OptionsStruct.ResidualTolerance = 1E-05;
OptionsStruct.NoiseScaleFactor = 0.05;
OptionsStruct.IncludeDDICutOff = false;
OptionsStruct.MaxIterations = 10;
OptionsStruct.VariationalWidth = 1.15;
OptionsStruct.WidthLowerBound = 0.01;
OptionsStruct.WidthUpperBound = 12;
OptionsStruct.WidthCutoff = 1e-2;
OptionsStruct.PlotLive = false;
OptionsStruct.JobNumber = 1;
OptionsStruct.JobNumber = 0;
OptionsStruct.RunOnGPU = true;
OptionsStruct.SaveData = true;
OptionsStruct.SaveDirectory = './Results/Data_TiltingOfDipoles/AdjustedSystemSize/Hz500';
OptionsStruct.SaveDirectory = sprintf('./Results/Data_3D/AspectRatio%s', strrep(num2str(AspectRatio), '.', '_'));
options = Helper.convertstruct2cell(OptionsStruct);
clear OptionsStruct
solver = VariationalSolver2D.DipolarGas(options{:});
pot = VariationalSolver2D.Potentials(options{:});
solver.Potential = pot.trap();
sim = Simulator.DipolarGas(options{:});
pot = Simulator.Potentials(options{:});
sim.Potential = pot.trap(); % + pot.repulsive_chopstick();
%-% Run Solver %-%
[Params, Transf, psi, V, VDk] = solver.run();
%-% Run Simulation %-%
[Params, Transf, psi, V, VDk] = sim.run();
%% v_z = 500, theta = 7.5: a_s = 75.00
a = 2.055;
scalingfactor = 5;
lx = scalingfactor*a;
ly = scalingfactor*sqrt(3)*a;
% Initialize OptionsStruct
OptionsStruct = struct;
% Assign values to OptionsStruct
OptionsStruct.NumberOfAtoms = ppum * (lx*ly);
OptionsStruct.DipolarPolarAngle = deg2rad(7.5);
OptionsStruct.DipolarAzimuthAngle = 0;
OptionsStruct.ScatteringLength = 75.00;
OptionsStruct.TrapFrequencies = [0, 0, 500];
OptionsStruct.TrapPotentialType = 'None';
OptionsStruct.NumberOfGridPoints = [128, 128];
OptionsStruct.Dimensions = [lx, ly];
OptionsStruct.TimeStepSize = 1E-3; % 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.15;
OptionsStruct.WidthLowerBound = 0.01;
OptionsStruct.WidthUpperBound = 12;
OptionsStruct.WidthCutoff = 1e-2;
OptionsStruct.PlotLive = false;
OptionsStruct.JobNumber = 2;
OptionsStruct.RunOnGPU = true;
OptionsStruct.SaveData = true;
OptionsStruct.SaveDirectory = './Results/Data_TiltingOfDipoles/AdjustedSystemSize/Hz500';
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 = 10: a_s = 75.00
a = 2.055;
scalingfactor = 5;
lx = scalingfactor*a;
ly = scalingfactor*sqrt(3)*a;
% Initialize OptionsStruct
OptionsStruct = struct;
% Assign values to OptionsStruct
OptionsStruct.NumberOfAtoms = ppum * (lx*ly);
OptionsStruct.DipolarPolarAngle = deg2rad(10);
OptionsStruct.DipolarAzimuthAngle = 0;
OptionsStruct.ScatteringLength = 75.00;
OptionsStruct.TrapFrequencies = [0, 0, 500];
OptionsStruct.TrapPotentialType = 'None';
OptionsStruct.NumberOfGridPoints = [128, 128];
OptionsStruct.Dimensions = [lx, ly];
OptionsStruct.TimeStepSize = 1E-3; % 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.15;
OptionsStruct.WidthLowerBound = 0.01;
OptionsStruct.WidthUpperBound = 12;
OptionsStruct.WidthCutoff = 1e-2;
OptionsStruct.PlotLive = false;
OptionsStruct.JobNumber = 3;
OptionsStruct.RunOnGPU = true;
OptionsStruct.SaveData = true;
OptionsStruct.SaveDirectory = './Results/Data_TiltingOfDipoles/AdjustedSystemSize/Hz500';
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 = 15: a_s = 75.00
a = 2.0875;
scalingfactor = 5;
lx = scalingfactor*a;
ly = scalingfactor*sqrt(3)*a;
% Initialize OptionsStruct
OptionsStruct = struct;
% Assign values to OptionsStruct
OptionsStruct.NumberOfAtoms = ppum * (lx*ly);
OptionsStruct.DipolarPolarAngle = deg2rad(15);
OptionsStruct.DipolarAzimuthAngle = 0;
OptionsStruct.ScatteringLength = 75.00;
OptionsStruct.TrapFrequencies = [0, 0, 500];
OptionsStruct.TrapPotentialType = 'None';
OptionsStruct.NumberOfGridPoints = [128, 128];
OptionsStruct.Dimensions = [lx, ly];
OptionsStruct.TimeStepSize = 1E-3; % 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.15;
OptionsStruct.WidthLowerBound = 0.01;
OptionsStruct.WidthUpperBound = 12;
OptionsStruct.WidthCutoff = 1e-2;
OptionsStruct.PlotLive = false;
OptionsStruct.JobNumber = 4;
OptionsStruct.RunOnGPU = true;
OptionsStruct.SaveData = true;
OptionsStruct.SaveDirectory = './Results/Data_TiltingOfDipoles/AdjustedSystemSize/Hz500';
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();

36
Dipolar-Gas-Simulator/+Simulator/@Calculator/calculateChemicalPotential.m
View File

@ -1,29 +1,29 @@
function muchem = calculateChemicalPotential(~,psi,Params,Transf,VDk,V)
%Parameters
% Parameters
normfac = Params.Lx*Params.Ly*Params.Lz/numel(psi);
KEop= 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
KEop = 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
% DDIs
frho=fftn(abs(psi).^2);
Phi=real(ifftn(frho.*VDk));
frho = fftn(abs(psi).^2);
Phi = real(ifftn(frho.*VDk));
Eddi = (Params.gdd*Phi.*abs(psi).^2);
Eddi = (Params.gdd*Phi.*abs(psi).^2);
% Kinetic energy
Ekin = KEop.*abs(fftn(psi)*normfac).^2;
Ekin = sum(Ekin(:))*Transf.dkx*Transf.dky*Transf.dkz/(2*pi)^3;
%Kinetic energy
Ekin = KEop.*abs(fftn(psi)*normfac).^2;
Ekin = trapz(Ekin(:))*Transf.dkx*Transf.dky*Transf.dkz/(2*pi)^3;
% Potential energy
Epot = V.*abs(psi).^2;
%Potential energy
Epot = V.*abs(psi).^2;
% Contact interactions
Eint = Params.gs*abs(psi).^4;
%Contact interactions
Eint = Params.gs*abs(psi).^4;
% Quantum fluctuations
Eqf = Params.gammaQF*abs(psi).^5;
%Quantum fluctuations
Eqf = Params.gammaQF*abs(psi).^5;
%Total energy
muchem = Ekin + trapz(Epot(:) + Eint(:) + Eddi(:) + Eqf(:))*Transf.dx*Transf.dy*Transf.dz; %
muchem = muchem / Params.N;
% Total energy
muchem = Ekin + sum(Epot(:) + Eint(:) + Eddi(:) + Eqf(:))*Transf.dx*Transf.dy*Transf.dz; %
muchem = muchem / Params.N;
end

35
Dipolar-Gas-Simulator/+Simulator/@Calculator/calculateEnergyComponents.m
View File

@ -1,35 +1,30 @@
function E = calculateEnergyComponents(~,psi,Params,Transf,VDk,V)
%Parameters
KEop= 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
% Parameters
KEop = 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
normfac = Params.Lx*Params.Ly*Params.Lz/numel(psi);
% DDIs
frho = fftn(abs(psi).^2);
Phi = real(ifftn(frho.*VDk));
frho = fftn(abs(psi).^2);
Phi = real(ifftn(frho.*VDk));
Eddi = 0.5*Params.gdd*Phi.*abs(psi).^2;
E.Eddi = trapz(Eddi(:))*Transf.dx*Transf.dy*Transf.dz;
Eddi = 0.5*Params.gdd*Phi.*abs(psi).^2;
E.Eddi = sum(Eddi(:))*Transf.dx*Transf.dy*Transf.dz;
% EddiTot = trapz(Eddi(:))*Transf.dx*Transf.dy*Transf.dz;
%Kinetic energy
% psik = ifftshift(fftn(fftshift(psi)))*normfac;
Ekin = KEop.*abs(fftn(psi)*normfac).^2;
E.Ekin = trapz(Ekin(:))*Transf.dkx*Transf.dky*Transf.dkz/(2*pi)^3;
% Kinetic energy
Ekin = KEop.*abs(fftn(psi)*normfac).^2;
E.Ekin = sum(Ekin(:))*Transf.dkx*Transf.dky*Transf.dkz/(2*pi)^3;
% Potential energy
Epot = V.*abs(psi).^2;
E.Epot = trapz(Epot(:))*Transf.dx*Transf.dy*Transf.dz;
Epot = V.*abs(psi).^2;
E.Epot = sum(Epot(:))*Transf.dx*Transf.dy*Transf.dz;
%Contact interactions
Eint = 0.5*Params.gs*abs(psi).^4;
E.Eint = trapz(Eint(:))*Transf.dx*Transf.dy*Transf.dz;
Eint = 0.5*Params.gs*abs(psi).^4;
E.Eint = sum(Eint(:))*Transf.dx*Transf.dy*Transf.dz;
%Quantum fluctuations
Eqf = 0.4*Params.gammaQF*abs(psi).^5;
E.Eqf = trapz(Eqf(:))*Transf.dx*Transf.dy*Transf.dz;
Eqf = 0.4*Params.gammaQF*abs(psi).^5;
E.Eqf = sum(Eqf(:))*Transf.dx*Transf.dy*Transf.dz;
end

21
Dipolar-Gas-Simulator/+Simulator/@Calculator/calculateNormalizedResiduals.m
View File

@ -1,25 +1,26 @@
function res = calculateNormalizedResiduals(~,psi,Params,Transf,VDk,V,muchem)
KEop= 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
KEop = 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
% DDIs
frho=fftn(abs(psi).^2);
Phi=real(ifftn(frho.*VDk));
frho = fftn(abs(psi).^2);
Phi = real(ifftn(frho.*VDk));
Eddi = Params.gdd*Phi.*psi;
%Kinetic energy
% Kinetic energy
Ekin = ifftn(KEop.*fftn(psi));
%Potential energy
% Potential energy
Epot = V.*psi;
%Contact interactions
% Contact interactions
Eint = Params.gs*abs(psi).^2.*psi;
%Quantum fluctuations
Eqf = Params.gammaQF*abs(psi).^3.*psi;
% Quantum fluctuations
Eqf = Params.gammaQF*abs(psi).^3.*psi;
%Total energy
res = trapz(abs(Ekin(:) + Epot(:) + Eint(:) + Eddi(:) + Eqf(:) - muchem*psi(:))*Transf.dx*Transf.dy*Transf.dz)/trapz(abs(muchem*psi(:))*Transf.dx*Transf.dy*Transf.dz);
% Total energy
res = sum(abs(Ekin(:) + Epot(:) + Eint(:) + Eddi(:) + Eqf(:) - muchem*psi(:))*Transf.dx*Transf.dy*Transf.dz)/sum(abs(muchem*psi(:))*Transf.dx*Transf.dy*Transf.dz);
end

22
Dipolar-Gas-Simulator/+Simulator/@Calculator/calculateNumericalHankelTransform.m
View File

@ -5,25 +5,26 @@ function VDkSemi = calculateNumericalHankelTransform(~,kr,kz,Rmax,Zmax,Nr)
Nr = 5e4;
end
Nz = 64;
Nz = 64;
farRmultiple = 2000;
% midpoint grids for the integration over 0<z<Zmax, Rmax<r<farRmultiple*Rmax (i.e. starts at Rmax)
dr=(farRmultiple-1)*Rmax/Nr;
r = ((1:Nr)'-0.5)*dr+Rmax;
dz=Zmax/Nz;
z = ((1:Nz)-0.5)*dz;
[R, Z] = ndgrid(r,z);
Rsq = R.^2 + Z.^2;
dr = (farRmultiple-1)*Rmax/Nr;
r = ((1:Nr)'-0.5)*dr+Rmax;
dz = Zmax/Nz;
z = ((1:Nz)-0.5)*dz;
[R, Z] = ndgrid(r,z);
Rsq = R.^2 + Z.^2;
% real space interaction to be transformed
igrandbase = (1 - 3*Z.^2./Rsq)./Rsq.^(3/2);
igrandbase = (1 - 3*Z.^2./Rsq)./Rsq.^(3/2);
% do the Hankel/Fourier-Bessel transform numerically
% prestore to ensure each besselj is only calculated once
% cell is faster than (:,:,krn) slicing
Nkr = numel(kr);
besselr = cell(Nkr,1);
Nkr = numel(kr);
besselr = cell(Nkr,1);
for krn = 1:Nkr
besselr{krn} = repmat(r.*besselj(0,kr(krn)*r),1,Nz);
end
@ -36,4 +37,5 @@ function VDkSemi = calculateNumericalHankelTransform(~,kr,kz,Rmax,Zmax,Nr)
VDkSemi(krn,kzn) = VDkSemi(krn,kzn) - sum(igrand(:))*dz*dr;
end
end
end

72
Dipolar-Gas-Simulator/+Simulator/@Calculator/calculateOrderParameter.m
View File

@ -2,57 +2,57 @@ function [m_Order] = calculateOrderParameter(~,psi,Transf,Params,VDk,V,T,muchem)
NumRealiz = 100;
Mx = numel(Transf.x);
My = numel(Transf.y);
Mz = numel(Transf.z);
Mx = numel(Transf.x);
My = numel(Transf.y);
Mz = numel(Transf.z);
r = normrnd(0,1,size(psi));
theta = rand(size(psi));
noise = r.*exp(2*pi*1i*theta);
r = normrnd(0,1,size(psi));
theta = rand(size(psi));
noise = r.*exp(2*pi*1i*theta);
KEop= 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
Gamma = 1-1i*Params.gamma_S;
dt = Params.dt;
KEop = 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
Gamma = 1-1i*Params.gamma_S;
dt = Params.dt;
avgpsi = 0;
avgpsi2 = 0;
avgpsi1 = 0;
avgpsi2 = 0;
for jj = 1:NumRealiz
%generate initial state
xi = sqrt(2*Params.gamma_S*Params.kbol*T*10^(-9)*dt/(Params.hbar*Params.w0*Transf.dx*Transf.dy*Transf.dz));
swapx = randi(length(Transf.x),1,length(Transf.x));
swapy = randi(length(Transf.y),1,length(Transf.y));
swapz = randi(length(Transf.z),1,length(Transf.z));
psi_j = psi + xi * noise(swapx,swapy,swapz);
% Generate initial state
xi = sqrt(2*Params.gamma_S*Params.kbol*T*10^(-9)*dt/(Params.hbar*Params.w0*Transf.dx*Transf.dy*Transf.dz));
swapx = randi(length(Transf.x),1,length(Transf.x));
swapy = randi(length(Transf.y),1,length(Transf.y));
swapz = randi(length(Transf.z),1,length(Transf.z));
psi_j = psi + xi * noise(swapx,swapy,swapz);
% --- % propagate forward in time 1 time step:
%kin
psi_j = fftn(psi_j);
psi_j = psi_j.*exp(-0.5*1i*Gamma*dt*KEop);
psi_j = ifftn(psi_j);
% Kin
psi_j = fftn(psi_j);
psi_j = psi_j.*exp(-0.5*1i*Gamma*dt*KEop);
psi_j = ifftn(psi_j);
%DDI
frho = fftn(abs(psi_j).^2);
Phi = real(ifftn(frho.*VDk));
% DDI
frho = fftn(abs(psi_j).^2);
Phi = real(ifftn(frho.*VDk));
%Real-space
psi_j = psi_j.*exp(-1i*Gamma*dt*(V + Params.gs*abs(psi_j).^2 + Params.gammaQF*abs(psi_j).^3 + Params.gdd*Phi - muchem));
% Real-space
psi_j = psi_j.*exp(-1i*Gamma*dt*(V + Params.gs*abs(psi_j).^2 + Params.gammaQF*abs(psi_j).^3 + Params.gdd*Phi - muchem));
%kin
psi_j = fftn(psi_j);
psi_j = psi_j.*exp(-0.5*1i*Gamma*dt*KEop);
psi_j = ifftn(psi_j);
% Kin
psi_j = fftn(psi_j);
psi_j = psi_j.*exp(-0.5*1i*Gamma*dt*KEop);
psi_j = ifftn(psi_j);
%Projection
kcut = sqrt(2*Params.e_cut);
K = (Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2)<kcut.^2;
psi_j = ifftn(K.*fftn(psi_j));
% Projection
kcut = sqrt(2*Params.e_cut);
K = (Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2)<kcut.^2;
psi_j = ifftn(K.*fftn(psi_j));
% --- %
avgpsi = avgpsi + abs(sum(psi_j(:)))/NumRealiz;
avgpsi1 = avgpsi1 + abs(sum(psi_j(:))) /NumRealiz;
avgpsi2 = avgpsi2 + sum(abs(psi_j(:)).^2)/NumRealiz;
end
m_Order = 1/sqrt(Mx*My*Mz)*avgpsi/sqrt(avgpsi2);
m_Order = 1/sqrt(Mx*My*Mz)*avgpsi1/sqrt(avgpsi2);
end

20
Dipolar-Gas-Simulator/+Simulator/@Calculator/calculatePhaseCoherence.m
View File

@ -1,19 +1,21 @@
function [PhaseC] = calculatePhaseCoherence(~,psi,Transf,Params)
norm = sum(sum(sum(abs(psi).^2,1),2),3)*Transf.dx*Transf.dy*Transf.dz;
psi = psi/sqrt(norm);
norm = sum(sum(sum(abs(psi).^2,1),2),3)*Transf.dx*Transf.dy*Transf.dz;
psi = psi/sqrt(norm);
NumGlobalShifts = 800;
betas = []; phishift = [];
NumGlobalShifts = 800;
betas = [];
phishift = [];
for jj = 1:NumGlobalShifts
phishift(jj) = -pi + 2*pi*(jj-1)/NumGlobalShifts;
betas(jj) = sum(sum(sum(abs(angle(psi*exp(-1i*phishift(jj)))).*abs(psi).^2)));
betas(jj) = sum(sum(sum(abs(angle(psi*exp(-1i*phishift(jj)))).*abs(psi).^2)));
end
[minbeta,minidx] = min(betas);
psi = psi*exp(-1i*phishift(minidx));
phi = angle(psi);
psi = psi*exp(-1i*phishift(minidx));
phi = angle(psi);
avgphi = sum(sum(sum(phi.*abs(psi).^2,1),2),3)*Transf.dx*Transf.dy*Transf.dz;
PhaseC = sum(sum(sum(abs(angle(psi)-avgphi).*abs(psi).^2)))*Transf.dx*Transf.dy*Transf.dz;
avgphi = sum(sum(sum(phi.*abs(psi).^2,1),2),3)*Transf.dx*Transf.dy*Transf.dz;
PhaseC = sum(sum(sum(abs(angle(psi)-avgphi).*abs(psi).^2)))*Transf.dx*Transf.dy*Transf.dz;
end

38
Dipolar-Gas-Simulator/+Simulator/@Calculator/calculateTotalEnergy.m
View File

@ -1,32 +1,28 @@
function E = calculateTotalEnergy(~,psi,Params,Transf,VDk,V)
function E = calculateTotalEnergy(~,psi,Params,Transf,VDk,V)
%Parameters
KEop= 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
% Parameters
normfac = Params.Lx*Params.Ly*Params.Lz/numel(psi);
KEop = 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
% DDIs
frho = fftn(abs(psi).^2);
Phi = real(ifftn(frho.*VDk));
frho = fftn(abs(psi).^2);
Phi = real(ifftn(frho.*VDk));
Eddi = 0.5*Params.gdd*Phi.*abs(psi).^2;
Eddi = 0.5*Params.gdd*Phi.*abs(psi).^2;
% EddiTot = trapz(Eddi(:))*Transf.dx*Transf.dy*Transf.dz;
%Kinetic energy
% psik = ifftshift(fftn(fftshift(psi)))*normfac;
Ekin = KEop.*abs(fftn(psi)*normfac).^2;
Ekin = trapz(Ekin(:))*Transf.dkx*Transf.dky*Transf.dkz/(2*pi)^3;
% Kinetic energy
Ekin = KEop.*abs(fftn(psi)*normfac).^2;
Ekin = sum(Ekin(:))*Transf.dkx*Transf.dky*Transf.dkz/(2*pi)^3;
% Potential energy
Epot = V.*abs(psi).^2;
Epot = V.*abs(psi).^2;
%Contact interactions
Eint = 0.5*Params.gs*abs(psi).^4;
% Contact interactions
Eint = 0.5*Params.gs*abs(psi).^4;
%Quantum fluctuations
Eqf = 0.4*Params.gammaQF*abs(psi).^5;
% Quantum fluctuations
Eqf = 0.4*Params.gammaQF*abs(psi).^5;
E = Ekin + trapz(Epot(:) + Eint(:) + Eddi(:) + Eqf(:))*Transf.dx*Transf.dy*Transf.dz;
% Total energy
E = Ekin + sum(Epot(:) + Eint(:) + Eddi(:) + Eqf(:))*Transf.dx*Transf.dy*Transf.dz;
end

73
Dipolar-Gas-Simulator/+Simulator/@Calculator/calculateVDk.m Normal file
View File

@ -0,0 +1,73 @@
function VDk = calculateVDk(this,Params,Transf,TransfRad,IncludeDDICutOff)
% makes the dipolar interaction matrix, size numel(Params.kr) * numel(Params.kz)
% Rmax and Zmax are the interaction cutoffs
% VDk needs to be multiplied by Cdd
% approach is that of Lu, PRA 82, 023622 (2010)
% == Calculating the DDI potential in Fourier space with appropriate cutoff == %
% Cylindrical (semianalytic)
% Cylindrical infinite Z, polarized along x (analytic)
% Spherical
if IncludeDDICutOff
switch this.CutoffType
case 'Cylindrical' % Cylindrical (semianalytic)
Zcutoff = Params.Lz/2;
alph = acos((Transf.KX*sin(Params.theta)*cos(Params.phi)+Transf.KY*sin(Params.theta)*sin(Params.phi)+Transf.KZ*cos(Params.theta))./sqrt(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2));
alph(1) = pi/2;
% Analytic part of cutoff for slice 0<z<Zmax, 0<r<Inf Ronen, PRL 98, 030406 (2007)
cossq = cos(alph).^2;
VDk = cossq-1/3;
sinsq = 1 - cossq;
VDk = VDk + exp(-Zcutoff*sqrt(Transf.KX.^2+Transf.KY.^2)).*( sinsq .* cos(Zcutoff * Transf.KZ) - sqrt(sinsq.*cossq).*sin(Zcutoff * Transf.KZ) );
% Nonanalytic part
% For a cylindrical cutoff, we need to construct a kr grid based on the 3D parameters using Bessel quadrature
VDkNon = this.calculateNumericalHankelTransform(TransfRad.kr, TransfRad.kz, TransfRad.Rmax, Zcutoff);
% Interpolating the nonanalytic part onto 3D grid
fullkr = [-flip(TransfRad.kr)',TransfRad.kr'];
[KR,KZ] = ndgrid(fullkr,TransfRad.kz);
[KX3D,KY3D,KZ3D] = ndgrid(ifftshift(Transf.kx),ifftshift(Transf.ky),ifftshift(Transf.kz));
KR3D = sqrt(KX3D.^2 + KY3D.^2);
fullVDK = [flip(VDkNon',2),VDkNon']';
VDkNon = interpn(KR,KZ,fullVDK,KR3D,KZ3D,'spline',0); %Last argument is -1/3 for full VDk. 0 for nonanalytic piece?
VDkNon = fftshift(VDkNon);
VDk = VDk + VDkNon;
case 'CylindricalInfiniteZ' % Cylindrical infinite Z, polarized along x -- PRA 107, 033301 (2023)
alph = acos((Transf.KX*sin(Params.theta)*cos(Params.phi)+Transf.KY*sin(Params.theta)*sin(Params.phi)+Transf.KZ*cos(Params.theta))./sqrt(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2));
alph(1) = pi/2;
rhoc = max([abs(Transf.x),abs(Transf.y)]);
KR = sqrt(Transf.KX.^2+Transf.KY.^2);
func = @(n,u,v) v.^2./(u.^2+v.^2).*(v.*besselj(n,u).*besselk(n+1,v) - u.*besselj(n+1,u).*besselk(n,v));
VDk = -0.5*func(0,KR*rhoc,abs(Transf.KZ)*rhoc) + (Transf.KX.^2./KR.^2 - 0.5).*func(2,KR*rhoc,abs(Transf.KZ)*rhoc);
VDk = (1/3)*(3*(cos(alph).^2)-1) - VDk;
VDk(KR==0) = -1/3 + 1/2*abs(Transf.KZ(KR==0))*rhoc.*besselk(1,abs(Transf.KZ(KR==0))*rhoc);
VDk(Transf.KZ==0) = 1/6 + (Transf.KX(Transf.KZ==0).^2-Transf.KY(Transf.KZ==0).^2)./(KR(Transf.KZ==0).^2).*(1/2 - besselj(1,KR(Transf.KZ==0)*rhoc)./(KR(Transf.KZ==0)*rhoc));
VDk(1,1,1) = 1/6;
case 'Spherical' % Spherical
Rcut = min(Params.Lx/2,Params.Ly/2,Params.Lz/2);
alph = acos((Transf.KX*sin(Params.theta)*cos(Params.phi)+Transf.KY*sin(Params.theta)*sin(Params.phi)+Transf.KZ*cos(Params.theta))./sqrt(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2));
alph(1) = pi/2;
K = sqrt(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
VDk = (cos(alph).^2-1/3).*(1 + 3*cos(Rcut*K)./(Rcut^2.*K.^2) - 3*sin(Rcut*K)./(Rcut^3.*K.^3));
otherwise
disp('Choose a valid DDI cutoff type!')
return
end
else
alph = acos((Transf.KX*sin(Params.theta)*cos(Params.phi)+Transf.KY*sin(Params.theta)*sin(Params.phi)+Transf.KZ*cos(Params.theta))./sqrt(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2));
alph(1) = pi/2;
VDk = cos(alph).^2-1/3;
end
VDk = 3 * VDk;
end

64
Dipolar-Gas-Simulator/+Simulator/@Calculator/calculateVDkWithCutoff.m
View File

@ -1,64 +0,0 @@
function VDk = calculateVDCutoff(this,Params,Transf,TransfRad)
% makes the dipolar interaction matrix, size numel(Params.kr) * numel(Params.kz)
% Rmax and Zmax are the interaction cutoffs
% VDk needs to be multiplied by Cdd
% approach is that of Lu, PRA 82, 023622 (2010)
% == Calculating the DDI potential in Fourier space with appropriate cutoff == %
% Cylindrical (semianalytic)
% Cylindrical infinite Z, polarized along x (analytic)
% Spherical
switch this.CutoffType
case 'Cylindrical' %Cylindrical (semianalytic)
Zcutoff = Params.Lz/2;
alph = acos((Transf.KX*sin(Params.theta)*cos(Params.phi)+Transf.KY*sin(Params.theta)*sin(Params.phi)+Transf.KZ*cos(Params.theta))./sqrt(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2));
alph(1) = pi/2;
% Analytic part of cutoff for slice 0<z<Zmax, 0<r<Inf Ronen, PRL 98, 030406 (2007)
cossq = cos(alph).^2;
VDk = cossq-1/3;
sinsq = 1 - cossq;
VDk = VDk + exp(-Zcutoff*sqrt(Transf.KX.^2+Transf.KY.^2)).*( sinsq .* cos(Zcutoff * Transf.KZ) - sqrt(sinsq.*cossq).*sin(Zcutoff * Transf.KZ) );
% Nonanalytic part
% For a cylindrical cutoff, we need to construct a kr grid based on the 3D parameters using Bessel quadrature
VDkNon = this.calculateNumericalHankelTransform(TransfRad.kr, TransfRad.kz, TransfRad.Rmax, Zcutoff);
% Interpolating the nonanalytic part onto 3D grid
fullkr = [-flip(TransfRad.kr)',TransfRad.kr'];
[KR,KZ] = ndgrid(fullkr,TransfRad.kz);
[KX3D,KY3D,KZ3D] = ndgrid(ifftshift(Transf.kx),ifftshift(Transf.ky),ifftshift(Transf.kz));
KR3D = sqrt(KX3D.^2 + KY3D.^2);
fullVDK = [flip(VDkNon',2),VDkNon']';
VDkNon = interpn(KR,KZ,fullVDK,KR3D,KZ3D,'spline',0); %Last argument is -1/3 for full VDk. 0 for nonanalytic piece?
VDkNon = fftshift(VDkNon);
VDk = VDk + VDkNon;
case 'CylindricalInfiniteZ' %Cylindrical infinite Z, polarized along x -- PRA 107, 033301 (2023)
alph = acos((Transf.KX*sin(Params.theta)*cos(Params.phi)+Transf.KY*sin(Params.theta)*sin(Params.phi)+Transf.KZ*cos(Params.theta))./sqrt(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2));
alph(1) = pi/2;
rhoc = max([abs(Transf.x),abs(Transf.y)]);
KR = sqrt(Transf.KX.^2+Transf.KY.^2);
func = @(n,u,v) v.^2./(u.^2+v.^2).*(v.*besselj(n,u).*besselk(n+1,v) - u.*besselj(n+1,u).*besselk(n,v));
VDk = -0.5*func(0,KR*rhoc,abs(Transf.KZ)*rhoc) + (Transf.KX.^2./KR.^2 - 0.5).*func(2,KR*rhoc,abs(Transf.KZ)*rhoc);
VDk = (1/3)*(3*(cos(alph).^2)-1) - VDk;
VDk(KR==0) = -1/3 + 1/2*abs(Transf.KZ(KR==0))*rhoc.*besselk(1,abs(Transf.KZ(KR==0))*rhoc);
VDk(Transf.KZ==0) = 1/6 + (Transf.KX(Transf.KZ==0).^2-Transf.KY(Transf.KZ==0).^2)./(KR(Transf.KZ==0).^2).*(1/2 - besselj(1,KR(Transf.KZ==0)*rhoc)./(KR(Transf.KZ==0)*rhoc));
VDk(1,1,1) = 1/6;
case 'Spherical' %Spherical
Rcut = min(Params.Lx/2,Params.Ly/2,Params.Lz/2);
alph = acos((Transf.KX*sin(Params.theta)*cos(Params.phi)+Transf.KY*sin(Params.theta)*sin(Params.phi)+Transf.KZ*cos(Params.theta))./sqrt(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2));
alph(1) = pi/2;
K = sqrt(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
VDk = (cos(alph).^2-1/3).*(1 + 3*cos(Rcut*K)./(Rcut^2.*K.^2) - 3*sin(Rcut*K)./(Rcut^3.*K.^3));
otherwise
disp('Choose a valid DDI cutoff type!')
return
end
end

36
Dipolar-Gas-Simulator/+Simulator/@DipolarGas/DipolarGas.m
View File

@ -12,15 +12,17 @@ classdef DipolarGas < handle & matlab.mixin.Copyable
SimulationMode;
TimeStepSize;
NumberOfTimeSteps;
MinimumTimeStepSize;
TimeCutOff;
EnergyTolerance;
ResidualTolerance;
MinimumTimeStepSize;
NoiseScaleFactor;
Calculator;
SimulationParameters;
IncludeDDICutOff;
PlotLive;
JobNumber;
RunOnGPU;
DebugMode;
@ -56,14 +58,21 @@ classdef DipolarGas < handle & matlab.mixin.Copyable
@(x) assert(any(strcmpi(x,{'Cylindrical','CylindricalInfiniteZ', 'Spherical'}))));
addParameter(p, 'TimeStepSize', 5E-4,...
@(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));
addParameter(p, 'NumberOfTimeSteps', 2e6,...
addParameter(p, 'MinimumTimeStepSize', 1e-6,...
@(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));
addParameter(p, 'TimeCutOff', 2e6,...
@(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));
addParameter(p, 'EnergyTolerance', 1e-10,...
@(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));
addParameter(p, 'ResidualTolerance', 1e-10,...
@(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));
addParameter(p, 'MinimumTimeStepSize', 1e-6,...
@(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));
addParameter(p, 'NoiseScaleFactor', 4,...
@(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));
addParameter(p, 'IncludeDDICutOff', true,...
@islogical);
addParameter(p, 'PlotLive', false,...
@islogical);
addParameter(p, 'JobNumber', 1,...
@(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));
addParameter(p, 'RunOnGPU', false,...
@ -87,11 +96,14 @@ classdef DipolarGas < handle & matlab.mixin.Copyable
this.Potential = NaN;
this.SimulationMode = p.Results.SimulationMode;
this.TimeStepSize = p.Results.TimeStepSize;
this.NumberOfTimeSteps = p.Results.NumberOfTimeSteps;
this.MinimumTimeStepSize = p.Results.MinimumTimeStepSize;
this.TimeCutOff = p.Results.TimeCutOff;
this.EnergyTolerance = p.Results.EnergyTolerance;
this.ResidualTolerance = p.Results.ResidualTolerance;
this.MinimumTimeStepSize = p.Results.MinimumTimeStepSize;
this.NoiseScaleFactor = p.Results.NoiseScaleFactor;
this.IncludeDDICutOff = p.Results.IncludeDDICutOff;
this.PlotLive = p.Results.PlotLive;
this.JobNumber = p.Results.JobNumber;
this.RunOnGPU = p.Results.RunOnGPU;
this.DebugMode = p.Results.DebugMode;
@ -113,12 +125,12 @@ classdef DipolarGas < handle & matlab.mixin.Copyable
function ret = get.TimeStepSize(this)
ret = this.TimeStepSize;
end
function set.NumberOfTimeSteps(this, val)
function set.TimeCutOff(this, val)
assert(val <= 2E6, 'Not time efficient to compute for time spans longer than 2E6 seconds!');
this.NumberOfTimeSteps = val;
this.TimeCutOff = val;
end
function ret = get.NumberOfTimeSteps(this)
ret = this.NumberOfTimeSteps;
function ret = get.TimeCutOff(this)
ret = this.TimeCutOff;
end
function set.NumberOfAtoms(this, val)
assert(val <= 1E6, '!!Not time efficient to compute for atom numbers larger than 1,000,000!!');

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

@ -8,8 +8,12 @@ function [psi,V,VDk] = initialize(this,Params,Transf,TransfRad)
if isfile(strcat(this.SaveDirectory, '/VDk_M.mat'))
VDk = load(sprintf(strcat(this.SaveDirectory, '/VDk_M.mat')));
VDk = VDk.VDk;
if ~isequal(size(VDk), this.NumberOfGridPoints)
VDk = this.Calculator.calculateVDk(Params,Transf,TransfRad, this.IncludeDDICutOff);
save(sprintf(strcat(this.SaveDirectory, '/VDk_M.mat')),'VDk');
end
else
VDk = this.Calculator.calculateVDkWithCutoff(Params,Transf,TransfRad);
VDk = this.Calculator.calculateVDk(Params,Transf,TransfRad, this.IncludeDDICutOff);
save(sprintf(strcat(this.SaveDirectory, '/VDk_M.mat')),'VDk');
end
fprintf('Computed and saved DDI potential in Fourier space with %s cutoff.\n', this.Calculator.CutoffType)

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

@ -4,63 +4,70 @@ function [psi] = propagateWavefunction(this,psi,Params,Transf,VDk,V,t_idx,Observ
switch this.SimulationMode
case 'ImaginaryTimeEvolution'
dt=-1j*abs(this.TimeStepSize);
dt = -1j*abs(this.TimeStepSize); % Imaginary Time
KEop= 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
Observ.residual = 1; Observ.res = 1;
KEop = 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
Observ.residual = 1;
Observ.res = 1;
muchem = this.Calculator.calculateChemicalPotential(psi,Params,Transf,VDk,V);
AdaptIdx = 0;
muchem = this.Calculator.calculateChemicalPotential(psi,Params,Transf,VDk,V);
AdaptIdx = 0;
pb = Helper.ProgressBar();
pb.run('Running simulation in imaginary time: ');
if this.PlotLive
Plotter.plotLive(psi,Params,Transf,Observ)
drawnow
end
while t_idx < Params.sim_time_cut_off
%kin
psi = fftn(psi);
psi = psi.*exp(-0.5*1i*dt*KEop);
psi = ifftn(psi);
% kin
psi = fftn(psi);
psi = psi.*exp(-0.5*1i*dt*KEop);
psi = ifftn(psi);
%DDI
frho = fftn(abs(psi).^2);
Phi = real(ifftn(frho.*VDk));
% DDI
frho = fftn(abs(psi).^2);
Phi = real(ifftn(frho.*VDk));
%Real-space
psi = psi.*exp(-1i*dt*(V + Params.gs*abs(psi).^2 + Params.gdd*Phi + Params.gammaQF*abs(psi).^3 - muchem));
% Real-space
psi = psi.*exp(-1i*dt*(V + Params.gs*abs(psi).^2 + Params.gdd*Phi + Params.gammaQF*abs(psi).^3 - muchem));
%kin
psi = fftn(psi);
psi = psi.*exp(-0.5*1i*dt*KEop);
psi = ifftn(psi);
% kin
psi = fftn(psi);
psi = psi.*exp(-0.5*1i*dt*KEop);
psi = ifftn(psi);
%Renorm
Norm = trapz(abs(psi(:)).^2)*Transf.dx*Transf.dy*Transf.dz;
psi = sqrt(Params.N)*psi/sqrt(Norm);
% Renorm
Norm = sum(abs(psi(:)).^2)*Transf.dx*Transf.dy*Transf.dz;
psi = sqrt(Params.N)*psi/sqrt(Norm);
muchem = this.Calculator.calculateChemicalPotential(psi,Params,Transf,VDk,V);
muchem = this.Calculator.calculateChemicalPotential(psi,Params,Transf,VDk,V);
if mod(t_idx,1000) == 0
if mod(t_idx,500) == 0
%Change in Energy
E = this.Calculator.calculateTotalEnergy(psi,Params,Transf,VDk,V);
E = E/Norm;
Observ.EVec = [Observ.EVec E];
% Change in Energy
E = this.Calculator.calculateTotalEnergy(psi,Params,Transf,VDk,V);
E = E/Norm;
Observ.EVec = [Observ.EVec E];
%Chemical potential
Observ.mucVec = [Observ.mucVec muchem];
% Chemical potential
Observ.mucVec = [Observ.mucVec muchem];
%Normalized residuals
res = this.Calculator.calculateNormalizedResiduals(psi,Params,Transf,VDk,V,muchem);
% Normalized residuals
res = this.Calculator.calculateNormalizedResiduals(psi,Params,Transf,VDk,V,muchem);
Observ.residual = [Observ.residual res];
Observ.res_idx = Observ.res_idx + 1;
Observ.res_idx = Observ.res_idx + 1;
save(sprintf('./Data/Run_%03i/psi_gs.mat',Params.njob),'psi','muchem','Observ','t_idx','Transf','Params','VDk','V');
if this.PlotLive
Plotter.plotLive(psi,Params,Transf,Observ)
drawnow
end
save(sprintf(strcat(this.SaveDirectory, '/Run_%03i/psi_gs.mat'),Params.njob),'psi','muchem','Observ','t_idx','Transf','Params','VDk','V');
%Adaptive time step -- Careful, this can quickly get out of control
relres = abs(Observ.residual(Observ.res_idx)-Observ.residual(Observ.res_idx-1))/Observ.residual(Observ.res_idx);
if relres <1e-5
if relres < 1e-5
if AdaptIdx > 4 && abs(dt) > Params.mindt
dt = dt / 2;
fprintf('Time step changed to '); disp(dt);
@ -79,23 +86,22 @@ function [psi] = propagateWavefunction(this,psi,Params,Transf,VDk,V,t_idx,Observ
break
end
t_idx=t_idx+1;
pb.run(100*t_idx/Params.sim_time_cut_off);
end
%Change in Energy
E = this.Calculator.calculateTotalEnergy(psi,Params,Transf,VDk,V);
E = E/Norm;
Observ.EVec = [Observ.EVec E];
% Change in Energy
E = this.Calculator.calculateTotalEnergy(psi,Params,Transf,VDk,V);
E = E/Norm;
Observ.EVec = [Observ.EVec E];
% Phase coherence
[PhaseC] = this.Calculator.calculatePhaseCoherence(psi,Transf,Params);
Observ.PCVec = [Observ.PCVec PhaseC];
[PhaseC] = this.Calculator.calculatePhaseCoherence(psi,Transf,Params);
Observ.PCVec = [Observ.PCVec PhaseC];
Observ.res_idx = Observ.res_idx + 1;
pb.run(' - Job Completed!');
disp('Run completed!');
disp('Saving data...');
save(sprintf('./Data/Run_%03i/psi_gs.mat',Params.njob),'psi','muchem','Observ','t_idx','Transf','Params','VDk','V');
save(sprintf(strcat(this.SaveDirectory, '/Run_%03i/psi_gs.mat'),Params.njob),'psi','muchem','Observ','t_idx','Transf','Params','VDk','V');
disp('Save complete!');
case 'RealTimeEvolution'
@ -106,9 +112,6 @@ function [psi] = propagateWavefunction(this,psi,Params,Transf,VDk,V,t_idx,Observ
muchem = this.Calculator.calculateChemicalPotential(psi,Params,Transf,VDk,V);
pb = Helper.ProgressBar();
pb.run('Running simulation in real time: ');
while t_idx < Params.sim_time_cut_off
% Parameters at time t
@ -134,7 +137,7 @@ function [psi] = propagateWavefunction(this,psi,Params,Transf,VDk,V,t_idx,Observ
if mod(t_idx,1000)==0
%Change in Normalization
Norm = trapz(abs(psi(:)).^2)*Transf.dx*Transf.dy*Transf.dz; % normalisation
Norm = sum(abs(psi(:)).^2)*Transf.dx*Transf.dy*Transf.dz; % normalisation
Observ.NormVec = [Observ.NormVec Norm];
%Change in Energy
@ -156,11 +159,10 @@ function [psi] = propagateWavefunction(this,psi,Params,Transf,VDk,V,t_idx,Observ
break
end
t_idx=t_idx+1;
pb.run(100*t_idx/Params.sim_time_cut_off);
end
%Change in Normalization
Norm = trapz(abs(psi(:)).^2)*Transf.dx*Transf.dy*Transf.dz; % normalisation
Norm = sum(abs(psi(:)).^2)*Transf.dx*Transf.dy*Transf.dz; % normalisation
Observ.NormVec = [Observ.NormVec Norm];
%Change in Energy
@ -175,7 +177,7 @@ function [psi] = propagateWavefunction(this,psi,Params,Transf,VDk,V,t_idx,Observ
Observ.tVecPlot = [Observ.tVecPlot tVal];
Observ.res_idx = Observ.res_idx + 1;
pb.run(' - Run Completed!');
disp('Run completed!');
disp('Saving data...');
save(sprintf('./Data/Run_%03i/TimeEvolution/psi_%i.mat',Params.njob,Observ.res_idx),'psi','muchem','Observ','t_idx');
disp('Save complete!');

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

@ -1,22 +1,23 @@
function [Params, Transf, psi,V,VDk] = run(this)
% --- Obtain simulation parameters ---
[Params] = this.setupParameters();
[Params] = this.setupParameters();
this.SimulationParameters = Params;
% --- Set up spatial grids and transforms ---
[Transf] = this.setupSpace(Params);
[TransfRad] = this.setupSpaceRadial(Params);
[Transf] = this.setupSpace(Params);
[TransfRad] = this.setupSpaceRadial(Params);
% --- Initialize ---
mkdir(sprintf(this.SaveDirectory))
[psi,V,VDk] = this.initialize(Params,Transf,TransfRad);
[psi,V,VDk] = this.initialize(Params,Transf,TransfRad);
Observ.EVec = []; Observ.NormVec = []; Observ.PCVec = []; Observ.tVecPlot = []; Observ.mucVec = [];
t_idx = 1; % Start at t = 0;
Observ.res_idx = 1;
t_idx = 1; % Start at t = 0;
Observ.res_idx = 1;
% --- Run Simulation ---
mkdir(sprintf('./Data/Run_%03i',Params.njob))
mkdir(sprintf(this.SaveDirectory))
mkdir(sprintf(strcat(this.SaveDirectory, '/Run_%03i'),Params.njob))
[psi] = this.propagateWavefunction(psi,Params,Transf,VDk,V,t_idx,Observ);
end

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

@ -6,7 +6,7 @@ function [Params] = setupParameters(this)
muB = CONSTANTS.BohrMagneton; % [J/T]
a0 = CONSTANTS.BohrRadius; % [m]
m0 = CONSTANTS.AtomicMassUnit; % [kg]
w0 = 2*pi*100; % 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
@ -21,7 +21,6 @@ function [Params] = setupParameters(this)
% Mass, length scale
Params.m = CONSTANTS.Dy164Mass;
l0 = sqrt(hbar/(Params.m*w0)); % Defining a harmonic oscillator length
% Atom numbers
Params.N = this.NumberOfAtoms;
@ -44,14 +43,18 @@ function [Params] = setupParameters(this)
% Tolerances
Params.Etol = this.EnergyTolerance;
Params.rtol = this.ResidualTolerance;
Params.sim_time_cut_off = this.NumberOfTimeSteps; % sometimes the imaginary time gets a little stuck
% even though the solution is good, this just stops it going on forever
Params.mindt = this.MinimumTimeStepSize; % Minimum size for a time step using adaptive dt
Params.sim_time_cut_off = this.TimeCutOff; % sometimes the imaginary time gets a little stuck
% even though the solution is good, this just stops it going on forever
Params.mindt = this.MinimumTimeStepSize; % Minimum size for a time step using adaptive dt
Params.njob = this.JobNumber;
Params.nsf = this.NoiseScaleFactor;
Params.njob = this.JobNumber;
% ================ Parameters defined by those above ================ %
% Length scale
l0 = sqrt(hbar/(Params.m*w0)); % Defining a harmonic oscillator length
% Contact interaction strength (units of l0/m)
Params.gs = 4*pi*Params.as/l0;
@ -59,7 +62,7 @@ function [Params] = setupParameters(this)
Params.add = mu0*Params.mu^2*Params.m/(12*pi*hbar^2);
% DDI strength
Params.gdd = 12*pi*Params.add/l0; %sometimes the 12 is a 4 --> depends on how Vdk (DDI) is defined
Params.gdd = 4*pi*Params.add/l0; %sometimes the 12 is a 4 --> depends on how Vdk (DDI) is defined
% Trap gamma
Params.gx = (Params.wx/w0)^2;

65
Dipolar-Gas-Simulator/+Simulator/@DipolarGas/setupSpace.m
View File

@ -1,33 +1,48 @@
function [Transf] = setupSpace(~,Params)
Transf.Xmax = 0.5*Params.Lx;
Transf.Ymax = 0.5*Params.Ly;
Transf.Zmax = 0.5*Params.Lz;
Nz = Params.Nz; Nx = Params.Nx; Ny = Params.Ny;
Transf.Xmax = 0.5*Params.Lx;
Transf.Ymax = 0.5*Params.Ly;
Transf.Zmax = 0.5*Params.Lz;
Nz = Params.Nz;
Nx = Params.Nx;
Ny = Params.Ny;
% Fourier grids
x = linspace(-0.5*Params.Lx,0.5*Params.Lx-Params.Lx/Params.Nx,Params.Nx);
Kmax = pi*Params.Nx/Params.Lx;
kx = linspace(-Kmax,Kmax,Nx+1);
kx = kx(1:end-1); dkx = kx(2)-kx(1);
kx = fftshift(kx);
x = linspace(-0.5*Params.Lx,0.5*Params.Lx-Params.Lx/Params.Nx,Params.Nx);
Kmax = pi*Params.Nx/Params.Lx;
kx = linspace(-Kmax,Kmax,Nx+1);
kx = kx(1:end-1);
dkx = kx(2)-kx(1);
kx = fftshift(kx);
y = linspace(-0.5*Params.Ly,0.5*Params.Ly-Params.Ly/Params.Ny,Params.Ny);
Kmax = pi*Params.Ny/Params.Ly;
ky = linspace(-Kmax,Kmax,Ny+1);
ky = ky(1:end-1); dky = ky(2)-ky(1);
ky = fftshift(ky);
y = linspace(-0.5*Params.Ly,0.5*Params.Ly-Params.Ly/Params.Ny,Params.Ny);
Kmax = pi*Params.Ny/Params.Ly;
ky = linspace(-Kmax,Kmax,Ny+1);
ky = ky(1:end-1);
dky = ky(2)-ky(1);
ky = fftshift(ky);
z = linspace(-0.5*Params.Lz,0.5*Params.Lz-Params.Lz/Params.Nz,Params.Nz);
Kmax = pi*Params.Nz/Params.Lz;
kz = linspace(-Kmax,Kmax,Nz+1);
kz = kz(1:end-1); dkz = kz(2)-kz(1);
kz = fftshift(kz);
z = linspace(-0.5*Params.Lz,0.5*Params.Lz-Params.Lz/Params.Nz,Params.Nz);
Kmax = pi*Params.Nz/Params.Lz;
kz = linspace(-Kmax,Kmax,Nz+1);
kz = kz(1:end-1);
dkz = kz(2)-kz(1);
kz = fftshift(kz);
[Transf.X,Transf.Y,Transf.Z] = ndgrid(x,y,z);
[Transf.KX,Transf.KY,Transf.KZ] = ndgrid(kx,ky,kz);
Transf.x = x;
Transf.y = y;
Transf.z = z;
Transf.kx = kx;
Transf.ky = ky;
Transf.kz = kz;
Transf.dx = x(2)-x(1);
Transf.dy = y(2)-y(1);
Transf.dz = z(2)-z(1);
Transf.dkx = dkx;
Transf.dky = dky;
Transf.dkz = dkz;
[Transf.X,Transf.Y,Transf.Z]=ndgrid(x,y,z);
[Transf.KX,Transf.KY,Transf.KZ]=ndgrid(kx,ky,kz);
Transf.x = x; Transf.y = y; Transf.z = z;
Transf.kx = kx; Transf.ky = ky; Transf.kz = kz;
Transf.dx = x(2)-x(1); Transf.dy = y(2)-y(1); Transf.dz = z(2)-z(1);
Transf.dkx = dkx; Transf.dky = dky; Transf.dkz = dkz;
end

23
Dipolar-Gas-Simulator/+Simulator/@DipolarGas/setupWavefunction.m
View File

@ -1,25 +1,32 @@
function [psi] = setupWavefunction(~,Params,Transf)
X = Transf.X; Y = Transf.Y; Z = Transf.Z;
X = Transf.X;
Y = Transf.Y;
Z = Transf.Z;
ellx = sqrt(Params.hbar/(Params.m*Params.wx))/Params.l0;
elly = sqrt(Params.hbar/(Params.m*Params.wy))/Params.l0;
ellz = sqrt(Params.hbar/(Params.m*Params.wz))/Params.l0;
Rx = 4*ellx; Ry = 4*elly; Rz = 4*ellz;
X0 = 0.0*Transf.Xmax; Y0 = 0.0*Transf.Ymax; Z0 = 0*Transf.Zmax;
Rx = 8*ellx;
Ry = 8*elly;
Rz = 8*ellz;
X0 = 0.0*Transf.Xmax;
Y0 = 0.0*Transf.Ymax;
Z0 = 0*Transf.Zmax;
psi = exp(-(X-X0).^2/Rx^2-(Y-Y0).^2/Ry^2-(Z-Z0).^2/Rz^2);
cur_norm = trapz(abs(psi(:)).^2)*Transf.dx*Transf.dy*Transf.dz;
cur_norm = sum(abs(psi(:)).^2)*Transf.dx*Transf.dy*Transf.dz;
psi = psi/sqrt(cur_norm);
% add some noise
% --- Adding some noise ---
r = normrnd(0,1,size(X));
theta = rand(size(X));
noise = r.*exp(2*pi*1i*theta);
psi = psi + 0.01*noise;
psi = psi + Params.nsf*noise;
Norm = trapz(abs(psi(:)).^2)*Transf.dx*Transf.dy*Transf.dz;
psi = sqrt(Params.N)*psi/sqrt(Norm);
% Renormalize wavefunction
Norm = sum(abs(psi(:)).^2)*Transf.dx*Transf.dy*Transf.dz;
psi = sqrt(Params.N)*psi/sqrt(Norm);
end