From 118406eed056fab019ead740e3570eb50ce75491 Mon Sep 17 00:00:00 2001 From: Karthik Chandrashekara Date: Wed, 13 Nov 2024 18:36:50 +0100 Subject: [PATCH] Adding a 2D GPE solver using the Variational method. --- .../+Variational2D/@Calculator/Calculator.m | 78 ++++++++ .../@Calculator/calculateChemicalPotential.m | 32 ++++ .../@Calculator/calculateEnergyComponents.m | 29 +++ .../calculateNormalizedResiduals.m | 29 +++ .../@Calculator/calculateTotalEnergy.m | 34 ++++ .../@Calculator/calculateVDkWithCutoff.m | 29 +++ .../@Calculator/calculateVariationalEnergy.m | 38 ++++ .../+Variational2D/@Calculator/find_nk_fit.m | 24 +++ .../+Variational2D/@DipolarGas/DipolarGas.m | 177 ++++++++++++++++++ .../+Variational2D/@DipolarGas/initialize.m | 26 +++ .../@DipolarGas/propagateWavefunction.m | 117 ++++++++++++ .../+Variational2D/@DipolarGas/run.m | 72 +++++++ .../@DipolarGas/setupParameters.m | 109 +++++++++++ .../+Variational2D/@DipolarGas/setupSpace.m | 33 ++++ .../@DipolarGas/setupWavefunction.m | 21 +++ 15 files changed, 848 insertions(+) create mode 100644 Dipolar-Gas-Simulator/+Variational2D/@Calculator/Calculator.m create mode 100644 Dipolar-Gas-Simulator/+Variational2D/@Calculator/calculateChemicalPotential.m create mode 100644 Dipolar-Gas-Simulator/+Variational2D/@Calculator/calculateEnergyComponents.m create mode 100644 Dipolar-Gas-Simulator/+Variational2D/@Calculator/calculateNormalizedResiduals.m create mode 100644 Dipolar-Gas-Simulator/+Variational2D/@Calculator/calculateTotalEnergy.m create mode 100644 Dipolar-Gas-Simulator/+Variational2D/@Calculator/calculateVDkWithCutoff.m create mode 100644 Dipolar-Gas-Simulator/+Variational2D/@Calculator/calculateVariationalEnergy.m create mode 100644 Dipolar-Gas-Simulator/+Variational2D/@Calculator/find_nk_fit.m create mode 100644 Dipolar-Gas-Simulator/+Variational2D/@DipolarGas/DipolarGas.m create mode 100644 Dipolar-Gas-Simulator/+Variational2D/@DipolarGas/initialize.m create mode 100644 Dipolar-Gas-Simulator/+Variational2D/@DipolarGas/propagateWavefunction.m create mode 100644 Dipolar-Gas-Simulator/+Variational2D/@DipolarGas/run.m create mode 100644 Dipolar-Gas-Simulator/+Variational2D/@DipolarGas/setupParameters.m create mode 100644 Dipolar-Gas-Simulator/+Variational2D/@DipolarGas/setupSpace.m create mode 100644 Dipolar-Gas-Simulator/+Variational2D/@DipolarGas/setupWavefunction.m diff --git a/Dipolar-Gas-Simulator/+Variational2D/@Calculator/Calculator.m b/Dipolar-Gas-Simulator/+Variational2D/@Calculator/Calculator.m new file mode 100644 index 0000000..009fe03 --- /dev/null +++ b/Dipolar-Gas-Simulator/+Variational2D/@Calculator/Calculator.m @@ -0,0 +1,78 @@ +classdef Calculator < handle & matlab.mixin.Copyable + + properties (Access = private) + + CalculatorDefaults = struct('ChemicalPotential', 1, ... + 'EnergyComponents', 1, ... + 'NormalizedResiduals', 1, ... + 'TotalEnergy', 1); + + end + + properties (Access = public) + ChemicalPotential; + EnergyComponents; + NormalizedResiduals; + TotalEnergy; + end + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %- Methods + + methods + function this = Calculator(varargin) + + this.ChemicalPotential = this.CalculatorDefaults.ChemicalPotential; + this.EnergyComponents = this.CalculatorDefaults.EnergyComponents; + this.NormalizedResiduals = this.CalculatorDefaults.NormalizedResiduals; + this.TotalEnergy = this.CalculatorDefaults.TotalEnergy; + end + + function restoreDefaults(this) + this.ChemicalPotential = this.CalculatorDefaults.ChemicalPotential; + this.EnergyComponents = this.CalculatorDefaults.EnergyComponents; + this.NormalizedResiduals = this.CalculatorDefaults.NormalizedResiduals; + this.TotalEnergy = this.CalculatorDefaults.TotalEnergy; + end + + end % + + % methods + + % end % - setters and getters + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %- Methods + + methods(Access = protected) + function cp = copyElement(this) + % Shallow copy object + cp = copyElement@matlab.mixin.Copyable(this); + + % Forces the setter to redefine the function handles to the new copied object + + pl = properties(this); + for k = 1:length(pl) + sc = superclasses(this.(pl{k})); + if any(contains(sc,{'matlab.mixin.Copyable'})) + cp.(pl{k}) = this.(pl{k}).copy(); + end + end + end + end + + methods (Static) + + % Creates an Instance of Class, ensures singleton behaviour (that there + % can only be one Instance of this class + function singleObj = getInstance(varargin) + % Creates an Instance of Class, ensures singleton behaviour + persistent localObj; + if isempty(localObj) || ~isvalid(localObj) + localObj = Simulator.Calculator(varargin{:}); + end + singleObj = localObj; + end + end + +end diff --git a/Dipolar-Gas-Simulator/+Variational2D/@Calculator/calculateChemicalPotential.m b/Dipolar-Gas-Simulator/+Variational2D/@Calculator/calculateChemicalPotential.m new file mode 100644 index 0000000..0a4cb81 --- /dev/null +++ b/Dipolar-Gas-Simulator/+Variational2D/@Calculator/calculateChemicalPotential.m @@ -0,0 +1,32 @@ +function muchem = calculateChemicalPotential(~,psi,Params,VParams,Transf,VDk,V) + + g_eff = Params.gs*VParams.nu/(2^(1+1/VParams.nu)*VParams.ell*gamma(1/VParams.nu)); + gamma_eff = Params.gammaQF*2^(1/VParams.nu-1.5)*5^(-1/VParams.nu)*VParams.ell*gamma(1+1/VParams.nu)*( VParams.nu/(VParams.ell*gamma(1/VParams.nu)) )^(5/2); + EVar = VParams.nu^2*gamma(2-1/VParams.nu)/(8*VParams.ell^2*gamma(1/VParams.nu)) + 0.5*Params.gz*VParams.ell^2*gamma(3/VParams.nu)/gamma(1/VParams.nu); + + % Parameters + normfac = Params.Lx*Params.Ly/numel(psi); + KEop = 0.5*(Transf.KX.^2+Transf.KY.^2); + + % DDIs + frho = fftn(abs(psi).^2); + Phi = real(ifftn(frho.*VDk)); + Eddi = (Params.gdd*Phi.*abs(psi).^2)/(sqrt(2*pi)*VParams.ell_eff); + + %Kinetic energy + Ekin = KEop.*abs(fftn(psi)*normfac).^2; + Ekin = trapz(Ekin(:))*Transf.dkx*Transf.dky/(2*pi)^2; + + %Potential energy + Epot = V.*abs(psi).^2; + + %Contact interactions + Eint = g_eff*abs(psi).^4; + + %Quantum fluctuations + Eqf = gamma_eff*abs(psi).^5; + + %Total energy + muchem = Ekin + EVar*Params.N + trapz(Epot(:) + Eint(:) + Eddi(:) + Eqf(:))*Transf.dx*Transf.dy; % + muchem = muchem / Params.N; %Only use if psi is normalized to N +end \ No newline at end of file diff --git a/Dipolar-Gas-Simulator/+Variational2D/@Calculator/calculateEnergyComponents.m b/Dipolar-Gas-Simulator/+Variational2D/@Calculator/calculateEnergyComponents.m new file mode 100644 index 0000000..1011665 --- /dev/null +++ b/Dipolar-Gas-Simulator/+Variational2D/@Calculator/calculateEnergyComponents.m @@ -0,0 +1,29 @@ +function E = calculateEnergyComponents(~,psi,Params,Transf,VDk,V,VParams) + + % Parameters + KEop = 0.5*(Transf.KX.^2+Transf.KY.^2); + normfac = Params.Lx*Params.Ly/numel(psi); + + % DDIs + frho = fftn(abs(psi).^2); + Phi = (ifftn(frho.*VDk)); + Eddi = 0.5*Params.gdd*Phi.*abs(psi).^2/(sqrt(2*pi)*VParams.ell);% + E.Eddi = trapz(Eddi(:))*Transf.dx*Transf.dy; + + % Kinetic energy + Ekin = KEop.*abs(fftn(psi)*normfac).^2; + E.Ekin = trapz(Ekin(:))*Transf.dkx*Transf.dky/(2*pi)^2 + 0.25*(1/VParams.ell^2)*Params.N; + + % Potential energy + Epot = V.*abs(psi).^2; + E.Epot = trapz(Epot(:))*Transf.dx*Transf.dy + 0.25*Params.gz*VParams.ell^2*Params.N; + + % Contact interactions + Eint = 0.5*Params.gs*abs(psi).^4/(sqrt(2*pi)*VParams.ell); + E.Eint = trapz(Eint(:))*Transf.dx*Transf.dy; + + % Quantum fluctuations + Eqf = 0.4*Params.gammaQF*abs(psi).^5*sqrt(2/(5*VParams.ell^3*pi^(3/2))); + E.Eqf = trapz(Eqf(:))*Transf.dx*Transf.dy; + +end diff --git a/Dipolar-Gas-Simulator/+Variational2D/@Calculator/calculateNormalizedResiduals.m b/Dipolar-Gas-Simulator/+Variational2D/@Calculator/calculateNormalizedResiduals.m new file mode 100644 index 0000000..97523fd --- /dev/null +++ b/Dipolar-Gas-Simulator/+Variational2D/@Calculator/calculateNormalizedResiduals.m @@ -0,0 +1,29 @@ +function res = calculateNormalizedResiduals(~,psi,Params,VParams,Transf,VDk,V,muchem) + + g_eff = Params.gs*VParams.nu/(2^(1+1/VParams.nu)*VParams.ell*gamma(1/VParams.nu)); + gamma_eff = Params.gammaQF*2^(1/VParams.nu-1.5)*5^(-1/VParams.nu)*VParams.ell*gamma(1+1/VParams.nu)*( VParams.nu/(VParams.ell*gamma(1/VParams.nu)) )^(5/2); + EVar = VParams.nu^2*gamma(2-1/VParams.nu)/(8*VParams.ell^2*gamma(1/VParams.nu)) + 0.5*Params.gz*VParams.ell^2*gamma(3/VParams.nu)/gamma(1/VParams.nu); + + KEop = 0.5*(Transf.KX.^2+Transf.KY.^2); + + % DDIs + frho = fftn(abs(psi).^2); + Phi = real(ifftn(frho.*VDk)); + Eddi = Params.gdd*Phi.*psi/(sqrt(2*pi)*VParams.ell_eff); + + % Kinetic energy + Ekin = ifftn(KEop.*fftn(psi)); + + % Potential energy + Epot = V.*psi; + + % Contact interactions + Eint = g_eff*abs(psi).^2.*psi; + + % Quantum fluctuations + Eqf = gamma_eff*abs(psi).^3.*psi; + + % Total energy + res = trapz(abs(Ekin(:) + EVar*psi(:) + Epot(:) + Eint(:) + Eddi(:) + Eqf(:) - muchem*psi(:))*Transf.dx*Transf.dy)/trapz(abs(muchem*psi(:))*Transf.dx*Transf.dy); + +end \ No newline at end of file diff --git a/Dipolar-Gas-Simulator/+Variational2D/@Calculator/calculateTotalEnergy.m b/Dipolar-Gas-Simulator/+Variational2D/@Calculator/calculateTotalEnergy.m new file mode 100644 index 0000000..a36640f --- /dev/null +++ b/Dipolar-Gas-Simulator/+Variational2D/@Calculator/calculateTotalEnergy.m @@ -0,0 +1,34 @@ +function E = calculateTotalEnergy(~,psi,Params,VParams,Transf,VDk,V) + + g_eff = Params.gs*VParams.nu/(2^(1+1/VParams.nu)*VParams.ell*gamma(1/VParams.nu)); + gamma_eff = Params.gammaQF*2^(1/VParams.nu-1.5)*5^(-1/VParams.nu)*VParams.ell*gamma(1+1/VParams.nu)*( VParams.nu/(VParams.ell*gamma(1/VParams.nu)) )^(5/2); + EVar = VParams.nu^2*gamma(2-1/VParams.nu)/(8*VParams.ell^2*gamma(1/VParams.nu)) + 0.5*Params.gz*VParams.ell^2*gamma(3/VParams.nu)/gamma(1/VParams.nu); + + % Parameters + KEop = 0.5*(Transf.KX.^2+Transf.KY.^2); + normfac = Params.Lx*Params.Ly/numel(psi); + + % DDIs + frho = fftn(abs(psi).^2); + Phi = real(ifftn(frho.*VDk)); + Eddi = 0.5*Params.gdd*Phi.*abs(psi).^2/(sqrt(2*pi)*VParams.ell_eff);% + Eddi = trapz(Eddi(:))*Transf.dx*Transf.dy; + + % Kinetic energy + Ekin = KEop.*abs(fftn(psi)*normfac).^2; + Ekin = trapz(Ekin(:))*Transf.dkx*Transf.dky/(2*pi)^2; + + % Potential energy + Epot = V.*abs(psi).^2; + Epot = trapz(Epot(:))*Transf.dx*Transf.dy; + + % Contact interactions + Eint = 0.5*g_eff*abs(psi).^4; + Eint = trapz(Eint(:))*Transf.dx*Transf.dy; + + % Quantum fluctuations + Eqf = 0.4*gamma_eff*abs(psi).^5; + Eqf = trapz(Eqf(:))*Transf.dx*Transf.dy; + + E = EVar*Params.N + Ekin + Epot + Eint + Eddi + Eqf; +end \ No newline at end of file diff --git a/Dipolar-Gas-Simulator/+Variational2D/@Calculator/calculateVDkWithCutoff.m b/Dipolar-Gas-Simulator/+Variational2D/@Calculator/calculateVDkWithCutoff.m new file mode 100644 index 0000000..e7573ad --- /dev/null +++ b/Dipolar-Gas-Simulator/+Variational2D/@Calculator/calculateVDkWithCutoff.m @@ -0,0 +1,29 @@ +function VDk = calculateVDkWithCutoff(~, rcut, Transf, Params, VParams) +% == Calculating the DDI potential in Fourier space with appropriate cutoff == % + + % Interaction in K space + QX = Transf.KX*VParams.ell_eff/sqrt(2); + QY = Transf.KY*VParams.ell_eff/sqrt(2); + + Qsq = QX.^2 + QY.^2; + absQ = sqrt(Qsq); + QDsq = QX.^2*cos(Params.eta)^2 + QY.^2*sin(Params.eta)^2; + + %Bare interaction + Fpar = -1 + 3*sqrt(pi)*QDsq.*erfcx(absQ)./absQ; + Fperp = 2 - 3*sqrt(pi).*absQ.*erfcx(absQ); + Fpar(absQ == 0) = -1; %Fpar(isinf(Fpar)) = -1; + + %Full DDI + VDk = (Fpar*sin(Params.alpha)^2 + Fperp*cos(Params.alpha)^2); + + %Implementing a cutoff: + VDr = ifftn(VDk); + VDr = fftshift(VDr); + + R = sqrt(Transf.X.^2+Transf.Y.^2) < 0.9*rcut; + VDr = VDr.*double(R); + VDr = ifftshift(VDr); + VDk = fftn(VDr); + +end \ No newline at end of file diff --git a/Dipolar-Gas-Simulator/+Variational2D/@Calculator/calculateVariationalEnergy.m b/Dipolar-Gas-Simulator/+Variational2D/@Calculator/calculateVariationalEnergy.m new file mode 100644 index 0000000..fa7f170 --- /dev/null +++ b/Dipolar-Gas-Simulator/+Variational2D/@Calculator/calculateVariationalEnergy.m @@ -0,0 +1,38 @@ +function E = calculateVariationalEnergy(~, psi, Params, VarArray, Transf, VDk, V) + + VParams.ell = VarArray(1); + VParams.nu = VarArray(2); + + g_eff = Params.gs*VParams.nu/(2^(1+1/VParams.nu)*VParams.ell*gamma(1/VParams.nu)); + gamma_eff = Params.gammaQF*2^(1/VParams.nu-1.5)*5^(-1/VParams.nu)*VParams.ell*gamma(1+1/VParams.nu)*( VParams.nu/(VParams.ell*gamma(1/VParams.nu)) )^(5/2); + EVar = VParams.nu^2*gamma(2-1/VParams.nu)/(8*VParams.ell^2*gamma(1/VParams.nu)) + 0.5*Params.gz*VParams.ell^2*gamma(3/VParams.nu)/gamma(1/VParams.nu); + + % Parameters + KEop = 0.5*(Transf.KX.^2+Transf.KY.^2); + normfac = Params.Lx*Params.Ly/numel(psi); + + % DDIs + [VParams] = find_nk_fit(VParams); % Not totally sure this should be here + frho = fftn(abs(psi).^2); + Phi = real(ifftn(frho.*VDk)); + Eddi = 0.5*Params.gdd*Phi.*abs(psi).^2/(sqrt(2*pi)*VParams.ell_eff);% + Eddi = sum(Eddi(:))*Transf.dx*Transf.dy; + + % Kinetic energy + Ekin = KEop.*abs(fftn(psi)*normfac).^2; + Ekin = trapz(Ekin(:))*Transf.dkx*Transf.dky/(2*pi)^2; + + % Potential energy + Epot = V.*abs(psi).^2; + Epot = trapz(Epot(:))*Transf.dx*Transf.dy; + + % Contact interactions + Eint = 0.5*g_eff*abs(psi).^4; + Eint = trapz(Eint(:))*Transf.dx*Transf.dy; + + % Quantum fluctuations + Eqf = 0.4*gamma_eff*abs(psi).^5; + Eqf = trapz(Eqf(:))*Transf.dx*Transf.dy; + + E = EVar*Params.N + Ekin + Epot + Eint + Eddi + Eqf; +end \ No newline at end of file diff --git a/Dipolar-Gas-Simulator/+Variational2D/@Calculator/find_nk_fit.m b/Dipolar-Gas-Simulator/+Variational2D/@Calculator/find_nk_fit.m new file mode 100644 index 0000000..086c94d --- /dev/null +++ b/Dipolar-Gas-Simulator/+Variational2D/@Calculator/find_nk_fit.m @@ -0,0 +1,24 @@ +function [VParams] = find_nk_fit(VParams) + + Lz = 20; + Nz = 4000; + normfac = Lz/Nz; + + z = linspace(-0.5*Lz,0.5*Lz-Lz/Nz,Nz); + dz = abs(z(2)-z(1)); + Kmax = pi*Nz/Lz; + kz = linspace(-Kmax,Kmax,Nz+1); + kz = kz(1:end-1); + + psiz = exp(-0.5*(abs(z)/VParams.ell).^VParams.nu);% VParams.nu + normz = sum(abs(psiz).^2)*dz; + psiz = psiz / sqrt(normz); + + psik = fftn(psiz)*normfac/(sqrt(2*pi)); % fftshifted + nk = abs(psik).^2; + nk_forfit = ifftshift(nk); + + f = fit(kz',nk_forfit','gauss1'); % f(x) = a1*exp(-((x-b1)/c1)^2) + + VParams.ell_eff = 1/f.c1; +end \ No newline at end of file diff --git a/Dipolar-Gas-Simulator/+Variational2D/@DipolarGas/DipolarGas.m b/Dipolar-Gas-Simulator/+Variational2D/@DipolarGas/DipolarGas.m new file mode 100644 index 0000000..3224d42 --- /dev/null +++ b/Dipolar-Gas-Simulator/+Variational2D/@DipolarGas/DipolarGas.m @@ -0,0 +1,177 @@ +classdef DipolarGas < handle & matlab.mixin.Copyable + + properties (Access = public) + NumberOfAtoms; + DipolarPolarAngle; + DipolarAzimuthAngle + ScatteringLength; + TrapFrequencies; + NumberOfGridPoints; + Dimensions; + + SimulationMode; + TimeStepSize; + NumberOfTimeSteps; + EnergyTolerance; + ResidualTolerance; + MinimumTimeStepSize; + + Calculator; + + SimulationParameters; + + JobNumber; + RunOnGPU; + DebugMode; + DoSave; + SaveDirectory; + end + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %- Methods + + methods + function this = DipolarGas(varargin) + + p = inputParser; + p.KeepUnmatched = true; + addParameter(p, 'NumberOfAtoms', 1E6,... + @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); + addParameter(p, 'DipolarPolarAngle', pi/2,... + @(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, 'ScatteringLength', 120,... + @(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,2),... + @(x) assert(isnumeric(x) && isvector(x) && all(x > 0))); + addParameter(p, 'Dimensions', 10 * ones(1,2),... + @(x) assert(isnumeric(x) && isvector(x) && all(x > 0))); + addParameter(p, 'TimeStepSize', 5E-4,... + @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); + addParameter(p, 'NumberOfTimeSteps', 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, 'JobNumber', 1,... + @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); + addParameter(p, 'RunOnGPU', false,... + @islogical); + addParameter(p, 'DebugMode', false,... + @islogical); + addParameter(p, 'SaveData', false,... + @islogical); + addParameter(p, 'SaveDirectory', './Data',... + @ischar); + + 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.TimeStepSize = p.Results.TimeStepSize; + this.NumberOfTimeSteps = p.Results.NumberOfTimeSteps; + this.EnergyTolerance = p.Results.EnergyTolerance; + this.ResidualTolerance = p.Results.ResidualTolerance; + this.MinimumTimeStepSize = p.Results.MinimumTimeStepSize; + + 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(); + + this.SimulationParameters = this.setupParameters(); + + end + end + + methods + function set.TimeStepSize(this, val) + assert(val > 1e-06, 'Not time efficient to compute for time steps smaller than 1 microsecond!'); + this.TimeStepSize = val; + end + function ret = get.TimeStepSize(this) + ret = this.TimeStepSize; + end + function set.NumberOfTimeSteps(this, val) + assert(val <= 2E6, 'Not time efficient to compute for time spans longer than 2E6 seconds!'); + this.NumberOfTimeSteps = val; + end + function ret = get.NumberOfTimeSteps(this) + ret = this.NumberOfTimeSteps; + end + function set.NumberOfAtoms(this, val) + assert(val <= 1E6, '!!Not time efficient to compute for atom numbers larger than 1,000,000!!'); + this.NumberOfAtoms = val; + end + function ret = get.NumberOfAtoms(this) + ret = this.NumberOfAtoms; + end + function set.DebugMode(this, val) + this.DebugMode = val; + end + function ret = get.DebugMode(this) + ret = this.DebugMode; + end + function set.DoSave(this, val) + this.DoSave = val; + end + function ret = get.DoSave(this) + ret = this.DoSave; + end + function set.SaveDirectory(this, val) + this.SaveDirectory = val; + end + function ret = get.SaveDirectory(this) + ret = this.SaveDirectory; + end + end % - setters and getters + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %- Methods + + methods(Access = protected) + function cp = copyElement(this) + % Shallow copy object + cp = copyElement@matlab.mixin.Copyable(this); + + % Forces the setter to redefine the function handles to the new copied object + + pl = properties(this); + for k = 1:length(pl) + sc = superclasses(this.(pl{k})); + if any(contains(sc,{'matlab.mixin.Copyable'})) + cp.(pl{k}) = this.(pl{k}).copy(); + end + end + end + end + + methods (Static) + + % Creates an Instance of Class, ensures singleton behaviour (that there + % can only be one Instance of this class + function singleObj = getInstance(varargin) + % Creates an Instance of Class, ensures singleton behaviour + persistent localObj; + if isempty(localObj) || ~isvalid(localObj) + localObj = Simulator.DipolarGas(varargin{:}); + end + singleObj = localObj; + end + end + +end diff --git a/Dipolar-Gas-Simulator/+Variational2D/@DipolarGas/initialize.m b/Dipolar-Gas-Simulator/+Variational2D/@DipolarGas/initialize.m new file mode 100644 index 0000000..5abe950 --- /dev/null +++ b/Dipolar-Gas-Simulator/+Variational2D/@DipolarGas/initialize.m @@ -0,0 +1,26 @@ +function [psi,V,VDk] = initialize(this,Params,VParams,Transf) + + format long + X = Transf.X; Y = Transf.Y; + rcut = min(Transf.Xmax,Transf.Ymax); + + % == Potential == % + V = 0.5*(Params.gx.*X.^2+Params.gy.*Y.^2); + + % == Calculating the DDIs == % + if isfile(strcat(this.SaveDirectory, '/VDk_M.mat')) + VDk = load(sprintf(strcat(this.SaveDirectory, '/VDk_M.mat'))); + VDk = VDk.VDk; + else + VDk = this.Calculator.calculateVDkWithCutoff(rcut, Transf, Params, VParams); + 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) + + % == Setting up the initial wavefunction == % + psi = this.setupWavefunction(Params,Transf); + + if this.RunOnGPU + psi = gpuArray(psi); + end +end \ No newline at end of file diff --git a/Dipolar-Gas-Simulator/+Variational2D/@DipolarGas/propagateWavefunction.m b/Dipolar-Gas-Simulator/+Variational2D/@DipolarGas/propagateWavefunction.m new file mode 100644 index 0000000..01f68b4 --- /dev/null +++ b/Dipolar-Gas-Simulator/+Variational2D/@DipolarGas/propagateWavefunction.m @@ -0,0 +1,117 @@ +function [psi] = propagateWavefunction(this, psi, Params, VParams, Transf, VDk, V, t_idx, Observ) + set(0,'defaulttextInterpreter','latex') + set(groot, 'defaultAxesTickLabelInterpreter','latex'); set(groot, 'defaultLegendInterpreter','latex'); + + dt =-1j*abs(this.TimeStepSize); % Imaginary Time + + KEop = 0.5*(Transf.KX.^2+Transf.KY.^2); + + Observ.res = 1; + AdaptIdx = 0; + + % Change in Energy + E = this.Calculator.calculateTotalEnergy(psi,Params,VParams,Transf,VDk,V); + E = E/Params.N; + Observ.EVec = [Observ.EVec E]; + + % Chemical Potential + muchem = this.Calculator.calculateChemicalPotential(psi,Params,VParams,Transf,VDk,V); + Observ.mucVec = [Observ.mucVec muchem]; + + % Normalized residuals + res = this.Calculator.calculateNormalizedResiduals(psi,Params,VParams,Transf,VDk,V,muchem); + Observ.residual = [Observ.residual res]; + + g_eff = Params.gs*VParams.nu/(2^(1+1/VParams.nu)*VParams.ell*gamma(1/VParams.nu)); + gamma_eff = Params.gammaQF*2^(1/VParams.nu-1.5)*5^(-1/VParams.nu)*VParams.ell*gamma(1+1/VParams.nu)*( VParams.nu/(VParams.ell*gamma(1/VParams.nu)) )^(5/2); + EVar = VParams.nu^2*gamma(2-1/VParams.nu)/(8*VParams.ell^2*gamma(1/VParams.nu)) + 0.5*Params.gz*VParams.ell^2*gamma(3/VParams.nu)/gamma(1/VParams.nu); + + pb = Helper.ProgressBar(); + pb.run('Running simulation in imaginary time: '); + + while t_idx < Params.sim_time_cut_off + + % 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)); + + % Real-space + psi = psi.*exp(-1i*dt*(V + EVar + g_eff*abs(psi).^2 + gamma_eff*abs(psi).^3 + Params.gdd*Phi/(sqrt(2*pi)*VParams.ell_eff) - muchem)); + + % 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; + psi = sqrt(Params.N)*psi/sqrt(Norm); + + muchem = this.Calculator.calculateChemicalPotential(psi,Params,VParams,Transf,VDk,V); + + if mod(t_idx,1000) == 0 + + % Change in Energy + E = this.Calculator.calculateTotalEnergy(psi,Params,VParams,Transf,VDk,V); + E = E/Params.N; + Observ.EVec = [Observ.EVec E]; + + % Chemical potential + Observ.mucVec = [Observ.mucVec muchem]; + + % Normalized residuals + res = this.Calculator.calculateNormalizedResiduals(psi,Params,VParams,Transf,VDk,V,muchem); + Observ.residual = [Observ.residual res]; + 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'); + + %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-4 + if AdaptIdx > 3 && abs(dt) > Params.mindt + dt = dt / 2; + fprintf('Time step changed to '); disp(dt); + AdaptIdx = 0; + elseif AdaptIdx > 3 && abs(dt) < Params.mindt + break + elseif Observ.residual(Observ.res_idx) < Params.rtol + break + else + AdaptIdx = AdaptIdx + 1; + end + else + AdaptIdx = 0; + end + end + if any(isnan(psi(:))) + disp('NaNs encountered!') + break + end + t_idx=t_idx+1; + end + + % Change in Energy + E = this.Calculator.calculateTotalEnergy(psi,Params,VParams,Transf,VDk,V); + E = E/Norm; + Observ.EVec = [Observ.EVec E]; + + % Normalized residuals + res = this.Calculator.calculateNormalizedResiduals(psi,Params,VParams,Transf,VDk,V,muchem); + Observ.residual = [Observ.residual res]; + Observ.res_idx = Observ.res_idx + 1; + + %Chemical potential + Observ.mucVec = [Observ.mucVec muchem]; + + pb.run(' - Job Completed!'); + disp('Saving data...'); + save(sprintf('./Data/Run_%03i/psi_gs.mat',Params.njob),'psi','muchem','Observ','t_idx','Transf','Params','VDk','V'); + disp('Save complete!'); +end \ No newline at end of file diff --git a/Dipolar-Gas-Simulator/+Variational2D/@DipolarGas/run.m b/Dipolar-Gas-Simulator/+Variational2D/@DipolarGas/run.m new file mode 100644 index 0000000..66e9017 --- /dev/null +++ b/Dipolar-Gas-Simulator/+Variational2D/@DipolarGas/run.m @@ -0,0 +1,72 @@ +function [Params, Transf, psi, V, VDk] = run(this) + format long; + + % --- Obtain simulation parameters --- + [Params] = this.setupParameters(); + this.SimulationParameters = Params; + + % --- Set up spatial grids and transforms --- + [Transf] = this.setupSpace(Params); + + % --- Set up for Variational method --- + VarArray = Params.VarArray; + VParams.ell = VarArray(1); + VParams.nu = VarArray(2); + [VParams] = find_nk_fit(VParams); + + % --- Initialize --- + mkdir(sprintf(this.SaveDirectory)) + mkdir(sprintf('./Data/Run_%03i',Params.njob)) + fminconoptions = optimoptions('fmincon','Display','off','StepTolerance',1e-8); + + + [psi,V,VDk,~] = this.initialize(Params,VParams,Transf); + ells(1) = VarArray(1); + nus(1) = VarArray(2); + E_Var = @(x) this.Calculator.calculateVariationalEnergy(psi,Params,x,Transf,VDk,V)/Params.N; + E_vs_iter(1) = E_Var([ells(1) nus(1)]); + + t_idx = 1; % Start at t = 0; + + for nn = 1:Params.SelfConIter + Observ.EVec = []; Observ.NormVec = []; Observ.PCVec = []; Observ.tVecPlot = []; Observ.mucVec = []; + Observ.res_idx = 1; + + VParams.ell = VarArray(1); + VParams.nu = VarArray(2); + [VParams] = find_nk_fit(VParams); + + [psi,V,VDk,~] = this.initialize(Params,VParams,Transf); + + % --- Adding some noise --- + Norm = trapz(abs(psi(:)).^2)*Transf.dx*Transf.dy; % normalisation + psi = sqrt(Params.N)*psi/sqrt(Norm); + r = normrnd(0,1,size(psi)); + theta = rand(size(psi)); + noise = r.*exp(2*pi*1i*theta); + psi = psi + 4*noise; + + % --- Run --- + [psi,Observ] = this.propagateWavefunction(psi,Params,Transf,VDk,V,t_idx,Observ); + psi = gather(psi); + + % --- Constrained minimization --- + E_Var = @(x) this.Calculator.calculateVariationalEnergy(psi,Params,x,Transf,V)/Params.N; + VarArray = fmincon(E_Var,VarArray,[],[],[],[],[Params.ell_lower, Params.nu_lower],[Params.ell_upper, Params.nu_upper],[],fminconoptions); + + % --- Convergence check --- + ells = [ells VarArray(1)]; + nus = [nus VarArray(2)]; + relelldiff = abs(ells(nn+1)-ells(nn))/ells(nn); + relnudiff = abs(nus(nn+1)-nus(nn))/nus(nn); + E_vs_iter = [E_vs_iter E_Var(VarArray)]; + + save(sprintf('./Data/Run_%03i/psi_gs_%i.mat',runIdx,nn),'psi','Observ','Transf','Params','VDk','V','VarArray'); + + if relelldiff < Params.ellcutoff && relnudiff < Params.nucutoff + break + end + end + save(sprintf('./Data/Run_%03i/psi_gs.mat',runIdx),'psi','Observ','Transf','Params','VDk','V','VarArray'); + delete(sprintf('./Data/Run_%03i/psi_gs_*.mat',runIdx)) +end diff --git a/Dipolar-Gas-Simulator/+Variational2D/@DipolarGas/setupParameters.m b/Dipolar-Gas-Simulator/+Variational2D/@DipolarGas/setupParameters.m new file mode 100644 index 0000000..4fedd24 --- /dev/null +++ b/Dipolar-Gas-Simulator/+Variational2D/@DipolarGas/setupParameters.m @@ -0,0 +1,109 @@ +function [Params] = setupParameters(this) + CONSTANTS = Helper.PhysicsConstants; + hbar = CONSTANTS.PlanckConstantReduced; % [J.s] + kbol = CONSTANTS.BoltzmannConstant; % [J/K] + mu0 = CONSTANTS.VacuumPermeability; % [N/A^2] + muB = CONSTANTS.BohrMagneton; % [J/T] + a0 = CONSTANTS.BohrRadius; % [m] + m0 = CONSTANTS.AtomicMassUnit; % [kg] + w0 = 2*pi*100; % 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 + Params.Nx = this.NumberOfGridPoints(1); + Params.Ny = this.NumberOfGridPoints(2); + + % Dimensions (in units of l0) + Params.Lx = this.Dimensions(1); + Params.Ly = this.Dimensions(2); + + % Mass, length scale + Params.m = CONSTANTS.Dy164Mass; + l0 = sqrt(hbar/(Params.m*w0)); % Defining a harmonic oscillator length + + % Atom numbers + Params.N = this.NumberOfAtoms; + + % Dipole angle + Params.theta = this.DipolarPolarAngle; % pi/2 dipoles along x, theta=0 dipoles along z + Params.phi = this.DipolarAzimuthAngle; + + % Dipole lengths (units of muB) + Params.mu = CONSTANTS.DyMagneticMoment; + + % Scattering lengths + Params.as = this.ScatteringLength*a0; + + % Trapping frequencies + Params.wx = 2*pi*this.TrapFrequencies(1); + Params.wy = 2*pi*this.TrapFrequencies(2); + Params.wz = 2*pi*this.TrapFrequencies(3); + + % 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.njob = this.JobNumber; + + + % ================ Variational method parameters ================ % + % FMinCon Settings + Params.SelfConIter = 20; % Max number of iterations to perform self-consistent calculation + Params.VarArray = [4 2.5]; % initial [ell nu], expand for other params + % ell is the "width" and nu is the exponent. psi~ e^(z/ell)^nu + + + % Window of optimization + Params.ell_lower = 0.2; + Params.ell_upper = 12; + + Params.nu_lower = 1; + Params.nu_upper = 4; + + % Relative cutoffs + Params.ellcutoff = 1e-2; + Params.nucutoff = 1e-2; + + % ================ Parameters defined by those above ================ % + + % Contact interaction strength (units of l0/m) + Params.gs = 4*pi*Params.as/l0; + + % Dipole lengths + 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 + + % Trap gamma + Params.gx = (Params.wx/w0)^2; + Params.gy = (Params.wy/w0)^2; + Params.gz = (Params.wz/w0)^2; + + % == Calculate LHY correction to account for quantum fluctuations == % + + eps_dd = Params.add/Params.as; + if eps_dd == 0 + Q5 = 1; + elseif eps_dd == 1 + Q5 = 3*sqrt(3)/2; + else + yeps = (1-eps_dd)/(3*eps_dd); + Q5 = (3*eps_dd)^(5/2)*( (8+26*yeps+33*yeps^2)*sqrt(1+yeps) + 15*yeps^3*log((1+sqrt(1+yeps))/sqrt(yeps)) )/48; + Q5 = real(Q5); + end + + Params.gammaQF = 128/3*sqrt(pi*(Params.as/l0)^5)*Q5; + + % Loading the rest into Params + Params.hbar = hbar; + Params.kbol = kbol; + Params.mu0 = mu0; + Params.muB = muB; + Params.a0 = a0; + Params.w0 = w0; + Params.l0 = l0; +end \ No newline at end of file diff --git a/Dipolar-Gas-Simulator/+Variational2D/@DipolarGas/setupSpace.m b/Dipolar-Gas-Simulator/+Variational2D/@DipolarGas/setupSpace.m new file mode 100644 index 0000000..6c9f632 --- /dev/null +++ b/Dipolar-Gas-Simulator/+Variational2D/@DipolarGas/setupSpace.m @@ -0,0 +1,33 @@ +function [Transf] = setupSpace(~,Params) + Transf.Xmax = 0.5*Params.Lx; + Transf.Ymax = 0.5*Params.Ly; + + 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); + + 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); + + [Transf.X,Transf.Y] = ndgrid(x,y); + [Transf.KX,Transf.KY] = ndgrid(kx,ky); + Transf.x = x; + Transf.y = y; + Transf.kx = kx; + Transf.ky = ky; + Transf.dx = x(2)-x(1); + Transf.dy = y(2)-y(1); + Transf.dkx = dkx; + Transf.dky = dky; +end \ No newline at end of file diff --git a/Dipolar-Gas-Simulator/+Variational2D/@DipolarGas/setupWavefunction.m b/Dipolar-Gas-Simulator/+Variational2D/@DipolarGas/setupWavefunction.m new file mode 100644 index 0000000..3e488a9 --- /dev/null +++ b/Dipolar-Gas-Simulator/+Variational2D/@DipolarGas/setupWavefunction.m @@ -0,0 +1,21 @@ +function [psi] = setupWavefunction(~,Params,Transf) + + X = Transf.X; + Y = Transf.Y; + + ellx = sqrt(Params.hbar/(Params.m*Params.wx))/Params.l0; + elly = sqrt(Params.hbar/(Params.m*Params.wy))/Params.l0; + + Rx = 6*sqrt(2)*ellx; + Ry = 6*sqrt(2)*elly; + X0 = 0.0*Transf.Xmax; + Y0 = 0.0*Transf.Ymax; + + psi = exp(-(X-X0).^2/Rx^2-(Y-Y0).^2/Ry^2); + cur_norm = trapz(abs(psi(:)).^2)*Transf.dx*Transf.dy; + psi = psi/sqrt(cur_norm); + + Norm = trapz(abs(psi(:)).^2)*Transf.dx*Transf.dy; + psi = sqrt(Params.N)*psi/sqrt(Norm); + +end \ No newline at end of file