function [psi]  = runGradientDescent(this,psi,Params,Transf,VDk,V,Observ)
    
    format long;

    AdditionalParams.GradientDescentMethod = this.GradientDescentMethod;
    SimulationMode                         = this.SimulationMode;

    switch this.GradientDescentMethod
        case 'HeavyBall'
            % Convergence Criteria:
            alpha           = 1E-5;
            beta            = 0.9;
            epsilon         = 1E-6;
        
            Observ.residual = 1; 
            Observ.res      = 1;
            psi_old         = psi;         % Previous psi value (for heavy-ball method)
            
            % Live Plotter
            if this.PlotLive
                Plotter.plotLive(psi,Params,Transf,Observ,SimulationMode,AdditionalParams)
                drawnow
            end
        
            % Minimization Loop
            for idx      = 1:this.MaxIterationsForGD
                
                % Compute gradient
                J            = compute_gradient(psi, Params, Transf, VDk, V);
        
                % Calculate chemical potential
                muchem       = real(inner_product(J, psi, Transf)) / inner_product(psi, psi, Transf);
                
                % Calculate residual and check convergence
                diff         = J - (muchem * psi);
                residual     = inner_product(diff, diff, Transf);
                
                if residual < epsilon
                    fprintf('Convergence reached at iteration %d\n', idx);
                    break;
                else
                
                    % Update psi using heavy-ball method
                    psi_new  = (1 + beta) * psi - alpha * J - beta * psi_old;
                    psi_old  = psi;
                    
                    % Normalize psi
                    Norm     = inner_product(psi_new, psi_new, Transf);
                    psi      = psi_new / sqrt(Norm);
                    
                    % Write output at specified intervals
                    if mod(idx,100) == 0        
                
                        % Collect change in energy
                        E               = this.Calculator.calculateTotalEnergy(psi,Params,Transf,VDk,V) / Norm;
                        Observ.EVec     = [Observ.EVec E];
                
                        % Collect chemical potentials
                        Observ.mucVec   = [Observ.mucVec muchem];
                
                        % Collect Normalized residuals
                        res             = this.Calculator.calculateNormalizedResiduals(psi,Params,Transf,VDk,V,muchem);
                        Observ.residual = [Observ.residual res];
                        Observ.res_idx  = Observ.res_idx + 1;
                        
                        if this.PlotLive
                            Plotter.plotLive(psi,Params,Transf,Observ,SimulationMode,AdditionalParams)
                            drawnow
                        end
            
                        save(sprintf(strcat(this.SaveDirectory, '/Run_%03i/psi_gs.mat'),Params.njob),'psi','muchem','Observ','Transf','Params','VDk','V');
                    end
                end
            end
            
            % Check if max iterations were hit without convergence
            last_iteration_number = idx;
            if last_iteration_number == this.MaxIterationsForGD
                fprintf('Max iterations reached without convergence. Final chemical potential: %.6f, Residual: %.6f\n', muchem, residual);
            else
                fprintf('Converged in %d iterations. Final chemical potential: %.6f\n', last_iteration_number, muchem);
            end
            
            % Change in Energy
            E              = this.Calculator.calculateTotalEnergy(psi,Params,Transf,VDk,V) / Norm;
            Observ.EVec    = [Observ.EVec E];
        
            save(sprintf(strcat(this.SaveDirectory, '/Run_%03i/psi_gs.mat'),Params.njob),'psi','muchem','Observ','Transf','Params','VDk','V');
            disp('Completed and saved!');
        
        case 'NonLinearCGD'
            
            % Convergence Criteria:
            epsilon           = 1E-13;
            
            % Iteration Counter:
            i                 = 1;
            Observ.residual   = 1; 
            Observ.res        = 1; 
            Observ.theta      = NaN;
        
            % Initialize the PrematureExitFlag to false
            PrematureExitFlag = false;
            
            % Live plotter
            if this.PlotLive
                Plotter.plotLive(psi,Params,Transf,Observ,SimulationMode,AdditionalParams)
                drawnow
            end
            
            % Minimization Loop
            while true
                % Compute gradient
                J                = compute_gradient(psi, Params, Transf, VDk, V);
                % Calculate chemical potential
                muchem           = real(inner_product(J, psi, Transf)) / inner_product(psi, psi, Transf);
                % Calculate residual
                residual         = J - (muchem * psi);
                
                % Energy convergence check every 100 steps
                if mod(i,100) == 0 && length(Observ.EVec) > 1
                    energydifference = abs(Observ.EVec(end) - Observ.EVec(end-1));
                    if energydifference <= epsilon
                        disp('Tolerance reached: Energy difference is below the specified epsilon.');
                        PrematureExitFlag = true;
                        break;
                    end
                end
            
                if i >= this.MaxIterationsForGD
                    disp('Maximum number of iterations for CGD reached.');
                    PrematureExitFlag = true;
                    break;
                end
                
                % Initialize search direction if first iteration
                if i == 1
                    d            = -residual;
                else % Compute beta via Polak-Ribiere and create new direction
                    residual_new = residual;
                    beta         = compute_beta(residual_new, residual_old, Transf);
                    d            = -residual_new + beta * p_old;
                end
                
                residual_old     = residual;
                proj             = inner_product(d, psi, Transf);
                p                = d - (proj * psi);
                p_old            = p;
                
                % Compute optimal theta to generate psi in the direction of minimum in the energy landscape
                theta            = compute_optimal_theta(p, muchem, psi, Params, Transf, VDk, V);
                
                % Update solution
                gamma            = 1 / sqrt(inner_product(p, p, Transf));
                psi              = (cos(theta).*psi) + (sin(theta).*(p*gamma));
                
                % Normalize psi
                Norm             = abs(inner_product(psi, psi, Transf));
                psi              = psi / sqrt(Norm);
                i                = i + 1;

                % Calculate chemical potential with new psi
                J                = compute_gradient(psi, Params, Transf, VDk, V);
                muchem           = real(inner_product(J, psi, Transf)) / Norm;
                
                if mod(i,100) == 0        
            
                    % Collect Energy value
                    E               = this.Calculator.calculateTotalEnergy(psi,Params,Transf,VDk,V) / Norm;
                    Observ.EVec     = [Observ.EVec E];
            
                    % Collect Chemical potential value
                    Observ.mucVec   = [Observ.mucVec muchem];
            
                    % Collect Normalized residuals
                    res             = this.Calculator.calculateNormalizedResiduals(psi,Params,Transf,VDk,V,muchem);
                    Observ.residual = [Observ.residual res];
                    Observ.res_idx  = Observ.res_idx + 1;

                     % Collect calculated thetas
                    Observ.theta    = [Observ.theta theta];
                    
                    % Live plotter
                    if this.PlotLive
                        Plotter.plotLive(psi,Params,Transf,Observ,SimulationMode,AdditionalParams)
                        drawnow
                    end
                    save(sprintf(strcat(this.SaveDirectory, '/Run_%03i/psi_gs.mat'),Params.njob),'psi','muchem','Observ','Transf','Params','VDk','V');
                end
            end
            
            % Check if loop ended prematurely
            if PrematureExitFlag
                disp('Optimizer ended prematurely without convergence to a minimum.');
            else
                fprintf('Minimum found! Number of Iterations for Convergence: %d\n\n', i);
            end
            
            % Change in Energy
            E              = this.Calculator.calculateTotalEnergy(psi,Params,Transf,VDk,V) / Norm;
            Observ.EVec    = [Observ.EVec E];
        
            save(sprintf(strcat(this.SaveDirectory, '/Run_%03i/psi_gs.mat'),Params.njob),'psi','muchem','Observ','Transf','Params','VDk','V');
            disp('Completed and saved!');
    end
end

%% Modules

function J      = compute_gradient(psi, Params, Transf, VDk, V)
    % === Precompute quantities ===
    absPsi2   = abs(psi).^2;        % Density |ψ|^2
    absPsi3   = abs(psi).^3;        % |ψ|^3
    N         = Params.N;

    % === Kinetic energy operator: -∇²/2 ===
    KEop      = 0.5 * (Transf.KX.^2 + Transf.KY.^2 + Transf.KZ.^2);
    HKin      = @(w) ifftn(KEop .* fftn(w));

    % === External potential (already scaled by ℏω₀) ===
    HV        = @(w) V .* w;

    % === Contact interaction: C|ψ|² ===
    C         = Params.gs * N;
    Hint      = @(w) (C * absPsi2) .* w;

    % === Dipole-dipole interaction: D U_dd * |ψ|² ===
    frho      = fftn(absPsi2);
    Phi_dd    = ifftn(frho .* VDk);
    D         = Params.gdd * N;
    Hddi      = @(w) (D * Phi_dd) .* w;

    % === Quantum fluctuations: Q|ψ|³ ===
    Q         = Params.gammaQF * N^(3/2);
    Hqf       = @(w) (Q * absPsi3) .* w;

    % === Total Hamiltonian operator ===
    H         = @(w) HKin(w) + HV(w) + Hint(w) + Hddi(w) + Hqf(w);

    % === Compute gradient ===
    J         = H(psi);
end

function g      = compute_g(psi, p, Params, VDk)
    % === Precompute necessary quantities ===
    rho       = real(psi .* conj(p));                  % Mixed density
    absPsi    = abs(psi);                              % Magnitude of wavefunction
    N         = Params.N;

    % === Contact interaction term ===
    C         = Params.gs * N;
    g_contact = @(w) (C * rho) .* w;

    % === Dipole-dipole interaction term ===
    rhotilde  = fftn(rho);
    Phi_dd    = ifftn(rhotilde .* VDk);                % Convolution with DDI kernel
    D         = Params.gdd * N;
    g_ddi     = @(w) (D * Phi_dd) .* w;

    % === Quantum fluctuation (LHY) correction ===
    Q         = Params.gammaQF * N^(3/2);              % LHY prefactor
    g_qf      = @(w) ((3/2) * Q * absPsi .* rho) .* w;

    % === Total effective operator applied to psi ===
    g_total   = @(w) g_contact(w) + g_ddi(w) + g_qf(w);

    % === Multiply by factor of 2, as per gradient expression ===
    g         = 2 * g_total(psi);
end

% Optimal direction
function theta  = compute_optimal_theta(p, muchem, psi, Params, Transf, VDk, V)
    
    Hpsi        = compute_gradient(psi, Params, Transf, VDk, V);
    Hp          = compute_gradient(p, Params, Transf, VDk, V);
    g           = compute_g(psi, p, Params, VDk);
    gamma       = 1 / sqrt(inner_product(p, p, Transf));
    numerator   = gamma * real(inner_product(p, Hpsi, Transf));
    denominator = muchem - (gamma^2 * (real(inner_product(p, Hp, Transf)) + real(inner_product(g, p, Transf)))); 

    theta       = numerator / denominator;
end

% Optimal step size via Polak-Ribiere
function beta   = compute_beta(residual_new, residual_old, Transf)
    beta        = inner_product(residual_new, (residual_new - residual_old), Transf) / inner_product(residual_old, residual_old, Transf);
    beta        = max(0, real(beta));
end

function s      = inner_product(u, v, Transf)
    s           = sum(conj(u(:)) .* v(:)) * Transf.dx * Transf.dy * Transf.dz;
end