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 integrate_outward (lam,jam,e,mesh,ndm,grid,f, &
11     b,y,beta,ddd,qq,nbeta,nwfx,lls,jjs,ikk,ik)
12  !---------------------------------------------------------------------
13  !
14  !    Integrate the wavefunction from 0 to r(ik)
15  !    generalized separable or US pseudopotentials are allowed
16  !    This routine assumes that y countains already the
17  !    correct values in the first two points
18  !
19  use kinds, only : DP
20  use radial_grids, only : radial_grid_type
21  implicit none
22  type(radial_grid_type):: grid
23  integer ::   &
24       lam,    &     ! l angular momentum
25       mesh,   &    ! size of radial mesh
26       ndm,    &   ! maximum radial mesh
27       nbeta,  &   ! number of beta function
28       nwfx,   &   ! maximum number of beta functions
29       lls(nbeta),&! for each beta the angular momentum
30       ikk(nbeta),&! for each beta the integration point
31       ik         ! the last integration point
32
33  real(DP) :: &
34       e,       &  ! output eigenvalue
35       jam,     &  ! j angular momentum
36       f(mesh), &  ! the f function
37       b(0:3), &   ! the taylor expansion of the potential
38       y(mesh), &  ! the output solution
39       jjs(nwfx), & ! the j angular momentum
40       beta(ndm,nwfx),& ! the beta functions
41       ddd(nwfx,nwfx),qq(nwfx,nwfx) ! parameters for computing B_ij
42
43  integer ::  &
44       nst, &      ! the exponential around the origin
45       n,    &     ! counter on mesh points
46       iib,jjb, &  ! counter on beta with correct lam
47       ierr,    &  ! used to control allocation
48       ib,jb,   &  ! counter on beta
49       info      ! info on exit of LAPACK subroutines
50
51  integer, allocatable :: iwork(:) ! auxiliary space
52
53  real(DP) :: &
54       b0e,     & ! the expansion of the known part
55       ddx12,   & ! the deltax entering the equations
56       x4l6,    & ! auxiliary for small r expansion
57       j1(4),d(4),& ! auxiliary for starting values of chi
58       delta,xc(4),& ! auxiliary for starting values of eta
59       int_0_inf_dr  ! the integral function
60
61  real(DP), allocatable :: &
62       el(:), &  ! auxiliary for integration
63       cm(:,:), &! the linear system
64       bm(:), & ! the known part of the linear system
65       c(:), &   ! the chi functions
66       coef(:), & ! the solution of the linear system
67       eta(:,:) ! the partial solution of the nonomogeneous
68
69
70  allocate(c(ik), stat=ierr)
71  allocate(el(ik), stat=ierr)
72  allocate(cm(nbeta,nbeta), stat=ierr)
73  allocate(bm(nbeta), stat=ierr)
74  allocate(coef(nbeta), stat=ierr)
75  allocate(iwork(nbeta), stat=ierr)
76  allocate(eta(ik,nbeta), stat=ierr)
77
78  if (mesh.ne.grid%mesh) call errore('integrate_outward','mesh dimension is not as expected',1)
79  ddx12=grid%dx*grid%dx/12.0_DP
80  b0e=b(0)-e
81  x4l6=4*lam+6
82  nst=(lam+1)*2
83  !
84  !  first solve the homogeneous equation
85  !
86  ! f is the original function
87  ! of the form 1 + h^2/12 * d2Rdr2
88  do n=2,ik-1
89     y(n+1)=( 12.0_DP*y(n)  - 10.0_DP* f(n)*y(n) - f(n-1)*y(n-1) )/f(n+1)
90  enddo
91  !
92  !     for each beta function with correct angular momentum
93  !     solve the inhomogeneous equation
94  !
95  iib=0
96  jjb=0
97  do ib=1,nbeta
98     !
99     !    set up the known part
100     !
101     if (lls(ib).eq.lam.and.jjs(ib).eq.jam) then
102        iib=iib+1
103        c=0.0_DP
104        do jb=1,nbeta
105           if (lls(jb).eq.lam.and.jjs(jb).eq.jam) then
106              do n=1,ikk(jb)
107                 c(n)= c(n)+(ddd(jb,ib) -e*qq(jb,ib))*beta(n,jb)
108              enddo
109           endif
110        enddo
111        !
112        !     compute the starting values of the solutions
113        !
114        do n=1,4
115           j1(n)=c(n)/grid%r(n)**(lam+1)
116        enddo
117        call seriesbes(j1,grid%r,grid%r2,4,d)
118        delta=b0e**2+x4l6*b(2)
119        xc(1)=(-d(1)*b0e-x4l6*d(3))/delta
120        xc(3)=(-b0e*d(3)+d(1)*b(2))/delta
121        xc(2)=0.0_DP
122        xc(4)=0.0_DP
123        do n=1,3
124           eta(n,iib)=grid%r(n)**(lam+1)*(xc(1)+grid%r2(n)*xc(3))/grid%sqr(n)
125        enddo
126
127        do n=1,ik
128           c(n)=c(n)*grid%r2(n)/grid%sqr(n)
129        enddo
130        !
131        !   solve the inhomogeneous equation
132        !
133        do n=3,ik-1
134           eta(n+1,iib)=((12.0_DP-10.0_DP*f(n))*eta(n,iib) &
135                -f(n-1)*eta(n-1,iib) &
136                + ddx12*(10.0_DP*c(n)+c(n-1)+c(n+1)) )/f(n+1)
137        enddo
138        !
139        !    compute the coefficents of the linear system
140        !
141        jjb=0
142        do jb=1,nbeta
143           if (lls(jb).eq.lam.and.jjs(jb).eq.jam) then
144              jjb=jjb+1
145              do n=1,min(ik,ikk(jb))
146                 el(n)=beta(n,jb)*eta(n,iib)*grid%sqr(n)
147              enddo
148              cm(jjb,iib)=-int_0_inf_dr(el,grid,min(ik,ikk(jb)),nst)
149           endif
150        enddo
151
152        do n=1,min(ik,ikk(ib))
153           el(n)=beta(n,ib)*y(n)*grid%sqr(n)
154        enddo
155        bm(iib)=int_0_inf_dr(el,grid,min(ik,ikk(ib)),nst)
156        cm(iib,iib)=1.0_DP+cm(iib,iib)
157     endif
158  enddo
159  if (iib.ne.jjb) call errore('integrate_outward','jjb.ne.iib',1)
160
161  if (iib.gt.0) then
162     call dcopy(iib,bm,1,coef,1)
163
164     call DGESV(iib,1,cm,nbeta,iwork,coef,nbeta,info)
165
166     if (info /= 0) call errore('integrate_outward', &
167                &  'problems solving the linear system',info)
168
169     do ib=1,iib
170        do n=1,ik
171           y(n)= y(n)+coef(ib)*eta(n,ib)
172        enddo
173     enddo
174  endif
175
176  deallocate(eta)
177  deallocate(iwork)
178  deallocate(coef)
179  deallocate(bm)
180  deallocate(cm)
181  deallocate(el)
182  deallocate(c)
183
184  return
185end subroutine integrate_outward
186