Line data Source code
1 : !--------------------------------------------------------------------------------
2 : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
3 : ! This file is part of FLEUR and available as free software under the conditions
4 : ! of the MIT license as expressed in the LICENSE file in more detail.
5 : !--------------------------------------------------------------------------------
6 :
7 : MODULE m_wann_plot_um_dat
8 : use m_juDFT
9 : c******************************************************************
10 : c plot wannierfunctions directly within fleur
11 : c based on wann_plot
12 : c FF, September 2006
13 : c******************************************************************
14 : CONTAINS
15 0 : SUBROUTINE wann_plot_um_dat(
16 : > DIMENSION,stars,vacuum,atoms,sphhar,input,sym,mpi,
17 : > lapw,oneD,noco,cell,vTot,enpara,eig_id,l_real,
18 : > mpi_comm,sortrule,band_min,band_max,l_soc,
19 : > l_dulo,l_noco,l_ss,lmaxd,
20 : > ntypd,
21 : > neigd,natd,nop,nvd,jspd,nbasfcn,llod,nlod,ntype,
22 0 : > omtil,nlo,llo,lapw_l,invtab,mrot,ngopr,neq,lmax,
23 0 : > invsat,invsatnr,nkpt,taual,rmt,amat,bmat,bbmat,alph,
24 0 : > beta,qss,sk2,phi2,odi,ods,irank,isize,n3d,nmzxyd,nmzd,
25 0 : > jmtd,nlhd,nq3,nvac,invs,invs2,film,nlh,jri,ntypsd,
26 0 : > ntypsy,jspins,nkptd,dx,n2d,rmsh,e1s,e2s,ulo_der,
27 0 : > ustep,ig,k1d,k2d,k3d,rgphs,slice,kk,nnne,
28 0 : > z1,nv2d,nmzxy,nmz,delz,ig2,area,tau,zatom,nq2,nop2,
29 0 : > volint,symor,pos,ef,l_bzsym,
30 : > l_proj_plot,wan90version)
31 :
32 : use m_wann_rw_eig
33 : use m_abcof
34 : use m_wann_2dvacabcof
35 : use m_wann_1dvacabcof
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_read_umatrix
44 : use m_wann_kptsrotate
45 : use m_wann_plot_vac
46 : use m_wann_plot_od_vac
47 : use m_wann_abinv
48 :
49 :
50 : implicit none
51 :
52 : #include "cpp_double.h"
53 : #ifdef CPP_MPI
54 : include 'mpif.h'
55 : integer mpiierr(3)
56 : integer cpu_index
57 : integer stt(MPI_STATUS_SIZE)
58 : #endif
59 :
60 : TYPE(t_dimension),INTENT(IN) :: DIMENSION
61 : TYPE(t_stars),INTENT(IN) :: stars
62 : TYPE(t_vacuum),INTENT(IN) :: vacuum
63 : TYPE(t_atoms),INTENT(IN) :: atoms
64 : TYPE(t_sphhar),INTENT(IN) :: sphhar
65 : TYPE(t_input),INTENT(IN) :: input
66 : TYPE(t_sym),INTENT(IN) :: sym
67 : TYPE(t_mpi),INTENT(IN) :: mpi
68 : TYPE(t_lapw),INTENT(IN) :: lapw
69 : TYPE(t_oneD),INTENT(IN) :: oneD
70 : TYPE(t_noco),INTENT(IN) :: noco
71 : TYPE(t_cell),INTENT(IN) :: cell
72 : TYPE(t_potden),INTENT(IN) :: vTot
73 : TYPE(t_enpara),INTENT(IN) :: enpara
74 :
75 : integer, intent (in) :: band_min(2),band_max(2),mpi_comm,eig_id
76 : logical, intent (in) :: l_soc,l_real
77 : logical, intent (in) :: invs,invs2,film,slice,symor
78 : integer, intent (in) :: lmaxd,ntypd,neigd,nkptd,kk,nnne
79 : integer, intent (in) :: natd,nop,nvd,jspd,nbasfcn,nq2,nop2
80 : integer, intent (in) :: llod,nlod,ntype,n3d,n2d
81 : integer, intent (in) :: nmzxyd,nmzd,jmtd,nlhd,nq3,nvac
82 : integer, intent (in) :: ntypsd,jspins,k1d,k2d,k3d
83 : real, intent (in) :: omtil,e1s,e2s,delz,area,z1,volint
84 : integer, intent (in) :: ig(-k1d:k1d,-k2d:k2d,-k3d:k3d)
85 : complex, intent (in) :: rgphs(-k1d:k1d,-k2d:k2d,-k3d:k3d)
86 : integer, intent (in) :: nlh(ntypsd),jri(ntypd),ntypsy(natd)
87 : integer, intent (in) :: nlo(ntypd),llo(nlod,ntypd),lapw_l(ntypd)
88 : integer, intent (in) :: invtab(nop),mrot(3,3,nop),ngopr(natd)
89 : integer, intent (in) :: neq(ntypd),lmax(ntypd)
90 : integer, intent (in) :: wan90version
91 : integer, intent (in) :: invsat(natd),invsatnr(natd),nkpt
92 : integer, intent (in) :: irank,isize,nv2d,nmzxy,nmz
93 : integer, intent (in) :: ulo_der(nlod,ntypd),ig2(n3d)
94 : real, intent (in) :: taual(3,natd),rmt(ntypd),dx(ntypd)
95 : real, intent (in) :: amat(3,3),bmat(3,3),bbmat(3,3)
96 : real, intent (in) :: rmsh(jmtd,ntypd),tau(3,nop),zatom(ntype)
97 : real, intent (in) :: alph(ntypd),beta(ntypd),qss(3)
98 : real, intent (in) :: pos(3,natd),ef
99 : complex, intent (in) :: ustep(n3d)
100 : logical, intent (in) :: l_dulo(nlod,ntypd),l_noco,l_ss,l_bzsym
101 : logical,intent(in)::l_proj_plot
102 : integer,intent(in)::sortrule
103 : c-odim
104 : real, intent (in) :: sk2(n2d),phi2(n2d)
105 : type (od_inp), intent (in) :: odi
106 : type (od_sym), intent (in) :: ods
107 : c+odim
108 : logical l_spreadcal
109 0 : complex, allocatable::spreads(:,:)
110 0 : real,allocatable:: centers(:,:)
111 : cccccccccccccccccc local variables cccccccccccccccccccc
112 : integer lmd,nlotot,n,nmat,nw,ispin,iter,ikpt,ilo
113 : integer noccbd,nn,nkpts,i,jspin,j,l,i_rec,m,nwf,nwfp
114 : integer jsp_start,jsp_end,nrec,nrec1,nbands
115 : integer nodeu,noded,n_size,na,n_rank,nbnd,nkbnd,idum,jdum,kdum
116 : integer i1,i2,i3,in,ikpt_k,lda
117 0 : integer n_bands(0:neigd),nslibd
118 : character*8 dop,iop
119 : real bkpt(3),sfp,tpi,wronk,wk,wk_b,phase
120 0 : real eig(neigd),cp_time(9)
121 : logical l_p0,l_bkpts,l_proj,l_amn,l_mmn,l_eig,wann,wann_plott
122 : !!! energy window boundaries
123 0 : integer, allocatable :: nv(:)
124 0 : integer, allocatable :: k1(:,:),k2(:,:),k3(:,:)
125 0 : real, allocatable :: we(:)
126 :
127 0 : real, allocatable :: eigg(:)
128 : real kpoints(nkptd)
129 : !!! a and b coeff. constructed for each k-point
130 0 : complex, allocatable :: acof(:,:,:)
131 0 : complex, allocatable :: bcof(:,:,:)
132 0 : complex, allocatable :: 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(:,:,:),vz(:,:,:)
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 :
145 0 : real :: uuilon(nlod,ntypd),duilon(nlod,ntypd)
146 0 : real :: ulouilopn(nlod,nlod,ntypd)
147 : !!! energy parameters
148 : character(len=3) :: spin12(2)
149 : data spin12/'WF1' , 'WF2'/
150 : character(len=2)::spinspin12(2)
151 : data spinspin12/'.1','.2'/
152 0 : complex,allocatable::wannierfunc(:,:)
153 0 : complex,allocatable::wannierfunc_temp(:,:)
154 : character(len=33):: header
155 : integer :: num_nnmax
156 : integer :: posi
157 : complex,allocatable::u_matrix_tmp(:,:,:)
158 0 : complex,allocatable::u_matrix(:,:,:)
159 : real :: tmp_omi
160 : integer :: kpt,oper
161 : real :: poinint(3)
162 : real :: phas,tmax
163 : real :: bkrot(3)
164 : integer :: j1,j2,j3
165 : logical :: um_format
166 : logical :: have_disentangled,l_chk,l_umdat
167 0 : integer,allocatable :: ndimwin(:)
168 0 : logical,allocatable :: lwindow(:,:)
169 : integer :: chk_unit,nkp,ntmp,ierr,fullnkpts
170 0 : integer,allocatable::irreduc(:),mapkoper(:)
171 : character(len=20)::checkpoint
172 : real :: tmp_latt(3,3)
173 : real,allocatable:: tmp_kpt_latt(:,:)
174 : real omega_invariant
175 : complex,allocatable::u_matrix_opt(:,:,:)
176 : logical l_file
177 0 : logical,allocatable::inc_band(:)
178 : integer :: num_inc,counter,kptibz,wannierspin
179 : logical :: l_byindex, l_byenergy, l_bynumber
180 : integer :: num_wann,num_bands,kpun,jspin2,jspin3
181 : complex :: d_wgn(-3:3,-3:3,3,nop)
182 : integer :: pos1,pos2,ato,loc,invop
183 : c ..basis wavefunctions in the vacuum
184 0 : complex, allocatable :: ac(:,:,:),bc(:,:,:)
185 0 : complex, allocatable :: ac_1(:,:,:),bc_1(:,:,:)
186 : real, allocatable :: dt(:),dte(:)
187 : real, allocatable :: t(:),te(:),tei(:)
188 0 : real, allocatable :: u(:,:,:),ue(:,:,:)
189 0 : real, allocatable :: u_1(:,:,:),ue_1(:,:,:)
190 : real :: vz0(2)
191 : integer :: ik,nv2,ivac,jvac,symvac,symvacvac
192 : real :: evacp,sign,arg
193 : complex :: c_1
194 : integer :: kvac1(nv2d),kvac2(nv2d),map2(nvd)
195 : real :: fas,zks
196 : integer :: mesh
197 : integer :: n2
198 : real :: v(3),scale,ev
199 : complex :: av,bv
200 : real :: volume
201 : c external dotirp !module now
202 : c real dotirp
203 : REAL :: s,const
204 : COMPLEX :: xdnout,factor
205 : INTEGER :: ii3,ix,iy,iz,nplo,nbn
206 : INTEGER :: nbmin,nbmax
207 : INTEGER :: nplot,nq,nt,jm,ii1,ii2
208 : LOGICAL :: twodim
209 0 : real,allocatable::knorm(:,:)
210 0 : real,allocatable::wfnorm(:)
211 : REAL :: pt(3),vec1(3),vec2(3),vec3(3),zero(3)
212 : INTEGER :: grid(3),k,addnoco
213 0 : integer,allocatable :: shiftkpt(:,:)
214 : LOGICAL :: cartesian,xsf
215 : REAL :: rhocc(jmtd),realpart,imagpart
216 : REAL :: point(3),post(3,natd)
217 : CHARACTER(len=30):: filename
218 : CHARACTER(len=20):: name1,name2,name3
219 : CHARACTER(len=10):: vandername
220 : NAMELIST /plot/twodim,cartesian,vec1,vec2,vec3,grid,zero,filename
221 : complex :: nsfactor
222 0 : integer :: ngopr1(natd)
223 :
224 0 : TYPE(t_usdus) :: usdus
225 0 : TYPE(t_mat) :: zzMat, zMat
226 :
227 0 : um_format=.false.
228 0 : l_byindex=.false.
229 0 : l_byenergy=.false.
230 0 : l_bynumber=.false.
231 0 : if(sortrule==1)l_byindex=.true.
232 0 : if(sortrule==2)l_bynumber=.true.
233 0 : if(sortrule==3)l_byenergy=.true.
234 0 : ngopr1(:)=1
235 :
236 :
237 :
238 : c read in plot_inp
239 :
240 0 : INQUIRE(file ="plot_inp",exist= twodim)
241 0 : IF (.NOT.twodim) THEN !no input file exists, create a template and
242 : !exit
243 0 : OPEN(20,file ="plot_inp")
244 0 : WRITE(20,'(i2,a5,l1)') 1,",xsf=",.false.
245 : c WRITE(20,*) "&PLOT twodim=t,cartesian=t"
246 : c WRITE(20,*) " vec1(1)=10.0 vec2(2)=10.0"
247 : c WRITE(20,*) " filename='plot1' /"
248 0 : WRITE(20,*) "&PLOT twodim=f,cartesian=f"
249 0 : WRITE(20,*) " vec1(1)=1.0 vec1(2)=0.0 vec1(3)=0.0 "
250 0 : WRITE(20,*) " vec2(1)=0.0 vec2(2)=1.0 vec2(3)=0.0 "
251 0 : WRITE(20,*) " vec3(1)=0.0 vec3(2)=0.0 vec3(3)=1.0 "
252 0 : WRITE(20,*) " grid(1)=30 grid(2)=30 grid(3)=30 "
253 0 : WRITE(20,*) " zero(1)=0.0 zero(2)=0.0 zero(3)=0.0 "
254 0 : WRITE(20,*) " filename ='plot2' /"
255 0 : CLOSE(20)
256 0 : WRITE(*,*) "No plot_inp file found. Created a template"
257 : CALL juDFT_error("Missing input for plot; modify plot_inp"
258 0 : + ,calledby ="wann_plot_um_dat")
259 : ENDIF
260 :
261 0 : OPEN (18,file='plot_inp')
262 0 : READ(18,'(i2,5x,l1)') nplot,xsf
263 : ! If xsf is specified we create an input file for xcrysden
264 0 : IF (nplot.ge.2)
265 : & CALL juDFT_error
266 : + ("plots one by one, please, this is not charge density"
267 0 : + ,calledby ="wann_plot_um_dat")
268 0 : twodim = .TRUE.;cartesian=.TRUE.;grid=(/100,100,100/)
269 0 : vec1 = (/0.,0.,0./);vec2=(/0.,0.,0./);vec3=(/0.,0.,0./)
270 0 : zero = (/0.,0.,0./);filename="default"
271 0 : READ(18,plot)
272 0 : IF (twodim.AND.ANY(grid(1:2)<1))
273 : + CALL juDFT_error("Illegal grid size in plot",calledby
274 0 : + ="wann_plot_um_dat")
275 0 : IF (.NOT.twodim.AND.ANY(grid<1))
276 : + CALL juDFT_error("Illegal grid size in plot",calledby
277 0 : + ="wann_plot_um_dat")
278 0 : IF (twodim) grid(3) = 1
279 : !calculate cartesian coordinates if needed
280 0 : IF (.NOT.cartesian) THEN
281 0 : vec1=matmul(amat,vec1)
282 0 : vec2=matmul(amat,vec2)
283 0 : vec3=matmul(amat,vec3)
284 0 : zero=matmul(amat,zero)
285 : ENDIF
286 0 : Close(18)
287 :
288 : !calculate volume
289 : volume = vec1(1)*vec2(2)*vec3(3) + vec2(1)*vec3(2)*vec1(3) +
290 : & vec3(1)*vec1(2)*vec2(3) - vec1(3)*vec2(2)*vec3(1) -
291 0 : & vec2(3)*vec3(2)*vec1(1) - vec3(3)*vec1(2)*vec2(1)
292 :
293 :
294 :
295 0 : sfp = 2* sqrt( pimach() )
296 0 : tpi = 2* pimach()
297 0 : lmd = lmaxd*(lmaxd+2)
298 0 : nkpts = nkpt
299 :
300 0 : nrec = 0
301 0 : nlotot = 0
302 0 : do n = 1, ntype
303 0 : do l = 1,nlo(n)
304 0 : nlotot = nlotot + neq(n) * ( 2*llo(l,n) + 1 )
305 : enddo
306 : enddo
307 :
308 : cccccccccccccccc initialize the potential cccccccccccc
309 :
310 0 : allocate (vz(nmzd,2,jspd))
311 0 : allocate (vr(jmtd,ntypd,jspd))
312 :
313 0 : vz = vTot%vacz
314 :
315 0 : do jspin = 1,jspins
316 0 : do n = 1, ntype
317 0 : do j = 1,jri(n)
318 0 : vr(j,n,jspin) = vTot%mt(j,0,n,jspin)
319 : enddo
320 : enddo
321 : enddo
322 :
323 : cccccccccccccccc end of the potential part ccccccccccc
324 0 : wannierspin=jspins
325 0 : if(l_soc) wannierspin=2
326 :
327 0 : allocate ( nv(wannierspin) )
328 0 : allocate ( k1(nvd,jspd),k2(nvd,jspd),k3(nvd,jspd) )
329 0 : allocate ( ff(ntypd,jmtd,2,0:lmaxd) )
330 0 : allocate ( gg(ntypd,jmtd,2,0:lmaxd) )
331 0 : allocate ( usdus%us(0:lmaxd,ntypd,jspins) )
332 0 : allocate ( usdus%uds(0:lmaxd,ntypd,jspins) )
333 0 : allocate ( usdus%dus(0:lmaxd,ntypd,jspins) )
334 0 : allocate ( usdus%duds(0:lmaxd,ntypd,jspins) )
335 0 : allocate ( usdus%ddn(0:lmaxd,ntypd,jspins) )
336 0 : allocate ( usdus%ulos(nlod,ntypd,jspins) )
337 0 : allocate ( usdus%dulos(nlod,ntypd,jspins) )
338 0 : allocate ( usdus%uulon(nlod,ntypd,jspins) )
339 0 : allocate ( usdus%dulon(nlod,ntypd,jspins) )
340 0 : allocate ( usdus%uloulopn(nlod,nlod,ntypd,jspins) )
341 :
342 0 : do 110 jspin=1,wannierspin ! cycle by spins
343 0 : print*,"spin=",jspin
344 0 : jspin2=jspin
345 0 : if(l_soc.and.jspins.eq.1)jspin2=1
346 0 : jspin3=jspin
347 0 : if(l_soc)jspin3=1
348 0 : jsp_start = jspin ; jsp_end = jspin
349 :
350 0 : addnoco=0
351 0 : if(l_noco.and.(jspin.eq.2))then
352 0 : addnoco=nv(1)+nlotot
353 : endif
354 :
355 : c*******************************************************
356 : c get num_bands and num_wann from WF1.amn (WF2.amn)
357 : c*******************************************************
358 0 : l_file=.false.
359 0 : inquire(file=spin12(jspin3)//'.amn',exist=l_file)
360 0 : open(355,file=spin12(jspin3)//'.amn')
361 0 : read(355,*)
362 0 : read(355,*)num_bands,kpun,num_wann
363 0 : close(355)
364 0 : if(l_byindex.and..not.((1+band_max(jspin)-
365 : & band_min(jspin)).eq.num_bands))
366 : & CALL juDFT_error("1+band_max-band_min /= num_bands",
367 0 : + calledby ="wann_plot_um_dat")
368 :
369 : c**************************************************************
370 : ! for bzsym = .true.: determine mapping between kpts and w90kpts
371 : c**************************************************************
372 0 : if (l_bzsym) then
373 0 : l_file=.false.
374 0 : inquire(file='w90kpts',exist=l_file)
375 0 : IF(.NOT.l_file) CALL juDFT_error
376 : + ("w90kpts not found, needed if bzsym",calledby
377 0 : + ="wann_plot_um_dat")
378 0 : open(412,file='w90kpts',form='formatted')
379 0 : read(412,*)fullnkpts
380 0 : close(412)
381 0 : print*,"fullnkpts=",fullnkpts
382 0 : IF(fullnkpts<=nkpts) CALL juDFT_error("fullnkpts.le.nkpts"
383 0 : + ,calledby ="wann_plot_um_dat")
384 0 : allocate(irreduc(fullnkpts),mapkoper(fullnkpts))
385 0 : allocate(shiftkpt(3,fullnkpts))
386 0 : l_file=.false.
387 0 : inquire(file='kptsmap',exist=l_file)
388 0 : IF(.NOT.l_file) CALL juDFT_error
389 : + ("kptsmap not found, needed if bzsym",calledby
390 0 : + ="wann_plot_um_dat")
391 0 : open(713,file='kptsmap')
392 0 : do i=1,fullnkpts
393 0 : read(713,*)kpt,irreduc(i),mapkoper(i),shiftkpt(:,i)
394 0 : IF(kpt/=i) CALL juDFT_error("kpt.ne.i",calledby
395 0 : + ="wann_plot_um_dat")
396 0 : print*,i,irreduc(i),mapkoper(i)
397 : enddo
398 0 : close(713)
399 0 : IF(MAXVAL(irreduc(:))/=nkpts) CALL juDFT_error
400 0 : + ("max(irreduc(:))/=nkpts",calledby ="wann_plot_um_dat")
401 : else
402 0 : fullnkpts=nkpts
403 : endif
404 :
405 0 : IF(kpun/=fullnkpts) CALL juDFT_error
406 : + ("mismatch in kpun and fullnkpts",calledby
407 0 : + ="wann_plot_um_dat")
408 :
409 0 : if(.not.l_proj_plot)then
410 : c**************************************************************
411 : c read in chk
412 : c*************************************************************
413 0 : allocate( u_matrix(num_bands,num_wann,fullnkpts) )
414 0 : allocate( lwindow(num_bands,fullnkpts) )
415 0 : allocate( ndimwin(fullnkpts) )
416 : call wann_read_umatrix(
417 : > fullnkpts,num_wann,num_bands,
418 : > um_format,jspin,wan90version,
419 : < have_disentangled,
420 : < lwindow,ndimwin,
421 0 : < u_matrix)
422 : else
423 : c**************************************************************
424 : c read WF1.umn (WF2.umn) (if projmethod)
425 : c**************************************************************
426 0 : have_disentangled=.false.
427 0 : l_file=.false.
428 0 : inquire(file=spin12(jspin)//'.umn',exist=l_file)
429 0 : IF(.NOT.l_file) CALL juDFT_error("no umn file foun
430 0 : +d",calledby ="wann_plot_um_dat")
431 0 : open(419,file=spin12(jspin)//'.umn')
432 0 : read(419,*) !num_wann,num_bands
433 0 : allocate(u_matrix(num_bands,num_wann,fullnkpts))
434 0 : do ikpt=1,fullnkpts
435 0 : do j=1,num_wann
436 0 : do i=1,num_bands
437 0 : read(419,*)idum,jdum,kdum,realpart,imagpart
438 0 : u_matrix(i,j,ikpt)=cmplx(realpart,imagpart)
439 : enddo
440 : enddo
441 : enddo
442 0 : close(419)
443 : endif
444 :
445 : c if(um_format)then
446 : c open(419,file='umatrix_formatted')
447 : c do ikpt=1,fullnkpts
448 : c do j=1,num_wann
449 : c do i=1,num_bands
450 : c write(419,*)u_matrix(i,j,ikpt)
451 : c enddo
452 : c enddo
453 : c enddo
454 : c close(419)
455 : c endif
456 :
457 :
458 : ***********************************************************
459 : ***********************************************************
460 :
461 0 : print*,"num_wann=",num_wann
462 0 : print*,"num_bands=",num_bands
463 : allocate(wannierfunc(num_wann,
464 0 : & (grid(1))*(grid(2))*(grid(3))))
465 :
466 :
467 :
468 0 : wannierfunc(:,:)=cmplx(0.0,0.0)
469 :
470 :
471 : cccccccccccc read in the eigenvalues and vectors cccccc
472 :
473 0 : l_p0 = .false.
474 0 : if (irank.eq.0) l_p0 = .true.
475 :
476 : call cdn_read0(eig_id,irank,isize,jspin,wannierspin,l_noco,
477 0 : < n_bands,n_size)
478 :
479 :
480 0 : allocate ( flo(ntypd,jmtd,2,nlod) )
481 :
482 0 : na = 1
483 0 : do 40 n = 1,ntype
484 0 : do 30 l = 0,lmax(n)
485 : c...compute the l-dependent, k-independent radial MT- basis functions
486 : call radfun(
487 : > l,n,jspin,enpara%el0(l,n,jspin),vr(1,n,jspin2),atoms,
488 : < ff(n,:,:,l),gg(n,:,:,l),usdus,
489 0 : < nodeu,noded,wronk)
490 0 : 30 continue
491 : c...and the local orbital radial functions
492 : c do ilo = 1, nlo(n)
493 : call radflo(
494 : > atoms,n,jspin,enpara%ello0(:,:,jspin),vr(1,n,jspin2),
495 : > ff(n,1:,1:,0:),gg(n,1:,1:,0:),mpi,
496 0 : < usdus,uuilon,duilon,ulouilopn,flo(n,:,:,:))
497 : c enddo
498 : c na = na + neq(n)
499 0 : 40 continue
500 0 : i_rec = 0 ; n_rank = 0
501 :
502 : c******************************************************************
503 : c beginning of k-point loop,each may be a separate task
504 : c******************************************************************
505 :
506 0 : allocate(knorm(fullnkpts,num_bands))
507 0 : print*,"num_bands=",num_bands
508 0 : knorm(:,:)=0.0
509 :
510 0 : do ikpt = 1,fullnkpts ! loop by k-points starts
511 :
512 0 : i_rec = i_rec + 1
513 0 : if (mod(i_rec-1,isize).eq.irank) then
514 0 : print*,"k-point=",ikpt
515 0 : kptibz=ikpt
516 0 : if(l_bzsym) kptibz=irreduc(ikpt)
517 0 : if(l_bzsym) oper=mapkoper(ikpt)
518 :
519 0 : if(have_disentangled) then
520 0 : if(.not.allocated(inc_band))
521 0 : > allocate(inc_band(size(lwindow,1)))
522 0 : inc_band(:)=lwindow(:,ikpt)
523 0 : num_inc=ndimwin(ikpt)
524 : end if
525 :
526 0 : allocate (we(neigd),eigg(neigd))
527 :
528 0 : zzMat%l_real = l_real
529 0 : zzMat%matsize1 = nbasfcn
530 0 : zzMat%matsize2 = neigd
531 0 : IF(l_real) THEN
532 0 : ALLOCATE (zzMat%data_r(zzMat%matsize1,zzMat%matsize2))
533 : ELSE
534 0 : ALLOCATE (zzMat%data_c(zzMat%matsize1,zzMat%matsize2))
535 : END IF
536 0 : CALL judft_error("TODO:adjust in wann_plot_um_dat")
537 : ! call wann_read_eig(eig_id, lmaxd,ntypd,nlod,neigd,nvd,wannierspin,
538 : ! > irank,isize,kptibz,jspin,nbasfcn,nlotot, l_ss,l_noco,nrec,
539 : ! < nmat,nv, k1,k2,k3,bkpt,wk,nbands,eigg,zzMat, .false.,1)
540 :
541 0 : zMat%l_real = zzMat%l_real
542 0 : zMat%matsize1 = zzMat%matsize1
543 0 : zMat%matsize2 = zzMat%matsize2
544 0 : IF (zzMat%l_real) THEN
545 0 : ALLOCATE (zMat%data_r(zMat%matsize1,zMat%matsize2))
546 0 : zMat%data_r = 0.0
547 : ELSE
548 0 : ALLOCATE (zMat%data_c(zMat%matsize1,zMat%matsize2))
549 0 : zMat%data_c = CMPLX(0.0,0.0)
550 : END IF
551 :
552 :
553 0 : nslibd = 0
554 :
555 : c...we work only within the energy window
556 :
557 0 : eig(:) = 0.
558 :
559 0 : print*,"bands used"
560 0 : do i = 1,nbands
561 : if ((eigg(i).ge.e1s .and. nslibd.lt.num_bands.and.l_bynumber)
562 : &.or.(eigg(i).ge.e1s.and.eigg(i).le.e2s.and.l_byenergy)
563 0 : &.or.(i.ge.band_min(jspin).and.i.le.band_max(jspin)
564 0 : &.and.l_byindex))then
565 0 : print*,i
566 0 : nslibd = nslibd + 1
567 0 : eig(nslibd) = eigg(i)
568 0 : we(nslibd) = we(i)
569 0 : IF(zzMat%l_real) THEN
570 0 : do j = 1,nv(jspin) + nlotot
571 0 : zMat%data_r(j,nslibd) = zzMat%data_r(j,i)
572 : end do
573 : ELSE
574 0 : do j = 1,nv(jspin) + nlotot
575 0 : zMat%data_c(j,nslibd) = zzMat%data_c(j,i)
576 : end do
577 : END IF
578 : endif
579 : enddo
580 : c***********************************************************
581 : c rotate the wavefunction
582 : c***********************************************************
583 :
584 0 : if (l_bzsym.and.oper.ne.1) then !rotate bkpt
585 : call wann_kptsrotate(
586 : > natd,nlod,llod,
587 : > ntypd,nlo,llo,invsat,
588 : > l_noco,l_soc,
589 : > ntype,neq,nlotot,
590 : > jspin,
591 : > oper,nop,mrot,nvd,nv,
592 : > shiftkpt(:,ikpt),
593 : > tau,
594 : x bkpt,k1(:,:),k2(:,:),k3(:,:),
595 0 : x zMat,nsfactor)
596 : endif
597 :
598 0 : print*,"bkpt=",bkpt(:)
599 : c******************************************************************
600 :
601 :
602 0 : noccbd = nslibd
603 :
604 0 : allocate ( acof(noccbd,0:lmd,natd),
605 0 : & bcof(noccbd,0:lmd,natd),
606 0 : & ccof(-llod:llod,noccbd,nlod,natd))
607 :
608 0 : acof(:,:,:) = cmplx(0.,0.) ; bcof(:,:,:) = cmplx(0.,0.)
609 0 : ccof(:,:,:,:) = cmplx(0.,0.)
610 :
611 : c...generation of the A,B,C coefficients in the spheres
612 : c...for the lapws and local orbitals, summed by the basis functions
613 0 : ngopr1(:)=1
614 :
615 : CALL abcof(input,atoms,sym,cell,lapw,noccbd,usdus,
616 0 : > noco,jspin,oneD,acof,bcof,ccof,zMat)
617 :
618 : call wann_abinv(atoms,
619 0 : X acof,bcof,ccof)
620 :
621 :
622 : c***********************************************************************
623 : c make preparations for plotting in vacuum
624 : c***********************************************************************
625 0 : if (film.and. .not.odi%d1)then
626 0 : allocate ( ac(nv2d,nslibd,2),bc(nv2d,nslibd,2),
627 0 : + u(nmzd,nv2d,nvac),ue(nmzd,nv2d,nvac))
628 : call wann_2dvacabcof(
629 : > nv2d,nslibd,nvac,nmzd,nmz,omtil,vz(:,:,jspin2),
630 : > nv(jspin),bkpt,z1,
631 : > nvd,k1(:,jspin),k2(:,jspin),k3(:,jspin),enpara%evac0(:,jspin),
632 : > bbmat,delz,bmat,nbasfcn,
633 : > neigd,zMat,
634 0 : < ac,bc,u,ue,addnoco,l_ss,qss,jspin)
635 : endif !preparations for vacuum
636 0 : if (odi%d1)then
637 0 : allocate ( ac_1(nv2d,-odi%mb:odi%mb,nslibd),
638 0 : & bc_1(nv2d,-odi%mb:odi%mb,nslibd),
639 0 : & u_1(nmzd,nv2d,-odi%mb:odi%mb),
640 0 : & ue_1(nmzd,nv2d,-odi%mb:odi%mb) )
641 : call wann_1dvacabcof(
642 : > DIMENSION,oneD,vacuum,stars,cell,
643 : > nv2d,nslibd,nmzd,nmz,omtil,vz(:,:,jspin2),
644 : > nv(jspin),bkpt,z1,odi,ods,
645 : > nvd,k1(:,jspin),k2(:,jspin),k3(:,jspin),enpara%evac0(:,jspin),
646 : > bbmat,delz,bmat,nbasfcn,neigd,zMat,
647 : > n2d,n3d,ig,nmzxy,nmzxyd,ig2,sk2,phi2,k1d,k2d,k3d,
648 0 : < ac_1,bc_1,u_1,ue_1,addnoco,l_ss,qss,jspin)
649 : endif !preparations for vacuum
650 : c**************************************************************************
651 : c**************************************************************************
652 :
653 0 : nbmin=1
654 0 : nbmax=nslibd
655 0 : counter=1
656 :
657 0 : band:DO nbn = nbmin,nbmax
658 0 : if(have_disentangled) then
659 0 : if(counter>num_inc) exit
660 0 : if(.not.inc_band(nbn))cycle band
661 : endif
662 :
663 :
664 0 : DO iz = 0,grid(3)-1
665 0 : DO iy = 0,grid(2)-1
666 0 : xloop:DO ix = 0,grid(1)-1
667 0 : posi=ix+1+iy*(grid(1))+iz*(grid(1))*(grid(2))
668 : point = zero+vec1*(ix+0.0)/(grid(1)-1)+vec2*(iy+0.0)
669 0 : $ /(grid(2)-1)
670 0 : IF (.NOT.twodim) point = point+vec3*(iz+0.0)/(grid(3)-1)
671 0 : poinint=matmul(bmat,point)/tpi_const
672 : phas=tpi*(bkpt(1)*poinint(1)
673 0 : + +bkpt(2)*poinint(2)+bkpt(3)*poinint(3))
674 0 : factor=cmplx(cos(phas),sin(phas))
675 : !Check if the point is in MT-sphere
676 0 : ii1 = 3
677 0 : ii2 = 3
678 0 : ii3 = 3
679 0 : IF (film .AND. .NOT.odi%d1) ii3 = 0
680 0 : IF (odi%d1) THEN
681 0 : ii1 = 0 ; ii2 = 0
682 : END IF
683 0 : DO i1 = -ii1,ii1
684 0 : DO i2 = -ii2,ii2
685 0 : DO i3 = -ii3,ii3
686 0 : pt = point+MATMUL(amat,(/i1,i2,i3/))
687 0 : na = 0
688 0 : DO nt = 1,ntype
689 0 : DO nq = 1,neq(nt)
690 0 : na = na + 1
691 0 : s = SQRT(dot_PRODUCT(pos(:,na)-pt,pos(:,na)-pt))
692 0 : IF (s<rmsh(jri(nt),nt)) THEN
693 : CALL wann_real(
694 : > pt,nt,na,0,1,bkpt,nbn,
695 : > n3d,nmzxyd,n2d,ntypsd,lmaxd,jmtd,
696 : > natd,ntypd,nmzd,nop,nop2,mrot,tau,invtab,
697 : > nq3,nvac,invs,z1,delz,nmz,nmzxy,nq2,
698 : > lmax,rmsh,jri,pos,ngopr,ntypsy,nvd,
699 : > omtil,amat,bmat,odi,ods,nlod,llod,nlo,llo,
700 : > ff,gg,flo,acof(nbn,:,:),bcof(nbn,:,:),
701 : > ccof(:,nbn,:,:),zMat,
702 : > nv(jspin),k1(:,jspin),k2(:,jspin),k3(:,jspin),
703 : > lmd,nbasfcn,l_ss,qss,jspin,0,
704 0 : < xdnout)
705 : wannierfunc(:,posi)=
706 0 : = wannierfunc(:,posi)+xdnout*u_matrix(counter,:,ikpt)*factor
707 0 : knorm(ikpt,nbn)=knorm(ikpt,nbn)+(abs(xdnout))**2
708 0 : CYCLE xloop
709 : ENDIF
710 : ENDDO
711 : ENDDO !nt
712 : ENDDO
713 : ENDDO
714 : ENDDO !i1
715 : !Check for point in vacuum
716 0 : IF (film.AND..NOT.odi%d1.AND.ABS(point(3))>=z1) THEN
717 0 : ivac=1
718 0 : if (point(3).lt. 0.0)ivac=2
719 0 : jvac=ivac
720 0 : if(nvac==1)jvac=1
721 : call wann_plot_vac(point,z1,nmzd,nv2d,n3d,nvac,
722 : > nmz,delz,bmat,bbmat,enpara%evac0(:,jspin),bkpt,vz,jspin,
723 : > k1(:,jspin),k2(:,jspin),k3(:,jspin),nvd,
724 : > nbasfcn,neigd,nv(jspin),omtil,nslibd,
725 : > ac(:,nbn,ivac),
726 : & bc(:,nbn,ivac),
727 0 : & u(:,:,jvac),ue(:,:,jvac),xdnout)
728 :
729 : wannierfunc(:,posi)=
730 : = wannierfunc(:,posi)+
731 0 : + xdnout*u_matrix(counter,:,ikpt)*factor
732 : CYCLE xloop
733 : END IF
734 :
735 0 : IF (odi%d1) THEN
736 0 : IF (SQRT((pt(1))**2 + (pt(2))**2)>=z1) THEN
737 : call wann_plot_od_vac(point,z1,nmzd,nv2d,odi,
738 : > nmz,delz,bmat,bkpt,nvd,nv(jspin),omtil,k3(:,jspin),
739 : > ac_1(:,:,nbn),bc_1(:,:,nbn),u_1,ue_1,
740 0 : < xdnout)
741 : wannierfunc(:,posi)=
742 0 : =wannierfunc(:,posi)+xdnout*u_matrix(counter,:,ikpt)*factor
743 : CYCLE xloop
744 : END IF
745 : END IF
746 : CALL wann_real(
747 : > point,0,0,0,2,bkpt,nbn,
748 : > n3d,nmzxyd,n2d,ntypsd,lmaxd,jmtd,
749 : > natd,ntypd,nmzd,nop,nop2,mrot,tau,invtab,
750 : > nq3,nvac,invs,z1,delz,nmz,nmzxy,nq2,
751 : > lmax,rmsh,jri,pos,ngopr,ntypsy,nvd,
752 : > omtil,amat,bmat,odi,ods,nlod,llod,nlo,llo,
753 : > ff,gg,flo,acof(nbn,:,:),bcof(nbn,:,:),
754 : > ccof(:,nbn,:,:),zMat,
755 : > nv(jspin),k1(:,jspin),k2(:,jspin),k3(:,jspin),
756 : > lmd,nbasfcn,l_ss,qss,jspin,0,
757 0 : < xdnout)
758 : wannierfunc(:,posi)=
759 0 : =wannierfunc(:,posi)+xdnout*u_matrix(counter,:,ikpt)*factor
760 0 : knorm(ikpt,nbn)=knorm(ikpt,nbn)+(abs(xdnout))**2
761 : ENDDO xloop
762 : ENDDO
763 : ENDDO !z-loop
764 :
765 : c..end of the loop by the bands
766 0 : counter=counter+1
767 :
768 :
769 : ENDDO band
770 :
771 0 : deallocate ( acof,bcof,ccof,we,eigg )
772 :
773 0 : write (*,*) 'nslibd=',nslibd
774 :
775 :
776 0 : if(film)then
777 0 : if (.not.odi%d1) deallocate(ac,bc,u,ue)
778 0 : if (odi%d1) deallocate(ac_1,bc_1,u_1,ue_1)
779 : endif
780 :
781 :
782 : endif!processors
783 :
784 : enddo !loop over k-points
785 0 : mesh=grid(1)*grid(2)*grid(3)
786 :
787 : #ifdef CPP_MPI
788 : c call MPI_BARRIER(mpi_comm,ierr)
789 0 : if(l_p0)then
790 0 : if(isize.ne.1)then
791 0 : allocate(wannierfunc_temp(num_wann,mesh))
792 0 : do cpu_index=1,isize-1
793 0 : do ikpt=1,fullnkpts
794 0 : if(mod(ikpt-1,isize).eq.cpu_index)then
795 : call MPI_RECV(knorm(ikpt,1:num_bands),num_bands,
796 : & CPP_MPI_REAL,cpu_index,
797 0 : & ikpt,mpi_comm,stt,mpiierr)
798 : endif !processors
799 : enddo !ikpt
800 : call MPI_RECV(wannierfunc_temp(1:num_wann,1:mesh),
801 : & num_wann*mesh,
802 : & CPP_MPI_COMPLEX,cpu_index,
803 0 : & cpu_index+fullnkpts,mpi_comm,stt,mpiierr)
804 : wannierfunc(:,:)=wannierfunc(:,:)+
805 0 : & wannierfunc_temp(:,:)
806 :
807 :
808 :
809 : enddo !cpu_index
810 0 : deallocate(wannierfunc_temp)
811 : endif !isize
812 : else
813 0 : do ikpt=1,fullnkpts
814 0 : if(mod(ikpt-1,isize).eq.irank)then
815 : call MPI_SEND(knorm(ikpt,1:num_bands),num_bands,
816 0 : & CPP_MPI_REAL,0,ikpt,mpi_comm,mpiierr)
817 : endif !processors
818 : enddo !ikpt
819 : call MPI_SEND(wannierfunc(1:num_wann,1:mesh),
820 : & num_wann*mesh,
821 0 : & CPP_MPI_COMPLEX,0,fullnkpts+irank,mpi_comm,mpiierr)
822 :
823 :
824 :
825 : endif ! l_p0
826 :
827 :
828 :
829 : #endif
830 :
831 :
832 0 : deallocate(flo)
833 0 : if(l_p0)then
834 0 : wannierfunc(:,:)=wannierfunc(:,:)/real(fullnkpts)
835 0 : DO nplo=1,num_wann
836 :
837 :
838 : c****************************************************************
839 : c make Wannier function real (as much as possible)
840 : c****************************************************************
841 0 : phas=0.0
842 0 : do iz=0,grid(3)-1
843 0 : do iy=0,grid(2)-1
844 0 : do ix=0,grid(1)-1
845 0 : posi=ix+1+iy*(grid(1))+iz*(grid(1))*(grid(2))
846 0 : tmax=wannierfunc(nplo,posi)*conjg(wannierfunc(nplo,posi))
847 :
848 0 : if (tmax>phas) then
849 0 : phas=tmax
850 0 : factor=wannierfunc(nplo,posi)
851 : end if
852 : end do
853 : end do
854 : end do
855 0 : factor=factor/sqrt(real(factor)**2+aimag(factor)**2)
856 0 : wannierfunc(nplo,:)=wannierfunc(nplo,:)/factor
857 :
858 :
859 : c***************************************************************
860 : c open files for plot and make headers
861 : c***************************************************************
862 0 : IF (xsf) THEN
863 0 : write (name1,22) nplo,jspin
864 : 22 format (i3.3,'.real.',i1,'.xsf')
865 0 : write (name2,23) nplo,jspin
866 : 23 format (i3.3,'.imag.',i1,'.xsf')
867 0 : write (name3,24) nplo,jspin
868 : 24 format (i3.3,'.absv.',i1,'.xsf')
869 0 : OPEN(55,file=name1)
870 0 : CALL xsf_WRITE_atoms(55,atoms,film,odi%d1,amat)
871 0 : OPEN(56,file=name2)
872 0 : CALL xsf_WRITE_atoms(56,atoms,film,odi%d1,amat)
873 0 : OPEN(57,file=name3)
874 0 : CALL xsf_WRITE_atoms(57,atoms,film,odi%d1,amat)
875 : CALL xsf_WRITE_header(55,twodim,filename,vec1,vec2,vec3,zero
876 0 : $ ,grid)
877 : CALL xsf_WRITE_header(56,twodim,filename,vec1,vec2,vec3,zero
878 0 : $ ,grid)
879 : CALL xsf_WRITE_header(57,twodim,filename,vec1,vec2,vec3,zero
880 0 : $ ,grid)
881 : ELSE
882 0 : WRITE (vandername,201) nplo,jspin
883 : 201 FORMAT (i5.5,'.',i1)
884 0 : OPEN(55,file=vandername)
885 0 : WRITE (55,7) grid(1),grid(2),grid(3),ikpt,nslibd
886 : 7 FORMAT (5i4)
887 : ENDIF
888 : c********************************************************************
889 : c write data to files
890 : c********************************************************************
891 0 : DO iz = 0,grid(3)-1
892 0 : DO iy = 0,grid(2)-1
893 0 : DO ix = 0,grid(1)-1
894 0 : posi=ix+1+iy*grid(1)+iz*grid(1)*grid(2)
895 0 : IF (xsf) THEN
896 0 : WRITE(55,*) real(wannierfunc(nplo,posi))
897 0 : WRITE(56,*) aimag(wannierfunc(nplo,posi))
898 0 : WRITE(57,*) abs(wannierfunc(nplo,posi))
899 : ELSE
900 0 : WRITE(55,8) real(wannierfunc(nplo,posi))
901 : ENDIF
902 : enddo
903 : enddo
904 : enddo
905 0 : IF (xsf) THEN
906 0 : CALL xsf_WRITE_endblock(55,twodim)
907 0 : CALL xsf_WRITE_endblock(56,twodim)
908 0 : CALL xsf_WRITE_endblock(57,twodim)
909 0 : CLOSE (55) ; CLOSE (56) ; CLOSE (57)
910 : ENDIF
911 :
912 : ENDDO !nplo
913 0 : IF (.not.xsf) CLOSE(55)
914 : 8 FORMAT (2f7.3)
915 :
916 : c*******************************************************************
917 : c determine spreads, centers, norms
918 : c*******************************************************************
919 :
920 0 : l_spreadcal=.true.
921 : if(l_spreadcal)then !calculate spreads and centers from real space grid
922 0 : print*,"calculate spreads and centers"
923 0 : allocate (spreads(num_wann,num_wann))
924 0 : allocate (centers(3,num_wann))
925 0 : allocate(wfnorm(num_wann))
926 0 : centers(:,:)=0.0
927 0 : wfnorm(:)=0.0
928 0 : spreads(:,:)=cmplx(0.0,0.0)
929 0 : do nplo=1,num_wann
930 0 : do iz=0,grid(3)-1
931 0 : do iy=0,grid(2)-1
932 0 : do ix=0,grid(1)-1
933 0 : posi=ix+1+iy*(grid(1))+iz*(grid(1))*(grid(2))
934 : point = zero+vec1*(ix+0.0)/(grid(1)-1)+vec2*(iy+0.0)
935 0 : $ /(grid(2)-1)
936 0 : IF (.NOT.twodim) point = point+vec3*(iz+0.0)/(grid(3)-1)
937 : centers(:,nplo)=centers(:,nplo)
938 0 : + +point(:)*(abs(wannierfunc(nplo,posi)))**2
939 : wfnorm(nplo)=wfnorm(nplo)+
940 0 : + (abs(wannierfunc(nplo,posi)))**2
941 0 : do ii1=1,num_wann
942 : spreads(nplo,ii1)=spreads(nplo,ii1)
943 : + +wannierfunc(nplo,posi)*conjg(wannierfunc(ii1,posi))*
944 0 : * dot_product(point,point)
945 : enddo
946 : enddo
947 : enddo
948 : enddo
949 : enddo
950 :
951 0 : do nplo=1,num_wann !normalize centers
952 0 : centers(:,nplo)=centers(:,nplo)/(mesh)*volume
953 : enddo
954 :
955 0 : do nplo=1,num_wann !normalize spreads
956 0 : do ii1=1,num_wann
957 0 : spreads(nplo,ii1)=spreads(nplo,ii1)/mesh*volume
958 : enddo
959 : spreads(nplo,nplo)=spreads(nplo,nplo)
960 0 : & -dot_product(centers(1:3,nplo),centers(1:3,nplo))
961 : enddo
962 :
963 0 : wfnorm(:)=wfnorm(:)/(mesh)*volume !normalize wfnorm
964 :
965 0 : knorm(:,:)=knorm(:,:)/mesh*volume !normalize knorm
966 :
967 : c***********************************************************
968 : c write spreads and so on to files
969 : c***********************************************************
970 :
971 0 : open(518,file=spin12(jspin)//'.centers')
972 0 : do nplo=1,num_wann
973 0 : write(518,*)centers(:,nplo)
974 : enddo
975 0 : close(518)
976 0 : open(519,file=spin12(jspin)//'.spreads')
977 0 : do nplo=1,num_wann
978 0 : do ii1=1,num_wann
979 0 : write(519,*)nplo,ii1,spreads(nplo,ii1)
980 : enddo
981 : enddo
982 0 : close(519)
983 0 : open(521,file=spin12(jspin)//'.norm')
984 0 : do ii1=1,num_wann
985 0 : write(521,*)wfnorm(ii1)
986 : enddo
987 0 : close(521)
988 0 : deallocate(centers)
989 0 : deallocate(spreads)
990 0 : deallocate(wfnorm)
991 : endif
992 :
993 0 : open(611,file=spin12(jspin)//'.knorm')
994 0 : do nplo=1,num_bands
995 0 : do ikpt=1,fullnkpts
996 0 : write(611,*)ikpt,nplo,knorm(ikpt,nplo)
997 : enddo
998 : enddo
999 0 : close(611)
1000 :
1001 : c*****************************************************************
1002 : c*****************************************************************
1003 :
1004 : endif !l_p0
1005 :
1006 0 : deallocate(knorm)
1007 0 : if(have_disentangled)deallocate(inc_band)
1008 0 : if(.not.l_proj_plot)then
1009 0 : deallocate(lwindow,ndimwin)
1010 : endif
1011 0 : deallocate(wannierfunc)
1012 0 : nrec=nrec+nkpts
1013 0 : deallocate(u_matrix)
1014 :
1015 : #ifdef CPP_MPI
1016 0 : call MPI_BARRIER(mpi_comm,mpiierr)
1017 : #endif
1018 :
1019 0 : 110 continue ! end of cycle by spins
1020 :
1021 :
1022 0 : deallocate ( vr,vz,nv,k1,k2,k3 )
1023 0 : deallocate ( ff,gg)
1024 :
1025 : #ifdef CPP_MPI
1026 0 : call MPI_BARRIER(mpi_comm,mpiierr)
1027 : #endif
1028 :
1029 0 : END SUBROUTINE wann_plot_um_dat
1030 : END MODULE m_wann_plot_um_dat
|
