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