1!
2! Copyright (C) 2007 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 chargeps(rho_i,phi_i,nwf_i,ll_i,jj_i,oc_i,iswf_i)
11  !---------------------------------------------------------------
12  !
13  !   calculate the (spherical) pseudo charge density
14  !
15  use kinds, only: dp
16  use ld1_parameters, only: nwfsx
17  use radial_grids, only: ndmx
18  use ld1inc, only: grid, pseudotype, qvan, nbeta, betas, lls, jjs, ikk,  &
19                    which_augfun, lpaw, qvanl, nspin
20  implicit none
21
22  integer :: &
23       nwf_i,        & ! input: the number of wavefunctions
24       iswf_i(nwfsx),& ! input: their spin
25       ll_i(nwfsx)     ! input: their angular momentum
26
27  real(DP) ::  &
28       jj_i(nwfsx), & ! input: their total angular momentum
29       oc_i(nwfsx), & ! input: the occupation
30       phi_i(ndmx,nwfsx), & ! input: the functions to add
31       rho_i(ndmx,2)   ! output: the (nspin) components of the charge
32
33  integer ::     &
34       is,     &   ! counter on spin
35       n,n1,n2,&   ! counters on beta and mesh function
36       ns,nst,ikl  ! counter on wavefunctions
37
38  real(DP) ::    &
39       work(nwfsx), & ! auxiliary variable for becp
40       int_0_inf_dr,& ! integration function
41       gi(ndmx)        ! used to compute the integrals
42
43
44  rho_i=0.0_dp
45  !
46  !    compute the square modulus of the eigenfunctions
47  !
48  do ns=1,nwf_i
49     if (oc_i(ns).gt.0.0_dp) then
50        is=iswf_i(ns)
51        do n=1,grid%mesh
52           rho_i(n,is)=rho_i(n,is)+oc_i(ns)*phi_i(n,ns)**2
53        end do
54     endif
55  enddo
56  !
57  !    if US pseudopotential compute the augmentation part
58  !
59  if (pseudotype.eq.3) then
60     do ns=1,nwf_i
61        if (oc_i(ns).gt.0.0_dp) then
62           is=iswf_i(ns)
63           do n1=1,nbeta
64              if (ll_i(ns).eq.lls(n1).and. &
65                   abs(jj_i(ns)-jjs(n1)).lt.1.e-7_dp) then
66                 nst=(ll_i(ns)+1)*2
67                 ikl=ikk(n1)
68                 do n=1,ikl
69                    gi(n)=betas(n,n1)*phi_i(n,ns)
70                 enddo
71                 work(n1)=int_0_inf_dr(gi,grid,ikl,nst)
72              else
73                 work(n1)=0.0_dp
74              endif
75           enddo
76           !
77           !   and adding to the charge density
78           !
79           do n1=1,nbeta
80              do n2=1,nbeta
81                 IF (which_augfun=='PSQ') then
82                    do n=1,grid%mesh
83                       rho_i(n,is)=rho_i(n,is)+qvanl(n,n1,n2,0)*oc_i(ns)* &
84                            work(n1)*work(n2)
85                    enddo
86                 ELSE
87                    do n=1,grid%mesh
88                       rho_i(n,is)=rho_i(n,is)+qvan(n,n1,n2)*oc_i(ns)* &
89                            work(n1)*work(n2)
90                    enddo
91                 ENDIF
92              enddo
93           enddo
94        endif
95     enddo
96  endif
97!
98!  Check for negative charge
99!
100  do is=1,nspin
101     do n=2,grid%mesh
102        if (rho_i(n,is)<-1.d-12) call errore('chargeps','negative rho',1)
103     enddo
104  enddo
105
106  return
107end subroutine chargeps
108