function muchem = calculateChemicalPotential(~,psi,Params,Transf,VDk,V)

    %Parameters
    normfac = Params.Lx*Params.Ly*Params.Lz/numel(psi);
    KEop= 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
    
    % DDIs
    frho=fftn(abs(psi).^2);
    Phi=real(ifftn(frho.*VDk));
    
    Eddi = (Params.gdd*Phi.*abs(psi).^2);
    
    %Kinetic energy
    Ekin = KEop.*abs(fftn(psi)*normfac).^2;
    Ekin = trapz(Ekin(:))*Transf.dkx*Transf.dky*Transf.dkz/(2*pi)^3;
    
    %Potential energy
    Epot = V.*abs(psi).^2;
    
    %Contact interactions
    Eint = Params.gs*abs(psi).^4;
    
    %Quantum fluctuations
    Eqf = Params.gammaQF*abs(psi).^5;
    
    %Total energy
    muchem = Ekin + trapz(Epot(:) + Eint(:) + Eddi(:) + Eqf(:))*Transf.dx*Transf.dy*Transf.dz; %
    muchem = muchem / Params.N;
end