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