function E = calculateTotalEnergy(~,psi,Params,VParams,Transf,VDk,V)
    
    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);
    
    % Parameters
    KEop      = 0.5*(Transf.KX.^2+Transf.KY.^2);
    normfac   = Params.Lx*Params.Ly/numel(psi);
    
    % DDIs
    frho      = fftn(abs(psi).^2);
    Phi       = real(ifftn(frho.*VDk));
    Eddi      = 0.5*Params.gdd*Phi.*abs(psi).^2/(sqrt(2*pi)*VParams.ell);%
    Eddi      = trapz(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         = Ez*Params.N + Ekin + Epot + Eint + Eddi + Eqf;
end