function [psi]  = minimizeEnergyFunctional(this,psi,Params,Transf,VDk,V,Observ)
    
    format long;
    
    % Define the function handle
    f                 = @(X) this.Calculator.calculateTotalEnergy(X, Params, Transf, VDk, V)/Params.N;

    % Convergence Criteria:
    Epsilon           = 1E-5;
    
    % Iteration Counter:
    i                 = 1;
    Observ.residual   = 1; 
    Observ.res        = 1; 

    % Initialize the PrematureExitFlag to false
    PrematureExitFlag = false;
    
    if this.PlotLive
        Plotter.plotLive(psi,Params,Transf,Observ)
        drawnow
    end

    % 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);

        % Update solution
        psi     = psi + alpha * S;

        % 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;
        i       = i + 1;
        muchem = this.Calculator.calculateChemicalPotential(psi,Params,Transf,VDk,V);
                
        if mod(i,500) == 0        
    
            % Change in Energy
            E               = this.Calculator.calculateTotalEnergy(psi,Params,Transf,VDk,V);
            E               = E/Norm;
            Observ.EVec     = [Observ.EVec E];
    
            % Chemical potential
            Observ.mucVec   = [Observ.mucVec muchem];
    
            % 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)
                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);
    E              = E/Norm;
    Observ.EVec    = [Observ.EVec E];

    disp('Saving data...');
    save(sprintf(strcat(this.SaveDirectory, '/Run_%03i/psi_gs.mat'),Params.njob),'psi','muchem','Observ','Transf','Params','VDk','V');
    disp('Save complete!');
    
end

%% Helper functions
% Numerical Gradient Calculation using the finite differences method
function J      = compute_gradient(psi, Params, Transf, VDk, V)

    % Operators
    
    % Kinetic energy
    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    = @(w) (Params.gdd*Phi).*w;

    % 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

% Line Search (Step Size Optimization)
function alpha  = optimize_step_size(f, X, S, Params, Transf, VDk, V)
    alpha       = 1;        % Initial step size
    rho         = 0.5;      % Step size reduction factor
    c           = 1E-4;     % Armijo condition constant
    max_iter    = 100;      % Max iterations for backtracking
    tol         = 1E-4;     % Tolerance for stopping

    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(:))
            break;
        else
            alpha = rho * alpha;  % Reduce the step size
        end
        
        % Early stopping if step size becomes too small
        if alpha < tol
            break;
        end
    end

end

% Update Search Direction 
function S_new  = update_search_direction(S, J_new, J_old)
    % (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