Skip to content

Symmetrize X0(q->0) #228

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
muralidhar-nalabothula wants to merge 8 commits into
base: 5.4
Choose a base branch
from
Open

Symmetrize X0(q->0) #228

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
4f3d2a3
symmetrize eps_inv(q->0)
muralidhar-nalabothula May 11, 2026
b0326b0
fix lbar issue for finite q and default
muralidhar-nalabothula May 11, 2026
7eb945e
reduce memory
muralidhar-nalabothula May 11, 2026
d2dfc63
set it to default-full
muralidhar-nalabothula May 12, 2026
c075626
set it to default-full
muralidhar-nalabothula May 12, 2026
7843442
write to log
muralidhar-nalabothula May 17, 2026
9253c27
fix bug when using g parallization, remove debug code
muralidhar-nalabothula May 17, 2026
dd507a6
cleanup
muralidhar-nalabothula May 17, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
22 changes: 6 additions & 16 deletions src/bse/K_driver_init.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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"
Copy link
Copy Markdown
Member

@sangallidavide sangallidavide May 12, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The code would do, by default, Lbar at finite q in 1D and 2D. Is this a good default?

Copy link
Copy Markdown
Contributor Author

@muralidhar-nalabothula muralidhar-nalabothula May 12, 2026 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I did not change any defaults in this PR. The code was simply too strict when I ran the BSE calculation without specifying the Lkind variable. I removed only a redundant block of code. I am also very happy that we were able to get rid of bar for finite q.

Copy link
Copy Markdown
Member

@sangallidavide sangallidavide May 12, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

From what I see. The lines above (kept from the original version)

 if (l_col_cut.and.what=="init") then
   ...
   if (STRING_match(BSE_L_kind,"default")) BSE_L_kind="default-bar"
   ...

are imposing Lbare by default for l_col_cut case.

Then I see this PR removes the subsequent lines

     if (iq/=1) then
       ...
       if (STRING_match(BSE_L_kind,"default")) BSE_L_kind="default-full"
       ...
     endif

Accordingly, with coulomb cutoff, i.e. in D<3 L-bar is used at any q.

Copy link
Copy Markdown
Member

@sangallidavide sangallidavide May 12, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Note. Check this in the "files changed section". It should be easier to read it.

Copy link
Copy Markdown
Contributor Author

@muralidhar-nalabothula muralidhar-nalabothula May 12, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There is this block right above

 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")) 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

so this part is not needed

    if (iq/=1) then
      ...
      if (STRING_match(BSE_L_kind,"default")) BSE_L_kind="default-full"
      ...
    endif

Copy link
Copy Markdown
Member

@sangallidavide sangallidavide May 12, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah, but then also this part needs to be removed, otherwise it overwirtes the above block (?)

 if (l_col_cut.and.what=="init") then
   ...
   if (STRING_match(BSE_L_kind,"default")) BSE_L_kind="default-bar"
   ...
 endif

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
Expand Down
65 changes: 65 additions & 0 deletions src/pol_function/X_dielectric_matrix.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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))
!
Expand Down Expand Up @@ -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