1!
2! Copyright (C) 2004 PWSCF 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 compute_chi(lam,ikk_in,phi_in,chi_out,xc,e,lbes4)
11  !--------------------------------------------------------------------------
12  !
13  !     This routine computes the chi functions:
14  !          |chi> = (\epsilon -T -V_{loc}) |psi>
15  !
16  use io_global, only : stdout, ionode
17  use kinds, only : DP
18  use radial_grids, only: ndmx, series
19  use ld1inc, only:  grid, vpsloc, rho0, verbosity
20
21  implicit none
22  integer :: &
23       ikk_in,& ! the point after which the chi should be zero
24       lam      ! the angular momentum
25  logical :: &
26       lbes4
27
28  real(DP) :: &
29       e,     &       ! input: the energy
30       xc(8), &       ! input: the coefficients of the bessel function
31       phi_in(ndmx), & ! input: pseudo wavefunction
32       chi_out(ndmx)   ! output: the chi function
33  !
34  real(DP) :: &
35       j1(ndmx), aux(ndmx), gi(ndmx),&
36       b(4),c(4), arow(ndmx),brow(ndmx),crow(ndmx),drow(ndmx), &
37       b0e, g0, g1, g2, &
38       r_clean,         &
39       ddx12,           &
40       x4l6,            &
41       x6l12,           &
42       int_0_inf_dr,    &
43       integral
44
45  integer :: &
46       n, nstart, nst
47  !
48  !   RRKJ: first expand in a taylor series the phis function
49  !   Since we know that the phis functions are a sum of Bessel
50  !   functions with coefficients xc, we compute analytically
51  !   the asymptotic expansion
52  !
53  !
54  ddx12=grid%dx*grid%dx/12.0_dp
55  x4l6=4*lam+6
56  x6l12=6*lam+12
57
58  do n=1,6
59     j1(n)=phi_in(n)/grid%r(n)**(lam+1)
60  enddo
61  call seriesbes(j1,grid%r,grid%r2,6,c)
62  !
63  if (lam == 0) then
64     if(lbes4.or.rho0.eq.0.0_dp)then
65        c(1)=xc(1)+xc(2)+xc(3)
66        c(2)=0.0_dp
67        c(3)=-xc(1)*(xc(4)**2/6.0_dp) &
68             -xc(2)*(xc(5)**2/6.0_dp) &
69             -xc(3)*(xc(6)**2/6.0_dp)
70        c(4)=0.0_dp
71     else
72        c(1)=xc(1)+xc(2)+xc(3)+xc(4)
73        c(2)=0.0_dp
74        c(3)=-xc(1)*(xc(5)**2/6.0_dp)  &
75             -xc(2)*(xc(6)**2/6.0_dp)  &
76             -xc(3)*(xc(7)**2/6.0_dp)  &
77             -xc(4)*(xc(8)**2/6.0_dp)
78        c(4)=0.0_dp
79     endif
80  elseif (lam == 3) then
81     c(1)=xc(1)*(48.0_dp*xc(4)**3/5040.0_dp)+   &
82          xc(2)*(48.0_dp*xc(5)**3/5040.0_dp)+   &
83          xc(3)*(48.0_dp*xc(6)**3/5040.0_dp)
84     c(2)=0.0_dp
85     c(3)=-xc(1)*(192.0_dp*xc(4)**5/362880.0_dp)  &
86          -xc(2)*(192.0_dp*xc(5)**5/362880.0_dp)  &
87          -xc(3)*(192.0_dp*xc(6)**5/362880.0_dp)
88     c(4)=0.0_dp
89  elseif (lam == 2) then
90     c(1)=xc(1)*(xc(4)**2/15.0_dp)+   &
91          xc(2)*(xc(5)**2/15.0_dp)+   &
92          xc(3)*(xc(6)**2/15.0_dp)
93     c(2)=0.0_dp
94     c(3)=-xc(1)*(xc(4)**4/210.0_dp)  &
95          -xc(2)*(xc(5)**4/210.0_dp)  &
96          -xc(3)*(xc(6)**4/210.0_dp)
97     c(4)=0.0_dp
98  elseif (lam == 1) then
99     c(1)=xc(1)*(xc(4)/3.0_dp)+  &
100          xc(2)*(xc(5)/3.0_dp)+  &
101          xc(3)*(xc(6)/3.0_dp)
102     c(2)=0.0_dp
103     c(3)=-xc(1)*(xc(4)**3/30.0_dp) &
104          -xc(2)*(xc(5)**3/30.0_dp) &
105          -xc(3)*(xc(6)**3/30.0_dp)
106     c(4)=0.0_dp
107  else
108     call errore('compute_chi','lam not programmed',1)
109  endif
110  !
111  !     and the potential
112  !
113  do n=1,4
114     j1(n)=vpsloc(n)
115  enddo
116  if (abs(j1(1)-j1(4))>1.d-12) then
117      call series(j1,grid%r,grid%r2,b)
118  else
119     b=0.0_dp
120     b(1)=j1(1)
121  endif
122  !
123  !   and compute the taylor expansion of the chis
124  !
125  b0e=(b(1)-e)
126
127  g0=x4l6*c(3)-b0e*c(1)
128  g1=x6l12*c(4)-c(1)*b(2)
129  g2=-(b0e*c(3)+b(3)*c(1))
130  nstart=5
131  do n=1,nstart-1
132     chi_out(n)= (g0+grid%r(n)*(g1+g2*grid%r(n)))*grid%r(n)**(lam+3)/grid%sqr(n)
133  enddo
134  do n=1,grid%mesh
135     aux(n)= (g0+grid%r(n)*(g1+g2*grid%r(n)))
136  enddo
137  !
138  !    set up the equation
139  !
140  do n=1,grid%mesh
141     gi(n)=phi_in(n)/grid%sqr(n)
142  enddo
143  do n=1,grid%mesh
144     j1(n)=grid%r2(n)*(vpsloc(n)-e)+(lam+0.5_dp)**2
145     j1(n)=1.0_dp-ddx12*j1(n)
146  enddo
147
148  do n=nstart,grid%mesh-3
149     drow(n)= gi(n+1)*j1(n+1)   &
150            + gi(n)*(-12.0_dp+10.0_dp*j1(n))+ &
151              gi(n-1)*j1(n-1)
152
153     brow(n)=10.0_dp*ddx12
154     crow(n)=ddx12
155     arow(n)=ddx12
156  enddo
157  drow(nstart)=drow(nstart)-ddx12*chi_out(nstart-1)
158  chi_out(grid%mesh-2)=0.0_dp
159  chi_out(grid%mesh-1)=0.0_dp
160  chi_out(grid%mesh)=0.0_dp
161  !
162  !    and solve it
163  !
164  call tridiag(arow(nstart),brow(nstart),crow(nstart), &
165       drow(nstart),chi_out(nstart),grid%mesh-3-nstart)
166  !
167  !   put the correct normalization and r dependence
168  !
169  do n=1,grid%mesh
170     chi_out(n)=chi_out(n)*grid%sqr(n)/grid%r2(n)
171     !         if(lam.eq.0)
172     !     +    write(stdout,'(5(e20.13,1x))')
173     !     +          r(n),chi_out(n),chi_out(n)/r(n)**(lam+1),
174     !     +          aux(n),aux(n)*r(n)**(lam+1)
175  enddo
176  !
177  !    smooth close to the origin with asymptotic expansion
178  !
179  do n=nstart,grid%mesh
180     if (abs(chi_out(n)/grid%r(n)**(lam+1)-aux(n))  &
181          .lt.1.e-3_dp*abs(aux(n)) ) goto 100
182     chi_out(n)=aux(n)*grid%r(n)**(lam+1)
183  enddo
184
185100 if (n.eq.grid%mesh+1.or.grid%r(min(n,grid%mesh)).gt.0.05_dp)then
186     write(stdout,*) lam,n,grid%mesh,grid%r(min(n,grid%mesh))
187     call errore('compute_chi','n is too large',1)
188  endif
189!
190!    When the input wavefunction is a diverging scattering state
191!    this routine might become numerically unstable at large r.
192!    Here chi_out should be 0.0 but it is not due to this instability.
193!    We now clean chi_out when phi_in > 20.0.
194!    Clean also after 7.0 a.u..
195!
196  r_clean=100.d0
197  do n=1,grid%mesh
198     if (abs(phi_in(n))>20.0_DP) then
199        r_clean=grid%r(n)
200        exit
201     endif
202  enddo
203  r_clean=min(r_clean,7.0_DP)
204  if (r_clean<grid%r(ikk_in)) &
205     call errore ('compute_chi ','phi_in too large before r_c', 1)
206
207  do n=grid%mesh,1,-1
208     if (grid%r(n).lt.r_clean) goto 200
209     chi_out(n)=0.0_dp
210  enddo
211200 continue
212     !    check that the chi are zero beyond ikk
213  nst=0
214  gi=0.0_dp
215  do n=ikk_in+1,grid%mesh
216     gi(n)=chi_out(n)**2
217  enddo
218  do n=min(ikk_in+20,grid%mesh),grid%mesh
219     chi_out(n)=0.0_dp
220  enddo
221  integral=int_0_inf_dr(gi,grid,grid%mesh,nst)
222  if (integral > 2.e-6_dp) then
223      write(stdout, '(5x,'' l='',i4, '' integral='',f15.9, &
224           & '' r(ikk) '',f15.9)') lam, integral, grid%r(ikk_in)
225      IF (verbosity=='high') THEN
226         do n=ikk_in,grid%mesh
227            write(stdout,*) grid%r(n),gi(n)
228         enddo
229      ENDIF
230     call errore ('compute_chi ','chi too large beyond r_c', 1)
231  endif
232  return
233end subroutine compute_chi
234
235
236subroutine tridiag(a,b,c,r,u,n)
237  !
238  !     See Numerical Recipes.
239  !
240  use kinds, only : DP
241  implicit none
242
243  integer :: n
244  real(DP) :: a(n),b(n),c(n),r(n),u(n)
245  real(DP) :: gam(n), bet
246
247  integer j
248
249  if (abs(b(1)).lt.1.e-10_DP)  &
250       call errore('tridiag','b(1) is too small',1)
251
252  bet=b(1)
253  u(1)=r(1)/bet
254  do j=2,n
255     gam(j)=c(j-1)/bet
256     bet=b(j)-a(j)*gam(j)
257     if (abs(bet) < 1.e-10_DP) &
258          call errore('tridiag','bet is too small',1)
259     u(j)=(r(j)-a(j)*u(j-1))/bet
260  enddo
261  do j=n-1,1,-1
262     u(j)=u(j)-gam(j+1)*u(j+1)
263  enddo
264  return
265end subroutine tridiag
266