diff --git a/src/bse/K_driver_init.F b/src/bse/K_driver_init.F
index 1750b595fe..a0ecb75a7e 100644
--- a/src/bse/K_driver_init.F
+++ b/src/bse/K_driver_init.F
@@ -59,11 +59,13 @@ subroutine K_driver_init(what,iq,Ken,Xk)
  endif
  !
  if (iq/=1) then
-   ! Default is Lfull as Lbar is NOT defined and gives WRONG results
+   ! For finite momentum Lbar is NOT defined and gives WRONG results
    ! BSE Hamiltonian is fully symmetric with Lfull
-   if(STRING_match(BSE_L_kind,"default")) BSE_L_kind="default-full"
-   if(STRING_match(BSE_L_kind,"bar")) call error("Lbar not defined at q/=0")
-endif
+   if (.not. STRING_match(BSE_L_kind, "full")) then
+    call warning('Lbar is defined only for q = 0. For q /= 0, it is always set to "full".')
+    BSE_L_kind="default-full"
+   endif
+ endif
  !
  if (l_col_cut.and.what=="init") then
    ! We are here in 0D, 1D or 2D
@@ -78,18 +80,6 @@ subroutine K_driver_init(what,iq,Ken,Xk)
      call warning('Coulomb cutoff at q=0. Default/Suggested Lkind= bar')
      if (STRING_match(BSE_L_kind,"default")) BSE_L_kind="default-bar"
    endif
-   if (iq/=1) then
-     ! We are here in 1D or 2D.
-     ! At q>0 there is L-T spliting in 1D or 2D
-     ! Lbar and Lfull give the same result only in the limit V->infinity
-     ! The above mentioned G=0 term is part of this effect.
-     ! Here we assume it would be better to include it, to be double checked
-     call warning('Coulomb cutoff at q>0. Default/Suggested Lkind= full')
-     if (STRING_match(BSE_L_kind,"default")) BSE_L_kind="default-full"
-     ! N.B.: One could compute the calculation without L-T slitting,
-     ! by removing such full line (2D) / plane (1D). Not coded at the moment.
-     ! I would define it as L="transverse"
-   endif
  endif
  !
  ! Is absorption (alpha or epsilon) proportional to chi
diff --git a/src/pol_function/X_dielectric_matrix.F b/src/pol_function/X_dielectric_matrix.F
index eae1511d25..238706c707 100644
--- a/src/pol_function/X_dielectric_matrix.F
+++ b/src/pol_function/X_dielectric_matrix.F
@@ -348,6 +348,8 @@ integer function X_dielectric_matrix(Xen,Xk,q,X,Xw,Dip,SILENT_MODE,CHILD)
      !
      if (.not.l_rpa_IP) call X_redux(iq,"X"//trim(X_what),X_par(iq_mem),Xw,X)
      !
+     if (iq==1) call X_symmetrize_q0(X_par(iq_mem),X%ng,Xw)
+     !
      if (n_OPTICAL_dir_to_eval>1.and.iq_dir==1) call MATRIX_copy(X_par(iq_mem),X_par_average,.TRUE.)
      if (n_OPTICAL_dir_to_eval>1)               call X_AVERAGE_do_it("ACCUMULATE",X_par(iq_mem))
      !
@@ -486,3 +488,66 @@ subroutine elemental_IO(iq_,this_is_Xo)
    end subroutine
    !
 end function
+
+subroutine X_symmetrize_q0(X_par, NG, Xw)
+  use pars,          ONLY: SP, DP, cZERO
+  use matrix,        ONLY: PAR_matrix
+  use frequency,     ONLY: w_samp
+  use D_lattice,     ONLY: nsym, i_time_rev
+  use R_lattice,     ONLY: g_rot
+  use parallel_int,  ONLY: PP_redux_wait
+  use com,           ONLY: msg
+
+  implicit none
+
+  type(PAR_matrix), intent(inout) :: X_par
+  integer,          intent(in)    :: NG
+  type(w_samp),     intent(in)    :: Xw
+
+  integer :: ig1, ig2, ig1_sym, ig2_sym, isym, iw
+  integer :: nsym_loop, ig1_loc
+  complex(SP), allocatable :: X_full(:,:)
+  !
+  if (nsym <= 1) return
+  !
+  call msg('s', 'Applying q->0 symmetries to the macroscopic dielectric matrix')
+  !
+  allocate(X_full(NG, NG))
+  ! 
+  do iw = 1, Xw%n_freqs  
+     X_full = cZERO
+     if (X_par%nrows>0) X_full(X_par%rows(1):X_par%rows(2), 1:NG) = X_par%blc(X_par%rows(1):X_par%rows(2), 1:NG, iw)
+     call PP_redux_wait(X_full, COMM=X_par%INTER_comm%COMM)
+     if (X_par%nrows>0) X_par%blc(X_par%rows(1):X_par%rows(2), 1:NG, iw) = cZERO
+     ! Note :
+     ! Time reversal maps (X(q,...,E) to X(q,...,-E)
+     ! so we only apply to E = 0 case
+     if (abs(Xw%p(iw)) < 1.e-5_SP) then
+        nsym_loop = nsym
+     else
+        nsym_loop = nsym / (i_time_rev + 1)
+     endif
+     do ig2 = 1, NG
+        do ig1_loc = 1, X_par%nrows
+           ig1 = ig1_loc + X_par%rows(1) - 1
+           do isym = 1, nsym_loop
+              ig1_sym = g_rot(ig1, isym)
+              ig2_sym = g_rot(ig2, isym)
+              !
+              if (ig1_sym > 0 .and. ig1_sym <= NG .and. ig2_sym > 0 .and. ig2_sym <= NG) then
+                 if (isym > nsym / (i_time_rev + 1)) then
+                    X_par%blc(ig1, ig2, iw) = X_par%blc(ig1, ig2, iw) + &
+                         conjg(X_full(ig1_sym, ig2_sym)) / real(nsym_loop, SP)
+                 else
+                    X_par%blc(ig1, ig2, iw) = X_par%blc(ig1, ig2, iw) + &
+                         X_full(ig1_sym, ig2_sym) / real(nsym_loop, SP)
+                 endif
+              endif
+           enddo
+        enddo
+     enddo
+  enddo
+  !
+  deallocate(X_full)
+!
+end subroutine X_symmetrize_q0