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