Problem about periodic alda TDDFT

[freeyourself-all]

I suspect there may be an issue with the periodic ALDA TDDFT implementation in Yambo(5.3.0) when using multiple k-points. I'm using Yambo as a reference to compare results from the periodic Casida equation, so I extracted the BS_mat matrix from the subroutine K_stored_in_a_big_matrix, as well as the preceding H_x variable from K_kernel.F (which represents half of the matrix).
To make the computation more efficient, Yambo calculates only half of the matrix and reconstructs the other half using symmetry — where matrix A is Hermitian and matrix B is symmetric.
However, when I run a calculation for silicon with a 2×2×2 k-point grid, I find that the exchange-correlation (xc) part of the B matrix is identical to that of the A matrix. This is extracted from H_x when iHxc = 2 (with mode_now = R and C) in K_kernel.F.
According to the Casida equation, the matrix element of A and B should be (ai|f_{xc}|bj) and (ai|f_{xc}|jb) respectively. These should generally not be the same.
Notably, this issue does not occur when I use a single k-point, where the A and B matrix xc components are correctly different.
Therefore, I'm not sure whether this behavior is expected, a result of incorrect usage on my part, or a bug in the program — so I'm submitting this issue for clarification.
This is the command I used to generate the Yambo input file.
yambo -F yambo.in_ALDA -J 3D_ALDA -o b -k alda -Ksolver d
Working directories:

Si_111_yambo_cp.zip
Si_222_yambo_cp.zip

[sangallidavide]

Dear @freeyourself-all thank you for the feedback.

You are right, the fxc imlpementation was recently re-written, but an error was introduced with the coupling term defined without swapping the indexes, as you mention. I'm not sure why for a single k-point you get them different.

One simple solution is to set in input to use the old implementation (much slower, but hopefully correct). Put in your input file
FxcMode= "R-def"

Instead, to use the new algorithm, you need to patch the code. Here a solution from the lumen fork of the yambo code:
https://gitlab.com/lumen-code/lumen/-/commit/2f46ca13e0c56d9784084c884b2202a3cbf103e2

Please let me know if the patch works. I did not test it.

[freeyourself-all]

Dear @freeyourself-all thank you for the feedback.

You are right, the fxc imlpementation was recently re-written, but an error was introduced with the coupling term defined without swapping the indexes, as you mention. I'm not sure why for a single k-point you get them different.

One simple solution is to set in input to use the old implementation (much slower, but hopefully correct). Put in your input file FxcMode= "R-def"

Instead, to use the new algorithm, you need to patch the code. Here a solution from the lumen fork of the yambo code: https://gitlab.com/lumen-code/lumen/-/commit/2f46ca13e0c56d9784084c884b2202a3cbf103e2

Please let me know if the patch works. I did not test it.

Dear @sangallidavide ,
Thank you very much for your prompt reply to my question.
I have tried both solutions you provided and obtained two different results. The old implementation (FxcMode="R-def") produces results that are closer to those I calculated myself, though there is still a discrepancy of about 10-5 Hartree. This may be due to the implementation of fxc (finite differences in Yambo) or some other factors, and I will continue investigating to find the cause.
As for the result from the patch, it appears to be different from the old implementation, or perhaps I am making the same mistake as in the old implementation (FxcMode="R-def") .

[sangallidavide]

Dear @freeyourself-all

by looking again into the code, the old implementation (FxcMode="R-def") was inverting the indexes in this subroutine
https://github.com/yambo-code/yambo/blob/master/src/tddft/TDDFT_ALDA_eh_space_R_kernel.F

in the new implementation (FxcMode="G-def" which is the default), without the patch, the inversion of the indexes was done here
https://github.com/yambo-code/yambo/blob/master/src/tddft/TDDFT_ALDA_eh_space_G_collisions_R.F

In both cases check the action of mode="C" where "C" means "Coupling".
I'd expect the two approaches to be equivalent. So, not sure why you get different results.

Indeed, my previous reply was triggered by the fact that I missed the operation done in TDDFT_ALDA_eh_space_G_collisions_R.F for the coupling case).

If you use the patch, you should remove such operation e.g. comment the lines

   if (mode=="R".or.mode=="Q".or.mode=="F") then
     is(1)=i_c
     os(1)=i_v
   else if (mode=="A".or.mode=="C") then
     ! Same momenta and spin indexes, inverted bands
     is(1)=i_v
     os(1)=i_c
   endif

This would adopt the procedure also used for the exchange term in the coupling. Not sure if that (patch+ commented lines) would give differeneces w.r.t. the original version for FxcMode="G-def".

[freeyourself-all]

Dear @freeyourself-all

by looking again into the code, the old implementation (FxcMode="R-def") was inverting the indexes in this subroutine https://github.com/yambo-code/yambo/blob/master/src/tddft/TDDFT_ALDA_eh_space_R_kernel.F

in the new implementation (FxcMode="G-def" which is the default), without the patch, the inversion of the indexes was done here https://github.com/yambo-code/yambo/blob/master/src/tddft/TDDFT_ALDA_eh_space_G_collisions_R.F

In both cases check the action of mode="C" where "C" means "Coupling". I'd expect the two approaches to be equivalent. So, not sure why you get different results.

Indeed, my previous reply was triggered by the fact that I missed the operation done in TDDFT_ALDA_eh_space_G_collisions_R.F for the coupling case).

If you use the patch, you should remove such operation e.g. comment the lines

   if (mode=="R".or.mode=="Q".or.mode=="F") then
     is(1)=i_c
     os(1)=i_v
   else if (mode=="A".or.mode=="C") then
     ! Same momenta and spin indexes, inverted bands
     is(1)=i_v
     os(1)=i_c
   endif

This would adopt the procedure also used for the exchange term in the coupling. Not sure if that (patch+ commented lines) would give differeneces w.r.t. the original version for FxcMode="G-def".

Dear @sangallidavide
Thank you for your detailed and patient reply. I have also looked into the code and confirmed that Yambo does swap the indexes as you described https://github.com/yambo-code/yambo/blob/master/src/tddft/TDDFT_ALDA_eh_space_G_collisions_R.F.
So the difference in results between the new implementation (FxcMode="G-def" which is the default) and the old implementation (FxcMode="R-def") is strange. To investigate further, I printed more intermediate variables and found that the variable mode in https://github.com/yambo-code/yambo/blob/master/src/tddft/TDDFT_ALDA_eh_space_G_collisions_R.F always set to "R", never to "C".
So I suspect that this condition in https://github.com/yambo-code/yambo/blob/master/src/bse/K_kernel.F does not reallocate memory when calculating the coupling.

     if (.not.allocated( BS_T_grp(i_Tgrp_k)%O_tddft_L) ) then
       call TDDFT_oscillators_alloc_L(i_Tgrp_k)
       call TDDFT_ALDA_eh_space_G_collisions_L(iq,Xk,i_Tgrp_k,NG,l_bs_tddft_wf_in_loop,tddft_wf,mode_now)
     endif
     if (.not.allocated( BS_T_grp(i_Tgrp_p)%O_tddft_R) ) then
       call TDDFT_oscillators_alloc_R(i_Tgrp_p)
       call TDDFT_ALDA_eh_space_G_collisions_R(iq,Xk,i_Tgrp_p,NG,l_bs_tddft_wf_in_loop,tddft_wf,mode_now)
     endif
   endif

To address this, I added an extra check to explicitly deallocate the memory before allocating it again. With this change, I found that the results from the new implementation (FxcMode="G-def" which is the default) and the old implementation (FxcMode="R-def") are now consistent. Here is the modified code:

if(iHxc==2 .and. l_tddft_gsum) then
     if (.not.allocated( BS_T_grp(i_Tgrp_k)%O_tddft_L) ) then
       call TDDFT_oscillators_alloc_L(i_Tgrp_k)
       call TDDFT_ALDA_eh_space_G_collisions_L(iq,Xk,i_Tgrp_k,NG,l_bs_tddft_wf_in_loop,tddft_wf,mode_now)
     else
       deallocate( BS_T_grp(i_Tgrp_k)%O_tddft_L)
       call TDDFT_oscillators_alloc_L(i_Tgrp_k)
       call TDDFT_ALDA_eh_space_G_collisions_L(iq,Xk,i_Tgrp_k,NG,l_bs_tddft_wf_in_loop,tddft_wf,mode_now)
     endif
     if (.not.allocated( BS_T_grp(i_Tgrp_p)%O_tddft_R) ) then
       call TDDFT_oscillators_alloc_R(i_Tgrp_p)
       call TDDFT_ALDA_eh_space_G_collisions_R(iq,Xk,i_Tgrp_p,NG,l_bs_tddft_wf_in_loop,tddft_wf,mode_now)
     else
       deallocate( BS_T_grp(i_Tgrp_p)%O_tddft_R)
       call TDDFT_oscillators_alloc_R(i_Tgrp_p)
       call TDDFT_ALDA_eh_space_G_collisions_R(iq,Xk,i_Tgrp_p,NG,l_bs_tddft_wf_in_loop,tddft_wf,mode_now)
     endif
   endif

[sangallidavide]

Dear @freeyourself-all you are right.

This also explains the discrepancy between the case of an isolated system (0D) and the case if a periodic system.

Indeed the allocation and de-allocation is done differently in the two cases. See the definition of l_BSE_minimize_memory here:
https://github.com/yambo-code/yambo/blob/master/src/bse/K_driver_init.F#L182

So, instead of your code, you can just set in input PAR_def_mode="memory" (it will act as your patch).
I'll release a fix for the case where the oscillators are kept in memory.

[sangallidavide]

Dear @freeyourself-all I released a general fix. No need to set PAR_def_mode="memory"

Check this branch:
https://gitlab.com/lumen-code/lumen/-/tree/fix-issue-215

Note 1: the new algorithm implemented in the lumen fork of the yambo code was almost correct. There was just a missing conjg(). See last line of function TDDFT_ALDA_eh_space_G_kernel_coupling inside TDDFT_ALDA_eh_space_G_kernel.F
However such code only works for q=0

Note 2: as you noticed the inversion of the indexes inside TDDFT_ALDA_eh_space_G_collisions_R.F / TDDFT_ALDA_eh_space_G_collisions_L.F was never activated. It has been removed.

Note 3: if you want to use inverted indexes at q=0 you can set in input ImposeAsymBSE. This activate the construction of the kernel with an explicit loop over both cv and vc indexes. Otherwise it loops over cv and uses the operation already mentioned in the function TDDFT_ALDA_eh_space_G_kernel_coupling to get the vc kernel from the cv oscilators
At finite q, instead, it always loops over both cv and vc. See also comments around line 475