function E = calculateVariationalEnergy(~, psi, Params, VarArray, Transf, VDk, V)
                               
    VParams.ell = VarArray(1);
    VParams.nu  = VarArray(2);
    
    g_eff       = Params.gs*VParams.nu/(2^(1+1/VParams.nu)*VParams.ell*gamma(1/VParams.nu));
    gamma_eff   = Params.gammaQF*2^(1/VParams.nu-1.5)*5^(-1/VParams.nu)*VParams.ell*gamma(1+1/VParams.nu)*( VParams.nu/(VParams.ell*gamma(1/VParams.nu)) )^(5/2);
    EVar        = VParams.nu^2*gamma(2-1/VParams.nu)/(8*VParams.ell^2*gamma(1/VParams.nu)) + 0.5*Params.gz*VParams.ell^2*gamma(3/VParams.nu)/gamma(1/VParams.nu);
    
    % Parameters
    KEop        = 0.5*(Transf.KX.^2+Transf.KY.^2);
    normfac     = Params.Lx*Params.Ly/numel(psi);
    
    % DDIs
    [VParams]   = find_nk_fit(VParams); % Not totally sure this should be here
    frho        = fftn(abs(psi).^2);
    Phi         = real(ifftn(frho.*VDk));
    Eddi        = 0.5*Params.gdd*Phi.*abs(psi).^2/(sqrt(2*pi)*VParams.ell_eff);%
    Eddi        = sum(Eddi(:))*Transf.dx*Transf.dy;
    
    % Kinetic energy
    Ekin        = KEop.*abs(fftn(psi)*normfac).^2;
    Ekin        = trapz(Ekin(:))*Transf.dkx*Transf.dky/(2*pi)^2;
    
    % Potential energy
    Epot        = V.*abs(psi).^2;
    Epot        = trapz(Epot(:))*Transf.dx*Transf.dy;
    
    % Contact interactions
    Eint        = 0.5*g_eff*abs(psi).^4;
    Eint        = trapz(Eint(:))*Transf.dx*Transf.dy;
    
    % Quantum fluctuations
    Eqf         = 0.4*gamma_eff*abs(psi).^5;
    Eqf         = trapz(Eqf(:))*Transf.dx*Transf.dy;
    
    E           = EVar*Params.N + Ekin + Epot + Eint + Eddi + Eqf;
end