1!
2! Copyright (C) 2001-2008 Quantum ESPRESSO group
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!
9!-----------------------------------------------------------------------
10subroutine compute_drhous (drhous, dbecsum, wgg, becq, alpq)
11  !-----------------------------------------------------------------------
12  !
13  !    This routine computes the part of the change of the charge density
14  !    which is due to the orthogonalization constraint on wavefunctions
15  !
16  !
17  USE kinds,      ONLY : DP
18  USE ions_base,  ONLY : nat
19  USE wavefunctions,  ONLY: evc
20  USE buffers,    ONLY : get_buffer
21  USE uspp,       ONLY : okvan, nkb, vkb
22  USE uspp_param, ONLY : nhm
23  USE lsda_mod,   ONLY : lsda, nspin, current_spin, isk
24  USE klist,      ONLY : xk, wk, ngk, igk_k
25  USE fft_base,   ONLY: dffts, dfftp
26  USE fft_interfaces, ONLY: invfft
27  USE wvfct,      ONLY : nbnd
28
29  USE qpoint,     ONLY : nksq, ikks, ikqs
30  USE eqv,        ONLY : evq
31  USE control_lr, ONLY : lgamma
32  USE units_lr,   ONLY : iuwfc, lrwfc
33  USE becmod,     ONLY : bec_type
34
35  implicit none
36  !
37  !     the dummy variables
38  !
39
40  complex(DP) :: dbecsum (nhm * (nhm + 1) / 2, nat, nspin, 3 * nat) &
41       , drhous (dfftp%nnr, nspin, 3 * nat)
42  !output:the derivative of becsum
43  ! output: add the orthogonality term
44  type (bec_type) :: becq(nksq), & ! (nkb, nbnd)
45                     alpq (3, nksq)
46  ! input: the becp with psi_{k+q}
47  ! input: the alphap with psi_{k+q}
48
49  real(DP) :: wgg (nbnd, nbnd, nksq)
50  ! input: the weights
51
52  integer :: npw, npwq, ik, ikq, ikk, ig, nu_i, ibnd, ios
53  ! counter on k points
54  ! the point k+q
55  ! record for wfcs at k point
56  ! counter on spin
57  ! counter on g vectors
58  ! counter on modes
59  ! counter on the bands
60  ! integer variable for I/O control
61
62  real(DP) :: weight
63  ! the weight of the k point
64
65  complex(DP), allocatable :: evcr (:,:)
66  ! the wavefunctions in real space
67
68  if (.not.okvan) return
69
70  call start_clock ('com_drhous')
71  allocate (evcr( dffts%nnr, nbnd))
72  !
73  drhous(:,:,:) = (0.d0, 0.d0)
74  dbecsum (:,:,:,:) = (0.d0, 0.d0)
75
76  do ik = 1, nksq
77     ikk = ikks(ik)
78     ikq = ikqs(ik)
79     npw = ngk(ikk)
80     npwq= ngk(ikq)
81     weight = wk (ikk)
82     if (lsda) current_spin = isk (ikk)
83     !
84     !   For each k point we construct the beta functions
85     !
86     call init_us_2 (npwq, igk_k(1,ikq), xk (1, ikq), vkb)
87     !
88     !   Read the wavefunctions at k and transform to real space
89     !
90     call get_buffer (evc, lrwfc, iuwfc, ikk)
91     evcr(:,:) = (0.d0, 0.d0)
92     do ibnd = 1, nbnd
93        do ig = 1, npw
94           evcr (dffts%nl (igk_k(ig,ikk) ), ibnd) = evc (ig, ibnd)
95        enddo
96        CALL invfft ('Wave', evcr (:, ibnd), dffts)
97     enddo
98     !
99     !   Read the wavefunctions at k+q
100     !
101     if (.not.lgamma.and.nksq.gt.1) call get_buffer (evq, lrwfc, iuwfc, ikq)
102     !
103     !   And compute the contribution of this k point to the change of
104     !   the charge density
105     !
106     do nu_i = 1, 3 * nat
107        call incdrhous (drhous (1, current_spin, nu_i), weight, ik, &
108             dbecsum (1, 1, current_spin, nu_i), evcr, wgg, becq, alpq, nu_i)
109     enddo
110
111  enddo
112
113  deallocate(evcr)
114  call stop_clock ('com_drhous')
115  return
116
117end subroutine compute_drhous
118