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

[Andrea Dal Corso]

I think that the phonon with constraint magnetization has never been
implemented. If your modification makes the two fermi levels case work,
please go on an modify the code.

Andrea

On Fri, 2017-10-13 at 17:04 +0000, Timrov Iurii wrote:
> 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
> 
> 
> _______________________________________________
> Q-e-developers mailing list
> Q-e-developers at qe-forge.org
> http://qe-forge.org/mailman/listinfo/q-e-developers

-- 
Andrea Dal Corso                    Tel. 0039-040-3787428
SISSA, Via Bonomea 265              Fax. 0039-040-3787249
I-34136 Trieste (Italy)             e-mail: dalcorso at sissa.it