src/writelsj.f90

   1
   2 ! Copyright (C) 2002-2007 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 subroutine writelsj
   7 use modmain
   8 implicit none
   9 ! local variables
  10 integer ik,ispn,jspn,i,j,kst
  11 integer is,ia,ias,ist,irc,n
  12 real(8) xl(3),xs(3),xj(3),wo
  13 complex(8) zt1
  14 ! allocatable arrays
  15 logical, allocatable :: done(:,:)
  16 real(8), allocatable :: xlt(:,:),xst(:,:),xjt(:,:)
  17 complex(8), allocatable :: apwalm(:,:,:,:,:)
  18 complex(8), allocatable :: evecfv(:,:,:)
  19 complex(8), allocatable :: evecsv(:,:)
  20 complex(8), allocatable :: wfmt1(:,:,:,:)
  21 complex(8), allocatable :: wfmt2(:,:,:)
  22 complex(8), allocatable :: wfmt3(:,:,:,:)
  23 complex(8), allocatable :: wfmt4(:,:,:,:)
  24 complex(8), allocatable :: zlflm(:,:)
  25 ! external functions
  26 complex(8) zfmtinp
  27 external zfmtinp
  28 ! initialise universal variables
  29 call init0
  30 call init1
  31 allocate(done(nstfv,nspnfv))
  32 allocate(xlt(3,natmtot),xst(3,natmtot),xjt(3,natmtot))
  33 allocate(apwalm(ngkmax,apwordmax,lmmaxapw,natmtot,nspnfv))
  34 allocate(evecfv(nmatmax,nstfv,nspnfv))
  35 allocate(evecsv(nstsv,nstsv))
  36 allocate(wfmt1(lmmaxapw,nrcmtmax,nstfv,nspnfv))
  37 allocate(wfmt2(lmmaxapw,nrcmtmax,2))
  38 allocate(wfmt3(lmmaxapw,nrcmtmax,2,3))
  39 allocate(wfmt4(lmmaxapw,nrcmtmax,2,3))
  40 allocate(zlflm(lmmaxapw,3))
  41 ! read density and potentials from file
  42 call readstate
  43 ! find the new linearisation energies
  44 call linengy
  45 ! generate the APW radial functions
  46 call genapwfr
  47 ! generate the local-orbital radial functions
  48 call genlofr
  49 if (task.eq.15) then
  50 ! get the occupancies from file
  51   do ik=1,nkpt
  52     call getoccsv(vkl(1,ik),occsv(1,ik))
  53   end do
  54 else
  55 ! flag the occupancies for the k-points and states in kstlist
  56   occsv(:,:)=0.d0
  57   do kst=1,nkstlist
  58     ik=kstlist(1,kst)
  59     j=kstlist(2,kst)
  60     if ((ik.le.0).or.(ik.gt.nkpt)) then
  61       write(*,*)
  62       write(*,'("Error(writelsj): k-point out of range : ",I8)') ik
  63       write(*,*)
  64       stop
  65     end if
  66     if ((j.le.0).or.(j.gt.nstsv)) then
  67       write(*,*)
  68       write(*,'("Error(writelsj): state out of range : ",I8)') j
  69       write(*,*)
  70       stop
  71     end if
  72     occsv(j,ik)=1.d0
  73   end do
  74 end if
  75 if (task.eq.15) then
  76   open(50,file='LSJ.OUT',action='WRITE',form='FORMATTED')
  77 else
  78   open(50,file='LSJ_KST.OUT',action='WRITE',form='FORMATTED')
  79 end if
  80 write(50,*)
  81 write(50,'("Expectation values are computed only over the muffin-tin")')
  82 ! zero the total expection values
  83 xlt(:,:)=0.d0
  84 xst(:,:)=0.d0
  85 xjt(:,:)=0.d0
  86 ! begin loop over k-point and state list
  87 do ik=1,nkpt
  88 ! get the eigenvectors from file
  89   call getevecfv(vkl(1,ik),vgkl(1,1,ik,1),evecfv)
  90   call getevecsv(vkl(1,ik),evecsv)
  91 ! find the matching coefficients
  92   do ispn=1,nspnfv
  93     call match(ngk(ik,ispn),gkc(1,ik,ispn),tpgkc(1,1,ik,ispn), &
  94      sfacgk(1,1,ik,ispn),apwalm(1,1,1,1,ispn))
  95   end do
  96 ! begin loop over atoms and species
  97   do is=1,nspecies
  98     n=lmmaxapw*nrcmt(is)
  99     do ia=1,natoms(is)
 100       ias=idxas(ia,is)
 101       done(:,:)=.false.
 102 ! begin loop over second-variational states
 103       do j=1,nstsv
 104         wo=wkpt(ik)*occsv(j,ik)
 105         if (abs(wo).gt.epsocc) then
 106           if (tevecsv) then
 107 ! generate spinor wavefunction from second-variational eigenvectors
 108             wfmt2(:,:,:)=0.d0
 109             i=0
 110             do ispn=1,nspinor
 111               if (spinsprl) then
 112                 jspn=ispn
 113               else
 114                 jspn=1
 115               end if
 116               do ist=1,nstfv
 117                 i=i+1
 118                 zt1=evecsv(i,j)
 119                 if (abs(dble(zt1))+abs(aimag(zt1)).gt.epsocc) then
 120                   if (.not.done(ist,jspn)) then
 121                     call wavefmt(lradstp,lmaxapw,is,ia,ngk(ik,jspn), &
 122                      apwalm(1,1,1,1,jspn),evecfv(1,ist,jspn),lmmaxapw, &
 123                      wfmt1(1,1,ist,jspn))
 124                     done(ist,jspn)=.true.
 125                   end if
 126 ! add to spinor wavefunction
 127                   call zaxpy(n,zt1,wfmt1(1,1,ist,jspn),1,wfmt2(1,1,ispn),1)
 128                 end if
 129               end do
 130             end do
 131           else
 132 ! spin-unpolarised wavefunction
 133             call wavefmt(lradstp,lmaxapw,is,ia,ngk(ik,1),apwalm,evecfv(1,j,1), &
 134              lmmaxapw,wfmt2)
 135           end if
 136 ! apply the L operator
 137           do ispn=1,nspinor
 138             do irc=1,nrcmt(is)
 139               call lopzflm(lmaxapw,wfmt2(1,irc,ispn),lmmaxapw,zlflm)
 140               do i=1,3
 141                 wfmt3(:,irc,ispn,i)=zlflm(:,i)
 142               end do
 143             end do
 144           end do
 145 ! compute the expectation value of L
 146           xl(:)=0.d0
 147           do i=1,3
 148             do ispn=1,nspinor
 149               zt1=zfmtinp(.true.,lmaxapw,nrcmt(is),rcmt(1,is),lmmaxapw, &
 150                wfmt2(1,1,ispn),wfmt3(1,1,ispn,i))
 151               xl(i)=xl(i)+dble(zt1)
 152             end do
 153           end do
 154 ! compute the expecations value of S
 155           if (spinpol) then
 156 ! apply the S operator
 157             do irc=1,nrcmt(is)
 158 ! S_x
 159               wfmt4(:,irc,1,1)=0.5d0*wfmt2(:,irc,2)
 160               wfmt4(:,irc,2,1)=0.5d0*wfmt2(:,irc,1)
 161 ! S_y
 162               wfmt4(:,irc,1,2)=-0.5d0*zi*wfmt2(:,irc,2)
 163               wfmt4(:,irc,2,2)=0.5d0*zi*wfmt2(:,irc,1)
 164 ! S_z
 165               wfmt4(:,irc,1,3)=0.5d0*wfmt2(:,irc,1)
 166               wfmt4(:,irc,2,3)=-0.5d0*wfmt2(:,irc,2)
 167             end do
 168             xs(:)=0.d0
 169             do i=1,3
 170               do ispn=1,nspinor
 171                 zt1=zfmtinp(.true.,lmaxapw,nrcmt(is),rcmt(1,is),lmmaxapw, &
 172                  wfmt2(1,1,ispn),wfmt4(1,1,ispn,i))
 173                 xs(i)=xs(i)+dble(zt1)
 174               end do
 175             end do
 176           else
 177 ! spin-unpolarised case
 178             xs(:)=0.d0
 179           end if
 180 ! expectation value of J
 181           xj(:)=xl(:)+xs(:)
 182 ! add to the total
 183           xlt(:,ias)=xlt(:,ias)+wo*xl(:)
 184           xst(:,ias)=xst(:,ias)+wo*xs(:)
 185           xjt(:,ias)=xjt(:,ias)+wo*xj(:)
 186 ! write the particular L, S and J to file
 187           if (task.eq.16) then
 188             write(50,*)
 189             write(50,'("k-point : ",I6,3G18.10)') ik,vkl(:,ik)
 190             write(50,'("state : ",I6)') j
 191             write(50,'("species : ",I4," (",A,"), atom : ",I4)') is, &
 192              trim(spsymb(is)),ia
 193             write(50,'(" L : ",3G18.10)') xl
 194             write(50,'(" S : ",3G18.10)') xs
 195             write(50,'(" J : ",3G18.10)') xj
 196           end if
 197         end if
 198 ! end loop over states
 199       end do
 200 ! end loops over atoms and species
 201     end do
 202   end do
 203 ! end loop over k-points
 204 end do
 205 ! write the total L, S and J to file if required
 206 if (task.eq.15) then
 207   do is=1,nspecies
 208     do ia=1,natoms(is)
 209       ias=idxas(ia,is)
 210       write(50,*)
 211       write(50,'("Species : ",I4," (",A,"), atom : ",I4)') is, &
 212        trim(spsymb(is)),ia
 213       write(50,'(" L : ",3G18.10)') xlt(:,ias)
 214       write(50,'(" S : ",3G18.10)') xst(:,ias)
 215       write(50,'(" J : ",3G18.10)') xjt(:,ias)
 216     end do
 217   end do
 218 end if
 219 close(50)
 220 write(*,*)
 221 write(*,'("Info(writelsj):")')
 222 if (task.eq.15) then
 223   write(*,'(" total L, S and J expectation values written to LSJ.OUT")')
 224 else
 225   write(*,'(" L, S and J expectation values for each k-point and state")')
 226   write(*,'("  in kstlist written to LSJ_KST.OUT")')
 227 end if
 228 write(*,*)
 229 deallocate(done,xlt,xst,xjt)
 230 deallocate(apwalm,evecfv,evecsv)
 231 deallocate(wfmt1,wfmt2,wfmt3,wfmt4,zlflm)
 232 return
 233 end subroutine
 234