Added functionality for real-time dynamics - quench of scattering length and rotation of dipoles.
This commit is contained in:
parent
6566f53559
commit
178350bdc2
@ -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,6 +14,8 @@ function plotLive(psi,Params,Transf,Observ)
|
||||
% Compute probability density |psi|^2
|
||||
n = abs(psi).^2;
|
||||
|
||||
switch SimulationMode
|
||||
case 'ImaginaryTimeEvolution'
|
||||
%Plotting
|
||||
figure(1);
|
||||
clf
|
||||
@ -73,4 +75,66 @@ function plotLive(psi,Params,Transf,Observ)
|
||||
ylabel('$\mu$', 'FontSize', 14); xlabel('Time steps', 'FontSize', 14);
|
||||
title('Chemical Potential', 'FontSize', 14);
|
||||
grid on
|
||||
|
||||
case 'RealTimeEvolution'
|
||||
|
||||
%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
|
||||
@ -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';
|
||||
|
||||
@ -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();
|
||||
@ -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)
|
||||
|
||||
@ -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
|
||||
@ -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,...
|
||||
@ -99,6 +123,12 @@ classdef DipolarGas < handle & matlab.mixin.Copyable
|
||||
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;
|
||||
@ -112,7 +142,8 @@ classdef DipolarGas < handle & matlab.mixin.Copyable
|
||||
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;
|
||||
@ -122,7 +153,7 @@ classdef DipolarGas < handle & matlab.mixin.Copyable
|
||||
this.DoSave = p.Results.SaveData;
|
||||
this.SaveDirectory = p.Results.SaveDirectory;
|
||||
|
||||
this.Calculator = Simulator.Calculator('CutOffType', p.Results.CutOffType);
|
||||
this.Calculator = Simulator.Calculator(varargin{:});
|
||||
|
||||
this.SimulationParameters = this.setupParameters();
|
||||
|
||||
|
||||
@ -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;
|
||||
|
||||
@ -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
|
||||
|
||||
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,19 +98,40 @@ function [psi] = propagateWavefunction(this,psi,Params,Transf,VDk,V,t_idx,Observ
|
||||
Observ.res_idx = Observ.res_idx + 1;
|
||||
|
||||
disp('Run completed!');
|
||||
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);
|
||||
|
||||
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
|
||||
|
||||
@ -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];
|
||||
@ -141,10 +164,21 @@ function [psi] = propagateWavefunction(this,psi,Params,Transf,VDk,V,t_idx,Observ
|
||||
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;
|
||||
|
||||
save(sprintf('./Data/Run_%03i/TimeEvolution/psi_%i.mat',Params.njob,Observ.res_idx),'psi','muchem','Observ','t_idx');
|
||||
if this.PlotLive
|
||||
Plotter.plotLive(psi,Params,Transf,Observ,this.SimulationMode)
|
||||
drawnow
|
||||
end
|
||||
|
||||
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!')
|
||||
@ -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!');
|
||||
if this.DoSave
|
||||
disp('Saving data...');
|
||||
save(sprintf('./Data/Run_%03i/TimeEvolution/psi_%i.mat',Params.njob,Observ.res_idx),'psi','muchem','Observ','t_idx');
|
||||
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
|
||||
@ -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);
|
||||
|
||||
@ -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;
|
||||
|
||||
@ -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
|
||||
Loading…
Reference in New Issue
Block a user