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