LCOV - code coverage report
Current view: top level - wannier - wann_plot_um_dat.F (source / functions) Coverage Total Hit
Test: FLEUR test coverage Lines: 0.0 % 422 0
Test Date: 2025-06-14 04:34:23 Functions: 0.0 % 1 0

            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
        

Generated by: LCOV version 2.0-1