diff --git a/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/runGradientDescent.m b/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/runGradientDescent.m
index e115ee0..562c07b 100644
--- a/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/runGradientDescent.m
+++ b/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/runGradientDescent.m
@@ -8,9 +8,9 @@ function [psi]  = runGradientDescent(this,psi,Params,Transf,VDk,V,Observ)
     switch this.GradientDescentMethod
         case 'HeavyBall'
             % Convergence Criteria:
-            epsilon         = 1E-6;
-            alpha           = 1E-3;
+            alpha           = 1E-6;
             beta            = 0.9;
+            epsilon         = 1E-8;
         
             Observ.residual = 1; 
             Observ.res      = 1;
@@ -46,7 +46,7 @@ function [psi]  = runGradientDescent(this,psi,Params,Transf,VDk,V,Observ)
                     
                     % Normalize psi
                     Norm     = inner_product(psi_new, psi_new, Transf);
-                    psi      = sqrt(Params.N) * psi_new / sqrt(Norm);
+                    psi      = psi_new / sqrt(Norm);
                     
                     % Write output at specified intervals
                     if mod(idx,100) == 0        
@@ -156,7 +156,7 @@ function [psi]  = runGradientDescent(this,psi,Params,Transf,VDk,V,Observ)
                 
                 % Normalize psi
                 Norm             = abs(inner_product(psi, psi, Transf));
-                psi              = sqrt(Params.N) * psi / sqrt(Norm);
+                psi              = psi / sqrt(Norm);
                 i                = i + 1;
 
                 % Calculate chemical potential with new psi
@@ -208,59 +208,67 @@ 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;
 
-    % Operators
-    
-    % Kinetic energy
-    KEop    = 0.5 * (Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
-    HKin    = @(w) ifft(KEop.*fft(w));
+    % === Kinetic energy operator: -∇²/2 ===
+    KEop      = 0.5 * (Transf.KX.^2 + Transf.KY.^2 + Transf.KZ.^2);
+    HKin      = @(w) ifftn(KEop .* fftn(w));
 
-    % Trap Potential
-    HV      = @(w) V.*w;
+    % === External potential (already scaled by ℏω₀) ===
+    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;
+    % === Contact interaction: C|ψ|² ===
+    C         = Params.gs * N;
+    Hint      = @(w) (C * absPsi2) .* w;
 
-    % Quantum fluctuations
-    Hqf     = @(w) (Params.gammaQF*abs(psi).^3).*w;
-    
-    H       = @(w) HKin(w) + HV(w) + Hint(w) + Hddi(w) + Hqf(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;
 
-    J       = H(psi);
+    % === 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)
-    
-    rho      = real(psi.*conj(p));
+    % === Precompute necessary quantities ===
+    rho       = real(psi .* conj(p));                  % Mixed density
+    absPsi    = abs(psi);                              % Magnitude of wavefunction
+    N         = Params.N;
 
-    % 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;
+    % === Contact interaction term ===
+    C         = Params.gs * N;
+    g_contact = @(w) (C * rho) .* 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);
+    % === 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;
 
-    g        = gop(psi);
+    % === 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 via line search
+% Optimal direction
 function theta  = compute_optimal_theta(p, muchem, psi, Params, Transf, VDk, V)
     
     Hpsi        = compute_gradient(psi, Params, Transf, VDk, V);