1!
2! Copyright (C) 2004-201 Quantum ESPRESSO 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 lderiv
11  !---------------------------------------------------------------
12  !
13  !  numerical integration of the radial schroedinger equation
14  !  computing logarithmic derivatives for Coulomb potential
15  !
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 ld1_parameters, ONLY : nwfsx
22  USE ld1inc,    ONLY : file_logder, grid, vpot, rel, nspin, nld, zed, &
23                        npte, deld, eminld, emaxld, rlderiv
24
25  IMPLICIT NONE
26
27  INTEGER ::        &
28       lam,    &   ! the angular momentum
29       ikrld,  &   ! index of matching radius
30       nc,     &   ! counter on logarithmic derivatives
31       nin,    &   ! integer variable for lschps
32       is,     &   ! counter on spin
33       nstop,  &   ! integer to monitor errors
34       ios,    &   ! used for I/O control
35       n,ie        ! generic counter
36
37  real(DP) ::           &
38       aux(ndmx),         & ! the square of the wavefunction
39       aux_dir(ndmx,2),   & ! the square of the wavefunction
40       ze2,              & ! the nuclear charge in Ry units
41       e,                & ! the eigenvalue
42       j,                & ! total angular momentum for log_der
43       thrdum = 0.0_dp     ! real variable (not used) for lschps
44
45  real(DP), EXTERNAL :: compute_log
46
47  real(DP), ALLOCATABLE ::        &
48       ene(:),        &    ! the energy grid
49       dlchi(:, :)         ! the logarithmic derivative
50
51  CHARACTER(len=256) :: flld
52
53  IF (nld == 0 .or. file_logder == ' ') RETURN
54  IF (nld > nwfsx) CALL errore('lderiv','nld is too large',1)
55
56  ze2=-zed*2.0_dp
57
58  DO n=1,grid%mesh
59     IF (grid%r(n) > rlderiv) GOTO 10
60  ENDDO
61  CALL errore('lderiv','wrong rlderiv?',1)
6210 ikrld = n-1
63  WRITE(stdout,'(5x,''Computing logarithmic derivative in'',f10.5)') &
64       (grid%r(ikrld)+grid%r(ikrld+1))*0.5_dp
65
66  npte= (emaxld-eminld)/deld + 1
67  ALLOCATE ( dlchi(npte, nld) )
68  ALLOCATE ( ene(npte) )
69  DO ie=1,npte
70     ene(ie)= eminld+deld*(ie-1)
71  ENDDO
72
73  DO is=1,nspin
74     DO nc=1,nld
75        IF (rel < 2) THEN
76           lam=nc-1
77           j=0.0_dp
78        ELSE
79           lam=nc/2
80           IF (mod(nc,2)==0) j=lam-0.5_dp
81           IF (mod(nc,2)==1) j=lam+0.5_dp
82        ENDIF
83        DO ie=1,npte
84           e=ene(ie)
85           !
86           !    integrate outward up to ikrld+1
87           !
88           IF (rel == 1) THEN
89              nin = ikrld+5
90              CALL lschps (3, zed, thrdum, grid, nin, 1, lam, e, &
91                   vpot(1,is), aux, nstop )
92           ELSEIF (rel == 2) THEN
93              CALL dir_outward(ndmx,ikrld+5,lam,j,e,grid%dx,&
94                   aux_dir,grid%r,grid%rab,vpot(1,is))
95              aux(:)=aux_dir(:,1)
96           ELSE
97              CALL intref(lam,e,ikrld+5,grid,vpot(1,is),ze2,aux)
98           ENDIF
99           !
100           !    compute the logarithmic derivative and save in dlchi
101           !
102           dlchi(ie, nc) = compute_log(aux(ikrld-3),grid%r(ikrld),grid%dx)
103        ENDDO
104     ENDDO
105
106     IF (nspin == 2 .and. is == 1) THEN
107        flld = trim(file_logder)//'up'
108     ELSEIF (nspin == 2 .and. is == 2) THEN
109        flld = trim(file_logder)//'dw'
110     ELSE
111        flld = trim(file_logder)
112     ENDIF
113
114     CALL write_efun(flld,dlchi,ene,npte,nld)
115     !
116  ENDDO
117
118  DEALLOCATE (ene)
119  DEALLOCATE (dlchi)
120  RETURN
121END SUBROUTINE lderiv
122
123