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 lderivps
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_logderps, &
24                        betas, ddd, qq, lls, jjs, pseudotype, vpstot, vnl, &
25                        rlderiv, npte, emaxld, eminld, deld
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       ze2,     &    ! the nuclear charge in Ry units
37       jam,     &    ! the total angular momentum
38       e,       &    ! the eigenvalue
39       lamsq,   &    ! combined angular momentum
40       ddx12,   &    !
41       b(0:3)          ! used for starting guess of the solution
42
43  real(DP),allocatable :: &
44       ene(:),      &  ! the energy mesh
45       dlchis(:,:), &  ! the logarithmic derivatives
46       vaux(:),     &  ! auxiliary: the potential
47       aux(:),      &  ! the square of the wavefunction
48       al(:)           ! the known part of the differential equation
49
50  real(DP), external :: compute_log
51  real(DP), external :: int_0_inf_dr
52
53  integer :: &
54       ikmin,         &  ! minimum value of ik
55       ios,           &  ! used for I/O control
56       is, ind           ! counters on index
57
58  character(len=256) :: flld
59
60
61  if (nld == 0 .or. file_logderps == ' ') return
62  if (nld > nwfsx) call errore('lderivps','nld is too large',1)
63
64  allocate( al(grid%mesh), aux(grid%mesh), vaux(grid%mesh) )
65
66  ze2=0.0_dp
67
68  do n=1,grid%mesh
69     if (grid%r(n) > rlderiv) go to 10
70  enddo
71  call errore('lderivps','wrong rlderiv?',1)
7210 ikrld = n-1
73  write(stdout,'(5x,''Computing logarithmic derivative in'',f10.5)') &
74       (grid%r(ikrld)+grid%r(ikrld+1))*0.5_dp
75  npte= (emaxld-eminld)/deld + 1
76  allocate ( dlchis(npte,nld) )
77  allocate ( ene(npte) )
78  do ie=1,npte
79     ene(ie)= eminld+deld*(ie-1)
80  enddo
81
82  ikmin=ikrld+5
83  if (nbeta>0) then
84     do nbf=1,nbeta
85        ikmin=max(ikmin,ikk(nbf))
86     enddo
87  endif
88  do is=1,nspin
89     do nc=1,nld
90        if (rel < 2) then
91           lam=nc-1
92           jam=0.0_dp
93        else
94           lam=nc/2
95           if (mod(nc,2)==0) jam=lam-0.5_dp
96           if (mod(nc,2)==1) jam=lam+0.5_dp
97        endif
98
99        ddx12=grid%dx*grid%dx/12.0_dp
100        nbf=nbeta
101        if (pseudotype == 1) then
102           if (rel < 2 .or. lam == 0 .or. abs(jam-lam+0.5_dp) < 0.001_dp) then
103              ind=1
104           else if (rel==2 .and. lam>0 .and. abs(jam-lam-0.5_dp)<0.001_dp) then
105              ind=2
106           endif
107           do n=1,grid%mesh
108              vaux(n) = vpstot(n,is) + vnl(n,lam,ind)
109           enddo
110           nbf=0
111        else
112           do n=1,grid%mesh
113              vaux(n) = vpstot(n,is)
114           enddo
115        endif
116
117        do n=1,4
118           al(n)=vaux(n)-ze2/grid%r(n)
119        enddo
120        call series(al,grid%r,grid%r2,b)
121
122        do ie=1,npte
123           e=ene(ie)
124           lamsq=(lam+0.5_dp)**2
125           !
126           !     b) find the value of solution s in the first two points
127           !
128           call start_scheq( lam, e, b, grid, ze2, aux )
129
130           do n=1,grid%mesh
131              al(n)=( (vaux(n)-e)*grid%r2(n) + lamsq )*ddx12
132              al(n)=1.0_dp-al(n)
133           enddo
134
135           call integrate_outward (lam,jam,e,grid%mesh,ndmx,grid,al,b,aux,betas,ddd,&
136                qq,nbf,nwfsx,lls,jjs,ikk,ikmin)
137
138           !
139           !    compute the logarithmic derivative and save in dlchi
140           !
141           do n=-3,3
142              aux(ikrld+n)= aux(ikrld+n)*grid%sqr(ikrld+n)
143           enddo
144
145           dlchis(ie,nc)=compute_log(aux(ikrld-3),grid%r(ikrld),grid%dx)
146        enddo
147     enddo
148
149     if (nspin == 2 .and. is == 1) then
150        flld = trim(file_logderps)//'up'
151     else if (nspin == 2 .and. is == 2) then
152        flld = trim(file_logderps)//'dw'
153     else
154        flld = trim(file_logderps)
155     end if
156
157     call write_efun(flld,dlchis,ene,npte,nld)
158     !
159  enddo
160
161  deallocate(ene)
162  deallocate(dlchis)
163  deallocate(vaux, aux, al)
164
165  return
166end subroutine lderivps
167