Line data Source code
1 : !--------------------------------------------------------------------------------
2 : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
3 : ! This file is part of FLEUR and available as free software under the conditions
4 : ! of the MIT license as expressed in the LICENSE file in more detail.
5 : !--------------------------------------------------------------------------------
6 :
7 : module m_wann_perpmag_rs
8 : use m_juDFT
9 : contains
10 0 : subroutine wann_perpmag_rs(
11 0 : > rvecnum,rvec,kpoints,
12 : > jspins_in,nkpts,l_bzsym,film,
13 : > l_soc,band_min,band_max,neigd,
14 0 : > l_socmmn0,l_ndegen,ndegen,wan90version,
15 : > l_unformatted)
16 : c*************************************************
17 : c Calculate the matrix elements of
18 : c <n|B_{eff}|m> in real space from the
19 : c files WF1.chk (and WF1_um.dat) produced
20 : c by wannier90.
21 : c
22 : c Frank Freimuth
23 : c*************************************************
24 : use m_constants
25 : use m_wann_read_umatrix
26 :
27 : implicit none
28 : integer, intent(in) :: rvecnum
29 : integer, intent(in) :: rvec(:,:)
30 : real, intent(in) :: kpoints(:,:)
31 : integer, intent(in) :: jspins_in
32 : integer, intent(in) :: nkpts
33 : logical, intent(in) :: l_bzsym,l_soc
34 : logical, intent(in) :: film
35 : integer, intent(in) :: band_min(2),band_max(2),neigd
36 : logical, intent(in) :: l_socmmn0
37 : logical, intent(in) :: l_ndegen
38 : integer, intent(in) :: ndegen(:)
39 : integer, intent(in) :: wan90version
40 : logical, intent(in) :: l_unformatted
41 :
42 : integer :: ikpt,jspins
43 : integer :: kpts
44 : logical :: l_file
45 : integer :: num_wann,num_kpts,num_nnmax,jspin
46 : integer :: kspin,kkspin
47 : integer :: wann_shift,num_wann2
48 : integer :: i,j,k,m,info,r1,r2,r3,dummy1
49 : integer :: dummy2,dummy3
50 : integer :: hopmin,hopmax,counter,m1,m2
51 : integer :: num_bands2
52 : integer,allocatable :: iwork(:)
53 : real,allocatable :: energy(:,:),ei(:)
54 : real,allocatable :: eigw(:,:),rwork(:)
55 : complex,allocatable :: work(:),vec(:,:)
56 0 : complex,allocatable :: u_matrix(:,:,:,:),hwann(:,:,:,:)
57 0 : complex,allocatable :: hwannmix(:,:,:,:)
58 0 : complex,allocatable :: hreal(:,:,:,:)
59 0 : complex,allocatable :: paulimat(:,:,:,:)
60 0 : complex,allocatable :: paulimatmix(:,:,:,:)
61 0 : complex,allocatable :: paulimat_opt(:,:,:,:)
62 0 : complex,allocatable :: paulimat2(:,:,:,:)
63 : complex :: fac,eulav,eulav1
64 : real :: tmp_omi,rdotk,tpi,minenerg,maxenerg
65 : real, allocatable :: minieni(:),maxieni(:)
66 : character :: jobz,uplo
67 : integer :: kpt,band,lee,lwork,lrwork,liwork,n,lda
68 : complex :: value(4)
69 : logical :: um_format
70 : logical :: repro_eig
71 : logical :: l_chk,l_proj
72 : logical :: have_disentangled
73 0 : integer,allocatable :: ndimwin(:,:)
74 0 : logical,allocatable :: lwindow(:,:,:)
75 : integer :: chk_unit,nkp,ntmp,ierr
76 : character(len=33) :: header
77 : character(len=20) :: checkpoint
78 : real :: tmp_latt(3,3), tmp_kpt_latt(3,nkpts)
79 : real :: omega_invariant
80 0 : complex,allocatable :: u_matrix_opt(:,:,:,:)
81 : integer :: num_bands,counter1,counter2
82 : logical :: l_umdat
83 : real,allocatable :: eigval2(:,:)
84 : real,allocatable :: eigval_opt(:,:)
85 : real :: scale,a,b
86 : character(len=2) :: spinspin12(0:2)
87 : character(len=3) :: spin12(2)
88 : character(len=6) :: filename
89 : integer :: jp,mp,kk
90 : integer :: hopmin_x,hopmin_y,hopmin_z
91 : integer :: hopmax_x,hopmax_y,hopmax_z
92 : integer :: rvecind
93 : integer :: ii,int_dummy
94 : integer :: dir,num(3)
95 : integer :: spin1,spin2
96 : logical :: l_addint
97 : integer :: nbnd,nwfs,fullnkpts
98 0 : complex,allocatable :: amn(:,:,:)
99 0 : complex,allocatable :: perpmag(:,:,:,:)
100 :
101 : data spinspin12/' ','.1' , '.2'/
102 : data spin12/'WF1','WF2'/
103 :
104 0 : call timestart("wann_perpmag_rs")
105 0 : tpi=2*pimach()
106 :
107 0 : jspins=jspins_in
108 0 : if(l_soc)jspins=1
109 :
110 0 : write(oUnit,*)"nkpts=",nkpts
111 :
112 :
113 : c*****************************************************
114 : c get num_bands and num_wann from the proj file
115 : c*****************************************************
116 0 : do j=jspins,0,-1
117 0 : inquire(file=trim('proj'//spinspin12(j)),exist=l_file)
118 0 : if(l_file)then
119 0 : filename='proj'//spinspin12(j)
120 0 : exit
121 : endif
122 : enddo
123 0 : if(l_file)then
124 0 : open (203,file=trim(filename),status='old')
125 0 : rewind (203)
126 : else
127 : CALL juDFT_error("no proj/proj.1/proj.2",calledby
128 0 : + ="wann_perpmag_rs")
129 : endif
130 0 : read (203,*) num_wann,num_bands
131 0 : close (203)
132 0 : write(oUnit,*)'According to proj there are ',num_bands,' bands'
133 0 : write(oUnit,*)"and ",num_wann," wannier functions."
134 :
135 0 : allocate( u_matrix_opt(num_bands,num_wann,nkpts,jspins) )
136 0 : allocate( u_matrix(num_wann,num_wann,nkpts,jspins) )
137 0 : allocate( lwindow(num_bands,nkpts,jspins) )
138 0 : allocate( ndimwin(nkpts,jspins) )
139 :
140 0 : do jspin=1,jspins !spin loop
141 : c****************************************************************
142 : c read in chk
143 : c****************************************************************
144 0 : num_kpts=nkpts
145 :
146 : call wann_read_umatrix2(
147 : > nkpts,num_wann,num_bands,
148 : > um_format,jspin,wan90version,
149 : < have_disentangled,
150 : < lwindow(:,:,jspin),
151 : < ndimwin(:,jspin),
152 : < u_matrix_opt(:,:,:,jspin),
153 0 : < u_matrix(:,:,:,jspin) )
154 :
155 :
156 :
157 : enddo !jspin
158 :
159 0 : jspin=1
160 :
161 0 : num_bands2=jspins*num_bands
162 0 : num_wann2=jspins*num_wann
163 :
164 0 : allocate( paulimat(3,num_bands2,num_bands2,nkpts) )
165 0 : paulimat=0.0
166 :
167 : c****************************************************
168 : c Read the file "updown.perpmag".
169 : c****************************************************
170 0 : if(l_unformatted)then
171 0 : inquire(file='updown.perpmagint_unf',exist=l_addint)
172 0 : open(304,file='updown.perpmag_unf',form='unformatted')
173 0 : if(l_addint)then
174 0 : open(307,file='updown.perpmagint_unf',form='unformatted')
175 : endif
176 0 : read(304)nbnd,fullnkpts,nwfs
177 0 : allocate(amn(nbnd,nwfs,fullnkpts))
178 0 : read(304)amn(1:nbnd,1:nwfs,1:fullnkpts)
179 0 : do nkp=1,num_kpts
180 0 : do i=1,num_bands
181 0 : do j=1,num_bands
182 0 : paulimat(1,j,i+num_bands*(jspins-1),nkp)=amn(j,i,nkp)
183 : enddo !j
184 : enddo !i
185 : enddo !nkp
186 0 : if(l_addint)then
187 0 : read(307)nbnd,fullnkpts,nwfs
188 0 : read(307)amn(1:nbnd,1:nwfs,1:fullnkpts)
189 0 : do nkp=1,num_kpts
190 0 : do i=1,num_bands
191 0 : do j=1,num_bands
192 : paulimat(1,j,i+num_bands*(jspins-1),nkp)=
193 0 : & paulimat(1,j,i+num_bands*(jspins-1),nkp)+amn(j,i,nkp)
194 : enddo !j
195 : enddo !i
196 : enddo !nkp
197 : endif
198 0 : deallocate(amn)
199 : else
200 0 : inquire(file='updown.perpmagint',exist=l_addint)
201 0 : open(304,file='updown.perpmag',form='formatted')
202 0 : read(304,*)
203 0 : read(304,*)
204 0 : if(l_addint)then
205 0 : open(307,file='updown.perpmagint',form='formatted')
206 0 : read(307,*)
207 0 : read(307,*)
208 : endif
209 0 : do nkp=1,num_kpts
210 0 : do i=1,num_bands
211 0 : do j=1,num_bands
212 0 : read(304,*)dummy1,dummy2,dummy3,a,b
213 0 : paulimat(1,j,i+num_bands*(jspins-1),nkp)=cmplx(a,b)
214 : enddo !j
215 : enddo !i
216 0 : if(l_addint)then
217 0 : do i=1,num_bands
218 0 : do j=1,num_bands
219 0 : read(307,*)dummy1,dummy2,dummy3,a,b
220 : paulimat(1,j,i+num_bands*(jspins-1),nkp)=
221 : = paulimat(1,j,i+num_bands*(jspins-1),nkp)+
222 0 : + cmplx(a,b)
223 : enddo !j
224 : enddo !i
225 : endif
226 : enddo !nkp
227 : endif
228 0 : close(304)
229 0 : if(l_addint)then
230 0 : close(307)
231 : endif
232 :
233 0 : do nkp=1,num_kpts
234 0 : paulimat(2,:,:,nkp)=ImagUnit*paulimat(1,:,:,nkp)
235 : paulimat(1,:,:,nkp)=paulimat(1,:,:,nkp)+
236 0 : & transpose(conjg( paulimat(1,:,:,nkp) ))
237 : paulimat(2,:,:,nkp)=paulimat(2,:,:,nkp)+
238 0 : & transpose(conjg( paulimat(2,:,:,nkp) ))
239 : enddo
240 :
241 : c****************************************************************
242 : c Calculate matrix elements of B_{eff} in the basis of
243 : c rotated Bloch functions.
244 : c****************************************************************
245 0 : allocate( paulimat2(3,num_wann2,num_wann2,nkpts) )
246 : write(oUnit,*)"calculate matrix elements of exchange field
247 0 : & between wannier orbitals"
248 :
249 0 : if(have_disentangled) then
250 :
251 0 : allocate( paulimat_opt(3,num_bands2,num_bands2,nkpts) )
252 0 : allocate( paulimatmix(3,num_wann2,num_bands2,nkpts))
253 :
254 0 : do nkp=1,num_kpts
255 : counter1=0
256 0 : do m=1,num_bands2
257 0 : spin1=(m-1)/num_bands
258 0 : if(lwindow(m-spin1*num_bands,nkp,spin1+1))then
259 0 : counter1=counter1+1
260 0 : counter2=0
261 0 : do mp=1,num_bands2
262 0 : spin2=(mp-1)/num_bands
263 0 : if(lwindow(mp-spin2*num_bands,nkp,spin2+1))then
264 0 : counter2=counter2+1
265 0 : do k=1,3
266 : paulimat_opt(k,counter2,counter1,nkp)=
267 0 : & paulimat(k,mp,m,nkp)
268 : enddo
269 : endif
270 : enddo !mp
271 : endif
272 : enddo !m
273 : enddo
274 :
275 0 : paulimatmix=0.0
276 0 : do spin1=0,jspins-1
277 0 : do spin2=0,jspins-1
278 0 : do nkp=1,num_kpts
279 0 : do jp=1,num_wann
280 0 : do m=1,ndimwin(nkp,spin2+1)
281 0 : do mp=1,ndimwin(nkp,spin1+1)
282 0 : do k=1,3
283 :
284 :
285 : paulimatmix(k,jp+spin1*num_wann,m+spin2*ndimwin(nkp,1),nkp)=
286 : = paulimatmix(k,jp+spin1*num_wann,m+spin2*ndimwin(nkp,1),nkp)+
287 : & conjg(u_matrix_opt(mp,jp,nkp,spin1+1))*
288 : & paulimat_opt(k,mp+spin1*ndimwin(nkp,1),
289 0 : & m+spin2*ndimwin(nkp,1),nkp)
290 : enddo !k
291 : enddo !mp
292 : enddo !m
293 : enddo !jp
294 : enddo !nkp
295 : enddo !spin2
296 : enddo !spin1
297 :
298 0 : paulimat2=0.0
299 0 : do spin1=0,jspins-1
300 0 : do spin2=0,jspins-1
301 0 : do nkp=1,num_kpts
302 0 : do j=1,num_wann
303 0 : do jp=1,num_wann
304 0 : do m=1,ndimwin(nkp,spin2+1)
305 :
306 0 : do k=1,3
307 :
308 :
309 : paulimat2(k,jp+spin1*num_wann,j+spin2*num_wann,nkp)=
310 : = paulimat2(k,jp+spin1*num_wann,j+spin2*num_wann,nkp)+
311 : & paulimatmix(k,jp+spin1*num_wann,
312 :
313 : & m+spin2*ndimwin(nkp,1),nkp)*
314 0 : & u_matrix_opt(m,j,nkp,spin2+1)
315 : enddo !k
316 : enddo !m
317 : enddo !jp
318 : enddo !j
319 : enddo !nkp
320 : enddo !spin2
321 : enddo !spin1
322 :
323 :
324 0 : deallocate(paulimat_opt)
325 0 : deallocate(paulimatmix)
326 :
327 : else
328 0 : paulimat2=paulimat
329 : end if !have_disentangled
330 :
331 0 : allocate(hwann(3,num_wann2,num_wann2,num_kpts))
332 0 : hwann=cmplx(0.0,0.0)
333 0 : wann_shift=0
334 :
335 :
336 0 : allocate(hwannmix(3,num_wann2,num_wann2,num_kpts))
337 0 : hwannmix=cmplx(0.0,0.0)
338 :
339 0 : do spin1=0,jspins-1
340 0 : do spin2=0,jspins-1
341 0 : do k=1,num_kpts
342 0 : do i=1,num_wann
343 0 : do mp=1,num_wann
344 0 : do j=1,num_wann
345 :
346 0 : do kk=1,3
347 : hwannmix(kk,mp+spin1*num_wann,i+spin2*num_wann,k)=
348 : = hwannmix(kk,mp+spin1*num_wann,i+spin2*num_wann,k)+
349 : * conjg(u_matrix(j,mp,k,spin1+1))*
350 : * paulimat2(kk,j+spin1*num_wann,
351 0 : * i+spin2*num_wann,k)
352 :
353 : enddo !kk
354 : enddo !j
355 : enddo !mp
356 :
357 : enddo !i
358 : enddo !k
359 : enddo !spin2
360 : enddo !spin1
361 :
362 0 : do spin1=0,jspins-1
363 0 : do spin2=0,jspins-1
364 0 : do k=1,num_kpts
365 0 : do m=1,num_wann
366 0 : do i=1,num_wann
367 0 : do mp=1,num_wann
368 0 : do kk=1,3
369 : hwann(kk,mp+spin1*num_wann,m+spin2*num_wann,k)=
370 : = hwann(kk,mp+spin1*num_wann,m+spin2*num_wann,k)+
371 : * hwannmix(kk,mp+spin1*num_wann,
372 : * i+spin2*num_wann,k)*
373 0 : * u_matrix(i,m,k,spin2+1)
374 : enddo !kk
375 : enddo !mp
376 : enddo !i
377 : enddo !m
378 : enddo !k
379 : enddo !spin2
380 : enddo !spin1
381 :
382 0 : deallocate(hwannmix)
383 :
384 :
385 :
386 : c************************************************************
387 : c Calculate matrix elements in real space.
388 : c***********************************************************
389 0 : write(oUnit,*)"calculate pauli-mat in rs"
390 :
391 :
392 0 : allocate( hreal(3,num_wann2,num_wann2,rvecnum) )
393 0 : hreal=cmplx(0.0,0.0)
394 0 : if(l_ndegen)then
395 0 : do rvecind=1,rvecnum
396 0 : do k=1,nkpts
397 : rdotk=tpi*( kpoints(1,k)*rvec(1,rvecind)+
398 : & kpoints(2,k)*rvec(2,rvecind)+
399 0 : & kpoints(3,k)*rvec(3,rvecind) )
400 0 : fac=cmplx(cos(rdotk),-sin(rdotk))/real(ndegen(rvecind))
401 0 : do m2=1,num_wann2
402 0 : do m1=1,num_wann2
403 0 : do dir=1,3
404 : hreal(dir,m1,m2,rvecind)=hreal(dir,m1,m2,rvecind)+
405 0 : & fac*hwann(dir,m1,m2,k)
406 : enddo !dir
407 : enddo !m1
408 : enddo !m2
409 : enddo !k
410 : enddo !rvecind
411 : else
412 0 : do rvecind=1,rvecnum
413 0 : do k=1,nkpts
414 : rdotk=tpi*( kpoints(1,k)*rvec(1,rvecind)+
415 : & kpoints(2,k)*rvec(2,rvecind)+
416 0 : & kpoints(3,k)*rvec(3,rvecind) )
417 0 : fac=cmplx(cos(rdotk),-sin(rdotk))
418 0 : do m2=1,num_wann2
419 0 : do m1=1,num_wann2
420 0 : do dir=1,3
421 : hreal(dir,m1,m2,rvecind)=hreal(dir,m1,m2,rvecind)+
422 0 : & fac*hwann(dir,m1,m2,k)
423 : enddo !dir
424 : enddo !m1
425 : enddo !m2
426 : enddo !k
427 : enddo !rvecind
428 : endif
429 0 : hreal=hreal/cmplx(real(nkpts),0.0)
430 :
431 : c open(321,file='rsperpmag'//spinspin12(1),form='formatted')
432 : c do rvecind=1,rvecnum
433 : c r3=rvec(3,rvecind)
434 : c r2=rvec(2,rvecind)
435 : c r1=rvec(1,rvecind)
436 : c do j=1,num_wann2
437 : c do i=1,num_wann2
438 : c do kk=1,3
439 : c write(321,'(i3,1x,i3,1x,i3,1x,i3,
440 : c & 1x,i3,1x,i3,1x,f20.8,1x,f20.8)')
441 : c & r1,r2,r3,i,j,kk,hreal(kk,i,j,rvecind)
442 : c enddo !kk
443 : c enddo !i
444 : c enddo !j
445 : c enddo !rvecnum
446 : c close(321)
447 0 : if(l_unformatted)then
448 0 : allocate(perpmag(num_wann2,num_wann2,3,rvecnum))
449 0 : if(l_ndegen)then
450 : open(321,file='tunfndegen',form='unformatted'
451 : #ifdef CPP_INTEL
452 : & ,convert='BIG_ENDIAN'
453 : #endif
454 0 : & )
455 : else
456 : open(321,file='tunf',form='unformatted'
457 : #ifdef CPP_INTEL
458 : & ,convert='BIG_ENDIAN'
459 : #endif
460 0 : & )
461 : endif
462 0 : do rvecind=1,rvecnum
463 0 : do j=1,num_wann2
464 0 : do i=1,num_wann2
465 0 : do kk=1,3
466 0 : perpmag(i,j,kk,rvecind)=hreal(kk,i,j,rvecind)
467 : enddo !kk
468 : enddo !i
469 : enddo !j
470 : enddo !rvecind
471 0 : write(321)perpmag
472 0 : deallocate(perpmag)
473 : else
474 0 : if(l_ndegen)then
475 :
476 : open(321,file='rsperpmag_ndegen'//spinspin12(1),
477 : & form='formatted',
478 0 : & recl=1000)
479 : else
480 : open(321,file='rsperpmag'//spinspin12(1),form='formatted',
481 0 : & recl=1000)
482 : endif
483 :
484 0 : do rvecind=1,rvecnum
485 0 : r3=rvec(3,rvecind)
486 0 : r2=rvec(2,rvecind)
487 0 : r1=rvec(1,rvecind)
488 0 : do j=1,num_wann2
489 0 : do i=1,num_wann2
490 0 : do kk=1,3
491 0 : write(321,*) r1,r2,r3,i,j,kk,real(hreal(kk,i,j,rvecind)),
492 0 : & aimag(hreal(kk,i,j,rvecind))
493 : enddo !kk
494 : enddo !i
495 : enddo !j
496 : enddo !rvecnum
497 : endif !l_unformatted
498 0 : close(321)
499 :
500 :
501 0 : deallocate(lwindow,u_matrix_opt,ndimwin)
502 0 : deallocate(u_matrix,hwann,hreal)
503 :
504 0 : call timestop("wann_perpmag_rs")
505 0 : end subroutine wann_perpmag_rs
506 : end module m_wann_perpmag_rs
|
