Symmetrize X0(q->0)

[muralidhar-nalabothula]

The code gives incorrect results when performing BSE calculations on systems with only a single Gamma point (e.g. molecules). The issue arises because, when computing $W(q \to 0)$, we choose a specific direction, which is no longer symmetric with respect to the point group. Therefore, this contribution should be symmetrized. The proper fix would be to perform a RIM treatment for W over several q-points, but that is considerably more demanding.

In this PR, I implement an alternative workaround for cases where the user enables symmetries (this should be a robust fix as long as symmetries are used). The idea is to rotate $\epsilon$ over all symmetry operations and average the result. In this approach, we explicitly avoid computing $\chi_0$ for more than a single $q \to 0$ point.

In addition, I slightly adjusted the Lkind defaults. With the current implementation, the fix can break down easily if Lbar is not specified correctly, and I think the existing way is too aggressive.

I also have one question:

1. I am not sure whether I am using the correct communicator to gather $X(G,G')$ onto a single CPU. Currently I am doing (it works, but please cross check):

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

Is this the correct approach, or should I instead use a different communicator / reduction routine here?

Tested on SiH4
old BSE energies at Q=0

    0.08776613       8.33745301e-03          1
      0.11589517       1.20850524e-03          2
      0.15104856       3.88729223e-03          3
      0.32560593       4.40265499e-02          4
      0.35400522       3.05387867e-03          5
      0.36274657       1.34110535e-02          6
      0.43202174       4.83342893e-02          7
      0.46968186       6.08044565e-02          8
      0.51835233       2.47406457e-02          9
      0.55277359       2.40801051e-01         10
      0.58717853       3.59762572e-02         11
      0.61871099       3.44642907e-01         12
      0.64019614       1.67792067e-01         13
      0.66484177       3.51507291e-02         14
      0.68922377       2.67852675e-02         15
      0.71332335       5.36881149e-01         16
      0.73797375       3.09134752e-01         17
      0.77709079       6.68502450e-02         18
      1.42262530       1.62761547e-02         19
      1.42498660       1.05558652e-02         20

new

      0.15306734       4.93323887e-06          1
      0.15328379       1.20675338e-07          2
      0.15377326       7.34474859e-04          3
      0.39103109       1.86077710e-02          4
      0.39110294       3.34145315e-03          5
      0.39134952       2.39068162e-04          6
      0.54828960       5.36223115e-06          7
      0.57273978       9.17885998e-08          8
      0.57323676       1.04717208e-06          9
      0.57616794       5.25700525e-05         10
      0.57637793       1.10940891e-05         11
      0.57665282       9.54457282e-05         12
      0.63472039       1.02440258e-02         13
      0.63489306       3.45902830e-01         14
      0.63537771       1.50725543e-01         15
      0.69839782       6.42127618e-02         16
      0.69868815       4.33492869e-01         17
      0.69879669       5.25059029e-02         18
      1.44143164       1.41905184e-06         19
      1.44206464       1.60859628e-07         20

the first three excitons must be degenerate

[sangallidavide]

Extra note. With Riccardo we went on discussing the average procedure also in other systems.

The improved implementation does two extra things:
a) for systems without spatial inversion (with SOC) the average is done over 6 points, and not just 3, since +dir and -dir are not equivalent in such case
See here: https://gitlab.com/lumen-code/lumen/-/work_items/277 ;
Implementation here: https://gitlab.com/lumen-code/lumen/-/blob/bug-fixes/src/pol_function/X_AVERAGE_setup.F?ref_type=heads#L58

b) instead of averaging over the Cartesian direction, the average is now done over the reciprocal lattice directions (1 0 0), (0 1 0) and (0 0 1) in reciprocal lattice units.
See here: https://gitlab.com/lumen-code/lumen/-/blob/bug-fixes/src/pol_function/X_AVERAGE_setup.F?ref_type=heads#L83
This should respect symmetries

Not sure if this helps with your issue, or if the alternative average you propose would be needed anyway. Now I'm in a call. As soon as I have some time, I'll think about this as well. Meanwhile, please let me know you opinion on points a) and b) above.

[sangallidavide]

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

[muralidhar-nalabothula]

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.

[sangallidavide]

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.

[sangallidavide]

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

[muralidhar-nalabothula]

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

[sangallidavide]

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

[muralidhar-nalabothula]

Hi Davide,

Thanks for the comments :)

1. First, I would strongly prefer to avoid computing the response function 6 times. Even performing a single $\chi_0$ calculation is already quite expensive for my systems.

2. For a 2D material, averaging over $\pm z$ is not really correct, since the concept of $q_z$ does not exist when there is no periodicity along that direction.

3. The star of $q$ consists of symmetry-equivalent $q$-points. Unless we explicitly probe a particular direction, we should not distinguish between them. If we arbitrarily choose one point from the star of $q$ and compute $W(q)$, we effectively end up in a broken-symmetry situation (which is what happens currently). In the present implementation we are not intentionally selecting any preferred direction, so in order to preserve the symmetry or equivalently, not distinguish one q-point from the others — we should consider all symmetry-related directions and average over them.

4. What you suggest is somewhat different: you pick 3/6 inequivalent points and average over them, but those points may themselves not be related by symmetry. Also, there is no guarantee that applying all symmetry operations to $W(q \to 0)$ leaves it invariant. In some cases the errors may partially cancel if we are lucky, but that is not something we can generally rely on. That said, your approach can still work if the user does not use symmetries, which is also closer in spirit to the RIM-like approach I mentioned earlier.

[sangallidavide]

1. First, I would strongly prefer to avoid computing the response function 6 times. Even performing a single χ 0 calculation is already quite expensive for my systems.

Clear. If a better option is available, it is welcome.

2. For a 2D material, averaging over ± z is not really correct, since the concept of q z does not exist when there is no periodicity along that direction.

In 3D the average is done with (without inversion) over 3 (6) directions, in 2D over 2 (4), in 1D over 1 (2).
In 0D I'm not sure why we should average. Shouldn't the head and the wings (only terms which are direction dependent) go to zero say with a spherical coulomb cut-off ?

3. The star of q consists of symmetry-equivalent q -points. Unless we explicitly probe a particular direction, we should not distinguish between them. If we arbitrarily choose one point from the star of q and compute W ( q ) , we effectively end up in a broken-symmetry situation (which is what happens currently). In the present implementation we are not intentionally selecting any preferred direction, so in order to preserve the symmetry or equivalently, not distinguish one q-point from the others — we should consider all symmetry-related directions and average over them.

4. What you suggest is somewhat different: you pick 3/6 inequivalent points and average over them, but those points may themselves not be related by symmetry. Also, there is no guarantee that applying all symmetry operations to W ( q → 0 ) leaves it invariant. In some cases the errors may partially cancel if we are lucky, but that is not something we can generally rely on. That said, your approach can still work if the user does not use symmetries, which is also closer in spirit to the RIM-like approach I mentioned earlier.

Yes, the phylosophy is similar to the RIM (random integration methon), but I'd rather call it SIM (Selected-points integration method. :-)

Thinking loud here. The starting idea (averaging over xyz) was thought for system with cubic/rectangular unit cell. The average over the rlu direction was meant as an improvement, thinking to hexagonal cells in 2D as an example. Re-thinking the hexagonal case, one would instead need to average over 3 (6) directions (1 0 0), (0 1 0) (-1 1 0) (plus same with -1). To define the number we could count the number of G vectors with the same modulus as (1 0 0), (0 1 0) in 2D, also considering (0 0 1) in 3D. That might imply something even larger than 6.

[daniele-varsano]

Hi Murali,
the test you are showing is done with bare coulomb 1/q^2 right? If a truncated coulomb cutoff is used, I agree with Davide that head and wings of Wc(q=0) should be zero in 0D (numerically order of 10^-10 and 10^-5), while the other components should not depend on the choice of the direction. Am I missing something?

[muralidhar-nalabothula]

Hi Daniele,

I am using the cutoff, and the results do not change with or without it, which is actually quite surprising. Maybe I am doing something wrong. I agree that when using a spherical cutoff, $W(q=0)$ should become direction independent.

Nevertheless, this issue affects all single-point calculations and is not limited to molecules.

Here is the input files (I changed few values like scisor and can be run on a laptop)
SiH4.zip

here is the updated plot with the input files given

[daniele-varsano]

Hi Murali,
I just had a look to your input. The problem is that the keyword of the runlevel "rim_cut" is missing and in this way the cutoff potential is ignored, and a bare potential is considered. Other two minor comments:

1. CUTGeo= "sphere xyz". # set the directions xyz. Even if in principle they are not needed, if I'm not wrong Yambo expects it.
2. CUTRadius= 8.0. # Should be ok, but you can set to 10 (half of the cell size). In this way you are sure to be safe.

These two issues actually could be hard-coded to make the user's life easier.

[muralidhar-nalabothula]

Hi Daniele,

My bad, rim_cut fixed it, and now both approaches give identical results.

Nevertheless, I still think we need to symmetrize it for 2D or 3D systems whenever we have to choose a direction for q \to ->0 (particularly when using a single k-point or only a few k-points).

[sangallidavide]

Ok, nice that we confirm no average is needed in 0D. Also my bad I suggested to try it in the past for the 0D case.

I agree, the new implementation is useful for 2D/3D (even 1D in some cases)

[daniele-varsano]

Hi,
Happy to know it worked. I agree that for 2d/3d it is useful and thanks for pointing this out.