File: eveqnss.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 (220 lines) | stat: -rw-r--r-- 6,059 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
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220

! Copyright (C) 2006 F. Bultmark, F. Cricchio, L. Nordstrom and J. K. Dewhurst.
! This file is distributed under the terms of the GNU General Public License.
! See the file COPYING for license details.

subroutine eveqnss(ngp,igpig,apwalm,evalfv,evecfv,evalsvp,evecsv)
use modmain
use moddftu
use modomp
implicit none
! arguments
integer, intent(in) :: ngp(nspnfv),igpig(ngkmax,nspnfv)
complex(8), intent(in) :: apwalm(ngkmax,apwordmax,lmmaxapw,natmtot,nspnfv)
real(8), intent(in) :: evalfv(nstfv,nspnfv)
complex(8), intent(in) :: evecfv(nmatmax,nstfv,nspnfv)
real(8), intent(out) :: evalsvp(nstsv)
complex(8), intent(out) :: evecsv(nstsv,nstsv)
! local variables
integer ist,jst,ispn,jspn
integer is,ia,ias,i,j,k
integer nrc,nrci,nrco
integer l,lm,nm,npc,npci
integer igp,ld,nthd
real(8) t1
real(8) ts0,ts1
complex(8) zq
! allocatable arrays
complex(8), allocatable :: wfmt1(:,:,:),wfmt2(:,:),wfmt3(:),wfmt4(:,:)
complex(8), allocatable :: wfir1(:,:),wfir2(:),wfgp(:,:)
! external functions
complex(8) zdotc,zfmtinp
external zdotc,zfmtinp
if (.not.spinpol) then
  write(*,*)
  write(*,'("Error(eveqnss): spin-unpolarised calculation")')
  write(*,*)
  stop
end if
call timesec(ts0)
ld=lmmaxdm*nspinor
! zero the second-variational Hamiltonian (stored in the eigenvector array)
evecsv(:,:)=0.d0
!-------------------------!
!     muffin-tin part     !
!-------------------------!
allocate(wfmt1(npcmtmax,nstfv,nspnfv))
allocate(wfmt2(npcmtmax,nspnfv))
allocate(wfmt3(npcmtmax),wfmt4(npcmtmax,3))
do ias=1,natmtot
  is=idxis(ias)
  ia=idxia(ias)
  nrc=nrcmt(is)
  nrci=nrcmti(is)
  nrco=nrc-nrci
  npc=npcmt(is)
  npci=npcmti(is)
! de-phasing factor (FC, FB & LN)
  t1=-0.5d0*dot_product(vqcss(:),atposc(:,ia,is))
  zq=cmplx(cos(t1),sin(t1),8)
! compute the first-variational wavefunctions
  do ispn=1,nspnfv
    if (ispn.eq.2) zq=conjg(zq)
    do ist=1,nstfv
      call wavefmt(lradstp,ias,ngp(ispn),apwalm(:,:,:,ias,ispn), &
       evecfv(:,ist,ispn),wfmt1(:,ist,ispn))
! de-phase if required
      if (ssdph) wfmt1(1:npc,ist,ispn)=zq*wfmt1(1:npc,ist,ispn)
    end do
  end do
  do jst=1,nstfv
! convert wavefunction to spherical coordinates
    do ispn=1,nspnfv
      call zbsht(nrc,nrci,wfmt1(:,jst,ispn),wfmt2(:,ispn))
    end do
! apply effective magnetic field and convert to spherical harmonics
    wfmt3(1:npc)=bsmt(1:npc,ias,3)*wfmt2(1:npc,1)
    call zfsht(nrc,nrci,wfmt3,wfmt4(:,1))
    wfmt3(1:npc)=-bsmt(1:npc,ias,3)*wfmt2(1:npc,2)
    call zfsht(nrc,nrci,wfmt3,wfmt4(:,2))
    wfmt3(1:npc)=cmplx(bsmt(1:npc,ias,1),-bsmt(1:npc,ias,2),8)*wfmt2(1:npc,2)
    call zfsht(nrc,nrci,wfmt3,wfmt4(:,3))
! apply muffin-tin potential matrix if required
    if (tvmatmt) then
      do l=0,lmaxdm
        if (tvmmt(l,ias)) then
          nm=2*l+1
          lm=idxlm(l,-l)
          do k=1,3
            if (k.eq.1) then
              ispn=1
              jspn=1
            else if (k.eq.2) then
              ispn=2
              jspn=2
            else
              ispn=1
              jspn=2
            end if
            if (l.le.lmaxi) then
              call zgemm('N','N',nm,nrci,nm,zone,vmatmt(lm,ispn,lm,jspn,ias), &
               ld,wfmt1(lm,jst,jspn),lmmaxi,zone,wfmt4(lm,k),lmmaxi)
            end if
            i=npci+lm
            call zgemm('N','N',nm,nrco,nm,zone,vmatmt(lm,ispn,lm,jspn,ias),ld, &
             wfmt1(i,jst,jspn),lmmaxo,zone,wfmt4(i,k),lmmaxo)
          end do
        end if
      end do
    end if
! second-variational Hamiltonian matrix
    call omp_hold(nstfv,nthd)
!$OMP PARALLEL DEFAULT(SHARED) &
!$OMP PRIVATE(i,j,k,ispn) &
!$OMP NUM_THREADS(nthd)
!$OMP DO
    do ist=1,nstfv
      do k=1,3
        if (k.eq.1) then
          ispn=1
          i=ist
          j=jst
        else if (k.eq.2) then
          ispn=2
          i=ist+nstfv
          j=jst+nstfv
        else
          ispn=1
          i=ist
          j=jst+nstfv
        end if
        if (i.le.j) then
          evecsv(i,j)=evecsv(i,j)+zfmtinp(nrc,nrci,rcmt(:,is),r2cmt(:,is), &
           wfmt1(:,ist,ispn),wfmt4(:,k))
        end if
      end do
    end do
!$OMP END DO
!$OMP END PARALLEL
    call omp_free(nthd)
  end do
end do
deallocate(wfmt1,wfmt2,wfmt3,wfmt4)
!---------------------------!
!     interstitial part     !
!---------------------------!
allocate(wfir1(ngtot,nspnfv),wfir2(ngtot),wfgp(ngkmax,3))
do jst=1,nstfv
  do ispn=1,nspnfv
    wfir1(:,ispn)=0.d0
    do igp=1,ngp(ispn)
      wfir1(igfft(igpig(igp,ispn)),ispn)=evecfv(igp,jst,ispn)
    end do
! Fourier transform wavefunction to real-space
    call zfftifc(3,ngridg,1,wfir1(:,ispn))
  end do
! multiply with magnetic field and transform to G-space
  wfir2(:)=bsir(:,3)*wfir1(:,1)
  call zfftifc(3,ngridg,-1,wfir2)
  do igp=1,ngp(1)
    wfgp(igp,1)=wfir2(igfft(igpig(igp,1)))
  end do
  wfir2(:)=-bsir(:,3)*wfir1(:,2)
  call zfftifc(3,ngridg,-1,wfir2)
  do igp=1,ngp(2)
    wfgp(igp,2)=wfir2(igfft(igpig(igp,2)))
  end do
  wfir2(:)=cmplx(bsir(:,1),-bsir(:,2),8)*wfir1(:,2)
  call zfftifc(3,ngridg,-1,wfir2)
  do igp=1,ngp(1)
    wfgp(igp,3)=wfir2(igfft(igpig(igp,1)))
  end do
! add to Hamiltonian matrix
  call omp_hold(nstfv,nthd)
!$OMP PARALLEL DEFAULT(SHARED) &
!$OMP PRIVATE(i,j,k,ispn) &
!$OMP NUM_THREADS(nthd)
!$OMP DO
  do ist=1,nstfv
    do k=1,3
      if (k.eq.1) then
        ispn=1
        i=ist
        j=jst
      else if (k.eq.2) then
        ispn=2
        i=ist+nstfv
        j=jst+nstfv
      else
        ispn=1
        i=ist
        j=jst+nstfv
      end if
      if (i.le.j) then
        evecsv(i,j)=evecsv(i,j)+zdotc(ngp(ispn),evecfv(:,ist,ispn),1, &
         wfgp(:,k),1)
      end if
    end do
  end do
!$OMP END DO
!$OMP END PARALLEL
  call omp_free(nthd)
end do
deallocate(wfir1,wfir2,wfgp)
! add the diagonal first-variational part
i=0
do ispn=1,nspinor
  do ist=1,nstfv
    i=i+1
    evecsv(i,i)=evecsv(i,i)+evalfv(ist,ispn)
  end do
end do
! diagonalise the second-variational Hamiltonian
call eveqnz(nstsv,nstsv,evecsv,evalsvp)
call timesec(ts1)
!$OMP CRITICAL(eveqnss_)
timesv=timesv+ts1-ts0
!$OMP END CRITICAL(eveqnss_)
return
end subroutine