[Q-e-developers] minor bug in new preconditioner of g_psi.f90

[Brian Austin]

Dear QE developers,

I believe I've discovered a minor bug in the released version of QE-5.0.
The problem is in the TEST_NEW_PRECONDITIONING block of PW/src/g_psi.f90.
Note that TEST_NEW_PRECONDITIONING is defined at the top of the file,
so this is an active block of code.

! g_psi.f90:42
#ifdef TEST_NEW_PRECONDITIONING
  scala = 1.d0
  do ipol=1,npol
     do k = 1, m
        do i = 1, n
           x = (h_diag(i,ipol) - e(k)*s_diag(i,ipol))*scala
           denm = (1.d0+x+sqrt(1.d0+(x-1)*(x-1.d0)))/scala
        !         denm = 1.d0 + 16*x*x*x*x/(27.d0+18*x+12*x*x+8*x*x*x)
           psi (i, ipol, k) = psi (i, ipol, k) / denm
        enddo
     enddo
  enddo
#else
! snip !
#endif

For large x, denm approaches 2*x, but I think it should approach 1*x,
so that g_psi approximates ((H-eS)^{-1}) psi
I suggest the following change:
           x = (h_diag(i,ipol) - e(k)*s_diag(i,ipol))
           denm = 0.5_dp * (1.d0+x+sqrt(1.d0+(x-1)*(x-1.d0)))


Fortunately, this bug is inconsequential because the extra factor of
1/2 gets propagated to both the Hamiltonian and overlap matrices, so
it will be accounted for during the diagonalization process, and
eventually corrected by refresh_evc().

Thanks for reading, and for your efforts with QE.
-Brian