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