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