function res  = calculateNormalizedResiduals(~,psi,Params,VParams,Transf,VDk,V,muchem)
    
    g_eff     = Params.gs      * (1/(sqrt(2*pi)*VParams.ell));
    gamma_eff = Params.gammaQF * (sqrt(2/5)/(pi^(3/4)*VParams.ell^(3/2)));
    Ez        = (0.25/VParams.ell^2) + (0.25*Params.gz*VParams.ell^2);
    
    KEop      = 0.5*(Transf.KX.^2+Transf.KY.^2);
    
    % DDIs
    frho      = fftn(abs(psi).^2);
    Phi       = real(ifftn(frho.*VDk));
    Eddi      = Params.gdd*Phi.*psi/(sqrt(2*pi)*VParams.ell);
    
    % Kinetic energy
    Ekin      = ifftn(KEop.*fftn(psi));
    
    % Potential energy
    Epot      = V.*psi;
    
    % Contact interactions
    Eint      = g_eff*abs(psi).^2.*psi;

    % Quantum fluctuations
    Eqf       = gamma_eff*abs(psi).^3.*psi;
    
    % Total energy
    res       = sum(abs(Ekin(:) + Ez*psi(:) + Epot(:) + Eint(:) + Eddi(:) + Eqf(:) - muchem*psi(:))*Transf.dx*Transf.dy)/sum(abs(muchem*psi(:))*Transf.dx*Transf.dy);

end