File: writelsj.f90

package info (click to toggle)
elkcode 5.4.24-2
  • links: PTS, VCS
  • area: main
  • in suites: buster, sid
  • size: 12,840 kB
  • sloc: f90: 48,415; fortran: 22,457; perl: 965; makefile: 384; sh: 369; python: 105; ansic: 67
file content (182 lines) | stat: -rw-r--r-- 5,884 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182

! Copyright (C) 2002-2007 J. K. Dewhurst, S. Sharma, C. Ambrosch-Draxl and
! F. Cricchio. This file is distributed under the terms of the GNU General
! Public License. See the file COPYING for license details.

subroutine writelsj
use modmain
implicit none
! local variables
integer kst,ik,ist,lm
integer ispn,is,ia,ias
real(8) xl(3),xs(3),t1
! allocatable arrays
complex(8), allocatable :: apwalm(:,:,:,:,:)
complex(8), allocatable :: evecfv(:,:,:),evecsv(:,:)
complex(8), allocatable :: dmat1(:,:,:,:,:),dmat2(:,:,:,:,:)
complex(8), allocatable :: zlflm(:,:)
! initialise universal variables
call init0
call init1
allocate(apwalm(ngkmax,apwordmax,lmmaxapw,natmtot,nspnfv))
allocate(evecfv(nmatmax,nstfv,nspnfv),evecsv(nstsv,nstsv))
allocate(dmat1(lmmaxo,nspinor,lmmaxo,nspinor,natmtot))
allocate(dmat2(lmmaxo,nspinor,lmmaxo,nspinor,nstsv))
allocate(zlflm(lmmaxo,3))
! read density and potentials from file
call readstate
! read Fermi energy from file
call readfermi
! find the new linearisation energies
call linengy
! generate the APW radial functions
call genapwfr
! generate the local-orbital radial functions
call genlofr
if (task.eq.15) then
! compute total L, S and J
  dmat1(:,:,:,:,:)=0.d0
  do ik=1,nkpt
! get the eigenvectors and occupancies from file
    call getevecfv(filext,ik,vkl(:,ik),vgkl(:,:,:,ik),evecfv)
    call getevecsv(filext,ik,vkl(:,ik),evecsv)
    call getoccsv(filext,ik,vkl(:,ik),occsv(:,ik))
! find the matching coefficients
    do ispn=1,nspnfv
      call match(ngk(ispn,ik),gkc(:,ispn,ik),tpgkc(:,:,ispn,ik), &
       sfacgk(:,:,ispn,ik),apwalm(:,:,:,:,ispn))
    end do
! loop over species and atoms
    do is=1,nspecies
      do ia=1,natoms(is)
        ias=idxas(ia,is)
! generate the density matrix
        call gendmat(.false.,.false.,0,lmaxo,ias,ngk(:,ik),apwalm,evecfv, &
         evecsv,lmmaxo,dmat2)
        do ist=1,nstsv
          t1=wkpt(ik)*occsv(ist,ik)
          dmat1(:,:,:,:,ias)=dmat1(:,:,:,:,ias)+t1*dmat2(:,:,:,:,ist)
        end do
      end do
    end do
! end loop over k-points
  end do
! symmetrise the density matrix
  call symdmat(lmaxo,lmmaxo,dmat1)
  open(50,file='LSJ.OUT',form='FORMATTED')
  write(50,*)
  write(50,'("Expectation values are computed only over the muffin-tin")')
! loop over species and atoms
  do is=1,nspecies
    write(50,*)
    write(50,'("Species : ",I4," (",A,")")') is,trim(spsymb(is))
    do ia=1,natoms(is)
      ias=idxas(ia,is)
! compute tr(LD)
      xl(:)=0.d0
      do ispn=1,nspinor
        do lm=1,lmmaxo
          call lopzflm(lmaxo,dmat1(:,ispn,lm,ispn,ias),lmmaxo,zlflm)
          xl(:)=xl(:)+dble(zlflm(lm,:))
        end do
      end do
! compute tr(sigma D)
      xs(:)=0.d0
      if (spinpol) then
        do lm=1,lmmaxo
          xs(1)=xs(1)+dble(dmat1(lm,2,lm,1,ias)+dmat1(lm,1,lm,2,ias))
          xs(2)=xs(2)+dble(-zi*dmat1(lm,2,lm,1,ias)+zi*dmat1(lm,1,lm,2,ias))
          xs(3)=xs(3)+dble(dmat1(lm,1,lm,1,ias)-dmat1(lm,2,lm,2,ias))
        end do
      end if
! S = 1/2 sigma
      xs(:)=0.5d0*xs(:)
      write(50,'(" atom : ",I4)') ia
      write(50,'("  L : ",3G18.10)') xl(:)
      write(50,'("  S : ",3G18.10)') xs(:)
      write(50,'("  J : ",3G18.10)') xl(:)+xs(:)
! end loop over atoms and species
    end do
  end do
  close(50)
  write(*,*)
  write(*,'("Info(writelsj):")')
  write(*,'(" total L, S and J expectation values written to LSJ.OUT")')
else
! compute L, S and J for all states in kstlist
  open(50,file='LSJ_KST.OUT',form='FORMATTED')
  write(50,*)
  write(50,'("Expectation values are computed only over the muffin-tin")')
  do kst=1,nkstlist
    ik=kstlist(1,kst)
    ist=kstlist(2,kst)
    if ((ik.le.0).or.(ik.gt.nkpt)) then
      write(*,*)
      write(*,'("Error(writelsj): k-point out of range : ",I8)') ik
      write(*,*)
      stop
    end if
    if ((ist.le.0).or.(ist.gt.nstsv)) then
      write(*,*)
      write(*,'("Error(writelsj): state out of range : ",I8)') ist
      write(*,*)
      stop
    end if
! get the eigenvectors and occupancies from file
    call getevecfv(filext,ik,vkl(:,ik),vgkl(:,:,:,ik),evecfv)
    call getevecsv(filext,ik,vkl(:,ik),evecsv)
    call getoccsv(filext,ik,vkl(:,ik),occsv(:,ik))
! find the matching coefficients
    do ispn=1,nspnfv
      call match(ngk(ispn,ik),gkc(:,ispn,ik),tpgkc(:,:,ispn,ik), &
       sfacgk(:,:,ispn,ik),apwalm(:,:,:,:,ispn))
    end do
! loop over species and atoms
    do is=1,nspecies
      do ia=1,natoms(is)
        ias=idxas(ia,is)
! generate the density matrix
        call gendmat(.false.,.false.,0,lmaxo,ias,ngk(:,ik),apwalm,evecfv, &
         evecsv,lmmaxo,dmat2)
! compute tr(LD)
        xl(:)=0.d0
        do ispn=1,nspinor
          do lm=1,lmmaxo
            call lopzflm(lmaxo,dmat2(:,ispn,lm,ispn,ist),lmmaxo,zlflm)
            xl(:)=xl(:)+dble(zlflm(lm,:))
          end do
        end do
! compute tr(sigma D)
        xs(:)=0.d0
        if (spinpol) then
          do lm=1,lmmaxo
            xs(1)=xs(1)+dble(dmat2(lm,2,lm,1,ist)+dmat2(lm,1,lm,2,ist))
            xs(2)=xs(2)+dble(-zi*dmat2(lm,2,lm,1,ist)+zi*dmat2(lm,1,lm,2,ist))
            xs(3)=xs(3)+dble(dmat2(lm,1,lm,1,ist)-dmat2(lm,2,lm,2,ist))
          end do
        else
          xs(3)=1.d0
        end if
! S = 1/2 sigma
        xs(:)=0.5d0*xs(:)
        write(50,*)
        write(50,'("k-point : ",I6,3G18.10)') ik,vkl(:,ik)
        write(50,'("state : ",I6)') ist
        write(50,'("species : ",I4," (",A,"), atom : ",I4)') is, &
         trim(spsymb(is)),ia
        write(50,'(" L : ",3G18.10)') xl(:)
        write(50,'(" S : ",3G18.10)') xs(:)
        write(50,'(" J : ",3G18.10)') xl(:)+xs(:)
      end do
    end do
  end do
  close(50)
  write(*,*)
  write(*,'("Info(writelsj):")')
  write(*,'(" L, S and J expectation values for each k-point and state")')
  write(*,'("  in kstlist written to LSJ_KST.OUT")')
end if
deallocate(apwalm,evecfv,evecsv,dmat1,dmat2,zlflm)
return
end subroutine