karthik/Calculations
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Modifications to have only width of Gaussian ansatz as the variational parameter and other changes to keep consistent with the source code as written originally by Wyatt.

Browse Source
This commit is contained in:
Karthik 2024-11-15 14:33:46 +01:00
parent 968467a186
commit b3ab383300
12 changed files with 63 additions and 72 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
Dipolar-Gas-Simulator/+BdGSolver2D/solveBogoliubovdeGennesIn2D.m
View File

@ -1,16 +1,15 @@
function [evals, modes] = solveBogoliubovdeGennesIn2D(psi, Params, VDk, VParams, Transf, muchem)
wz_tilde = Params.wz / Params.w0;
gs = Params.gs;
gdd = Params.gdd;
gammaQF = Params.gammaQF;
KEop = 0.5*(Transf.KX.^2+Transf.KY.^2);
Ez = (0.25*VParams.sigma^2) + (0.25*wz_tilde^2*VParams.sigma^2);
Ez = (0.25*VParams.ell^2) + (0.25*Params.gz*VParams.ell^2);
muchem_tilde = muchem - Ez;
g_pf_2D = 1/(sqrt(2*pi)*VParams.sigma);
gQF_pf_2D = sqrt(2/5)/(pi^(3/4)*VParams.sigma^(3/2));
g_pf_2D = 1/(sqrt(2*pi)*VParams.ell);
gQF_pf_2D = sqrt(2/5)/(pi^(3/4)*VParams.ell^(3/2));
% eigs only works with column vectors
psi = psi.';

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

@ -59,7 +59,7 @@ OptionsStruct.ScatteringLength = 102.515;
OptionsStruct.TrapFrequencies = [10, 10, 72.4];
OptionsStruct.NumberOfGridPoints = [1024, 1024];
OptionsStruct.NumberOfGridPoints = [128, 128];
OptionsStruct.Dimensions = [100, 100];
OptionsStruct.TimeStepSize = 1E-3; % in s
OptionsStruct.TimeCutOff = 2E6; % in s
@ -73,7 +73,7 @@ OptionsStruct.SaveDirectory = './Data';
options = Helper.convertstruct2cell(OptionsStruct);
clear OptionsStruct
solver = Variational2D.DipolarGas(options{:});
solver = VariationalSolver2D.DipolarGas(options{:});
%-% Run Solver %-%
[Params, Transf, psi, V, VDk] = solver.run();

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

@ -9,7 +9,7 @@ OptionsStruct.ScatteringLength = 102.515;
OptionsStruct.TrapFrequencies = [10, 10, 72.4];
OptionsStruct.NumberOfGridPoints = [1024, 1024];
OptionsStruct.NumberOfGridPoints = [128, 128];
OptionsStruct.Dimensions = [100, 100];
OptionsStruct.TimeStepSize = 1E-3; % in s
OptionsStruct.TimeCutOff = 2E6; % in s
@ -23,7 +23,7 @@ OptionsStruct.SaveDirectory = './Data';
options = Helper.convertstruct2cell(OptionsStruct);
clear OptionsStruct
solver = Variational2D.DipolarGas(options{:});
solver = VariationalSolver2D.DipolarGas(options{:});
%-% Run Solver %-%
[Params, Transf, psi, V, VDk] = solver.run();

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

@ -1,8 +1,8 @@
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);
g_eff = Params.gs * (1/(sqrt(2*pi)*VParams.ell));
gamma_eff = Params.gammaQF * (sqrt(2/5)/(pi^(3/4)*VParams.ell^(3/2)));
Ez = (0.25*VParams.ell^2) + (0.25*Params.gz*VParams.ell^2);
% Parameters
normfac = Params.Lx*Params.Ly/numel(psi);
@ -11,7 +11,7 @@ function muchem = calculateChemicalPotential(~,psi,Params,VParams,Transf,VDk,V)
% DDIs
frho = fftn(abs(psi).^2);
Phi = real(ifftn(frho.*VDk));
Eddi = (Params.gdd*Phi.*abs(psi).^2)/(sqrt(2*pi)*VParams.ell_eff);
Eddi = (Params.gdd*Phi.*abs(psi).^2)/(sqrt(2*pi)*VParams.ell);
%Kinetic energy
Ekin = KEop.*abs(fftn(psi)*normfac).^2;
@ -27,6 +27,6 @@ function muchem = calculateChemicalPotential(~,psi,Params,VParams,Transf,VDk,V)
Eqf = gamma_eff*abs(psi).^5;
%Total energy
muchem = Ekin + EVar*Params.N + trapz(Epot(:) + Eint(:) + Eddi(:) + Eqf(:))*Transf.dx*Transf.dy; %
muchem = Ekin + Ez*Params.N + trapz(Epot(:) + Eint(:) + Eddi(:) + Eqf(:))*Transf.dx*Transf.dy; %
muchem = muchem / Params.N; %Only use if psi is normalized to N
end

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

@ -1,15 +1,15 @@
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);
g_eff = Params.gs * (1/(sqrt(2*pi)*VParams.ell));
gamma_eff = Params.gammaQF * (sqrt(2/5)/(pi^(3/4)*VParams.ell^(3/2)));
Ez = (0.25*VParams.ell^2) + (0.25*Params.gz*VParams.ell^2);
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);
Eddi = Params.gdd*Phi.*psi/(sqrt(2*pi)*VParams.ell);
% Kinetic energy
Ekin = ifftn(KEop.*fftn(psi));
@ -24,6 +24,6 @@ function res = calculateNormalizedResiduals(~,psi,Params,VParams,Transf,VDk,V,m
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);
res = trapz(abs(Ekin(:) + Ez*psi(:) + Epot(:) + Eint(:) + Eddi(:) + Eqf(:) - muchem*psi(:))*Transf.dx*Transf.dy)/trapz(abs(muchem*psi(:))*Transf.dx*Transf.dy);
end

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

@ -1,8 +1,8 @@
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);
g_eff = Params.gs * (1/(sqrt(2*pi)*VParams.ell));
gamma_eff = Params.gammaQF * (sqrt(2/5)/(pi^(3/4)*VParams.ell^(3/2)));
Ez = (0.25*VParams.ell^2) + (0.25*Params.gz*VParams.ell^2);
% Parameters
KEop = 0.5*(Transf.KX.^2+Transf.KY.^2);
@ -11,7 +11,7 @@ function E = calculateTotalEnergy(~,psi,Params,VParams,Transf,VDk,V)
% 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 = 0.5*Params.gdd*Phi.*abs(psi).^2/(sqrt(2*pi)*VParams.ell);%
Eddi = trapz(Eddi(:))*Transf.dx*Transf.dy;
% Kinetic energy
@ -30,5 +30,5 @@ function E = calculateTotalEnergy(~,psi,Params,VParams,Transf,VDk,V)
Eqf = 0.4*gamma_eff*abs(psi).^5;
Eqf = trapz(Eqf(:))*Transf.dx*Transf.dy;
E = EVar*Params.N + Ekin + Epot + Eint + Eddi + Eqf;
E = Ez*Params.N + Ekin + Epot + Eint + Eddi + Eqf;
end

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

@ -2,8 +2,8 @@ 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);
QX = Transf.KX*VParams.ell/sqrt(2);
QY = Transf.KY*VParams.ell/sqrt(2);
Qsq = QX.^2 + QY.^2;
absQ = sqrt(Qsq);

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

@ -1,21 +1,17 @@
function E = calculateVariationalEnergy(this, psi, Params, VarArray, Transf, VDk, V)
function E = calculateVariationalEnergy(~, psi, Params, ell, 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);
g_eff = Params.gs * (1/(sqrt(2*pi)*ell));
gamma_eff = Params.gammaQF * (sqrt(2/5)/(pi^(3/4)*ell^(3/2)));
Ez = (0.25*ell^2) + (0.25*Params.gz*ell^2);
% Parameters
KEop = 0.5*(Transf.KX.^2+Transf.KY.^2);
normfac = Params.Lx*Params.Ly/numel(psi);
% DDIs
[VParams] = this.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 = 0.5*Params.gdd*Phi.*abs(psi).^2/(sqrt(2*pi)*ell);%
Eddi = sum(Eddi(:))*Transf.dx*Transf.dy;
% Kinetic energy
@ -35,5 +31,6 @@ function E = calculateVariationalEnergy(this, psi, Params, VarArray, Transf, VDk
Eqf = trapz(Eqf(:))*Transf.dx*Transf.dy;
% Total
E = EVar*Params.N + Ekin + Epot + Eint + Eddi + Eqf;
E = Ez*Params.N + Ekin + Epot + Eint + Eddi + Eqf;
end

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

@ -90,7 +90,7 @@ classdef DipolarGas < handle & matlab.mixin.Copyable
this.DoSave = p.Results.SaveData;
this.SaveDirectory = p.Results.SaveDirectory;
this.Calculator = Variational2D.Calculator();
this.Calculator = VariationalSolver2D.Calculator();
this.SimulationParameters = this.setupParameters();

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

@ -22,9 +22,9 @@ function [psi] = propagateWavefunction(this, psi, Params, VParams, Transf, VDk,
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);
g_eff = Params.gs * (1/(sqrt(2*pi)*VParams.ell));
gamma_eff = Params.gammaQF * (sqrt(2/5)/(pi^(3/4)*VParams.ell^(3/2)));
Ez = (0.25*VParams.ell^2) + (0.25*Params.gz*VParams.ell^2);
pb = Helper.ProgressBar();
pb.run('Propagating in imaginary time: ');
@ -41,7 +41,7 @@ function [psi] = propagateWavefunction(this, psi, Params, VParams, Transf, VDk,
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));
psi = psi.*exp(-1i*dt*(V + Ez + g_eff*abs(psi).^2 + gamma_eff*abs(psi).^3 + Params.gdd*Phi/(sqrt(2*pi)*VParams.ell) - muchem));
% kin
psi = fftn(psi);

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

@ -9,10 +9,15 @@ function [Params, Transf, psi, V, VDk] = run(this)
[Transf] = this.setupSpace(Params);
% --- Set up for Variational method ---
VarArray = Params.VarArray;
VParams.ell = VarArray(1);
VParams.nu = VarArray(2);
[VParams] = this.Calculator.find_nk_fit(VParams);
VParams.SelfConIter = Params.SelfConIter; % Max number of iterations to perform self-consistent calculation
VParams.ell = Params.ell; % initial [ell], ell is the "width" - psi ~ e^(z^2/ell^2)
% Window of optimization
VParams.ell_lower = Params.ell_lower;
VParams.ell_upper = Params.ell_upper;
% Relative cutoffs
VParams.ellcutoff = Params.ellcutoff;
% --- Initialize ---
mkdir(sprintf(this.SaveDirectory))
@ -21,24 +26,22 @@ function [Params, Transf, psi, V, VDk] = run(this)
[psi,V,VDk] = this.initialize(Params,VParams,Transf);
ells(1) = VarArray(1);
nus(1) = VarArray(2);
ells(1) = VParams.ell;
E_Var = @(x) this.Calculator.calculateVariationalEnergy(psi, Params, x, Transf, VDk, V)/Params.N;
E_vs_iter(1) = E_Var([ells(1) nus(1)]);
E_vs_iter(1) = E_Var(ells(1));
t_idx = 1; % Start at t = 0;
for nn = 1:Params.SelfConIter
Observ.EVec = []; Observ.NormVec = []; Observ.PCVec = []; Observ.tVecPlot = []; Observ.mucVec = []; Observ.residual = [];
Observ.res_idx = 1;
Observ.res_idx = 1;
VParams.ell = VarArray(1);
VParams.nu = VarArray(2);
[VParams] = this.Calculator.find_nk_fit(VParams);
[psi,V,VDk] = this.initialize(Params,VParams,Transf);
[~,V,VDk] = this.initialize(Params,VParams,Transf); % Do not overwrite psi, susbequent iterations should use psi generated at the end of the loop to converge faster
% This is still needed however to recalculate VDk with new value
% for the variational parameter
% --- Adding some noise ---
% 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
psi = sqrt(Params.N)*psi/sqrt(Norm);
r = normrnd(0,1,size(psi));
@ -52,24 +55,22 @@ function [Params, Transf, psi, V, VDk] = run(this)
% --- Constrained minimization ---
E_Var = @(x) this.Calculator.calculateVariationalEnergy(psi, Params, x, Transf, VDk, V)/Params.N;
VarArray = fmincon(E_Var,VarArray,[],[],[],[],[Params.ell_lower, Params.nu_lower],[Params.ell_upper, Params.nu_upper],[],fminconoptions);
VParams.ell = fmincon(E_Var,VParams.ell,[],[],[],[],Params.ell_lower,Params.ell_upper,[],fminconoptions);
% --- Convergence check ---
ells = [ells VarArray(1)];
nus = [nus VarArray(2)];
ells = [ells VParams.ell];
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)];
E_vs_iter = [E_vs_iter E_Var(VParams.ell)];
save(sprintf('./Data/Run_%03i/psi_gs_%i.mat',Params.njob),'psi','Observ','Transf','Params','VDk','V','VarArray');
save(sprintf('./Data/Run_%03i/psi_gs_%i.mat',Params.njob),'psi','Observ','Transf','Params','VDk','V','VParams');
if relelldiff < Params.ellcutoff && relnudiff < Params.nucutoff
if relelldiff < Params.ellcutoff
break
end
end
disp('Saving data...');
save(sprintf('./Data/Run_%03i/psi_gs.mat',Params.njob),'psi','Observ','Transf','Params','VDk','V','VarArray');
save(sprintf('./Data/Run_%03i/psi_gs.mat',Params.njob),'psi','Observ','Transf','Params','VDk','V','VParams');
disp('Save complete!');
delete(sprintf('./Data/Run_%03i/psi_gs_*.mat',Params.njob))

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

@ -52,20 +52,14 @@ function [Params] = setupParameters(this)
% ================ 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
Params.ell = 4; % initial [ell], ell is the "width" - psi ~ e^(z^2/ell^2)
% 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 ================ %