File: eveqnsv.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 (381 lines) | stat: -rw-r--r-- 10,524 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
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381

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

subroutine eveqnsv(ngp,igpig,vgpc,apwalm,evalfv,evecfv,evalsvp,evecsv)
use modmain
use moddftu
use modomp
implicit none
! arguments
integer, intent(in) :: ngp,igpig(ngkmax)
real(8), intent(in) :: vgpc(3,ngkmax)
complex(8), intent(in) :: apwalm(ngkmax,apwordmax,lmmaxapw,natmtot)
real(8), intent(in) :: evalfv(nstfv)
complex(8), intent(in) :: evecfv(nmatmax,nstfv)
real(8), intent(out) :: evalsvp(nstsv)
complex(8), intent(out) :: evecsv(nstsv,nstsv)
! local variables
logical socz
integer nsc,nsd,ld,ist,jst
integer ispn,jspn,is,ias
integer nrc,nrci,nrco,irco,irc
integer l,lm,nm,npc,npci
integer igp,i,j,k,nthd
real(8) ca,t1
real(8) ts0,ts1
complex(8) zsum,z1
! automatic arrays
complex(8) zlflm(lmmaxo,3)
! allocatable arrays
complex(8), allocatable :: wfmt1(:,:),wfmt2(:),wfmt3(:),wfmt4(:,:)
complex(8), allocatable :: gwfmt(:,:,:),gzfmt(:,:)
complex(8), allocatable :: wfir1(:),wfir2(:),wfgp(:,:),gwfgp(:,:)
! external functions
complex(8) zdotc,zfmtinp
external zdotc,zfmtinp
! no calculation of second-variational eigenvectors
if (.not.tevecsv) then
  do i=1,nstsv
    evalsvp(i)=evalfv(i)
  end do
  evecsv(:,:)=0.d0
  do i=1,nstsv
    evecsv(i,i)=1.d0
  end do
  return
end if
call timesec(ts0)
! coupling constant of the external A-field (1/c)
ca=1.d0/solsc
! number of spin combinations after application of Hamiltonian
if (spinpol) then
  if (ncmag.or.spinorb) then
    nsc=3
  else
    nsc=2
  end if
  nsd=2
else
  nsc=1
  nsd=1
end if
! special case of spin-orbit coupling and collinear magnetism
if (spinorb.and.cmagz) then
  socz=.true.
else
  socz=.false.
end if
ld=lmmaxdm*nspinor
! zero the second-variational Hamiltonian (stored in the eigenvector array)
evecsv(:,:)=0.d0
!-------------------------!
!     muffin-tin part     !
!-------------------------!
allocate(wfmt1(npcmtmax,nstfv),wfmt2(npcmtmax))
allocate(wfmt3(npcmtmax),wfmt4(npcmtmax,nsc))
if (tafield.or.(xcgrad.eq.4)) allocate(gzfmt(npcmtmax,3))
if (xcgrad.eq.4) allocate(gwfmt(npcmtmax,3,nstfv))
! begin loop over atoms
do ias=1,natmtot
  is=idxis(ias)
  nrc=nrcmt(is)
  nrci=nrcmti(is)
  nrco=nrc-nrci
  irco=nrci+1
  npc=npcmt(is)
  npci=npcmti(is)
! compute the first-variational wavefunctions
  do ist=1,nstfv
    call wavefmt(lradstp,ias,ngp,apwalm(:,:,:,ias),evecfv(:,ist),wfmt1(:,ist))
  end do
! begin loop over states
  do jst=1,nstfv
    if (spinpol) then
! convert wavefunction to spherical coordinates
      call zbsht(nrc,nrci,wfmt1(:,jst),wfmt2)
! apply Kohn-Sham effective magnetic field
      wfmt3(1:npc)=bsmt(1:npc,ias,ndmag)*wfmt2(1:npc)
! convert to spherical harmonics and store in wfmt4
      call zfsht(nrc,nrci,wfmt3,wfmt4(:,1))
      wfmt4(1:npc,2)=-wfmt4(1:npc,1)
! non-collinear magnetic field
      if (ncmag) then
        wfmt3(1:npc)=cmplx(bsmt(1:npc,ias,1),-bsmt(1:npc,ias,2),8)*wfmt2(1:npc)
        call zfsht(nrc,nrci,wfmt3,wfmt4(:,3))
      end if
      if (socz) wfmt4(1:npc,3)=0.d0
! apply spin-orbit coupling if required
      if (spinorb) then
! inner part of muffin-tin
        i=1
        do irc=1,nrci
          call lopzflm(lmaxi,wfmt1(i,jst),lmmaxo,zlflm)
          t1=socfr(irc,ias)
          do j=1,lmmaxi
            wfmt4(i,1)=wfmt4(i,1)+t1*zlflm(j,3)
            wfmt4(i,2)=wfmt4(i,2)-t1*zlflm(j,3)
            wfmt4(i,3)=wfmt4(i,3)+t1*(zlflm(j,1) &
             +cmplx(aimag(zlflm(j,2)),-dble(zlflm(j,2)),8))
            i=i+1
          end do
        end do
! outer part of muffin-tin
        do irc=irco,nrc
          call lopzflm(lmaxo,wfmt1(i,jst),lmmaxo,zlflm)
          t1=socfr(irc,ias)
          do j=1,lmmaxo
            wfmt4(i,1)=wfmt4(i,1)+t1*zlflm(j,3)
            wfmt4(i,2)=wfmt4(i,2)-t1*zlflm(j,3)
            wfmt4(i,3)=wfmt4(i,3)+t1*(zlflm(j,1) &
             +cmplx(aimag(zlflm(j,2)),-dble(zlflm(j,2)),8))
            i=i+1
          end do
        end do
      end if
    else
      do k=1,nsc
        wfmt4(1:npc,k)=0.d0
      end do
    end if
! 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,nsc
            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),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),lmmaxo,zone,wfmt4(i,k),lmmaxo)
          end do
        end if
      end do
    end if
! apply vector potential if required
    if (tafield) then
      call gradzfmt(nrc,nrci,rcmt(:,is),wfmt1(:,jst),npcmtmax,gzfmt)
      do i=1,npc
        z1=afieldc(1)*gzfmt(i,1)+afieldc(2)*gzfmt(i,2)+afieldc(3)*gzfmt(i,3)
        z1=ca*cmplx(-aimag(z1),dble(z1),8)
        wfmt4(i,1:nsd)=wfmt4(i,1:nsd)+z1
      end do
    end if
! second-variational Hamiltonian matrix
    call omp_hold(nstfv,nthd)
!$OMP PARALLEL DEFAULT(SHARED) &
!$OMP PRIVATE(i,j,k) &
!$OMP NUM_THREADS(nthd)
!$OMP DO
    do ist=1,nstfv
      do k=1,nsc
        if (k.eq.1) then
          i=ist
          j=jst
        else if (k.eq.2) then
          i=ist+nstfv
          j=jst+nstfv
        else
          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),wfmt4(:,k))
        end if
      end do
    end do
!$OMP END DO
!$OMP END PARALLEL
    call omp_free(nthd)
! end loop over states
  end do
! apply tau-DFT potential if required
  if (xcgrad.eq.4) then
    do ist=1,nstfv
      call gradzfmt(nrc,nrci,rcmt(:,is),wfmt1(:,ist),npcmtmax,gwfmt(:,:,ist))
    end do
    do jst=1,nstfv
      do k=1,3
        call zbsht(nrc,nrci,gwfmt(:,k,jst),wfmt2)
        wfmt2(1:npc)=wsmt(1:npc,ias)*wfmt2(1:npc)
        call zfsht(nrc,nrci,wfmt2,gzfmt(:,k))
      end do
      do ist=1,nstfv
        zsum=0.d0
        do k=1,3
          z1=zfmtinp(nrc,nrci,rcmt(:,is),r2cmt(:,is),gwfmt(:,k,ist),gzfmt(:,k))
          zsum=zsum+z1
        end do
        zsum=0.5d0*zsum
        evecsv(ist,jst)=evecsv(ist,jst)+zsum
        if (spinpol) then
          i=ist+nstfv
          j=jst+nstfv
          evecsv(i,j)=evecsv(i,j)+zsum
        end if
      end do
    end do
  end if
! end loop over atoms
end do
deallocate(wfmt1,wfmt2,wfmt3,wfmt4)
if (tafield.or.(xcgrad.eq.4)) deallocate(gzfmt)
if (xcgrad.eq.4) deallocate(gwfmt)
!---------------------------!
!     interstitial part     !
!---------------------------!
if (spinpol.or.tafield.or.(xcgrad.eq.4)) then
  if (socz) nsc=2
  allocate(wfir1(ngtot),wfir2(ngtot),wfgp(ngkmax,nsc))
  if (xcgrad.eq.4) allocate(gwfgp(ngkmax,nstfv))
! begin loop over states
  do jst=1,nstfv
    wfir1(:)=0.d0
    do igp=1,ngp
      wfir1(igfft(igpig(igp)))=evecfv(igp,jst)
    end do
! Fourier transform wavefunction to real-space
    call zfftifc(3,ngridg,1,wfir1)
! multiply with magnetic field and transform to G-space
    if (spinpol) then
      wfir2(:)=bsir(:,ndmag)*wfir1(:)
      call zfftifc(3,ngridg,-1,wfir2)
      do igp=1,ngp
        wfgp(igp,1)=wfir2(igfft(igpig(igp)))
      end do
      wfgp(1:ngp,2)=-wfgp(1:ngp,1)
      if (ncmag) then
        wfir2(:)=cmplx(bsir(:,1),-bsir(:,2),8)*wfir1(:)
        call zfftifc(3,ngridg,-1,wfir2)
        do igp=1,ngp
          wfgp(igp,3)=wfir2(igfft(igpig(igp)))
        end do
      end if
    else
      wfgp(1:ngp,1:nsd)=0.d0
    end if
! apply vector potential if required
    if (tafield) then
      wfir1(:)=0.d0
      do igp=1,ngp
        t1=-ca*dot_product(afieldc(:),vgpc(:,igp))
        wfir1(igfft(igpig(igp)))=t1*evecfv(igp,jst)
      end do
      call zfftifc(3,ngridg,1,wfir1)
      wfir1(:)=wfir1(:)*cfunir(:)
      call zfftifc(3,ngridg,-1,wfir1)
      do igp=1,ngp
        z1=wfir1(igfft(igpig(igp)))
        wfgp(igp,1:nsd)=wfgp(igp,1:nsd)+z1
      end do
    end if
! add to Hamiltonian matrix
    call omp_hold(nstfv,nthd)
!$OMP PARALLEL DEFAULT(SHARED) &
!$OMP PRIVATE(i,j,k) &
!$OMP NUM_THREADS(nthd)
!$OMP DO
    do ist=1,nstfv
      do k=1,nsc
        if (k.eq.1) then
          i=ist
          j=jst
        else if (k.eq.2) then
          i=ist+nstfv
          j=jst+nstfv
        else
          i=ist
          j=jst+nstfv
        end if
        if (i.le.j) then
          evecsv(i,j)=evecsv(i,j)+zdotc(ngp,evecfv(:,ist),1,wfgp(:,k),1)
        end if
      end do
    end do
!$OMP END DO
!$OMP END PARALLEL
    call omp_free(nthd)
! end loop over states
  end do
! apply tau-DFT potential if required
  if (xcgrad.eq.4) then
    do k=1,3
      do ist=1,nstfv
        do igp=1,ngp
          z1=evecfv(igp,ist)
          gwfgp(igp,ist)=vgpc(k,igp)*cmplx(-aimag(z1),dble(z1),8)
        end do
      end do
      do jst=1,nstfv
        wfir1(:)=0.d0
        do igp=1,ngp
          wfir1(igfft(igpig(igp)))=gwfgp(igp,jst)
        end do
        call zfftifc(3,ngridg,1,wfir1)
        wfir1(:)=wsir(:)*wfir1(:)
        call zfftifc(3,ngridg,-1,wfir1)
        do igp=1,ngp
          wfgp(igp,1)=wfir1(igfft(igpig(igp)))
        end do
        do ist=1,nstfv
          z1=0.5d0*zdotc(ngp,gwfgp(:,ist),1,wfgp,1)
          evecsv(ist,jst)=evecsv(ist,jst)+z1
          if (spinpol) then
            i=ist+nstfv
            j=jst+nstfv
            evecsv(i,j)=evecsv(i,j)+z1
          end if
        end do
      end do
    end do
  end if
  deallocate(wfir1,wfir2,wfgp)
  if (xcgrad.eq.4) deallocate(gwfgp)
end if
! 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)
  end do
end do
if (spcpl.or.(.not.spinpol)) then
! spins are coupled; or spin-unpolarised: full diagonalisation
  call eveqnz(nstsv,nstsv,evecsv,evalsvp)
else
! spins not coupled: block diagonalise H
  call eveqnz(nstfv,nstsv,evecsv,evalsvp)
  i=nstfv+1
  call eveqnz(nstfv,nstsv,evecsv(i,i),evalsvp(i))
  do i=1,nstfv
    do j=1,nstfv
      evecsv(i,j+nstfv)=0.d0
      evecsv(i+nstfv,j)=0.d0
    end do
  end do
end if
call timesec(ts1)
!$OMP CRITICAL(eveqnsv_)
timesv=timesv+ts1-ts0
!$OMP END CRITICAL(eveqnsv_)
return
end subroutine