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