Line data Source code
1 : !--------------------------------------------------------------------------------
2 : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
3 : ! This file is part of FLEUR and available as free software under the conditions
4 : ! of the MIT license as expressed in the LICENSE file in more detail.
5 : !--------------------------------------------------------------------------------
6 :
7 : c******************************************************************
8 : c Calculate lapw-representation of wannierfunctions.
9 : c Lapw means here: lapw-like, but not identical to
10 : c fleur-lapw. Wannierfunctions can not be expanded
11 : c in fleur-lapw-set of a single k-point in general.
12 : c
13 : c Frank Freimuth, November 2006
14 : c******************************************************************
15 : module m_wannier_to_lapw
16 : use m_juDFT
17 : #ifdef CPP_MPI
18 : use mpi
19 : #endif
20 : contains
21 0 : subroutine wannier_to_lapw(
22 : > mpi_communicatior,eig_id,l_real,
23 : > input,lapw ,noco,nococonv,sym,cell,atoms,stars,vacuum,
24 : > sphhar,
25 : > vTot,
26 : > l_soc,unigrid,sortrule,band_min,band_max,
27 : > 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 : > 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 : > ustep,ig,k1d,k2d,k3d,rgphs,slice,kk,nnne,
35 0 : > z1,nv2d,nmzxy,nmz,delz,ig2,area,tau,zatom,nq2,nop2,
36 0 : > volint,symor,pos,ef,l_bzsym,l_proj_plot,wan90version)
37 : use m_wann_rw_eig
38 : use m_wann_read_umatrix
39 : use m_abcof
40 : use m_radfun
41 : use m_radflo
42 : use m_cdnread, only : cdn_read0
43 : use m_types
44 : use m_constants
45 : use m_wann_real
46 : use m_xsf_io
47 : use m_wann_plot_vac
48 : USE m_abcrot
49 :
50 :
51 : implicit none
52 :
53 : TYPE(t_input),INTENT(IN) :: input
54 : TYPE(t_lapw),INTENT(IN) :: lapw
55 :
56 : TYPE(t_noco),INTENT(IN) :: noco
57 : TYPE(t_nococonv),INTENT(IN):: nococonv
58 : TYPE(t_sym),INTENT(IN) :: sym
59 : TYPE(t_cell),INTENT(IN) :: cell
60 : TYPE(t_atoms),INTENT(IN) :: atoms
61 : TYPE(t_stars),INTENT(IN) :: stars
62 : TYPE(t_vacuum),INTENT(IN) :: vacuum
63 : TYPE(t_sphhar),INTENT(IN) :: sphhar
64 : TYPE(t_potden),INTENT(IN) :: vTot
65 :
66 : #ifdef CPP_MPI
67 : integer mpiierr(3)
68 : integer cpu_index
69 : integer stt(MPI_STATUS_SIZE)
70 : #endif
71 : logical,intent(in):: l_soc, l_real
72 : integer,intent(in)::unigrid(4),mpi_communicatior,eig_id
73 : integer,intent(in)::band_min(2),band_max(2)
74 : logical, intent (in) :: invs,invs2,film,slice,symor
75 : integer, intent (in) :: lmaxd,ntypd,neigd,nkptd,kk,nnne
76 : integer, intent (in) :: natd,nop,nvd,jspd,nbasfcn,nq2,nop2
77 : integer, intent (in) :: llod,nlod,ntype,n3d,n2d
78 : integer, intent (in) :: nmzxyd,nmzd,jmtd,nlhd,nq3,nvac
79 : integer, intent (in) :: ntypsd,jspins,k1d,k2d,k3d
80 : real, intent (in) :: omtil,e1s,e2s,delz,area,z1,volint
81 : integer, intent (in) :: ig(-k1d:k1d,-k2d:k2d,-k3d:k3d)
82 : complex, intent (in) :: rgphs(-k1d:k1d,-k2d:k2d,-k3d:k3d)
83 : integer, intent (in) :: nlh(ntypsd),jri(ntypd),ntypsy(natd)
84 : integer, intent (in) :: nlo(ntypd),llo(nlod,ntypd),lapw_l(ntypd)
85 : integer, intent (in) :: invtab(nop),mrot(3,3,nop),ngopr(natd)
86 : integer, intent (in) :: neq(ntypd),lmax(ntypd)
87 : integer, intent (in) :: invsat(natd),invsatnr(natd),nkpt
88 : integer, intent (in) :: irank,isize,nv2d,nmzxy,nmz
89 : integer, intent (in) :: ulo_der(nlod,ntypd),ig2(n3d)
90 : real, intent (in) :: taual(3,natd),rmt(ntypd),dx(ntypd)
91 : real, intent (in) :: amat(3,3),bmat(3,3),bbmat(3,3)
92 : real, intent (in) :: rmsh(jmtd,ntypd),tau(3,nop),zatom(ntype)
93 : real, intent (in) :: alph(ntypd),beta(ntypd),qss(3)
94 : real, intent (in) :: pos(3,natd),ef
95 : complex, intent (in) :: ustep(n3d)
96 : logical, intent (in) :: l_dulo(nlod,ntypd),l_noco,l_ss,l_bzsym
97 : logical,intent(in)::l_proj_plot
98 : integer,intent(in)::sortrule
99 : integer, intent(in):: wan90version
100 : c-odim
101 : real, intent (in) :: sk2(n2d),phi2(n2d)
102 :
103 : c+odim
104 : logical l_spreadcal
105 : complex, allocatable::spreads(:,:)
106 : real,allocatable:: centers(:,:)
107 : cccccccccccccccccc local variables cccccccccccccccccccc
108 : integer lmd,nlotot,n,nmat,nw,ispin,iter,ikpt,ilo
109 : integer :: wannierspin,jspin2
110 : integer noccbd,nn,nkpts,i,jspin,j,l,i_rec,m,nwf,nwfp
111 : integer jsp_start,jsp_end,nrec,nrec1,nbands
112 : integer nodeu,noded,n_size,na,n_rank,nbnd,nkbnd
113 : integer i1,i2,i3,in,ikpt_k,lda
114 : integer n_bands(0:neigd),nslibd
115 : character*8 dop,iop,name(10)
116 : real bkpt(3),wronk,wk,wk_b,phase
117 0 : real eig(neigd),cp_time(9)
118 : logical l_p0,l_bkpts,l_proj,l_amn,l_mmn,l_eig,wann,wann_plott
119 : !!! energy window boundaries
120 0 : integer, allocatable :: nv(:)
121 0 : integer, allocatable :: k1(:,:),k2(:,:),k3(:,:)
122 0 : real, allocatable :: we(:)
123 :
124 0 : real, allocatable :: eigg(:)
125 : real kpoints(nkptd)
126 : !!! a and b coeff. constructed for each k-point
127 0 : complex, allocatable :: acof(:,:,:)
128 0 : complex, allocatable :: bcof(:,:,:)
129 0 : complex, allocatable :: ccof(:,:,:,:)
130 0 : complex, allocatable :: wann_acof(:,:,:,:,:,:)
131 0 : complex, allocatable :: wann_bcof(:,:,:,:,:,:)
132 0 : complex, allocatable :: wann_ccof(:,:,:,:,:,:,:)
133 : !!! the parameters for the number of wfs
134 : integer :: nwfs
135 : !!! the potential in the spheres and the vacuum
136 0 : real, allocatable :: vr(:,:,:)
137 : !!! bkpts data
138 : integer nntot,ikpt_help
139 : integer, allocatable :: gb(:,:,:),bpt(:,:)
140 : !!! radial wavefunctions in the muffin-tins and more ...
141 0 : real, allocatable :: flo(:,:,:,:)
142 0 : real, allocatable :: ff(:,:,:,:),gg(:,:,:,:)
143 :
144 : real uuilon(nlod,ntypd),duilon(nlod,ntypd)
145 : real ulouilopn(nlod,nlod,ntypd)
146 : !!! energy parameters
147 : real ello(nlod,ntypd,max(2,jspd)),evac(2,max(2,jspd))
148 : real epar(0:lmaxd,ntypd,max(2,jspd)),evdu(2,max(jspd,2))
149 : character(len=3) :: spin12(2)
150 : data spin12/'WF1' , 'WF2'/
151 : complex,allocatable::wannierfunc(:,:)
152 : complex,allocatable::wannierfunc_temp(:,:)
153 : integer posi
154 0 : complex,allocatable::u_matrix(:,:,:)
155 : integer kpt,oper
156 : real poinint(3)
157 : real phas,tmax
158 : real bkrot(3)
159 : integer j1,j2,j3
160 : logical um_format
161 : logical have_disentangled
162 0 : integer,allocatable :: ndimwin(:)
163 0 : logical,allocatable :: lwindow(:,:)
164 : integer :: chk_unit,nkp,ntmp,ierr,fullnkpts
165 0 : integer,allocatable::irreduc(:),mapkoper(:)
166 : logical l_file
167 0 : logical,allocatable::inc_band(:)
168 : integer num_inc,counter,kptibz
169 : logical l_byindex, l_byenergy, l_bynumber
170 : integer num_wann,num_bands,kpun
171 : complex d_wgn(-3:3,-3:3,3,nop)
172 : integer pos1,pos2,ato,loc,invop
173 : real vz0(2)
174 : integer ik,nv2,ivac,jvac,symvac,symvacvac
175 : real evacp,sign,arg
176 : complex c_1
177 : integer kvac1(nv2d),kvac2(nv2d),map2(nvd)
178 : real fas,zks
179 : integer mesh
180 : integer n2
181 : real v(3),scale,ev
182 : complex av,bv
183 : real volume
184 : REAL :: s,const
185 : COMPLEX :: xdnout,factor
186 : INTEGER :: ii3,ix,iy,iz,nplo,nbn
187 : INTEGER :: nbmin,nbmax
188 : INTEGER :: nplot,nq,nt,jm,ii1,ii2
189 : LOGICAL :: twodim
190 0 : real,allocatable::knorm(:,:)
191 : real,allocatable::wfnorm(:)
192 : REAL :: pt(3),vec1(3),vec2(3),vec3(3),zero(3)
193 : INTEGER :: grid(3),k
194 : LOGICAL :: cartesian,xsf
195 : REAL :: rhocc(jmtd)
196 : REAL :: point(3),post(3,natd)
197 : CHARACTER(len=30):: filename
198 : CHARACTER(len=20):: name1,name2,name3
199 : CHARACTER(len=10):: vandername
200 : NAMELIST /plot/twodim,cartesian,vec1,vec2,vec3,grid,zero,filename
201 : integer cell1,cell2,cell3,pw1,pw2,pw3
202 0 : complex,allocatable::wannint(:,:,:,:)
203 0 : complex,allocatable::wannz(:,:),wannz2(:,:)
204 : real denom
205 :
206 0 : TYPE(t_mat) :: zzMat, zMat
207 0 : TYPE(t_usdus) :: usdus
208 :
209 0 : call timestart("wannier_to_lapw")
210 : c specify number of unit-cells that are calculated
211 0 : cell1=3
212 0 : cell2=3
213 0 : cell3=3
214 :
215 0 : um_format=.true.
216 0 : l_byindex=.false.
217 0 : l_byenergy=.false.
218 0 : l_bynumber=.false.
219 0 : if(sortrule==1)l_byindex=.true.
220 0 : if(sortrule==2)l_bynumber=.true.
221 0 : if(sortrule==3)l_byenergy=.true.
222 :
223 :
224 0 : lmd = lmaxd*(lmaxd+2)
225 0 : nkpts = nkpt
226 :
227 0 : nrec = 0
228 0 : nlotot = 0
229 0 : do n = 1, ntype
230 0 : do l = 1,nlo(n)
231 0 : nlotot = nlotot + neq(n) * ( 2*llo(l,n) + 1 )
232 : enddo
233 : enddo
234 :
235 :
236 : cccccccccccccccc initialize the potential cccccccccccc
237 :
238 0 : allocate (vr(jmtd,ntypd,jspd))
239 :
240 0 : do jspin = 1,jspins
241 0 : do n = 1, ntype
242 0 : do j = 1,jri(n)
243 0 : vr(j,n,jspin) = vTot%mt(j,0,n,jspin)
244 : enddo
245 : enddo
246 : enddo
247 :
248 : cccccccccccccccc end of the potential part ccccccccccc
249 0 : wannierspin=jspd
250 0 : if(l_soc) wannierspin=2
251 :
252 0 : allocate ( nv(jspd) )
253 0 : allocate ( k1(nvd,jspd),k2(nvd,jspd),k3(nvd,jspd) )
254 0 : allocate ( ff(ntypd,jmtd,2,0:lmaxd) )
255 0 : allocate ( gg(ntypd,jmtd,2,0:lmaxd) )
256 0 : allocate ( usdus%us(0:lmaxd,ntypd,jspins) )
257 0 : allocate ( usdus%uds(0:lmaxd,ntypd,jspins) )
258 0 : allocate ( usdus%dus(0:lmaxd,ntypd,jspins) )
259 0 : allocate ( usdus%duds(0:lmaxd,ntypd,jspins) )
260 0 : allocate ( usdus%ddn(0:lmaxd,ntypd,jspins) )
261 0 : allocate ( usdus%ulos(nlod,ntypd,jspins) )
262 0 : allocate ( usdus%dulos(nlod,ntypd,jspins) )
263 0 : allocate ( usdus%uulon(nlod,ntypd,jspins) )
264 0 : allocate ( usdus%dulon(nlod,ntypd,jspins) )
265 0 : allocate ( usdus%uloulopn(nlod,nlod,ntypd,jspins) )
266 : c****************************************************
267 : c cycle by spins starts!
268 : c****************************************************
269 0 : do 110 jspin=1,wannierspin ! cycle by spins
270 :
271 0 : jspin2=jspin
272 : if(l_soc .and. jspins.eq.1)jspin2=1
273 : jsp_start = jspin ; jsp_end = jspin
274 :
275 0 : jsp_start = jspin ; jsp_end = jspin
276 :
277 : c*******************************************************
278 : c get num_bands and num_wann from WF1.amn (WF2.amn)
279 : c*******************************************************
280 0 : l_file=.false.
281 0 : inquire(file=spin12(jspin)//'.amn',exist=l_file)
282 0 : open(355,file=spin12(jspin)//'.amn')
283 0 : read(355,*)
284 0 : read(355,*)num_bands,kpun,num_wann
285 0 : close(355)
286 0 : if(l_byindex.and.
287 : & .not.((1+band_max(jspin)-band_min(jspin)).eq.num_bands))
288 : & CALL juDFT_error("1+band_max-band_min/=num_bands",calledby
289 0 : + ="wannier_to_lapw")
290 :
291 : c**************************************************************
292 : ! for bzsym = .true.: determine mapping between kpts and w90kpts
293 : c**************************************************************
294 0 : if (l_bzsym) then
295 0 : l_file=.false.
296 0 : inquire(file='w90kpts',exist=l_file)
297 0 : IF(.NOT.l_file) CALL juDFT_error
298 : + ("w90kpts not found, needed if bzsym",calledby
299 0 : + ="wannier_to_lapw")
300 0 : open(412,file='w90kpts',form='formatted')
301 0 : read(412,*)fullnkpts
302 0 : close(412)
303 0 : print*,"fullnkpts=",fullnkpts
304 0 : IF(fullnkpts<=nkpts) CALL juDFT_error("fullnkpts.le.nkpts"
305 0 : + ,calledby ="wannier_to_lapw")
306 0 : allocate(irreduc(fullnkpts),mapkoper(fullnkpts))
307 0 : l_file=.false.
308 0 : inquire(file='kptsmap',exist=l_file)
309 0 : IF(.NOT.l_file) CALL juDFT_error
310 : + ("kptsmap not found, needed if bzsym",calledby
311 0 : + ="wannier_to_lapw")
312 0 : open(713,file='kptsmap')
313 0 : do i=1,fullnkpts
314 0 : read(713,*)kpt,irreduc(i),mapkoper(i)
315 0 : IF(kpt/=i) CALL juDFT_error("kpt.ne.i",calledby
316 0 : + ="wannier_to_lapw")
317 0 : print*,i,irreduc(i),mapkoper(i)
318 : enddo
319 0 : close(713)
320 0 : IF(MAXVAL(irreduc(:))/=nkpts) CALL juDFT_error
321 0 : + ("max(irreduc(:))/=nkpts",calledby ="wannier_to_lapw")
322 : else
323 0 : fullnkpts=nkpts
324 :
325 : endif
326 :
327 0 : IF(kpun/=fullnkpts) CALL juDFT_error
328 : + ("mismatch in kpun and fullnkpts",calledby ="wannier_to_lapw"
329 0 : + )
330 :
331 : c**************************************************************
332 : c read in u_matrix
333 : c**************************************************************
334 :
335 0 : if(.not.l_proj_plot)then
336 0 : allocate(lwindow(num_bands,fullnkpts))
337 0 : allocate(ndimwin(fullnkpts))
338 0 : allocate(u_matrix(num_bands,num_wann,fullnkpts))
339 : call wann_read_umatrix(fullnkpts,num_wann,num_bands,
340 : > um_format,jspin,wan90version,
341 : < have_disentangled,lwindow,ndimwin,
342 0 : < u_matrix)
343 0 : if(.not.have_disentangled)
344 0 : & deallocate(lwindow,ndimwin)
345 0 : if(have_disentangled)allocate(inc_band(num_bands))
346 :
347 : else
348 : c**************************************************************
349 : c read WF1.umn (WF2.umn) (if projmethod)
350 : c**************************************************************
351 0 : l_file=.false.
352 0 : inquire(file=spin12(jspin)//'.umn',exist=l_file)
353 0 : IF(.NOT.l_file) CALL juDFT_error("no umn file found" ,calledby
354 0 : + ="wannier_to_lapw")
355 0 : open(419,file=spin12(jspin)//'.umn')
356 0 : read(419,*) !num_wann,num_bands
357 0 : allocate(u_matrix(num_bands,num_wann,fullnkpts))
358 0 : do ikpt=1,fullnkpts
359 0 : do j=1,num_wann
360 0 : do i=1,num_bands
361 0 : read(419,*)u_matrix(i,j,ikpt)
362 : enddo
363 : enddo
364 : enddo
365 0 : close(419)
366 : endif
367 :
368 : ***********************************************************
369 : ***********************************************************
370 :
371 0 : print*,"num_wann=",num_wann
372 0 : print*,"num_bands=",num_bands
373 :
374 : cccccccccccc read in the eigenvalues and vectors cccccc
375 :
376 : l_p0 = .false.
377 0 : if (irank.eq.0) l_p0 = .true.
378 0 : allocate ( flo(ntypd,jmtd,2,nlod) )
379 : #ifdef CPP_NEVER
380 : call cdn_read0(
381 : > lmaxd,ntypd,nlod,neigd,wannierspin,
382 : > irank,isize,jspin,jsp_start,jsp_end,
383 : > l_noco,nrec,66,
384 : < ello,evac,epar,bkpt,wk,n_bands,nrec1,n_size)
385 :
386 :
387 : na = 1
388 : do 40 n = 1,ntype
389 : do 30 l = 0,lmax(n)
390 : c...compute the l-dependent, k-independent radial MT- basis functions
391 : call radfun(
392 : > l,epar(l,n,jspin),vr(1,n,jspin),jri(n),rmsh(1,n),
393 : > dx(n),jmtd,
394 : < ff(n,:,:,l),gg(n,:,:,l),us(l,n),
395 : < dus(l,n),uds(l,n),duds(l,n),
396 : < ddn(l,n),nodeu,noded,wronk)
397 : 30 continue
398 : c...and the local orbital radial functions
399 : c do ilo = 1, nlo(n)
400 : call radflo(
401 : > ntypd,nlod,jspd,jmtd,lmaxd,n,jspin,
402 : > ello(1,1,jspin),vr(1,n,jspin),
403 : > jri(n),rmsh(1,n),dx(n),ff(n,1:,1:,0:),
404 : > gg(n,1:,1:,0:),llo,nlo,l_dulo(1,n),irank,ulo_der,
405 : < ulos(1,1),dulos(1,1),uulon(1,1),dulon(1,1),
406 : < uloulopn(1,1,1),uuilon,duilon,ulouilopn,flo(n,:,:,:))
407 : c enddo
408 : c na = na + neq(n)
409 : 40 continue
410 : #else
411 0 : call judft_error("NOT implemented")
412 : #endif
413 :
414 :
415 0 : allocate(knorm(fullnkpts,num_bands))
416 0 : allocate(wann_acof(num_wann,0:lmd,natd,-cell1:cell1,-cell2:cell2,
417 0 : & -cell3:cell3))
418 0 : allocate(wann_bcof(num_wann,0:lmd,natd,-cell1:cell1,-cell2:cell2,
419 0 : & -cell3:cell3))
420 0 : allocate(wann_ccof(num_wann,-llod:llod,nlod,natd,-cell1:cell1,
421 0 : & -cell2:cell2,-cell3:cell3))
422 0 : allocate(wannint(-unigrid(4):unigrid(4),-unigrid(4):unigrid(4),
423 0 : , -unigrid(4):unigrid(4),num_wann))
424 0 : wann_acof(:,:,:,:,:,:)=0.0
425 0 : wann_bcof(:,:,:,:,:,:)=0.0
426 0 : wann_ccof(:,:,:,:,:,:,:)=0.0
427 0 : wannint(:,:,:,:)=0.0
428 :
429 0 : print*,"num_bands=",num_bands
430 0 : print*,"num_wann=",num_wann
431 0 : knorm(:,:)=0.0
432 :
433 : c******************************************************************
434 : c beginning of k-point loop,each may be a separate task
435 : c******************************************************************
436 0 : i_rec = 0 ; n_rank = 0
437 0 : do ikpt = 1,fullnkpts ! loop by k-points starts
438 :
439 0 : i_rec = i_rec + 1
440 0 : if (mod(i_rec-1,isize).eq.irank) then
441 0 : print*,"k-point=",ikpt
442 0 : kptibz=ikpt
443 0 : if(l_bzsym) kptibz=irreduc(ikpt)
444 0 : if(l_bzsym) oper=mapkoper(ikpt)
445 :
446 0 : if(have_disentangled) then
447 0 : inc_band(:)=lwindow(:,ikpt)
448 0 : num_inc=ndimwin(ikpt)
449 : end if
450 :
451 0 : allocate (we(neigd),eigg(neigd))
452 :
453 0 : zzMat%l_real = l_real
454 0 : zzMat%matsize1 = nbasfcn
455 0 : zzMat%matsize2 = neigd
456 0 : IF(l_real) THEN
457 0 : ALLOCATE (zzMat%data_r(zzMat%matsize1,zzMat%matsize2))
458 : ELSE
459 0 : ALLOCATE (zzMat%data_c(zzMat%matsize1,zzMat%matsize2))
460 : END IF
461 :
462 : call wann_read_eig(
463 : > eig_id,
464 : > ntypd,neigd,nvd,jspd,
465 : > irank,isize,kptibz,jspin,nbasfcn,
466 : > l_ss,l_noco,nrec,
467 : < nmat,nbands,eigg,zzMat,
468 0 : > .false.,1)
469 :
470 :
471 0 : zMat%l_real = zzMat%l_real
472 0 : zMat%matsize1 = zzMat%matsize1
473 0 : zMat%matsize2 = zzMat%matsize2
474 0 : IF (zzMat%l_real) THEN
475 0 : ALLOCATE (zMat%data_r(zMat%matsize1,zMat%matsize2))
476 0 : zMat%data_r = 0.0
477 : ELSE
478 0 : ALLOCATE (zMat%data_c(zMat%matsize1,zMat%matsize2))
479 0 : zMat%data_c = CMPLX(0.0,0.0)
480 : END IF
481 :
482 0 : nslibd = 0
483 :
484 : c...we work only within the energy window
485 :
486 0 : eig(:) = 0.
487 :
488 0 : print*,"bands used"
489 0 : do i = 1,nbands
490 : if ((eigg(i).ge.e1s .and. nslibd.lt.num_bands.and.l_bynumber)
491 : &.or.(eigg(i).ge.e1s.and.eigg(i).le.e2s.and.l_byenergy)
492 0 : &.or.(i.ge.band_min(jspin)
493 0 : & .and.i.le.band_max(jspin).and.l_byindex))then
494 0 : print*,i
495 0 : nslibd = nslibd + 1
496 0 : eig(nslibd) = eigg(i)
497 0 : we(nslibd) = we(i)
498 0 : IF(zzMat%l_real) THEN
499 0 : do j = 1,nv(jspin) + nlotot
500 0 : zMat%data_r(j,nslibd) = zzMat%data_r(j,i)
501 : end do
502 : ELSE
503 0 : do j = 1,nv(jspin) + nlotot
504 0 : zMat%data_c(j,nslibd) = zzMat%data_c(j,i)
505 : end do
506 : END IF
507 : endif
508 : enddo
509 :
510 : c***********************************************************
511 : c rotate the wavefunction
512 : c***********************************************************
513 0 : if (l_bzsym.and.oper.ne.1) then !rotate bkpt
514 0 : bkrot(:)=0.0
515 0 : do k=1,3
516 0 : bkrot(:)=bkrot(:)+mrot(k,:,oper)*bkpt(k)
517 : enddo
518 0 : bkpt(:)=bkrot(:)
519 :
520 0 : jloop:do j=1,nv(jspin)
521 : j1=mrot(1,1,oper)*k1(j,jspin)+
522 0 : + mrot(2,1,oper)*k2(j,jspin)+mrot(3,1,oper)*k3(j,jspin)
523 : j2=mrot(1,2,oper)*k1(j,jspin)+
524 0 : + mrot(2,2,oper)*k2(j,jspin)+mrot(3,2,oper)*k3(j,jspin)
525 : j3=mrot(1,3,oper)*k1(j,jspin)+
526 0 : + mrot(2,3,oper)*k2(j,jspin)+mrot(3,3,oper)*k3(j,jspin)
527 0 : k1(j,jspin)=j1
528 0 : k2(j,jspin)=j2
529 0 : k3(j,jspin)=j3
530 : enddo jloop
531 :
532 : endif
533 0 : print*,"bkpt=",bkpt(:)
534 : c******************************************************************
535 : c calculate a-, b-, and c-coefficients
536 : c******************************************************************
537 :
538 0 : noccbd = nslibd
539 :
540 0 : allocate ( acof(noccbd,0:lmd,natd),
541 0 : & bcof(noccbd,0:lmd,natd),
542 0 : & ccof(-llod:llod,noccbd,nlod,natd))
543 :
544 0 : acof(:,:,:) = cmplx(0.,0.) ; bcof(:,:,:) = cmplx(0.,0.)
545 0 : ccof(:,:,:,:) = cmplx(0.,0.)
546 :
547 : c...generation of the A,B,C coefficients in the spheres
548 : c...for the lapws and local orbitals, summed by the basis functions
549 :
550 : CALL abcof(input,atoms,sym,cell,lapw,noccbd,usdus,
551 0 : + noco,nococonv,jspin ,acof,bcof,ccof,zMat)
552 :
553 : call abcrot(
554 : > ntypd,natd,noccbd,lmaxd,lmd,llod,nlod,ntype,neq,
555 : > noccbd,lmax,nlo,llo,nop,ngopr,mrot,invsat,invsatnr,
556 : > bmat,
557 0 : X acof,bcof,ccof)
558 :
559 : c***************************************************************
560 : c calculate wannierfunctions' a-,b-, and c-coefficients
561 : c***************************************************************
562 0 : do i1=-cell1,cell1
563 0 : do i2=-cell2,cell2
564 0 : do i3=-cell3,cell3
565 : factor=exp(ImagUnit*tpi_const*
566 0 : * (i1*bkpt(1)+i2*bkpt(2)+i3*bkpt(3)))
567 0 : do i=1,num_wann
568 0 : do n=1,noccbd
569 : wann_acof(i,0:lmd,1:natd,i1,i2,i3)=
570 : = wann_acof(i,0:lmd,1:natd,i1,i2,i3)+
571 0 : + u_matrix(n,i,ikpt)*acof(n,0:lmd,1:natd)*factor
572 :
573 : wann_bcof(i,0:lmd,1:natd,i1,i2,i3)=
574 : = wann_bcof(i,0:lmd,1:natd,i1,i2,i3)+
575 0 : + u_matrix(n,i,ikpt)*bcof(n,0:lmd,1:natd)*factor
576 :
577 : wann_ccof(i,-llod:llod,1:nlod,1:natd,i1,i2,i3)=
578 : = wann_ccof(i,-llod:llod,1:nlod,1:natd,i1,i2,i3)+
579 0 : + u_matrix(n,i,ikpt)*ccof(-llod:llod,n,1:nlod,1:natd)*factor
580 : enddo
581 : enddo
582 : enddo
583 : enddo
584 : enddo
585 : c***************************************************************
586 : c calculate wannierfunctions' planewave-expansion
587 : c***************************************************************
588 0 : allocate(wannz(nv(jspin),num_wann))
589 0 : allocate(wannz2(nv(jspin),num_wann))
590 0 : wannz(:,:)=cmplx(0.0,0.0)
591 0 : do n=1,noccbd
592 0 : do m=1,num_wann
593 0 : IF(zMat%l_real) THEN
594 0 : do j=1, nv(jspin)
595 : wannz(j,m)=wannz(j,m)+zMat%data_r(j,n)*
596 0 : + u_matrix(n,m,ikpt)/sqrt(omtil)
597 : enddo
598 : ELSE
599 0 : do j=1, nv(jspin)
600 : wannz(j,m)=wannz(j,m)+zMat%data_c(j,n)*
601 0 : + u_matrix(n,m,ikpt)/sqrt(omtil)
602 : enddo
603 : END IF
604 :
605 : enddo
606 : enddo
607 0 : print*,"unigrid=",unigrid(:)
608 0 : do j=1,nv(jspin)
609 0 : do pw1=-unigrid(4),unigrid(4)
610 0 : do pw2=-unigrid(4),unigrid(4)
611 0 : do pw3=-unigrid(4),unigrid(4)
612 0 : denom=-pw1+unigrid(1)*(k1(j,jspin)+bkpt(1))
613 0 : wannz2(j,:)=wannz(j,:)
614 0 : if(abs(denom).gt.1e-5)then
615 0 : denom=denom*tpi_const/2
616 0 : factor=cmplx(cos(denom),sin(denom))
617 : wannz2(j,:)=wannz2(j,:)*(factor-conjg(factor))/
618 0 : / (ImagUnit*denom*2)
619 : endif
620 0 : denom=-pw2+unigrid(2)*(k2(j,jspin)+bkpt(2))
621 0 : if(abs(denom).gt.1e-5)then
622 0 : denom=denom*tpi_const/2
623 0 : factor=cmplx(cos(denom),sin(denom))
624 : wannz2(j,:)=wannz2(j,:)*(factor-conjg(factor))/
625 0 : / (ImagUnit*denom*2)
626 : endif
627 0 : denom=-pw3+unigrid(3)*(k3(j,jspin)+bkpt(3))
628 0 : if(abs(denom).gt.1e-5)then
629 0 : denom=denom*tpi_const/2
630 0 : factor=cmplx(cos(denom),sin(denom))
631 : wannz2(j,:)=wannz2(j,:)*(factor-conjg(factor))/
632 0 : / (ImagUnit*denom*2)
633 : endif
634 : wannint(pw1,pw2,pw3,:)=wannint(pw1,pw2,pw3,:)+
635 0 : + wannz2(j,:)
636 : enddo
637 : enddo
638 : enddo
639 : enddo
640 0 : deallocate(wannz,wannz2)
641 :
642 0 : deallocate ( acof,bcof,ccof,we,eigg )
643 :
644 0 : write (*,*) 'nslibd=',nslibd
645 :
646 : endif!processors
647 :
648 : enddo !loop over k-points
649 : c****************************************************************
650 : c end of k-loop
651 : c****************************************************************
652 :
653 : c****************************************************************
654 : c write radial wavefunctions to file
655 : c****************************************************************
656 0 : open(344,file=spin12(jspin)//'.lapw',form='unformatted')
657 0 : write(344)num_wann,cell1,cell2,cell3
658 0 : write(344)amat(1:3,1:3)
659 0 : write(344)ntype,jmtd,lmaxd,nlod,natd,lmd,llod
660 0 : write(344)pos(1:3,1:natd)
661 0 : write(344)neq(1:ntype)
662 0 : write(344)zatom(1:ntype)
663 0 : write(344)dx(1:ntype)
664 0 : write(344)rmt(1:ntype)
665 0 : write(344)jri(1:ntype)
666 0 : write(344)rmsh(1:jmtd,1:ntype)
667 0 : write(344)ff(1:ntype,1:jmtd,1:2,0:lmaxd)
668 0 : write(344)gg(1:ntype,1:jmtd,1:2,0:lmaxd)
669 0 : write(344)flo(1:ntype,1:jmtd,1:2,1:nlod)
670 : c****************************************************************
671 : c write a-,b-, and c-coefficients to file
672 : c****************************************************************
673 :
674 : write(344)wann_acof(1:num_wann,0:lmd,1:natd,-cell1:cell1,
675 0 : & -cell2:cell2,-cell3:cell3)
676 : write(344)wann_bcof(1:num_wann,0:lmd,1:natd,-cell1:cell1,
677 0 : & -cell2:cell2,-cell3:cell3)
678 : write(344)wann_ccof(1:num_wann,-llod:llod,1:nlod,1:natd,
679 0 : & -cell1:cell1,-cell2:cell2,-cell3:cell3)
680 : c****************************************************************
681 : c write planewave expansion to file
682 : c****************************************************************
683 0 : write(344)unigrid(1:4)
684 0 : write(344)wannint(:,:,:,:)
685 0 : close(344)
686 :
687 :
688 :
689 :
690 :
691 :
692 0 : deallocate(flo)
693 :
694 :
695 0 : if(have_disentangled)deallocate(lwindow,ndimwin,inc_band)
696 :
697 :
698 0 : deallocate(u_matrix)
699 :
700 : #ifdef CPP_MPI
701 0 : call MPI_BARRIER(mpi_communicatior,mpiierr(1))
702 : #endif
703 :
704 0 : nrec=nrec+nkpts
705 0 : 110 continue ! end of cycle by spins
706 :
707 :
708 0 : deallocate ( vr,nv,k1,k2,k3 )
709 0 : deallocate ( ff,gg)
710 :
711 : #ifdef CPP_MPI
712 0 : call MPI_BARRIER(mpi_communicatior,mpiierr(2))
713 : #endif
714 :
715 0 : call timestop("wannier_to_lapw")
716 0 : end subroutine wannier_to_lapw
717 : end module m_wannier_to_lapw
|
