Latest scripts for testing of energy minimization via CGD.

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

@@ -23,6 +23,7 @@ OptionsStruct.UseApproximationForLHY  = true;
 OptionsStruct.IncludeDDICutOff        = true;
 OptionsStruct.CutoffType              = 'Cylindrical';
 OptionsStruct.SimulationMode          = 'EnergyMinimization'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution' | 'EnergyMinimization'
+OptionsStruct.MaxIterationsForCGD     = 10;
 OptionsStruct.TimeStepSize            = 1E-4;            % in s
 OptionsStruct.MinimumTimeStepSize     = 2E-10;           % in s
 OptionsStruct.TimeCutOff              = 2E6;             % in s

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

@@ -45,10 +45,10 @@ sim.Potential                         = pot.trap();
 
 OptionsStruct = struct;
 
-OptionsStruct.NumberOfAtoms           = 4E5;
+OptionsStruct.NumberOfAtoms           = 8E5;
 OptionsStruct.DipolarPolarAngle       = deg2rad(0);
 OptionsStruct.DipolarAzimuthAngle     = 0;
-OptionsStruct.ScatteringLength        = 87;
+OptionsStruct.ScatteringLength        = 85;
 
 AspectRatio                           = 2.0;
 HorizontalTrapFrequency               = 125;
@@ -63,7 +63,7 @@ OptionsStruct.IncludeDDICutOff        = true;
 OptionsStruct.CutoffType              = 'Cylindrical';
 OptionsStruct.SimulationMode          = 'ImaginaryTimeEvolution'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution'
 OptionsStruct.TimeStepSize            = 1E-4;            % in s
-OptionsStruct.MinimumTimeStepSize     = 2E-10;           % in s
+OptionsStruct.MinimumTimeStepSize     = 1E-10;            % in s
 OptionsStruct.TimeCutOff              = 2E6;             % in s
 OptionsStruct.EnergyTolerance         = 5E-10;
 OptionsStruct.ResidualTolerance       = 1E-08;
@@ -73,7 +73,7 @@ OptionsStruct.PlotLive                = false;
 OptionsStruct.JobNumber               = 0;
 OptionsStruct.RunOnGPU                = true;
 OptionsStruct.SaveData                = true;
-OptionsStruct.SaveDirectory           = './Results/Data_3D/BeyondSSD_Labyrinth';
+OptionsStruct.SaveDirectory           = './Results/Data_3D/ApproximateLHY/BeyondSSD_Labyrinth';
 options                               = Helper.convertstruct2cell(OptionsStruct);
 clear OptionsStruct
 

Dipolar-Gas-Simulator/+Simulator/@DipolarGas/DipolarGas.m

@@ -17,6 +17,7 @@ classdef DipolarGas < handle & matlab.mixin.Copyable
         EnergyTolerance;
         ResidualTolerance;
         NoiseScaleFactor;
+        MaxIterationsForCGD;
         
         Calculator;
         
@@ -69,6 +70,8 @@ classdef DipolarGas < handle & matlab.mixin.Copyable
                 @(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));
             addParameter(p, 'NoiseScaleFactor',       4,...
                 @(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));
+            addParameter(p, 'MaxIterationsForCGD',  100,...
+                @(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));
             
             addParameter(p, 'IncludeDDICutOff',       true,...
                 @islogical);
@@ -104,6 +107,7 @@ classdef DipolarGas < handle & matlab.mixin.Copyable
             this.EnergyTolerance      = p.Results.EnergyTolerance;
             this.ResidualTolerance    = p.Results.ResidualTolerance;
             this.NoiseScaleFactor     = p.Results.NoiseScaleFactor;
+            this.MaxIterationsForCGD  = p.Results.MaxIterationsForCGD;
             
             this.IncludeDDICutOff       = p.Results.IncludeDDICutOff;
             this.UseApproximationForLHY = p.Results.UseApproximationForLHY;

Dipolar-Gas-Simulator/+Simulator/@DipolarGas/minimizeEnergyFunctional.m

@@ -1,4 +1,4 @@
-function [psi]      = minimizeEnergyFunctional(this,psi,Params,Transf,VDk,V)
+function [psi]  = minimizeEnergyFunctional(this,psi,Params,Transf,VDk,V)
     
     format long;
     
@@ -8,35 +8,58 @@ function [psi]      = minimizeEnergyFunctional(this,psi,Params,Transf,VDk,V)
     % Convergence Criteria:
     Epsilon     = 1E-5;
     
-    % Gradient Calculation
-    J           = compute_gradient(psi, Params, Transf, VDk, V);  % Gradient at initial point
-    S           = -J;                     % Initial search direction
-    
     % Iteration Counter:
     i           = 1;
     
-    % Minimization
+    % Initialize the PrematureExitFlag to false
-    while norm(S(:)) > Epsilon % Halting Criterion
+    PrematureExitFlag = false;
+
+    % Minimization Loop
+    while true
+        % Compute gradient
+        J       = compute_gradient(psi, Params, Transf, VDk, V);
+
+        % Check stopping criterion (Gradient norm)
+        if norm(J(:)) < Epsilon
+            disp('Tolerance reached: Gradient norm is below the specified epsilon.');
+            PrematureExitFlag = true;  % Set flag to indicate premature exit
+            break;
+        elseif i >= this.MaxIterationsForCGD
+            disp('Maximum number of iterations for CGD reached.');
+            PrematureExitFlag = true;  % Set flag to indicate premature exit
+            break;
+        end
+
+        % Initialize search direction if first iteration
+        if i == 1
+            S   = -J;
+        else
+            % Update search direction
+            S   = update_search_direction(S, J, J_old);
+        end
+
         % Step Size Optimization (Line Search)
-        alpha  = optimize_step_size(f, psi, S, Params, Transf, VDk, V);
+        alpha   = optimize_step_size(f, psi, S, Params, Transf, VDk, V);
-    
+
-        % Update the solution
+        % Update solution
-        x_o_new = psi + alpha * S;
+        psi     = psi + alpha * S;
-    
+
-        % Update the gradient and search direction
+        % Normalize psi
+        Norm    = sum(abs(psi(:)).^2) * Transf.dx * Transf.dy * Transf.dz;
+        psi     = sqrt(Params.N) * psi / sqrt(Norm);
+
+        % Store old gradient
         J_old   = J;
-        J       = compute_gradient(x_o_new, Params, Transf, VDk, V);
-        S       = update_search_direction(S, J, J_old);
-    
-        % Update for next iteration
-        psi     = x_o_new;
-        Norm    = sum(abs(psi(:)).^2)*Transf.dx*Transf.dy*Transf.dz;
-        psi     = sqrt(Params.N)*psi/sqrt(Norm);
         i       = i + 1;
     end
     
-    disp('Minimum found!');
+    % Check if loop ended prematurely
-    fprintf('Number of Iterations for Convergence: %d\n\n', i);
+    if PrematureExitFlag
+        disp('Optimizer ended prematurely without convergence to a minimum.');
+    else
+        disp('Minimum found!');
+        fprintf('Number of Iterations for Convergence: %d\n\n', i);
+    end
     
     muchem      = this.Calculator.calculateChemicalPotential(psi,Params,Transf,VDk,V);
 
@@ -48,26 +71,31 @@ end
 
 %% Helper functions
 % Numerical Gradient Calculation using the finite differences method
-function J = compute_gradient(psi, Params, Transf, VDk, V)
+function J      = compute_gradient(psi, Params, Transf, VDk, V)
 
-    % Kinetic energy
+    % Operators
-    KEop    = 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
-    H_o     = KEop + V;
     
-    %Contact interactions
+    % Kinetic energy
-    Hint    = Params.gs*abs(psi).^2;
+    KEop    = 0.5 * (Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
+    HKin    = @(w) real(ifft(KEop.*fft(w)));
+
+    % Trap Potential
+    HV      = @(w) V.*w;
+
+    % Contact interactions
+    Hint    = @(w) (Params.gs*abs(psi).^2).*w;
     
     % DDIs
     frho    = fftn(abs(psi).^2);
     Phi     = real(ifftn(frho.*VDk));
-    Hddi    = Params.gdd*Phi;
+    Hddi    = @(w) (Params.gdd*Phi).*w;
-    
-    %Quantum fluctuations
-    Hqf     = Params.gammaQF*abs(psi).^3;
-    
-    H = H_o + Hint + Hddi + Hqf;
 
-    J = H.*psi;
+    % Quantum fluctuations
+    Hqf     = @(w) (Params.gammaQF*abs(psi).^3).*w;
+    
+    H       = @(w) HKin(w) + HV(w) + Hint(w) + Hddi(w) + Hqf(w);
+
+    J = H(psi);
 
 end
 
@@ -81,7 +109,7 @@ function alpha  = optimize_step_size(f, X, S, Params, Transf, VDk, V)
 
     grad        = compute_gradient(X, Params, Transf, VDk, V);  % Compute gradient once
     f_X         = f(X);                                         % Evaluate f(X) once
-    
+
     for k = 1:max_iter
         % Evaluate Armijo condition with precomputed f(X) and grad
         if f(X + alpha * S) <= f_X + c * alpha * (S(:)' * grad(:))
@@ -98,8 +126,13 @@ function alpha  = optimize_step_size(f, X, S, Params, Transf, VDk, V)
 
 end
 
-% Update Search Direction (Polak-Ribiere method)
+% Update Search Direction 
 function S_new  = update_search_direction(S, J_new, J_old)
-    beta        = (norm(J_new(:))^2) / (norm(J_old(:))^2);
+    % (Fletcher-Reeves method)
+    % beta        = (norm(J_new(:))^2) / (norm(J_old(:))^2);
+    % S_new       = -J_new + beta * S;
+
+    % (Polak-Ribiere method)
+    beta        = max(0, (J_new(:)' * (J_new(:) - J_old(:))) / (norm(J_old(:))^2));
     S_new       = -J_new + beta * S;
 end