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

    % 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
    while norm(S(:)) > Epsilon % Halting Criterion
        % Step Size Optimization (Line Search)
        alpha  = optimize_step_size(f, psi, S, Params, Transf, VDk, V);
    
        % Update the solution
        x_o_new = psi + alpha * S;
    
        % Update the gradient and search direction
        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!');
    fprintf('Number of Iterations for Convergence: %d\n\n', i);
    
    muchem      = this.Calculator.calculateChemicalPotential(psi,Params,Transf,VDk,V);

    disp('Saving data...');
    save(sprintf(strcat(this.SaveDirectory, '/Run_%03i/psi_gs.mat'),Params.njob),'psi','muchem','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)

    % Kinetic energy
    KEop    = 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
    H_o     = KEop + V;
    
    %Contact interactions
    Hint    = Params.gs*abs(psi).^2;
    
    % DDIs
    frho    = fftn(abs(psi).^2);
    Phi     = real(ifftn(frho.*VDk));
    Hddi    = Params.gdd*Phi;
    
    %Quantum fluctuations
    Hqf     = Params.gammaQF*abs(psi).^3;
    
    H = H_o + Hint + Hddi + Hqf;

    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-8;     % 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 (Polak-Ribiere method)
function S_new  = update_search_direction(S, J_new, J_old)
    beta        = (norm(J_new(:))^2) / (norm(J_old(:))^2);
    S_new       = -J_new + beta * S;
end