diff --git a/Dipolar-Gas-Simulator/+Plotter/plotLive.m b/Dipolar-Gas-Simulator/+Plotter/plotLive.m index bec2d08..3f8e5ba 100644 --- a/Dipolar-Gas-Simulator/+Plotter/plotLive.m +++ b/Dipolar-Gas-Simulator/+Plotter/plotLive.m @@ -1,4 +1,4 @@ -function plotLive(psi,Params,Transf,Observ) +function plotLive(psi,Params,Transf,Observ,SimulationMode) set(0,'defaulttextInterpreter','latex') set(groot, 'defaultAxesTickLabelInterpreter','latex'); set(groot, 'defaultLegendInterpreter','latex'); @@ -14,63 +14,127 @@ function plotLive(psi,Params,Transf,Observ) % 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 + switch SimulationMode + case 'ImaginaryTimeEvolution' + %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'); + 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,z)|^2$', 'Interpreter', 'latex', 'FontSize', 14) + + nexttile; + plotyz = pcolor(y,z,nyz'); + set(plotyz, 'EdgeColor', 'none'); + 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(y,z)|^2$', 'Interpreter', 'latex', 'FontSize', 14) + + nexttile; + plotxy = pcolor(x,y,nxy'); + set(plotxy, 'EdgeColor', 'none'); + 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) + + nexttile; + plot(-log10(Observ.residual),'-b') + 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_{tot}$', 'FontSize', 14); xlabel('Time steps', 'FontSize', 14); + title('Total Energy', 'FontSize', 14); + grid on + + nexttile; + plot(Observ.mucVec,'-b') + ylabel('$\mu$', 'FontSize', 14); xlabel('Time steps', 'FontSize', 14); + title('Chemical Potential', 'FontSize', 14); + grid on + + case 'RealTimeEvolution' - 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'); - 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,z)|^2$', 'Interpreter', 'latex', 'FontSize', 14) - - nexttile; - plotyz = pcolor(y,z,nyz'); - set(plotyz, 'EdgeColor', 'none'); - 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(y,z)|^2$', 'Interpreter', 'latex', 'FontSize', 14) - - nexttile; - plotxy = pcolor(x,y,nxy'); - set(plotxy, 'EdgeColor', 'none'); - 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) - - nexttile; - plot(-log10(Observ.residual),'-b') - 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_{tot}$', 'FontSize', 14); xlabel('Time steps', 'FontSize', 14); - title('Total Energy', 'FontSize', 14); - grid on - - nexttile; - plot(Observ.mucVec,'-b') - ylabel('$\mu$', 'FontSize', 14); xlabel('Time steps', 'FontSize', 14); - title('Chemical Potential', 'FontSize', 14); - grid on + %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'); + 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,z)|^2$', 'Interpreter', 'latex', 'FontSize', 14) + + nexttile; + plotyz = pcolor(y,z,nyz'); + set(plotyz, 'EdgeColor', 'none'); + 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(y,z)|^2$', 'Interpreter', 'latex', 'FontSize', 14) + + nexttile; + plotxy = pcolor(x,y,nxy'); + set(plotxy, 'EdgeColor', 'none'); + 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) + + nexttile; + plot(Observ.tVecPlot,Observ.NormVec,'-b') + ylabel('$N$', 'FontSize', 14); xlabel('Time steps', 'FontSize', 14); + title('Normalization', 'FontSize', 14); + grid on + + nexttile; + plot(Observ.tVecPlot,1-2*Observ.PCVec/pi,'-b') + ylabel('$C$', 'FontSize', 14); xlabel('Time steps', 'FontSize', 14); + title('Coherence', 'FontSize', 14); + grid on + + nexttile; + plot(Observ.tVecPlot,Observ.EVec,'-b') + ylabel('$E_{tot}$', 'FontSize', 14); xlabel('Time steps', 'FontSize', 14); + title('Energy', 'FontSize', 14); + grid on end \ No newline at end of file diff --git a/Dipolar-Gas-Simulator/+Scripts/run_locally.m b/Dipolar-Gas-Simulator/+Scripts/run_locally.m index 45c9e76..f359596 100644 --- a/Dipolar-Gas-Simulator/+Scripts/run_locally.m +++ b/Dipolar-Gas-Simulator/+Scripts/run_locally.m @@ -3,7 +3,7 @@ % Important: Run only sectionwise!! %% - Create Simulator, Potential and Calculator object with specified options - +% - Imaginary-Time OptionsStruct = struct; OptionsStruct.NumberOfAtoms = 40000; @@ -50,8 +50,79 @@ sim.Potential = pot.trap(); % + pot.repulsive_chopstick( Plotter.visualizeTrapPotential(sim.Potential,Params,Transf) %% - Plot initial wavefunction Plotter.visualizeWavefunction(psi,Params,Transf) -%% +%% Imaginary-Time followed by Real-time +% - Imaginary-Time +OptionsStruct = struct; + +OptionsStruct.NumberOfAtoms = 80000; +OptionsStruct.DipolarPolarAngle = deg2rad(0); +OptionsStruct.DipolarAzimuthAngle = 0; +OptionsStruct.ScatteringLength = 110; + +OptionsStruct.TrapFrequencies = [50, 20, 150]; +OptionsStruct.TrapPotentialType = 'Harmonic'; + +OptionsStruct.NumberOfGridPoints = [64, 128, 64]; +OptionsStruct.Dimensions = [50, 50, 20]; +OptionsStruct.UseApproximationForLHY = true; +OptionsStruct.IncludeDDICutOff = true; +OptionsStruct.CutoffType = 'CustomCylindrical'; +OptionsStruct.CustomCylindricalCutOffRadius = 20.0; +OptionsStruct.CustomCylindricalCutOffHeight = 5.0; + +OptionsStruct.SimulationMode = 'ImaginaryTimeEvolution'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution' | 'EnergyMinimization' +OptionsStruct.TimeStepSize = 5E-3; % in s +OptionsStruct.MinimumTimeStepSize = 1E-6; % in s +OptionsStruct.TimeCutOff = 1E5; % in s +OptionsStruct.EnergyTolerance = 5E-10; +OptionsStruct.ResidualTolerance = 1E-08; +OptionsStruct.NoiseScaleFactor = 0.010; +OptionsStruct.PlotLive = true; +OptionsStruct.JobNumber = 0; +OptionsStruct.RunOnGPU = false; +OptionsStruct.SaveData = false; +options = Helper.convertstruct2cell(OptionsStruct); + +sim = Simulator.DipolarGas(options{:}); +pot = Simulator.Potentials(options{:}); +sim.Potential = pot.trap(); + +%-% Run Simulation %-% +[Params, Transf, psi, V, VDk] = sim.run(); + +SaveDirectory = './Results/Data_3D/RealTimeDynamics'; +save(sprintf(strcat(SaveDirectory, '/psi_init.mat'),Params.njob),'psi','Transf','Params','VDk','V'); +close all +clear + +OptionsStruct.SimulationMode = 'RealTimeEvolution'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution' | 'EnergyMinimization' +OptionsStruct.EquilibrationTime = 20E-3; +OptionsStruct.QuenchTime = 30E-3; +OptionsStruct.HoldTime = 20E-3; +OptionsStruct.QuenchScatteringLength = true; +OptionsStruct.RotateDipoles = false; +OptionsStruct.FinalScatteringLength = 85; +OptionsStruct.FinalDipolarPolarAngle = deg2rad(90); +OptionsStruct.FinalDipolarAzimuthAngle = 0; +OptionsStruct.TimeStepSize = 1E-3; % in s +OptionsStruct.NoiseScaleFactor = 0.010; + +OptionsStruct.PlotLive = true; +OptionsStruct.JobNumber = 0; +OptionsStruct.RunOnGPU = false; +OptionsStruct.SaveData = true; +OptionsStruct.SaveDirectory = SaveDirectory; +options = Helper.convertstruct2cell(OptionsStruct); +clear OptionsStruct + +sim = Simulator.DipolarGas(options{:}); +pot = Simulator.Potentials(options{:}); +sim.Potential = pot.trap(); +%-% Run Simulation %-% +[Params, Transf, psi, V, VDk] = sim.run(); + +%% % 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) @@ -435,12 +506,16 @@ SaveDirectory = './Results/Data_3D/ApproximateLHY/BeyondSSD_Honeycomb'; JobNumber = 0; Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber) %% -SaveDirectory = './Results/Data_3D/TiltedDipoles15'; +SaveDirectory = './Results/Data_3D/TiltedDipoles0'; JobNumber = 0; Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber) %% +SaveDirectory = './Results/Data_3D/TiltedDipoles15'; +JobNumber = 1; +Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber) +%% SaveDirectory = './Results/Data_3D/TiltedDipoles30'; -JobNumber = 0; +JobNumber = 1; Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber) %% SaveDirectory = './Results/Data_3D/TiltedDipoles45'; diff --git a/Dipolar-Gas-Simulator/+Scripts/run_on_cluster.m b/Dipolar-Gas-Simulator/+Scripts/run_on_cluster.m index 1608c1e..1745af3 100644 --- a/Dipolar-Gas-Simulator/+Scripts/run_on_cluster.m +++ b/Dipolar-Gas-Simulator/+Scripts/run_on_cluster.m @@ -1,3 +1,4 @@ +%{ theta = 0; phi = 0; @@ -44,4 +45,52 @@ sim.Potential = pot.trap(); %-% Run Simulation %-% [Params, Transf, psi, V, VDk] = sim.run(); +%} + +%% - Aspect Ratio: 2.5 + +theta = 45; + +OptionsStruct = struct; + +OptionsStruct.NumberOfAtoms = 1E5; +OptionsStruct.DipolarPolarAngle = deg2rad(theta); +OptionsStruct.DipolarAzimuthAngle = 0; +OptionsStruct.ScatteringLength = 85; + +AspectRatio = 2.5; +HorizontalTrapFrequency = 125; +VerticalTrapFrequency = AspectRatio * HorizontalTrapFrequency; +OptionsStruct.TrapFrequencies = [HorizontalTrapFrequency, HorizontalTrapFrequency, VerticalTrapFrequency]; +OptionsStruct.TrapPotentialType = 'Harmonic'; + +OptionsStruct.NumberOfGridPoints = [128, 128, 64]; +OptionsStruct.Dimensions = [12, 12, 12]; +OptionsStruct.UseApproximationForLHY = true; +OptionsStruct.IncludeDDICutOff = true; +OptionsStruct.CutoffType = 'CustomCylindrical'; +OptionsStruct.CustomCylindricalCutOffRadius = 4.5; +OptionsStruct.CustomCylindricalCutOffHeight = 4.5; +OptionsStruct.SimulationMode = 'ImaginaryTimeEvolution'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution' +OptionsStruct.TimeStepSize = 1E-3; % in s +OptionsStruct.MinimumTimeStepSize = 2E-6; % in s +OptionsStruct.TimeCutOff = 2E6; % in s +OptionsStruct.EnergyTolerance = 5E-10; +OptionsStruct.ResidualTolerance = 1E-05; +OptionsStruct.NoiseScaleFactor = 0.01; + +OptionsStruct.PlotLive = false; +OptionsStruct.JobNumber = 2; +OptionsStruct.RunOnGPU = true; +OptionsStruct.SaveData = true; +OptionsStruct.SaveDirectory = sprintf('./Results/Data_3D/TiltedDipoles%s', strrep(num2str(round(rad2deg(OptionsStruct.DipolarPolarAngle),2)), '.', '_')); +options = Helper.convertstruct2cell(OptionsStruct); +clear OptionsStruct + +sim = Simulator.DipolarGas(options{:}); +pot = Simulator.Potentials(options{:}); +sim.Potential = pot.trap(); + +%-% Run Simulation %-% +[Params, Transf, psi, V, VDk] = sim.run(); \ No newline at end of file diff --git a/Dipolar-Gas-Simulator/+Simulator/@Calculator/Calculator.m b/Dipolar-Gas-Simulator/+Simulator/@Calculator/Calculator.m index f352bd4..824cc50 100644 --- a/Dipolar-Gas-Simulator/+Simulator/@Calculator/Calculator.m +++ b/Dipolar-Gas-Simulator/+Simulator/@Calculator/Calculator.m @@ -36,7 +36,10 @@ classdef Calculator < handle & matlab.mixin.Copyable p.KeepUnmatched = true; addParameter(p, 'CutOffType', this.CalculatorDefaults.CutOffType,... @(x) any(strcmpi(x,{'None', 'Cylindrical','CylindricalInfiniteZ', 'Spherical', 'CustomCylindrical'}))); - + addParameter(p, 'CustomCylindricalCutOffRadius', 1,... + @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); + addParameter(p, 'CustomCylindricalCutOffHeight', 1,... + @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); p.parse(varargin{:}); this.ChemicalPotential = this.CalculatorDefaults.ChemicalPotential; @@ -46,8 +49,8 @@ classdef Calculator < handle & matlab.mixin.Copyable this.PhaseCoherence = this.CalculatorDefaults.PhaseCoherence; this.TotalEnergy = this.CalculatorDefaults.TotalEnergy; this.CutOffType = p.Results.CutOffType; - this.CustomCylindricalCutOffRadius = this.CalculatorDefaults.CustomCylindricalCutOffRadius; - this.CustomCylindricalCutOffHeight = this.CalculatorDefaults.CustomCylindricalCutOffHeight; + this.CustomCylindricalCutOffRadius = p.Results.CustomCylindricalCutOffRadius; + this.CustomCylindricalCutOffHeight = p.Results.CustomCylindricalCutOffHeight; end function restoreDefaults(this) diff --git a/Dipolar-Gas-Simulator/+Simulator/@Calculator/calculatePhaseCoherence.m b/Dipolar-Gas-Simulator/+Simulator/@Calculator/calculatePhaseCoherence.m new file mode 100644 index 0000000..6ff5edf --- /dev/null +++ b/Dipolar-Gas-Simulator/+Simulator/@Calculator/calculatePhaseCoherence.m @@ -0,0 +1,19 @@ +function [PhaseC] = calculatePhaseCoherence(~,psi,Transf) + + norm = sum(sum(sum(abs(psi).^2,1),2),3)*Transf.dx*Transf.dy*Transf.dz; + psi = psi/sqrt(norm); + + 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))); + end + [~,minidx] = min(betas); + + 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; +end \ No newline at end of file diff --git a/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/DipolarGas.m b/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/DipolarGas.m index 288d30d..0f78ade 100644 --- a/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/DipolarGas.m +++ b/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/DipolarGas.m @@ -5,12 +5,20 @@ classdef DipolarGas < handle & matlab.mixin.Copyable DipolarPolarAngle; DipolarAzimuthAngle ScatteringLength; + FinalDipolarPolarAngle; + FinalDipolarAzimuthAngle + FinalScatteringLength; + EquilibrationTime; + QuenchTime; + HoldTime; TrapFrequencies; NumberOfGridPoints; Dimensions; Potential; SimulationMode; + QuenchScatteringLength; + RotateDipoles; TimeStepSize; MinimumTimeStepSize; TimeCutOff; @@ -23,6 +31,7 @@ classdef DipolarGas < handle & matlab.mixin.Copyable Calculator; SimulationParameters; + QuenchSettings; IncludeDDICutOff; UseApproximationForLHY; PlotLive; @@ -47,8 +56,14 @@ classdef DipolarGas < handle & matlab.mixin.Copyable @(x) assert(isnumeric(x) && isscalar(x) && (x > -2*pi) && (x < 2*pi))); addParameter(p, 'DipolarAzimuthAngle', pi/2,... @(x) assert(isnumeric(x) && isscalar(x) && (x > -2*pi) && (x < 2*pi))); + addParameter(p, 'FinalDipolarPolarAngle', 0,... + @(x) assert(isnumeric(x) && isscalar(x) && (x > -2*pi) && (x < 2*pi))); + addParameter(p, 'FinalDipolarAzimuthAngle', 0,... + @(x) assert(isnumeric(x) && isscalar(x) && (x > -2*pi) && (x < 2*pi))); addParameter(p, 'ScatteringLength', 120,... @(x) assert(isnumeric(x) && isscalar(x) && (x > -150) && (x < 150))); + addParameter(p, 'FinalScatteringLength', 88,... + @(x) assert(isnumeric(x) && isscalar(x) && (x > -150) && (x < 150))); addParameter(p, 'TrapFrequencies', 100 * ones(1,3),... @(x) assert(isnumeric(x) && isvector(x) && all(x > 0))); addParameter(p, 'NumberOfGridPoints', 128 * ones(1,3),... @@ -59,6 +74,12 @@ classdef DipolarGas < handle & matlab.mixin.Copyable @(x) assert(any(strcmpi(x,{'ImaginaryTimeEvolution','RealTimeEvolution', 'EnergyMinimization'})))); addParameter(p, 'CutOffType', 'Cylindrical',... @(x) assert(any(strcmpi(x,{'None', 'Cylindrical','CylindricalInfiniteZ', 'Spherical', 'CustomCylindrical'})))); + addParameter(p, 'EquilibrationTime', 20E-3,... + @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); + addParameter(p, 'QuenchTime', 30E-3,... + @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); + addParameter(p, 'HoldTime', 20E-3,... + @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); addParameter(p, 'TimeStepSize', 5E-4,... @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); addParameter(p, 'MinimumTimeStepSize', 1e-6,... @@ -75,7 +96,10 @@ classdef DipolarGas < handle & matlab.mixin.Copyable @(x) assert(any(strcmpi(x,{'HeavyBall','NonLinearCGD'})))); addParameter(p, 'MaxIterationsForGD', 100,... @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); - + addParameter(p, 'QuenchScatteringLength', false,... + @islogical); + addParameter(p, 'RotateDipoles', false,... + @islogical); addParameter(p, 'IncludeDDICutOff', true,... @islogical); addParameter(p, 'UseApproximationForLHY', false,... @@ -95,37 +119,44 @@ classdef DipolarGas < handle & matlab.mixin.Copyable p.parse(varargin{:}); - this.NumberOfAtoms = p.Results.NumberOfAtoms; - this.DipolarPolarAngle = p.Results.DipolarPolarAngle; - this.DipolarAzimuthAngle = p.Results.DipolarAzimuthAngle; - this.ScatteringLength = p.Results.ScatteringLength; - this.TrapFrequencies = p.Results.TrapFrequencies; - this.NumberOfGridPoints = p.Results.NumberOfGridPoints; - this.Dimensions = p.Results.Dimensions; - this.Potential = NaN; - this.SimulationMode = p.Results.SimulationMode; - this.TimeStepSize = p.Results.TimeStepSize; - this.MinimumTimeStepSize = p.Results.MinimumTimeStepSize; - this.TimeCutOff = p.Results.TimeCutOff; - this.EnergyTolerance = p.Results.EnergyTolerance; - this.ResidualTolerance = p.Results.ResidualTolerance; - this.NoiseScaleFactor = p.Results.NoiseScaleFactor; - this.GradientDescentMethod = p.Results.GradientDescentMethod; - this.MaxIterationsForGD = p.Results.MaxIterationsForGD; + this.NumberOfAtoms = p.Results.NumberOfAtoms; + this.DipolarPolarAngle = p.Results.DipolarPolarAngle; + this.DipolarAzimuthAngle = p.Results.DipolarAzimuthAngle; + this.ScatteringLength = p.Results.ScatteringLength; + this.FinalDipolarPolarAngle = p.Results.FinalDipolarPolarAngle; + this.FinalDipolarAzimuthAngle = p.Results.FinalDipolarAzimuthAngle; + this.FinalScatteringLength = p.Results.FinalScatteringLength; + this.EquilibrationTime = p.Results.EquilibrationTime; + this.QuenchTime = p.Results.QuenchTime; + this.HoldTime = p.Results.HoldTime; + this.TrapFrequencies = p.Results.TrapFrequencies; + this.NumberOfGridPoints = p.Results.NumberOfGridPoints; + this.Dimensions = p.Results.Dimensions; + this.Potential = NaN; + this.SimulationMode = p.Results.SimulationMode; + this.TimeStepSize = p.Results.TimeStepSize; + this.MinimumTimeStepSize = p.Results.MinimumTimeStepSize; + this.TimeCutOff = p.Results.TimeCutOff; + this.EnergyTolerance = p.Results.EnergyTolerance; + this.ResidualTolerance = p.Results.ResidualTolerance; + this.NoiseScaleFactor = p.Results.NoiseScaleFactor; + this.GradientDescentMethod = p.Results.GradientDescentMethod; + this.MaxIterationsForGD = p.Results.MaxIterationsForGD; + this.QuenchScatteringLength = p.Results.QuenchScatteringLength; + this.RotateDipoles = p.Results.RotateDipoles; + this.IncludeDDICutOff = p.Results.IncludeDDICutOff; + this.UseApproximationForLHY = p.Results.UseApproximationForLHY; + this.PlotLive = p.Results.PlotLive; + this.JobNumber = p.Results.JobNumber; + this.RunOnGPU = p.Results.RunOnGPU; + this.DebugMode = p.Results.DebugMode; + this.DoSave = p.Results.SaveData; + this.SaveDirectory = p.Results.SaveDirectory; - this.IncludeDDICutOff = p.Results.IncludeDDICutOff; - this.UseApproximationForLHY = p.Results.UseApproximationForLHY; - this.PlotLive = p.Results.PlotLive; - this.JobNumber = p.Results.JobNumber; - this.RunOnGPU = p.Results.RunOnGPU; - this.DebugMode = p.Results.DebugMode; - this.DoSave = p.Results.SaveData; - this.SaveDirectory = p.Results.SaveDirectory; - - this.Calculator = Simulator.Calculator('CutOffType', p.Results.CutOffType); - - this.SimulationParameters = this.setupParameters(); + this.Calculator = Simulator.Calculator(varargin{:}); + this.SimulationParameters = this.setupParameters(); + end end diff --git a/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/initialize.m b/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/initialize.m index 0d16c30..16743fe 100644 --- a/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/initialize.m +++ b/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/initialize.m @@ -15,7 +15,7 @@ function [psi,V,VDk] = initialize(this,Params,Transf,TransfRad) [VDk, DDICutOffIncluded, CutOffType] = deal(loadedData.VDk, loadedData.DDICutOffIncluded, loadedData.CutOffType); end - if isempty(VDk) || ~isequal(size(VDk), this.NumberOfGridPoints) || this.IncludeDDICutOff ~= DDICutOffIncluded + if isempty(VDk) || ~isequal(size(VDk), this.NumberOfGridPoints) || this.IncludeDDICutOff ~= DDICutOffIncluded || strcmp(this.Calculator.CutOffType, 'CustomCylindrical') % Calculate VDk if necessary VDk = this.Calculator.calculateVDk(Params, Transf, TransfRad, this.IncludeDDICutOff); DDICutOffIncluded = this.IncludeDDICutOff; diff --git a/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/propagateWavefunction.m b/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/propagateWavefunction.m index 87198bf..0f4c9f6 100644 --- a/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/propagateWavefunction.m +++ b/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/propagateWavefunction.m @@ -3,7 +3,7 @@ function [psi] = propagateWavefunction(this,psi,Params,Transf,VDk,V,t_idx,Observ switch this.SimulationMode case 'ImaginaryTimeEvolution' - dt = -1j*abs(this.TimeStepSize); % Imaginary Time + dt = -1j*abs(Params.dt); % Imaginary Time KEop = 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2); Observ.residual = 1; @@ -12,7 +12,7 @@ function [psi] = propagateWavefunction(this,psi,Params,Transf,VDk,V,t_idx,Observ muchem = this.Calculator.calculateChemicalPotential(psi,Params,Transf,VDk,V); if this.PlotLive - Plotter.plotLive(psi,Params,Transf,Observ) + Plotter.plotLive(psi,Params,Transf,Observ,this.SimulationMode) drawnow end @@ -59,11 +59,13 @@ function [psi] = propagateWavefunction(this,psi,Params,Transf,VDk,V,t_idx,Observ Observ.res_idx = Observ.res_idx + 1; if this.PlotLive - Plotter.plotLive(psi,Params,Transf,Observ) + Plotter.plotLive(psi,Params,Transf,Observ,this.SimulationMode) drawnow end - - save(sprintf(strcat(this.SaveDirectory, '/Run_%03i/psi_gs.mat'),Params.njob),'psi','muchem','Observ','t_idx','Transf','Params','VDk','V'); + + if this.DoSave + save(sprintf(strcat(this.SaveDirectory, '/Run_%03i/psi_gs.mat'),Params.njob),'psi','muchem','Observ','t_idx','Transf','Params','VDk','V'); + end %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); @@ -96,20 +98,41 @@ function [psi] = propagateWavefunction(this,psi,Params,Transf,VDk,V,t_idx,Observ Observ.res_idx = Observ.res_idx + 1; disp('Run completed!'); - disp('Saving data...'); - save(sprintf(strcat(this.SaveDirectory, '/Run_%03i/psi_gs.mat'),Params.njob),'psi','muchem','Observ','t_idx','Transf','Params','VDk','V'); - disp('Save complete!'); + if this.DoSave + disp('Saving data...'); + save(sprintf(strcat(this.SaveDirectory, '/Run_%03i/psi_gs.mat'),Params.njob),'psi','muchem','Observ','t_idx','Transf','Params','VDk','V'); + disp('Save complete!'); + end case 'RealTimeEvolution' - dt = abs(this.TimeStepSize); + dt = abs(Params.dt); - KEop= 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2); + KEop = 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2); muchem = this.Calculator.calculateChemicalPotential(psi,Params,Transf,VDk,V); - while t_idx < Params.sim_time_cut_off + if this.PlotLive + Plotter.plotLive(psi,Params,Transf,Observ,this.SimulationMode) + drawnow + end + while t_idx < this.QuenchSettings.tSteps + + %Time + tVal = this.QuenchSettings.tVec(t_idx); + + if this.QuenchScatteringLength + Params.gs = this.QuenchSettings.gs_vec(t_idx); + Params.gammaQF = this.QuenchSettings.gammaQF_vec(t_idx); + end + + if this.RotateDipoles + Params.theta = this.QuenchSettings.theta_vec(t_idx); + Params.phi = this.QuenchSettings.phi_vec(t_idx); + VDk = this.Calculator.calculateVDk(Params, Transf, TransfRad, this.IncludeDDICutOff); + end + % Parameters at time t %kin @@ -119,7 +142,7 @@ function [psi] = propagateWavefunction(this,psi,Params,Transf,VDk,V,t_idx,Observ %DDI frho = fftn(abs(psi).^2); - Phi = real(ifftn(frho.*VDk)); + Phi = real(ifftn(frho.*VDk)); %Real-space psi = psi.*exp(-1i*dt*(V + Params.gs*abs(psi).^2 + Params.gammaQF*abs(psi).^3 + Params.gdd*Phi - muchem)); @@ -131,7 +154,7 @@ function [psi] = propagateWavefunction(this,psi,Params,Transf,VDk,V,t_idx,Observ muchem = this.Calculator.calculateChemicalPotential(psi,Params,Transf,VDk,V); - if mod(t_idx,1000)==0 + if mod(t_idx,this.QuenchSettings.saveinterval)==0 %Change in Normalization Norm = sum(abs(psi(:)).^2)*Transf.dx*Transf.dy*Transf.dz; % normalisation Observ.NormVec = [Observ.NormVec Norm]; @@ -140,17 +163,28 @@ function [psi] = propagateWavefunction(this,psi,Params,Transf,VDk,V,t_idx,Observ 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); + Observ.PCVec = [Observ.PCVec PhaseC]; + Observ.tVecPlot = [Observ.tVecPlot tVal]; Observ.res_idx = Observ.res_idx + 1; + + if this.PlotLive + Plotter.plotLive(psi,Params,Transf,Observ,this.SimulationMode) + drawnow + end - save(sprintf('./Data/Run_%03i/TimeEvolution/psi_%i.mat',Params.njob,Observ.res_idx),'psi','muchem','Observ','t_idx'); + if this.DoSave + save(sprintf(strcat(this.SaveDirectory, '/Run_%03i/psi_%i.mat'), Params.njob, Observ.res_idx),'psi','muchem','Observ','t_idx'); + end end if any(isnan(psi(:))) disp('NaNs encountered!') break end - t_idx=t_idx+1; + t_idx = t_idx+1; end %Change in Normalization @@ -170,12 +204,15 @@ function [psi] = propagateWavefunction(this,psi,Params,Transf,VDk,V,t_idx,Observ Observ.res_idx = Observ.res_idx + 1; 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!'); + if this.DoSave + disp('Saving data...'); + save(sprintf(strcat(this.SaveDirectory, '/Run_%03i/psi_%i.mat'),Params.njob,Observ.res_idx),'psi','muchem','Observ','t_idx'); + disp('Save complete!'); + end + otherwise - disp('Choose a valid DDI cutoff type!') + disp('Choose a valid method!') return end end \ No newline at end of file diff --git a/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/run.m b/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/run.m index 0393059..5d58c0e 100644 --- a/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/run.m +++ b/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/run.m @@ -3,6 +3,10 @@ function [Params, Transf, psi,V,VDk] = run(this) [Params] = this.setupParameters(); this.SimulationParameters = Params; + if strcmp(this.SimulationMode, 'RealTimeEvolution') + this.QuenchSettings = this.setupQuenchSettings(Params); + end + % --- Set up spatial grids and transforms --- [Transf] = this.setupSpace(Params); [TransfRad] = this.setupSpaceRadial(Params); diff --git a/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/setupParameters.m b/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/setupParameters.m index 6c3e49a..741fb8d 100644 --- a/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/setupParameters.m +++ b/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/setupParameters.m @@ -46,7 +46,7 @@ function [Params] = setupParameters(this) 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.dt = this.TimeStepSize; Params.nsf = this.NoiseScaleFactor; Params.njob = this.JobNumber; diff --git a/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/setupQuenchSettings.m b/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/setupQuenchSettings.m new file mode 100644 index 0000000..35f35fd --- /dev/null +++ b/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/setupQuenchSettings.m @@ -0,0 +1,43 @@ +function Quench = setupQuenchSettings(this, Params) + + Quench.eq_time = this.EquilibrationTime; % Equilibration time in s + Quench.quench_time = this.QuenchTime; % Quench time in s + Quench.hold_time = this.HoldTime; % Hold time in s + + Quench.as_final = this.FinalScatteringLength*Params.a0; % Final scattering length in a0 (initial is Params.as) + Quench.theta_final = this.FinalDipolarPolarAngle; % pi/2 dipoles along x, theta=0 dipoles along z + Quench.phi_final = this.FinalDipolarAzimuthAngle; + + % --- Save steps + Quench.saveinterval = 200; + + % --- Constructing quench vectors + eq_time = Quench.eq_time*Params.w0; + quench_time = Quench.quench_time*Params.w0; + hold_time = Quench.hold_time*Params.w0; + + tVec_Eq = 0:Params.dt:eq_time; + tVec_Q = (eq_time+Params.dt):Params.dt:(eq_time+quench_time); + tVec_H = (eq_time+quench_time+Params.dt):Params.dt:(eq_time+quench_time+hold_time); + + Quench.tSteps = length(tVec_Eq) + length(tVec_Q) + length(tVec_H); + Quench.as_vec = [Params.as*ones(1, length(tVec_Eq)), linspace(Params.as,Quench.as_final,length(tVec_Q)), Quench.as_final*ones(1, length(tVec_H))]; + Quench.gs_vec = 4*pi*Quench.as_vec/Params.l0; + Quench.theta_vec = [Params.theta*ones(1, length(tVec_Eq)), linspace(Params.theta,Quench.theta_final,length(tVec_Q)), Quench.theta_final*ones(1, length(tVec_H))]; + Quench.phi_vec = [Params.phi*ones(1, length(tVec_Eq)), linspace(Params.phi,Quench.phi_final,length(tVec_Q)), Quench.phi_final*ones(1, length(tVec_H))]; + Quench.tVec = [tVec_Eq, tVec_Q, tVec_H]; + + % Quantum fluctuations need to be calculated over the quench + eps_dd_vec = Params.add./Quench.as_vec; + + if this.UseApproximationForLHY + Q5_vec = 1 + ((3*eps_dd_vec.^2)/2); + else + yeps_vec = (1-eps_dd_vec)./(3*eps_dd_vec); + Q5_vec = (3*eps_dd_vec).^(5/2).*( (8+26*yeps_vec+33*yeps_vec.^2).*sqrt(1+yeps_vec) + 15*yeps_vec.^3.*log((1+sqrt(1+yeps_vec))./sqrt(yeps_vec)) )/48; + Q5_vec = real(Q5_vec); + Q5_vec(eps_dd_vec == 0) = 1; + Q5_vec(eps_dd_vec == 1) = 3*sqrt(3)/2; + end + Quench.gammaQF_vec = 128/3*sqrt(pi*(Quench.as_vec/Params.l0).^5).*Q5_vec; +end \ No newline at end of file