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

Generated by: LCOV version 2.0-1