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