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 partial_wave_expansion
11  !---------------------------------------------------------------
12  !
13  !  numerical integration of the radial schroedinger equation
14  !  computing logarithmic derivatives for pseudo-potentials
15  !  multiple nonlocal projectors are allowed
16  !
17  use kinds,     only : DP
18  use radial_grids, only : ndmx
19  use io_global, only : stdout
20  use mp,        only : mp_bcast
21  use radial_grids, only: series
22  use ld1_parameters, only : nwfsx
23  use ld1inc,    only : grid, nld, nbeta, nspin, rel, ikk, file_pawexp, &
24                        betas, ddd, qq, lls, jjs, pseudotype, vpstot, vnl, &
25                        rpwe, npte, emaxld, eminld, deld, phis, rcutus, rcloc
26  implicit none
27
28  integer  ::       &
29       lam,   &      ! the angular momentum
30       ikrld, &      ! index of matching radius
31       nc,    &      ! counter on logarithmic derivatives
32       nbf,   &      ! number of b functions
33       n,ie          ! generic counters
34
35  real(DP) ::  &
36       norm, norm2,& ! auxiliary variables
37       ze2,     &    ! the nuclear charge in Ry units
38       jam,     &    ! the total angular momentum
39       e,       &    ! the eigenvalue
40       lamsq,   &    ! combined angular momentum
41       b(0:3),  &    ! used for starting guess of the solution
42       ddx12
43
44  real(DP),allocatable :: &
45       ene(:),      &  ! the energy grid
46       nnn(:,:),    &  ! the expansion fraction
47       vaux(:),     &  ! auxiliary: the potential
48       aux(:),      &  ! the square of the wavefunction
49       aux2(:),     &  ! auxiliary wavefunction
50       al(:)           ! the known part of the differential equation
51
52  real(DP), external :: compute_log
53  real(DP), external :: int_0_inf_dr
54
55  integer :: &
56       ik, jb,        &  ! auxiliary variables
57       ikmin,         &  ! minimum value of ik
58       nst,           &  ! auxiliary for integrals
59       ios,           &  ! used for I/O control
60       is, ind           ! counters on index
61
62  logical :: found
63
64  character(len=256) :: flld
65
66
67  if (nld == 0 .or. file_pawexp == ' ') return
68  if (nld > nwfsx) call errore('partial_wave_expansion','nld is too large',1)
69
70  allocate( al(grid%mesh), aux(grid%mesh), aux2(grid%mesh),vaux(grid%mesh) )
71
72  ze2=0.0_dp
73
74  do n=1,grid%mesh
75     if (grid%r(n) > rpwe) go to 10
76  enddo
77  call errore('partial_wave_expansion','wrong rpwe?',1)
7810 ikrld = n-1
79  write(stdout,'(5x,''Computing the partial wave expansion '')')
80  npte= (emaxld-eminld)/deld + 1
81  allocate ( nnn(npte,nld) )
82  allocate ( ene(npte) )
83  do ie=1,npte
84     ene(ie)=eminld+deld*(ie-1)
85  enddo
86  ikmin=ikrld+5
87  if (nbeta>0) then
88     do nbf=1,nbeta
89        ikmin=max(ikmin,ikk(nbf))
90     enddo
91  endif
92  spinloop : &
93  do is=1,nspin
94     nlogdloop : &
95     do nc=1,nld
96        if (rel < 2) then
97           lam=nc-1
98           jam=0.0_dp
99        else
100           lam=nc/2
101           if (mod(nc,2)==0) jam=lam-0.5_dp
102           if (mod(nc,2)==1) jam=lam+0.5_dp
103        endif
104        ddx12=grid%dx*grid%dx/12.0_dp
105        nst=(lam+1)*2
106        nbf=nbeta
107        if (pseudotype == 1) then
108           if (rel < 2 .or. lam == 0 .or. abs(jam-lam+0.5_dp) < 0.001_dp) then
109              ind=1
110           else if (rel==2 .and. lam>0 .and. abs(jam-lam-0.5_dp)<0.001_dp) then
111              ind=2
112           endif
113           do n=1,grid%mesh
114              vaux(n) = vpstot(n,is) + vnl(n,lam,ind)
115           enddo
116           nbf=0
117        else
118           do n=1,grid%mesh
119              vaux(n) = vpstot(n,is)
120           enddo
121        endif
122
123        do n=1,4
124           al(n)=vaux(n)-ze2/grid%r(n)
125        enddo
126        call series(al,grid%r,grid%r2,b)
127
128        energiesloop : &
129        do ie=1,npte
130           e=eminld+deld*(ie-1.0_dp)
131           lamsq=(lam+0.5_dp)**2
132           !
133           !     b) find the value of solution s in the first two points
134           !
135           call start_scheq( lam, e, b, grid, ze2, aux )
136
137           do n=1,grid%mesh
138              al(n)=( (vaux(n)-e)*grid%r2(n) + lamsq )*ddx12
139              al(n)=1.0_dp-al(n)
140           enddo
141
142           !
143           ! integrate outward at fixed energy
144           !
145           call integrate_outward (lam,jam,e,grid%mesh,ndmx,grid,al,b,aux,betas,ddd,&
146                qq,nbf,nwfsx,lls,jjs,ikk,ikmin)
147           !
148           aux(1:ikmin) = aux(1:ikmin)*grid%sqr(1:ikmin)
149           !
150           ! calculate the accuracy of the expasion of the solution at
151           ! a given energy in terms of the partial waves at the chosen
152           ! reference energies.
153           ! a vanishing value of nnn indicates a perfect expansion
154           !
155           aux2(:) = 0.d0
156           ik=min(ikrld,ikmin)
157           if (mod(ik,2)==0) ik=ik+1
158           found = .false.
159           do jb=1,nbeta
160              if (lls(jb).eq.lam.and.jjs(jb).eq.jam) then
161                 found = .true.
162                 IF (grid%r(ik)>max(rcutus(jb),rcloc).and.ie==1) &
163                    write(stdout, &
164                 '(5x,"R is outside the sphere for l=",i5," j=",f5.2)') lam,jam
165                 al(1:ikk(jb))=betas(1:ikk(jb),jb)*aux(1:ikk(jb))
166                 norm = int_0_inf_dr(al,grid,ikk(jb),nst)
167                 aux2(1:ik) = aux2(1:ik) + phis(1:ik,jb)*norm
168              endif
169           enddo
170           if( .not. found) then
171               nnn(:,nc) = 0._dp
172               write(stdout, '(7x,a,i3)') "no projector for channel: ", lam
173               cycle nlogdloop
174           endif
175
176           aux2(ik+1:ikmin) = 0._dp
177           al(1:ik)= aux(1:ik)**2
178           !
179           norm=int_0_inf_dr(al,grid,ik,nst)
180           !
181           al(1:ik)= (aux2(1:ik)-aux(1:ik))**2
182           norm2=int_0_inf_dr(al,grid,ik,nst)
183           nnn(ie,nc) = norm2/norm
184           !
185           !
186        enddo energiesloop
187     enddo nlogdloop
188
189     flld = trim(file_pawexp)
190     nnn=100.0_DP*nnn
191     call write_efun(flld,nnn,ene,npte,nld)
192  enddo spinloop
193
194  deallocate(ene)
195  deallocate(nnn)
196  deallocate(vaux, aux, al)
197
198  return
199end subroutine partial_wave_expansion
200