[Q-e-developers] [SUSPECT ATTACHMENT REMOVED] Bug in PHonon when nspin=2 and tot_magnetization is used

[Timrov Iurii]

Dear developers,


It seems there is a bug in the PHonon code when nspin=2 and tot_magnetization is used in PWscf. The problem is in the routine LR_Modules/setup_nbnd_occ.f90 where the code computes the number of occupied bands nbnd_occ for every k point. So, when tot_magnetization is used in PWscf, two_fermi_energies=.TRUE., and nbnd_occ is zero for k+q points. It is so because the weights for k+q points are set to zero:


LR_Modules/setup_nscf.f90 -> PW/src/set_kplusq.f90


and as a consequence wg_up and wg_dw are also zero for k+q points (because the routine iweights is used), and this in turn leads to the fact that nbnd_occ is also zero for k+q points. The thing is that nbnd_occ(k+q) is used in several places in PHonon, and so the code gives wrong results:


      Representation #  1 mode #   1

     Self-consistent Calculation
     kpoint   2 ibnd  17 solve_linter: root not converged  3.846E+07
     kpoint   6 ibnd  17 solve_linter: root not converged  4.992E+07
     kpoint   8 ibnd  17 solve_linter: root not converged  1.526E+06
     kpoint  12 ibnd  17 solve_linter: root not converged  2.761E+01


The solution to this problem is simple: to replace the way nbnd_occ is computed when two_fermi_energies=.TRUE. by these lines:


           DO ik = 1, nks
              IF (isk(ik)==1) THEN
                 nbnd_occ(ik) = NINT(nelup)
              ELSE
                 nbnd_occ(ik) = NINT(neldw)
              ENDIF
           ENDDO


Please note that setting wk for k+q points different from zero [e.g. equal to the weights of k points (which is actually not true, I think)] is not a good solution, because nonzero weights for k+q points will mess up the occupations in the case of metals at the NSCF stage of the PHonon calculation.


I did tests for NiO and this fix seems to work fine. The test script is attached. Actually this bug is in QE since a very long time. I did my tests with QE 5.1. I also wanted to make tests with the latest version of QE (the SVN version), but it seems there is yet another problem in the PWscf calculation (2nd SCF calculation in the script which I attached crashes). The error message is:




Self-consistent Calculation

     iteration #  1     ecut=    50.00 Ry     beta= 0.70
     Davidson diagonalization with overlap

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     Error in routine cdiaghg (58):
     S matrix not positive definite
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

     stopping ...


Please let me know if you agree that there is a bug in the calculation of nbnd_occ. And concerning the second problem, it does not appear in QE 5.1, so I guess it may be related to some recent changes in QE.


Greetings,

Iurii



--
Iurii Timrov
Swiss Federal Institute of Technology Lausanne (EPFL)
Laboratory of Theory and Simulation of Materials (THEOS)
CH-1015 Lausanne, Switzerland


-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.quantum-espresso.org/pipermail/developers/attachments/20171013/5d9d57b2/attachment.html>