1!
2! Copyright (C) 2001-2016 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!---------------------------------------------------------------------------
9SUBROUTINE lr_calc_dens_eels (drhoscf, dpsi)
10  !--------------------------------------------------------------------------
11  !
12  ! Calculates response charge-density from linear-response
13  ! orbitals and ground-state orbitals.
14  ! See Eq.(36) in B. Walker and R. Gebauer, J. Chem. Phys. 127, 164106 (2007)
15  !
16  ! It does:
17  ! 1) Calculates the response charge-density
18  ! 2) Calculates the contribution to the charge-density due to US PP's
19  ! 3) Sums up the normal and ultrasoft terms
20  ! 4) Symmetrizes the total response charge-density
21  !
22  ! Written by Iurii Timrov (2013)
23  !
24  USE kinds,                 ONLY : DP
25  USE ions_base,             ONLY : nat
26  USE gvect,                 ONLY : ngm, g
27  USE fft_base,              ONLY : dfftp, dffts
28  USE klist,                 ONLY : xk, wk, igk_k, ngk
29  USE wvfct,                 ONLY : nbnd, npwx
30  USE gvecw,                 ONLY : gcutw
31  USE qpoint,                ONLY : nksq, ikks, ikqs
32  USE wavefunctions,  ONLY : evc
33  USE uspp_param,            ONLY : nhm
34  USE uspp,                  ONLY : okvan, vkb
35  USE mp_global,             ONLY : inter_pool_comm, intra_bgrp_comm
36  USE mp,                    ONLY : mp_sum
37  USE io_files,              ONLY : iunwfc, nwordwfc
38  USE buffers,               ONLY : get_buffer
39  USE fft_interfaces,        ONLY : fft_interpolate
40  !
41  IMPLICIT NONE
42  !
43  COMPLEX(DP), INTENT(in) :: dpsi(npwx,nbnd,nksq)
44  ! input: the perturbed wavefunctions
45  COMPLEX(DP), INTENT(out) :: drhoscf(dfftp%nnr)
46  ! input/output: the accumulated change to the charge density on the dense mesh
47  !
48  COMPLEX(DP), ALLOCATABLE :: drhoscfh(:), dbecsum(:,:)
49  ! the accumulated change to the charge density on the smooth mesh
50  ! the ultrasoft term
51  INTEGER ::  ik,  &
52              ikk, & ! index of the point k
53              ikq, & ! index of the point k+q
54              npwq   ! number of the plane-waves at point k+q
55  REAL(DP) :: weight ! weight of the k point
56  !
57  CALL start_clock('lr_calc_dens')
58  !
59  ALLOCATE ( drhoscfh(dffts%nnr) )
60  drhoscfh(:) = (0.0d0, 0.0d0)
61  !
62  IF (okvan) THEN
63     ALLOCATE (dbecsum(nhm*(nhm+1)/2,nat))
64     dbecsum(:,:) = (0.0d0, 0.0d0)
65  ENDIF
66  !
67  DO ik = 1, nksq
68     !
69     ikk  = ikks(ik)
70     ikq  = ikqs(ik)
71     npwq = ngk(ikq)
72     !
73     ! Read the unperturbed wavefuctions evc at k
74     !
75     IF (nksq > 1) CALL get_buffer (evc, nwordwfc, iunwfc, ikk)
76     !
77     ! The weight of the k point
78     !
79     weight = wk(ikk)
80     !
81     ! Calculate beta-functions vkb at point k+q
82     !
83     IF (okvan) CALL init_us_2 (npwq, igk_k(1,ikq), xk(1,ikq), vkb)
84     !
85     ! Calculation of the response charge density
86     !
87     CALL incdrhoscf(drhoscfh(:), weight, ik, dbecsum(:,:), dpsi(:,:,ik))
88     !
89  ENDDO
90  !
91  ! Interpolate the charge density from the smooth mesh
92  ! to a thicker mesh (if doublegrid=.true.)
93  ! drhoscfh -> drhoscf
94  !
95  CALL fft_interpolate(dffts, drhoscfh, dfftp, drhoscf)
96  !
97  IF (okvan) THEN
98     !
99     ! The calculation of dbecsum is distributed across processors (see addusdbec).
100     ! Sum over processors the contributions coming from each slice of bands.
101     !
102#if defined(__MPI)
103     CALL mp_sum (dbecsum, intra_bgrp_comm)
104#endif
105     !
106     ! Calculate the total charge density response
107     ! (sum up the normal and ultrasoft terms)
108     !
109     CALL lr_addusddens (drhoscf, dbecsum)
110     !
111  ENDIF
112  !
113  ! Reduce the response charge density across pools.
114  !
115#if defined(__MPI)
116  CALL mp_sum (drhoscf, inter_pool_comm)
117#endif
118  !
119  ! Symmetrization of the charge density response
120  ! wrt the small group of q.
121  !
122#if defined(__MPI)
123  CALL lr_psym_eels(drhoscf)
124#else
125  CALL lr_sym_eels(drhoscf)
126#endif
127  !
128  DEALLOCATE (drhoscfh)
129  IF (okvan) DEALLOCATE (dbecsum)
130  !
131  CALL stop_clock('lr_calc_dens')
132  !
133  RETURN
134  !
135END SUBROUTINE lr_calc_dens_eels
136