1!
2! Copyright (C) 2016 Quantum ESPRESSO Foundation
3! This file is distributed under the terms of the
4! GNU General Public License. See the file `License'
5! in the root directory of the present distribution,
6! or http://www.gnu.org/copyleft/gpl.txt .
7!
8!-----------------------------------------------------------------------
9subroutine h_prec (ik, evq, h_diag)
10  !-----------------------------------------------------------------------
11  !
12  ! ... Compute the precondition vector h_diag used in the solution of the
13  ! ... linear system - On input:
14  ! ... ik     index of k-point
15  ! ... evq    wavefunction at k+q point
16  !...  h_diag must be allocated
17  ! ... g2kin  contains kinetic energy (k+q+G)^2 for the current k+q point
18  !
19  USE kinds,      ONLY : dp
20  USE klist,      ONLY : ngk
21  USE qpoint,     ONLY : ikqs, ikks
22  USE wvfct,      ONLY : g2kin, npwx, nbnd
23  USE gvect,      ONLY : gstart
24  USE control_lr, ONLY : nbnd_occ
25  USE mp,         ONLY : mp_sum
26  USE mp_bands,   ONLY : intra_bgrp_comm
27  USE control_flags,    ONLY : gamma_only
28  USE noncollin_module, ONLY : noncolin, npol
29  !
30  IMPLICIT NONE
31  INTEGER, INTENT(in) :: ik
32  COMPLEX(dp), INTENT(in) :: evq(npwx*npol, nbnd)
33  REAL(dp), INTENT(out) :: h_diag(npwx*npol, nbnd)
34  !
35  REAL(dp), ALLOCATABLE :: eprec(:)
36  COMPLEX(dp), ALLOCATABLE :: aux(:)
37  INTEGER :: ibnd, nbnd_, ig, ikk, ikq, npwq
38  REAL(dp), EXTERNAL :: DDOT
39  !
40  ikk = ikks(ik)
41  ikq = ikqs(ik)
42  npwq = ngk(ikq)
43  nbnd_=nbnd_occ(ikk)
44
45  ALLOCATE ( eprec(nbnd_) )
46  ALLOCATE ( aux(npol*npwx) )
47  DO ibnd = 1, nbnd_
48     aux=(0.d0,0.d0)
49     DO ig = 1, npwq
50        aux (ig) = g2kin (ig) * evq (ig, ibnd)
51     END DO
52     ! NOTE: eprec(i) = 1.35*<\psi_i|Ek|\psi_i> is always real
53     IF (noncolin) THEN
54        DO ig = 1, npwq
55           aux (ig+npwx) = g2kin (ig)* evq (ig+npwx, ibnd)
56        END DO
57        eprec(ibnd) = DDOT(2*npwx*npol,evq(1,ibnd),1,aux(1),1)
58     ELSE IF ( gamma_only) THEN
59        eprec(ibnd) = 2.0_dp*DDOT(2*npwq,evq(1,ibnd),1,aux(1),1)
60        ! the following line is actually not needed
61        ! because q=0 in gamma-only case, so |k+q+G|=0 for G=0
62        IF (gstart==2) eprec(ibnd) = eprec(ibnd)-DBLE(evq(1,ibnd))*DBLE(aux(1))
63     ELSE
64        eprec(ibnd) = DDOT(2*npwq,evq(1,ibnd),1,aux(1),1)
65     END IF
66     eprec(ibnd) = 1.35_dp * eprec(ibnd)
67     !
68  END DO
69  DEALLOCATE (aux)
70  CALL mp_sum(eprec, intra_bgrp_comm)
71  !
72  h_diag=0.0_dp
73  DO ibnd = 1, nbnd_
74     DO ig = 1, npwq
75        ! Diagonal preconditining:
76        ! h_diag(G) = <Ek>/|k+q+G|^2 if |k+q+G|^2 .gt. <Ek>
77        ! h_diag(G) = 1 otherwise
78        ! written in this funny way because g2kin may be zero
79        h_diag(ig,ibnd)=1.0_dp/max(1.0_dp,g2kin(ig)/eprec(ibnd))
80     END DO
81     IF (noncolin) THEN
82        DO ig = 1, npwq
83           h_diag(ig+npwx,ibnd)=h_diag(ig,ibnd)
84        END DO
85     END IF
86  END DO
87  DEALLOCATE (eprec)
88
89END SUBROUTINE h_prec
90