Line data Source code
1 : module m_wann_torque_rs
2 : use m_juDFT
3 : contains
4 0 : subroutine wann_torque_rs(
5 0 : > ntype,neq,rvecnum,rvec,kpoints,
6 : > jspins_in,nkpts,l_bzsym,film,
7 : > l_soc,band_min,band_max,neigd,
8 0 : > l_ndegen,ndegen,wan90version,
9 : > l_unformatted)
10 : c*************************************************
11 : c Calculate the matrix elements of the
12 : c muffin-tin spin current operator
13 : c in real space from the
14 : c files WF1.chk (and WF1_um.dat) produced
15 : c by wannier90.
16 : c
17 : c Frank Freimuth, January 2011
18 : c*************************************************
19 : use m_constants, only:pimach
20 : use m_wann_read_umatrix
21 :
22 : implicit none
23 : integer, intent(in) :: ntype
24 : integer, intent(in) :: neq(:)
25 : integer, intent(in) :: rvecnum
26 : integer, intent(in) :: rvec(:,:)
27 : real, intent(in) :: kpoints(:,:)
28 :
29 : integer, intent(in) :: jspins_in
30 : integer, intent(in) :: nkpts
31 : logical, intent(in) :: l_bzsym,l_soc
32 : logical, intent(in) :: film
33 :
34 :
35 : integer, intent(in) :: band_min(2),band_max(2),neigd
36 :
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 :: nbnd,fullnkpts
43 : integer :: ikpt,jspins
44 : integer :: kpts
45 : logical :: l_file
46 : integer :: num_wann,num_kpts,num_nnmax,jspin
47 : integer :: kspin,kkspin
48 : integer :: num_wann2
49 : integer :: i,j,k,m,info,r1,r2,r3,dummy1
50 : integer :: dummy2,dummy3
51 : integer :: counter,m1,m2
52 : integer :: num_bands2
53 : integer,allocatable :: iwork(:)
54 : real,allocatable :: energy(:,:),ei(:)
55 : real,allocatable :: eigw(:,:),rwork(:)
56 : complex,allocatable :: work(:),vec(:,:)
57 0 : complex,allocatable :: u_matrix(:,:,:,:),hwann(:,:,:,:)
58 0 : complex,allocatable :: hwannmix(:,:,:,:)
59 0 : complex,allocatable :: hreal(:,:,:,:)
60 : complex,allocatable :: hrealsoc(:,:,:,:,:,:,:)
61 : complex,allocatable :: hwannsoc(:,:,:,:,:)
62 0 : complex,allocatable :: paulimat_opt(:,:,:,:)
63 0 : complex,allocatable :: nablamat(:,:,:,:)
64 0 : complex,allocatable :: nablamat2(:,:,:,:)
65 0 : complex,allocatable :: nablamatmix(:,:,:,:)
66 : complex :: fac,eulav,eulav1
67 : real :: tmp_omi,rdotk,tpi,minenerg,maxenerg
68 : real, allocatable :: minieni(:),maxieni(:)
69 : character :: jobz,uplo
70 : integer :: kpt,band,lee,lwork,lrwork,liwork,n,lda
71 : complex :: value(4)
72 : logical :: um_format
73 : logical :: repro_eig
74 : logical :: l_chk,l_proj
75 : logical :: have_disentangled
76 0 : integer,allocatable :: ndimwin(:,:)
77 0 : logical,allocatable :: lwindow(:,:,:)
78 : integer :: chk_unit,nkp,ntmp,ierr
79 : character(len=33) :: header
80 : character(len=20) :: checkpoint
81 : real :: tmp_latt(3,3), tmp_kpt_latt(3,nkpts)
82 : real :: omega_invariant
83 0 : complex,allocatable :: u_matrix_opt(:,:,:,:)
84 : integer :: num_bands,counter1,counter2
85 : logical :: l_umdat
86 : real,allocatable :: eigval2(:,:)
87 : real,allocatable :: eigval_opt(:,:)
88 : real :: scale,a,b
89 : character(len=2) :: spinspin12(0:2)
90 : character(len=3) :: spin12(2)
91 : character(len=6) :: filename
92 : integer :: jp,mp,kk
93 : integer :: rvecind
94 : integer :: spin1,spin2
95 : character(len=10) :: torquename
96 : character(len=9) :: torquename2
97 : character(len=15) :: torquename3
98 : integer :: at,nn,itype
99 0 : complex,allocatable :: perpmag(:,:,:,:)
100 : logical :: l_yesss
101 :
102 : data spinspin12/' ','.1' , '.2'/
103 : data spin12/'WF1','WF2'/
104 :
105 0 : call timestart("wann_torque_rs")
106 :
107 0 : tpi=2*pimach()
108 :
109 0 : jspins=jspins_in
110 0 : if(l_soc)jspins=1
111 :
112 0 : write(6,*)"nkpts=",nkpts
113 :
114 : c*****************************************************
115 : c get num_bands and num_wann from the proj file
116 : c*****************************************************
117 0 : do j=jspins,0,-1
118 0 : inquire(file=trim('proj'//spinspin12(j)),exist=l_file)
119 0 : if(l_file)then
120 0 : filename='proj'//spinspin12(j)
121 0 : exit
122 : endif
123 : enddo
124 0 : if(l_file)then
125 0 : open (203,file=trim(filename),status='old')
126 0 : rewind (203)
127 : else
128 0 : stop 'no proj/proj.1/proj.2'
129 : endif
130 0 : read (203,*) num_wann,num_bands
131 0 : close (203)
132 0 : write(6,*)'According to proj there are ',num_bands,' bands'
133 0 : write(6,*)"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 : enddo !jspin
156 :
157 0 : jspin=1
158 :
159 0 : num_bands2=jspins*num_bands
160 0 : num_wann2=jspins*num_wann
161 :
162 : c************************************************
163 : c Read the file torque_...
164 : c************************************************
165 0 : allocate( nablamat(3,num_bands2,num_bands2,nkpts) )
166 0 : at=0
167 0 : do itype=1,ntype
168 0 : do nn=1,neq(itype)
169 0 : at=at+1
170 0 : write(torquename,fmt='("torque_",i3.3)')at
171 0 : if(l_unformatted)then
172 :
173 : inquire(file=trim(torquename)//'_unf',
174 0 : & exist=l_yesss)
175 0 : if(.not.l_yesss)cycle
176 :
177 : open(304,file=trim(torquename)//'_unf',
178 0 : & form='unformatted')
179 0 : read(304)nbnd,nbnd,fullnkpts
180 0 : if(nbnd.ne.num_bands2)stop'problem1'
181 0 : if(fullnkpts.ne.nkpts)stop'problem2'
182 0 : read(304)nablamat(1:3,1:nbnd,1:nbnd,1:fullnkpts)
183 : else
184 0 : open(304,file=torquename,form='formatted')
185 0 : read(304,*)
186 0 : read(304,*)
187 0 : do nkp=1,num_kpts
188 0 : do i=1,num_bands2
189 0 : do j=1,num_bands2
190 0 : do k=1,3
191 0 : read(304,*)dummy1,dummy2,dummy3,a,b
192 0 : nablamat(k,j,i,nkp)=cmplx(a,b)
193 : enddo !k
194 : enddo !j
195 : enddo !i
196 : enddo !nkp
197 : endif !l_unformatted
198 0 : close(304)
199 :
200 : c********************************************************************
201 : c Calculate matrix elements of MT spin current in the basis of
202 : c rotated Bloch functions.
203 : c********************************************************************
204 0 : allocate( nablamat2(3,num_wann2,num_wann2,nkpts) )
205 : write(6,*)"calculate matrix elements of MT spin current
206 0 : & operator between wannier orbitals"
207 :
208 0 : if(have_disentangled) then
209 :
210 0 : allocate( paulimat_opt(3,num_bands2,num_bands2,nkpts) )
211 0 : allocate( nablamatmix(3,num_bands2,num_wann2,nkpts) )
212 :
213 0 : do nkp=1,num_kpts
214 : counter1=0
215 0 : do m=1,num_bands2
216 0 : spin1=(m-1)/num_bands
217 0 : if(lwindow(m-spin1*num_bands,nkp,spin1+1))then
218 0 : counter1=counter1+1
219 0 : counter2=0
220 0 : do mp=1,num_bands2
221 0 : spin2=(mp-1)/num_bands
222 0 : if(lwindow(mp-spin2*num_bands,nkp,spin2+1))then
223 0 : counter2=counter2+1
224 0 : do k=1,3
225 : paulimat_opt(k,counter2,counter1,nkp)=
226 0 : & nablamat(k,mp,m,nkp)
227 : enddo
228 : endif
229 : enddo !mp
230 : endif
231 : enddo !m
232 : enddo
233 :
234 0 : nablamatmix=0.0
235 0 : do spin1=0,jspins-1
236 0 : do spin2=0,jspins-1
237 0 : do nkp=1,num_kpts
238 0 : do j=1,num_wann
239 0 : do m=1,ndimwin(nkp,spin2+1)
240 0 : do mp=1,ndimwin(nkp,spin1+1)
241 0 : do k=1,3
242 : nablamatmix(k,mp+spin1*ndimwin(nkp,1),j+spin2*num_wann,nkp)=
243 : &nablamatmix(k,mp+spin1*ndimwin(nkp,1),j+spin2*num_wann,nkp)+
244 : & paulimat_opt(k,mp+spin1*ndimwin(nkp,1),
245 : & m +spin2*ndimwin(nkp,1),nkp)*
246 0 : & u_matrix_opt(m,j,nkp,spin2+1)
247 : enddo !k
248 : enddo !mp
249 : enddo !m
250 : enddo !j
251 : enddo !nkp
252 : enddo !spin2
253 : enddo !spin1
254 :
255 0 : deallocate(paulimat_opt)
256 :
257 0 : nablamat2=0.0
258 0 : do spin1=0,jspins-1
259 0 : do spin2=0,jspins-1
260 0 : do nkp=1,num_kpts
261 0 : do j=1,num_wann
262 0 : do jp=1,num_wann
263 0 : do mp=1,ndimwin(nkp,spin1+1)
264 0 : do k=1,3
265 : nablamat2(k,jp+spin1*num_wann,j+spin2*num_wann,nkp)=
266 : & nablamat2(k,jp+spin1*num_wann,j+spin2*num_wann,nkp)+
267 : & conjg(u_matrix_opt(mp,jp,nkp,spin1+1))*
268 0 : & nablamatmix(k,mp+spin1*ndimwin(nkp,1),j+spin2*num_wann,nkp)
269 : enddo !k
270 : enddo !mp
271 : enddo !jp
272 : enddo !j
273 : enddo !nkp
274 : enddo !spin2
275 : enddo !spin1
276 0 : deallocate( nablamatmix )
277 : else
278 0 : nablamat2(:,:,:,:)=nablamat(:,:,:,:)
279 : end if !have_disentangled
280 :
281 0 : allocate(hwann(3,num_wann2,num_wann2,num_kpts))
282 0 : hwann=cmplx(0.0,0.0)
283 :
284 0 : allocate(hwannmix(3,num_wann2,num_wann2,num_kpts))
285 0 : hwannmix=cmplx(0.0,0.0)
286 :
287 0 : do spin2=0,jspins-1
288 0 : do spin1=0,jspins-1
289 0 : do k=1,num_kpts
290 0 : do m=1,num_wann
291 0 : do j=1,num_wann
292 0 : do mp=1,num_wann
293 0 : do kk=1,3
294 : hwannmix(kk,mp+spin1*num_wann,m+spin2*num_wann,k)=
295 : & hwannmix(kk,mp+spin1*num_wann,m+spin2*num_wann,k)+
296 : * conjg(u_matrix(j,mp,k,spin1+1))*
297 : * nablamat2(kk,j+spin1*num_wann,
298 0 : * m+spin2*num_wann,k)
299 : enddo !kk
300 : enddo !mp
301 : enddo !j
302 : enddo !m
303 : enddo !k
304 : enddo !spin2
305 : enddo !spin1
306 :
307 :
308 0 : do spin2=0,jspins-1
309 0 : do spin1=0,jspins-1
310 0 : do k=1,num_kpts
311 0 : do m=1,num_wann
312 0 : do i=1,num_wann
313 0 : do mp=1,num_wann
314 0 : do kk=1,3
315 : hwann(kk,mp+spin1*num_wann,m+spin2*num_wann,k)=
316 : * hwann(kk,mp+spin1*num_wann,m+spin2*num_wann,k)+
317 : * hwannmix(kk,mp+spin1*num_wann,
318 : * i+spin2*num_wann,k)*
319 0 : * u_matrix(i,m,k,spin2+1)
320 : enddo !kk
321 : enddo !mp
322 : enddo !i
323 : enddo !m
324 : enddo !k
325 : enddo !spin1
326 : enddo !spin2
327 0 : deallocate(hwannmix)
328 0 : deallocate(nablamat2)
329 : c************************************************************
330 : c Calculate matrix elements in real space.
331 : c***********************************************************
332 0 : write(6,*)"calculate nabla-mat in rs"
333 :
334 0 : allocate(hreal(3,num_wann2,num_wann2,rvecnum))
335 0 : hreal=cmplx(0.0,0.0)
336 :
337 0 : if(l_ndegen)then
338 0 : do rvecind=1,rvecnum
339 0 : do k=1,nkpts
340 : rdotk=tpi*(kpoints(1,k)*rvec(1,rvecind)+
341 : & kpoints(2,k)*rvec(2,rvecind)+
342 0 : & kpoints(3,k)*rvec(3,rvecind) )
343 : fac=cmplx(cos(rdotk),-sin(rdotk))/
344 0 : / cmplx(real(ndegen(rvecind)),0.0)
345 0 : do i=1,num_wann2
346 0 : do j=1,num_wann2
347 0 : do kk=1,3
348 : hreal(kk,j,i,rvecind)=
349 : & hreal(kk,j,i,rvecind)+
350 0 : & fac*hwann(kk,j,i,k)
351 : enddo !kk
352 : enddo !j
353 : enddo !i
354 : enddo !k
355 : enddo !rvecind
356 : else
357 0 : do rvecind=1,rvecnum
358 0 : do k=1,nkpts
359 : rdotk=tpi*(kpoints(1,k)*rvec(1,rvecind)+
360 : & kpoints(2,k)*rvec(2,rvecind)+
361 0 : & kpoints(3,k)*rvec(3,rvecind) )
362 0 : fac=cmplx(cos(rdotk),-sin(rdotk))
363 0 : do i=1,num_wann2
364 0 : do j=1,num_wann2
365 0 : do kk=1,3
366 : hreal(kk,j,i,rvecind)=
367 : & hreal(kk,j,i,rvecind)+
368 0 : & fac*hwann(kk,j,i,k)
369 : enddo !kk
370 : enddo !j
371 : enddo !i
372 : enddo !k
373 : enddo !rvecind
374 : endif
375 0 : hreal=hreal/cmplx(real(nkpts),0.0)
376 0 : deallocate(hwann)
377 0 : if(l_unformatted)then
378 0 : if(l_ndegen)then
379 0 : write(torquename3,fmt='("rstorndegen_",i3.3)')at
380 : open(321,file=trim(torquename3)//'_unf',
381 : & form='unformatted'
382 : #ifdef CPP_INTEL
383 : & ,convert='BIG_ENDIAN'
384 : #endif
385 0 : & )
386 : else ! l_ndegen
387 0 : write(torquename2,fmt='("rstor_",i3.3)')at
388 : open(321,file=trim(torquename2)//'_unf',
389 : & form='unformatted'
390 : #ifdef CPP_INTEL
391 : & ,convert='BIG_ENDIAN'
392 : #endif
393 0 : & )
394 : endif ! l_ndegen
395 0 : allocate(perpmag(num_wann2,num_wann2,3,rvecnum))
396 0 : do rvecind=1,rvecnum
397 0 : do j=1,num_wann2
398 0 : do i=1,num_wann2
399 0 : do kk=1,3
400 0 : perpmag(i,j,kk,rvecind)=hreal(kk,i,j,rvecind)
401 : enddo !kk
402 : enddo !i
403 : enddo !j
404 : enddo !rvecind
405 0 : write(321)perpmag
406 0 : deallocate(perpmag)
407 : else
408 0 : if(l_ndegen)then
409 0 : write(torquename3,fmt='("rstorndegen_",i3.3)')at
410 0 : open(321,file=torquename3)
411 : else
412 0 : write(torquename2,fmt='("rstor_",i3.3)')at
413 0 : open(321,file=torquename2)
414 : endif
415 0 : do rvecind=1,rvecnum
416 0 : r3=rvec(3,rvecind)
417 0 : r2=rvec(2,rvecind)
418 0 : r1=rvec(1,rvecind)
419 0 : do j=1,num_wann2
420 0 : do i=1,num_wann2
421 0 : do kk=1,3
422 : write(321,'(i3,1x,i3,1x,i3,1x,i3,
423 : & 1x,i3,1x,i3,1x,f20.8,1x,f20.8)')
424 0 : & r1,r2,r3,i,j,kk,hreal(kk,i,j,rvecind)
425 : enddo !kk
426 : enddo !j
427 : enddo !i
428 : enddo !rvecind
429 : endif !l_unformatted
430 0 : close(321)
431 0 : deallocate(hreal)
432 : enddo !nn
433 : enddo !itype
434 :
435 0 : call timestop("wann_torque_rs")
436 0 : end subroutine wann_torque_rs
437 : end module m_wann_torque_rs
|
