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_updown
8 : use m_juDFT
9 : CONTAINS
10 0 : SUBROUTINE wann_updown(
11 : > mpi,input,sym,atoms,stars,vacuum,sphhar,oneD,noco,cell,vTot,
12 : > enpara,
13 0 : > eig_id,l_real,mpi_comm,l_dulo,l_noco,l_ss,lmaxd,ntypd,
14 : > neigd,natd,nop,nvd,jspd,nbasfcn,llod,nlod,ntype,
15 0 : > omtil,nlo,llo,lapw_l,invtab,mrot,ngopr,neq,lmax,
16 0 : > invsat,invsatnr,nkpt,taual,rmt,amat,bmat,bbmat,alph,
17 0 : > beta,qss,sk2,phi2,odi,ods,irank,isize,n3d,nmzxyd,nmzd,
18 0 : > jmtd,nlhd,nq3,nvac,invs,invs2,film,nlh,jri,ntypsd,
19 0 : > ntypsy,jspins,nkptd,dx,n2d,rmsh,e1s,e2s,ulo_der,
20 0 : > ustep,ig,k1d,k2d,k3d,rgphs,slice,kk,nnne,
21 0 : > z1,nv2d,nmzxy,nmz,delz,ig2,area,tau,zatom,nq2,nop2,
22 : > volint,symor,pos,ef,l_soc,
23 : > memd,lnonsph,clnu,lmplmd,mlh,nmem,llh,lo1l,
24 : > theta,phi)
25 : c****************************************************************************
26 : c Calculate spin-up-spin-down matrix elements.
27 : c Useful in the case of spin-orbit coupling and noncollinear magnetism.
28 : c Implemented matrix elements:
29 : c (i) Overlaps => Pauli matrix elements.
30 : c (ii) Momentum operator => Spin current.
31 : c (iii) H_{SOC} => Include SOC in the basis set of Wannier functions.
32 : c (iv) Spin-current into the MT-spheres => Torque.
33 : c (v) Matrix elements of exchange field => Torque.
34 : c Frank Freimuth
35 : c****************************************************************************
36 : use m_types
37 : use m_juDFT
38 : use m_wann_rw_eig
39 : use m_abcof
40 : use m_radfun
41 : use m_radflo
42 : use m_cdnread, only : cdn_read0, cdn_read
43 : use m_constants, only : pimach
44 : ! use m_wann_mmk0_od_vac
45 : use m_wann_mmk0_vac
46 : use m_wann_mmk0_updown_sph
47 : use m_wann_abinv
48 : use m_wann_kptsrotate
49 : use m_wann_read_inp
50 : use m_wann_radovlp_integrals
51 : use m_matmul,only : matmul3,matmul3r
52 : use m_wann_socmat
53 : use m_wann_write_hsomtx
54 : use m_wann_write_nabla
55 : use m_wann_write_amn
56 : use m_wann_write_nabla
57 : use m_wann_mmkb_int
58 : #ifdef CPP_TOPO
59 : use m_wann_perpmag
60 : use m_wann_nabla_updown
61 : use m_wann_surfcurr_updown
62 : use m_wann_surfcurr_int_updown
63 : use m_wann_gabeffgagaunt
64 : #endif
65 : IMPLICIT NONE
66 : #include "cpp_double.h"
67 : #ifdef CPP_MPI
68 : include 'mpif.h'
69 : integer ierr(3)
70 : integer cpu_index
71 : integer stt(MPI_STATUS_SIZE)
72 :
73 : #endif
74 :
75 : TYPE(t_mpi),INTENT(IN) :: mpi
76 : TYPE(t_input),INTENT(IN) :: input
77 : TYPE(t_sym),INTENT(IN) :: sym
78 : TYPE(t_atoms),INTENT(IN) :: atoms
79 : TYPE(t_stars),INTENT(IN) :: stars
80 : TYPE(t_vacuum),INTENT(IN) :: vacuum
81 : TYPE(t_sphhar),INTENT(IN) :: sphhar
82 : TYPE(t_oneD),INTENT(IN) :: oneD
83 : TYPE(t_noco),INTENT(IN) :: noco
84 : TYPE(t_cell),INTENT(IN) :: cell
85 : TYPE(t_potden),INTENT(IN) :: vTot
86 : TYPE(t_enpara),INTENT(IN) :: enpara
87 :
88 : logical, intent (in) :: invs,invs2,film,slice,symor,l_real
89 : integer, intent (in) :: lmaxd,ntypd,neigd,nkptd,kk,nnne,mpi_comm
90 : integer, intent (in) :: natd,nop,nvd,jspd,nbasfcn,nq2,nop2,eig_id
91 : integer, intent (in) :: llod,nlod,ntype,n3d,n2d
92 : integer, intent (in) :: nmzxyd,nmzd,jmtd,nlhd,nq3,nvac
93 : integer, intent (in) :: ntypsd,jspins,k1d,k2d,k3d
94 : real, intent (in) :: omtil,e1s,e2s,delz,area,z1,volint
95 : integer, intent (in) :: ig(-k1d:k1d,-k2d:k2d,-k3d:k3d)
96 : complex, intent (in) :: rgphs(-k1d:k1d,-k2d:k2d,-k3d:k3d)
97 : integer, intent (in) :: nlh(ntypsd),jri(ntypd),ntypsy(natd)
98 : integer, intent (in) :: nlo(ntypd),llo(nlod,ntypd),lapw_l(ntypd)
99 : integer, intent (in) :: invtab(nop),mrot(3,3,nop),ngopr(natd)
100 : integer, intent (in) :: neq(ntypd),lmax(ntypd)
101 : integer, intent (in) :: invsat(natd),invsatnr(natd),nkpt
102 : integer, intent (in) :: irank,isize,nv2d,nmzxy,nmz
103 : integer, intent (in) :: ulo_der(nlod,ntypd),ig2(n3d)
104 : real, intent (in) :: taual(3,natd),rmt(ntypd),dx(ntypd)
105 : real, intent (in) :: amat(3,3),bmat(3,3),bbmat(3,3)
106 : real, intent (in) :: rmsh(jmtd,ntypd),tau(3,nop),zatom(ntype)
107 : real, intent (in) :: alph(ntypd),beta(ntypd),qss(3)
108 : real, intent (in) :: pos(3,natd)
109 : real, intent(in) :: ef
110 : complex, intent (in) :: ustep(n3d)
111 : logical, intent (in) :: l_dulo(nlod,ntypd),l_noco,l_ss
112 : logical, intent (in) :: l_soc
113 : integer, intent (in) :: memd
114 : integer, intent (in) :: lnonsph(ntypd)
115 : complex, intent (in) :: clnu(memd,0:nlhd,ntypsd)
116 : integer, intent (in) :: lmplmd
117 : integer, intent (in) :: mlh(memd,0:nlhd,ntypsd)
118 : integer, intent (in) :: nmem(0:nlhd,ntypsd)
119 : integer, intent (in) :: llh(0:nlhd,ntypsd)
120 : integer, intent (in) :: lo1l(0:llod,ntypd)
121 : real, intent (in) :: sk2(n2d),phi2(n2d)
122 : type (od_inp), intent (in) :: odi
123 : type (od_sym), intent (in) :: ods
124 : real, intent(in) :: theta,phi
125 :
126 : cccccccccccccccccc local variables cccccccccccccccccccc
127 0 : type(t_wann) :: wann
128 : integer :: lmd,nlotot,n,nmat,iter,ikpt,ikpt_b,nmat_b
129 : integer :: addnoco,addnoco2,funbas,loplod
130 : integer :: noccbd,noccbd_b,nn,nkpts,i,jspin,j,l,i_rec,m,nwf,nwfp
131 : integer :: jsp_start,jsp_end,nrec,nrec1,nbands,nbands_b
132 : integer :: nodeu,noded,n_size,na,n_rank,nbnd,numbands
133 : integer :: i1,i2,i3,in,lda
134 : integer :: n_bands(0:neigd),nslibd,nslibd_b
135 : character(len=8) :: dop,iop,name(10)
136 : real :: bkpt(3),bkpt_b(3),sfp,tpi,wronk,wk,wk_b,phase
137 : complex :: c_phase, zzConjgTemp
138 0 : real :: eig(neigd),cp_time(9),efermi
139 : logical :: l_p0,l_bkpts,l_proj,l_amn,l_mmn
140 : !!! energy window boundaries
141 0 : integer, allocatable :: nv(:)
142 : integer, allocatable :: nv_b(:)
143 0 : integer, allocatable :: k1(:,:),k2(:,:),k3(:,:)
144 : integer, allocatable :: k1_b(:,:),k2_b(:,:),k3_b(:,:)
145 :
146 0 : real, allocatable :: eigg(:,:)
147 0 : real kpoints(nkptd)
148 : !!! a and b coeff. constructed for each k-point
149 0 : complex, allocatable :: acof(:,:,:,:)
150 0 : complex, allocatable :: bcof(:,:,:,:)
151 0 : complex, allocatable :: ccof(:,:,:,:,:)
152 : !!! the parameters for the number of wfs
153 : integer :: nwfs
154 : !!! the potential in the spheres and the vacuum
155 0 : real, allocatable :: vr(:,:,:),vz(:,:,:)
156 : !!! bkpts data
157 : integer nntot,ikpt_help
158 : integer, allocatable :: gb(:,:,:),bpt(:,:)
159 : !!! radial wavefunctions in the muffin-tins and more ...
160 0 : real, allocatable :: flo(:,:,:,:,:)
161 0 : real, allocatable :: ff(:,:,:,:,:),gg(:,:,:,:,:)
162 :
163 0 : real :: uuilon(nlod,ntypd),duilon(nlod,ntypd)
164 0 : real :: ulouilopn(nlod,nlod,ntypd)
165 : !!! energy parameters
166 0 : real :: ello(nlod,ntypd,max(2,jspd)),evac(2,max(2,jspd))
167 0 : real :: epar(0:lmaxd,ntypd,max(2,jspd)),evdu(2,max(jspd,2))
168 : !!! the Mmn matrices
169 0 : complex, allocatable :: surfcurr(:,:,:,:)
170 0 : complex, allocatable :: hsomtx(:,:,:,:,:),hsomtx_l(:,:,:,:)
171 0 : complex, allocatable :: mmnk(:,:,:,:),mmn(:,:,:)
172 0 : complex, allocatable :: perpmag(:,:,:)
173 0 : complex, allocatable :: amn(:,:,:),nablamat(:,:,:,:)
174 : complex, allocatable :: socamn(:,:,:)
175 0 : complex, allocatable :: socmmn(:,:,:)
176 : complex, allocatable :: socmmnk(:,:,:,:)
177 : complex, allocatable :: omn(:,:,:),vec(:,:),a(:)
178 : complex, allocatable :: smn(:,:,:)
179 : complex, allocatable :: psiw(:,:,:)
180 : c..wf-hamiltonian in reciprocal space Hwf(k)
181 : complex, allocatable :: hwfk(:,:,:)
182 : c..egenvalues of the wf-hamiltonian Hwf(k)
183 : real, allocatable :: eigw(:,:)
184 : c..wf-hamiltonian in real space (hopping in the same unit cell)
185 : complex, allocatable :: hwfr(:,:)
186 : real, allocatable :: ei(:)
187 : complex, allocatable :: work(:)
188 : real,allocatable::centers(:,:,:)
189 : integer, allocatable :: iwork(:)
190 : integer :: info,lwork,lrwork,liwork,lee
191 : character :: jobz, uplo
192 : logical :: l_file
193 : logical :: l_amn2, l_conjugate
194 : character(len=3) :: spin12(2)
195 : data spin12/'WF1' , 'WF2'/
196 : character(len=30) :: task
197 0 : integer,allocatable::irreduc(:)
198 0 : integer,allocatable::mapkoper(:)
199 : integer :: fullnkpts,kpt,kptibz,kptibz_b,j1,j2,j3,oper,oper_b,k
200 : real :: bkrot(3)
201 : integer :: ios,kplot,kplotoper,plotoper,gfcut
202 : integer :: funbas_start,funbas_stop
203 : complex :: phasust
204 : integer,allocatable::pair_to_do(:,:)
205 0 : integer :: ngopr1(natd)
206 : integer,allocatable::maptopair(:,:,:)
207 : integer :: ldu,ldv,wannierspin,jspin2
208 : complex, allocatable :: amn_copy(:,:)
209 : complex, allocatable :: sunit(:,:)
210 : real, allocatable :: sdiag(:)
211 : complex, allocatable :: umn(:,:),vmn(:,:)
212 : real, allocatable :: rwork(:)
213 : character :: jobu*1,jobv*1
214 : real,allocatable::kdiff(:,:)
215 0 : integer,allocatable :: shiftkpt(:,:)
216 : integer :: unigrid(6),gfthick
217 : complex,allocatable::ujug(:,:,:,:),ujdg(:,:,:,:)
218 : complex,allocatable::djug(:,:,:,:),djdg(:,:,:,:)
219 : complex,allocatable::ujulog(:,:,:,:)
220 : complex,allocatable::djulog(:,:,:,:)
221 : complex,allocatable::ulojug(:,:,:,:)
222 : complex,allocatable::ulojdg(:,:,:,:)
223 : complex,allocatable::ulojulog(:,:,:,:)
224 : real delta,delta1,time_interstitial,time_mmn
225 : real time_total,delta2,delta3
226 : real time_lapw_expand,time_rw,time_symm,time_film
227 : real time_lapw_plot,time_ujugaunt,time_abcof
228 : integer :: n_start,n_end,mlotot,mlolotot,err
229 : integer :: mlot_d,mlolot_d,ilo,dir,l1,l2,ii,jj
230 : character(len=2) :: spin012(0:2)
231 : data spin012/' ', '.1', '.2'/
232 : character(len=6) :: filename
233 : real :: arg,hescale
234 : complex :: nsfactor,nsfactor_b,value
235 : real :: b1(3),b2(3),b3(3)
236 0 : real,allocatable :: radial1_ff(:,:,:,:,:)
237 0 : real,allocatable :: radial1_gg(:,:,:,:,:)
238 0 : real,allocatable :: radial1_fg(:,:,:,:,:)
239 0 : real,allocatable :: radial1_gf(:,:,:,:,:)
240 0 : real,allocatable :: radial1_flo(:,:,:,:,:)
241 0 : real,allocatable :: radial1_glo(:,:,:,:,:)
242 0 : real,allocatable :: radial1_lof(:,:,:,:,:)
243 0 : real,allocatable :: radial1_log(:,:,:,:,:)
244 0 : real,allocatable :: radial1_lolo(:,:,:,:,:)
245 :
246 : real,parameter :: bohrtocm=0.529177e-8
247 : real,parameter :: condquant=7.7480917e-5
248 : real :: dirfacs(3)
249 : complex,allocatable :: ubeffug(:,:,:)
250 : complex,allocatable :: ubeffdg(:,:,:)
251 : complex,allocatable :: dbeffug(:,:,:)
252 : complex,allocatable :: dbeffdg(:,:,:)
253 : complex,allocatable :: ubeffulog(:,:,:)
254 : complex,allocatable :: dbeffulog(:,:,:)
255 : complex,allocatable :: ulobeffug(:,:,:)
256 : complex,allocatable :: ulobeffdg(:,:,:)
257 : complex,allocatable :: ulobeffulog(:,:,:)
258 : integer :: spinloop1,spinloop2,ilop
259 :
260 0 : TYPE(t_usdus) :: usdus
261 0 : TYPE(t_mat) :: zMat(2),zzMat(2)
262 0 : TYPE(t_lapw) :: lapw
263 :
264 0 : IF(l_ss) CALL juDFT_error("spin-spiral not yet in this version"
265 0 : + ,calledby ="wann_updown")
266 :
267 : c-----initializations
268 0 : time_interstitial=0.0
269 0 : time_mmn=0.0
270 0 : time_total=0.0
271 0 : time_ujugaunt=0.0
272 0 : time_abcof=0.0
273 0 : time_rw=0.0
274 0 : time_symm=0.0
275 0 : time_film=0.0
276 0 : ngopr1(:)=1
277 :
278 0 : call cpu_time(time_total)
279 :
280 0 : l_p0 = .false.
281 0 : if (irank.eq.0) l_p0 = .true.
282 :
283 0 : sfp = 2* sqrt( pimach() )
284 0 : tpi = 2* pimach()
285 0 : lmd = lmaxd*(lmaxd+2)
286 :
287 : !!! should be changed in case the windows are really used
288 0 : nkpts = nkpt
289 :
290 : c-----read the input file to determine what to do
291 0 : call wann_read_inp(input,l_p0,wann)
292 :
293 0 : IF(wann%l_byenergy.AND.wann%l_byindex) CALL juDFT_error
294 0 : + ("byenergy.and.byindex",calledby ="wann_updown")
295 0 : IF(wann%l_byenergy.AND.wann%l_bynumber) CALL juDFT_error
296 0 : + ("byenergy.and.bynumber",calledby ="wann_updown")
297 0 : IF(wann%l_bynumber.AND.wann%l_byindex) CALL juDFT_error
298 0 : + ("bynumber.and.byindex",calledby ="wann_updown")
299 0 : IF(.NOT.(wann%l_bynumber.OR.wann%l_byindex.OR.wann%l_byenergy))
300 : & CALL juDFT_error("no rule to sort bands",calledby
301 0 : + ="wann_updown")
302 :
303 0 : efermi=ef
304 0 : if(.not.wann%l_fermi)efermi=0.0
305 :
306 : #ifdef CPP_MPI
307 0 : call MPI_BARRIER(mpi_comm,ierr)
308 : #endif
309 :
310 : c**************************************************************
311 : c for bzsym=.true.: determine mapping between kpts and w90kpts
312 : c**************************************************************
313 0 : if (wann%l_bzsym) then
314 0 : l_file=.false.
315 0 : inquire(file='w90kpts',exist=l_file)
316 0 : IF(.NOT.l_file) CALL juDFT_error
317 : + ("w90kpts not found, needed if bzsym",calledby
318 0 : + ="wann_updown")
319 0 : open(412,file='w90kpts',form='formatted')
320 0 : read(412,*)fullnkpts
321 0 : close(412)
322 0 : if(l_p0)print*,"fullnkpts=",fullnkpts
323 0 : IF(fullnkpts<nkpts) CALL juDFT_error("fullnkpts.lt.nkpts"
324 0 : + ,calledby ="wann_updown")
325 0 : allocate(irreduc(fullnkpts),mapkoper(fullnkpts))
326 0 : allocate(shiftkpt(3,fullnkpts))
327 0 : l_file=.false.
328 0 : inquire(file='kptsmap',exist=l_file)
329 0 : IF(.NOT.l_file) CALL juDFT_error
330 : + ("kptsmap not found, needed if bzsym",calledby
331 0 : + ="wann_updown")
332 0 : open(713,file='kptsmap')
333 0 : do i=1,fullnkpts
334 0 : read(713,*)kpt,irreduc(i),mapkoper(i),shiftkpt(:,i)
335 0 : IF(kpt/=i) CALL juDFT_error("kpt.ne.i",
336 0 : + calledby ="wann_updown")
337 0 : if(l_p0)print*,i,irreduc(i),mapkoper(i)
338 : enddo
339 0 : close(713)
340 0 : IF(MAXVAL(irreduc(:))/=nkpts) CALL juDFT_error
341 0 : + ("max(irreduc(:))/=nkpts",calledby ="wann_updown")
342 : else
343 0 : fullnkpts=nkpts
344 : endif
345 :
346 0 : nrec = 0
347 0 : if(l_p0)then
348 0 : write (*,*) 'fermi energy:',efermi
349 0 : write (*,*) 'emin,emax=',e1s,e2s
350 0 : write (*,*) 'nbasfcn =',nbasfcn
351 : endif
352 0 : nlotot = 0
353 0 : do n = 1, ntype
354 0 : do l = 1,nlo(n)
355 0 : nlotot = nlotot + neq(n) * ( 2*llo(l,n) + 1 )
356 : enddo
357 : enddo
358 :
359 : c*********************************************************
360 : cccccccccccccccc initialize the potential cccccccccccc
361 : c*********************************************************
362 :
363 0 : allocate (vz(nmzd,2,jspd))
364 0 : allocate (vr(jmtd,ntypd,jspd))
365 :
366 0 : vz = vTot%vacz
367 :
368 0 : do jspin = 1,jspins
369 0 : do n = 1, ntype
370 0 : do j = 1,jri(n)
371 0 : vr(j,n,jspin) = vTot%mt(j,0,n,jspin)
372 : enddo
373 : enddo
374 : enddo
375 :
376 : cccccccccccccccc end of the potential part ccccccccccc
377 0 : wannierspin=jspd
378 0 : if(l_soc) wannierspin=2
379 :
380 0 : allocate ( nv(wannierspin) )
381 0 : allocate ( k1(nvd,wannierspin),k2(nvd,wannierspin),
382 0 : & k3(nvd,wannierspin) )
383 0 : allocate ( ff(ntypd,jmtd,2,0:lmaxd,2) )
384 0 : allocate ( gg(ntypd,jmtd,2,0:lmaxd,2) )
385 0 : allocate ( usdus%us(0:lmaxd,ntypd,2) )
386 0 : allocate ( usdus%uds(0:lmaxd,ntypd,2) )
387 0 : allocate ( usdus%dus(0:lmaxd,ntypd,2) )
388 0 : allocate ( usdus%duds(0:lmaxd,ntypd,2) )
389 0 : allocate ( usdus%ddn(0:lmaxd,ntypd,2) )
390 0 : allocate ( usdus%ulos(nlod,ntypd,2) )
391 0 : allocate ( usdus%dulos(nlod,ntypd,2) )
392 0 : allocate ( usdus%uulon(nlod,ntypd,2) )
393 0 : allocate ( usdus%dulon(nlod,ntypd,2) )
394 0 : allocate ( usdus%uloulopn(nlod,nlod,ntypd,2) )
395 : c****************************************************
396 : c jspin=1
397 : c****************************************************
398 0 : jspin=1
399 :
400 : c...read number of bands and wannier functions from file proj
401 :
402 : c..reading the proj.1 / proj.2 / proj file
403 0 : l_proj=.false.
404 0 : do j=jspin,0,-1
405 0 : inquire(file=trim('proj'//spin012(j)),exist=l_proj)
406 0 : if(l_proj)then
407 0 : filename='proj'//spin012(j)
408 0 : exit
409 : endif
410 : enddo
411 :
412 0 : if(l_proj)then
413 0 : open (203,file=trim(filename),status='old')
414 0 : rewind (203)
415 0 : read (203,*) nwfs,numbands
416 0 : rewind (203)
417 0 : close (203)
418 : else
419 : CALL juDFT_error("no proj/proj.1/proj.2",
420 0 : + calledby ="wann_updown")
421 : endif
422 :
423 :
424 0 : jspin2=jspin
425 : if(l_soc .and. jspins.eq.1)jspin2=1
426 0 : jsp_start = jspin ; jsp_end = jspin
427 :
428 : cccccccccccc read in the eigenvalues and vectors cccccc
429 :
430 :
431 0 : do jspin=1,jspins
432 0 : CALL judft_error("TODO:wann_updown")
433 : ! call cdn_read0(eig_id,irank,isize,jspin,jspd,l_noco,
434 : ! < ello,evac,epar,wk,n_bands,n_size)
435 : enddo
436 :
437 0 : jspin=1
438 :
439 : c.. now we want to define the maximum number of the bands by all kpts
440 :
441 0 : nbnd = 0
442 :
443 0 : i_rec = 0 ; n_rank = 0
444 : c*************************************************************
445 : c..writing down the eig.1 and/or eig.2 files
446 : c*************************************************************
447 0 : if(l_p0)then
448 : call wann_write_eig(
449 : > eig_id,l_real,
450 : > lmaxd,ntypd,nlod,neigd,nvd,wannierspin,
451 : > isize,jspin,nbasfcn,nlotot,
452 : > l_ss,l_noco,nrec,fullnkpts,
453 : > wann%l_bzsym,wann%l_byindex,wann%l_bynumber,
454 : > wann%l_byenergy,
455 : > irreduc,odi,wann%band_min(jspin),
456 : > wann%band_max(jspin),
457 : > numbands,
458 : > e1s,e2s,efermi,wann%l_pauli,nkpts,
459 0 : < nbnd,kpoints,.false.,-1)
460 0 : if(odi%d1)then
461 0 : kpoints(:)=kpoints(:)*bmat(3,3)
462 : endif
463 : endif!l_p0
464 :
465 : ! nbnd is calculated for process zero and is sent here to the others
466 : #ifdef CPP_MPI
467 0 : if(l_p0)then
468 0 : do cpu_index=1,isize-1
469 0 : call MPI_SEND(nbnd,1,MPI_INTEGER,cpu_index,1,mpi_comm,ierr)
470 : enddo
471 : else
472 0 : call MPI_RECV(nbnd,1,MPI_INTEGER,0,1,mpi_comm,stt,ierr)
473 : endif
474 : #endif
475 :
476 0 : print*,"process: ",irank," nbnd= ",nbnd
477 : c##################################################################
478 :
479 0 : allocate ( mmn(nbnd,nbnd,fullnkpts) )
480 0 : mmn(:,:,:) = cmplx(0.,0.)
481 0 : if((l_soc.or.l_noco) .and. (jspin.eq.1))
482 0 : & allocate( socmmn(nbnd,nbnd,fullnkpts) )
483 :
484 0 : if(wann%l_nabla)then
485 0 : allocate ( nablamat(3,nbnd,nbnd,fullnkpts) )
486 0 : nablamat = cmplx(0.,0.)
487 : endif
488 :
489 0 : if(wann%l_surfcurr)then
490 0 : allocate ( surfcurr(3,nbnd,nbnd,fullnkpts) )
491 0 : surfcurr = cmplx(0.,0.)
492 : endif
493 :
494 0 : if(wann%l_socmat)then
495 0 : allocate ( hsomtx(2,2,nbnd,nbnd,fullnkpts) )
496 0 : allocate ( hsomtx_l(nbnd,nbnd,2,2) )
497 : endif
498 :
499 0 : write (*,*) 'nwfs=',nwfs
500 :
501 0 : allocate ( flo(ntypd,jmtd,2,nlod,2) )
502 :
503 : c***********************************************************
504 : c Calculate the radial wave functions for the two spins.
505 : c***********************************************************
506 0 : do jspin=1,wannierspin
507 0 : jspin2=jspin
508 0 : if(l_soc .and. jspins.eq.1)jspin2=1
509 0 : do n = 1,ntype
510 0 : do l = 0,lmax(n)
511 : call radfun(
512 : > l,n,jspin,epar(l,n,jspin2),vr(1,n,jspin2),atoms,
513 : < ff(n,:,:,l,jspin),gg(n,:,:,l,jspin),usdus,
514 0 : < nodeu,noded,wronk)
515 : enddo !l
516 0 : do ilo = 1, nlo(n)
517 :
518 : call radflo(
519 : > atoms,n,jspin,ello(:,:,jspin2),vr(1,n,jspin2),
520 : > ff(n,1:,1:,0:,jspin),gg(n,1:,1:,0:,jspin),mpi,
521 0 : < usdus,uuilon,duilon,ulouilopn,flo(n,:,:,:,jspin))
522 :
523 : enddo !ilo
524 : enddo !n
525 : enddo !jspin
526 :
527 : c****************************************************************
528 : c Calculate the radial overlap integrals.
529 : c****************************************************************
530 0 : if(wann%l_mmn0)then
531 0 : allocate( radial1_ff(2,2,0:lmaxd,0:lmaxd,ntype) )
532 0 : allocate( radial1_fg(2,2,0:lmaxd,0:lmaxd,ntype) )
533 0 : allocate( radial1_gf(2,2,0:lmaxd,0:lmaxd,ntype) )
534 0 : allocate( radial1_gg(2,2,0:lmaxd,0:lmaxd,ntype) )
535 :
536 0 : allocate( radial1_flo(2,2,0:lmaxd,1:nlod,ntype) )
537 0 : allocate( radial1_glo(2,2,0:lmaxd,1:nlod,ntype) )
538 0 : allocate( radial1_lof(2,2,1:nlod,0:lmaxd,ntype) )
539 0 : allocate( radial1_log(2,2,1:nlod,0:lmaxd,ntype) )
540 0 : allocate( radial1_lolo(2,2,1:nlod,1:nlod,ntype) )
541 :
542 0 : do n=1,ntype
543 0 : do l1=0,lmaxd
544 0 : do l2=0,lmaxd
545 0 : do spinloop1=1,2
546 0 : do spinloop2=1,2
547 : call wann_radovlp_integrals(
548 : > rmsh(:,n),jri(n),dx(n),
549 : > ff(n,:,:,l1,spinloop2),
550 : > ff(n,:,:,l2,spinloop1),
551 0 : < radial1_ff(spinloop2,spinloop1,l1,l2,n))
552 : call wann_radovlp_integrals(
553 : > rmsh(:,n),jri(n),dx(n),
554 : > ff(n,:,:,l1,spinloop2),
555 : > gg(n,:,:,l2,spinloop1),
556 0 : < radial1_fg(spinloop2,spinloop1,l1,l2,n))
557 : call wann_radovlp_integrals(
558 : > rmsh(:,n),jri(n),dx(n),
559 : > gg(n,:,:,l1,spinloop2),
560 : > ff(n,:,:,l2,spinloop1),
561 0 : < radial1_gf(spinloop2,spinloop1,l1,l2,n))
562 : call wann_radovlp_integrals(
563 : > rmsh(:,n),jri(n),dx(n),
564 : > gg(n,:,:,l1,spinloop2),
565 : > gg(n,:,:,l2,spinloop1),
566 0 : < radial1_gg(spinloop2,spinloop1,l1,l2,n))
567 : enddo !spinloop2
568 : enddo !spinloop1
569 : enddo !l2
570 : enddo !l1
571 : enddo !n
572 :
573 0 : do n=1,ntype
574 0 : do ilo=1,nlo(n)
575 0 : do ilop=1,nlo(n)
576 0 : do spinloop1=1,2
577 0 : do spinloop2=1,2
578 : call wann_radovlp_integrals(
579 : > rmsh(:,n),jri(n),dx(n),
580 : > flo(n,:,:,ilo,spinloop2),
581 : > flo(n,:,:,ilop,spinloop1),
582 0 : < radial1_lolo(spinloop2,spinloop1,ilo,ilop,n))
583 : enddo !spinloop2
584 : enddo !spinloop1
585 : enddo!ilop
586 : enddo !ilo
587 : enddo !n
588 :
589 0 : do n=1,ntype
590 0 : do ilo=1,nlo(n)
591 0 : do l2=0,lmaxd
592 0 : do spinloop1=1,2
593 0 : do spinloop2=1,2
594 : call wann_radovlp_integrals(
595 : > rmsh(:,n),jri(n),dx(n),
596 : > flo(n,:,:,ilo,spinloop2),
597 : > ff(n,:,:,l2,spinloop1),
598 0 : < radial1_lof(spinloop2,spinloop1,ilo,l2,n))
599 : enddo !spinloop2
600 : enddo !spinloop1
601 : enddo!l2
602 : enddo !ilo
603 : enddo !n
604 :
605 0 : do n=1,ntype
606 0 : do ilo=1,nlo(n)
607 0 : do l2=0,lmaxd
608 0 : do spinloop1=1,2
609 0 : do spinloop2=1,2
610 : call wann_radovlp_integrals(
611 : > rmsh(:,n),jri(n),dx(n),
612 : > flo(n,:,:,ilo,spinloop2),
613 : > gg(n,:,:,l2,spinloop1),
614 0 : < radial1_log(spinloop2,spinloop1,ilo,l2,n))
615 : enddo !spinloop2
616 : enddo !spinloop1
617 : enddo!l2
618 : enddo !ilo
619 : enddo !n
620 :
621 0 : do n=1,ntype
622 0 : do ilo=1,nlo(n)
623 0 : do l2=0,lmaxd
624 0 : do spinloop1=1,2
625 0 : do spinloop2=1,2
626 : call wann_radovlp_integrals(
627 : > rmsh(:,n),jri(n),dx(n),
628 : > ff(n,:,:,l2,spinloop2),
629 : > flo(n,:,:,ilo,spinloop1),
630 0 : < radial1_flo(spinloop2,spinloop1,l2,ilo,n))
631 : enddo !spinloop2
632 : enddo !spinloop1
633 : enddo!l2
634 : enddo !ilo
635 : enddo !n
636 :
637 0 : do n=1,ntype
638 0 : do ilo=1,nlo(n)
639 0 : do l2=0,lmaxd
640 0 : do spinloop1=1,2
641 0 : do spinloop2=1,2
642 : call wann_radovlp_integrals(
643 : > rmsh(:,n),jri(n),dx(n),
644 : > gg(n,:,:,l2,spinloop2),
645 : > flo(n,:,:,ilo,spinloop1),
646 0 : < radial1_glo(spinloop2,spinloop1,l2,ilo,n))
647 : enddo !spinloop2
648 : enddo !spinloop1
649 : enddo!l2
650 : enddo !ilo
651 : enddo !n
652 :
653 : endif
654 :
655 :
656 0 : jspin=1
657 :
658 : c****************************************************************
659 : c compute matrix elements of exchange field
660 : c****************************************************************
661 : #ifdef CPP_TOPO
662 : if(wann%l_perpmag)then
663 : allocate( perpmag(nbnd,nbnd,fullnkpts) )
664 : perpmag=0.0
665 : allocate( ubeffug(0:lmd,0:lmd,1:ntype) )
666 : allocate( ubeffdg(0:lmd,0:lmd,1:ntype) )
667 : allocate( dbeffug(0:lmd,0:lmd,1:ntype) )
668 : allocate( dbeffdg(0:lmd,0:lmd,1:ntype) )
669 : allocate( ubeffulog(0:lmd,0:15,1:ntype) )
670 : allocate( dbeffulog(0:lmd,0:15,1:ntype) )
671 : allocate( ulobeffug(0:lmd,0:15,1:ntype) )
672 : allocate( ulobeffdg(0:lmd,0:15,1:ntype) )
673 : allocate( ulobeffulog(0:15,0:15,1:ntype) )
674 : call wann_gabeffgagaunt(
675 : > neq,mlh,nlhd,nlh,ntypsy,llh,lmax,
676 : > nmem,ntype,ntypd,bbmat,bmat,
677 : > nlod,nlo,llo,flo,ff,gg,jri,rmsh,dx,jmtd,
678 : > lmaxd,lmd,clnu,
679 : < ubeffug,ubeffdg,dbeffug,dbeffdg,ubeffulog,
680 : < dbeffulog,
681 : < ulobeffug,ulobeffdg,ulobeffulog)
682 : endif
683 : #endif
684 :
685 :
686 0 : i_rec = 0 ; n_rank = 0
687 :
688 : c****************************************************************
689 : c.. loop by kpoints starts! each may be a separate task
690 : c****************************************************************
691 0 : do 10 ikpt = wann%ikptstart,fullnkpts ! loop by k-points starts
692 0 : kptibz=ikpt
693 0 : if(wann%l_bzsym) kptibz=irreduc(ikpt)
694 0 : if(wann%l_bzsym) oper=mapkoper(ikpt)
695 :
696 0 : i_rec = i_rec + 1
697 0 : if (mod(i_rec-1,isize).eq.irank) then
698 :
699 : c****************************************************************
700 : c Obtain eigenvectors for the two spin channels.
701 : c****************************************************************
702 :
703 0 : allocate ( eigg(neigd,2) )
704 0 : call cpu_time(delta)
705 0 : do jspin=1,wannierspin
706 0 : if(jspin.eq.1)nrec=0
707 0 : if(jspin.eq.2)nrec=nkpts
708 0 : if(l_noco)nrec=0
709 0 : n_start=1
710 0 : n_end=neigd
711 :
712 0 : zzMat(jspin)%l_real = l_real
713 0 : zzMat(jspin)%matsize1 = nbasfcn
714 0 : zzMat(jspin)%matsize2 = neigd
715 0 : IF(l_real) THEN
716 0 : ALLOCATE (zzMat(jspin)%data_r(zzMat(jspin)%matsize1,
717 0 : + zzMat(jspin)%matsize2))
718 : ELSE
719 0 : ALLOCATE (zzMat(jspin)%data_c(zzMat(jspin)%matsize1,
720 0 : + zzMat(jspin)%matsize2))
721 : END IF
722 :
723 : ! CALL cdn_read(
724 : ! > eig_id,
725 : ! > nvd,jspd,irank,isize,kptibz,jspin,nbasfcn,
726 : ! > l_ss,l_noco,neigd,n_start,n_end,
727 : ! < nmat,nv,ello,evdu,epar,
728 : ! < k1,k2,k3,bkpt,wk,nbands,eigg(:,jspin),zzMat(jspin))
729 0 : call judft_error("cdn_read in wann_updown")
730 0 : nrec=0
731 : enddo !jspin
732 0 : call cpu_time(delta1)
733 0 : time_rw=time_rw+delta1-delta
734 0 : nslibd = 0
735 :
736 0 : jspin=1
737 :
738 : c...we work only within the energy window
739 :
740 0 : DO jspin = 1, wannierspin
741 0 : zMat(jspin)%l_real = zzMat(jspin)%l_real
742 0 : zMat(jspin)%matsize1 = zzMat(jspin)%matsize1
743 0 : zMat(jspin)%matsize2 = zzMat(jspin)%matsize2
744 0 : IF (zzMat(jspin)%l_real) THEN
745 0 : IF(.not.allocated(zmat(jspin)%data_r))
746 0 : + ALLOCATE (zMat(jspin)%data_r(zMat(jspin)%matsize1,
747 0 : + zMat(jspin)%matsize2))
748 0 : zMat(jspin)%data_r = 0.0
749 : ELSE
750 0 : IF(.not.allocated(zMat(jspin)%data_c))
751 0 : + ALLOCATE (zMat(jspin)%data_c(zMat(jspin)%matsize1,
752 0 : + zMat(jspin)%matsize2))
753 0 : zMat(jspin)%data_c = CMPLX(0.0,0.0)
754 : END IF
755 : END DO
756 :
757 0 : eig(:) = 0.
758 :
759 0 : do jspin=1,wannierspin
760 0 : nslibd=0
761 : !print*,"bands used:"
762 0 : do i = 1,nbands
763 : if ((eigg(i,1).ge.e1s.and.nslibd.lt.numbands.
764 : & and.wann%l_bynumber)
765 : & .or.(eigg(i,1).ge.e1s.and.eigg(i,1).le.e2s.
766 : & and.wann%l_byenergy)
767 0 : & .or.(i.ge.wann%band_min(jspin).and.
768 : & (i.le.wann%band_max(jspin))
769 0 : & .and.wann%l_byindex))then
770 :
771 : !print*,i
772 0 : nslibd = nslibd + 1
773 0 : eig(nslibd) = eigg(i,1)
774 0 : if(l_noco)then
775 0 : funbas=nv(1)+nlotot
776 0 : funbas=funbas+nv(2)+nlotot
777 : else
778 0 : funbas=nv(jspin)+nlotot
779 : endif
780 0 : IF (zzMat(jspin)%l_real) THEN
781 0 : do j = 1,funbas
782 0 : zMat(jspin)%data_r(j,nslibd) = zzMat(jspin)%data_r(j,i)
783 : enddo
784 : ELSE
785 0 : do j = 1,funbas
786 0 : zMat(jspin)%data_c(j,nslibd) = zzMat(jspin)%data_c(j,i)
787 : enddo
788 : END IF
789 : endif
790 : enddo
791 : enddo
792 0 : jspin=1
793 :
794 : c***********************************************************
795 : c rotate the wavefunction
796 : c***********************************************************
797 0 : do jspin=1,wannierspin
798 0 : if (wann%l_bzsym.and.oper.ne.1) then !rotate bkpt
799 : call wann_kptsrotate(
800 : > natd,nlod,llod,
801 : > ntypd,nlo,llo,invsat,
802 : > l_noco,l_soc,
803 : > ntype,neq,nlotot,
804 : > jspin,
805 : > oper,nop,mrot,nvd,nv,
806 : > shiftkpt(:,ikpt),
807 : > tau,
808 : x bkpt,k1(:,:),k2(:,:),k3(:,:),
809 0 : x zMat(jspin),nsfactor)
810 : else
811 0 : nsfactor=cmplx(1.0,0.0)
812 : endif
813 : enddo !jspin
814 : c print*,"bkpt1=",bkpt
815 :
816 : c******************************************************************
817 :
818 : c...the overlap matrix Mmn which is computed for each k- and b-point
819 :
820 0 : noccbd = nslibd
821 :
822 0 : allocate ( acof(noccbd,0:lmd,natd,wannierspin),
823 0 : & bcof(noccbd,0:lmd,natd,wannierspin),
824 0 : & ccof(-llod:llod,noccbd,nlod,natd,wannierspin))
825 :
826 0 : acof(:,:,:,:) = cmplx(0.,0.) ; bcof(:,:,:,:) = cmplx(0.,0.)
827 0 : ccof(:,:,:,:,:) = cmplx(0.,0.)
828 :
829 : c...generation the A,B,C coefficients in the spheres
830 : c...for the lapws and local orbitals, summed by the basis functions
831 :
832 0 : ALLOCATE(lapw%k1(SIZE(k1,1),SIZE(k1,2)))
833 0 : ALLOCATE(lapw%k2(SIZE(k1,1),SIZE(k1,2)))
834 0 : ALLOCATE(lapw%k3(SIZE(k1,1),SIZE(k1,2)))
835 0 : lapw%nv = nv
836 0 : lapw%nmat = nmat
837 0 : lapw%k1 = k1
838 0 : lapw%k2 = k2
839 0 : lapw%k3 = k3
840 : ! I think the other variables of lapw are not needed here.
841 :
842 0 : call cpu_time(delta)
843 0 : do jspin=1,wannierspin
844 :
845 : CALL abcof(input,atoms,sym,cell,lapw,noccbd,usdus,
846 : + noco,jspin,oneD,acof(1:,0:,1:,jspin),
847 : + bcof(1:,0:,1:,jspin),ccof(-llod:,1:,1:,1:,jspin),
848 0 : + zMat(jspin))
849 :
850 : call wann_abinv(atoms,
851 : X acof(:,0:,:,jspin),bcof(:,0:,:,jspin),
852 0 : < ccof(-llod:,:,:,:,jspin))
853 : enddo !jspin
854 :
855 0 : DEALLOCATE(lapw%k1,lapw%k2,lapw%k3)
856 :
857 0 : jspin=1
858 0 : call cpu_time(delta1)
859 0 : time_abcof=time_abcof+delta1-delta
860 :
861 0 : if(wann%l_socmat)then
862 : call wann_socmat(
863 : > mpi,enpara,input,noco,atoms,
864 : > lmaxd,ntypd,nlod,natd,noccbd,
865 : > llod,jmtd,jspd,nlhd,neq,
866 : > ntype,theta,phi,jspins,irank,
867 : > jri,lmax,dx,rmsh,epar,ello,nlo,llo,
868 : > l_dulo,ulo_der,vTot%mt,
869 : > acof,bcof,ccof,
870 0 : < hsomtx_l)
871 0 : do i1 = 1,2
872 0 : do i2 = 1,2
873 0 : do i = 1,nbnd
874 0 : do j = 1,nbnd
875 0 : hsomtx(i1,i2,i,j,ikpt) = hsomtx_l(i,j,i1,i2)
876 : enddo
877 : enddo
878 : enddo
879 : enddo
880 : endif
881 :
882 : #ifdef CPP_TOPO
883 :
884 : if(wann%l_perpmag)then
885 : call wann_perpmag(
886 : > nslibd,neq,mlh,nlhd,nlh,ntypsy,llh,lmax,
887 : > nmem,ntype,ntypd,bbmat,bmat,
888 : > nlod,nlo,llo,flo,ff,gg,jri,rmsh,dx,jmtd,
889 : > lmaxd,lmd,clnu,
890 : > ubeffug,ubeffdg,dbeffug,dbeffdg,ubeffulog,
891 : > dbeffulog,
892 : > ulobeffug,ulobeffdg,ulobeffulog,
893 : > acof,bcof,
894 : < perpmag(:,:,ikpt))
895 : endif
896 :
897 : if(wann%l_surfcurr)then
898 : call wann_surfcurr_int_updown(
899 : > nv2d,jspin,odi,ods,n3d,nmzxyd,n2d,ntypsd,
900 : > ntype,lmaxd,jmtd,ntypd,natd,nmzd,neq,nq3,nvac,
901 : > nmz,nmzxy,nq2,nop,nop2,volint,film,slice,symor,
902 : > invs,invs2,z1,delz,ngopr,ntypsy,jri,
903 : > pos,taual,zatom,rmt,
904 : > lmax,mrot,tau,rmsh,invtab,amat,bmat,bbmat,ikpt,
905 : > nvd,nv(jspin),bkpt,omtil,
906 : > k1(:,jspin),k2(:,jspin),k3(:,jspin),
907 : > nslibd,nbasfcn,neigd,zMat,
908 : > dirfacs,
909 : > surfcurr(:,:,:,ikpt))
910 : call wann_surfcurr_updown(
911 : > dirfacs,amat,
912 : > ntype,lmaxd,lmax,natd,
913 : > neq,noccbd,lmd,natd,llod,nlod,
914 : > nlo,llo,
915 : > acof,bcof,ccof,
916 : > us,dus,duds,uds,
917 : > ulos,dulos,
918 : > rmt,pos,
919 : & surfcurr(:,:,:,ikpt))
920 : write(6,*)"dirfacs=",dirfacs
921 : endif
922 :
923 : if(wann%l_nabla)then
924 : call wann_nabla_updown(
925 : > jspin,ntype,lmaxd,lmax,natd,
926 : > jmtd,jri,rmsh,dx,neq,
927 : > noccbd,lmd,natd,llod,nlod,
928 : > ff(:,:,:,:,1),gg(:,:,:,:,1),
929 : > ff(:,:,:,:,2),gg(:,:,:,:,2),
930 : > acof(:,:,:,1),bcof(:,:,:,1),
931 : > ccof(:,:,:,:,1),
932 : > acof(:,:,:,2),bcof(:,:,:,2),
933 : > ccof(:,:,:,:,2),
934 : & nablamat(:,:,:,ikpt))
935 :
936 : addnoco=0
937 : do 41 i = n_rank+1,nv(jspin),n_size
938 : b1(1)=bkpt(1)+k1(i,jspin)
939 : b1(2)=bkpt(2)+k2(i,jspin)
940 : b1(3)=bkpt(3)+k3(i,jspin)
941 : b2(1)=b1(1)*bmat(1,1)+b1(2)*bmat(2,1)+b1(3)*bmat(3,1)
942 : b2(2)=b1(1)*bmat(1,2)+b1(2)*bmat(2,2)+b1(3)*bmat(3,2)
943 : b2(3)=b1(1)*bmat(1,3)+b1(2)*bmat(2,3)+b1(3)*bmat(3,3)
944 : do 42 j = n_rank+1,nv(jspin),n_size
945 : c b1(1)=bkpt(1)+k1(j,jspin)
946 : c b1(2)=bkpt(2)+k2(j,jspin)
947 : c b1(3)=bkpt(3)+k3(j,jspin)
948 : c b3(1)=b1(1)*bmat(1,1)+b1(2)*bmat(2,1)+b1(3)*bmat(3,1)
949 : c b3(2)=b1(1)*bmat(1,2)+b1(2)*bmat(2,2)+b1(3)*bmat(3,2)
950 : c b3(3)=b1(1)*bmat(1,3)+b1(2)*bmat(2,3)+b1(3)*bmat(3,3)
951 : c--> determine index and phase factor
952 : i1 = k1(j,jspin) - k1(i,jspin)
953 : i2 = k2(j,jspin) - k2(i,jspin)
954 : i3 = k3(j,jspin) - k3(i,jspin)
955 : in = ig(i1,i2,i3)
956 : if (in.eq.0) goto 42
957 : phase = rgphs(i1,i2,i3)
958 : phasust = cmplx(phase,0.0)*ustep(in)
959 :
960 : do m = 1,nslibd
961 : do n = 1,nslibd
962 : do dir=1,3
963 : IF (zMat%l_real) THEN
964 : zzConjgTemp = cmplx(zMat(1)%data_r(i+addnoco,m) *
965 : + zMat(2)%data_r(j+addnoco,n), 0.0)
966 : ELSE
967 : zzConjgTemp = zMat(1)%data_c(i+addnoco,m) *
968 : + CONJG(zMat(2)%data_c(j+addnoco,n))
969 : END IF
970 : value=phasust*zzConjgTemp
971 : nablamat(dir,m,n,ikpt) = nablamat(dir,m,n,ikpt) -
972 : + value*b2(dir)
973 : enddo
974 : enddo
975 : enddo
976 :
977 : 42 continue
978 : 41 continue
979 : endif
980 : #endif
981 :
982 : c***************************************************
983 : c Calculate matrix elements of the Pauli matrix.
984 : c***************************************************
985 0 : if(wann%l_mmn0)then
986 :
987 : c------> interstitial contribution to updown.mmn0-matrix
988 0 : addnoco=0
989 0 : if(l_noco)then
990 0 : addnoco=nv(1)+nlotot
991 : endif
992 :
993 : c$$$ do i = n_rank+1,nv(jspin),n_size
994 : c$$$ do 22 j = n_rank+1,nv(jspin),n_size
995 : c$$$
996 : c$$$c--> determine index and phase factor
997 : c$$$ i1 = k1(j,jspin) - k1(i,jspin)
998 : c$$$ i2 = k2(j,jspin) - k2(i,jspin)
999 : c$$$ i3 = k3(j,jspin) - k3(i,jspin)
1000 : c$$$ in = ig(i1,i2,i3)
1001 : c$$$ if (in.eq.0) goto 22
1002 : c$$$ phase = rgphs(i1,i2,i3)
1003 : c$$$ phasust = cmplx(phase,0.0)*ustep(in)
1004 : c$$$ do m = 1,nslibd
1005 : c$$$ do n = 1,nslibd
1006 : c$$$#if ( !defined(CPP_INVERSION) || defined(CPP_SOC) )
1007 : c$$$ mmn(m,n,ikpt) =
1008 : c$$$ = mmn(m,n,ikpt) +
1009 : c$$$ + phasust*z(i,m,1)*conjg(z(j+addnoco,n,2))
1010 : c$$$#else
1011 : c$$$ mmn(m,n,ikpt) =
1012 : c$$$ = mmn(m,n,ikpt) +
1013 : c$$$ + phasust*cmplx(z(i,m,1)*z(j+addnoco,n,2),0.0)
1014 : c$$$#endif
1015 : c$$$ enddo !n
1016 : c$$$ enddo !m
1017 : c$$$
1018 : c$$$ 22 continue
1019 : c$$$ enddo !i
1020 :
1021 0 : addnoco2 = addnoco ! TODO: correct for addnoco2
1022 : call wann_mmkb_int(
1023 : > cmplx(1.,0.),addnoco,addnoco2,
1024 : > nvd,k1d,k2d,k3d,
1025 : > n3d,k1(:,jspin),k2(:,jspin),k3(:,jspin),
1026 : > nv(jspin),neigd,nbasfcn,zMat(1),nslibd,
1027 : > k1(:,jspin),k2(:,jspin),k3(:,jspin),
1028 : > nv(jspin),zMat(2),nslibd,
1029 : > nbnd,
1030 : > rgphs,ustep,ig,(/ 0,0,0 /),
1031 0 : < mmn(:,:,ikpt))
1032 :
1033 :
1034 : c---> spherical contribution to updown.mmn0-matrix
1035 : call wann_mmk0_updown_sph(
1036 : > l_noco,alph,beta,
1037 : > llod,noccbd,nlod,natd,ntypd,lmaxd,lmax,lmd,
1038 : > ntype,neq,nlo,llo,
1039 : > radial1_ff,radial1_gg,radial1_fg,radial1_gf,
1040 : > radial1_flo,radial1_glo,
1041 : > radial1_lof,radial1_log,
1042 : > radial1_lolo,
1043 : > acof(1:noccbd,:,:,:),
1044 : > bcof(1:noccbd,:,:,:),ccof(:,1:noccbd,:,:,:),
1045 : > usdus%ddn,usdus%uulon,usdus%dulon,usdus%uloulopn,
1046 0 : = mmn(1:noccbd,1:noccbd,ikpt))
1047 :
1048 : c---> vacuum contribution to updown.mmn0-matrix
1049 :
1050 : if (film .and. .not.odi%d1) then
1051 :
1052 : c call wann_mmk0_vac(
1053 : c > z1,nmzd,nv2d,k1d,k2d,k3d,n3d,nvac,
1054 : c > ig,nmz,delz,ig2,area,bmat,
1055 : c > bbmat,evac(:,jspin),bkpt,vz(:,:,jspin2),
1056 : c > nslibd,jspin,k1,k2,k3,wannierspin,nvd,
1057 : c > nbasfcn,neigd,z,nv,omtil,
1058 : c < mmn(1:noccbd,1:noccbd,ikpt))
1059 :
1060 : elseif (odi%d1) then
1061 :
1062 : c call wann_mmk0_od_vac(
1063 : c > z1,nmzxyd,nmzd,nv2d,k1d,k2d,k3d,n2d,n3d,
1064 : c > ig,nmzxy,nmz,delz,ig2,odi%n2d,
1065 : c > bbmat,evac(1,jspin),bkpt,odi%M,odi%mb,
1066 : c > vz(:,1,jspin2),odi,
1067 : c > nslibd,jspin,k1,k2,k3,wannierspin,nvd,area,
1068 : c > nbasfcn,neigd,z,nv,sk2,phi2,omtil,qss,
1069 : c < mmn(1:noccbd,1:noccbd,ikpt))
1070 :
1071 : endif
1072 : endif !l_mmn0
1073 :
1074 0 : deallocate ( acof,bcof,ccof,eigg )
1075 :
1076 0 : write (*,*) 'nslibd=',nslibd
1077 :
1078 :
1079 : endif ! loop by processors
1080 :
1081 0 : 10 continue ! end of cycle by the k-points
1082 :
1083 : #ifdef CPP_MPI
1084 0 : call MPI_BARRIER(mpi_comm,ierr)
1085 : #endif
1086 : c******************************************************
1087 : c Write down the projections.
1088 : c******************************************************
1089 : 5 continue
1090 :
1091 0 : if(wann%l_nabla)then
1092 0 : hescale=tpi*condquant/bohrtocm/omtil
1093 0 : hescale=sqrt(hescale)
1094 0 : nablamat=nablamat*hescale
1095 : call wann_write_nabla(
1096 : > mpi_comm,l_p0,'nabl.updown',
1097 : > 'Matrix elements of nabla operator',
1098 : > nbnd,fullnkpts,nbnd,
1099 : > irank,isize,
1100 0 : < nablamat)
1101 : endif
1102 :
1103 0 : if(wann%l_socmat)then
1104 : call wann_write_hsomtx(
1105 : > mpi_comm,l_p0,'WF1.hsomtx',
1106 : > 'Matrix elements of spin-orbit coupling',
1107 : > 2,2,nbnd,nbnd,fullnkpts,
1108 : > irank,isize,
1109 0 : < hsomtx)
1110 : endif
1111 :
1112 0 : if(wann%l_mmn0)then
1113 : call wann_write_amn(
1114 : > mpi_comm,l_p0,'updown.mmn0',
1115 : > 'Overlaps of the wavefunct. at the same kpoint',
1116 : > nbnd,fullnkpts,nbnd,
1117 : > irank,isize,.false.,
1118 0 : < mmn,wann%l_unformatted)
1119 : endif !l_soc and l_mmn0
1120 :
1121 0 : if(wann%l_perpmag)then
1122 : call wann_write_amn(
1123 : > mpi_comm,l_p0,'updown.perpmag',
1124 : > 'Overlaps of the wavefunct. with Beff',
1125 : > nbnd,fullnkpts,nbnd,
1126 : > irank,isize,.true.,
1127 0 : < perpmag,wann%l_unformatted)
1128 : endif !l_soc and l_mmn0
1129 :
1130 0 : if(wann%l_surfcurr)then
1131 0 : surfcurr=conjg(surfcurr/cmplx(0.0,2.0))
1132 0 : hescale=tpi*condquant/bohrtocm/omtil
1133 0 : hescale=sqrt(hescale)
1134 0 : surfcurr=surfcurr*hescale
1135 : call wann_write_nabla(
1136 : > mpi_comm,l_p0,'surfcurr.updown',
1137 : > 'Surface currents',
1138 : > nbnd,fullnkpts,nbnd,
1139 : > irank,isize,
1140 0 : < surfcurr)
1141 : endif
1142 :
1143 0 : if( allocated (nablamat) ) deallocate( nablamat )
1144 0 : deallocate ( mmn )
1145 :
1146 0 : deallocate(flo)
1147 :
1148 0 : deallocate ( vr,vz,nv,k1,k2,k3 )
1149 0 : deallocate ( ff,gg )
1150 0 : if (wann%l_bzsym) deallocate(irreduc,mapkoper)
1151 :
1152 0 : call cpu_time(delta)
1153 0 : time_total=delta-time_total
1154 0 : write(6,*)"time_total=",time_total
1155 0 : write(6,*)"time_mmn=",time_mmn
1156 0 : write(6,*)"time_interstitial=",time_interstitial
1157 0 : write(6,*)"time_ujugaunt=",time_ujugaunt
1158 0 : write(6,*)"time_abcof=",time_abcof
1159 0 : write(6,*)"time_symm=",time_symm
1160 0 : write(6,*)"time_rw=",time_rw
1161 0 : write(6,*)"time_film=",time_film
1162 :
1163 : #ifdef CPP_MPI
1164 0 : call MPI_BARRIER(mpi_comm,ierr)
1165 : #endif
1166 :
1167 : #ifdef CPP_MPI
1168 0 : call MPI_BARRIER(mpi_comm,ierr)
1169 : #endif
1170 :
1171 :
1172 0 : END SUBROUTINE wann_updown
1173 : END MODULE m_wann_updown
|
