File: zpotcoul.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 (296 lines) | stat: -rw-r--r-- 10,260 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
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296

! 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: zpotcoul
! !INTERFACE:
subroutine zpotcoul(nr,nri,np,npi,ld1,r,ngdg,igf,ngp,gpc,gclgp,ld2,jlgprmt, &
 ylmgp,sfacgp,zrhoir,ld3,zvclmt,zvclir)
! !USES:
use modmain
use modphonon
! !INPUT/OUTPUT PARAMETERS:
!   nr      : number of radial points for each species (in,integer(nspecies))
!   nri     : number of radial points on inner part (in,integer(nspecies))
!   np      : total number of points in muffin-tins (in,integer(nspecies))
!   npi     : number of points on inner part (in,integer(nspecies))
!   ld1     : leading dimension (in,integer)
!   r       : radial mesh for each species (in,real(ld1,nspecies))
!   ngdg    : G-vector grid sizes (in,integer(3))
!   igf     : map from G-vector index to FFT array (in,integer(*))
!   ngp     : number of G+p-vectors (in,integer)
!   gpc     : G+p-vector lengths (in,real(ngp))
!   gclgp   : Coulomb Green's function in G+p-space (in,real(ngp))
!   ld2     : leading dimension (in,integer)
!   jlgprmt : spherical Bessel functions for evergy G+p-vector and muffin-tin
!             radius (in,real(0:lnpsd,ld2,nspecies))
!   ylmgp   : spherical harmonics of the G+p-vectors (in,complex(lmmaxo,ngp))
!   sfacgp  : structure factors of the G+p-vectors (in,complex(ld2,natmtot))
!   zrhoir  : interstitial charge density (in,complex(*))
!   ld3     : leading dimension (in,integer)
!   zvclmt  : muffin-tin Coulomb potential, with the contribution from the
!             isolated muffin-tin density precalculated and passed in
!             (inout,complex(ld3,natmtot))
!   zvclir  : interstitial Coulomb potential (out,complex(*))
! !DESCRIPTION:
!   Calculates the Coulomb potential of a complex charge density by solving
!   Poisson's equation using the method of M. Weinert, {\it J. Math. Phys.}
!   {\bf 22}, 2433 (1981). First, the multipole moments of the muffin-tin charge
!   are determined for the $j$th atom of the $i$th species by
!   $$ q_{ij;lm}^{\rm MT}=\int_0^{R_i}r^{l+2}\rho_{ij;lm}(r)dr+z_{ij}Y_{00}
!   \,\delta_{l,0}\;, $$
!   where $R_i$ is the muffin-tin radius and $z_{ij}$ is a point charge located
!   at the atom center (usually the nuclear charge, which should be taken as
!   {\bf negative}). Next, the multipole moments of the continuation of the
!   interstitial density, $\rho^{\rm I}$, into the muffin-tin are found with
!   $$ q_{ij;lm}^{\rm I}=4\pi i^l R_i^{l+3}\sum_{\bf G}\frac{j_{l+1}(GR_i)}
!    {GR_i}\rho^{\rm I}({\bf G})\exp(i{\bf G}\cdot{\bf r}_{ij})Y_{lm}^*
!    (\hat{\bf G}), $$
!   remembering that
!   $$ \lim_{x\rightarrow 0}\frac{j_{l+n}(x)}{x^n}=\frac{1}{(2n+1)!!}
!    \delta_{l,0} $$
!   should be used for the case ${\bf G}=0$. A pseudocharge is now constructed
!   which is equal to the real density in the interstitial region and whose
!   multipoles are the difference between the real and interstitial muffin-tin
!   multipoles. This pseudocharge density is smooth in the sense that it can be
!   expanded in terms of the finite set of ${\bf G}$-vectors. In each muffin-tin
!   the pseudocharge has the form
!   $$ \rho_{ij}^{\rm P}({\bf r})=\rho^{\rm I}({\bf r}-{\bf r}_{ij})+\sum_{lm}
!    \rho_{ij;lm}^{\rm P}\frac{1}{R_i^{l+3}}\left(\frac{r}{R_i}\right)^l\left(1-
!    \frac{r^2}{R_i^2}\right)^{N_i}Y_{lm}(\hat{\bf r}) $$
!   where
!   $$ \rho_{ij;lm}^{\rm P}=\frac{(2l+2N_i+3)!!}{2^N_iN_i!(2l+1)!!}\left(
!    q_{ij;lm}^{\rm MT}-q_{ij;lm}^{\rm I}\right) $$
!   and $N_i\approx\frac{1}{2}R_iG_{\rm max}$ is generally a good choice.
!   The pseudocharge in reciprocal space is given by
!   $$ \rho^{\rm P}({\bf G})=\rho^{\rm I}({\bf G})+\sum_{ij;lm}2^{N_i}N_i!
!    \frac{4\pi(-i)^l}{\Omega R_i^l}\frac{j_{l+N_i+1}(GR_i)}{(GR_i)^{N_i+1}}
!    \rho_{ij;lm}^{\rm P}\exp(-i{\bf G}\cdot{\bf r}_{ij})Y_{lm}(\hat{\bf G}) $$
!   which may be used for solving Poisson's equation directly
!   $$ V^{\rm P}({\bf G})=\begin{cases}
!     4\pi\frac{\rho^{\rm P}({\bf G})}{G^2} & G>0 \\
!     0 & G=0 \end{cases}\;. $$
!   The usual Green's function approach is then employed to determine the
!   potential in the muffin-tin sphere due to charge in the sphere. In other
!   words
!   $$ V_{ij;lm}^{\rm MT}(r)=\frac{4\pi}{2l+1}\left(\frac{1}{r^{l+1}}\int_0^r
!    \rho_{ij;lm}^{\rm MT}(r'){r'}^{l+2}dr'+r^l\int_r^{R_i}\frac{
!    \rho_{ij;lm}^{\rm MT}(r')}{{r'}^{l-1}}dr'\right)+\frac{1}{Y_{00}}
!    \frac{z_{ij}}{r}\delta_{l,0} $$
!   where the last term is the monopole arising from the point charge. All that
!   remains is to add the homogenous solution of Poisson's equation,
!   $$ V_{ij}^{\rm H}({\bf r})=\sum_{lm}V_{ij;lm}^{\rm H}\left(\frac{r}
!    {R_i}\right)^lY_{lm}(\hat{\bf r}), $$
!   to the muffin-tin potential so that it is continuous at the muffin-tin
!   boundary. Therefore the coefficients, $\rho_{ij;lm}^{\rm H}$, are given by
!   $$ V_{ij;lm}^{\rm H}=4\pi i^l\sum_{\bf G}j_{l}(Gr)V^{\rm P}({\bf G})
!    \exp(i{\bf G}\cdot{\bf r}_{ij})Y_{lm}^*(\hat{\bf G})-V_{ij;lm}^{\rm MT}
!    (R_i). $$
!   Finally note that the ${\bf G}$-vectors passed to the routine can represent
!   vectors with a non-zero offset, ${\bf G}+{\bf p}$ say, which is required for
!   calculating Coulomb matrix elements.
!
! !REVISION HISTORY:
!   Created April 2003 (JKD)
!EOP
!BOC
implicit none
! arguments
integer, intent(in) :: nr(nspecies),nri(nspecies)
integer, intent(in) :: np(nspecies),npi(nspecies)
integer, intent(in) :: ld1
real(8), intent(in) :: r(ld1,nspecies)
integer, intent(in) :: ngdg(3),igf(*),ngp
real(8), intent(in) :: gpc(ngp),gclgp(ngp)
integer, intent(in) :: ld2
real(8), intent(in) :: jlgprmt(0:lnpsd,ld2,nspecies)
complex(8), intent(in) :: ylmgp(lmmaxo,ngp),sfacgp(ld2,natmtot)
complex(8), intent(in) :: zrhoir(*)
integer, intent(in) :: ld3
complex(8), intent(inout) :: zvclmt(ld3,natmtot)
complex(8), intent(out) :: zvclir(*)
! local variables
integer is,ia,ias
integer ir,l,m,lm
integer ig,jg,i,j
real(8) t0,t1,t2,t3
complex(8) z1,z2
! automatic arrays
real(8) rmtl(0:lmaxo+3,nspecies),rl(ld1,0:lmaxo)
complex(8) qlm(lmmaxo,natmtot)
complex(8) zl(0:lmaxo),zlm(lmmaxo)
complex(8) zhmt(ld3)
! external functions
real(8) factnm
external factnm
! compute (R_mt)^l
do is=1,nspecies
  rmtl(0,is)=1.d0
  do l=1,lmaxo+3
    rmtl(l,is)=rmtl(l-1,is)*rmt(is)
  end do
end do
! compute the multipole moments from the muffin-tin potentials
do ias=1,natmtot
  is=idxis(ias)
  i=np(is)-lmmaxo
  lm=0
  do l=0,lmaxo
    t1=dble(2*l+1)*rmtl(l+1,is)/fourpi
    do m=-l,l
      lm=lm+1
      i=i+1
      qlm(lm,ias)=t1*zvclmt(i,ias)
    end do
  end do
end do
! Fourier transform density to G-space and store in zvclir
call zcopy(ngdg(1)*ngdg(2)*ngdg(3),zrhoir,1,zvclir,1)
call zfftifc(3,ngdg,-1,zvclir)
! subtract the multipole moments of the interstitial charge density
do is=1,nspecies
  do l=0,lmaxo
    zl(l)=fourpi*zil(l)*rmtl(l+2,is)
  end do
  do ia=1,natoms(is)
    ias=idxas(ia,is)
    do ig=1,ngp
      jg=igf(ig)
      if (gpc(ig).gt.epslat) then
        z1=zvclir(jg)*sfacgp(ig,ias)/gpc(ig)
        lm=0
        do l=0,lmaxo
          z2=jlgprmt(l+1,ig,is)*z1*zl(l)
          do m=-l,l
            lm=lm+1
            qlm(lm,ias)=qlm(lm,ias)-z2*conjg(ylmgp(lm,ig))
          end do
        end do
      else
        t1=(fourpi/3.d0)*rmtl(3,is)*y00
        qlm(1,ias)=qlm(1,ias)-t1*zvclir(jg)
      end if
    end do
  end do
end do
! find the smooth pseudocharge within the muffin-tin whose multipoles are the
! difference between the real muffin-tin and interstitial multipoles
t0=fourpi/omega
t1=factnm(2*lnpsd+1,2)
do ias=1,natmtot
  is=idxis(ias)
  lm=0
  do l=0,lmaxo
    t2=t1/(factnm(2*l+1,2)*rmtl(l,is))
    z1=t2*zilc(l)
    do m=-l,l
      lm=lm+1
      zlm(lm)=z1*qlm(lm,ias)
    end do
  end do
! add the pseudocharge and real interstitial densities in G-space
  do ig=1,ngp
    jg=igf(ig)
    if (gpc(ig).gt.epslat) then
      t2=gpc(ig)*rmt(is)
      t3=1.d0/t2**lnpsd
      z1=t3*zlm(1)*ylmgp(1,ig)
      lm=1
      do l=1,lmaxo
        lm=lm+1
        z2=zlm(lm)*ylmgp(lm,ig)
        do m=1-l,l
          lm=lm+1
          z2=z2+zlm(lm)*ylmgp(lm,ig)
        end do
        t3=t3*t2
        z1=z1+t3*z2
      end do
      z2=t0*jlgprmt(lnpsd,ig,is)*conjg(sfacgp(ig,ias))
      zvclir(jg)=zvclir(jg)+z1*z2
    else
      t2=t0*y00/factnm(2*lnpsd+1,2)
      zvclir(jg)=zvclir(jg)+t2*zlm(1)
    end if
  end do
end do
! solve Poisson's equation in G+p-space for the pseudocharge
do ig=1,ngp
  jg=igf(ig)
  zvclir(jg)=gclgp(ig)*zvclir(jg)
end do
! match potentials at muffin-tin boundary by adding homogeneous solution
do is=1,nspecies
! compute (r/R_mt)^l
  t1=1.d0/rmt(is)
  do ir=1,nr(is)
    t2=t1*r(ir,is)
    rl(ir,0)=1.d0
    do l=1,lmaxo
      rl(ir,l)=rl(ir,l-1)*t2
    end do
  end do
  do ia=1,natoms(is)
    ias=idxas(ia,is)
! find the spherical harmonic expansion of the interstitial potential at the
! muffin-tin radius
    zlm(:)=0.d0
    do ig=1,ngp
      jg=igf(ig)
      z1=fourpi*zvclir(jg)*sfacgp(ig,ias)
      lm=0
      do l=0,lmaxo
        z2=jlgprmt(l,ig,is)*z1*zil(l)
        do m=-l,l
          lm=lm+1
          zlm(lm)=zlm(lm)+z2*conjg(ylmgp(lm,ig))
        end do
      end do
    end do
! calculate the homogenous solution
    i=np(is)-lmmaxo
    lm=0
    do l=0,lmaxi
      do m=-l,l
        lm=lm+1
        i=i+1
        z1=zlm(lm)-zvclmt(i,ias)
        j=lm
        do ir=1,nri(is)
          zhmt(j)=z1*rl(ir,l)
          j=j+lmmaxi
        end do
        do ir=nri(is)+1,nr(is)
          zhmt(j)=z1*rl(ir,l)
          j=j+lmmaxo
        end do
      end do
    end do
    do l=lmaxi+1,lmaxo
      do m=-l,l
        lm=lm+1
        i=i+1
        z1=zlm(lm)-zvclmt(i,ias)
        j=npi(is)+lm
        do ir=nri(is)+1,nr(is)
          zhmt(j)=z1*rl(ir,l)
          j=j+lmmaxo
        end do
      end do
    end do
    zvclmt(1:np(is),ias)=zvclmt(1:np(is),ias)+zhmt(1:np(is))
! store the nuclear potential without the self-term for the phonon dynamical
! matrix calculation
    if (tphdyn) then
      if (ias.eq.iasph) zvnmt(1:np(is))=zhmt(1:np(is))
    end if
  end do
end do
! Fourier transform interstitial potential to real-space
call zfftifc(3,ngdg,1,zvclir)
return
end subroutine
!EOC