function [psi]  = runGradientDescent(this,psi,Params,Transf,VDk,V,Observ)
    
    format long;
    switch this.GradientDescentMethod
        case 'HeavyBall'
            % Convergence Criteria:
            epsilon         = 1E-6;
            alpha           = 1E-3;
            beta            = 0.9;
        
            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,this.SimulationMode)
                drawnow
            end
        
            % Minimization Loop
            for idx      = 1:this.MaxIterationsForGD
                
                % Compute gradient
                J        = compute_gradient(psi, Params, Transf, VDk, V);
        
                % Calculate chemical potential
                muchem   = sum(real(psi(:)' .* J(:))) / sum(abs(psi(:)).^2);
                
                % Calculate residual and check convergence
                residual = sum(abs(J(:) - (muchem * psi(:))).^2) * Transf.dx * Transf.dy * Transf.dz;
                
                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     = sum(abs(psi_new(:)).^2) * Transf.dx * Transf.dy * Transf.dz;
                    psi      = sqrt(Params.N) * 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);
                        E               = E/Norm;
                        Observ.EVec     = [Observ.EVec E];
                
                        % Collect chemical potentials
                        Observ.mucVec   = [Observ.mucVec muchem];
                
                        % Collect residuals
                        Observ.residual = [Observ.residual residual];
                        Observ.res_idx  = Observ.res_idx + 1;
                        
                        if this.PlotLive
                            Plotter.plotLive(psi,Params,Transf,Observ,this.SimulationMode)
                            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);
            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!');
        
        case 'NonLinearCGD'
            
            % Convergence Criteria:
            epsilon           = 1E-14;
            
            % Iteration Counter:
            i                 = 1;
            Observ.residual   = 1; 
            Observ.res        = 1; 
        
            % Initialize the PrematureExitFlag to false
            PrematureExitFlag = false;
            
            % Live plotter
            if this.PlotLive
                Plotter.plotLive(psi,Params,Transf,Observ,this.SimulationMode)
                drawnow
            end
            
            % Minimization Loop
            while true
                % Compute gradient
                J                = compute_gradient(psi, Params, Transf, VDk, V);
                % Calculate chemical potential
                muchem           = real(psi(:)' * J(:)) / sum(abs(psi(:)).^2);
                % Calculate residual
                residual         = J - (muchem * psi);
                
                % Compute energy difference between the last two saved energy values
                if i == 1
                    energydifference = NaN;
                elseif mod(i,100) == 0 && length(Observ.EVec) > 1
                    energydifference = abs(Observ.EVec(end) - Observ.EVec(end-1));
                end
                % Convergence check - if energy difference is below set tolerance, then exit
                if energydifference <= epsilon
                    disp('Tolerance reached: Energy difference is below the specified epsilon.');
                    PrematureExitFlag = true;  % Set flag to indicate premature exit
                    break;
                elseif i >= this.MaxIterationsForGD % If set maximum number of iterations reached, then exit
                    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
                    d            = -residual;
                else % Compute beta via Polak-Ribiere and create new direction
                    residual_new = residual;
                    beta         = compute_beta(residual_new, residual_old);
                    d            = -residual_new + beta * p_old;
                end
                
                residual_old     = residual;
                p                = d - (real(d(:)' * psi(:) * Transf.dx * Transf.dy * Transf.dz) .* 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(sum(abs(p(:)).^2));
                psi              = (cos(theta).*psi) + (gamma*sin(theta).*p);
                % Normalize psi
                Norm             = sum(abs(psi(:)).^2) * Transf.dx * Transf.dy * Transf.dz;
                psi              = sqrt(Params.N) * psi / sqrt(Norm);
                i                = i + 1;

                % Calculate chemical potential with new psi
                muchem           = real(psi(:)' * J(:)) / sum(abs(psi(:)).^2);
                
                if mod(i,100) == 0        
            
                    % Collect Energy value
                    E               = this.Calculator.calculateTotalEnergy(psi,Params,Transf,VDk,V);
                    E               = E/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;
                    
                    % Live plotter
                    if this.PlotLive
                        Plotter.plotLive(psi,Params,Transf,Observ,this.SimulationMode)
                        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
end

%% Modules

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) 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     = 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

function g      = compute_g(psi, p, Params, VDk)
    
    rho      = real(psi.*conj(p));

    % Contact interactions
    C        = Params.gs*Params.N;
    gint     = @(w)(C.*rho).*w;
    
    % DDIs
    D        = Params.gdd*Params.N;
    rhotilde = fftn(rho);
    Phi      = ifftn(rhotilde.*VDk);
    gaddi    = @(w)(D.*Phi).*w;

    % Quantum fluctuations
    eps_dd   = Params.add/Params.as;
    Q        = (4/(3*pi^2)) * (C^(5/2)/Params.N) * (1 + ((3*eps_dd^2)/2));
    gqf      = @(w)(((3/2)*Q).*abs(psi).*rho).*w;
    
    gop      = @(w) gint(w) + gaddi(w) + gqf(w);

    g        = gop(psi);

end

% Optimal direction via line search
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(sum(abs(p(:)).^2));
    numerator   = real(p(:)' * Hpsi(:) * Transf.dx * Transf.dy * Transf.dz)/gamma;
    denominator = muchem - (gamma^2 * (real(p(:)' * Hp(:) * Transf.dx * Transf.dy * Transf.dz) + real(g(:)' * p(:) * Transf.dx * Transf.dy * Transf.dz))); 

    theta       = numerator / denominator;
end

% Optimal step size via Polak-Ribiere
function beta   = compute_beta(residual_new, residual_old)
    beta        = (residual_new(:)' * (residual_new(:) - residual_old(:))) / sum(abs(residual_old(:)).^2);
    beta        = max(0,beta);
end