Line data Source code
1 : c*******************************************c
2 : c Set up uHu matrix necessary for c
3 : c Wannier-based calc. of orbital moment c
4 : c*******************************************c
5 : c keyword is 'matrixuhu' in wann_inp c
6 : c*******************************************c
7 : c uHu = < u_{k+b1} | H_{k} | u_{k+b2} > c
8 : c c
9 : c Contributions to Hamiltonian: c
10 : c (i) interstitial c
11 : c (ii) muffin tin (a) spherical c
12 : c (b) non-sph. c
13 : c (c) SOC c
14 : c (iii) vacuum c
15 : c*******************************************c
16 : c J.-P. Hanke, Dec. 2015 c
17 : c*******************************************c
18 : MODULE m_wann_uHu
19 : USE m_juDFT
20 : #ifdef CPP_MPI
21 : use mpi
22 : #endif
23 : CONTAINS
24 0 : SUBROUTINE wann_uHu(
25 : > stars,vacuum,atoms,sphhar,input,kpts,sym,fmpi,
26 : > banddos ,noco,nococonv,cell,vTot,wann,enpara,
27 0 : > eig_idList,l_real,l_dulo,l_noco,l_ss,lmaxd,ntypd,
28 : > neigd,natd,nop,nvd,jspd,llod,nlod,ntype,
29 0 : > omtil,nlo,llo,lapw_l,invtab,mrot,ngopr,neq,lmax,
30 0 : > invsat,invsatnr,nkpt,taual,rmt,amat,bmat,bbmat,alph,
31 0 : > beta,qss,sk2,phi2,irank,isize,n3d,nmzxyd,nmzd,
32 0 : > jmtd,nlhd,nq3,nvac,invs,invs2,film,nlh,jri,ntypsd,
33 0 : > ntypsy,jspins,nkptd,dx,n2d,rmsh,e1s,e2s,ulo_der,
34 0 : > ustep,ig,k1d,k2d,k3d,rgphs,slice,kk,nnne,
35 0 : > z1,nv2d,nmzxy,nmz,delz,ig2,area,tau,zatom,nq2,kv2,nop2,
36 : > volint,symor,pos,ef,l_soc,
37 0 : > memd,lnonsph,clnu,lmplmd,mlh,nmem,llh,lo1l,
38 : > theta,phi,
39 : > l_ms,l_sgwf,l_socgwf,aux_latt_const,
40 0 : > param_file,param_vec,nparampts,param_alpha,l_dim)
41 :
42 : use m_types
43 : use m_wann_mmnk_symm
44 : use m_wann_rw_eig
45 : use m_abcof
46 : use m_radfun
47 : use m_radflo
48 : use m_cdnread
49 : use m_constants, only : pimach
50 : use m_wann_projmethod
51 : use m_wann_abinv
52 : use m_wann_kptsrotate
53 : ! use m_wann_read_inp !Call wann_read_inp in fleur_init
54 : use m_wann_maxbnd
55 : use m_wann_uHu_tlmplm
56 : use m_wann_uHu_sph
57 : use m_wann_uHu_int
58 : use m_wann_uHu_soc
59 : use m_wann_uHu_vac
60 : use m_wann_uHu_util
61 : use m_wann_uHu_commat
62 : use m_wann_write_uHu
63 : USE m_eig66_io
64 :
65 : IMPLICIT NONE
66 :
67 : #ifdef CPP_MPI
68 : integer ierr(3)
69 : integer cpu_index
70 : integer stt(MPI_STATUS_SIZE)
71 : #endif
72 :
73 :
74 : TYPE(t_stars),INTENT(IN) :: stars
75 : TYPE(t_vacuum),INTENT(IN) :: vacuum
76 : TYPE(t_atoms),INTENT(IN) :: atoms
77 : TYPE(t_sphhar),INTENT(IN) :: sphhar
78 : TYPE(t_input),INTENT(IN) :: input
79 : TYPE(t_kpts),INTENT(IN) :: kpts
80 : TYPE(t_sym),INTENT(IN) :: sym
81 : TYPE(t_mpi),INTENT(IN) :: fmpi
82 : TYPE(t_banddos),INTENT(IN) :: banddos
83 :
84 : TYPE(t_noco),INTENT(IN) :: noco
85 : TYPE(t_nococonv),INTENT(IN) :: nococonv
86 : TYPE(t_cell),INTENT(IN) :: cell
87 : TYPE(t_potden),INTENT(IN) :: vTot
88 : TYPE(t_wann),INTENT(INOUT) :: wann
89 : TYPE(t_enpara),INTENT(IN) :: enpara
90 :
91 : c ..scalar arguments..
92 : character(len=20),intent(in) :: param_file
93 :
94 : INTEGER, INTENT (IN) :: eig_idList(wann%nparampts)
95 : logical, intent (in) :: invs,invs2,film,slice,symor
96 : logical, intent (in) :: l_real,l_noco,l_ss,l_soc
97 : logical, intent (in) :: l_ms,l_sgwf,l_socgwf
98 : integer, intent (in) :: lmaxd,ntypd,neigd,nkptd,kk,nnne
99 : integer, intent (in) :: natd,nop,nvd,jspd,nq2,nop2
100 : integer, intent (in) :: llod,nlod,ntype,n3d,n2d
101 : integer, intent (in) :: nmzxyd,nmzd,jmtd,nlhd,nq3,nvac
102 : integer, intent (in) :: ntypsd,jspins,k1d,k2d,k3d
103 : integer, intent (in) :: irank,isize,nv2d,nmzxy,nmz
104 : integer, intent (in) :: memd,lmplmd,nparampts
105 : real, intent (in) :: omtil,e1s,e2s,delz,area,z1,volint
106 : real, intent (in) :: ef,theta,phi,aux_latt_const
107 :
108 : c ..array arguments..
109 : logical, intent (in) :: l_dulo(nlod,ntypd)
110 : logical, intent (in) :: l_dim(3)
111 : integer, intent (in) :: ig(-k1d:k1d,-k2d:k2d,-k3d:k3d)
112 : integer, intent (in) :: nlh(ntypsd),jri(ntypd),ntypsy(natd)
113 : integer, intent (in) :: nlo(ntypd),llo(nlod,ntypd),lapw_l(ntypd)
114 : integer, intent (in) :: invtab(nop),mrot(3,3,nop),ngopr(natd)
115 : integer, intent (in) :: neq(ntypd),lmax(ntypd)
116 : integer, intent (in) :: invsat(natd),invsatnr(natd),nkpt
117 : integer, intent (in) :: ulo_der(nlod,ntypd),ig2(n3d),kv2(2,n2d)
118 : integer, intent (in) :: mlh(memd,0:nlhd,ntypsd)
119 : integer, intent (in) :: nmem(0:nlhd,ntypsd)
120 : integer, intent (in) :: llh(0:nlhd,ntypsd),lnonsph(ntypd)
121 : integer, intent (in) :: lo1l(0:llod,ntypd)
122 : complex, intent (in) :: rgphs(-k1d:k1d,-k2d:k2d,-k3d:k3d)
123 : real, intent (in) :: taual(3,natd),rmt(ntypd),dx(ntypd)
124 : real, intent (in) :: amat(3,3),bmat(3,3),bbmat(3,3)
125 : real, intent (in) :: rmsh(jmtd,ntypd),tau(3,nop),zatom(ntype)
126 : real, intent (in) :: alph(ntypd),beta(ntypd),qss(3)
127 : real, intent (in) :: pos(3,natd),sk2(n2d),phi2(n2d)
128 : real, intent (in) :: param_vec(3,nparampts)
129 : real, intent (in) :: param_alpha(ntypd,nparampts)
130 : complex, intent (in) :: ustep(n3d)
131 : complex, intent (in) :: clnu(memd,0:nlhd,ntypsd)
132 :
133 : c ..allocatable arrays..
134 : ! integer, allocatable :: kveclo(:) , nv(:)
135 : ! integer, allocatable :: kveclo_b(:) , nv_b(:)
136 : ! integer, allocatable :: kveclo_b2(:), nv_b2(:)
137 : ! integer, allocatable :: k1(:,:) , k2(:,:) , k3(:,:)
138 : ! integer, allocatable :: k1_b(:,:) , k2_b(:,:) , k3_b(:,:)
139 : ! integer, allocatable :: k1_b2(:,:), k2_b2(:,:), k3_b2(:,:)
140 0 : integer, allocatable :: irreduc(:),mapkoper(:)
141 0 : integer, allocatable :: irreduc_q(:),mapqoper(:)
142 0 : integer, allocatable :: shiftkpt(:,:),pair_to_do(:,:)
143 0 : integer, allocatable :: shiftqpt(:,:),pair_to_do_q(:,:)
144 0 : integer, allocatable :: maptopair(:,:,:)
145 0 : integer, allocatable :: maptopair_q(:,:,:)
146 0 : integer, allocatable :: counts(:),displs(:)
147 0 : integer, allocatable :: gb(:,:,:),bpt(:,:)
148 0 : integer, allocatable :: gb_q(:,:,:),bpt_q(:,:)
149 0 : INTEGER, ALLOCATABLE :: innerEig_idList(:)
150 0 : real, allocatable :: we(:),we_b(:),we_b2(:)
151 0 : real, allocatable :: eigg(:)
152 0 : real, allocatable :: vr(:,:,:),vz(:,:,:)
153 0 : real, allocatable :: flo(:,:,:,:,:)
154 0 : real, allocatable :: ff(:,:,:,:,:),gg(:,:,:,:,:)
155 0 : real, allocatable :: kdiff(:,:),qdiff(:,:),zero_qdiff(:,:)
156 0 : complex, allocatable :: vpw(:,:),vzxy(:,:,:,:)
157 0 : complex, allocatable :: uHu(:,:,:,:,:)
158 : c complex, allocatable :: uHuold(:,:)
159 0 : complex, allocatable :: acof_b(:,:,:),acof_b2(:,:,:)
160 0 : complex, allocatable :: bcof_b(:,:,:),bcof_b2(:,:,:)
161 0 : complex, allocatable :: ccof_b(:,:,:,:),ccof_b2(:,:,:,:)
162 0 : complex, allocatable :: tdd(:,:,:,:,:),tdu(:,:,:,:,:)
163 0 : complex, allocatable :: tud(:,:,:,:,:),tuu(:,:,:,:,:)
164 0 : complex, allocatable :: tdulo(:,:,:,:,:,:),tuulo(:,:,:,:,:,:)
165 0 : complex, allocatable :: tulod(:,:,:,:,:,:),tulou(:,:,:,:,:,:)
166 0 : complex, allocatable :: tuloulo(:,:,:,:,:,:,:)
167 0 : complex, allocatable :: tdd_soc(:,:,:,:),tdu_soc(:,:,:,:)
168 0 : complex, allocatable :: tud_soc(:,:,:,:),tuu_soc(:,:,:,:)
169 0 : complex, allocatable :: tdulo_soc(:,:,:,:,:)
170 0 : complex, allocatable :: tuulo_soc(:,:,:,:,:)
171 0 : complex, allocatable :: tulod_soc(:,:,:,:,:)
172 0 : complex, allocatable :: tulou_soc(:,:,:,:,:)
173 0 : complex, allocatable :: tuloulo_soc(:,:,:,:,:,:)
174 :
175 : c ..local arrays..
176 : character(len=2) :: spin012(0:2)
177 : data spin012/' ', '.1', '.2'/
178 : character(len=3) :: spin12(2)
179 : data spin12/'WF1' , 'WF2'/
180 : character(len=8) :: name(10)
181 0 : integer :: n_bands(0:neigd),ngopr1(natd)
182 : real :: bkpt(3),bkpt_b(3),bkpt_b2(3),bkrot(3)
183 0 : real :: eig(neigd),eig_b(neigd),eig_b2(neigd)
184 0 : real :: uuilon(nlod,ntypd),duilon(nlod,ntypd)
185 0 : real :: ulouilopn(nlod,nlod,ntypd)
186 0 : real :: evac(2,max(2,jspd))
187 : real :: evdu(2,max(jspd,2))
188 : real :: qpt_i(3),qptb_i(3)
189 0 : real :: alph_i(ntypd),alphb_i(ntypd)
190 0 : real :: beta_i(ntypd),betab_i(ntypd)
191 : real :: cp_time(9)
192 :
193 : c ..local scalars..
194 : character(len=6) :: filename
195 : character(len=8) :: dop,iop
196 : character(len=12) fending
197 : character(len=30) fstart
198 : logical :: l_p0,l_bkpts,l_proj,l_file
199 : logical :: l_bqpts,l_gwf,l_exist,l_nocosoc,l_symcheck
200 : logical :: l_skip_sph,l_skip_non,l_skip_soc
201 : logical :: l_skip_int,l_skip_vac,l_skip_loc
202 : integer :: lmd,nlotot,n,iter,ikpt,ikpt_b,ikpt_b2,iqpt,iqpt_b
203 : integer :: addnoco,addnoco2,funbas,loplod,igvm2
204 : integer :: nn,nkpts,i,j,l,i_rec,m,nwf,nwfp
205 : integer :: jsp_start,jsp_end,nrec,nrec_b,nrec1
206 : integer :: nodeu,noded,n_size,na,n_rank,nbnd,numbands
207 : integer :: i1,i2,i3,in,lda
208 : integer :: nmat,nmat_b,nmat_b2,nmat_qb
209 : integer :: nbands,nbands_b,nbands_b2,nbands_qb
210 : integer :: nslibd,nslibd_b,nslibd_b2,nslibd_qb
211 : integer :: noccbd,noccbd_b,noccbd_b2,noccbd_qb
212 : integer :: kptibz,kptibz_b,kptibz_b2
213 : integer :: qptibz, qptibz_b
214 : integer :: oper,oper_b,oper_b2,oper_q, oper_qb
215 : integer :: nwfs,nntot,nntot_q,fullnkpts,fullnqpts
216 : integer :: kpt,qpt,j1,j2,j3,k,ikpt_help,iqpt_help
217 : integer :: wannierspin,jspin,jspin_b,jspin2
218 : integer :: jspin3,jspin4_b,jspin4,jspin5,tspin,tspin2
219 : integer :: n_start,n_end,mlotot,mlolotot,err
220 : integer :: mlot_d,mlolot_d,ilo,dir,length
221 : integer :: npotmatfile,ig3,maxvac,irec,imz,ivac,ipot
222 : integer :: funit_start,band_help,sign2
223 : integer :: doublespin,doublespin_max,nrec5
224 : integer :: aoff,d1,d10,d100
225 : integer :: eig_id
226 : real :: tpi,wronk,wk,wk_b,wk_b2
227 : real :: t0,t00,t1,t_myTlmplm,t_init
228 : real :: t_int,t_sph,t_vac,t_abcof,t_eig,t_total
229 : real :: efermi,htr2ev
230 : real :: theta_i, thetab_i, phi_i, phib_i
231 : complex :: nsfactor,nsfactor_b,nsfactor_b2
232 : integer :: nbasfcn,nbasfcn_b,nbasfcn_b2
233 :
234 0 : TYPE(t_usdus) :: usdus
235 0 : TYPE(t_mat) :: zMat, zzMat, zMat_b, zMat_b2
236 0 : TYPE(t_lapw) :: lapw, lapw_b, lapw_b2
237 :
238 0 : eig_id = eig_idList(1)
239 :
240 :
241 : c ..initializations..
242 0 : call cpu_time(t00)
243 :
244 0 : t_init = 0.
245 0 : t_myTlmplm = 0.
246 0 : t_eig = 0.
247 0 : t_abcof = 0.
248 0 : t_int = 0.
249 0 : t_sph = 0.
250 0 : t_vac = 0.
251 0 : t_total = 0.
252 0 : htr2ev = 27.2
253 0 : nntot_q = 1
254 0 : fullnqpts = 1
255 0 : funit_start = 5000
256 :
257 : ! In some modes the code skips the Hamiltonian setup and
258 : ! jumps directly into the Wannier branch.
259 : if(.not.allocated (usdus%us))then
260 0 : call usdus%init(atoms,input%jspins)
261 : endif
262 :
263 0 : aoff = iachar('1')-1
264 0 : d1 = mod(irank,10)
265 0 : IF (irank < 100) THEN
266 0 : d10 = int( (irank + 0.5)/10 )
267 0 : fstart = 'eig'//achar(d10+aoff)//achar(d1+aoff)
268 : ELSE
269 0 : d10 = mod((irank-d1)/10,10)
270 0 : d100 = (irank-10*d10-d1)/100
271 0 : IF ( d100.GE.10 ) d100 = d100 + iachar('7')
272 : fstart =
273 0 : + 'eig'//achar(d100+aoff)//achar(d10+aoff)//achar(d1+aoff)
274 : ENDIF
275 :
276 :
277 0 : ngopr1(:)=1
278 :
279 0 : l_p0 = .false.
280 0 : if (irank.eq.0) l_p0 = .true.
281 0 : l_nocosoc = l_noco.or.l_soc
282 :
283 0 : tpi = 2* pimach()
284 0 : lmd = lmaxd*(lmaxd+2)
285 :
286 : !!! should be changed in case the windows are really used
287 0 : nkpts = nkpt
288 :
289 : ! do we have to construct GWF ?
290 : l_gwf = .false.
291 0 : l_gwf = l_sgwf.or.l_socgwf
292 :
293 :
294 : c-----read the input file to determine what to do
295 : ! call wann_read_inp(input,l_p0,wann) !Call wann_read_inp in fleur_init
296 :
297 0 : if(wann%l_byenergy.and.wann%l_byindex) CALL juDFT_error
298 0 : + ("byenergy.and.byindex",calledby ="wannier")
299 0 : if(wann%l_byenergy.and.wann%l_bynumber) CALL juDFT_error
300 0 : + ("byenergy.and.bynumber",calledby ="wannier")
301 0 : if(wann%l_bynumber.and.wann%l_byindex) CALL juDFT_error
302 0 : + ("bynumber.and.byindex",calledby ="wannier")
303 0 : if(.not.(wann%l_bynumber.or.wann%l_byindex.or.wann%l_byenergy))
304 0 : & CALL juDFT_error("no rule to sort bands",calledby ="wannier")
305 :
306 :
307 0 : efermi=ef
308 0 : if(.not.wann%l_fermi)efermi=0.0
309 :
310 : #ifdef CPP_MPI
311 0 : call MPI_BARRIER(MPI_COMM_WORLD,ierr(1))
312 : #endif
313 :
314 : c**************************************************************
315 : c for bzsym=.true.: determine mapping between kpts and w90kpts
316 : c**************************************************************
317 0 : if (wann%l_bzsym) then
318 0 : l_file=.false.
319 0 : inquire(file='w90kpts',exist=l_file)
320 0 : if(.not.l_file) CALL juDFT_error
321 0 : + ("w90kpts not found, needed if bzsym",calledby ="wannier")
322 0 : open(412,file='w90kpts',form='formatted')
323 0 : read(412,*)fullnkpts
324 0 : close(412)
325 0 : if(l_p0)print*,"fullnkpts=",fullnkpts
326 0 : if(fullnkpts<nkpts) CALL juDFT_error("fullnkpts.lt.nkpts"
327 0 : + ,calledby ="wannier")
328 0 : allocate(irreduc(fullnkpts),mapkoper(fullnkpts))
329 0 : allocate(shiftkpt(3,fullnkpts))
330 0 : l_file=.false.
331 0 : inquire(file='kptsmap',exist=l_file)
332 0 : if(.not.l_file) CALL juDFT_error
333 0 : + ("kptsmap not found, needed if bzsym",calledby ="wannier")
334 0 : open(713,file='kptsmap')
335 0 : do i=1,fullnkpts
336 0 : read(713,*)kpt,irreduc(i),mapkoper(i),shiftkpt(:,i)
337 0 : if(kpt/=i) CALL juDFT_error("kpt.ne.i",calledby ="wannier")
338 0 : if(l_p0)print*,i,irreduc(i),mapkoper(i)
339 : enddo
340 0 : close(713)
341 0 : if(maxval(irreduc(:))/=nkpts) CALL juDFT_error
342 0 : + ("max(irreduc(:))/=nkpts",calledby ="wannier")
343 : else
344 0 : fullnkpts=nkpts
345 : endif
346 :
347 :
348 0 : if(l_gwf) fullnqpts = nparampts
349 :
350 :
351 0 : nrec = 0
352 0 : if(l_p0)then
353 0 : write (*,*) 'fermi energy:',efermi
354 0 : write (*,*) 'emin,emax=',e1s,e2s
355 0 : write (*,*) 'nbasfcn =',nbasfcn
356 : endif
357 0 : nlotot = 0
358 0 : mlotot = 0
359 0 : mlolotot = 0
360 0 : do n = 1, ntype
361 0 : mlotot = mlotot + nlo(n)
362 0 : mlolotot = mlolotot + nlo(n)*(nlo(n)+1)/2
363 0 : do l = 1,nlo(n)
364 0 : nlotot = nlotot + neq(n) * ( 2*llo(l,n) + 1 )
365 : enddo
366 : enddo
367 :
368 :
369 0 : allocate(counts(0:isize-1),displs(0:isize-1))
370 0 : call array_split(fullnkpts,isize,counts,displs)
371 :
372 : c**********************************************************
373 : ccccccccccccccc read in the bkpts file ccccccccccccccccc
374 : c**********************************************************
375 0 : l_bkpts = .false.
376 0 : inquire (file='bkpts',exist=l_bkpts)
377 0 : if (.not.l_bkpts) CALL juDFT_error("need bkpts for matrixmmn"
378 0 : + ,calledby ="wannier")
379 0 : open (202,file='bkpts',form='formatted',status='old')
380 0 : rewind (202)
381 0 : read (202,'(i4)') nntot
382 0 : if(l_p0)then
383 0 : write (*,*) 'nntot=',nntot
384 0 : write(*,*) 'fullnkpts=',fullnkpts
385 0 : write(*,*) 'nkpts=',nkpts
386 : endif
387 0 : allocate ( gb(1:3,1:nntot,1:fullnkpts),bpt(1:nntot,1:fullnkpts))
388 0 : do ikpt=1,fullnkpts
389 0 : do nn=1,nntot
390 : read (202,'(2i6,3x,3i4)')
391 0 : & ikpt_help,bpt(nn,ikpt),(gb(i,nn,ikpt),i=1,3)
392 0 : if (ikpt/=ikpt_help) CALL juDFT_error("ikpt.ne.ikpt_help"
393 0 : + ,calledby ="wannier")
394 0 : if (bpt(nn,ikpt)>fullnkpts) CALL juDFT_error("bpt.gt.fullnkpts"
395 0 : + ,calledby ="wannier")
396 : enddo
397 : enddo
398 0 : close (202)
399 0 : allocate(kdiff(3,nntot))
400 :
401 : c**********************************************************
402 : ccccccccccccccc read in the bqpts file ccccccccccccccccc
403 : c**********************************************************
404 0 : if (l_gwf.or.l_ms) then ! for Omega functional minimization
405 0 : l_bqpts = .false.
406 0 : inquire (file='bqpts',exist=l_bqpts)
407 0 : if (.not.l_bqpts) CALL juDFT_error("need bqpts for matrixmmn"
408 0 : + ,calledby ="wannier")
409 0 : open (202,file='bqpts',form='formatted',status='old')
410 0 : rewind (202)
411 0 : read (202,'(i4)') nntot_q
412 0 : if(l_p0)then
413 0 : write (*,*) 'nntot_q=',nntot_q
414 0 : write(*,*) 'fullnqpts=',fullnqpts
415 : endif
416 : allocate ( gb_q(1:3,1:nntot_q,1:fullnqpts),
417 0 : & bpt_q(1:nntot_q,1:fullnqpts))
418 0 : do iqpt=1,fullnqpts
419 0 : do nn=1,nntot_q
420 : read (202,'(2i6,3x,3i4)')
421 0 : & iqpt_help,bpt_q(nn,iqpt),(gb_q(i,nn,iqpt),i=1,3)
422 0 : if (iqpt/=iqpt_help) CALL juDFT_error("iqpt.ne.iqpt_help"
423 0 : + ,calledby ="wannier")
424 0 : if (bpt_q(nn,iqpt)>fullnqpts)
425 0 : & CALL juDFT_error("bpt_q.gt.fullnqpts",calledby ="wannier")
426 : enddo
427 : enddo
428 0 : close (202)
429 0 : allocate(qdiff(3,nntot_q))
430 0 : allocate(zero_qdiff(3,nntot_q))
431 0 : zero_qdiff=0.0
432 : endif
433 :
434 :
435 : ! when treating gen. WF for spin spirals, the Brillouin zone
436 : ! of q-points is twice as large compared to k-BZ. Thus,
437 : ! the G-vectors connecting neighbors across the boundary
438 : ! need to be doubled
439 0 : if(l_sgwf) gb_q = 2*gb_q
440 0 : if(l_socgwf) gb_q = 2*gb_q
441 :
442 0 : if(wann%l_finishgwf) goto 9110
443 : c********************************************************
444 : c find symmetry-related elements in mmkb
445 : c********************************************************
446 0 : allocate(maptopair(3,fullnkpts,nntot))
447 0 : allocate(pair_to_do(fullnkpts,nntot))
448 : call wann_mmnk_symm(input,kpts,
449 : > fullnkpts,nntot,bpt,gb,wann%l_bzsym,
450 : > irreduc,mapkoper,l_p0,film,nop,invtab,mrot,
451 : > tau,
452 0 : < pair_to_do,maptopair,kdiff,.false.,param_file)
453 :
454 : ! do the same for q-points to construct GWFs
455 0 : if(l_gwf)then
456 0 : allocate(maptopair_q(3,fullnqpts,nntot_q))
457 0 : allocate(pair_to_do_q(fullnqpts,nntot_q))
458 : call wann_mmnk_symm(input,kpts,
459 : > fullnqpts,nntot_q,bpt_q,gb_q,wann%l_bzsym,
460 : > irreduc_q,mapqoper,l_p0,.false.,1,invtab(1),mrot(:,:,1),
461 : > tau,
462 0 : < pair_to_do_q,maptopair_q,qdiff,.true.,param_file)
463 : endif
464 :
465 :
466 : c*********************************************************
467 : cccccccccccccccc initialize the potential cccccccccccc
468 : c*********************************************************
469 :
470 0 : if(.not. l_noco) then
471 0 : allocate ( vpw(n3d,jspd),vzxy(nmzxyd,stars%ng2-1,2,jspd) )
472 : else
473 0 : allocate ( vpw(n3d,4),vzxy(nmzxyd,stars%ng2-1,2,4) )
474 : endif
475 :
476 0 : allocate (vz(nmzd,2,4))
477 0 : allocate (vr(jmtd,ntypd,jspd))
478 :
479 0 : vpw(:,1:SIZE(vTot%pw,2)) = vTot%pw(:,1:SIZE(vTot%pw,2))
480 0 : IF(film) THEN
481 : vz(:,:,1:SIZE(vTot%vac,4)) =
482 0 : & REAL(vTot%vac(:,1,:,1:SIZE(vTot%vac,4)))
483 0 : IF (SIZE(vTot%vac,4)==3) vz(:,:,4) = AIMAG(vTot%vac(:,1,:,3))
484 : vzxy(:,:,:,1:SIZE(vTot%vac,4)) =
485 0 : & vTot%vac(:nmzxyd,2:,:,1:SIZE(vTot%vac,4))
486 : END IF
487 :
488 0 : do jspin = 1,jspins
489 0 : do n = 1, ntype
490 0 : do j = 1,jri(n)
491 0 : vr(j,n,jspin) = vTot%mt(j,0,n,jspin)
492 : enddo
493 : enddo
494 : enddo
495 :
496 0 : if(.not.film) deallocate(vz,vzxy)
497 :
498 0 : if(l_noco)then
499 0 : npotmatfile=25
500 :
501 : OPEN (npotmatfile,FILE='potmat',FORM='unformatted',
502 0 : + STATUS='old')
503 : c---> load the interstitial potential
504 0 : if(l_p0 .and. nq3.ne.n3d) then
505 0 : write(*,*)'WARNING for reading potmat: nq3.ne.n3d'
506 : endif
507 0 : vpw = cmplx(0.0, 0.0)
508 0 : READ (npotmatfile) (vpw(ig3,1),ig3=1,n3d)
509 0 : READ (npotmatfile) (vpw(ig3,2),ig3=1,n3d)
510 0 : READ (npotmatfile) (vpw(ig3,3),ig3=1,n3d)
511 0 : vpw(:,4) = conjg(vpw(:,3))
512 0 : if(film) then
513 : maxvac=2
514 0 : DO ivac = 1,maxvac
515 : c---> if the two vacuua are equivalent, the potential file has to
516 : c---> be backspaced, because the potential is the same at both
517 : c---> surfaces of the film
518 0 : IF ((ivac.EQ.2) .AND. (nvac.EQ.1)) THEN
519 0 : DO irec = 1,4
520 0 : BACKSPACE (npotmatfile)
521 : ENDDO
522 : ENDIF
523 : c---> load the non-warping part of the potential
524 0 : READ (npotmatfile)((vz(imz,ivac,ipot),imz=1,nmzd),ipot=1,4)
525 :
526 : c---> load the warping part of the potential
527 0 : DO ipot = 1,3
528 0 : READ (npotmatfile)((vzxy(imz,igvm2,ivac,ipot),
529 0 : + imz=1,nmzxy),igvm2=1,nq2-1)
530 : ENDDO
531 0 : vzxy(:,:,:,4) = conjg(vzxy(:,:,:,3))
532 : enddo
533 : endif
534 0 : CLOSE (npotmatfile)
535 : endif
536 :
537 :
538 0 : if(film .and. l_p0) write(*,*)'nvac',nvac
539 :
540 : cccccccccccccccc end of the potential part ccccccccccc
541 0 : wannierspin=jspd
542 0 : if(l_soc) wannierspin=2
543 :
544 0 : allocate(flo(ntypd,jmtd,2,nlod,2))
545 0 : allocate ( ff(ntypd,jmtd,2,0:lmaxd,2) )
546 0 : allocate ( gg(ntypd,jmtd,2,0:lmaxd,2) )
547 : ! allocate ( kveclo(nlotot),nv(wannierspin) )
548 : ! allocate ( kveclo_b(nlotot),nv_b(wannierspin) )
549 : ! allocate ( kveclo_b2(nlotot),nv_b2(wannierspin) )
550 : ! allocate ( k1(nvd,wannierspin),k2(nvd,wannierspin),
551 : ! & k3(nvd,wannierspin) )
552 : ! allocate ( k1_b(nvd,wannierspin),k2_b(nvd,wannierspin),
553 : ! & k3_b(nvd,wannierspin) )
554 : ! allocate ( k1_b2(nvd,wannierspin),k2_b2(nvd,wannierspin),
555 : ! & k3_b2(nvd,wannierspin) )
556 :
557 0 : if(l_nocosoc) then
558 0 : doublespin_max=4
559 : else
560 0 : doublespin_max=wannierspin
561 : endif
562 :
563 : c if(l_soc.and.(jspins.eq.1)) doublespin_max=wannierspin
564 :
565 0 : l_skip_int = .false.; l_skip_soc = .false.; l_skip_vac = .false.
566 0 : l_skip_sph = .false.; l_skip_non = .false.; l_skip_loc = .false.
567 0 : inquire(file='debug_uHu',exist=l_exist)
568 0 : if(l_exist) then
569 0 : open(888,file='debug_uHu')
570 0 : read(888,*)l_skip_int
571 0 : read(888,*)l_skip_sph
572 0 : read(888,*)l_skip_non
573 0 : read(888,*)l_skip_soc
574 0 : read(888,*)l_skip_loc
575 0 : read(888,*)l_skip_vac
576 0 : read(888,*)doublespin_max
577 0 : close(888)
578 0 : if(l_p0) then
579 0 : write(*,*)'skip INT :',l_skip_int
580 0 : write(*,*)'skip SPH :',l_skip_sph
581 0 : write(*,*)'skip NON :',l_skip_non
582 0 : write(*,*)'skip SOC :',l_skip_soc
583 0 : write(*,*)'skip LOC :',l_skip_loc
584 0 : write(*,*)'skip VAC :',l_skip_vac
585 0 : write(*,*)'doublespin_max:',doublespin_max
586 : endif
587 : endif
588 :
589 0 : tspin = doublespin_max
590 0 : if(fullnqpts.eq.1) tspin=1
591 0 : if(l_p0) write(*,*)'tspin',tspin
592 0 : allocate( tdd(0:lmd,0:lmd,ntypd,nntot*nntot,tspin) )
593 0 : allocate( tdu(0:lmd,0:lmd,ntypd,nntot*nntot,tspin) )
594 0 : allocate( tud(0:lmd,0:lmd,ntypd,nntot*nntot,tspin) )
595 0 : allocate( tuu(0:lmd,0:lmd,ntypd,nntot*nntot,tspin) )
596 0 : allocate( tdulo(0:lmd,nlod,-llod:llod,ntypd,nntot*nntot,tspin) )
597 0 : allocate( tuulo(0:lmd,nlod,-llod:llod,ntypd,nntot*nntot,tspin) )
598 0 : allocate( tulou(0:lmd,nlod,-llod:llod,ntypd,nntot*nntot,tspin) )
599 0 : allocate( tulod(0:lmd,nlod,-llod:llod,ntypd,nntot*nntot,tspin) )
600 : allocate( tuloulo(nlod,-llod:llod,nlod,-llod:llod,
601 0 : > ntypd,nntot*nntot,tspin) )
602 0 : allocate( tdd_soc(0:lmd,0:lmd,ntypd,nntot*nntot) )
603 0 : allocate( tdu_soc(0:lmd,0:lmd,ntypd,nntot*nntot) )
604 0 : allocate( tud_soc(0:lmd,0:lmd,ntypd,nntot*nntot) )
605 0 : allocate( tuu_soc(0:lmd,0:lmd,ntypd,nntot*nntot) )
606 0 : allocate( tdulo_soc(0:lmd,nlod,-llod:llod,ntypd,nntot*nntot) )
607 0 : allocate( tuulo_soc(0:lmd,nlod,-llod:llod,ntypd,nntot*nntot) )
608 0 : allocate( tulou_soc(0:lmd,nlod,-llod:llod,ntypd,nntot*nntot) )
609 0 : allocate( tulod_soc(0:lmd,nlod,-llod:llod,ntypd,nntot*nntot) )
610 0 : allocate( tuloulo_soc(nlod,-llod:llod,nlod,-llod:llod,
611 0 : > ntypd,nntot*nntot) )
612 :
613 0 : tuu = cmplx(0.,0.); tdu = cmplx(0.,0.)
614 0 : tud = cmplx(0.,0.); tdd = cmplx(0.,0.)
615 0 : tuulo = cmplx(0.,0.); tdulo = cmplx(0.,0.)
616 0 : tulou = cmplx(0.,0.); tulod = cmplx(0.,0.)
617 0 : tuloulo = cmplx(0.,0.)
618 0 : tuu_soc = cmplx(0.,0.); tdu_soc = cmplx(0.,0.)
619 0 : tud_soc = cmplx(0.,0.); tdd_soc = cmplx(0.,0.)
620 0 : tuulo_soc = cmplx(0.,0.); tdulo_soc = cmplx(0.,0.)
621 0 : tulou_soc = cmplx(0.,0.); tulod_soc = cmplx(0.,0.)
622 0 : tuloulo_soc = cmplx(0.,0.)
623 :
624 0 : call cpu_time(t1)
625 0 : t_init = t1-t00
626 :
627 : c*****************************************************************c
628 : c START Q LOOP c
629 : c*****************************************************************c
630 0 : do 314 iqpt = 1,fullnqpts ! loop by q-points starts
631 :
632 0 : ALLOCATE(innerEig_idList(nntot_q))
633 :
634 0 : qptibz=iqpt
635 0 : if(wann%l_bzsym .AND. l_gwf) qptibz=irreduc_q(iqpt)
636 : if(wann%l_bzsym .AND. l_gwf) oper_q=mapqoper(iqpt)
637 :
638 0 : qpt_i = qss
639 0 : alph_i = alph
640 0 : beta_i = beta
641 0 : theta_i = theta
642 0 : phi_i = phi
643 0 : if(l_sgwf.or.l_ms) then
644 0 : qpt_i(:) = param_vec(:,qptibz)
645 0 : alph_i(:) = param_alpha(:,qptibz)
646 0 : elseif(l_socgwf) then
647 0 : if(l_dim(2)) phi_i = tpi*param_vec(2,qptibz)
648 0 : if(l_dim(3)) theta_i = tpi*param_vec(3,qptibz)
649 : endif
650 :
651 0 : if(l_p0) then
652 0 : write(*,*)'qpt_i',qpt_i
653 0 : do n=1,ntype
654 0 : write(*,*)'n=',n,alph_i(n),beta_i(n)
655 : enddo
656 0 : write(*,*)'theta_i=',theta_i
657 0 : write(*,*)'phi_i=',phi_i
658 : endif
659 :
660 0 : IF (l_gwf) THEN
661 0 : do iqpt_b=1,nntot_q
662 :
663 0 : innerEig_idList(iqpt_b) = eig_idList(bpt_q(iqpt_b,iqpt))
664 :
665 : ! WRITE(fending,'("_",i4.4)')bpt_q(iqpt_b,iqpt)
666 : ! innerEig_idList(iqpt_b)=open_eig(fmpi%mpi_comm,
667 : ! + lapw%dim_nbasfcn(),input%neig,
668 : ! + nkpts,wannierspin,atoms%lmaxd,
669 : ! + atoms%nlod,atoms%ntype,atoms%nlotot,
670 : ! + l_noco,.FALSE.,l_real,l_soc,.FALSE.,
671 : ! + mpi%n_size,filename=trim(fstart)//fending,
672 : ! + layers=banddos%layers,nstars=banddos%nstars,
673 : ! + ncored=29,nsld=atoms%nat,
674 : ! + nat=atoms%nat,l_dos=banddos%dos.OR.input%cdinf,
675 : ! + l_mcd=banddos%l_mcd,l_orb=banddos%l_orb)
676 :
677 : enddo
678 :
679 0 : eig_id = eig_idList(qptibz)
680 :
681 : ! WRITE(fending,'("_",i4.4)')qptibz
682 : ! eig_id=open_eig(fmpi%mpi_comm,lapw%dim_nbasfcn(),input%neig,
683 : ! + nkpts,wannierspin,atoms%lmaxd,
684 : ! + atoms%nlod,atoms%ntype,atoms%nlotot,
685 : ! + l_noco,.FALSE.,l_real,l_soc,.FALSE.,
686 : ! + mpi%n_size,filename=trim(fstart)//fending,
687 : ! + layers=banddos%layers,nstars=banddos%nstars,
688 : ! + ncored=29,nsld=atoms%nat,
689 : ! + nat=atoms%nat,l_dos=banddos%dos.OR.input%cdinf,
690 : ! + l_mcd=banddos%l_mcd,l_orb=banddos%l_orb)
691 :
692 : ! OPEN(66,file=trim(fstart)//fending,access='direct',
693 : ! > form='unformatted',recl=irecl,status='old')
694 0 : ELSEIF(l_ms) THEN
695 :
696 0 : eig_id = eig_idList(qptibz)
697 :
698 : ! WRITE(fending,'("_",i4.4)')qptibz
699 : ! eig_id=open_eig(fmpi%mpi_comm,lapw%dim_nbasfcn(),input%neig,
700 : ! + nkpts,wannierspin,atoms%lmaxd,
701 : ! + atoms%nlod,atoms%ntype,atoms%nlotot,
702 : ! + l_noco,.FALSE.,l_real,l_soc,.FALSE.,
703 : ! + mpi%n_size,filename=trim(fstart)//fending,
704 : ! + layers=banddos%layers,nstars=banddos%nstars,
705 : ! + ncored=29,nsld=atoms%nat,
706 : ! + nat=atoms%nat,l_dos=banddos%dos.OR.input%cdinf,
707 : ! + l_mcd=banddos%l_mcd,l_orb=banddos%l_orb)
708 :
709 : ! OPEN(66,file=trim(fstart)//fending,access='direct',
710 : ! > form='unformatted',recl=irecl,status='old')
711 : ELSE
712 0 : fending=''
713 : ENDIF ! l_gwf.or.l_ms
714 0 : nrec=0
715 0 : nrec_b=0
716 :
717 : c****************************************************
718 : c cycle by spins starts!
719 : c****************************************************
720 0 : do 110 doublespin=1,doublespin_max ! cycle by spins
721 0 : if(l_p0) write(*,*)'spin loop:',doublespin
722 :
723 : c jspin=mod(doublespin+1,2)+1
724 : c jspin_b=jspin
725 : c if(doublespin.eq.3) jspin_b=2
726 : c if(doublespin.eq.4) jspin_b=1
727 0 : jspin_b=mod(doublespin+1,2)+1
728 0 : jspin=jspin_b
729 0 : if(doublespin.eq.3) jspin=2
730 0 : if(doublespin.eq.4) jspin=1
731 :
732 0 : tspin2 = doublespin
733 0 : if(fullnqpts.eq.1) tspin2 = 1
734 :
735 0 : nrec_b = nrec
736 :
737 0 : if(.not.l_noco) then
738 0 : nrec = (jspin-1)*nkpts
739 0 : nrec_b = (jspin_b-1)*nkpts
740 : endif
741 :
742 : c...read number of bands and wannier functions from file proj
743 :
744 : c..reading the proj.1 / proj.2 / proj file
745 0 : l_proj=.false.
746 0 : do j=jspin,0,-1
747 0 : inquire(file=trim('proj'//spin012(j)),exist=l_proj)
748 0 : if(l_proj)then
749 0 : filename='proj'//spin012(j)
750 0 : exit
751 : endif
752 : enddo
753 :
754 0 : if(l_proj)then
755 0 : open (203,file=trim(filename),status='old')
756 0 : rewind (203)
757 0 : read (203,*) nwfs,numbands
758 0 : rewind (203)
759 0 : close (203)
760 : elseif(wann%l_projmethod.or.wann%l_bestproj
761 0 : & .or.wann%l_matrixamn)then
762 0 : CALL juDFT_error("no proj/proj.1/proj.2",calledby ="wannier")
763 : endif
764 :
765 :
766 0 : jspin2=jspin
767 0 : if(l_soc .and. jspins.eq.1)jspin2=1
768 0 : jsp_start = jspin ; jsp_end = jspin
769 :
770 : cccccccccccc read in the eigenvalues and vectors cccccc
771 0 : do jspin5=1,wannierspin!1,2
772 0 : jsp_start=jspin5; jsp_end=jspin5
773 0 : nrec5=0
774 : if(.not.l_noco) nrec5 = (jspin5-1)*nkpts
775 : ! CALL judft_error("TODO:wann_uHu")
776 : ! call cdn_read0(eig_id,irank,isize,jspin5,wannierspin,l_noco,
777 : ! < ello,evac,epar,wk,n_bands,n_size)
778 :
779 : CALL cdn_read0(eig_id,fmpi%irank,fmpi%isize,jspin5,
780 : & input%jspins, !wannierspin instead of DIMENSION%jspd?&
781 0 : & noco%l_noco, n_bands,n_size)
782 :
783 :
784 : enddo
785 :
786 : c.. now we want to define the maximum number of the bands by all kpts
787 0 : nbnd = 0
788 0 : i_rec = 0 ; n_rank = 0
789 :
790 0 : if(l_p0)then
791 : call wann_maxbnd(
792 : > eig_id,l_real,
793 : > lmaxd,ntypd,nlod,neigd,nvd,wannierspin,
794 : > isize,jspin,lapw%dim_nbasfcn(),nlotot,
795 : > l_ss,l_noco,nrec,fullnkpts,
796 : > wann%l_bzsym,wann%l_byindex,wann%l_bynumber,
797 : > wann%l_byenergy,
798 : > irreduc,wann%band_min(jspin),
799 : > wann%band_max(jspin),
800 : > numbands,e1s,e2s,efermi,nkpts,
801 0 : < nbnd,l_gwf,iqpt)
802 : endif!l_p0
803 :
804 : ! nbnd is calculated for process zero and is sent here to the others
805 : #ifdef CPP_MPI
806 0 : if(l_p0)then
807 0 : do cpu_index=1,isize-1
808 : call MPI_SEND(nbnd,1,MPI_INTEGER,cpu_index,1,MPI_COMM_WORLD,
809 0 : & ierr(1))
810 : enddo
811 : else
812 0 : call MPI_RECV(nbnd,1,MPI_INTEGER,0,1,MPI_COMM_WORLD,stt,ierr(1))
813 : endif
814 : #endif
815 :
816 : c##################################################################
817 0 : if(.not.allocated(uHu)) then
818 0 : allocate(uHu(nbnd,nbnd,nntot,nntot,counts(irank)))
819 : c allocate(uHuold(nbnd,nbnd))
820 0 : uHu = cmplx(0.,0.)
821 : c uHuold = cmplx(0.,0.)
822 : endif
823 :
824 :
825 :
826 : ! set up T-matrix for spherical and non-spherical part
827 : ! only once since it is k- and q-independent !
828 0 : if(iqpt.eq.1) then
829 0 : if(l_p0) write(*,*)'tspin2=',tspin2
830 0 : tuu(:,:,:,:,tspin2) = cmplx(0.,0.)
831 0 : tdu(:,:,:,:,tspin2) = cmplx(0.,0.)
832 0 : tud(:,:,:,:,tspin2) = cmplx(0.,0.)
833 0 : tdd(:,:,:,:,tspin2) = cmplx(0.,0.)
834 0 : tuulo(:,:,:,:,:,tspin2) = cmplx(0.,0.)
835 0 : tdulo(:,:,:,:,:,tspin2) = cmplx(0.,0.)
836 0 : tulou(:,:,:,:,:,tspin2) = cmplx(0.,0.)
837 0 : tulod(:,:,:,:,:,tspin2) = cmplx(0.,0.)
838 0 : tuloulo(:,:,:,:,:,:,tspin2) = cmplx(0.,0.)
839 :
840 0 : do jspin4=1,2
841 0 : jspin3=jspin4
842 0 : if(jspins.eq.1) jspin3=1
843 0 : na = 1
844 0 : do 40 n = 1,ntype
845 0 : do 30 l = 0,lmax(n)
846 : c...compute the l-dependent, k-independent radial MT- basis functions
847 :
848 : call radfun(
849 : > l,n,jspin4,enpara%el0(l,n,jspin3),vr(1,n,jspin3),atoms,
850 : < ff(n,:,:,l,jspin4),gg(n,:,:,l,jspin4),usdus,
851 0 : < nodeu,noded,wronk)
852 :
853 0 : 30 continue
854 : c...and the local orbital radial functions
855 0 : do ilo = 1, nlo(n)
856 :
857 : call radflo(
858 : > atoms,n,jspin4,enpara%ello0(:,:,jspin3),vr(1,n,jspin3),
859 : > ff(n,1:,1:,0:,jspin4),gg(n,1:,1:,0:,jspin4),fmpi,
860 0 : < usdus,uuilon,duilon,ulouilopn,flo(n,:,:,:,jspin4))
861 :
862 : enddo
863 : c na = na + neq(n)
864 0 : 40 continue
865 : enddo!jspin3
866 :
867 :
868 0 : if(jspin.eq.jspin_b) then
869 0 : if(l_p0) write(*,*)'wann_uHu_tlmplm'
870 0 : jspin3=jspin
871 0 : if(jspins.eq.1) jspin3=1
872 0 : call cpu_time(t0)
873 : call wann_uHu_tlmplm(
874 : > memd,nlhd,ntypsd,ntypd,jmtd,lmaxd,jspd,ntype,dx,
875 : > rmsh,jri,lmax,ntypsy,natd,lnonsph,lmd,lmplmd,clnu,
876 : > mlh,nmem,llh,nlh,neq,irank,mlotot,mlolotot,
877 : > vTot%mt(:,:,:,jspin3),nlod,llod,loplod,
878 : > enpara%ello0(1,1,jspin3),
879 : > llo,nlo,lo1l,l_dulo,ulo_der,ff(:,:,:,:,jspin),
880 : > gg(:,:,:,:,jspin),flo(:,:,:,:,jspin),
881 : > ff(:,:,:,:,jspin_b),gg(:,:,:,:,jspin_b),
882 : > flo(:,:,:,:,jspin_b),kdiff,kdiff,nntot,nntot,bmat,bbmat,
883 : > vr(1,1,jspin3),enpara%el0(0,1,jspin3),invsat,
884 : > l_skip_sph,l_skip_non,l_skip_loc,
885 : < tuu(:,:,:,:,tspin2),tud(:,:,:,:,tspin2),
886 : > tdu(:,:,:,:,tspin2),tdd(:,:,:,:,tspin2),
887 : > tuulo(:,:,:,:,:,tspin2),tulou(:,:,:,:,:,tspin2),
888 : > tdulo(:,:,:,:,:,tspin2),tulod(:,:,:,:,:,tspin2),
889 0 : > tuloulo(:,:,:,:,:,:,tspin2))
890 0 : call cpu_time(t1)
891 0 : t_myTlmplm = t_myTlmplm + t1-t0
892 : endif
893 :
894 : endif!iqpt.eq.1
895 :
896 :
897 : ! compute SOC-contribution for each theta_i !
898 0 : if(l_soc.and. (.not.l_skip_soc)) then
899 0 : if(l_p0) write(*,*)'wann_uHu_soc'
900 0 : tuu_soc = cmplx(0.,0.); tdu_soc = cmplx(0.,0.)
901 0 : tud_soc = cmplx(0.,0.); tdd_soc = cmplx(0.,0.)
902 0 : tuulo_soc = cmplx(0.,0.); tdulo_soc = cmplx(0.,0.)
903 0 : tulou_soc = cmplx(0.,0.); tulod_soc = cmplx(0.,0.)
904 0 : tuloulo_soc = cmplx(0.,0.)
905 0 : call cpu_time(t0)
906 : call wann_uHu_soc(
907 : > input,atoms,
908 : > ntypd,jmtd,lmaxd,jspd,
909 : > ntype,dx,rmsh,jri,lmax,natd,
910 : > lmd,lmplmd,neq,irank,
911 : > nlod,llod,loplod,enpara%ello0,llo,nlo,lo1l,l_dulo,ulo_der,
912 : > ff(:,:,:,:,jspin),gg(:,:,:,:,jspin),flo(:,:,:,:,jspin),
913 : > ff(:,:,:,:,jspin_b),gg(:,:,:,:,jspin_b),flo(:,:,:,:,jspin_b),
914 : > kdiff,kdiff,nntot,nntot,bmat,bbmat,
915 : > vr,enpara%el0,jspin,jspin_b,jspins,
916 : > .true.,theta_i,phi_i,alph_i,beta_i,l_noco,l_skip_loc,
917 : < tuu_soc,tud_soc,
918 : > tdu_soc,tdd_soc,
919 : > tuulo_soc,tulou_soc,
920 : > tdulo_soc,tulod_soc,
921 0 : > tuloulo_soc)
922 0 : call cpu_time(t1)
923 0 : t_myTlmplm = t_myTlmplm + t1-t0
924 : endif!l_soc
925 :
926 :
927 0 : i_rec = 0 ; n_rank = 0
928 :
929 : c****************************************************************
930 : c.. loop by kpoints starts! each may be a separate task
931 : c****************************************************************
932 0 : if(l_p0) write(*,*)'start k-loop'
933 0 : do 10 ikpt = wann%ikptstart,fullnkpts ! loop by k-points starts
934 0 : kptibz=ikpt
935 0 : if(wann%l_bzsym) kptibz=irreduc(ikpt)
936 : if(wann%l_bzsym) oper=mapkoper(ikpt)
937 :
938 0 : if(kpt_on_node(ikpt,isize,counts,displs).eq.irank) then
939 : c write(*,*)irank,ikpt
940 0 : i_rec = i_rec + 1
941 : c if (mod(i_rec-1,isize).eq.irank) then
942 :
943 0 : allocate ( eigg(neigd) )
944 :
945 0 : n_start=1
946 0 : n_end=input%neig
947 :
948 0 : call cpu_time(t0)
949 : ! get current bkpt vector
950 :
951 : CALL lapw%init(input,noco,nococonv,kpts,atoms,sym,kptibz,cell,
952 0 : & fmpi)
953 :
954 : nbasfcn = MERGE(lapw%nv(1)+lapw%nv(2)+2*atoms%nlotot,
955 0 : & lapw%nv(1)+atoms%nlotot,noco%l_noco)
956 0 : CALL zzMat%init(l_real,nbasfcn,input%neig)
957 0 : CALL zMat%init(l_real,nbasfcn,input%neig)
958 :
959 :
960 : CALL cdn_read(
961 : & eig_id,
962 : & lapw%dim_nvd(),input%jspins,fmpi%irank,fmpi%isize, !wannierspin instead of DIMENSION%jspd?&
963 : & kptibz,jspin,lapw%dim_nbasfcn(),
964 : & noco%l_ss,noco%l_noco,input%neig,
965 : & n_start,n_end,
966 0 : & nbands,eigg,zzMat)
967 :
968 :
969 : c$$$ nslibd=0
970 : ! do i = 1,nbands
971 : ! if((eigg(i).ge.e1s.and.nslibd.lt.numbands
972 : ! & .and.wann%l_bynumber)
973 : ! &.or.(eigg(i).ge.e1s.and.eigg(i).le.e2s.and.wann%l_byenergy)
974 : ! &.or.(i.ge.wann%band_min(jspin)
975 : ! & .and.
976 : ! & (i.le.wann%band_max(jspin))
977 : ! & .and.
978 : ! & wann%l_byindex))then
979 : ! nslibd = nslibd + 1
980 : ! eig(nslibd) = eigg(i)
981 : ! we(nslibd) = we(i)
982 : ! if(l_noco)then
983 : ! ! funbas = nv(1) + nlotot
984 : ! ! funbas = funbas+nv(2) + nlotot
985 : ! else
986 : ! ! funbas = nv(jspin) + nlotot
987 : ! endif
988 : ! IF (zzMat%l_real) THEN
989 : ! ! do j = 1,funbas
990 : ! zMat%data_r(:,nslibd) = zzMat%data_r(:,i)
991 : ! ! enddo
992 : ! ELSE
993 : ! ! do j = 1,funbas
994 : ! zMat%data_c(:,nslibd) = zzMat%data_c(:,i)
995 : ! ! enddo
996 : ! END IF
997 : ! endif
998 : ! enddo
999 :
1000 :
1001 :
1002 :
1003 0 : call cpu_time(t1)
1004 0 : t_eig = t_eig + t1 - t0
1005 :
1006 :
1007 0 : allocate (we_b(neigd), we_b2(neigd))
1008 :
1009 : !!! the cycle by the nearest neighbors (nntot) for each kpoint
1010 :
1011 0 : do 15 ikpt_b = 1,nntot
1012 0 : kptibz_b=bpt(ikpt_b,ikpt)
1013 :
1014 0 : if(wann%l_bzsym) oper_b=mapkoper(kptibz_b)
1015 0 : if (wann%l_bzsym) kptibz_b=irreduc(kptibz_b)
1016 :
1017 : n_start=1
1018 0 : n_end=input%neig
1019 :
1020 0 : eigg = 0.
1021 0 : call cpu_time(t0)
1022 : ! call judft_error("BUG not implement")
1023 : ! CALL cdn_read(
1024 : ! > eig_id,
1025 : ! > nvd,jspd,irank,isize,kptibz_b,jspin,nbasfcn, !wannierspin instead of jspd?
1026 : ! > l_ss,l_noco,neigd,n_start,n_end,
1027 : ! < nmat_b,nv_b,ello,evdu,epar,
1028 : ! < wk_b,nbands_b,eigg,zzMat)
1029 :
1030 :
1031 :
1032 : call lapw_b%init(input,noco,nococonv,kpts,atoms,sym,kptibz_b,
1033 0 : & cell,fmpi)
1034 :
1035 :
1036 : nbasfcn_b = MERGE(lapw_b%nv(1)+lapw_b%nv(2)+2*atoms%nlotot,
1037 0 : & lapw_b%nv(1)+atoms%nlotot,noco%l_noco)
1038 0 : CALL zMat_b%init(l_real,nbasfcn_b,input%neig)
1039 0 : CALL zzMat%init(l_real,nbasfcn_b,input%neig)
1040 :
1041 :
1042 : CALL cdn_read(
1043 : & eig_id,
1044 : & lapw%dim_nvd(),input%jspins,fmpi%irank,fmpi%isize, !wannierspin instead of DIMENSION%jspd?
1045 : & kptibz_b,jspin,lapw%dim_nbasfcn(),
1046 : & noco%l_ss,noco%l_noco,input%neig,n_start,n_end,
1047 0 : & nbands_b,eigg,zzMat)
1048 :
1049 :
1050 :
1051 :
1052 0 : nslibd_b = 0
1053 :
1054 0 : IF(zzMat%l_real) THEN
1055 0 : zMat_b%data_r = 0.0
1056 : ELSE
1057 0 : zMat_b%data_c = CMPLX(0.0,0.0)
1058 : END IF
1059 :
1060 0 : eig_b(:) = 0.
1061 :
1062 0 : do i = 1,nbands_b
1063 : if((eigg(i).ge.e1s.and.nslibd_b.lt.numbands
1064 : & .and.wann%l_bynumber)
1065 : &.or.(eigg(i).ge.e1s.and.eigg(i).le.e2s.and.wann%l_byenergy)
1066 0 : &.or.(i.ge.wann%band_min(jspin)
1067 : & .and.
1068 : & (i.le.wann%band_max(jspin))
1069 : & .and.
1070 0 : & wann%l_byindex))then
1071 0 : nslibd_b = nslibd_b + 1
1072 0 : eig_b(nslibd_b) = eigg(i)
1073 0 : we_b(nslibd_b) = we_b(i)
1074 : if(l_noco)then
1075 : ! funbas = nv_b(1) + nlotot
1076 : ! funbas = funbas+nv_b(2) + nlotot
1077 : else
1078 : ! funbas = nv_b(jspin) + nlotot
1079 : endif
1080 0 : IF (zzMat%l_real) THEN
1081 : ! do j = 1,funbas
1082 0 : zMat_b%data_r(:,nslibd_b) = zzMat%data_r(:,i)
1083 : ! enddo
1084 : ELSE
1085 : ! do j = 1,funbas
1086 : ! zMat_b%data_c(j,nslibd_b) = zzMat%data_c(j,i)
1087 0 : zMat_b%data_c(:,nslibd_b) = zzMat%data_c(:,i)
1088 : ! enddo
1089 : END IF
1090 : endif
1091 : enddo
1092 :
1093 : c***********************************************************
1094 : c Rotate the wavefunction of next neighbor.
1095 : c***********************************************************
1096 : if (wann%l_bzsym .and. (oper_b.ne.1) ) then
1097 : ! TODO
1098 : ! call wann_kptsrotate(
1099 : ! > natd,nlod,llod,
1100 : ! > ntypd,nlo,llo,invsat,
1101 : ! > l_noco,l_soc,
1102 : ! > ntype,neq,nlotot,
1103 : ! > kveclo_b,jspin,
1104 : ! > oper_b,nop,mrot,nvd,
1105 : ! > nv_b,
1106 : ! > shiftkpt(:,bpt(ikpt_b,ikpt)),
1107 : ! > tau,
1108 : ! x bkpt_b,k1_b(:,:),
1109 : ! x k2_b(:,:),k3_b(:,:),
1110 : ! x zMat_b,nsfactor_b)
1111 : else
1112 : nsfactor_b=cmplx(1.0,0.0)
1113 : endif
1114 0 : noccbd_b = nslibd_b
1115 0 : call cpu_time(t1)
1116 0 : t_eig = t_eig + t1 - t0
1117 :
1118 0 : addnoco = 0
1119 : ! if(l_noco.and.jspin.eq.2) addnoco=nv_b(1)+nlotot
1120 0 : if(l_noco.and.jspin.eq.2) addnoco=lapw_b%nv(1)+atoms%nlotot
1121 :
1122 : ! set up a(k+b1),b(k+b1),c(k+b1)
1123 0 : allocate( acof_b(noccbd_b,0:lmd,natd),
1124 0 : > bcof_b(noccbd_b,0:lmd,natd),
1125 0 : > ccof_b(-llod:llod,noccbd_b,nlod,natd) )
1126 :
1127 0 : call cpu_time(t0)
1128 :
1129 : ! ALLOCATE(lapw_b%k1(SIZE(k1_b,1),SIZE(k1_b,2)))
1130 : ! ALLOCATE(lapw_b%k2(SIZE(k1_b,1),SIZE(k1_b,2)))
1131 : ! ALLOCATE(lapw_b%k3(SIZE(k1_b,1),SIZE(k1_b,2)))
1132 : ! lapw_b%k1 = k1_b
1133 : ! lapw_b%k2 = k2_b
1134 : ! lapw_b%k3 = k3_b
1135 : ! lapw_b%nmat = nmat_b
1136 : ! lapw_b%nv = nv_b
1137 : ! I think the other variables of lapw are not needed here.
1138 : !TODO
1139 : ! CALL abcof(input,atoms,noccbd_b,sym,cell,bkpt_b,lapw_b,
1140 : ! + noccbd_b,usdus,noco,jspin,kveclo_b ,
1141 : ! + acof_b,bcof_b,ccof_b,zMat_b)
1142 :
1143 : CALL abcof(input,atoms,sym,cell,lapw_b,
1144 : & noccbd_b,usdus,noco,nococonv,jspin ,
1145 0 : & acof_b,bcof_b,ccof_b,zMat_b)
1146 :
1147 :
1148 : ! DEALLOCATE(lapw_b%k1,lapw_b%k2,lapw_b%k3)
1149 :
1150 : call wann_abinv(atoms,sym,
1151 0 : X acof_b,bcof_b,ccof_b)
1152 0 : call cpu_time(t1)
1153 0 : t_abcof = t_abcof + t1 - t0
1154 :
1155 0 : do 25 ikpt_b2 = 1,nntot!ikpt_b!nntot
1156 :
1157 0 : kptibz_b2=bpt(ikpt_b2,ikpt)
1158 0 : if(wann%l_bzsym) oper_b2=mapkoper(kptibz_b2)
1159 0 : if (wann%l_bzsym) kptibz_b2=irreduc(kptibz_b2)
1160 :
1161 : n_start=1
1162 0 : n_end=input%neig
1163 :
1164 0 : eigg = 0.
1165 0 : call cpu_time(t0)
1166 :
1167 : ! WRITE(*,*) 'Here probably the wrong record is read in'
1168 : ! WRITE(*,*) 'Should eig_id not be dependent on iqpt_b?'
1169 : ! WRITE(*,*) '(in wann_uHu)'
1170 : ! CALL cdn_read(
1171 : ! > eig_id,
1172 : ! > nvd,jspd,irank,isize,kptibz_b2,jspin,nbasfcn, !wannierspin instead of jspd?
1173 : ! > l_ss,l_noco,neigd,n_start,n_end,
1174 : ! < nmat_b2,nv_b2,ello,evdu,epar,kveclo_b2,
1175 : ! < k1_b2,k2_b2,k3_b2,bkpt_b2,wk_b2,nbands_b2,
1176 : ! < eigg,zzMat)
1177 :
1178 : call lapw_b2%init(input,noco,nococonv,kpts,atoms,sym,kptibz_b2,
1179 0 : & cell,fmpi)
1180 :
1181 : nbasfcn_b2 = MERGE(lapw_b2%nv(1)+lapw_b2%nv(2)+2*atoms%nlotot,
1182 0 : & lapw_b2%nv(1)+atoms%nlotot,noco%l_noco)
1183 0 : CALL zMat_b2%init(l_real,nbasfcn_b2,input%neig)
1184 0 : CALL zzMat%init(l_real,nbasfcn_b2,input%neig)
1185 :
1186 :
1187 :
1188 : CALL cdn_read(
1189 : & eig_id,
1190 : & lapw%dim_nvd(),input%jspins,fmpi%irank,fmpi%isize, !wannierspin instead of DIMENSION%jspd?&
1191 : & kptibz_b2,jspin,lapw%dim_nbasfcn(),
1192 : & noco%l_ss,noco%l_noco,input%neig,n_start,n_end,
1193 0 : & nbands_b2,eigg,zzMat)
1194 :
1195 :
1196 0 : nslibd_b2 = 0
1197 :
1198 :
1199 0 : eig_b2(:) = 0.
1200 :
1201 0 : do i = 1,nbands_b2
1202 : if((eigg(i).ge.e1s.and.nslibd_b2.lt.numbands
1203 : & .and.wann%l_bynumber)
1204 : &.or.(eigg(i).ge.e1s.and.eigg(i).le.e2s.and.wann%l_byenergy)
1205 0 : &.or.(i.ge.wann%band_min(jspin_b)
1206 : & .and.
1207 : & (i.le.wann%band_max(jspin_b))
1208 : & .and.
1209 0 : & wann%l_byindex))then
1210 0 : nslibd_b2 = nslibd_b2 + 1
1211 0 : eig_b2(nslibd_b2) = eigg(i)
1212 0 : we_b2(nslibd_b2) = we_b2(i)
1213 : if(l_noco)then
1214 : ! funbas = nv_b2(1) + nlotot
1215 : ! funbas = funbas+nv_b2(2) + nlotot
1216 : else
1217 : ! funbas = nv_b2(jspin_b) + nlotot
1218 : endif
1219 0 : IF (zzMat%l_real) THEN
1220 : ! do j = 1,funbas
1221 0 : zMat_b2%data_r(:,nslibd_b2) = zzMat%data_r(:,i)
1222 : ! enddo
1223 : ELSE
1224 : ! do j = 1,funbas
1225 0 : zMat_b2%data_c(:,nslibd_b2) = zzMat%data_c(:,i)
1226 : ! enddo
1227 : END IF
1228 : endif
1229 : enddo
1230 :
1231 : c***********************************************************
1232 : c Rotate the wavefunction of next neighbor.
1233 : c***********************************************************
1234 : if (wann%l_bzsym .and. (oper_b2.ne.1) ) then
1235 : ! TODO
1236 : ! call wann_kptsrotate(
1237 : ! > natd,nlod,llod,
1238 : ! > ntypd,nlo,llo,invsat,
1239 : ! > l_noco,l_soc,
1240 : ! > ntype,neq,nlotot,
1241 : ! > kveclo_b2,jspin_b,
1242 : ! > oper_b2,nop,mrot,nvd,
1243 : ! > nv_b2,
1244 : ! > shiftkpt(:,bpt(ikpt_b2,ikpt)),
1245 : ! > tau,
1246 : ! x bkpt_b2,k1_b2(:,:),
1247 : ! x k2_b2(:,:),k3_b2(:,:),
1248 : ! x zMat_b2,nsfactor_b2)
1249 : else
1250 : nsfactor_b2=cmplx(1.0,0.0)
1251 : endif
1252 0 : noccbd_b2 = nslibd_b2
1253 0 : call cpu_time(t1)
1254 0 : t_eig = t_eig + t1 - t0
1255 :
1256 0 : addnoco2 = 0
1257 : ! if(l_noco.and.jspin_b.eq.2) addnoco2=nv_b2(1)+nlotot
1258 0 : if(l_noco.and.jspin_b.eq.2) addnoco2=lapw_b2%nv(1)+atoms%nlotot
1259 :
1260 : ! set up a(k+b2),b(k+b2),c(k+b2)
1261 0 : allocate( acof_b2(noccbd_b2,0:lmd,natd),
1262 0 : > bcof_b2(noccbd_b2,0:lmd,natd),
1263 0 : > ccof_b2(-llod:llod,noccbd_b2,nlod,natd) )
1264 :
1265 0 : call cpu_time(t0)
1266 :
1267 :
1268 : ! ALLOCATE(lapw_b2%k1(SIZE(k1_b2,1),SIZE(k1_b2,2)))
1269 : ! ALLOCATE(lapw_b2%k2(SIZE(k1_b2,1),SIZE(k1_b2,2)))
1270 : ! ALLOCATE(lapw_b2%k3(SIZE(k1_b2,1),SIZE(k1_b2,2)))
1271 : ! lapw_b2%k1 = k1_b2
1272 : ! lapw_b2%k2 = k2_b2
1273 : ! lapw_b2%k3 = k3_b2
1274 : ! lapw_b2%nmat = nmat_b2
1275 : ! lapw_b2%nv = nv_b2
1276 : ! I think the other variables of lapw are not needed here.
1277 : !TODO
1278 : ! CALL abcof(input,atoms,noccbd_b2,sym,cell,bkpt_b2,lapw_b2,
1279 : ! + noccbd_b2,usdus,noco,jspin_b,kveclo_b2 ,
1280 : ! + acof_b2,bcof_b2,ccof_b2,zMat_b2)
1281 :
1282 : ! DEALLOCATE(lapw_b2%k1,lapw_b2%k2,lapw_b2%k3)
1283 :
1284 : CALL abcof(input,atoms,sym,cell,lapw_b2,
1285 : & noccbd_b2,usdus,noco,nococonv,jspin ,
1286 0 : & acof_b2,bcof_b2,ccof_b2,zMat_b2)
1287 :
1288 : call wann_abinv(atoms,sym,
1289 0 : X acof_b2,bcof_b2,ccof_b2)
1290 0 : call cpu_time(t1)
1291 0 : t_abcof = t_abcof + t1 - t0
1292 :
1293 :
1294 : c**************************************c
1295 : c calculate uHu matrix due to: c
1296 : c (i) interstitial c
1297 : c (ii) muffin tins c
1298 : c (iii) vacuum c
1299 : c**************************************c
1300 : c if(ANY(gb(3:3,ikpt_b,ikpt).ne.0) .or.
1301 : c > ANY(gb(3:3,ikpt_b2,ikpt).ne.0)) then
1302 : c write(*,*)ikpt,ikpt_b,ikpt_b2
1303 : c write(*,*)gb(:,ikpt_b,ikpt)
1304 : c write(*,*)gb(:,ikpt_b2,ikpt)
1305 : c uHu(:,:,ikpt_b2,ikpt_b,i_rec) = cmplx(0.,0.)
1306 : c goto 848
1307 : c endif
1308 :
1309 : ! MT (SPH, NON)
1310 0 : if(.not.(l_skip_sph.and.l_skip_non)) then
1311 0 : call cpu_time(t0)
1312 : call wann_uHu_sph(
1313 : > cmplx(1.,0.),nbnd,llod,nslibd_b,nslibd_b2,nlod,natd,ntypd,
1314 : > lmd,jmtd,taual,nop,lmax,ntype,neq,nlo,llo,
1315 : > acof_b,bcof_b,ccof_b,lapw_b2%bkpt,
1316 : > acof_b2,bcof_b2,ccof_b2,lapw_b%bkpt,lapw%bkpt,
1317 : > gb(:,ikpt_b,ikpt),gb(:,ikpt_b2,ikpt),
1318 : < tuu(:,:,:,:,tspin2),tud(:,:,:,:,tspin2),
1319 : > tdu(:,:,:,:,tspin2),tdd(:,:,:,:,tspin2),
1320 : > tuulo(:,:,:,:,:,tspin2),tulou(:,:,:,:,:,tspin2),
1321 : > tdulo(:,:,:,:,:,tspin2),tulod(:,:,:,:,:,tspin2),
1322 : > tuloulo(:,:,:,:,:,:,tspin2),
1323 : > kdiff,kdiff,nntot,nntot,
1324 0 : > uHu(:,:,ikpt_b2,ikpt_b,i_rec))
1325 0 : call cpu_time(t1)
1326 0 : t_sph = t_sph + t1 - t0
1327 : endif
1328 : c if(ikpt.eq.1) then
1329 : c write(*,*)'SPH,NON',ikpt_b2,ikpt_b
1330 : c write(*,*)uHu(1,3,ikpt_b2,ikpt_b,i_rec)
1331 : c write(*,*)uHu(3,1,ikpt_b2,ikpt_b,i_rec)
1332 : c uHuold(:,:) = uHu(:,:,ikpt_b2,ikpt_b,i_rec)
1333 : c endif
1334 :
1335 : ! MT (SOC)
1336 0 : if(l_soc.and.(.not.l_skip_soc)) then
1337 0 : call cpu_time(t0)
1338 : call wann_uHu_sph(
1339 : > cmplx(1.,0.),nbnd,llod,nslibd_b,nslibd_b2,nlod,natd,ntypd,
1340 : > lmd,jmtd,taual,nop,lmax,ntype,neq,nlo,llo,
1341 : > acof_b,bcof_b,ccof_b,lapw_b2%bkpt,
1342 : > acof_b2,bcof_b2,ccof_b2,lapw_b%bkpt,lapw%bkpt,
1343 : > gb(:,ikpt_b,ikpt),gb(:,ikpt_b2,ikpt),
1344 : < tuu_soc,tud_soc,
1345 : > tdu_soc,tdd_soc,
1346 : > tuulo_soc,tulou_soc,
1347 : > tdulo_soc,tulod_soc,
1348 : > tuloulo_soc,
1349 : > kdiff,kdiff,nntot,nntot,
1350 0 : > uHu(:,:,ikpt_b2,ikpt_b,i_rec))
1351 0 : call cpu_time(t1)
1352 0 : t_sph = t_sph + t1 - t0
1353 : c if(ikpt.eq.1) then
1354 : c write(*,*)'SOC',ikpt_b2,ikpt_b
1355 : c write(*,*)uHu(1,3,ikpt_b2,ikpt_b,i_rec)-uHuold(1,3)
1356 : c write(*,*)uHu(3,1,ikpt_b2,ikpt_b,i_rec)-uHuold(3,1)
1357 : c uHuold(:,:) = uHu(:,:,ikpt_b2,ikpt_b,i_rec)
1358 : c endif
1359 : endif
1360 :
1361 :
1362 : ! In collinear calculation, the potential in interstital and
1363 : ! vacuum region is always diagonal with respect to the spin
1364 0 : IF( l_noco .OR. ((.NOT. l_noco) .AND. (doublespin.LT.3)) ) THEN
1365 :
1366 0 : jspin3 = doublespin
1367 0 : if(jspins.eq.1) jspin3=1
1368 :
1369 0 : jspin4 = jspin
1370 0 : jspin4_b=jspin_b
1371 0 : if(jspins.eq.1) then
1372 0 : jspin4=1
1373 0 : jspin4_b=1
1374 : endif
1375 :
1376 0 : sign2 = 1
1377 0 : if(doublespin.EQ.4) sign2=-1
1378 :
1379 : ! INT
1380 0 : if(.not.l_skip_int) then
1381 0 : call cpu_time(t0)
1382 : call wann_uHu_int(cmplx(1.,0.),nvd,k1d,k2d,k3d,n3d,
1383 : > lapw_b%nv(jspin),!nv_b(jspin),
1384 : > lapw_b2%nv(jspin),!nv_b2(jspin_b),
1385 : > nbnd,neigd,
1386 : > nslibd_b,nslibd_b2,nbasfcn_b,nbasfcn_b2,addnoco,addnoco2,
1387 : ! > k1_b(:,jspin), k2_b(:,jspin), k3_b(:,jspin),
1388 : > lapw_b%k1(:,jspin),lapw_b%k2(:,jspin),lapw_b%k3(:,jspin),
1389 : > gb(:,ikpt_b,ikpt),
1390 : ! > k1_b2(:,jspin_b),k2_b2(:,jspin_b),k3_b2(:,jspin_b),
1391 : > lapw_b2%k1(:,jspin),lapw_b2%k2(:,jspin),lapw_b2%k3(:,jspin),
1392 : > gb(:,ikpt_b2,ikpt),
1393 : > lapw%bkpt,bbmat,vpw(:,jspin3),zMat_b,zMat_b2,rgphs,
1394 : > ustep,ig,jspin.eq.jspin_b,sign2,
1395 0 : > uHu(:,:,ikpt_b2,ikpt_b,i_rec))
1396 0 : call cpu_time(t1)
1397 0 : t_int = t_int + t1 - t0
1398 : endif
1399 : c if(ikpt.eq.1) then
1400 : c write(*,*)'INT',ikpt_b2,ikpt_b
1401 : c write(*,*)uHu(1,3,ikpt_b2,ikpt_b,i_rec)-uHuold(1,3)
1402 : c write(*,*)uHu(3,1,ikpt_b2,ikpt_b,i_rec)-uHuold(3,1)
1403 : c uHuold(:,:) = uHu(:,:,ikpt_b2,ikpt_b,i_rec)
1404 : c endif
1405 :
1406 :
1407 : ! VAC
1408 0 : if ((.not.l_skip_vac) .and. film ) then
1409 :
1410 0 : call cpu_time(t0)
1411 : call wann_uHu_vac(
1412 : > cmplx(1.,0.),l_noco,l_soc,jspins,nlotot,qpt_i,
1413 : > nbnd,z1,nmzxyd,nmzd,n2d,nv2d,k1d,k2d,k3d,n3d,nvac,ig,
1414 : > rgphs,nmzxy,nmz,delz,ig2,nq2,kv2,area,bmat,bbmat,
1415 : > evac(:,jspin4),evac(:,jspin4_b),lapw_b%bkpt,lapw_b2%bkpt,
1416 : > vzxy(:,:,:,jspin3),vz,nslibd_b,nslibd_b2,
1417 : > jspin,jspin_b,doublespin,
1418 : ! > k1_b,k2_b,k3_b,
1419 : > lapw_b%k1,lapw_b%k2,lapw_b%k3,
1420 : ! > k1_b2,k2_b2,k3_b2,
1421 : > lapw_b2%k1,lapw_b2%k2,lapw_b2%k3,
1422 : > wannierspin,nvd,nbasfcn,neigd,
1423 : > zMat_b,zMat_b2,
1424 : > lapw_b%nv,!nv_b,
1425 : > lapw_b2%nv,! nv_b2,
1426 : > omtil,gb(:,ikpt_b,ikpt),
1427 : > gb(:,ikpt_b2,ikpt),sign2,
1428 0 : > uHu(:,:,ikpt_b2,ikpt_b,i_rec))
1429 0 : call cpu_time(t1)
1430 0 : t_vac = t_vac + t1 - t0
1431 :
1432 :
1433 : endif
1434 : c if(ikpt.eq.1) then
1435 : c write(*,*)'VAC',ikpt_b2,ikpt_b
1436 : c write(*,*)uHu(1,3,ikpt_b2,ikpt_b,i_rec)-uHuold(1,3)
1437 : c write(*,*)uHu(3,1,ikpt_b2,ikpt_b,i_rec)-uHuold(3,1)
1438 : c uHuold(:,:) = uHu(:,:,ikpt_b2,ikpt_b,i_rec)
1439 : c endif
1440 : ELSE
1441 0 : IF(ikpt.EQ.1 .AND. ikpt_b.EQ.1 .AND. ikpt_b2.EQ.1) THEN
1442 0 : if(l_p0) write(*,*)'skip INT and VAC'
1443 : ENDIF
1444 : ENDIF
1445 :
1446 :
1447 : 848 continue
1448 :
1449 0 : deallocate (acof_b2,bcof_b2,ccof_b2)
1450 0 : 25 continue ! end of loop by the nearest k-neighbors
1451 :
1452 0 : deallocate (acof_b,bcof_b,ccof_b)
1453 0 : 15 continue ! end of loop by the nearest k-neighbors
1454 0 : deallocate ( we_b,we_b2 )
1455 0 : deallocate ( eigg )
1456 :
1457 :
1458 :
1459 :
1460 :
1461 : endif ! loop by processors
1462 :
1463 :
1464 :
1465 :
1466 :
1467 :
1468 0 : 10 continue ! end of cycle by the k-points
1469 :
1470 : #ifdef CPP_MPI
1471 0 : call MPI_BARRIER(MPI_COMM_WORLD,ierr(1))
1472 : #endif
1473 : 5 continue
1474 :
1475 0 : l_symcheck = (.not.l_nocosoc).or.(doublespin.eq.doublespin_max)
1476 0 : jspin4=jspin
1477 0 : if(l_nocosoc) jspin4=1
1478 0 : if(l_p0) write(*,*)'write uHu file'
1479 0 : uHu = uHu * htr2ev
1480 : call wann_write_uHu(jspin4,l_p0,fullnkpts,nntot,nntot,wann,
1481 : > nbnd,bpt,gb,isize,irank,fending,'',uHu,
1482 : > counts(irank),counts,displs,isize,
1483 0 : > wann%l_unformatted,.true.,l_symcheck)
1484 0 : uHu = uHu / htr2ev
1485 0 : if(.not.l_nocosoc) deallocate(uHu)
1486 :
1487 0 : if(.not.l_noco)nrec=nrec+nkpts
1488 :
1489 0 : 110 continue ! end of cycle by spins
1490 :
1491 : #ifdef CPP_MPI
1492 0 : call MPI_BARRIER(MPI_COMM_WORLD,ierr(1))
1493 : #endif
1494 :
1495 : ! close eig files
1496 : IF (l_gwf) THEN
1497 :
1498 : ! CALL close_eig(eig_id)
1499 :
1500 : ! CLOSE(66)
1501 : DO iqpt_b=1,nntot_q
1502 : ! CALL close_eig(innerEig_idList(iqpt_b))
1503 : ENDDO
1504 : ENDIF
1505 :
1506 :
1507 : c if(l_p0) write(*,*)'write uHu file'
1508 : c uHu = htr2ev * uHu
1509 : c call wann_write_uHu(1,l_p0,fullnkpts,nntot,nntot,wann,
1510 : c > nbnd,bpt,gb,isize,irank,fending,'',uHu,
1511 : c > counts(irank),counts,displs,isize,
1512 : c > wann%l_unformatted,.true.)
1513 0 : if(allocated(uHu)) deallocate(uHu)
1514 : c if(allocated(uHuold)) deallocate(uHuold)
1515 :
1516 0 : DEALLOCATE(innerEig_idList)
1517 :
1518 0 : 314 continue ! iqpt, q-points
1519 : c************************************************c
1520 : c END Q LOOP c
1521 : c************************************************c
1522 :
1523 :
1524 : ! deallocate ( kveclo,nv,k1,k2,k3 )
1525 0 : deallocate(flo)
1526 0 : deallocate ( ff,gg)
1527 : ! if (allocated(nv_b))deallocate(kveclo_b,nv_b,k1_b,k2_b,k3_b)
1528 : ! if (allocated(nv_b2))deallocate(kveclo_b2,nv_b2,k1_b2,k2_b2,k3_b2)
1529 0 : if (wann%l_bzsym)deallocate(irreduc,mapkoper,shiftkpt)
1530 0 : if (wann%l_bzsym.AND.l_gwf)deallocate(irreduc_q,mapqoper,shiftqpt)
1531 0 : if (allocated(pair_to_do)) deallocate(pair_to_do,maptopair)
1532 0 : if (allocated(pair_to_do_q)) deallocate(pair_to_do_q,maptopair_q)
1533 0 : if (allocated(kdiff)) deallocate ( kdiff )
1534 0 : if (allocated(qdiff)) deallocate(qdiff,zero_qdiff)
1535 0 : if(allocated(vpw)) deallocate(vpw)
1536 0 : if(allocated(vzxy)) deallocate(vzxy)
1537 0 : if(allocated(vz)) deallocate(vz)
1538 0 : deallocate(vr)
1539 0 : deallocate(tdd,tdu,tud,tuu)
1540 0 : deallocate(tdulo,tuulo)
1541 0 : deallocate(tulou,tulod)
1542 0 : deallocate(tuloulo)
1543 0 : deallocate(tdd_soc,tdu_soc,tud_soc,tuu_soc)
1544 0 : deallocate(tdulo_soc,tuulo_soc)
1545 0 : deallocate(tulou_soc,tulod_soc)
1546 0 : deallocate(tuloulo_soc)
1547 0 : deallocate(counts,displs)
1548 :
1549 : 9110 continue
1550 :
1551 0 : if(l_sgwf .or. l_socgwf) gb_q = gb_q/2
1552 :
1553 0 : if(l_p0.and.l_gwf) then
1554 : call wann_uHu_commat(
1555 : > fullnkpts,nntot,bpt,fullnqpts,nntot_q,bpt_q,
1556 : > gb,gb_q,aux_latt_const,wann%l_unformatted,l_dim,
1557 0 : > nparampts,param_vec/2.0)
1558 : endif
1559 :
1560 0 : if(allocated(gb)) deallocate(gb,bpt)
1561 0 : if(allocated(gb_q)) deallocate(gb_q,bpt_q)
1562 :
1563 : #ifdef CPP_MPI
1564 0 : call MPI_BARRIER(MPI_COMM_WORLD,ierr(1))
1565 : #endif
1566 0 : call cpu_time(t1)
1567 0 : t_total = t1-t00
1568 0 : if(l_p0) then
1569 0 : write(*,900)'t_init =',t_init
1570 0 : write(*,900)'t_myTlmplm=',t_myTlmplm
1571 0 : write(*,900)'t_eig =',t_eig
1572 0 : write(*,900)'t_abcof =',t_abcof
1573 0 : write(*,900)'t_int =',t_int
1574 0 : write(*,900)'t_sph =',t_sph
1575 0 : write(*,900)'t_vac =',t_vac
1576 0 : write(*,900)'t_total =',t_total
1577 : endif
1578 :
1579 : 900 FORMAT(a,f14.4)
1580 :
1581 0 : END SUBROUTINE wann_uHu
1582 :
1583 :
1584 : END MODULE m_wann_uHu
|
