1
2! Copyright (C) 2002-2005 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
6!BOP
7! !ROUTINE: elfplot
8! !INTERFACE:
9subroutine elfplot
10! !USES:
11use modmain
12! !DESCRIPTION:
13!   Outputs the electron localisation function (ELF) for 1D, 2D or 3D plotting.
14!   The spin-averaged ELF is given by
15!   $$ f_{\rm ELF}({\bf r})=\frac{1}{1+[D({\bf r})/D^0({\bf r})]^2}, $$
16!   where
17!   $$ D({\bf r})=\frac{1}{2}\left(\tau({\bf r})-\frac{1}{4}
18!    \frac{[\nabla n({\bf r})]^2}{n({\bf r})}\right) $$
19!   and
20!   $$ \tau({\bf r})=\sum_{i=1}^N \left|\nabla\Psi_i({\bf r})
21!    \right|^2 $$
22!   is the spin-averaged kinetic energy density from the spinor wavefunctions.
23!   The function $D^0$ is the kinetic energy density for the homogeneous
24!   electron gas evaluated for $n({\bf r})$:
25!   $$ D^0({\bf r})=\frac{3}{5}(6\pi^2)^{2/3}\left(\frac{n({\bf r})}{2}
26!    \right)^{5/3}. $$
27!   The ELF is useful for the topological classification of bonding. See for
28!   example T. Burnus, M. A. L. Marques and E. K. U. Gross [Phys. Rev. A 71,
29!   10501 (2005)].
30!
31! !REVISION HISTORY:
32!   Created September 2003 (JKD)
33!   Fixed bug found by F. Wagner (JKD)
34!EOP
35!BOC
36implicit none
37! local variables
38integer ik,is,ias
39integer nr,nri,ir
40integer ig,ifg,i
41real(8) r,t1,t2
42complex(8) z1
43! allocatable arrays
44real(8), allocatable :: gwf2mt(:,:),gwf2ir(:)
45real(8), allocatable :: rfmt1(:),rfmt2(:),grfir(:)
46real(8), allocatable :: grfmt1(:,:),grfmt2(:,:)
47real(8), allocatable :: elfmt(:,:),elfir(:)
48complex(8), allocatable :: evecfv(:,:),evecsv(:,:)
49complex(8), allocatable :: zfft1(:),zfft2(:)
50! initialise universal variables
51call init0
52call init1
53! allocate local arrays
54allocate(gwf2mt(npmtmax,natmtot),gwf2ir(ngtot))
55allocate(rfmt1(npmtmax),rfmt2(npmtmax),grfir(ngtot))
56allocate(grfmt1(npmtmax,3),grfmt2(npmtmax,3))
57allocate(elfmt(npmtmax,natmtot),elfir(ngtot))
58allocate(evecfv(nmatmax,nstfv),evecsv(nstsv,nstsv))
59allocate(zfft1(ngtot),zfft2(ngtot))
60! read density and potentials from file
61call readstate
62! generate the core wavefunctions and densities
63call gencore
64! read Fermi energy from file
65call readfermi
66! find the new linearisation energies
67call linengy
68! generate the APW radial functions
69call genapwfr
70! generate the local-orbital radial functions
71call genlofr
72! set the gradient squared to zero
73gwf2mt(:,:)=0.d0
74gwf2ir(:)=0.d0
75do ik=1,nkpt
76! get the eigenvectors and occupancies from file
77  call getoccsv(filext,ik,vkl(:,ik),occsv(:,ik))
78  call getevecfv(filext,ik,vkl(:,ik),vgkl(:,:,:,ik),evecfv)
79  call getevecsv(filext,ik,vkl(:,ik),evecsv)
80! add the valence wavefunction gradient squared
81  call gradwf2(ik,evecfv,evecsv,gwf2mt,gwf2ir)
82end do
83! add core wavefunction gradient squared
84call gradwfcr2(gwf2mt)
85!------------------------!
86!     muffin-tin ELF     !
87!------------------------!
88do ias=1,natmtot
89  is=idxis(ias)
90  nr=nrmt(is)
91  nri=nrmti(is)
92! convert rho from spherical harmonics to spherical coordinates
93  call rbsht(nr,nri,rhomt(:,ias),rfmt1)
94! compute the gradient of the density
95  call gradrfmt(nr,nri,rlmt(:,-1,is),wcrmt(:,:,is),rhomt(:,ias),npmtmax,grfmt1)
96! convert gradient to spherical coordinates
97  do i=1,3
98    call rbsht(nr,nri,grfmt1(:,i),grfmt2(:,i))
99  end do
100  do i=1,npmt(is)
101    r=abs(rfmt1(i))
102! square of gradient of rho
103    t1=grfmt2(i,1)**2+grfmt2(i,2)**2+grfmt2(i,3)**2
104! D for inhomogeneous density
105    t1=(1.d0/2.d0)*(gwf2mt(i,ias)-(1.d0/4.d0)*t1/r)
106! D0 for uniform electron gas
107    t2=(3.d0/5.d0)*((6.d0*pi**2)**(2.d0/3.d0))*(r/2.d0)**(5.d0/3.d0)
108! ELF function
109    rfmt2(i)=1.d0/(1.d0+(t1/t2)**2)
110  end do
111! convert ELF from spherical coordinates to spherical harmonics
112  call rfsht(nr,nri,rfmt2,elfmt(:,ias))
113end do
114!--------------------------!
115!     interstitial ELF     !
116!--------------------------!
117! Fourier transform density to G-space
118zfft1(:)=rhoir(:)
119call zfftifc(3,ngridg,-1,zfft1)
120! calculate the square of gradient of rho
121grfir(:)=0.d0
122do i=1,3
123  zfft2(:)=0.d0
124  do ig=1,ngvec
125    ifg=igfft(ig)
126! take the gradient
127    z1=zfft1(ifg)
128    zfft2(ifg)=vgc(i,ig)*cmplx(-aimag(z1),dble(z1),8)
129  end do
130! Fourier transform gradient to real-space
131  call zfftifc(3,ngridg,1,zfft2)
132  do ir=1,ngtot
133    grfir(ir)=grfir(ir)+dble(zfft2(ir))**2
134  end do
135end do
136do ir=1,ngtot
137  r=abs(rhoir(ir))
138! D for inhomogeneous density
139  t1=(1.d0/2.d0)*(gwf2ir(ir)-(1.d0/4.d0)*grfir(ir)/r)
140! D0 for homogeneous electron gas
141  t2=(3.d0/5.d0)*((6.d0*pi**2)**(2.d0/3.d0))*(r/2.d0)**(5.d0/3.d0)
142! ELF function
143  elfir(ir)=1.d0/(1.d0+(t1/t2)**2)
144end do
145! symmetrise the ELF
146call symrf(nrmt,nrmti,npmt,ngridg,ngtot,ngvec,igfft,npmtmax,elfmt,elfir)
147! plot the ELF to file
148select case(task)
149case(51)
150  open(50,file='ELF1D.OUT',form='FORMATTED')
151  open(51,file='ELFLINES.OUT',form='FORMATTED')
152  call plot1d(50,51,1,elfmt,elfir)
153  close(50)
154  close(51)
155  write(*,*)
156  write(*,'("Info(elfplot):")')
157  write(*,'(" 1D ELF plot written to ELF1D.OUT")')
158  write(*,'(" vertex location lines written to ELFLINES.OUT")')
159case(52)
160  open(50,file='ELF2D.OUT',form='FORMATTED')
161  call plot2d(.false.,50,1,elfmt,elfir)
162  close(50)
163  write(*,*)
164  write(*,'("Info(elfplot): 2D ELF plot written to ELF2D.OUT")')
165case(53)
166  open(50,file='ELF3D.OUT',form='FORMATTED')
167  call plot3d(50,1,elfmt,elfir)
168  close(50)
169  write(*,*)
170  write(*,'("Info(elfplot): 3D ELF plot written to ELF3D.OUT")')
171end select
172deallocate(gwf2mt,gwf2ir,rfmt1,rfmt2,grfir)
173deallocate(grfmt1,grfmt2,elfmt,elfir)
174deallocate(evecfv,evecsv,zfft1,zfft2)
175end subroutine
176!EOC
177
178