function muchem = calculateChemicalPotential(~,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
    normfac     = Params.Lx*Params.Ly/numel(psi);
    KEop        = 0.5*(Transf.KX.^2+Transf.KY.^2);
    
    % DDIs
    frho        = fftn(abs(psi).^2);
    Phi         = real(ifftn(frho.*VDk));
    Eddi        = (Params.gdd*Phi.*abs(psi).^2)/(sqrt(2*pi)*VParams.ell);
    
    % 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;
    
    % Contact interactions
    Eint        = g_eff*abs(psi).^4;
    
    % Quantum fluctuations
    Eqf         = gamma_eff*abs(psi).^5;
    
    % Total energy
    muchem      = Ekin + Ez*Params.N +  trapz(Epot(:) + Eint(:) + Eddi(:) + Eqf(:))*Transf.dx*Transf.dy; %
    muchem      = muchem / Params.N; %Only use if psi is normalized to N
end