File: gradzfmt.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 (221 lines) | stat: -rw-r--r-- 6,126 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

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

!BOP
! !ROUTINE: gradzfmt
! !INTERFACE:
subroutine gradzfmt(nr,nri,r,zfmt,ld,gzfmt)
! !USES:
use modmain
! !INPUT/OUTPUT PARAMETERS:
!   nr    : number of radial mesh points (in,integer)
!   nri   : number of points on inner part of muffin-tin (in,integer)
!   r     : radial mesh (in,real(nr))
!   zfmt  : complex muffin-tin function (in,complex(*))
!   ld    : leading dimension (in,integer)
!   gzfmt : gradient of zfmt (out,complex(ld,3))
! !DESCRIPTION:
!   Calculates the gradient of a complex muffin-tin function. In other words,
!   given the spherical harmonic expansion coefficients, $f_{lm}(r)$, of a
!   function $f({\bf r})$, the routine returns ${\bf F}_{lm}$ where
!   $$ \sum_{lm}{\bf F}_{lm}(r)Y_{lm}(\hat{\bf r})=\nabla f({\bf r}). $$
!   This is done using the formula (see, for example, V. Devanathan,
!   {\em Angular Momentum Techniques In Quantum Mechanics})
!   \begin{align*}
!    \nabla_{\mu}^s f_{lm}(r)Y_{lm}(\hat{\bf r})&=\sqrt{\frac{l+1}{2l+3}}
!    C(l,1,l+1|m,\mu,m+\mu)Y_{l+1m+\mu}(\hat{\bf r})\left(\frac{d}{dr}
!    -\frac{l}{r}\right)f_{lm}(r)\\
!    &-\sqrt{\frac{l}{2l-1}}C(l,1,l-1|m,\mu,m+\mu)Y_{l-1,m+\mu}(\hat{\bf r})
!    \left(\frac{d}{dr}+\frac{l+1}{r}\right)f_{lm}(r),
!   \end{align*}
!   where $C$ are Clebsch-Gordan coefficients and the gradient $\nabla_{\mu}^s$
!   is in terms of the spherical unit vectors $\hat{\bf e}_{\mu}$:
!   $$ \hat{\bf e}_{+1}=-\frac{\hat{\bf x}+i\hat{\bf y}}{\sqrt{2}},
!    \qquad\hat{\bf e}_0=\hat{\bf z},\qquad
!    \hat{\bf e}_{-1}=\frac{\hat{\bf x}-i\hat{\bf y}}{\sqrt{2}}. $$
!   Note that the gradient returned is in terms of the global
!   $(\hat{\bf x},\hat{\bf y},\hat{\bf z})$ coordinate system.
!
! !REVISION HISTORY:
!   Rewritten May 2009 (JKD)
!EOP
!BOC
implicit none
! arguments
integer, intent(in) :: nr,nri
real(8), intent(in) :: r(nr)
complex(8), intent(in) :: zfmt(*)
integer, intent(in) :: ld
complex(8), intent(out) :: gzfmt(ld,3)
! local variables
integer nro,iro,ir,i,j
integer np,npi,i1,i2
integer l,m,lm1,lm2
! real constant 1/sqrt(2)
real(8), parameter :: c1=0.7071067811865475244d0
real(8) t1,t2,t3
complex(8) z1
! automatic arrays
real(8) ri(nr),f1(nr),f2(nr),g1(nr),g2(nr)
! allocatable arrays
complex(8), allocatable :: drmt(:)
! external functions
real(8) clebgor
external clebgor
nro=nr-nri
iro=nri+1
npi=lmmaxi*nri
np=npi+lmmaxo*nro
! store 1/r
ri(:)=1.d0/r(:)
!----------------------------------------!
!     compute the radial derivatives     !
!----------------------------------------!
allocate(drmt(np))
do lm1=1,lmmaxi
  i1=lm1
  do ir=1,nri
    f1(ir)=dble(zfmt(i1))
    f2(ir)=aimag(zfmt(i1))
    i1=i1+lmmaxi
  end do
  do ir=iro,nr
    f1(ir)=dble(zfmt(i1))
    f2(ir)=aimag(zfmt(i1))
    i1=i1+lmmaxo
  end do
  call fderiv(1,nr,r,f1,g1)
  call fderiv(1,nr,r,f2,g2)
  i1=lm1
  do ir=1,nri
    drmt(i1)=cmplx(g1(ir),g2(ir),8)
    i1=i1+lmmaxi
  end do
  do ir=iro,nr
    drmt(i1)=cmplx(g1(ir),g2(ir),8)
    i1=i1+lmmaxo
  end do
end do
do lm1=lmmaxi+1,lmmaxo
  i1=npi+lm1
  do ir=iro,nr
    f1(ir)=dble(zfmt(i1))
    f2(ir)=aimag(zfmt(i1))
    i1=i1+lmmaxo
  end do
  call fderiv(1,nro,r(iro),f1(iro),g1(iro))
  call fderiv(1,nro,r(iro),f2(iro),g2(iro))
  i1=npi+lm1
  do ir=iro,nr
    drmt(i1)=cmplx(g1(ir),g2(ir),8)
    i1=i1+lmmaxo
  end do
end do
!-------------------------------------------------------!
!     compute the gradient in spherical coordinates     !
!-------------------------------------------------------!
! zero the gradient array
gzfmt(1:np,:)=0.d0
! inner part of muffin-tin
lm1=0
do l=0,lmaxi
  t1=sqrt(dble(l+1)/dble(2*l+3))
  if (l.gt.0) then
    t2=-sqrt(dble(l)/dble(2*l-1))
  else
    t2=0.d0
  end if
  do m=-l,l
    lm1=lm1+1
    j=1
    do i=-1,1
      if (i.eq.0) j=3
      if (i.eq.1) j=2
      if ((l+1.le.lmaxi).and.(abs(m+i).le.l+1)) then
! index to (l,m) is l*(l+1)+m+1, therefore index to (l+1,m+i) is
        lm2=(l+1)*(l+2)+(m+i)+1
        t3=t1*clebgor(l,1,l+1,m,i,m+i)
        i1=lm1; i2=lm2
        do ir=1,nri
          gzfmt(i2,j)=gzfmt(i2,j)+t3*(drmt(i1)-dble(l)*ri(ir)*zfmt(i1))
          i1=i1+lmmaxi; i2=i2+lmmaxi
        end do
      end if
      if (abs(m+i).le.l-1) then
! index to (l-1,m+i)
        lm2=(l-1)*l+(m+i)+1
        t3=t2*clebgor(l,1,l-1,m,i,m+i)
        i1=lm1; i2=lm2
        do ir=1,nri
          gzfmt(i2,j)=gzfmt(i2,j)+t3*(drmt(i1)+dble(l+1)*ri(ir)*zfmt(i1))
          i1=i1+lmmaxi; i2=i2+lmmaxi
        end do
      end if
    end do
  end do
end do
! outer part of muffin-tin
lm1=0
do l=0,lmaxo
  t1=sqrt(dble(l+1)/dble(2*l+3))
  if (l.gt.0) then
    t2=-sqrt(dble(l)/dble(2*l-1))
  else
    t2=0.d0
  end if
  do m=-l,l
    lm1=lm1+1
    j=1
    do i=-1,1
      if (i.eq.0) j=3
      if (i.eq.1) j=2
      if ((l+1.le.lmaxo).and.(abs(m+i).le.l+1)) then
        lm2=(l+1)*(l+2)+(m+i)+1
        t3=t1*clebgor(l,1,l+1,m,i,m+i)
        i1=npi+lm1; i2=npi+lm2
        do ir=iro,nr
          gzfmt(i2,j)=gzfmt(i2,j)+t3*(drmt(i1)-dble(l)*ri(ir)*zfmt(i1))
          i1=i1+lmmaxo; i2=i2+lmmaxo
        end do
      end if
      if (abs(m+i).le.l-1) then
        lm2=(l-1)*l+(m+i)+1
        t3=t2*clebgor(l,1,l-1,m,i,m+i)
        i1=npi+lm1; i2=npi+lm2
        do ir=iro,nr
          gzfmt(i2,j)=gzfmt(i2,j)+t3*(drmt(i1)+dble(l+1)*ri(ir)*zfmt(i1))
          i1=i1+lmmaxo; i2=i2+lmmaxo
        end do
      end if
    end do
  end do
end do
deallocate(drmt)
!--------------------------------------------------------!
!     convert from spherical components to Cartesian     !
!--------------------------------------------------------!
i1=0
do ir=1,nri
  do lm1=1,lmmaxi
    i1=i1+1
    z1=gzfmt(i1,1)
    gzfmt(i1,1)=c1*(z1-gzfmt(i1,2))
    z1=c1*(z1+gzfmt(i1,2))
    gzfmt(i1,2)=cmplx(-aimag(z1),dble(z1),8)
  end do
end do
do ir=iro,nr
  do lm1=1,lmmaxo
    i1=i1+1
    z1=gzfmt(i1,1)
    gzfmt(i1,1)=c1*(z1-gzfmt(i1,2))
    z1=c1*(z1+gzfmt(i1,2))
    gzfmt(i1,2)=cmplx(-aimag(z1),dble(z1),8)
  end do
end do
return
end subroutine
!EOC