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

Generated by: LCOV version 1.13