236       real*8  eh,ph  local
283              eh = dble(coulomb_screened_s_e(real_mb(dn(1))))
294               eh = dble(0.5*seh*dv)
298            eh = eh*hfx_parameter
300            if (ispin.eq.1) eh = eh + eh
301            ph = 2.0d0*eh
352       real*8  eh,ph,ss,teh  local
526                        eh = dble(icoulomb_screened_small_s_e(
532                        eh = dble(icoulomb_screened_s_e(real_mb(dn(1))))
538                     eh = eh*hfx_parameter*ss
539                     if (ispin0.eq.1) eh = eh + eh
540                     ph = 2.0d0*eh
542                     ehfx = ehfx - eh
545      >               = dbl_mb(ehfx_orb(1,ms)+i2-1) - eh
549                        ehfx = ehfx - eh
552      >                  = dbl_mb(ehfx_orb(1,ms)+i2-1) - eh
673               eh = dble(0.5*seh*dv)
676               eh = eh*hfx_parameter*ss
679               if (ispin0.eq.1) eh = eh + eh
680               ph = 2.0d0*eh
684               ehfx = ehfx - eh
686               dbl_mb(ehfx_orb(1,ms)+i-1) =dbl_mb(ehfx_orb(1,ms)+i-1)-eh
702                  ehfx = ehfx - eh
705      >                                      - eh
767       real*8  eh,ph,ss  local
839                 eh = dble(coulomb_screened_s_e(real_mb(dn(1))))
846                  eh = dble(0.5*seh*dv)
850               eh = eh*hfx_parameter*ss
851               if (ispin0.eq.1) eh = eh + eh
852               ph = 2.0d0*eh
855               ehfx = ehfx - eh
857               dbl_mb(ehfx_orb(1,ms)+i-1) = dbl_mb(ehfx_orb(1,ms)+i-1)-eh
863                  ehfx = ehfx - eh
866      >                                      - eh
924       real*8  eh,ph  local
975                 eh = dble(coulomb_screened_s_e(real_mb(dn(1))))
988                  eh = dble(0.5*seh*dv)
992               eh = eh*hfx_parameter
994               if (ispin.eq.1) eh = eh + eh
995               ph = 2.0d0*eh
1081       real*8  eh,ph  local
1176                  eh = dble(icoulomb_screened_s_e(real_mb(dn(1))))
1183                  ehfx = ehfx - eh
1189                     ehfx = ehfx - eh
1280                eh = dble(0.5*seh*dv)
1284                ehfx = ehfx - eh
1299                   ehfx = ehfx - eh
1481       real*8  eh,ph  local
1556                 eh = dble(0.5*seh*lattice_omega())
1563                  eh = dble(0.5*seh*dv)
1570               eh = eh*hfx_parameter
1571               if (ispin0.eq.1) eh = eh + eh
1572               ph = 2.0d0*eh
1575               ehfx = ehfx - eh
1577               dbl_mb(ehfx_orb(1,ms)+i-1) = dbl_mb(ehfx_orb(1,ms)+i-1)-eh
1583                  ehfx = ehfx - eh
1586      >                                      - eh