1 2! Copyright (C) 2008 J. K. Dewhurst, S. Sharma and C. Ambrosch-Draxl. 3! This file is distributed under the terms of the GNU General Public License. 4! See the file COPYING for license details. 5 6subroutine elnes 7use modmain 8use modomp 9use modtest 10implicit none 11! local variables 12integer ik,jk,ikq,isym,nsk(3) 13integer ist,jst,iw,n,nthd 14real(8) vgqc(3),gqc 15real(8) vkql(3),v(3) 16real(8) q,wd,dw,w,t1 17! allocatable arrays 18real(8), allocatable :: jlgqr(:,:),ddcs(:) 19real(8), allocatable :: e(:,:,:),f(:,:,:) 20complex(8), allocatable :: ylmgq(:),sfacgq(:) 21complex(8), allocatable :: expmt(:,:),emat(:,:) 22! initialise universal variables 23call init0 24call init1 25call init2 26! check q-vector is commensurate with k-point grid 27v(:)=dble(ngridk(:))*vecql(:) 28v(:)=abs(v(:)-nint(v(:))) 29if ((v(1).gt.epslat).or.(v(2).gt.epslat).or.(v(3).gt.epslat)) then 30 write(*,*) 31 write(*,'("Error(elnes): q-vector incommensurate with k-point grid")') 32 write(*,'(" ngridk : ",3I6)') ngridk 33 write(*,'(" vecql : ",3G18.10)') vecql 34 write(*,*) 35 stop 36end if 37! read in the density and potentials from file 38call readstate 39! read Fermi energy from file 40call readfermi 41! find the new linearisation energies 42call linengy 43! generate the APW radial functions 44call genapwfr 45! generate the local-orbital radial functions 46call genlofr 47! get the second-variational eigenvalues and occupancies from file 48do ik=1,nkpt 49 call getevalsv(filext,ik,vkl(:,ik),evalsv(:,ik)) 50 call getoccsv(filext,ik,vkl(:,ik),occsv(:,ik)) 51end do 52! generate the phase factor function exp(iq.r) in the muffin-tins 53allocate(jlgqr(njcmax,nspecies)) 54allocate(ylmgq(lmmaxo),sfacgq(natmtot)) 55allocate(expmt(npcmtmax,natmtot)) 56call gengqf(1,vecqc,vgqc,gqc,jlgqr,ylmgq,sfacgq) 57call genexpmt(1,jlgqr,ylmgq,1,sfacgq,expmt) 58deallocate(jlgqr,ylmgq,sfacgq) 59allocate(e(nstsv,nstsv,nkptnr),f(nstsv,nstsv,nkptnr)) 60e(:,:,:)=0.d0 61f(:,:,:)=0.d0 62! begin parallel loop over non-reduced k-points 63call holdthd(nkptnr,nthd) 64!$OMP PARALLEL DEFAULT(SHARED) & 65!$OMP PRIVATE(emat,jk,vkql,isym) & 66!$OMP PRIVATE(ikq,ist,jst,t1) & 67!$OMP NUM_THREADS(nthd) 68allocate(emat(nstsv,nstsv)) 69!$OMP DO 70do ik=1,nkptnr 71!$OMP CRITICAL(elnes_) 72 write(*,'("Info(elnes): ",I6," of ",I6," k-points")') ik,nkptnr 73!$OMP END CRITICAL(elnes_) 74! equivalent reduced k-point 75 jk=ivkik(ivk(1,ik),ivk(2,ik),ivk(3,ik)) 76! k+q-vector in lattice coordinates 77 vkql(:)=vkl(:,ik)+vecql(:) 78! index to k+q-vector 79 call findkpt(vkql,isym,ikq) 80! compute < i,k+q | exp(iq.r) | j,k > matrix elements 81 call genexpmat(vkl(:,ik),expmt,emat) 82! add to the double differential scattering cross-section 83 do jst=1,nstsv 84 if (evalsv(jst,jk).lt.emaxelnes) then 85 do ist=1,nstsv 86 e(ist,jst,ik)=evalsv(ist,ikq)-evalsv(jst,jk) 87 t1=dble(emat(ist,jst))**2+aimag(emat(ist,jst))**2 88 f(ist,jst,ik)=t1*occsv(jst,jk)*(occmax-occsv(ist,ikq)) 89 end do 90 end if 91 end do 92end do 93!$OMP END DO 94deallocate(emat) 95!$OMP END PARALLEL 96call freethd(nthd) 97! number of subdivisions used for interpolation 98nsk(:)=max(ngrkf/ngridk(:),1) 99n=nstsv*nstsv 100! integrate over the Brillouin zone 101allocate(ddcs(nwplot)) 102call brzint(nswplot,ngridk,nsk,ivkiknr,nwplot,wplot,n,n,e,f,ddcs) 103q=sqrt(vecqc(1)**2+vecqc(2)**2+vecqc(3)**2) 104t1=2.d0/(omega*occmax) 105if (q.gt.epslat) t1=t1/q**4 106ddcs(:)=t1*ddcs(:) 107open(50,file='ELNES.OUT',form='FORMATTED') 108wd=wplot(2)-wplot(1) 109dw=wd/dble(nwplot) 110do iw=1,nwplot 111 w=dw*dble(iw-1)+wplot(1) 112 write(50,'(2G18.10)') w,ddcs(iw) 113end do 114close(50) 115write(*,*) 116write(*,'("Info(elnes):")') 117write(*,'(" ELNES double differential cross-section written to ELNES.OUT")') 118! write ELNES distribution to test file 119call writetest(140,'ELNES cross-section',nv=nwplot,tol=1.d-2,rva=ddcs) 120deallocate(e,f,ddcs,expmt) 121end subroutine 122 123