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