SIRIUS LR results discrepancies for non-zero q

[gsavva]

Using bulk Ni (ferromagnetic metal) from HP examples as use case, with NC-PP (Ni_inputs.zip).

The scf step (pw.x) is run by native QE 7.2.

For q-grid = 1x1x1 (nq1 = 1, nq2 = 1, nq3 = 1), the calculated U value from native QE 7.2 and q-e-sirius, agree perfectly:

                                 Hubbard U parameters:

       site n.  type  label  spin  new_type  new_label  manifold  Hubbard U (eV)
         1        1    Ni      1      1         Ni         3d       2.7925

For non-zero q, the "RESPONSE OCCUPATION MATRICES" (enabled by iverbosity = 4), even in the 1st iteration differ for spin 2

native QE

      atom #  1   q point #   2   iter #   1               0.180E-12
      RESPONSE OCCUPATION MATRICES:
      Hubbard atom  1  spin  1
     -0.0140399990    0.0002194360   -0.0008829688    0.0229230653    0.0001443888
      0.0002194360   -0.0740814893   -0.0459746036   -0.0009443695    0.0460902817
     -0.0008829688   -0.0459746036   -0.1069103422    0.0015041946    0.1041923848
      0.0229230653   -0.0009443695    0.0015041946   -0.0403943551   -0.0006496459
      0.0001443888    0.0460902817    0.1041923848   -0.0006496459   -0.1047424709
      Hubbard atom  1  spin  2
     -7.9025880489   -2.8511248445   -0.1191074931   -4.3180679758    4.0357382722
     -2.8511248445   -7.8404436305    1.4706187977   -2.7743955636    3.0094071347
     -0.1191074931    1.4706187977   -0.1780607553    0.1578916184    1.4496396617
     -4.3180679758   -2.7743955636    0.1578916184   -3.0448601351    1.1061530681
      4.0357382722    3.0094071347    1.4496396617    1.1061530681   -7.6327640468

q-e-sirius

      atom #  1   q point #   2   iter #   1               0.180E-12
      RESPONSE OCCUPATION MATRICES:
      Hubbard atom  1  spin  1
     -0.0140399791    0.0002194353   -0.0008829670    0.0229230321    0.0001443881
      0.0002194353   -0.0740813124   -0.0459744642   -0.0009443677    0.0460901421
     -0.0008829670   -0.0459744642   -0.1069100088    0.0015041915    0.1041920539
      0.0229230321   -0.0009443677    0.0015041915   -0.0403942969   -0.0006496444
      0.0001443881    0.0460901421    0.1041920539   -0.0006496444   -0.1047421394
      Hubbard atom  1  spin  2
     -2.1820462688   -2.7824075388   -3.3309471976   -7.4166891395   -0.5254798096
     -2.7824075388   -8.5839242230    1.4589533516   -1.9899165453    2.9814711806
     -3.3309471976    1.4589533516    1.6888641235    1.9193620090    3.8262795116
     -7.4166891395   -1.9899165453    1.9193620090   -2.2566624204    3.5171287057
     -0.5254798096    2.9814711806    3.8262795116    3.5171287057   -4.0354559194

absolute difference

      RESPONSE OCCUPATION MATRICES:
      Hubbard atom  1  spin  1
          1.99E-08      7.00E-10        1.80E-09       3.32E-08       7.00E-10
          7.00E-10      1.77E-07        1.39E-07       1.80E-09       1.40E-07
          1.80E-09      1.39E-07        3.33E-07       3.10E-09       3.31E-07
          3.32E-08      1.80E-09        3.10E-09       5.82E-08       1.50E-09
          7.00E-10      1.40E-07        3.31E-07       1.50E-09       3.31E-07
      Hubbard atom  1  spin  2
      5.7205417801   0.0687173057   3.2118397045   3.0986211637   4.5612180818
      0.0687173057   0.7434805925   0.0116654461   0.7844790183   0.0279359541
      3.2118397045   0.0116654461   1.8669248788   1.7614703906   2.3766398499
      3.0986211637   0.7844790183   1.7614703906   0.7881977147   2.4109756376
      4.5612180818   0.0279359541   2.3766398499   2.4109756376   3.5973081274

The "RESPONSE OCCUPATION MATRICES" are calculated by hp_dnsq, based on dspi, the output of the sirius_linear_solver

[gsavva]

A simpler example that reproduces the above issue: CTi, a non-magnetic metallic system (CTi_inputs.zip)

Using the following:

• Norm-conserving pseudopotentials (NC PP)
• U and V values in the scf step = 1e-08
• For the Response run via hp.x:

  iverbosity = 4
  nq1 = 2
  nq2 = 2
  nq3 = 2
  perturb_only_atom(1) = .true.  ! perturb only Ti
  start_q = 2
  last_q = 2

• very tight threshold for chi: conv_thr_chi = 1.0000000000d-12

For q!=0 (hence start_q = 2, last_q = 2), the "RESPONSE OCCUPATION MATRICES", even in the 1st iteration, differ for Hubbard atom 1 spin 1

[gsavva]

After investigation with @gcistaro, and inspecting ch_psi_all_k from QE vs. Linear_response_operator::multiply() from SIRIUS, we found out the following.

In QE, the index for the number of occupied bands is ikqs (i.e. k+q): code line.
On the other hand, SIRIUS, only knows about the number of occupied bands at k: code line

Even more strangely, the native QE cgsolve_all (ref) gets nbnd_occ(ikk) as input argument, but uses instead nbnd_occ (ikqs(ik)) from global scope ! (for which SIRIUS has no information about)

[gsavva]

Commits by @gcistaro eee2fa4 on the SIRIUS side and electronic-structure/q-e-sirius@fa45193 on the q-e-sirius side should fix the problem. Testing underway...

[gsavva]

issue fixed via PRs:

• SIRIUS: Fix linear solver for metals #997
• q-e-sirius:Fix linear solver for metals q-e-sirius#59