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

Adapted elements from a different working version of the solver from Wyatt - minor changes such as replacing trapz with sum and a few others to make the versions similar.

Browse Source
This commit is contained in:
Karthik 2024-11-30 08:26:12 +01:00
parent 5049970246
commit ac47adcce3
12 changed files with 57 additions and 75 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -69,7 +69,7 @@ OptionsStruct.DipolarPolarAngle = 0;
OptionsStruct.DipolarAzimuthAngle = 0; OptionsStruct.DipolarAzimuthAngle = 0;
OptionsStruct.ScatteringLength = 98.0676; OptionsStruct.ScatteringLength = 98.0676;
OptionsStruct.TrapFrequencies = [10, 10, 72.4]; OptionsStruct.TrapFrequencies = [0, 0, 72.4];
OptionsStruct.TrapPotentialType = 'None'; OptionsStruct.TrapPotentialType = 'None';
OptionsStruct.NumberOfGridPoints = [128, 128]; OptionsStruct.NumberOfGridPoints = [128, 128];
@ -78,14 +78,14 @@ OptionsStruct.TimeStepSize = 500E-6; % in s
OptionsStruct.MinimumTimeStepSize = 1E-5; % in s OptionsStruct.MinimumTimeStepSize = 1E-5; % in s
OptionsStruct.TimeCutOff = 1E6; % in s OptionsStruct.TimeCutOff = 1E6; % in s
OptionsStruct.EnergyTolerance = 5E-10; OptionsStruct.EnergyTolerance = 5E-10;
OptionsStruct.ResidualTolerance = 1E-04; OptionsStruct.ResidualTolerance = 1E-05;
OptionsStruct.NoiseScaleFactor = 4; OptionsStruct.NoiseScaleFactor = 4;
OptionsStruct.MaxIterations = 1; OptionsStruct.MaxIterations = 10;
OptionsStruct.VariationalWidth = 5.0; OptionsStruct.VariationalWidth = 5.0;
OptionsStruct.WidthLowerBound = 0.2; OptionsStruct.WidthLowerBound = 1;
OptionsStruct.WidthUpperBound = 12; OptionsStruct.WidthUpperBound = 12;
OptionsStruct.WidthCutoff = 1e-2; OptionsStruct.WidthCutoff = 5e-3;
OptionsStruct.PlotLive = true; OptionsStruct.PlotLive = true;
OptionsStruct.JobNumber = 1; OptionsStruct.JobNumber = 1;
@ -101,7 +101,7 @@ solver.Potential = pot.trap();
%-% Run Solver %-% %-% Run Solver %-%
[Params, Transf, psi, V, VDk] = solver.run(); [Params, Transf, psi, V, VDk] = solver.run();
%%
% Solve BdG equations % Solve BdG equations
% Load data % Load data
Data = load(sprintf(horzcat(solver.SaveDirectory, '/Run_%03i/psi_gs.mat'),solver.JobNumber),'psi','Transf','Params','VParams','V'); Data = load(sprintf(horzcat(solver.SaveDirectory, '/Run_%03i/psi_gs.mat'),solver.JobNumber),'psi','Transf','Params','VParams','V');

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

@ -15,7 +15,7 @@ function muchem = calculateChemicalPotential(~,psi,Params,VParams,Transf,VDk,V)
% Kinetic energy % Kinetic energy
Ekin = KEop.*abs(fftn(psi)*normfac).^2; Ekin = KEop.*abs(fftn(psi)*normfac).^2;
Ekin = trapz(Ekin(:))*Transf.dkx*Transf.dky/(2*pi)^2; Ekin = sum(Ekin(:))*Transf.dkx*Transf.dky/(2*pi)^2;
% Potential energy % Potential energy
Epot = V.*abs(psi).^2; Epot = V.*abs(psi).^2;
@ -27,6 +27,7 @@ function muchem = calculateChemicalPotential(~,psi,Params,VParams,Transf,VDk,V)
Eqf = gamma_eff*abs(psi).^5; Eqf = gamma_eff*abs(psi).^5;
% Total energy % Total energy
muchem = Ekin + Ez*Params.N + trapz(Epot(:) + Eint(:) + Eddi(:) + Eqf(:))*Transf.dx*Transf.dy; % muchem = Ekin + Ez*Params.N + sum(Epot(:) + Eint(:) + Eddi(:) + Eqf(:))*Transf.dx*Transf.dy;
muchem = muchem / Params.N; % Only use if psi is normalized to N muchem = muchem / Params.N; % Only use if psi is normalized to N
end end

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

@ -1,29 +0,0 @@
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

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

@ -24,6 +24,6 @@ function res = calculateNormalizedResiduals(~,psi,Params,VParams,Transf,VDk,V,m
Eqf = gamma_eff*abs(psi).^3.*psi; Eqf = gamma_eff*abs(psi).^3.*psi;
% Total energy % Total energy
res = trapz(abs(Ekin(:) + Ez*psi(:) + Epot(:) + Eint(:) + Eddi(:) + Eqf(:) - muchem*psi(:))*Transf.dx*Transf.dy)/trapz(abs(muchem*psi(:))*Transf.dx*Transf.dy); res = sum(abs(Ekin(:) + Ez*psi(:) + Epot(:) + Eint(:) + Eddi(:) + Eqf(:) - muchem*psi(:))*Transf.dx*Transf.dy)/sum(abs(muchem*psi(:))*Transf.dx*Transf.dy);
end end

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

@ -12,23 +12,24 @@ function E = calculateTotalEnergy(~,psi,Params,VParams,Transf,VDk,V)
frho = fftn(abs(psi).^2); frho = fftn(abs(psi).^2);
Phi = real(ifftn(frho.*VDk)); Phi = real(ifftn(frho.*VDk));
Eddi = 0.5*Params.gdd*Phi.*abs(psi).^2/(sqrt(2*pi)*VParams.ell); Eddi = 0.5*Params.gdd*Phi.*abs(psi).^2/(sqrt(2*pi)*VParams.ell);
Eddi = trapz(Eddi(:))*Transf.dx*Transf.dy; Eddi = sum(Eddi(:))*Transf.dx*Transf.dy;
% Kinetic energy % Kinetic energy
Ekin = KEop.*abs(fftn(psi)*normfac).^2; Ekin = KEop.*abs(fftn(psi)*normfac).^2;
Ekin = trapz(Ekin(:))*Transf.dkx*Transf.dky/(2*pi)^2; Ekin = sum(Ekin(:))*Transf.dkx*Transf.dky/(2*pi)^2;
% Potential energy % Potential energy
Epot = V.*abs(psi).^2; Epot = V.*abs(psi).^2;
Epot = trapz(Epot(:))*Transf.dx*Transf.dy; Epot = sum(Epot(:))*Transf.dx*Transf.dy;
% Contact interactions % Contact interactions
Eint = 0.5*g_eff*abs(psi).^4; Eint = 0.5*g_eff*abs(psi).^4;
Eint = trapz(Eint(:))*Transf.dx*Transf.dy; Eint = sum(Eint(:))*Transf.dx*Transf.dy;
% Quantum fluctuations % Quantum fluctuations
Eqf = 0.4*gamma_eff*abs(psi).^5; Eqf = 0.4*gamma_eff*abs(psi).^5;
Eqf = trapz(Eqf(:))*Transf.dx*Transf.dy; Eqf = sum(Eqf(:))*Transf.dx*Transf.dy;
E = Ez*Params.N + Ekin + Epot + Eint + Eddi + Eqf; E = Ez*Params.N + Ekin + Epot + Eint + Eddi + Eqf;
end end

10
Dipolar-Gas-Simulator/+VariationalSolver2D/@Calculator/calculateVariationalEnergy.m
View File

@ -12,24 +12,24 @@ function E = calculateVariationalEnergy(this, psi, Params, ell, Transf, V)
VDk = this.calculateVDkWithCutoff(Transf, Params, ell); % VDk depends on the variational parameters. THIS HAS TO BE RECALCULATED HERE FOR THE CONSTRAINED OPTIMIZATION TO WORK! VDk = this.calculateVDkWithCutoff(Transf, Params, ell); % VDk depends on the variational parameters. THIS HAS TO BE RECALCULATED HERE FOR THE CONSTRAINED OPTIMIZATION TO WORK!
frho = fftn(abs(psi).^2); frho = fftn(abs(psi).^2);
Phi = real(ifftn(frho.*VDk)); Phi = real(ifftn(frho.*VDk));
Eddi = 0.5*Params.gdd*Phi.*abs(psi).^2/(sqrt(2*pi)*ell);% Eddi = 0.5*Params.gdd*Phi.*abs(psi).^2/(sqrt(2*pi)*ell);
Eddi = sum(Eddi(:))*Transf.dx*Transf.dy; Eddi = sum(Eddi(:))*Transf.dx*Transf.dy;
% Kinetic energy % Kinetic energy
Ekin = KEop.*abs(fftn(psi)*normfac).^2; Ekin = KEop.*abs(fftn(psi)*normfac).^2;
Ekin = trapz(Ekin(:))*Transf.dkx*Transf.dky/(2*pi)^2; Ekin = sum(Ekin(:))*Transf.dkx*Transf.dky/(2*pi)^2;
% Potential energy % Potential energy
Epot = V.*abs(psi).^2; Epot = V.*abs(psi).^2;
Epot = trapz(Epot(:))*Transf.dx*Transf.dy; Epot = sum(Epot(:))*Transf.dx*Transf.dy;
% Contact interactions % Contact interactions
Eint = 0.5*g_eff*abs(psi).^4; Eint = 0.5*g_eff*abs(psi).^4;
Eint = trapz(Eint(:))*Transf.dx*Transf.dy; Eint = sum(Eint(:))*Transf.dx*Transf.dy;
% Quantum fluctuations % Quantum fluctuations
Eqf = 0.4*gamma_eff*abs(psi).^5; Eqf = 0.4*gamma_eff*abs(psi).^5;
Eqf = trapz(Eqf(:))*Transf.dx*Transf.dy; Eqf = sum(Eqf(:))*Transf.dx*Transf.dy;
% Total % Total
E = Ez*Params.N + Ekin + Epot + Eint + Eddi + Eqf; E = Ez*Params.N + Ekin + Epot + Eint + Eddi + Eqf;

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

@ -53,7 +53,7 @@ classdef DipolarGas < handle & matlab.mixin.Copyable
addParameter(p, 'ScatteringLength', 120,... addParameter(p, 'ScatteringLength', 120,...
@(x) assert(isnumeric(x) && isscalar(x) && (x > -150) && (x < 150))); @(x) assert(isnumeric(x) && isscalar(x) && (x > -150) && (x < 150)));
addParameter(p, 'TrapFrequencies', 100 * ones(1,3),... addParameter(p, 'TrapFrequencies', 100 * ones(1,3),...
@(x) assert(isnumeric(x) && isvector(x) && all(x > 0))); @(x) assert(isnumeric(x) && isvector(x) && all(x >= 0)));
addParameter(p, 'NumberOfGridPoints', 128 * ones(1,2),... addParameter(p, 'NumberOfGridPoints', 128 * ones(1,2),...
@(x) assert(isnumeric(x) && isvector(x) && all(x > 0))); @(x) assert(isnumeric(x) && isvector(x) && all(x > 0)));
addParameter(p, 'Dimensions', 10 * ones(1,2),... addParameter(p, 'Dimensions', 10 * ones(1,2),...

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

@ -51,7 +51,7 @@ function [psi, Observ] = propagateWavefunction(this, psi, Params, VParams, Trans
psi = ifftn(psi); psi = ifftn(psi);
% Renorm % Renorm
Norm = trapz(abs(psi(:)).^2)*Transf.dx*Transf.dy; Norm = sum(abs(psi(:)).^2)*Transf.dx*Transf.dy;
psi = sqrt(Params.N)*psi/sqrt(Norm); psi = sqrt(Params.N)*psi/sqrt(Norm);
muchem = this.Calculator.calculateChemicalPotential(psi,Params,VParams,Transf,VDk,V); muchem = this.Calculator.calculateChemicalPotential(psi,Params,VParams,Transf,VDk,V);

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

@ -23,29 +23,38 @@ function [Params, Transf, psi, V, VDk] = run(this)
% --- Initialize --- % --- Initialize ---
mkdir(sprintf(this.SaveDirectory)) mkdir(sprintf(this.SaveDirectory))
mkdir(sprintf(strcat(this.SaveDirectory, '/Run_%03i'),Params.njob)) mkdir(sprintf(strcat(this.SaveDirectory, '/Run_%03i'),Params.njob))
fminconoptions = optimoptions('fmincon','Display','off','StepTolerance',1e-8); fminconoptions = optimoptions('fmincon','Display','off','StepTolerance',1e-10,'Algorithm','sqp','MaxIterations',2000,'MaxFunctionEvaluations',4000);
[psi,V,VDk] = this.initialize(Params,VParams,Transf); [psi,V,VDk] = this.initialize(Params,VParams,Transf);
ells(1) = VParams.ell; ells(1) = VParams.ell;
E_Var = @(x) this.Calculator.calculateVariationalEnergy(psi, Params, x, Transf, V)/Params.N; E_Var = @(x) this.Calculator.calculateVariationalEnergy(psi, Params, x, Transf, V)/Params.N;
E_vs_iter(1) = E_Var(ells(1)); E_vs_iter(1) = E_Var(ells(1));
t_idx = 1; firstpoint = 1;
for nn = 1:Params.SelfConIter for nn = 1:Params.SelfConIter
Observ.EVec = []; Observ.NormVec = []; Observ.PCVec = []; Observ.tVecPlot = []; Observ.mucVec = []; Observ.residual = []; Observ.res_idx = 1; Observ.EVec = []; Observ.NormVec = []; Observ.PCVec = []; Observ.tVecPlot = []; Observ.mucVec = []; Observ.residual = []; Observ.res_idx = 1;
t_idx = 1; % Start at t = 0; if firstpoint == 1
[psi,V,VDk] = this.initialize(Params,VParams,Transf); % Recalculate Psi,VDk with new value for the variational parameter
firstpoint = 0;
else
[~,V,VDk] = this.initialize(Params,VParams,Transf); % Recalculate Psi,VDk with new value for the variational parameter [~,V,VDk] = this.initialize(Params,VParams,Transf); % Recalculate Psi,VDk with new value for the variational parameter
end
% --- Adding some noise --- % --- Adding some noise ---
% Noise added in every iteration to ensure it is not stuck in some local minimum % Noise added in every iteration to ensure it is not stuck in some local minimum
Norm = trapz(abs(psi(:)).^2)*Transf.dx*Transf.dy; % normalisation Norm = sum(abs(psi(:)).^2)*Transf.dx*Transf.dy; % normalisation
psi = sqrt(Params.N)*psi/sqrt(Norm); psi = psi/sqrt(Norm);
r = normrnd(0,1,size(psi)); r = normrnd(0,1,size(psi));
theta = rand(size(psi)); theta = rand(size(psi));
noise = r.*exp(2*pi*1i*theta); noise = r.*exp(2*pi*1i*theta);
psi = psi + Params.nsf*noise; psi = psi + Params.nsf*noise;
% Renormalize wavefunction
Norm = sum(abs(psi(:)).^2)*Transf.dx*Transf.dy; % normalisation
psi = sqrt(Params.N)*psi/sqrt(Norm);
% --- Run --- % --- Run ---
[psi, Observ] = this.propagateWavefunction(psi, Params, VParams, Transf, VDk, V, t_idx, Observ, nn); [psi, Observ] = this.propagateWavefunction(psi, Params, VParams, Transf, VDk, V, t_idx, Observ, nn);
psi = gather(psi); psi = gather(psi);

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

@ -17,9 +17,8 @@ function [Params] = setupParameters(this)
Params.Lx = this.Dimensions(1); Params.Lx = this.Dimensions(1);
Params.Ly = this.Dimensions(2); Params.Ly = this.Dimensions(2);
% Mass, length scale % Mass
Params.m = CONSTANTS.Dy164Mass; Params.m = CONSTANTS.Dy164Mass;
l0 = sqrt(hbar/(Params.m*w0)); % Defining a harmonic oscillator length
% Atom numbers % Atom numbers
Params.N = this.NumberOfAtoms; Params.N = this.NumberOfAtoms;
@ -65,12 +64,16 @@ function [Params] = setupParameters(this)
Params.evarcutoff = this.VariationalEnergyTolerance; Params.evarcutoff = this.VariationalEnergyTolerance;
% ================ Parameters defined by those above ================ % % ================ 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) % Contact interaction strength (units of l0/m)
Params.gs = 4*pi*Params.as/l0; Params.gs = 4*pi*Params.as/l0;
% Dipole lengths % Dipole lengths
Params.add = mu0*Params.mu^2*Params.m/(12*pi*hbar^2); Params.add = mu0*Params.mu^2*Params.m/(12*pi*hbar^2);
Params.ppum2 = Params.N/(Params.Lx*Params.Ly*(l0*1e6)^2);% Particles per squared micron
Params.ppadd2 = Params.ppum2*(Params.add*1e6)^2; % Particles per squared add
% DDI strength % DDI strength
Params.gdd = 4*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

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

@ -8,16 +8,13 @@ function [psi] = setupWavefunction(~,Params,Transf)
ellx = sqrt(Params.hbar/(Params.m*Params.wx))/Params.l0; ellx = sqrt(Params.hbar/(Params.m*Params.wx))/Params.l0;
elly = sqrt(Params.hbar/(Params.m*Params.wy))/Params.l0; elly = sqrt(Params.hbar/(Params.m*Params.wy))/Params.l0;
Rx = 6*sqrt(2)*ellx; Rx = 8*sqrt(2)*ellx;
Ry = 6*sqrt(2)*elly; Ry = 8*sqrt(2)*elly;
X0 = 0.0*Transf.Xmax; X0 = 0.0*Transf.Xmax;
Y0 = 0.0*Transf.Ymax; Y0 = 0.0*Transf.Ymax;
psi = exp(-(X-X0).^2/Rx^2-(Y-Y0).^2/Ry^2); psi = exp(-(X-X0).^2/Rx^2-(Y-Y0).^2/Ry^2);
cur_norm = trapz(abs(psi(:)).^2)*Transf.dx*Transf.dy; Norm = sum(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); psi = sqrt(Params.N)*psi/sqrt(Norm);
end end

4
Dipolar-Gas-Simulator/+VariationalSolver2D/@Potentials/Potentials.m
View File

@ -26,7 +26,7 @@ classdef Potentials < handle & matlab.mixin.Copyable
addParameter(p, 'TrapPotentialType', this.PotentialDefaults.TrapPotentialType,... addParameter(p, 'TrapPotentialType', this.PotentialDefaults.TrapPotentialType,...
@(x) assert(any(strcmpi(x,{'None','Harmonic','SquareBox','RoundBox'})))); @(x) assert(any(strcmpi(x,{'None','Harmonic','SquareBox','RoundBox'}))));
addParameter(p, 'TrapFrequencies', this.PotentialDefaults.TrapFrequencies,... addParameter(p, 'TrapFrequencies', this.PotentialDefaults.TrapFrequencies,...
@(x) assert(isnumeric(x) && isvector(x) && all(x > 0))); @(x) assert(isnumeric(x) && isvector(x) && all(x >= 0)));
addParameter(p, 'NumberOfGridPoints', 128 * ones(1,3),... addParameter(p, 'NumberOfGridPoints', 128 * ones(1,3),...
@(x) assert(isnumeric(x) && isvector(x) && all(x > 0))); @(x) assert(isnumeric(x) && isvector(x) && all(x > 0)));
addParameter(p, 'Dimensions', 10 * ones(1,3),... addParameter(p, 'Dimensions', 10 * ones(1,3),...
@ -68,7 +68,7 @@ classdef Potentials < handle & matlab.mixin.Copyable
methods methods
function set.TrapFrequencies(this, val) function set.TrapFrequencies(this, val)
assert(isnumeric(val) && isvector(val) && all(val > 0), 'Incorrectly defined trap frequencies!'); assert(isnumeric(val) && isvector(val) && all(val >= 0), 'Incorrectly defined trap frequencies!');
this.TrapFrequencies = val; this.TrapFrequencies = val;
end end
function ret = get.TrapFrequencies(this) function ret = get.TrapFrequencies(this)