Added and tested 2D BdG solver - works until Matlab starts solving for the eigen-values and vectors.

Dipolar-Gas-Simulator/+BdGSolver2D/BdGMatrix.m

@@ -0,0 +1,35 @@
+function BdGVec = BdGMatrix(gVec, psi, Params, VDk, VParams, Transf, muchem)
+    
+    gs           = Params.gs; 
+    gdd          = Params.gdd;
+    gammaQF      = Params.gammaQF; 
+    
+    KEop         = 0.5*(Transf.KX.^2+Transf.KY.^2);
+    
+    g_pf_2D      = 1/(sqrt(2*pi)*VParams.ell);
+    gQF_pf_2D    = sqrt(2/5)/(pi^(3/4)*VParams.ell^(3/2));
+    Ez           = (0.25/VParams.ell^2) + (0.25*Params.gz*VParams.ell^2);
+    muchem_tilde = muchem - Ez;
+    
+    % eigs only works with column vectors
+    psi          = psi.';
+    KEop         = KEop.';
+    VDk          = VDk.';
+    
+    % Interaction Potential
+    frho         = fftn(abs(psi).^2);
+    Phi          = real(ifftn(frho.*VDk));
+    % Operators
+    H            = @(w) real(ifft(KEop.*fft(w)));
+    C            = @(w) (((g_pf_2D*gs*abs(psi).^2) + (g_pf_2D*gdd*Phi)).*w) + (gQF_pf_2D*gammaQF*abs(psi).^3.*w);
+    muHC         = @(w) (-muchem_tilde * w) + H(w) + C(w);
+    X            = @(w) (psi.*real(ifft(VDk.*fft(psi.*w)))) + (3/2)*(gQF_pf_2D*gammaQF*abs(psi).^3).*w;
+    
+    g            = reshape(gVec, size(psi));
+
+    % Operate in order on g
+    BdG          = muHC(muHC(g) + (2.*X(g)));
+
+    BdGVec       = reshape(BdG, [], 1);
+
+end

Dipolar-Gas-Simulator/+BdGSolver2D/solveBogoliubovdeGennesIn2D.m

@@ -1,44 +1,16 @@
 function [evals, modes] = solveBogoliubovdeGennesIn2D(psi, Params, VDk, VParams, Transf, muchem)
 
-gs           = Params.gs; 
+Size         = length(psi(:));
-gdd          = Params.gdd;
+Neigs        = length(psi(:));
-gammaQF      = Params.gammaQF; 
-
-KEop         = 0.5*(Transf.KX.^2+Transf.KY.^2);
-
-g_pf_2D      = 1/(sqrt(2*pi)*VParams.ell);
-gQF_pf_2D    = sqrt(2/5)/(pi^(3/4)*VParams.ell^(3/2));
-Ez           = (0.25/VParams.ell^2) + (0.25*Params.gz*VParams.ell^2);
-muchem_tilde = muchem - Ez;
-
-% eigs only works with column vectors
-psi          = psi.';
-KEop         = KEop.';
-VDk          = VDk.';
-
-% Interaction Potential
-frho         = fftn(abs(psi).^2);
-Phi          = real(ifftn(frho.*VDk));
-
-% Operators
-H            = @(w)     real(ifft(KEop.*fft(w)));
-C            = @(w)     (((g_pf_2D*gs*abs(psi).^2) + (g_pf_2D*gdd*Phi)).*w) + (gQF_pf_2D*gammaQF*abs(psi).^3.*w);
-muHC         = @(w)     (-muchem_tilde * w) + H(w) + C(w);
-X            = @(w,psi) (psi.*real(ifft(VDk.*fft(psi.*w)))) + (3/2)*(gQF_pf_2D*gammaQF*abs(psi).^3).*w;
-
-% Operate in order on g
-BdG          = @(g) muHC(muHC(g) + (2.*X(g)));
-
-syssize      = size(psi);
-opts.v0      = psi(:);
 opts.tol     = 1e-16;
 opts.disp    = 1;
 opts.issym   = 0;
 opts.isreal  = 1;
 opts.maxit   = 1e4;
-Neigs        = syssize;
 
-[g,D]        = eigs(BdG,syssize,Neigs,'sr',opts);
+BdGVec       = @(g) BdGSolver2D.BdGMatrix(g, psi, Params, VDk, VParams, Transf, muchem);
+
+[g,D]        = eigs(BdGVec,Size,Neigs,'sr',opts);
 evals        = diag(D); 
 clear D;
 

Dipolar-Gas-Simulator/+Scripts/run_locally.m

@@ -88,7 +88,7 @@ OptionsStruct.WidthUpperBound         = 12;
 OptionsStruct.WidthCutoff             = 1e-2;
 
 OptionsStruct.PlotLive                = true;
-OptionsStruct.JobNumber               = 2;
+OptionsStruct.JobNumber               = 1;
 OptionsStruct.RunOnGPU                = false;
 OptionsStruct.SaveData                = true;
 OptionsStruct.SaveDirectory           = './Data_TriangularPhase';
@@ -100,7 +100,33 @@ pot                                   = VariationalSolver2D.Potentials(options{:
 solver.Potential                      = pot.trap();
 
 %-% Run Solver %-%
-[Params, Transf, psi, V, VDk]         = solver.run();
+% [Params, Transf, psi, V, VDk]         = solver.run();
+
+% Solve BdG equations
+% Load data
+Data = load(sprintf(horzcat(solver.SaveDirectory, '/Run_%03i/psi_gs.mat'),solver.JobNumber),'psi','Params','Transf');
+Params = Data.Params;
+Transf = Data.Transf; 
+ 
+if isgpuarray(Data.psi)
+    psi = gather(Data.psi);
+else
+    psi = Data.psi;
+end
+
+VParams.ell    = Params.ell;
+
+% == DDI potential == %
+VDk = solver.Calculator.calculateVDkWithCutoff(Transf, Params, VParams.ell);
+
+% == Trap potential == %
+X = Transf.X; Y = Transf.Y;
+V = 0.0*(Params.gx.*X.^2+Params.gy.*Y.^2); 
+
+% == Chemical potential == %
+muchem = solver.Calculator.calculateChemicalPotential(psi,Params,VParams,Transf,VDk,V);
+
+[evals, modes] = BdGSolver2D.solveBogoliubovdeGennesIn2D(psi, Params, VDk, VParams, Transf, muchem);
 
 %% - Create Variational2D and Calculator object with specified options
 
@@ -130,7 +156,7 @@ OptionsStruct.WidthUpperBound         = 12;
 OptionsStruct.WidthCutoff             = 1e-2;
 
 OptionsStruct.PlotLive                = true;
-OptionsStruct.JobNumber               = 2;
+OptionsStruct.JobNumber               = 1;
 OptionsStruct.RunOnGPU                = false;
 OptionsStruct.SaveData                = true;
 OptionsStruct.SaveDirectory           = './Data_StripePhase';
@@ -172,7 +198,7 @@ OptionsStruct.WidthUpperBound         = 12;
 OptionsStruct.WidthCutoff             = 1e-2;
 
 OptionsStruct.PlotLive                = true;
-OptionsStruct.JobNumber               = 2;
+OptionsStruct.JobNumber               = 1;
 OptionsStruct.RunOnGPU                = false;
 OptionsStruct.SaveData                = true;
 OptionsStruct.SaveDirectory           = './Data_HoneycombPhase';
@@ -212,3 +238,4 @@ SaveDirectory = './Data_HoneycombPhase';
 JobNumber     = 3;
 % Plotter.visualizeGSWavefunction2D(SaveDirectory, JobNumber)
 Scripts.analyzeRun2D(SaveDirectory, JobNumber)
+

Dipolar-Gas-Simulator/+Scripts/run_on_cluster.m

@@ -10,6 +10,7 @@
 % as = ((as/add)*Params.add)/Params.a0
 % Critical point: 102.5133; Triangular phase: 98.0676; Stripe phase: 100.0289; Honeycomb phase: 101.9903
 
+%{
 %% - Create Variational2D and Calculator object with specified options
 
 OptionsStruct = struct;
@@ -135,3 +136,71 @@ solver.Potential                      = pot.trap();
 
 %-% Run Solver %-%
 [Params, Transf, psi, V, VDk]         = solver.run();
+%}
+%% - Create Variational2D and Calculator object with specified options
+
+OptionsStruct = struct;
+
+OptionsStruct.NumberOfAtoms           = 2.0030e+07;     
+OptionsStruct.DipolarPolarAngle       = 0;
+OptionsStruct.DipolarAzimuthAngle     = 0;
+OptionsStruct.ScatteringLength        = 98.0676;        
+
+OptionsStruct.TrapFrequencies         = [10, 10, 72.4];
+OptionsStruct.TrapPotentialType       = 'None';
+
+OptionsStruct.NumberOfGridPoints      = [128, 128];
+OptionsStruct.Dimensions              = [100, 100];      
+OptionsStruct.TimeStepSize            = 500E-6;          % in s
+OptionsStruct.MinimumTimeStepSize     = 1E-5;            % in s
+OptionsStruct.TimeCutOff              = 1E6;             % in s
+OptionsStruct.EnergyTolerance         = 5E-10;
+OptionsStruct.ResidualTolerance       = 1E-04;
+OptionsStruct.NoiseScaleFactor        = 4;
+
+OptionsStruct.MaxIterations           = 1;      
+OptionsStruct.VariationalWidth        = 5.0;       
+OptionsStruct.WidthLowerBound         = 0.2;
+OptionsStruct.WidthUpperBound         = 12;
+OptionsStruct.WidthCutoff             = 1e-2;
+
+OptionsStruct.PlotLive                = true;
+OptionsStruct.JobNumber               = 1;
+OptionsStruct.RunOnGPU                = false;
+OptionsStruct.SaveData                = true;
+OptionsStruct.SaveDirectory           = './Data_TriangularPhase';
+options                               = Helper.convertstruct2cell(OptionsStruct);
+clear OptionsStruct
+
+solver                                = VariationalSolver2D.DipolarGas(options{:});
+pot                                   = VariationalSolver2D.Potentials(options{:});
+solver.Potential                      = pot.trap();
+
+%-% Run Solver %-%
+% [Params, Transf, psi, V, VDk]         = solver.run();
+
+% Solve BdG equations
+% Load data
+Data = load(sprintf(horzcat(solver.SaveDirectory, '/Run_%03i/psi_gs.mat'),solver.JobNumber),'psi','Params','Transf');
+Params = Data.Params;
+Transf = Data.Transf; 
+ 
+if isgpuarray(Data.psi)
+    psi = gather(Data.psi);
+else
+    psi = Data.psi;
+end
+
+VParams.ell    = Params.ell;
+
+% == DDI potential == %
+VDk = solver.Calculator.calculateVDkWithCutoff(Transf, Params, VParams.ell);
+
+% == Trap potential == %
+X = Transf.X; Y = Transf.Y;
+V = 0.0*(Params.gx.*X.^2+Params.gy.*Y.^2); 
+
+% == Chemical potential == %
+muchem = solver.Calculator.calculateChemicalPotential(psi,Params,VParams,Transf,VDk,V);
+
+[evals, modes] = BdGSolver2D.solveBogoliubovdeGennesIn2D(psi, Params, VDk, VParams, Transf, muchem);