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

Generated by: LCOV version 1.13