File: rhoinit.f90

package info (click to toggle)
elkcode 5.4.24-2
  • links: PTS, VCS
  • area: main
  • in suites: buster, sid
  • size: 12,840 kB
  • sloc: f90: 48,415; fortran: 22,457; perl: 965; makefile: 384; sh: 369; python: 105; ansic: 67
file content (161 lines) | stat: -rw-r--r-- 3,844 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161

! Copyright (C) 2002-2005 J. K. Dewhurst, S. Sharma and C. Ambrosch-Draxl.
! This file is distributed under the terms of the GNU General Public License.
! See the file COPYING for license details.

!BOP
! !ROUTINE: rhoinit
! !INTERFACE:
subroutine rhoinit
! !USES:
use modmain
use modomp
! !DESCRIPTION:
!   Initialises the crystal charge density. Inside the muffin-tins it is set to
!   the spherical atomic density. In the interstitial region it is taken to be
!   constant such that the total charge is correct. Requires that the atomic
!   densities have already been calculated.
!
! !REVISION HISTORY:
!   Created January 2003 (JKD)
!EOP
!BOC
implicit none
! local variables
integer lmax,is,ia,ias
integer nr,nri,nro,ir,i
integer nrc,nrci,irco,irc
integer l,m,lm,ig,ifg,nthd
real(8) x,t1,t2
complex(8) z1,z2,z3
! allocatable arrays
real(8), allocatable :: jl(:,:),ffg(:),fr(:)
complex(8), allocatable :: zfmt(:),zfft(:)
! external functions
real(8) splint
external splint
lmax=min(lmaxi,1)
! zero the charge density arrays
rhomt(:,:)=0.d0
rhoir(:)=0.d0
! compute the superposition of all the atomic density tails
allocate(zfft(ngtot))
zfft(:)=0.d0
call omp_hold(nspecies,nthd)
!$OMP PARALLEL DEFAULT(SHARED) &
!$OMP PRIVATE(ffg,fr,nr,nro,ig,ir) &
!$OMP PRIVATE(x,t1,ia,ias,ifg) &
!$OMP NUM_THREADS(nthd)
!$OMP DO
do is=1,nspecies
  allocate(ffg(ngvec),fr(nrspmax))
  nr=nrmt(is)
  nro=nrsp(is)-nr+1
  do ig=1,ngvec
    do ir=nr,nrsp(is)
! spherical bessel function j_0(x)
      x=gc(ig)*rsp(ir,is)
      if (x.gt.1.d-8) then
        t1=sin(x)/x
      else
        t1=1.d0
      end if
      fr(ir)=t1*rhosp(ir,is)*r2sp(ir,is)
    end do
    t1=splint(nro,rsp(nr,is),fr(nr))
! apply low-pass filter
    t1=t1*exp(-4.d0*(gc(ig)/gmaxvr)**2)
    ffg(ig)=(fourpi/omega)*t1
  end do
  do ia=1,natoms(is)
    ias=idxas(ia,is)
!$OMP CRITICAL(rhoinit_)
    do ig=1,ngvec
      ifg=igfft(ig)
      zfft(ifg)=zfft(ifg)+ffg(ig)*conjg(sfacg(ig,ias))
    end do
!$OMP END CRITICAL(rhoinit_)
  end do
  deallocate(fr,ffg)
end do
!$OMP END DO
!$OMP END PARALLEL
call omp_free(nthd)
! compute the tails in each muffin-tin
call omp_hold(natmtot,nthd)
!$OMP PARALLEL DEFAULT(SHARED) &
!$OMP PRIVATE(jl,zfmt,is,nrc,nrci) &
!$OMP PRIVATE(irco,ig,ifg,irc,x) &
!$OMP PRIVATE(z1,z2,z3,lm,l,m,i) &
!$OMP NUM_THREADS(nthd)
!$OMP DO
do ias=1,natmtot
  allocate(jl(0:lmax,nrcmtmax),zfmt(npcmtmax))
  is=idxis(ias)
  nrc=nrcmt(is)
  nrci=nrcmti(is)
  irco=nrci+1
  zfmt(1:npcmt(is))=0.d0
  do ig=1,ngvec
    ifg=igfft(ig)
    do irc=1,nrc
      x=gc(ig)*rcmt(irc,is)
      call sbessel(lmax,x,jl(:,irc))
    end do
    z1=fourpi*zfft(ifg)*sfacg(ig,ias)
    lm=0
    do l=0,lmax
      z2=z1*zil(l)
      do m=-l,l
        lm=lm+1
        z3=z2*conjg(ylmg(lm,ig))
        i=lm
        do irc=1,nrci
          zfmt(i)=zfmt(i)+jl(l,irc)*z3
          i=i+lmmaxi
        end do
        do irc=irco,nrc
          zfmt(i)=zfmt(i)+jl(l,irc)*z3
          i=i+lmmaxo
        end do
      end do
    end do
  end do
  call ztorfmt(nrc,nrci,zfmt,rhomt(:,ias))
  deallocate(jl,zfmt)
end do
!$OMP END DO
!$OMP END PARALLEL
call omp_free(nthd)
! convert the density from a coarse to a fine radial mesh
call rfmtctof(rhomt)
! add the atomic charge density and the excess charge in each muffin-tin
t1=chgexs/omega
do ias=1,natmtot
  is=idxis(ias)
  nr=nrmt(is)
  nri=nrmti(is)
  i=1
  do ir=1,nri
    t2=(t1+rhosp(ir,is))/y00
    rhomt(i,ias)=rhomt(i,ias)+t2
    i=i+lmmaxi
  end do
  do ir=nri+1,nr
    t2=(t1+rhosp(ir,is))/y00
    rhomt(i,ias)=rhomt(i,ias)+t2
    i=i+lmmaxo
  end do
end do
! interstitial density determined from the atomic tails and excess charge
call zfftifc(3,ngridg,1,zfft)
do ir=1,ngtot
  rhoir(ir)=dble(zfft(ir))+t1
! make sure that the density is always positive
  if (rhoir(ir).lt.1.d-10) rhoir(ir)=1.d-10
end do
deallocate(zfft)
return
end subroutine
!EOC