IffGit has a new shared runner for building Docker images in GitLab CI. Visit https://iffgit.fz-juelich.de/examples/ci-docker-in-docker for more details.

wann_mmkb_vac.F 11.2 KB
Newer Older
1
2
3
4
5
6
!--------------------------------------------------------------------------------
! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------

7
8
9
10
11
12
13
14
15
16
17
      MODULE m_wann_mmkb_vac
      use m_juDFT
c**************************************************************
c      Determines the overlap matrix Mmn(k) in the vacuum
c      in the film case for the wannier functions.
c      For more details see routine wannier.F 
c
c      Y. Mokrousov, F. Freimuth
c*************************************************************** 
      CONTAINS
      SUBROUTINE wann_mmkb_vac(
18
     >     vacchi,l_noco,nlotot,qss,
19
20
     >     nbnd,z1,nmzd,nv2d,k1d,k2d,k3d,n3d,nvac,
     >     ig,nmz,delz,ig2,area,bmat,
21
22
     >     bbmat,evac,evac_b,bkpt,bkpt_b,vz,vz_b,
     >     nslibd,nslibd_b,jspin,jspin_b,
23
     >     k1,k2,k3,k1_b,k2_b,k3_b,jspd,nvd,
24
     >     nbasfcn,neigd,zMat,zMat_b,nv,nv_b,omtil,gb,
25
26
     <     mmn)

27
      USE m_types
28
      use m_constants
29
      use m_intgr, only : intgz0
30
31
32
      USE m_vacuz
      USE m_vacudz

33
      implicit none
34

35
      TYPE(t_mat), INTENT(IN) :: zMat, zMat_b
36

37
38
c     .. scalar Arguments..
      logical, intent (in) :: l_noco
39
      integer, intent (in) :: nlotot,jspin_b
40
41
42
43
44
45
46
      real,    intent (in) :: qss(3)
      integer, intent (in) :: nmzd,nv2d,k1d,k2d,k3d,n3d,nbnd
      integer, intent (in) :: nmz,nslibd,nslibd_b,nvac
      integer, intent (in) :: jspin,jspd,nvd
      integer, intent (in) :: nbasfcn,neigd
      integer, intent (in) :: gb(3)
      real,    intent (in) :: delz,z1,omtil,area
47
      complex, intent (in) :: vacchi
48
49

c     ..array arguments..
50
      real,    intent (in) :: bkpt(3),bkpt_b(3),evac(2),evac_b(2)
51
52
      integer, intent (in) :: ig(-k1d:k1d,-k2d:k2d,-k3d:k3d)
      integer, intent (in) :: ig2(n3d),nv(jspd),nv_b(jspd)
53
54
      real,    intent (in) :: vz(nmzd,2),vz_b(nmzd,2)
      real,    intent (in) :: bbmat(3,3),bmat(3,3)
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
      integer, intent (in) :: k1(nvd,jspd),k2(nvd,jspd),k3(nvd,jspd)
      integer,intent(in)::k1_b(nvd,jspd),k2_b(nvd,jspd),k3_b(nvd,jspd)
      complex, intent (inout) :: mmn(nbnd,nbnd)

c     ..basis wavefunctions in the vacuum
      complex, allocatable :: ac(:,:),bc(:,:),ac_b(:,:),bc_b(:,:)
      real,    allocatable :: dt(:),dte(:)
      real,    allocatable :: t(:),te(:),tei(:)
      real,    allocatable :: u(:,:),ue(:,:),v(:)
      real,    allocatable :: dt_b(:),dte_b(:)
      real,    allocatable :: t_b(:),te_b(:),tei_b(:)
      real,    allocatable :: u_b(:,:),ue_b(:,:)

c     ..local scalars..
      logical tail
70
      real wronk,arg,zks,tpi,vz0(2),vz0_b(2),scale,evacp,ev,const
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
      real uu,ud,du,dd,xx(nmz),xximag(nmz)
      real :: uuimag,udimag,duimag,ddimag,qss1,qss2
      integer i,m,l,j,k,n,nv2,nv2_b,ivac,n2,n2_b,sign,ik
      integer :: lprime,np1,addnoco,addnoco2
      complex :: av,bv,ic,c_1
      integer, allocatable :: kvac1(:),kvac2(:),map2(:)
      integer, allocatable :: kvac1_b(:),kvac2_b(:),map2_b(:)
      integer symvac,symvacvac

      allocate ( ac(nv2d,nslibd),bc(nv2d,nslibd),
     +           ac_b(nv2d,nslibd_b),bc_b(nv2d,nslibd_b),
     +           dt(nv2d),dte(nv2d),t(nv2d),te(nv2d),
     +           tei(nv2d),u(nmzd,nv2d),ue(nmzd,nv2d),
     +           dt_b(nv2d),dte_b(nv2d),t_b(nv2d),te_b(nv2d),
     +           tei_b(nv2d),u_b(nmzd,nv2d),ue_b(nmzd,nv2d),
     +           v(3),kvac1(nv2d),kvac2(nv2d),map2(nvd),
     +           kvac1_b(nv2d),kvac2_b(nv2d),map2_b(nvd) )

      tpi = 2 * pimach() ; ic = cmplx(0.,1.)

      tail = .true.
      np1 = nmz + 1

c.. determining the indexing array (in-plane stars)
c.. for the k-point

      wronk = 2.0
      const = 1.0 / ( sqrt(omtil)*wronk )

      do ivac = 1,2
         vz0(ivac) = vz(nmz,ivac)
102
         vz0_b(ivac) = vz_b(nmz,ivac)
103
104
105
106
107
108
109
110
      enddo

      n2 = 0 ; n2_b = 0

      addnoco=0
      addnoco2=0
      if(l_noco.and.jspin.eq.2)then
        addnoco=nv(1)+nlotot
111
112
      endif
      if(l_noco.and.jspin_b.eq.2)then
113
114
115
116
117
118
119
120
121
122
123
124
        addnoco2=nv_b(1)+nlotot
      endif

      do 40 k = 1,nv(jspin)
         do 30 j = 1,n2
            if ( k1(k,jspin).eq.kvac1(j) .and.
     +          k2(k,jspin).eq.kvac2(j) ) then
                map2(k) = j
                goto 40
             endif 
 30      continue
         n2 = n2 + 1
125
126
127
128
129
         
         IF(n2>nv2d) then
            write(*,*)n2,nv2d,'jspin',jspin
         endif
 
130
131
132
133
134
135
136
         IF (n2>nv2d)  CALL juDFT_error("wannier Mmn vac",calledby
     +        ="wann_mmkb_vac")

         kvac1(n2) = k1(k,jspin)
         kvac2(n2) = k2(k,jspin)
         map2(k) = n2
 40   continue
137
         !write(*,*)'ok',n2,nv2d,'jspin',jspin
138
139
140

c.. and for the b-point
 
141
      do 41 k = 1,nv_b(jspin_b)
142
         do 31 j = 1,n2_b
143
144
            if ( k1_b(k,jspin_b).eq.kvac1_b(j) .and.
     +          k2_b(k,jspin_b).eq.kvac2_b(j) ) then
145
146
147
148
149
150
                map2_b(k) = j
                goto 41
             endif
 31      continue
         n2_b = n2_b + 1
     
151
152
153
154
         IF(n2_b>nv2d) then
            write(*,*)n2_b,nv2d,'jspin_b',jspin_b
         endif

155
156
157
         IF (n2_b>nv2d)  CALL juDFT_error("wannier Mmn vac",calledby
     +        ="wann_mmkb_vac")

158
159
         kvac1_b(n2_b) = k1_b(k,jspin_b)
         kvac2_b(n2_b) = k2_b(k,jspin_b)
160
161
         map2_b(k) = n2_b
 41   continue
162
         !write(*,*)'ok',n2_b,nv2d,'jspin_b',jspin_b
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

c...cycle by the vacua
      do 140 ivac = 1,nvac



      sign = 3. - 2.*ivac
      evacp = evac(ivac)

      nv2 = n2 ; nv2_b = n2_b

c.. the body of the routine

      qss1=0.0
      qss2=0.0
      if(l_noco.and.jspin.eq.1)then
        qss1=-qss(1)/2.0
        qss2=-qss(2)/2.0
      elseif(l_noco.and.jspin.eq.2)then
        qss1=qss(1)/2.0
        qss2=qss(2)/2.0
      endif

      do ik = 1,nv2
         v(1) = bkpt(1) + kvac1(ik) + qss1
         v(2) = bkpt(2) + kvac2(ik) + qss2
         v(3) = 0.
Daniel Wortmann's avatar
Daniel Wortmann committed
190
         ev = evacp - 0.5*dot_product(v,matmul(bbmat,v))
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
         call vacuz(ev,vz(1,ivac),vz0(ivac),nmz,delz,t(ik),dt(ik),
     +        u(1,ik))
         call vacudz(ev,vz(1,ivac),vz0(ivac),nmz,delz,te(ik),
     +        dte(ik),tei(ik),ue(1,ik),dt(ik),
     +        u(1,ik))
         scale = wronk/ (te(ik)*dt(ik)-dte(ik)*t(ik))
         te(ik) = scale*te(ik)
         dte(ik) = scale*dte(ik)
         tei(ik) = scale*tei(ik)
         do j = 1,nmz
            ue(j,ik) = scale*ue(j,ik)
         enddo
      enddo
c-----> construct a and b coefficients for the k-point
       symvacvac=1
       if (nvac==1) symvacvac=2
       do symvac=1,symvacvac
         do i = 1,nv2d
            do n = 1,nslibd
               ac(i,n) = cmplx(0.0,0.0)
               bc(i,n) = cmplx(0.0,0.0)
            enddo   
            do n = 1,nslibd_b
               ac_b(i,n) = cmplx(0.0,0.0)
               bc_b(i,n) = cmplx(0.0,0.0)
            enddo   
         enddo   

        if (symvac==2) sign=-1.0   

      do k = 1,nv(jspin)
         l = map2(k)
         zks = k3(k,jspin)*bmat(3,3)*sign
         arg = zks*z1
         c_1 = cmplx(cos(arg),sin(arg)) * const
         av = -c_1 * cmplx( dte(l),zks*te(l) )
         bv =  c_1 * cmplx(  dt(l),zks* t(l) )
c-----> loop over basis functions
229
230
         IF(zMat%l_real) THEN
            do n = 1,nslibd
231
232
               ac(l,n) = ac(l,n) + zMat%data_r(k+addnoco,n)*av
               bc(l,n) = bc(l,n) + zMat%data_r(k+addnoco,n)*bv
233
234
235
            enddo
         ELSE
            do n = 1,nslibd
236
237
               ac(l,n) = ac(l,n) + zMat%data_c(k+addnoco,n)*av
               bc(l,n) = bc(l,n) + zMat%data_c(k+addnoco,n)*bv
238
239
            enddo
         END IF
240
241
242
243
      enddo

c...now for the k+b point

244
      evacp = evac_b(ivac)
245
246
247
248
      do ik = 1,nv2_b
         v(1) = bkpt_b(1) + kvac1_b(ik) + qss1
         v(2) = bkpt_b(2) + kvac2_b(ik) + qss2
         v(3) = 0.
Daniel Wortmann's avatar
Daniel Wortmann committed
249
         ev = evacp - 0.5*dot_product(v,matmul(bbmat,v))
250
251
252
         call vacuz(ev,vz_b(1,ivac),vz0_b(ivac),nmz,delz,t_b(ik),
     +        dt_b(ik),u_b(1,ik))
         call vacudz(ev,vz_b(1,ivac),vz0_b(ivac),nmz,delz,te_b(ik),
253
254
255
256
257
258
259
260
261
262
263
264
     +        dte_b(ik),tei_b(ik),ue_b(1,ik),dt_b(ik),
     +        u_b(1,ik))
         scale = wronk/ (te_b(ik)*dt_b(ik)-dte_b(ik)*t_b(ik))
         te_b(ik) = scale*te_b(ik)
         dte_b(ik) = scale*dte_b(ik)
         tei_b(ik) = scale*tei_b(ik)
         do j = 1,nmz
            ue_b(j,ik) = scale*ue_b(j,ik)
         enddo
      enddo
c-----> construct a and b coefficients for the k+b point

265
      do k = 1,nv_b(jspin_b)
266
         l = map2_b(k)
267
         zks = k3_b(k,jspin_b)*bmat(3,3)*sign
268
269
270
271
272
         arg = zks*z1
         c_1 = cmplx(cos(arg),sin(arg)) * const
         av = -c_1 * cmplx( dte_b(l),zks*te_b(l) )
         bv =  c_1 * cmplx( dt_b(l),zks*t_b(l) )
c-----> loop over basis functions
273
274
         IF(zMat_b%l_real) THEN
            do n = 1,nslibd_b
275
276
               ac_b(l,n) = ac_b(l,n) + zMat_b%data_r(k+addnoco,n)*av
               bc_b(l,n) = bc_b(l,n) + zMat_b%data_r(k+addnoco,n)*bv
277
278
279
            enddo
         ELSE
            do n = 1,nslibd_b
280
281
               ac_b(l,n) = ac_b(l,n) + zMat_b%data_c(k+addnoco,n)*av
               bc_b(l,n) = bc_b(l,n) + zMat_b%data_c(k+addnoco,n)*bv
282
283
            enddo
         END IF
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
      enddo


      do l = 1,nv2
      do lprime = 1,nv2_b
      if (kvac1(l).eq.(kvac1_b(lprime)-gb(1))
     & .and. kvac2(l).eq.(kvac2_b(lprime)-gb(2)))then
         zks = gb(3)*bmat(3,3)*sign

         do i = 1,nmz
             xx(np1-i) = u(i,l)*u_b(i,lprime)*
     *          cos((z1+(i-1)*delz)*zks)
            xximag(np1-i) = u(i,l)*u_b(i,lprime)*
     *          sin((z1+(i-1)*delz)*zks)
         enddo   
         call intgz0(xx,delz,nmz,uu,tail)
         call intgz0(xximag,delz,nmz,uuimag,tail)

         do i = 1,nmz
            xx(np1-i) = u(i,l)*ue_b(i,lprime)*
     *   cos((z1+(i-1)*delz)*zks)
            xximag(np1-i) = u(i,l)*ue_b(i,lprime)*
     *          sin((z1+(i-1)*delz)*zks)
         enddo   
         call intgz0(xx,delz,nmz,ud,tail)
         call intgz0(xximag,delz,nmz,udimag,tail)

         do i = 1,nmz
            xx(np1-i) = ue(i,l)*u_b(i,lprime)*
     *   cos((z1+(i-1)*delz)*zks)
            xximag(np1-i) = ue(i,l)*u_b(i,lprime)*
     *          sin((z1+(i-1)*delz)*zks)
         enddo   
         call intgz0(xx,delz,nmz,du,tail)
         call intgz0(xximag,delz,nmz,duimag,tail)
         do i = 1,nmz
            xx(np1-i) = ue(i,l)*ue_b(i,lprime)*
     *   cos((z1+(i-1)*delz)*zks)
            xximag(np1-i) = ue(i,l)*ue_b(i,lprime)*
     *          sin((z1+(i-1)*delz)*zks)
         enddo   
         call intgz0(xx,delz,nmz,dd,tail)
         call intgz0(xximag,delz,nmz,ddimag,tail)

         do i = 1,nslibd
            do j = 1,nslibd_b
               mmn(i,j) = mmn(i,j) + area*(
     *  ac(l,i)*conjg(ac_b(lprime,j))*cmplx(uu,uuimag) +
     +  ac(l,i)*conjg(bc_b(lprime,j))*cmplx(ud,udimag) +
     *  bc(l,i)*conjg(ac_b(lprime,j))*cmplx(du,duimag) +
334
     +  bc(l,i)*conjg(bc_b(lprime,j))*cmplx(dd,ddimag) )*vacchi
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
            enddo
         enddo





      endif ! kvac1=kvac1_b and kvac2=kvac2_b
      enddo ! lprime
      enddo ! l

      enddo !symvac
c... cycle by the vacua finishes
 140  enddo      

      deallocate ( ac,bc,dt,dte,t,te,tei,u,ue,
     +             v,kvac1,kvac2,map2 )
      deallocate ( ac_b,bc_b,dt_b,dte_b,t_b,te_b,tei_b,u_b,ue_b,
     +             kvac1_b,kvac2_b,map2_b )

      END SUBROUTINE wann_mmkb_vac
      END MODULE m_wann_mmkb_vac