1 2! Copyright (C) 2002-2005 J. K. Dewhurst, S. Sharma and C. Ambrosch-Draxl. 3! This file is distributed under the terms of the GNU General Public License. 4! See the file COPYING for license details. 5 6!BOP 7! !ROUTINE: rhoinit 8! !INTERFACE: 9subroutine rhoinit 10! !USES: 11use modmain 12use modomp 13! !DESCRIPTION: 14! Initialises the crystal charge density. Inside the muffin-tins it is set to 15! the spherical atomic density. In the interstitial region it is taken to be 16! constant such that the total charge is correct. Requires that the atomic 17! densities have already been calculated. 18! 19! !REVISION HISTORY: 20! Created January 2003 (JKD) 21!EOP 22!BOC 23implicit none 24! local variables 25integer lmax,is,ia,ias,i 26integer nr,nri,nro,nrs,ir 27integer nrc,nrci,irco,irc 28integer l,m,lm,ig,ifg,nthd 29real(8) x,sm,t1,t2 30complex(8) z1,z2,z3 31! automatic arrays 32real(8) ffg(ngvc),wr(nrspmax),jl(0:lmaxi,nrcmtmax) 33complex(8) zfmt(npcmtmax) 34! allocatable arrays 35complex(8), allocatable :: zfft(:) 36lmax=min(lmaxi,1) 37! compute the superposition of all the atomic density tails 38allocate(zfft(ngtot)) 39zfft(:)=0.d0 40call holdthd(nspecies,nthd) 41!$OMP PARALLEL DO DEFAULT(SHARED) & 42!$OMP PRIVATE(ffg,wr,nr,nrs,nro,ig) & 43!$OMP PRIVATE(t1,t2,sm,ir,x,ia,ias,ifg) & 44!$OMP NUM_THREADS(nthd) 45do is=1,nspecies 46 nr=nrmt(is) 47 nrs=nrsp(is) 48 nro=nrs-nr+1 49! determine the weights for the radial integral 50 call wsplint(nro,rsp(nr,is),wr(nr)) 51 do ig=1,ngvc 52 t1=gc(ig) 53! spherical bessel function j_0(x) times the atomic density tail 54 if (t1.gt.epslat) then 55 t2=1.d0/t1 56 sm=0.d0 57 do ir=nr,nrs 58 x=t1*rsp(ir,is) 59 sm=sm+t2*sin(x)*rhosp(ir,is)*rsp(ir,is)*wr(ir) 60 end do 61 else 62 sm=sum(rhosp(nr:nrs,is)*(rsp(nr:nrs,is)**2)*wr(nr:nrs)) 63 end if 64! apply low-pass filter 65 t1=sm*exp(-4.d0*(gc(ig)/gmaxvr)**2) 66 ffg(ig)=(fourpi/omega)*t1 67 end do 68 do ia=1,natoms(is) 69 ias=idxas(ia,is) 70!$OMP CRITICAL(rhoinit_) 71 do ig=1,ngvc 72 ifg=igfft(ig) 73 zfft(ifg)=zfft(ifg)+ffg(ig)*conjg(sfacg(ig,ias)) 74 end do 75!$OMP END CRITICAL(rhoinit_) 76 end do 77end do 78!$OMP END PARALLEL DO 79call freethd(nthd) 80! compute the tails in each muffin-tin 81call holdthd(natmtot,nthd) 82!$OMP PARALLEL DO DEFAULT(SHARED) & 83!$OMP PRIVATE(jl,zfmt,is,nrc,nrci) & 84!$OMP PRIVATE(irco,ig,ifg,irc,x) & 85!$OMP PRIVATE(z1,z2,z3,lm,l,m,i) & 86!$OMP NUM_THREADS(nthd) 87do ias=1,natmtot 88 is=idxis(ias) 89 nrc=nrcmt(is) 90 nrci=nrcmti(is) 91 irco=nrci+1 92 zfmt(1:npcmt(is))=0.d0 93 do ig=1,ngvc 94 ifg=igfft(ig) 95 do irc=1,nrc 96 x=gc(ig)*rcmt(irc,is) 97 call sbessel(lmax,x,jl(:,irc)) 98 end do 99 z1=fourpi*zfft(ifg)*sfacg(ig,ias) 100 lm=0 101 do l=0,lmax 102 z2=z1*zil(l) 103 do m=-l,l 104 lm=lm+1 105 z3=z2*conjg(ylmg(lm,ig)) 106 i=lm 107 do irc=1,nrci 108 zfmt(i)=zfmt(i)+jl(l,irc)*z3 109 i=i+lmmaxi 110 end do 111 do irc=irco,nrc 112 zfmt(i)=zfmt(i)+jl(l,irc)*z3 113 i=i+lmmaxo 114 end do 115 end do 116 end do 117 end do 118 call ztorfmt(nrc,nrci,zfmt,rhomt(:,ias)) 119end do 120!$OMP END PARALLEL DO 121call freethd(nthd) 122! convert the density from a coarse to a fine radial mesh 123call rfmtctof(rhomt) 124! add the atomic charge density and the excess charge in each muffin-tin 125t1=chgexs/omega 126do ias=1,natmtot 127 is=idxis(ias) 128 nr=nrmt(is) 129 nri=nrmti(is) 130 i=1 131 do ir=1,nri 132 t2=(t1+rhosp(ir,is))/y00 133 rhomt(i,ias)=rhomt(i,ias)+t2 134 i=i+lmmaxi 135 end do 136 do ir=nri+1,nr 137 t2=(t1+rhosp(ir,is))/y00 138 rhomt(i,ias)=rhomt(i,ias)+t2 139 i=i+lmmaxo 140 end do 141end do 142! interstitial density determined from the atomic tails and excess charge 143call zfftifc(3,ngridg,1,zfft) 144do ir=1,ngtot 145 rhoir(ir)=dble(zfft(ir))+t1 146! make sure that the density is always positive 147 if (rhoir(ir).lt.1.d-10) rhoir(ir)=1.d-10 148end do 149deallocate(zfft) 150end subroutine 151!EOC 152 153