From 4f3d2a30ef27f1a8906d36ee160088a3dd372855 Mon Sep 17 00:00:00 2001 From: Murali Date: Mon, 11 May 2026 22:18:53 +0200 Subject: [PATCH 1/8] symmetrize eps_inv(q->0) --- src/pol_function/X_dielectric_matrix.F | 66 ++++++++++++++++++++++++++ 1 file changed, 66 insertions(+) diff --git a/src/pol_function/X_dielectric_matrix.F b/src/pol_function/X_dielectric_matrix.F index eae1511d25..c9b571ec54 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,67 @@ 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('r', 'Applying q->0 symmetries to the macroscopic dielectric matrix') + + allocate(X_full(NG, NG, 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(1:X_par%nrows, 1:NG, :) + + call PP_redux_wait(X_full, COMM=X_par%INTER_comm%COMM) + + X_par%blc = cZERO + + do iw = 1, Xw%n_freqs + 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_loc, ig2, iw) = X_par%blc(ig1_loc, ig2, iw) + & + conjg(X_full(ig1_sym, ig2_sym, iw)) / real(nsym_loop, SP) + else + X_par%blc(ig1_loc, ig2, iw) = X_par%blc(ig1_loc, ig2, iw) + & + X_full(ig1_sym, ig2_sym, iw) / real(nsym_loop, SP) + endif + endif + enddo + enddo + enddo + enddo + + deallocate(X_full) + +end subroutine X_symmetrize_q0 From b0326b0a430bff7b6edf122f72c5d88964036320 Mon Sep 17 00:00:00 2001 From: Murali Date: Mon, 11 May 2026 22:33:41 +0200 Subject: [PATCH 2/8] fix lbar issue for finite q and default --- src/bse/K_driver_init.F | 21 +++++---------------- 1 file changed, 5 insertions(+), 16 deletions(-) diff --git a/src/bse/K_driver_init.F b/src/bse/K_driver_init.F index 1750b595fe..45e166e8bf 100644 --- a/src/bse/K_driver_init.F +++ b/src/bse/K_driver_init.F @@ -59,11 +59,12 @@ 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")) & +& call warning('Lbar is defined only for q = 0. For q /= 0, it is always set to "full".') + BSE_L_kind="full" + endif ! if (l_col_cut.and.what=="init") then ! We are here in 0D, 1D or 2D @@ -78,18 +79,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 From 7eb945ec83b99b0d87a0a9eb9a9433ab0e8476c5 Mon Sep 17 00:00:00 2001 From: Murali Date: Mon, 11 May 2026 23:12:48 +0200 Subject: [PATCH 3/8] reduce memory --- src/pol_function/X_dielectric_matrix.F | 39 +++++++++++++------------- 1 file changed, 19 insertions(+), 20 deletions(-) diff --git a/src/pol_function/X_dielectric_matrix.F b/src/pol_function/X_dielectric_matrix.F index c9b571ec54..1487e06274 100644 --- a/src/pol_function/X_dielectric_matrix.F +++ b/src/pol_function/X_dielectric_matrix.F @@ -506,49 +506,48 @@ subroutine X_symmetrize_q0(X_par, NG, Xw) integer :: ig1, ig2, ig1_sym, ig2_sym, isym, iw integer :: nsym_loop, ig1_loc - complex(SP), allocatable :: X_full(:,:,:) - + complex(SP), allocatable :: X_full(:,:) + ! if (nsym <= 1) return - + ! call msg('r', 'Applying q->0 symmetries to the macroscopic dielectric matrix') - - allocate(X_full(NG, NG, 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(1:X_par%nrows, 1:NG, :) - - call PP_redux_wait(X_full, COMM=X_par%INTER_comm%COMM) - - X_par%blc = cZERO - - do iw = 1, Xw%n_freqs + ! + 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(1:X_par%nrows, 1:NG, iw) + call PP_redux_wait(X_full, COMM=X_par%INTER_comm%COMM) + if (X_par%nrows>0) X_par%blc(1:X_par%nrows, 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_loc, ig2, iw) = X_par%blc(ig1_loc, ig2, iw) + & - conjg(X_full(ig1_sym, ig2_sym, iw)) / real(nsym_loop, SP) + conjg(X_full(ig1_sym, ig2_sym)) / real(nsym_loop, SP) else X_par%blc(ig1_loc, ig2, iw) = X_par%blc(ig1_loc, ig2, iw) + & - X_full(ig1_sym, ig2_sym, iw) / real(nsym_loop, SP) + X_full(ig1_sym, ig2_sym) / real(nsym_loop, SP) endif endif enddo enddo enddo enddo - + ! deallocate(X_full) - +! end subroutine X_symmetrize_q0 From d2dfc63ed6b789baf60158ac1a2c1c2e01ca0a72 Mon Sep 17 00:00:00 2001 From: Murali Date: Tue, 12 May 2026 12:00:38 +0200 Subject: [PATCH 4/8] set it to default-full --- src/bse/K_driver_init.F | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/bse/K_driver_init.F b/src/bse/K_driver_init.F index 45e166e8bf..01964cce84 100644 --- a/src/bse/K_driver_init.F +++ b/src/bse/K_driver_init.F @@ -63,7 +63,7 @@ subroutine K_driver_init(what,iq,Ken,Xk) ! BSE Hamiltonian is fully symmetric with Lfull if (.not. STRING_match(BSE_L_kind, "full")) & & call warning('Lbar is defined only for q = 0. For q /= 0, it is always set to "full".') - BSE_L_kind="full" + BSE_L_kind="default-full" endif ! if (l_col_cut.and.what=="init") then From c07562671d2fd5e81272c16f8be845817f6b8f75 Mon Sep 17 00:00:00 2001 From: Murali Date: Tue, 12 May 2026 13:15:19 +0200 Subject: [PATCH 5/8] set it to default-full --- src/bse/K_driver_init.F | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/bse/K_driver_init.F b/src/bse/K_driver_init.F index 01964cce84..a0ecb75a7e 100644 --- a/src/bse/K_driver_init.F +++ b/src/bse/K_driver_init.F @@ -61,9 +61,10 @@ subroutine K_driver_init(what,iq,Ken,Xk) if (iq/=1) then ! For finite momentum Lbar is NOT defined and gives WRONG results ! BSE Hamiltonian is fully symmetric with Lfull - if (.not. STRING_match(BSE_L_kind, "full")) & -& call warning('Lbar is defined only for q = 0. For q /= 0, it is always set to "full".') - BSE_L_kind="default-full" + 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 From 7843442c0f50bda020aa01caa1f11538205ebb17 Mon Sep 17 00:00:00 2001 From: Murali Date: Sun, 17 May 2026 11:04:34 +0200 Subject: [PATCH 6/8] write to log --- src/pol_function/X_dielectric_matrix.F | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/src/pol_function/X_dielectric_matrix.F b/src/pol_function/X_dielectric_matrix.F index 1487e06274..1e61b274f4 100644 --- a/src/pol_function/X_dielectric_matrix.F +++ b/src/pol_function/X_dielectric_matrix.F @@ -497,6 +497,7 @@ subroutine X_symmetrize_q0(X_par, NG, Xw) use R_lattice, ONLY: g_rot use parallel_int, ONLY: PP_redux_wait use com, ONLY: msg + use stderr, ONLY: write_to_log implicit none @@ -507,18 +508,22 @@ subroutine X_symmetrize_q0(X_par, NG, Xw) integer :: ig1, ig2, ig1_sym, ig2_sym, isym, iw integer :: nsym_loop, ig1_loc complex(SP), allocatable :: X_full(:,:) + logical :: saved_log ! if (nsym <= 1) return ! - call msg('r', 'Applying q->0 symmetries to the macroscopic dielectric matrix') + saved_log = write_to_log + if (X_par%INTER_comm%CPU_id == 0) write_to_log = .TRUE. + call msg('s', 'Applying q->0 symmetries to the macroscopic dielectric matrix') + write_to_log = saved_log ! 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(1:X_par%nrows, 1:NG, iw) + 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(1:X_par%nrows, 1:NG, iw) = cZERO + 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 @@ -536,10 +541,10 @@ subroutine X_symmetrize_q0(X_par, NG, Xw) ! 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_loc, ig2, iw) = X_par%blc(ig1_loc, ig2, iw) + & + 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_loc, ig2, iw) = X_par%blc(ig1_loc, ig2, iw) + & + X_par%blc(ig1, ig2, iw) = X_par%blc(ig1, ig2, iw) + & X_full(ig1_sym, ig2_sym) / real(nsym_loop, SP) endif endif From 9253c272c01090cd1ef58b3e04c4957a1a619858 Mon Sep 17 00:00:00 2001 From: Murali Date: Sun, 17 May 2026 12:48:09 +0200 Subject: [PATCH 7/8] fix bug when using g parallization, remove debug code --- src/pol_function/X_dielectric_matrix.F | 4 ---- 1 file changed, 4 deletions(-) diff --git a/src/pol_function/X_dielectric_matrix.F b/src/pol_function/X_dielectric_matrix.F index 1e61b274f4..888a5d3768 100644 --- a/src/pol_function/X_dielectric_matrix.F +++ b/src/pol_function/X_dielectric_matrix.F @@ -497,7 +497,6 @@ subroutine X_symmetrize_q0(X_par, NG, Xw) use R_lattice, ONLY: g_rot use parallel_int, ONLY: PP_redux_wait use com, ONLY: msg - use stderr, ONLY: write_to_log implicit none @@ -512,10 +511,7 @@ subroutine X_symmetrize_q0(X_par, NG, Xw) ! if (nsym <= 1) return ! - saved_log = write_to_log - if (X_par%INTER_comm%CPU_id == 0) write_to_log = .TRUE. call msg('s', 'Applying q->0 symmetries to the macroscopic dielectric matrix') - write_to_log = saved_log ! allocate(X_full(NG, NG)) ! From dd507a6e5b75aa9872dcd74ff8ddfc1c0673cca5 Mon Sep 17 00:00:00 2001 From: Murali Date: Sun, 17 May 2026 12:50:22 +0200 Subject: [PATCH 8/8] cleanup --- src/pol_function/X_dielectric_matrix.F | 1 - 1 file changed, 1 deletion(-) diff --git a/src/pol_function/X_dielectric_matrix.F b/src/pol_function/X_dielectric_matrix.F index 888a5d3768..238706c707 100644 --- a/src/pol_function/X_dielectric_matrix.F +++ b/src/pol_function/X_dielectric_matrix.F @@ -507,7 +507,6 @@ subroutine X_symmetrize_q0(X_par, NG, Xw) integer :: ig1, ig2, ig1_sym, ig2_sym, isym, iw integer :: nsym_loop, ig1_loc complex(SP), allocatable :: X_full(:,:) - logical :: saved_log ! if (nsym <= 1) return !