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