LCOV - code coverage report
Current view: top level - wannier/uhu - wann_uHu.F (source / functions) Hit Total Coverage
Test: combined.info Lines: 0 658 0.0 %
Date: 2019-09-08 04:53:50 Functions: 0 1 0.0 %

          Line data    Source code
       1             : c*******************************************c
       2             : c   Set up uHu matrix necessary for         c
       3             : c   Wannier-based calc. of orbital moment   c
       4             : c*******************************************c
       5             : c   keyword is 'matrixuhu' in wann_inp      c
       6             : c*******************************************c
       7             : c   uHu = < u_{k+b1} | H_{k} | u_{k+b2} >   c
       8             : c                                           c
       9             : c   Contributions to Hamiltonian:           c
      10             : c       (i)   interstitial                  c
      11             : c       (ii)  muffin tin  (a) spherical     c
      12             : c                         (b) non-sph.      c
      13             : c                         (c) SOC           c
      14             : c       (iii) vacuum                        c
      15             : c*******************************************c
      16             : c                  J.-P. Hanke, Dec. 2015   c
      17             : c*******************************************c
      18             :       MODULE m_wann_uHu
      19             :       USE m_juDFT
      20             :       CONTAINS
      21           0 :       SUBROUTINE wann_uHu(
      22             :      >      DIMENSION,stars,vacuum,atoms,sphhar,input,kpts,sym,mpi,
      23             :      >      banddos,oneD,noco,cell,vTot,wann,
      24           0 :      >      eig_idList,l_real,l_dulo,l_noco,l_ss,lmaxd,ntypd,
      25             :      >      neigd,natd,nop,nvd,jspd,nbasfcn,llod,nlod,ntype,
      26           0 :      >      omtil,nlo,llo,lapw_l,invtab,mrot,ngopr,neq,lmax,
      27           0 :      >      invsat,invsatnr,nkpt,taual,rmt,amat,bmat,bbmat,alph,
      28           0 :      >      beta,qss,sk2,phi2,odi,ods,irank,isize,n3d,nmzxyd,nmzd,
      29           0 :      >      jmtd,nlhd,nq3,nvac,invs,invs2,film,nlh,jri,ntypsd,
      30           0 :      >      ntypsy,jspins,nkptd,dx,n2d,rmsh,e1s,e2s,ulo_der,
      31           0 :      >      ustep,ig,k1d,k2d,k3d,rgphs,slice,kk,nnne,
      32           0 :      >      z1,nv2d,nmzxy,nmz,delz,zrfs,ig2,area,tau,zatom,nq2,kv2,nop2,
      33             :      >      volint,symor,pos,ef,l_soc,
      34           0 :      >      memd,lnonsph,clnu,lmplmd,mlh,nmem,llh,lo1l,
      35             :      >      theta,phi,
      36             :      >      l_ms,l_sgwf,l_socgwf,aux_latt_const,
      37           0 :      >      param_file,param_vec,nparampts,param_alpha,l_dim)
      38             : 
      39             :       use m_types
      40             :       use m_wann_mmnk_symm
      41             :       use m_wann_rw_eig
      42             :       use m_abcof
      43             :       use m_radfun
      44             :       use m_radflo
      45             :       use m_cdnread
      46             :       use m_constants, only : pimach
      47             :       use m_wann_projmethod
      48             :       use m_wann_abinv
      49             :       use m_wann_kptsrotate
      50             :       use m_wann_read_inp
      51             :       use m_matmul,only : matmul3,matmul3r
      52             :       use m_wann_maxbnd
      53             :       use m_wann_uHu_tlmplm
      54             :       use m_wann_uHu_sph
      55             :       use m_wann_uHu_int
      56             :       use m_wann_uHu_soc
      57             :       use m_wann_uHu_vac
      58             :       use m_wann_uHu_od_vac
      59             :       use m_wann_uHu_util
      60             :       use m_wann_uHu_commat
      61             :       use m_wann_write_uHu
      62             :       USE m_eig66_io
      63             : 
      64             :       IMPLICIT NONE
      65             : #include "cpp_double.h"
      66             : #ifdef CPP_MPI
      67             :       include 'mpif.h'
      68             :       integer ierr(3)
      69             :       integer cpu_index
      70             :       integer stt(MPI_STATUS_SIZE)
      71             : #endif
      72             : 
      73             :       TYPE(t_dimension),INTENT(IN) :: DIMENSION
      74             :       TYPE(t_stars),INTENT(IN)     :: stars
      75             :       TYPE(t_vacuum),INTENT(IN)    :: vacuum
      76             :       TYPE(t_atoms),INTENT(IN)     :: atoms
      77             :       TYPE(t_sphhar),INTENT(IN)    :: sphhar
      78             :       TYPE(t_input),INTENT(IN)     :: input
      79             :       TYPE(t_kpts),INTENT(IN)      :: kpts
      80             :       TYPE(t_sym),INTENT(IN)       :: sym
      81             :       TYPE(t_mpi),INTENT(IN)       :: mpi
      82             :       TYPE(t_banddos),INTENT(IN)   :: banddos
      83             :       TYPE(t_oneD),INTENT(IN)      :: oneD
      84             :       TYPE(t_noco),INTENT(IN)      :: noco
      85             :       TYPE(t_cell),INTENT(IN)      :: cell
      86             :       TYPE(t_potden),INTENT(IN)    :: vTot
      87             :       TYPE(t_wann),INTENT(INOUT)   :: wann
      88             : 
      89             : c     ..scalar arguments..
      90             :       character(len=20),intent(in) :: param_file
      91             :       type (od_inp), intent (in) :: odi
      92             :       type (od_sym), intent (in) :: ods
      93             :       INTEGER, INTENT (IN) :: eig_idList(wann%nparampts)
      94             :       logical, intent (in) :: invs,invs2,film,slice,symor,zrfs
      95             :       logical, intent (in) :: l_real,l_noco,l_ss,l_soc
      96             :       logical, intent (in) :: l_ms,l_sgwf,l_socgwf
      97             :       integer, intent (in) :: lmaxd,ntypd,neigd,nkptd,kk,nnne
      98             :       integer, intent (in) :: natd,nop,nvd,jspd,nbasfcn,nq2,nop2
      99             :       integer, intent (in) :: llod,nlod,ntype,n3d,n2d
     100             :       integer, intent (in) :: nmzxyd,nmzd,jmtd,nlhd,nq3,nvac
     101             :       integer, intent (in) :: ntypsd,jspins,k1d,k2d,k3d
     102             :       integer, intent (in) :: irank,isize,nv2d,nmzxy,nmz
     103             :       integer, intent (in) :: memd,lmplmd,nparampts
     104             :       real,    intent (in) :: omtil,e1s,e2s,delz,area,z1,volint
     105             :       real,    intent (in) :: ef,theta,phi,aux_latt_const
     106             : 
     107             : c     ..array arguments..
     108             :       logical, intent (in) :: l_dulo(nlod,ntypd)
     109             :       logical, intent (in) :: l_dim(3)
     110             :       integer, intent (in) :: ig(-k1d:k1d,-k2d:k2d,-k3d:k3d)
     111             :       integer, intent (in) :: nlh(ntypsd),jri(ntypd),ntypsy(natd)
     112             :       integer, intent (in) :: nlo(ntypd),llo(nlod,ntypd),lapw_l(ntypd)
     113             :       integer, intent (in) :: invtab(nop),mrot(3,3,nop),ngopr(natd)
     114             :       integer, intent (in) :: neq(ntypd),lmax(ntypd)
     115             :       integer, intent (in) :: invsat(natd),invsatnr(natd),nkpt
     116             :       integer, intent (in) :: ulo_der(nlod,ntypd),ig2(n3d),kv2(2,n2d)
     117             :       integer, intent (in) :: mlh(memd,0:nlhd,ntypsd)
     118             :       integer, intent (in) :: nmem(0:nlhd,ntypsd)
     119             :       integer, intent (in) :: llh(0:nlhd,ntypsd),lnonsph(ntypd)
     120             :       integer, intent (in) :: lo1l(0:llod,ntypd)
     121             :       complex, intent (in) :: rgphs(-k1d:k1d,-k2d:k2d,-k3d:k3d)
     122             :       real,    intent (in) :: taual(3,natd),rmt(ntypd),dx(ntypd)
     123             :       real,    intent (in) :: amat(3,3),bmat(3,3),bbmat(3,3)
     124             :       real,    intent (in) :: rmsh(jmtd,ntypd),tau(3,nop),zatom(ntype)
     125             :       real,    intent (in) :: alph(ntypd),beta(ntypd),qss(3)
     126             :       real,    intent (in) :: pos(3,natd),sk2(n2d),phi2(n2d)
     127             :       real,    intent (in) :: param_vec(3,nparampts)
     128             :       real,    intent (in) :: param_alpha(ntypd,nparampts)
     129             :       complex, intent (in) :: ustep(n3d)
     130             :       complex, intent (in) :: clnu(memd,0:nlhd,ntypsd)
     131             : 
     132             : c     ..allocatable arrays..
     133           0 :       integer, allocatable :: kveclo(:)   , nv(:)
     134           0 :       integer, allocatable :: kveclo_b(:) , nv_b(:)
     135           0 :       integer, allocatable :: kveclo_b2(:), nv_b2(:)
     136           0 :       integer, allocatable :: k1(:,:)   , k2(:,:)   , k3(:,:)
     137           0 :       integer, allocatable :: k1_b(:,:) , k2_b(:,:) , k3_b(:,:)
     138           0 :       integer, allocatable :: k1_b2(:,:), k2_b2(:,:), k3_b2(:,:)
     139           0 :       integer, allocatable :: irreduc(:),mapkoper(:)
     140           0 :       integer, allocatable :: irreduc_q(:),mapqoper(:)        
     141           0 :       integer, allocatable :: shiftkpt(:,:),pair_to_do(:,:)
     142           0 :       integer, allocatable :: shiftqpt(:,:),pair_to_do_q(:,:)  
     143           0 :       integer, allocatable :: maptopair(:,:,:)
     144           0 :       integer, allocatable :: maptopair_q(:,:,:)              
     145           0 :       integer, allocatable :: counts(:),displs(:)
     146           0 :       integer, allocatable :: gb(:,:,:),bpt(:,:)
     147           0 :       integer, allocatable :: gb_q(:,:,:),bpt_q(:,:)
     148           0 :       INTEGER, ALLOCATABLE :: innerEig_idList(:)
     149           0 :       real,    allocatable :: we(:),we_b(:),we_b2(:)
     150           0 :       real,    allocatable :: eigg(:)
     151           0 :       real,    allocatable :: vr(:,:,:),vz(:,:,:)
     152           0 :       real,    allocatable :: flo(:,:,:,:,:)
     153           0 :       real,    allocatable :: ff(:,:,:,:,:),gg(:,:,:,:,:)
     154           0 :       real,    allocatable :: kdiff(:,:),qdiff(:,:),zero_qdiff(:,:)
     155           0 :       complex, allocatable :: vpw(:,:),vzxy(:,:,:,:)
     156           0 :       complex, allocatable :: uHu(:,:,:,:,:)
     157             : c      complex, allocatable :: uHuold(:,:)
     158           0 :       complex, allocatable :: acof_b(:,:,:),acof_b2(:,:,:)
     159           0 :       complex, allocatable :: bcof_b(:,:,:),bcof_b2(:,:,:)
     160           0 :       complex, allocatable :: ccof_b(:,:,:,:),ccof_b2(:,:,:,:)
     161           0 :       complex, allocatable :: tdd(:,:,:,:,:),tdu(:,:,:,:,:)
     162           0 :       complex, allocatable :: tud(:,:,:,:,:),tuu(:,:,:,:,:)
     163           0 :       complex, allocatable :: tdulo(:,:,:,:,:,:),tuulo(:,:,:,:,:,:)
     164           0 :       complex, allocatable :: tulod(:,:,:,:,:,:),tulou(:,:,:,:,:,:)
     165           0 :       complex, allocatable :: tuloulo(:,:,:,:,:,:,:)
     166           0 :       complex, allocatable :: tdd_soc(:,:,:,:),tdu_soc(:,:,:,:)
     167           0 :       complex, allocatable :: tud_soc(:,:,:,:),tuu_soc(:,:,:,:)
     168           0 :       complex, allocatable :: tdulo_soc(:,:,:,:,:)
     169           0 :       complex, allocatable :: tuulo_soc(:,:,:,:,:)
     170           0 :       complex, allocatable :: tulod_soc(:,:,:,:,:)
     171           0 :       complex, allocatable :: tulou_soc(:,:,:,:,:)
     172           0 :       complex, allocatable :: tuloulo_soc(:,:,:,:,:,:)
     173             : 
     174             : c     ..local arrays..
     175             :       character(len=2) :: spin012(0:2)
     176             :       data spin012/'  ', '.1', '.2'/
     177             :       character(len=3) :: spin12(2)
     178             :       data   spin12/'WF1' , 'WF2'/
     179             :       character(len=8) :: name(10)
     180           0 :       integer :: n_bands(0:neigd),ngopr1(natd)
     181             :       real    :: bkpt(3),bkpt_b(3),bkpt_b2(3),bkrot(3)
     182           0 :       real    :: eig(neigd),eig_b(neigd),eig_b2(neigd)
     183           0 :       real    :: uuilon(nlod,ntypd),duilon(nlod,ntypd)
     184           0 :       real    :: ulouilopn(nlod,nlod,ntypd)
     185           0 :       real    :: ello(nlod,ntypd,max(2,jspd)),evac(2,max(2,jspd))
     186           0 :       real    :: epar(0:lmaxd,ntypd,max(2,jspd)),evdu(2,max(jspd,2))
     187             :       real    :: qpt_i(3),qptb_i(3)
     188           0 :       real    :: alph_i(ntypd),alphb_i(ntypd)
     189           0 :       real    :: beta_i(ntypd),betab_i(ntypd)
     190             :       real    :: cp_time(9)
     191             : 
     192             : c     ..local scalars..
     193             :       character(len=6) :: filename
     194             :       character(len=8) :: dop,iop
     195             :       character(len=12) fending
     196             :       character(len=30) fstart
     197             :       logical :: l_p0,l_bkpts,l_proj,l_file
     198             :       logical :: l_bqpts,l_gwf,l_exist,l_nocosoc,l_symcheck
     199             :       logical :: l_skip_sph,l_skip_non,l_skip_soc
     200             :       logical :: l_skip_int,l_skip_vac,l_skip_loc
     201             :       integer :: lmd,nlotot,n,iter,ikpt,ikpt_b,ikpt_b2,iqpt,iqpt_b
     202             :       integer :: addnoco,addnoco2,funbas,loplod,igvm2
     203             :       integer :: nn,nkpts,i,j,l,i_rec,m,nwf,nwfp
     204             :       integer :: jsp_start,jsp_end,nrec,nrec_b,nrec1
     205             :       integer :: nodeu,noded,n_size,na,n_rank,nbnd,numbands
     206             :       integer :: i1,i2,i3,in,lda
     207             :       integer :: nmat,nmat_b,nmat_b2,nmat_qb
     208             :       integer :: nbands,nbands_b,nbands_b2,nbands_qb
     209             :       integer :: nslibd,nslibd_b,nslibd_b2,nslibd_qb
     210             :       integer :: noccbd,noccbd_b,noccbd_b2,noccbd_qb
     211             :       integer :: kptibz,kptibz_b,kptibz_b2
     212             :       integer :: qptibz, qptibz_b
     213             :       integer :: oper,oper_b,oper_b2,oper_q, oper_qb
     214             :       integer :: nwfs,nntot,nntot_q,fullnkpts,fullnqpts
     215             :       integer :: kpt,qpt,j1,j2,j3,k,ikpt_help,iqpt_help
     216             :       integer :: wannierspin,jspin,jspin_b,jspin2
     217             :       integer :: jspin3,jspin4_b,jspin4,jspin5,tspin,tspin2
     218             :       integer :: n_start,n_end,mlotot,mlolotot,err
     219             :       integer :: mlot_d,mlolot_d,ilo,dir,length
     220             :       integer :: npotmatfile,ig3,maxvac,irec,imz,ivac,ipot
     221             :       integer :: funit_start,band_help,sign2
     222             :       integer :: doublespin,doublespin_max,nrec5
     223             :       integer :: aoff,d1,d10,d100
     224             :       integer :: eig_id
     225             :       real    :: tpi,wronk,wk,wk_b,wk_b2
     226             :       real    :: t0,t00,t1,t_myTlmplm,t_init
     227             :       real    :: t_int,t_sph,t_vac,t_abcof,t_eig,t_total
     228             :       real    :: efermi,htr2ev
     229             :       real    :: theta_i, thetab_i, phi_i, phib_i
     230             :       complex :: nsfactor,nsfactor_b,nsfactor_b2
     231             : 
     232           0 :       TYPE(t_usdus) :: usdus
     233           0 :       TYPE(t_zmat)  :: zMat, zzMat, zMat_b, zMat_b2
     234           0 :       TYPE(t_lapw)  :: lapw_b, lapw_b2
     235             : 
     236             : 
     237             : c     ..initializations..
     238           0 :       call cpu_time(t00)
     239             : 
     240           0 :       t_init = 0.
     241           0 :       t_myTlmplm = 0.
     242           0 :       t_eig = 0.
     243           0 :       t_abcof = 0.
     244           0 :       t_int = 0.
     245           0 :       t_sph = 0.
     246           0 :       t_vac = 0.
     247           0 :       t_total = 0.
     248           0 :       htr2ev = 27.2
     249           0 :       nntot_q = 1
     250           0 :       fullnqpts = 1
     251           0 :       funit_start = 5000
     252             :       
     253             :     ! In some modes the code skips the Hamiltonian setup and
     254             :     ! jumps directly into the Wannier branch.
     255             :       if(.not.allocated (usdus%us))then
     256           0 :          call usdus%init(atoms,input%jspins)
     257             :       endif
     258             : 
     259           0 :       aoff = iachar('1')-1
     260           0 :       d1  = mod(irank,10)
     261           0 :       IF (irank < 100) THEN
     262           0 :         d10 = int( (irank + 0.5)/10 )
     263           0 :         fstart = 'eig'//achar(d10+aoff)//achar(d1+aoff)
     264             :       ELSE
     265           0 :         d10 = mod((irank-d1)/10,10)
     266           0 :         d100 = (irank-10*d10-d1)/100
     267           0 :         IF ( d100.GE.10 ) d100 = d100 + iachar('7')
     268             :         fstart =
     269           0 :      +  'eig'//achar(d100+aoff)//achar(d10+aoff)//achar(d1+aoff)
     270             :       ENDIF
     271             : 
     272             : 
     273           0 :       ngopr1(:)=1
     274             : 
     275           0 :       l_p0 = .false.
     276           0 :       if (irank.eq.0) l_p0 = .true.
     277           0 :       l_nocosoc = l_noco.or.l_soc
     278             : 
     279           0 :       tpi = 2* pimach()
     280           0 :       lmd = lmaxd*(lmaxd+2)
     281             : 
     282             : !!!   should be changed in case the windows are really used
     283           0 :       nkpts = nkpt
     284             : 
     285             :       ! do we have to construct GWF ?
     286             :       l_gwf = .false.
     287           0 :       l_gwf = l_sgwf.or.l_socgwf 
     288             :       
     289             : 
     290             : c-----read the input file to determine what to do
     291           0 :       call wann_read_inp(input,l_p0,wann)
     292             : 
     293           0 :       if(wann%l_byenergy.and.wann%l_byindex) CALL juDFT_error
     294           0 :      +     ("byenergy.and.byindex",calledby ="wannier")
     295           0 :       if(wann%l_byenergy.and.wann%l_bynumber) CALL juDFT_error
     296           0 :      +     ("byenergy.and.bynumber",calledby ="wannier")
     297           0 :       if(wann%l_bynumber.and.wann%l_byindex) CALL juDFT_error
     298           0 :      +     ("bynumber.and.byindex",calledby ="wannier")
     299           0 :       if(.not.(wann%l_bynumber.or.wann%l_byindex.or.wann%l_byenergy))
     300           0 :      &     CALL juDFT_error("no rule to sort bands",calledby ="wannier")
     301             : 
     302             : 
     303           0 :       efermi=ef
     304           0 :       if(.not.wann%l_fermi)efermi=0.0
     305             : 
     306             : #ifdef CPP_MPI
     307           0 :          call MPI_BARRIER(MPI_COMM_WORLD,ierr)
     308             : #endif
     309             : 
     310             : c**************************************************************
     311             : c   for bzsym=.true.: determine mapping between kpts and w90kpts
     312             : c**************************************************************
     313           0 :       if (wann%l_bzsym) then
     314           0 :          l_file=.false.
     315           0 :          inquire(file='w90kpts',exist=l_file)
     316           0 :          if(.not.l_file)  CALL juDFT_error
     317           0 :      +        ("w90kpts not found, needed if bzsym",calledby ="wannier")
     318           0 :          open(412,file='w90kpts',form='formatted')
     319           0 :          read(412,*)fullnkpts
     320           0 :          close(412)
     321           0 :          if(l_p0)print*,"fullnkpts=",fullnkpts
     322           0 :          if(fullnkpts<nkpts) CALL juDFT_error("fullnkpts.lt.nkpts"
     323           0 :      +        ,calledby ="wannier")
     324           0 :          allocate(irreduc(fullnkpts),mapkoper(fullnkpts))
     325           0 :          allocate(shiftkpt(3,fullnkpts))
     326           0 :          l_file=.false.
     327           0 :          inquire(file='kptsmap',exist=l_file)
     328           0 :          if(.not.l_file)  CALL juDFT_error
     329           0 :      +        ("kptsmap not found, needed if bzsym",calledby ="wannier")
     330           0 :          open(713,file='kptsmap')
     331           0 :          do i=1,fullnkpts
     332           0 :             read(713,*)kpt,irreduc(i),mapkoper(i),shiftkpt(:,i)
     333           0 :             if(kpt/=i) CALL juDFT_error("kpt.ne.i",calledby ="wannier")
     334           0 :             if(l_p0)print*,i,irreduc(i),mapkoper(i)
     335             :          enddo
     336           0 :          close(713)
     337           0 :          if(maxval(irreduc(:))/=nkpts) CALL juDFT_error
     338           0 :      +        ("max(irreduc(:))/=nkpts",calledby ="wannier")
     339             :       else
     340           0 :          fullnkpts=nkpts
     341             :       endif
     342             : 
     343             : 
     344           0 :       if(l_gwf) fullnqpts = nparampts
     345             : 
     346             : 
     347           0 :       nrec = 0
     348           0 :       if(l_p0)then
     349           0 :       write (*,*) 'fermi energy:',efermi
     350           0 :       write (*,*) 'emin,emax=',e1s,e2s
     351           0 :       write (*,*) 'nbasfcn =',nbasfcn
     352             :       endif
     353           0 :       nlotot = 0
     354           0 :       mlotot = 0
     355           0 :       mlolotot = 0
     356           0 :       do n = 1, ntype
     357           0 :         mlotot = mlotot + nlo(n)
     358           0 :         mlolotot = mlolotot + nlo(n)*(nlo(n)+1)/2
     359           0 :         do l = 1,nlo(n)
     360           0 :           nlotot = nlotot + neq(n) * ( 2*llo(l,n) + 1 )
     361             :         enddo
     362             :       enddo
     363             : 
     364             : 
     365           0 :       allocate(counts(0:isize-1),displs(0:isize-1))
     366           0 :       call array_split(fullnkpts,isize,counts,displs)
     367             : 
     368             : c**********************************************************
     369             : ccccccccccccccc   read in the bkpts file  ccccccccccccccccc
     370             : c**********************************************************
     371           0 :        l_bkpts = .false.
     372           0 :        inquire (file='bkpts',exist=l_bkpts)
     373           0 :        if (.not.l_bkpts)  CALL juDFT_error("need bkpts for matrixmmn"
     374           0 :      +      ,calledby ="wannier")
     375           0 :        open (202,file='bkpts',form='formatted',status='old')
     376           0 :        rewind (202)
     377           0 :        read (202,'(i4)') nntot
     378           0 :        if(l_p0)then
     379           0 :        write (*,*) 'nntot=',nntot
     380           0 :        write(*,*) 'fullnkpts=',fullnkpts
     381           0 :        write(*,*) 'nkpts=',nkpts
     382             :        endif
     383           0 :        allocate ( gb(1:3,1:nntot,1:fullnkpts),bpt(1:nntot,1:fullnkpts))
     384           0 :        do ikpt=1,fullnkpts
     385           0 :         do nn=1,nntot
     386             :          read (202,'(2i6,3x,3i4)')
     387           0 :      &     ikpt_help,bpt(nn,ikpt),(gb(i,nn,ikpt),i=1,3)
     388           0 :          if (ikpt/=ikpt_help)  CALL juDFT_error("ikpt.ne.ikpt_help"
     389           0 :      +        ,calledby ="wannier")       
     390           0 :          if (bpt(nn,ikpt)>fullnkpts) CALL juDFT_error("bpt.gt.fullnkpts"
     391           0 :      +        ,calledby ="wannier")
     392             :         enddo
     393             :        enddo
     394           0 :        close (202)
     395           0 :         allocate(kdiff(3,nntot))
     396             : 
     397             : c**********************************************************
     398             : ccccccccccccccc   read in the bqpts file  ccccccccccccccccc         
     399             : c**********************************************************
     400           0 :       if (l_gwf.or.l_ms) then ! for Omega functional minimization
     401           0 :        l_bqpts = .false.
     402           0 :        inquire (file='bqpts',exist=l_bqpts)
     403           0 :        if (.not.l_bqpts)  CALL juDFT_error("need bqpts for matrixmmn"
     404           0 :      +     ,calledby ="wannier")
     405           0 :        open (202,file='bqpts',form='formatted',status='old')
     406           0 :        rewind (202)
     407           0 :        read (202,'(i4)') nntot_q
     408           0 :        if(l_p0)then
     409           0 :        write (*,*) 'nntot_q=',nntot_q
     410           0 :        write(*,*) 'fullnqpts=',fullnqpts
     411             :        endif
     412             :        allocate ( gb_q(1:3,1:nntot_q,1:fullnqpts),
     413           0 :      &            bpt_q(1:nntot_q,1:fullnqpts))
     414           0 :        do iqpt=1,fullnqpts
     415           0 :         do nn=1,nntot_q
     416             :          read (202,'(2i6,3x,3i4)')
     417           0 :      &     iqpt_help,bpt_q(nn,iqpt),(gb_q(i,nn,iqpt),i=1,3)
     418           0 :          if (iqpt/=iqpt_help)  CALL juDFT_error("iqpt.ne.iqpt_help"
     419           0 :      +        ,calledby ="wannier")
     420           0 :          if (bpt_q(nn,iqpt)>fullnqpts)
     421           0 :      &        CALL juDFT_error("bpt_q.gt.fullnqpts",calledby ="wannier")
     422             :         enddo
     423             :        enddo
     424           0 :        close (202)
     425           0 :         allocate(qdiff(3,nntot_q))
     426           0 :         allocate(zero_qdiff(3,nntot_q))
     427           0 :         zero_qdiff=0.0
     428             :       endif                                                        
     429             : 
     430             : 
     431             : ! when treating gen. WF for spin spirals, the Brillouin zone
     432             : ! of q-points is twice as large compared to k-BZ. Thus,
     433             : ! the G-vectors connecting neighbors across the boundary
     434             : ! need to be doubled
     435           0 :       if(l_sgwf) gb_q = 2*gb_q    
     436           0 :       if(l_socgwf) gb_q = 2*gb_q 
     437             : 
     438           0 :       if(wann%l_finishgwf) goto 9110
     439             : c********************************************************
     440             : c      find symmetry-related elements in mmkb
     441             : c********************************************************
     442           0 :          allocate(maptopair(3,fullnkpts,nntot))
     443           0 :          allocate(pair_to_do(fullnkpts,nntot))
     444             :          call wann_mmnk_symm(input,kpts,
     445             :      >     fullnkpts,nntot,bpt,gb,wann%l_bzsym,
     446             :      >     irreduc,mapkoper,l_p0,film,nop,invtab,mrot,odi%d1,
     447             :      >     tau,
     448           0 :      <     pair_to_do,maptopair,kdiff,.false.,param_file)
     449             : 
     450             :       ! do the same for q-points to construct GWFs
     451           0 :       if(l_gwf)then 
     452           0 :          allocate(maptopair_q(3,fullnqpts,nntot_q))
     453           0 :          allocate(pair_to_do_q(fullnqpts,nntot_q))
     454             :          call wann_mmnk_symm(input,kpts,
     455             :      >     fullnqpts,nntot_q,bpt_q,gb_q,wann%l_bzsym,
     456             :      >     irreduc_q,mapqoper,l_p0,.false.,1,invtab(1),mrot(:,:,1),
     457             :      >     .false.,tau,
     458           0 :      <     pair_to_do_q,maptopair_q,qdiff,.true.,param_file)
     459             :       endif
     460             : 
     461             : 
     462             : c*********************************************************
     463             : cccccccccccccccc   initialize the potential   cccccccccccc
     464             : c*********************************************************
     465             : 
     466           0 :       if(.not. l_noco) then
     467           0 :        allocate ( vpw(n3d,jspd),vzxy(nmzxyd,odi%n2d-1,2,jspd) )
     468             :       else
     469           0 :        allocate ( vpw(n3d,4),vzxy(nmzxyd,odi%n2d-1,2,4) )
     470             :       endif
     471             : 
     472           0 :       allocate (vz(nmzd,2,4))
     473           0 :       allocate (vr(jmtd,ntypd,jspd))
     474             : 
     475           0 :       vpw(:,1:SIZE(vTot%pw,2)) = vTot%pw(:,1:SIZE(vTot%pw,2))
     476           0 :       IF(film) THEN
     477             :          vz(:,:,1:SIZE(vTot%vacz,3)) =
     478           0 :      &      vTot%vacz(:,:,1:SIZE(vTot%vacz,3))
     479             :          vzxy(:,:,:,1:SIZE(vTot%vacxy,4)) =
     480           0 :      &      vTot%vacxy(:,:,:,1:SIZE(vTot%vacxy,4))
     481             :       END IF
     482             : 
     483           0 :       do jspin = 1,jspins
     484           0 :         do n = 1, ntype
     485           0 :           do j = 1,jri(n)
     486           0 :             vr(j,n,jspin) = vTot%mt(j,0,n,jspin)
     487             :           enddo
     488             :         enddo
     489             :       enddo
     490             : 
     491           0 :       if(.not.film) deallocate(vz,vzxy)
     492             : 
     493           0 :       if(l_noco)then
     494           0 :          npotmatfile=25
     495             : 
     496             :          OPEN (npotmatfile,FILE='potmat',FORM='unformatted',
     497           0 :      +                  STATUS='old')
     498             : c--->    load the interstitial potential
     499           0 :          if(l_p0 .and. nq3.ne.n3d) then
     500           0 :           write(*,*)'WARNING for reading potmat: nq3.ne.n3d'
     501             :          endif
     502           0 :          vpw = cmplx(0.0, 0.0)
     503           0 :          READ (npotmatfile) (vpw(ig3,1),ig3=1,n3d)
     504           0 :          READ (npotmatfile) (vpw(ig3,2),ig3=1,n3d)
     505           0 :          READ (npotmatfile) (vpw(ig3,3),ig3=1,n3d)
     506           0 :          vpw(:,4) = conjg(vpw(:,3))
     507           0 :          if(film) then
     508           0 :           maxvac=2
     509           0 :           if(odi%d1)maxvac=1
     510           0 :           DO ivac = 1,maxvac
     511             : c--->       if the two vacuua are equivalent, the potential file has to
     512             : c--->       be backspaced, because the potential is the same at both
     513             : c--->       surfaces of the film
     514           0 :             IF ((ivac.EQ.2) .AND. (nvac.EQ.1)) THEN
     515           0 :                DO irec = 1,4
     516           0 :                   BACKSPACE (npotmatfile)
     517             :                ENDDO
     518             :             ENDIF
     519             : c--->       load the non-warping part of the potential
     520           0 :             READ (npotmatfile)((vz(imz,ivac,ipot),imz=1,nmzd),ipot=1,4)
     521             : 
     522             : c--->       load the warping part of the potential
     523           0 :             if(.not.odi%d1)then
     524           0 :                DO ipot = 1,3
     525           0 :                   READ (npotmatfile)((vzxy(imz,igvm2,ivac,ipot),
     526           0 :      +                           imz=1,nmzxy),igvm2=1,nq2-1)
     527             :                ENDDO   
     528             :             else
     529           0 :                DO ipot = 1,3
     530           0 :                   READ (npotmatfile)((vzxy(imz,igvm2,ivac,ipot),
     531           0 :      +                           imz=1,nmzxy),igvm2=1,odi%n2d-1)
     532             :                ENDDO   
     533             :             endif   
     534           0 :             vzxy(:,:,:,4) = conjg(vzxy(:,:,:,3))
     535             :           enddo   
     536             :          endif
     537           0 :          CLOSE (npotmatfile)
     538             :       endif   
     539             :    
     540             : 
     541           0 :       if(film .and. l_p0) write(*,*)'nvac',nvac
     542             : 
     543             : cccccccccccccccc   end of the potential part  ccccccccccc
     544           0 :       wannierspin=jspd
     545           0 :       if(l_soc) wannierspin=2
     546             :      
     547           0 :       allocate(flo(ntypd,jmtd,2,nlod,2))
     548           0 :       allocate ( ff(ntypd,jmtd,2,0:lmaxd,2) )
     549           0 :       allocate ( gg(ntypd,jmtd,2,0:lmaxd,2) )
     550           0 :       allocate ( kveclo(nlotot),nv(wannierspin) )
     551           0 :       allocate ( kveclo_b(nlotot),nv_b(wannierspin) )
     552           0 :       allocate ( kveclo_b2(nlotot),nv_b2(wannierspin) )
     553             :       allocate ( k1(nvd,wannierspin),k2(nvd,wannierspin),
     554           0 :      &           k3(nvd,wannierspin) )
     555           0 :       allocate ( k1_b(nvd,wannierspin),k2_b(nvd,wannierspin),
     556           0 :      &           k3_b(nvd,wannierspin) )
     557           0 :       allocate ( k1_b2(nvd,wannierspin),k2_b2(nvd,wannierspin),
     558           0 :      &           k3_b2(nvd,wannierspin) )
     559             : 
     560           0 :       if(l_nocosoc) then
     561           0 :          doublespin_max=4
     562             :       else
     563           0 :          doublespin_max=wannierspin
     564             :       endif
     565             :       
     566             : c      if(l_soc.and.(jspins.eq.1)) doublespin_max=wannierspin
     567             : 
     568           0 :       l_skip_int = .false.; l_skip_soc = .false.; l_skip_vac = .false.
     569           0 :       l_skip_sph = .false.; l_skip_non = .false.; l_skip_loc = .false.
     570           0 :       inquire(file='debug_uHu',exist=l_exist)
     571           0 :       if(l_exist) then
     572           0 :        open(888,file='debug_uHu')
     573           0 :        read(888,*)l_skip_int
     574           0 :        read(888,*)l_skip_sph
     575           0 :        read(888,*)l_skip_non
     576           0 :        read(888,*)l_skip_soc
     577           0 :        read(888,*)l_skip_loc
     578           0 :        read(888,*)l_skip_vac
     579           0 :        read(888,*)doublespin_max 
     580           0 :        close(888)
     581           0 :        if(l_p0) then
     582           0 :         write(*,*)'skip INT :',l_skip_int
     583           0 :         write(*,*)'skip SPH :',l_skip_sph
     584           0 :         write(*,*)'skip NON :',l_skip_non
     585           0 :         write(*,*)'skip SOC :',l_skip_soc
     586           0 :         write(*,*)'skip LOC :',l_skip_loc
     587           0 :         write(*,*)'skip VAC :',l_skip_vac
     588           0 :         write(*,*)'doublespin_max:',doublespin_max
     589             :        endif
     590             :       endif
     591             : 
     592           0 :       tspin = doublespin_max
     593           0 :       if(fullnqpts.eq.1) tspin=1
     594           0 :       if(l_p0) write(*,*)'tspin',tspin
     595           0 :       allocate( tdd(0:lmd,0:lmd,ntypd,nntot*nntot,tspin) )
     596           0 :       allocate( tdu(0:lmd,0:lmd,ntypd,nntot*nntot,tspin) )
     597           0 :       allocate( tud(0:lmd,0:lmd,ntypd,nntot*nntot,tspin) )
     598           0 :       allocate( tuu(0:lmd,0:lmd,ntypd,nntot*nntot,tspin) )
     599           0 :       allocate( tdulo(0:lmd,nlod,-llod:llod,ntypd,nntot*nntot,tspin) )
     600           0 :       allocate( tuulo(0:lmd,nlod,-llod:llod,ntypd,nntot*nntot,tspin) )
     601           0 :       allocate( tulou(0:lmd,nlod,-llod:llod,ntypd,nntot*nntot,tspin) )
     602           0 :       allocate( tulod(0:lmd,nlod,-llod:llod,ntypd,nntot*nntot,tspin) )
     603             :       allocate( tuloulo(nlod,-llod:llod,nlod,-llod:llod,
     604           0 :      >                  ntypd,nntot*nntot,tspin) )
     605           0 :       allocate( tdd_soc(0:lmd,0:lmd,ntypd,nntot*nntot) )
     606           0 :       allocate( tdu_soc(0:lmd,0:lmd,ntypd,nntot*nntot) )
     607           0 :       allocate( tud_soc(0:lmd,0:lmd,ntypd,nntot*nntot) )
     608           0 :       allocate( tuu_soc(0:lmd,0:lmd,ntypd,nntot*nntot) )
     609           0 :       allocate( tdulo_soc(0:lmd,nlod,-llod:llod,ntypd,nntot*nntot) )
     610           0 :       allocate( tuulo_soc(0:lmd,nlod,-llod:llod,ntypd,nntot*nntot) )
     611           0 :       allocate( tulou_soc(0:lmd,nlod,-llod:llod,ntypd,nntot*nntot) )
     612           0 :       allocate( tulod_soc(0:lmd,nlod,-llod:llod,ntypd,nntot*nntot) )
     613           0 :       allocate( tuloulo_soc(nlod,-llod:llod,nlod,-llod:llod,
     614           0 :      >                  ntypd,nntot*nntot) )
     615             : 
     616           0 :       tuu = cmplx(0.,0.); tdu = cmplx(0.,0.)
     617           0 :       tud = cmplx(0.,0.); tdd = cmplx(0.,0.)
     618           0 :       tuulo = cmplx(0.,0.); tdulo = cmplx(0.,0.)
     619           0 :       tulou = cmplx(0.,0.); tulod = cmplx(0.,0.)
     620           0 :       tuloulo = cmplx(0.,0.)
     621           0 :       tuu_soc = cmplx(0.,0.); tdu_soc = cmplx(0.,0.)
     622           0 :       tud_soc = cmplx(0.,0.); tdd_soc = cmplx(0.,0.)
     623           0 :       tuulo_soc = cmplx(0.,0.); tdulo_soc = cmplx(0.,0.)
     624           0 :       tulou_soc = cmplx(0.,0.); tulod_soc = cmplx(0.,0.)
     625           0 :       tuloulo_soc = cmplx(0.,0.)
     626             : 
     627           0 :       call cpu_time(t1)
     628           0 :       t_init = t1-t00
     629             : 
     630             : c*****************************************************************c
     631             : c                         START Q LOOP                            c
     632             : c*****************************************************************c
     633           0 :       do 314 iqpt = 1,fullnqpts  ! loop by q-points starts
     634             : 
     635           0 :        ALLOCATE(innerEig_idList(nntot_q))
     636             : 
     637           0 :         qptibz=iqpt                          
     638           0 :         if(wann%l_bzsym .AND. l_gwf) qptibz=irreduc_q(iqpt)
     639             :         if(wann%l_bzsym .AND. l_gwf) oper_q=mapqoper(iqpt)
     640             : 
     641           0 :        qpt_i = qss
     642           0 :        alph_i = alph
     643           0 :        beta_i = beta
     644           0 :        theta_i = theta
     645           0 :        phi_i = phi
     646           0 :        if(l_sgwf.or.l_ms) then
     647           0 :           qpt_i(:) = param_vec(:,qptibz)
     648           0 :           alph_i(:) = param_alpha(:,qptibz)
     649           0 :        elseif(l_socgwf) then 
     650           0 :           if(l_dim(2)) phi_i = tpi*param_vec(2,qptibz)
     651           0 :           if(l_dim(3)) theta_i = tpi*param_vec(3,qptibz)
     652             :        endif
     653             : 
     654           0 :        if(l_p0) then
     655           0 :         write(*,*)'qpt_i',qpt_i
     656           0 :         do n=1,ntype
     657           0 :          write(*,*)'n=',n,alph_i(n),beta_i(n)
     658             :         enddo
     659           0 :         write(*,*)'theta_i=',theta_i
     660           0 :         write(*,*)'phi_i=',phi_i
     661             :        endif
     662             : 
     663           0 :        IF (l_gwf) THEN
     664           0 :           do iqpt_b=1,nntot_q
     665             : 
     666           0 :             innerEig_idList(iqpt_b) = eig_idList(bpt_q(iqpt_b,iqpt))
     667             : 
     668             : !            WRITE(fending,'("_",i4.4)')bpt_q(iqpt_b,iqpt)
     669             : !            innerEig_idList(iqpt_b)=open_eig(mpi%mpi_comm,
     670             : !     +                  DIMENSION%nbasfcn,DIMENSION%neigd,
     671             : !     +                  nkpts,wannierspin,atoms%lmaxd,
     672             : !     +                  atoms%nlod,atoms%ntype,atoms%nlotot,
     673             : !     +                  l_noco,.FALSE.,l_real,l_soc,.FALSE.,
     674             : !     +                  mpi%n_size,filename=trim(fstart)//fending,
     675             : !     +                  layers=vacuum%layers,nstars=vacuum%nstars,
     676             : !     +                  ncored=DIMENSION%nstd,nsld=atoms%nat,
     677             : !     +                  nat=atoms%nat,l_dos=banddos%dos.OR.input%cdinf,
     678             : !     +                  l_mcd=banddos%l_mcd,l_orb=banddos%l_orb)
     679             : 
     680             :           enddo
     681             : 
     682           0 :           eig_id = eig_idList(qptibz)
     683             : 
     684             : !          WRITE(fending,'("_",i4.4)')qptibz
     685             : !         eig_id=open_eig(mpi%mpi_comm,DIMENSION%nbasfcn,DIMENSION%neigd,
     686             : !     +                   nkpts,wannierspin,atoms%lmaxd,
     687             : !     +                   atoms%nlod,atoms%ntype,atoms%nlotot,
     688             : !     +                   l_noco,.FALSE.,l_real,l_soc,.FALSE.,
     689             : !     +                   mpi%n_size,filename=trim(fstart)//fending,
     690             : !     +                   layers=vacuum%layers,nstars=vacuum%nstars,
     691             : !     +                   ncored=DIMENSION%nstd,nsld=atoms%nat,
     692             : !     +                   nat=atoms%nat,l_dos=banddos%dos.OR.input%cdinf,
     693             : !     +                   l_mcd=banddos%l_mcd,l_orb=banddos%l_orb)
     694             : 
     695             : !        OPEN(66,file=trim(fstart)//fending,access='direct',
     696             : !     >          form='unformatted',recl=irecl,status='old')
     697           0 :        ELSEIF(l_ms) THEN
     698             : 
     699           0 :           eig_id = eig_idList(qptibz)
     700             : 
     701             : !        WRITE(fending,'("_",i4.4)')qptibz
     702             : !        eig_id=open_eig(mpi%mpi_comm,DIMENSION%nbasfcn,DIMENSION%neigd,
     703             : !     +                  nkpts,wannierspin,atoms%lmaxd,
     704             : !     +                  atoms%nlod,atoms%ntype,atoms%nlotot,
     705             : !     +                  l_noco,.FALSE.,l_real,l_soc,.FALSE.,
     706             : !     +                  mpi%n_size,filename=trim(fstart)//fending,
     707             : !     +                  layers=vacuum%layers,nstars=vacuum%nstars,
     708             : !     +                  ncored=DIMENSION%nstd,nsld=atoms%nat,
     709             : !     +                  nat=atoms%nat,l_dos=banddos%dos.OR.input%cdinf,
     710             : !     +                  l_mcd=banddos%l_mcd,l_orb=banddos%l_orb)
     711             : 
     712             : !        OPEN(66,file=trim(fstart)//fending,access='direct',
     713             : !     >          form='unformatted',recl=irecl,status='old')
     714             :        ELSE
     715           0 :           fending=''
     716             :        ENDIF ! l_gwf.or.l_ms
     717           0 :        nrec=0
     718           0 :        nrec_b=0
     719             : 
     720             : c****************************************************
     721             : c cycle by spins starts! 
     722             : c****************************************************
     723           0 :       do 110 doublespin=1,doublespin_max   ! cycle by spins
     724           0 :          if(l_p0) write(*,*)'spin loop:',doublespin
     725             : 
     726             : c         jspin=mod(doublespin+1,2)+1
     727             : c         jspin_b=jspin
     728             : c         if(doublespin.eq.3) jspin_b=2
     729             : c         if(doublespin.eq.4) jspin_b=1
     730           0 :          jspin_b=mod(doublespin+1,2)+1
     731           0 :          jspin=jspin_b
     732           0 :          if(doublespin.eq.3) jspin=2
     733           0 :          if(doublespin.eq.4) jspin=1
     734             : 
     735           0 :          tspin2 = doublespin
     736           0 :          if(fullnqpts.eq.1) tspin2 = 1
     737             : 
     738           0 :          nrec_b = nrec
     739             : 
     740           0 :          if(.not.l_noco) then
     741           0 :             nrec = (jspin-1)*nkpts
     742           0 :             nrec_b = (jspin_b-1)*nkpts
     743             :          endif
     744             : 
     745             : c...read number of bands and wannier functions from file proj
     746             : 
     747             : c..reading the proj.1 / proj.2 / proj file
     748           0 :        l_proj=.false.  
     749           0 :        do j=jspin,0,-1
     750           0 :          inquire(file=trim('proj'//spin012(j)),exist=l_proj)
     751           0 :          if(l_proj)then
     752           0 :             filename='proj'//spin012(j)
     753           0 :             exit
     754             :          endif
     755             :        enddo
     756             : 
     757           0 :        if(l_proj)then
     758           0 :          open (203,file=trim(filename),status='old')
     759           0 :          rewind (203)
     760           0 :          read (203,*) nwfs,numbands
     761           0 :          rewind (203)
     762           0 :          close (203)
     763             :        elseif(wann%l_projmethod.or.wann%l_bestproj
     764           0 :      &                         .or.wann%l_matrixamn)then
     765           0 :           CALL juDFT_error("no proj/proj.1/proj.2",calledby ="wannier")
     766             :        endif  
     767             : 
     768             : 
     769           0 :        jspin2=jspin
     770             :        if(l_soc .and. jspins.eq.1)jspin2=1
     771           0 :        jsp_start = jspin ; jsp_end = jspin
     772             : 
     773             : cccccccccccc   read in the eigenvalues and vectors   cccccc
     774           0 :        do jspin5=1,wannierspin!1,2
     775           0 :        jsp_start=jspin5; jsp_end=jspin5
     776           0 :        nrec5=0
     777             :        if(.not.l_noco) nrec5 = (jspin5-1)*nkpts
     778           0 :        CALL judft_error("TODO:wann_uHu")
     779             : !       call cdn_read0(eig_id,irank,isize,jspin5,wannierspin,l_noco,
     780             : !     <                ello,evac,epar,wk,n_bands,n_size)
     781             : 
     782             :        enddo
     783             : 
     784             : c..   now we want to define the maximum number of the bands by all kpts
     785           0 :       nbnd = 0
     786           0 :       i_rec = 0 ; n_rank = 0
     787             : 
     788           0 :       if(l_p0)then         
     789             :          call wann_maxbnd(
     790             :      >            eig_id,l_real,
     791             :      >            lmaxd,ntypd,nlod,neigd,nvd,wannierspin,
     792             :      >            isize,jspin,nbasfcn,nlotot,
     793             :      >            l_ss,l_noco,nrec,fullnkpts,
     794             :      >            wann%l_bzsym,wann%l_byindex,wann%l_bynumber,
     795             :      >            wann%l_byenergy,
     796             :      >            irreduc,odi,wann%band_min(jspin),
     797             :      >            wann%band_max(jspin),
     798             :      >            numbands,e1s,e2s,efermi,nkpts,
     799           0 :      <            nbnd,l_gwf,iqpt)       
     800             :       endif!l_p0
     801             : 
     802             : ! nbnd is calculated for process zero and is sent here to the others
     803             : #ifdef CPP_MPI
     804           0 :       if(l_p0)then
     805           0 :          do cpu_index=1,isize-1
     806           0 :       call MPI_SEND(nbnd,1,MPI_INTEGER,cpu_index,1,MPI_COMM_WORLD,ierr)
     807             :          enddo
     808             :       else
     809           0 :        call MPI_RECV(nbnd,1,MPI_INTEGER,0,1,MPI_COMM_WORLD,stt,ierr)
     810             :       endif
     811             : #endif
     812             :      
     813             : c##################################################################
     814           0 :          if(.not.allocated(uHu)) then
     815           0 :             allocate(uHu(nbnd,nbnd,nntot,nntot,counts(irank)))
     816             : c            allocate(uHuold(nbnd,nbnd))
     817           0 :             uHu = cmplx(0.,0.)
     818             : c            uHuold = cmplx(0.,0.)
     819             :          endif
     820             : 
     821             : 
     822             : 
     823             : ! set up T-matrix for spherical and non-spherical part
     824             : ! only once since it is k- and q-independent !
     825           0 :       if(iqpt.eq.1) then
     826           0 :          if(l_p0) write(*,*)'tspin2=',tspin2
     827           0 :          tuu(:,:,:,:,tspin2) = cmplx(0.,0.)
     828           0 :          tdu(:,:,:,:,tspin2) = cmplx(0.,0.)
     829           0 :          tud(:,:,:,:,tspin2) = cmplx(0.,0.)
     830           0 :          tdd(:,:,:,:,tspin2) = cmplx(0.,0.)
     831           0 :          tuulo(:,:,:,:,:,tspin2) = cmplx(0.,0.)
     832           0 :          tdulo(:,:,:,:,:,tspin2) = cmplx(0.,0.)
     833           0 :          tulou(:,:,:,:,:,tspin2) = cmplx(0.,0.)
     834           0 :          tulod(:,:,:,:,:,tspin2) = cmplx(0.,0.)
     835           0 :          tuloulo(:,:,:,:,:,:,tspin2) = cmplx(0.,0.)
     836             : 
     837           0 :       do jspin4=1,2
     838           0 :          jspin3=jspin4
     839           0 :          if(jspins.eq.1) jspin3=1
     840           0 :       na = 1
     841           0 :       do 40 n = 1,ntype
     842           0 :        do 30 l = 0,lmax(n)
     843             : c...compute the l-dependent, k-independent radial MT- basis functions
     844             : 
     845             :          call radfun(
     846             :      >              l,n,jspin4,epar(l,n,jspin3),vr(1,n,jspin3),atoms,
     847             :      <              ff(n,:,:,l,jspin4),gg(n,:,:,l,jspin4),usdus,
     848           0 :      <              nodeu,noded,wronk)
     849             : 
     850           0 :    30  continue
     851             : c...and the local orbital radial functions
     852           0 :        do ilo = 1, nlo(n)
     853             : 
     854             :          call radflo(
     855             :      >             atoms,n,jspin4,ello(:,:,jspin3),vr(1,n,jspin3),
     856             :      >             ff(n,1:,1:,0:,jspin4),gg(n,1:,1:,0:,jspin4),mpi,
     857           0 :      <             usdus,uuilon,duilon,ulouilopn,flo(n,:,:,:,jspin4))
     858             : 
     859             :        enddo
     860             : c       na = na + neq(n)
     861           0 :    40 continue
     862             :       enddo!jspin3
     863             : 
     864             : 
     865           0 :       if(jspin.eq.jspin_b) then
     866           0 :          if(l_p0) write(*,*)'wann_uHu_tlmplm'
     867           0 :          jspin3=jspin
     868           0 :          if(jspins.eq.1) jspin3=1
     869           0 :          call cpu_time(t0)
     870             :          call wann_uHu_tlmplm(
     871             :      >        memd,nlhd,ntypsd,ntypd,jmtd,lmaxd,jspd,ntype,dx,
     872             :      >        rmsh,jri,lmax,ntypsy,natd,lnonsph,lmd,lmplmd,clnu,
     873             :      >        mlh,nmem,llh,nlh,neq,irank,mlotot,mlolotot,
     874             :      >        vTot%mt(:,:,:,jspin3),nlod,llod,loplod,ello(1,1,jspin3),
     875             :      >        llo,nlo,lo1l,l_dulo,ulo_der,ff(:,:,:,:,jspin),
     876             :      >        gg(:,:,:,:,jspin),flo(:,:,:,:,jspin),
     877             :      >        ff(:,:,:,:,jspin_b),gg(:,:,:,:,jspin_b),
     878             :      >        flo(:,:,:,:,jspin_b),kdiff,kdiff,nntot,nntot,bmat,bbmat,
     879             :      >        vr(1,1,jspin3),epar(0,1,jspin3),invsat,
     880             :      >        l_skip_sph,l_skip_non,l_skip_loc,
     881             :      <        tuu(:,:,:,:,tspin2),tud(:,:,:,:,tspin2),
     882             :      >        tdu(:,:,:,:,tspin2),tdd(:,:,:,:,tspin2),
     883             :      >        tuulo(:,:,:,:,:,tspin2),tulou(:,:,:,:,:,tspin2),
     884             :      >        tdulo(:,:,:,:,:,tspin2),tulod(:,:,:,:,:,tspin2),
     885           0 :      >        tuloulo(:,:,:,:,:,:,tspin2))
     886           0 :          call cpu_time(t1)
     887           0 :          t_myTlmplm = t_myTlmplm + t1-t0
     888             :       endif
     889             : 
     890             :       endif!iqpt.eq.1
     891             : 
     892             : 
     893             : ! compute SOC-contribution for each theta_i !
     894           0 :       if(l_soc.and. (.not.l_skip_soc)) then
     895           0 :          if(l_p0) write(*,*)'wann_uHu_soc'
     896           0 :          tuu_soc = cmplx(0.,0.); tdu_soc = cmplx(0.,0.)
     897           0 :          tud_soc = cmplx(0.,0.); tdd_soc = cmplx(0.,0.)
     898           0 :          tuulo_soc = cmplx(0.,0.); tdulo_soc = cmplx(0.,0.)
     899           0 :          tulou_soc = cmplx(0.,0.); tulod_soc = cmplx(0.,0.)
     900           0 :          tuloulo_soc = cmplx(0.,0.)
     901           0 :          call cpu_time(t0)
     902             :          call wann_uHu_soc(
     903             :      >     input,atoms,
     904             :      >     ntypd,jmtd,lmaxd,jspd,
     905             :      >     ntype,dx,rmsh,jri,lmax,natd,
     906             :      >     lmd,lmplmd,neq,irank,
     907             :      >     nlod,llod,loplod,ello,llo,nlo,lo1l,l_dulo,ulo_der,
     908             :      >     ff(:,:,:,:,jspin),gg(:,:,:,:,jspin),flo(:,:,:,:,jspin),
     909             :      >     ff(:,:,:,:,jspin_b),gg(:,:,:,:,jspin_b),flo(:,:,:,:,jspin_b),
     910             :      >     kdiff,kdiff,nntot,nntot,bmat,bbmat,
     911             :      >     vr,epar,jspin,jspin_b,jspins,
     912             :      >     .true.,theta_i,phi_i,alph_i,beta_i,l_noco,l_skip_loc,
     913             :      <     tuu_soc,tud_soc,
     914             :      >     tdu_soc,tdd_soc,
     915             :      >     tuulo_soc,tulou_soc,
     916             :      >     tdulo_soc,tulod_soc,
     917           0 :      >     tuloulo_soc)
     918           0 :          call cpu_time(t1)
     919           0 :          t_myTlmplm = t_myTlmplm + t1-t0
     920             :       endif!l_soc
     921             : 
     922             : 
     923           0 :       i_rec = 0 ; n_rank = 0
     924             :       
     925             : c****************************************************************
     926             : c.. loop by kpoints starts!      each may be a separate task
     927             : c****************************************************************
     928           0 :       if(l_p0) write(*,*)'start k-loop'
     929           0 :       do 10 ikpt = wann%ikptstart,fullnkpts  ! loop by k-points starts
     930           0 :         kptibz=ikpt
     931             :         if(wann%l_bzsym) kptibz=irreduc(ikpt)
     932             :         if(wann%l_bzsym) oper=mapkoper(ikpt)
     933             : 
     934           0 :         if(kpt_on_node(ikpt,isize,counts,displs).eq.irank) then
     935             : c        write(*,*)irank,ikpt
     936           0 :         i_rec = i_rec + 1
     937             : c      if (mod(i_rec-1,isize).eq.irank) then
     938             : 
     939           0 :       allocate ( eigg(neigd) )
     940             : 
     941           0 :           n_start=1
     942           0 :           n_end=neigd
     943             :  
     944           0 :       call cpu_time(t0)
     945             :       ! get current bkpt vector
     946             : 
     947           0 :       zzMat%l_real = l_real
     948           0 :       zzMat%nbasfcn = nbasfcn
     949           0 :       zzMat%nbands = neigd
     950           0 :       IF(l_real) THEN
     951           0 :          ALLOCATE (zzMat%z_r(zzMat%nbasfcn,zzMat%nbands))
     952             :       ELSE
     953           0 :          ALLOCATE (zzMat%z_c(zzMat%nbasfcn,zzMat%nbands))
     954             :       END IF
     955             : 
     956             : !      CALL cdn_read(
     957             : !     >              eig_id,
     958             : !     >              nvd,jspd,irank,isize,kptibz,jspin,nbasfcn, !wannierspin instead of jspd?
     959             : !     >              l_ss,l_noco,neigd,n_start,n_end,
     960             : !     <              nmat,nv,ello,evdu,epar,
     961             : !     <              wk,nbands,eigg,zzMat)
     962             : 
     963           0 :       call cpu_time(t1)
     964           0 :       t_eig = t_eig + t1 - t0
     965             : 
     966           0 :       zMat_b%l_real = zzMat%l_real
     967           0 :       zMat_b2%l_real = zzMat%l_real
     968           0 :       zMat_b%nbasfcn = zzMat%nbasfcn
     969           0 :       zMat_b2%nbasfcn = zzMat%nbasfcn
     970           0 :       zMat_b%nbands = zzMat%nbands
     971           0 :       zMat_b2%nbands = zzMat%nbands
     972             :     !  IF(zzMat%l_real) THEN
     973             :     !     ALLOCATE (zMat_b%z_r(zMat%nbasfcn,zMat%nbands))
     974             :     !     ALLOCATE (zMat_b2%z_r(zMat%nbasfcn,zMat%nbands))
     975             :     !  ELSE
     976             :     !     ALLOCATE (zMat_b%z_c(zMat%nbasfcn,zMat%nbands))
     977             :     !     ALLOCATE (zMat_b2%z_c(zMat%nbasfcn,zMat%nbands))
     978             :     !  END IF
     979             :       
     980             :       
     981           0 :       IF(zzMat%l_real) THEN
     982           0 :          ALLOCATE (zMat_b%z_r(zMat_b%nbasfcn,zMat_b%nbands))
     983           0 :          ALLOCATE (zMat_b2%z_r(zMat_b2%nbasfcn,zMat_b2%nbands))
     984             :       ELSE
     985           0 :          ALLOCATE (zMat_b%z_c(zMat_b%nbasfcn,zMat_b%nbands))
     986           0 :          ALLOCATE (zMat_b2%z_c(zMat_b2%nbasfcn,zMat_b2%nbands))
     987             :       END IF
     988             :       
     989             :       
     990             : 
     991           0 :       allocate (we_b(neigd), we_b2(neigd))
     992             : 
     993             :   !!! the cycle by the nearest neighbors (nntot) for each kpoint
     994             : 
     995           0 :        do 15  ikpt_b = 1,nntot
     996           0 :           kptibz_b=bpt(ikpt_b,ikpt)
     997             : 
     998             :           if(wann%l_bzsym) oper_b=mapkoper(kptibz_b)
     999             :           if (wann%l_bzsym) kptibz_b=irreduc(kptibz_b)
    1000             : 
    1001             :           n_start=1
    1002             :           n_end=neigd
    1003             : 
    1004           0 :         eigg = 0.
    1005           0 :         call cpu_time(t0)
    1006           0 :         call judft_error("BUG not implement")
    1007             : !        CALL cdn_read(
    1008             : !     >                eig_id,
    1009             : !     >                nvd,jspd,irank,isize,kptibz_b,jspin,nbasfcn, !wannierspin instead of jspd?
    1010             : !     >                l_ss,l_noco,neigd,n_start,n_end,
    1011             : !     <                nmat_b,nv_b,ello,evdu,epar,
    1012             : !     <                wk_b,nbands_b,eigg,zzMat)
    1013             : 
    1014           0 :         nslibd_b = 0
    1015             : 
    1016           0 :       IF(zzMat%l_real) THEN
    1017           0 :          zMat_b%z_r = 0.0
    1018             :       ELSE
    1019           0 :          zMat_b%z_c = CMPLX(0.0,0.0)
    1020             :       END IF
    1021             : 
    1022           0 :       eig_b(:) = 0.
    1023             : 
    1024           0 :         do i = 1,nbands_b
    1025             :           if((eigg(i).ge.e1s.and.nslibd_b.lt.numbands
    1026             :      & .and.wann%l_bynumber)
    1027             :      &.or.(eigg(i).ge.e1s.and.eigg(i).le.e2s.and.wann%l_byenergy)
    1028           0 :      &.or.(i.ge.wann%band_min(jspin)
    1029             :      &                 .and.
    1030             :      &       (i.le.wann%band_max(jspin))
    1031             :      &                 .and.
    1032           0 :      &       wann%l_byindex))then
    1033           0 :             nslibd_b = nslibd_b + 1
    1034           0 :             eig_b(nslibd_b) = eigg(i)
    1035           0 :             we_b(nslibd_b) = we_b(i)
    1036           0 :             if(l_noco)then
    1037           0 :                funbas =        nv_b(1) + nlotot
    1038           0 :                funbas = funbas+nv_b(2) + nlotot
    1039             :             else
    1040           0 :                funbas = nv_b(jspin) + nlotot
    1041             :             endif
    1042           0 :             IF (zzMat%l_real) THEN
    1043           0 :                do j = 1,funbas
    1044           0 :                   zMat_b%z_r(j,nslibd_b) = zzMat%z_r(j,i)
    1045             :                enddo
    1046             :             ELSE
    1047           0 :                do j = 1,funbas
    1048           0 :                   zMat_b%z_c(j,nslibd_b) = zzMat%z_c(j,i)
    1049             :                enddo
    1050             :             END IF
    1051             :           endif
    1052             :         enddo
    1053             : 
    1054             : c***********************************************************
    1055             : c              Rotate the wavefunction of next neighbor.
    1056             : c***********************************************************
    1057             :       if (wann%l_bzsym .and. (oper_b.ne.1)  ) then
    1058             : !         TODO
    1059             : !         call wann_kptsrotate(
    1060             : !     >            natd,nlod,llod,
    1061             : !     >            ntypd,nlo,llo,invsat,
    1062             : !     >            l_noco,l_soc,
    1063             : !     >            ntype,neq,nlotot,
    1064             : !     >            kveclo_b,jspin,
    1065             : !     >            oper_b,nop,mrot,nvd,
    1066             : !     >            nv_b,
    1067             : !     >            shiftkpt(:,bpt(ikpt_b,ikpt)),
    1068             : !     >            tau,
    1069             : !     x            bkpt_b,k1_b(:,:),
    1070             : !     x            k2_b(:,:),k3_b(:,:),
    1071             : !     x            zMat_b,nsfactor_b)
    1072             :       else
    1073             :          nsfactor_b=cmplx(1.0,0.0)
    1074             :       endif
    1075           0 :         noccbd_b = nslibd_b
    1076           0 :       call cpu_time(t1)
    1077           0 :       t_eig = t_eig + t1 - t0
    1078             : 
    1079           0 :       addnoco = 0
    1080           0 :       if(l_noco.and.jspin.eq.2) addnoco=nv_b(1)+nlotot
    1081             : 
    1082             :       ! set up a(k+b1),b(k+b1),c(k+b1)
    1083           0 :       allocate( acof_b(noccbd_b,0:lmd,natd),
    1084           0 :      >          bcof_b(noccbd_b,0:lmd,natd),
    1085           0 :      >          ccof_b(-llod:llod,noccbd_b,nlod,natd) )
    1086             : 
    1087           0 :       call cpu_time(t0)
    1088             : 
    1089           0 :       ALLOCATE(lapw_b%k1(SIZE(k1_b,1),SIZE(k1_b,2)))
    1090           0 :       ALLOCATE(lapw_b%k2(SIZE(k1_b,1),SIZE(k1_b,2)))
    1091           0 :       ALLOCATE(lapw_b%k3(SIZE(k1_b,1),SIZE(k1_b,2)))
    1092           0 :       lapw_b%k1 = k1_b
    1093           0 :       lapw_b%k2 = k2_b
    1094           0 :       lapw_b%k3 = k3_b
    1095             :       lapw_b%nmat = nmat_b
    1096           0 :       lapw_b%nv = nv_b
    1097             :       ! I think the other variables of lapw are not needed here.
    1098             : !TODO
    1099             : !      CALL abcof(input,atoms,noccbd_b,sym,cell,bkpt_b,lapw_b,
    1100             : !     +           noccbd_b,usdus,noco,jspin,kveclo_b,oneD,
    1101             : !     +           acof_b,bcof_b,ccof_b,zMat_b)
    1102             : 
    1103           0 :       DEALLOCATE(lapw_b%k1,lapw_b%k2,lapw_b%k3)
    1104             : 
    1105             :       call wann_abinv(atoms,
    1106           0 :      X        acof_b,bcof_b,ccof_b)
    1107           0 :       call cpu_time(t1)
    1108           0 :       t_abcof = t_abcof + t1 - t0
    1109             : 
    1110           0 :        do 25  ikpt_b2 = 1,nntot!ikpt_b!nntot
    1111             : 
    1112           0 :           kptibz_b2=bpt(ikpt_b2,ikpt)
    1113             :           if(wann%l_bzsym) oper_b2=mapkoper(kptibz_b2)
    1114             :           if (wann%l_bzsym) kptibz_b2=irreduc(kptibz_b2)
    1115             : 
    1116             :           n_start=1
    1117             :           n_end=neigd
    1118             : 
    1119           0 :           eigg = 0.
    1120           0 :           call cpu_time(t0)
    1121             : 
    1122           0 :           WRITE(*,*) 'Here probably the wrong record is read in'
    1123           0 :           WRITE(*,*) 'Should eig_id not be dependent on iqpt_b?'
    1124           0 :           WRITE(*,*) '(in wann_uHu)'
    1125             : !          CALL cdn_read(
    1126             : !     >                  eig_id,
    1127             : !     >                  nvd,jspd,irank,isize,kptibz_b2,jspin,nbasfcn, !wannierspin instead of jspd?
    1128             : !     >                  l_ss,l_noco,neigd,n_start,n_end,
    1129             : !     <                  nmat_b2,nv_b2,ello,evdu,epar,kveclo_b2,
    1130             : !     <                  k1_b2,k2_b2,k3_b2,bkpt_b2,wk_b2,nbands_b2,
    1131             : !     <                  eigg,zzMat)
    1132             : 
    1133           0 :         nslibd_b2 = 0
    1134             : 
    1135           0 :       IF(zzMat%l_real) THEN
    1136           0 :          zMat_b2%z_r = 0.0
    1137             :       ELSE
    1138           0 :          zMat_b2%z_c = CMPLX(0.0,0.0)
    1139             :       END IF
    1140             : 
    1141           0 :       eig_b2(:) = 0.
    1142             : 
    1143           0 :         do i = 1,nbands_b2
    1144             :           if((eigg(i).ge.e1s.and.nslibd_b2.lt.numbands
    1145             :      & .and.wann%l_bynumber)
    1146             :      &.or.(eigg(i).ge.e1s.and.eigg(i).le.e2s.and.wann%l_byenergy)
    1147           0 :      &.or.(i.ge.wann%band_min(jspin_b)
    1148             :      &                 .and.
    1149             :      &       (i.le.wann%band_max(jspin_b))
    1150             :      &                 .and.
    1151           0 :      &       wann%l_byindex))then
    1152           0 :             nslibd_b2 = nslibd_b2 + 1
    1153           0 :             eig_b2(nslibd_b2) = eigg(i)
    1154           0 :             we_b2(nslibd_b2) = we_b2(i)
    1155           0 :             if(l_noco)then
    1156           0 :                funbas =        nv_b2(1) + nlotot
    1157           0 :                funbas = funbas+nv_b2(2) + nlotot
    1158             :             else
    1159           0 :                funbas = nv_b2(jspin_b) + nlotot
    1160             :             endif
    1161           0 :             IF (zzMat%l_real) THEN
    1162           0 :                do j = 1,funbas
    1163           0 :                   zMat_b2%z_r(j,nslibd_b2) = zzMat%z_r(j,i)
    1164             :                enddo
    1165             :             ELSE
    1166           0 :                do j = 1,funbas
    1167           0 :                   zMat_b2%z_c(j,nslibd_b2) = zzMat%z_c(j,i)
    1168             :                enddo
    1169             :             END IF
    1170             :           endif
    1171             :         enddo
    1172             : 
    1173             : c***********************************************************
    1174             : c              Rotate the wavefunction of next neighbor.
    1175             : c***********************************************************
    1176             :       if (wann%l_bzsym .and. (oper_b2.ne.1)  ) then
    1177             : !  TODO
    1178             : !         call wann_kptsrotate(
    1179             : !     >            natd,nlod,llod,
    1180             : !     >            ntypd,nlo,llo,invsat,
    1181             : !     >            l_noco,l_soc,
    1182             : !     >            ntype,neq,nlotot,
    1183             : !     >            kveclo_b2,jspin_b,
    1184             : !     >            oper_b2,nop,mrot,nvd,
    1185             : !     >            nv_b2,
    1186             : !     >            shiftkpt(:,bpt(ikpt_b2,ikpt)),
    1187             : !     >            tau,
    1188             : !     x            bkpt_b2,k1_b2(:,:),
    1189             : !     x            k2_b2(:,:),k3_b2(:,:),
    1190             : !     x            zMat_b2,nsfactor_b2)
    1191             :       else
    1192             :          nsfactor_b2=cmplx(1.0,0.0)
    1193             :       endif
    1194           0 :         noccbd_b2 = nslibd_b2
    1195           0 :       call cpu_time(t1)
    1196           0 :       t_eig = t_eig + t1 - t0
    1197             : 
    1198           0 :       addnoco2 = 0
    1199           0 :       if(l_noco.and.jspin_b.eq.2) addnoco2=nv_b2(1)+nlotot
    1200             : 
    1201             :       ! set up a(k+b2),b(k+b2),c(k+b2)
    1202           0 :       allocate( acof_b2(noccbd_b2,0:lmd,natd),
    1203           0 :      >          bcof_b2(noccbd_b2,0:lmd,natd),
    1204           0 :      >          ccof_b2(-llod:llod,noccbd_b2,nlod,natd) )
    1205             : 
    1206           0 :       call cpu_time(t0)
    1207             : 
    1208             : 
    1209           0 :       ALLOCATE(lapw_b2%k1(SIZE(k1_b2,1),SIZE(k1_b2,2)))
    1210           0 :       ALLOCATE(lapw_b2%k2(SIZE(k1_b2,1),SIZE(k1_b2,2)))
    1211           0 :       ALLOCATE(lapw_b2%k3(SIZE(k1_b2,1),SIZE(k1_b2,2)))
    1212           0 :       lapw_b2%k1 = k1_b2
    1213           0 :       lapw_b2%k2 = k2_b2
    1214           0 :       lapw_b2%k3 = k3_b2
    1215             :       lapw_b2%nmat = nmat_b2
    1216           0 :       lapw_b2%nv = nv_b2
    1217             :       ! I think the other variables of lapw are not needed here.
    1218             : !TODO
    1219             : !      CALL abcof(input,atoms,noccbd_b2,sym,cell,bkpt_b2,lapw_b2,
    1220             : !     +           noccbd_b2,usdus,noco,jspin_b,kveclo_b2,oneD,
    1221             : !     +           acof_b2,bcof_b2,ccof_b2,zMat_b2)
    1222             : 
    1223           0 :       DEALLOCATE(lapw_b2%k1,lapw_b2%k2,lapw_b2%k3)
    1224             : 
    1225             :       call wann_abinv(atoms,
    1226           0 :      X        acof_b2,bcof_b2,ccof_b2)
    1227           0 :       call cpu_time(t1)
    1228           0 :       t_abcof = t_abcof + t1 - t0
    1229             : 
    1230             : 
    1231             : c**************************************c
    1232             : c    calculate uHu matrix due to:      c
    1233             : c             (i)    interstitial      c
    1234             : c             (ii)   muffin tins       c
    1235             : c             (iii)  vacuum            c 
    1236             : c**************************************c
    1237             : c      if(ANY(gb(3:3,ikpt_b,ikpt).ne.0) .or.
    1238             : c     >   ANY(gb(3:3,ikpt_b2,ikpt).ne.0)) then
    1239             : c       write(*,*)ikpt,ikpt_b,ikpt_b2
    1240             : c       write(*,*)gb(:,ikpt_b,ikpt)
    1241             : c       write(*,*)gb(:,ikpt_b2,ikpt)
    1242             : c       uHu(:,:,ikpt_b2,ikpt_b,i_rec) = cmplx(0.,0.)
    1243             : c       goto 848
    1244             : c      endif
    1245             : 
    1246             :       ! MT (SPH, NON)
    1247           0 :       if(.not.(l_skip_sph.and.l_skip_non)) then
    1248           0 :          call cpu_time(t0)
    1249             :          call wann_uHu_sph(
    1250             :      >        cmplx(1.,0.),nbnd,llod,nslibd_b,nslibd_b2,nlod,natd,ntypd,
    1251             :      >        lmd,jmtd,taual,nop,lmax,ntype,neq,nlo,llo,
    1252             :      >        acof_b,bcof_b,ccof_b,bkpt_b2,
    1253             :      >        acof_b2,bcof_b2,ccof_b2,bkpt_b,bkpt,
    1254             :      >        gb(:,ikpt_b,ikpt),gb(:,ikpt_b2,ikpt),
    1255             :      <        tuu(:,:,:,:,tspin2),tud(:,:,:,:,tspin2),
    1256             :      >        tdu(:,:,:,:,tspin2),tdd(:,:,:,:,tspin2),
    1257             :      >        tuulo(:,:,:,:,:,tspin2),tulou(:,:,:,:,:,tspin2),
    1258             :      >        tdulo(:,:,:,:,:,tspin2),tulod(:,:,:,:,:,tspin2),
    1259             :      >        tuloulo(:,:,:,:,:,:,tspin2),
    1260             :      >        kdiff,kdiff,nntot,nntot,
    1261           0 :      >        uHu(:,:,ikpt_b2,ikpt_b,i_rec))
    1262           0 :          call cpu_time(t1)
    1263           0 :          t_sph = t_sph + t1 - t0
    1264             :       endif
    1265             : c      if(ikpt.eq.1) then
    1266             : c       write(*,*)'SPH,NON',ikpt_b2,ikpt_b
    1267             : c       write(*,*)uHu(1,3,ikpt_b2,ikpt_b,i_rec)
    1268             : c       write(*,*)uHu(3,1,ikpt_b2,ikpt_b,i_rec)
    1269             : c       uHuold(:,:) = uHu(:,:,ikpt_b2,ikpt_b,i_rec)
    1270             : c      endif
    1271             : 
    1272             :       ! MT (SOC)
    1273           0 :       if(l_soc.and.(.not.l_skip_soc)) then
    1274           0 :          call cpu_time(t0)
    1275             :          call wann_uHu_sph(
    1276             :      >        cmplx(1.,0.),nbnd,llod,nslibd_b,nslibd_b2,nlod,natd,ntypd,
    1277             :      >        lmd,jmtd,taual,nop,lmax,ntype,neq,nlo,llo,
    1278             :      >        acof_b,bcof_b,ccof_b,bkpt_b2,
    1279             :      >        acof_b2,bcof_b2,ccof_b2,bkpt_b,bkpt,
    1280             :      >        gb(:,ikpt_b,ikpt),gb(:,ikpt_b2,ikpt),
    1281             :      <        tuu_soc,tud_soc,
    1282             :      >        tdu_soc,tdd_soc,
    1283             :      >        tuulo_soc,tulou_soc,
    1284             :      >        tdulo_soc,tulod_soc,
    1285             :      >        tuloulo_soc,
    1286             :      >        kdiff,kdiff,nntot,nntot,
    1287           0 :      >        uHu(:,:,ikpt_b2,ikpt_b,i_rec))
    1288           0 :          call cpu_time(t1)
    1289           0 :          t_sph = t_sph + t1 - t0
    1290             : c      if(ikpt.eq.1) then
    1291             : c       write(*,*)'SOC',ikpt_b2,ikpt_b
    1292             : c       write(*,*)uHu(1,3,ikpt_b2,ikpt_b,i_rec)-uHuold(1,3)
    1293             : c       write(*,*)uHu(3,1,ikpt_b2,ikpt_b,i_rec)-uHuold(3,1)
    1294             : c       uHuold(:,:) = uHu(:,:,ikpt_b2,ikpt_b,i_rec)
    1295             : c      endif
    1296             :       endif
    1297             : 
    1298             : 
    1299             :       ! In collinear calculation, the potential in interstital and
    1300             :       ! vacuum region is always diagonal with respect to the spin
    1301           0 :       IF( l_noco .OR. ((.NOT. l_noco) .AND. (doublespin.LT.3)) ) THEN
    1302             : 
    1303           0 :          jspin3 = doublespin
    1304           0 :          if(jspins.eq.1) jspin3=1
    1305             : 
    1306           0 :          jspin4 = jspin
    1307           0 :          jspin4_b=jspin_b
    1308           0 :          if(jspins.eq.1) then
    1309           0 :           jspin4=1
    1310           0 :           jspin4_b=1
    1311             :          endif
    1312             : 
    1313           0 :          sign2 = 1
    1314           0 :          if(doublespin.EQ.4) sign2=-1
    1315             : 
    1316             :          ! INT
    1317           0 :          if(.not.l_skip_int) then
    1318           0 :             call cpu_time(t0)
    1319             :             call wann_uHu_int(cmplx(1.,0.),nvd,k1d,k2d,k3d,n3d,
    1320             :      >            nv_b(jspin),nv_b2(jspin_b),nbnd,neigd,
    1321             :      >            nslibd_b,nslibd_b2,nbasfcn,addnoco,addnoco2,
    1322             :      >            k1_b(:,jspin), k2_b(:,jspin), k3_b(:,jspin),
    1323             :      >            gb(:,ikpt_b,ikpt),
    1324             :      >            k1_b2(:,jspin_b),k2_b2(:,jspin_b),k3_b2(:,jspin_b),
    1325             :      >            gb(:,ikpt_b2,ikpt),
    1326             :      >            bkpt,bbmat,vpw(:,jspin3),zMat_b,zMat_b2,rgphs,
    1327             :      >            ustep,ig,jspin.eq.jspin_b,sign2,
    1328           0 :      >            uHu(:,:,ikpt_b2,ikpt_b,i_rec))
    1329           0 :             call cpu_time(t1)
    1330           0 :             t_int = t_int + t1 - t0
    1331             :          endif
    1332             : c      if(ikpt.eq.1) then
    1333             : c       write(*,*)'INT',ikpt_b2,ikpt_b
    1334             : c       write(*,*)uHu(1,3,ikpt_b2,ikpt_b,i_rec)-uHuold(1,3)
    1335             : c       write(*,*)uHu(3,1,ikpt_b2,ikpt_b,i_rec)-uHuold(3,1)
    1336             : c       uHuold(:,:) = uHu(:,:,ikpt_b2,ikpt_b,i_rec)
    1337             : c      endif
    1338             : 
    1339             : 
    1340             :          ! VAC
    1341           0 :          if ((.not.l_skip_vac) .and. film .and. (.not.odi%d1)) then
    1342             : 
    1343           0 :             call cpu_time(t0)
    1344             :             call wann_uHu_vac(
    1345             :      >            cmplx(1.,0.),l_noco,l_soc,zrfs,jspins,nlotot,qpt_i,
    1346             :      >            nbnd,z1,nmzxyd,nmzd,n2d,nv2d,k1d,k2d,k3d,n3d,nvac,ig,
    1347             :      >            rgphs,nmzxy,nmz,delz,ig2,nq2,kv2,area,bmat,bbmat,
    1348             :      >            evac(:,jspin4),evac(:,jspin4_b),bkpt_b,bkpt_b2,
    1349             :      >            vzxy(:,:,:,jspin3),vz,nslibd_b,nslibd_b2,
    1350             :      >            jspin,jspin_b,doublespin,k1_b,k2_b,k3_b,
    1351             :      >            k1_b2,k2_b2,k3_b2,wannierspin,nvd,nbasfcn,neigd,
    1352             :      >            zMat_b,zMat_b2,nv_b,nv_b2,omtil,gb(:,ikpt_b,ikpt),
    1353             :      >            gb(:,ikpt_b2,ikpt),sign2,
    1354           0 :      >            uHu(:,:,ikpt_b2,ikpt_b,i_rec))
    1355           0 :             call cpu_time(t1)
    1356           0 :             t_vac = t_vac + t1 - t0
    1357             : 
    1358           0 :          elseif ((.not.l_skip_vac) .and. odi%d1) then
    1359             : 
    1360           0 :             call cpu_time(t0)
    1361             :             call wann_uHu_od_vac(
    1362             :      >            DIMENSION,oneD,vacuum,stars,cell,
    1363             :      >            cmplx(1.,0.),l_noco,l_soc,jspins,nlotot,nbnd,z1,
    1364             :      >            nmzxyd,nmzd,nv2d,k1d,k2d,k3d,n2d,n3d,ig,nmzxy,nmz,
    1365             :      >            delz,ig2,bbmat,evac(1,jspin4),evac(1,jspin4_b),bkpt_b,
    1366             :      >            bkpt_b2,odi,vzxy(:,:,:,jspin3),vz,nslibd_b,nslibd_b2,
    1367             :      >            jspin,jspin_b,doublespin,k1_b,k2_b,k3_b,k1_b2,k2_b2,
    1368             :      >          k3_b2,wannierspin,nvd,area,nbasfcn,neigd,zMat_b,zMat_b2,
    1369             :      >            nv_b,nv_b2,sk2,phi2,omtil,gb(:,ikpt_b,ikpt),
    1370             :      >            gb(:,ikpt_b2,ikpt),qpt_i,sign2,
    1371           0 :      >            uHu(:,:,ikpt_b2,ikpt_b,i_rec))
    1372           0 :             call cpu_time(t1)
    1373           0 :             t_vac = t_vac + t1 - t0
    1374             : 
    1375             :          endif
    1376             : c      if(ikpt.eq.1) then
    1377             : c       write(*,*)'VAC',ikpt_b2,ikpt_b
    1378             : c       write(*,*)uHu(1,3,ikpt_b2,ikpt_b,i_rec)-uHuold(1,3)
    1379             : c       write(*,*)uHu(3,1,ikpt_b2,ikpt_b,i_rec)-uHuold(3,1)
    1380             : c       uHuold(:,:) = uHu(:,:,ikpt_b2,ikpt_b,i_rec)
    1381             : c      endif
    1382             :       ELSE
    1383           0 :        IF(ikpt.EQ.1 .AND. ikpt_b.EQ.1 .AND. ikpt_b2.EQ.1) THEN
    1384           0 :           if(l_p0) write(*,*)'skip INT and VAC'
    1385             :        ENDIF
    1386             :       ENDIF
    1387             : 
    1388             : 
    1389             : 848     continue
    1390             : 
    1391           0 :         deallocate (acof_b2,bcof_b2,ccof_b2)
    1392           0 : 25      continue ! end of loop by the nearest k-neighbors
    1393             : 
    1394           0 :         deallocate (acof_b,bcof_b,ccof_b)
    1395           0 : 15      continue ! end of loop by the nearest k-neighbors
    1396           0 :         deallocate ( we_b,we_b2 )
    1397           0 :       deallocate ( eigg )
    1398             : 
    1399             : 
    1400             : 
    1401           0 :       IF(l_real) THEN
    1402           0 :          DEALLOCATE ( zzMat%z_r )
    1403             :       ELSE
    1404           0 :          DEALLOCATE ( zzMat%z_c )
    1405             :       END IF
    1406             : 
    1407             : 
    1408             : 
    1409           0 :       IF(zzMat%l_real) THEN
    1410           0 :          DEALLOCATE (zMat_b%z_r)
    1411           0 :          DEALLOCATE (zMat_b2%z_r)
    1412             :       ELSE
    1413           0 :          DEALLOCATE (zMat_b%z_c)
    1414           0 :          DEALLOCATE (zMat_b2%z_c)
    1415             :       END IF
    1416             : 
    1417             : 
    1418             : 
    1419             :       endif   ! loop by processors
    1420             : 
    1421             : 
    1422             : 
    1423             : 
    1424             : 
    1425             : 
    1426           0 : 10    continue ! end of cycle by the k-points
    1427             : 
    1428             : #ifdef CPP_MPI
    1429           0 :       call MPI_BARRIER(MPI_COMM_WORLD,ierr)
    1430             : #endif
    1431             :   5   continue
    1432             : 
    1433           0 :       l_symcheck = (.not.l_nocosoc).or.(doublespin.eq.doublespin_max)
    1434           0 :       jspin4=jspin
    1435           0 :       if(l_nocosoc) jspin4=1
    1436           0 :       if(l_p0) write(*,*)'write uHu file'
    1437           0 :       uHu = uHu * htr2ev
    1438             :       call wann_write_uHu(jspin4,l_p0,fullnkpts,nntot,nntot,wann,
    1439             :      >          nbnd,bpt,gb,isize,irank,fending,'',uHu,
    1440             :      >          counts(irank),counts,displs,isize,
    1441           0 :      >          wann%l_unformatted,.true.,l_symcheck)
    1442           0 :       uHu = uHu / htr2ev
    1443           0 :       if(.not.l_nocosoc) deallocate(uHu)
    1444             : 
    1445           0 :       if(.not.l_noco)nrec=nrec+nkpts
    1446             : 
    1447           0 : 110   continue ! end of cycle by spins
    1448             : 
    1449             : #ifdef CPP_MPI
    1450           0 :       call MPI_BARRIER(MPI_COMM_WORLD,ierr)
    1451             : #endif
    1452             : 
    1453             :       ! close eig files
    1454             :       IF (l_gwf) THEN
    1455             : 
    1456             : !         CALL close_eig(eig_id)
    1457             : 
    1458             : !         CLOSE(66)
    1459             :          DO iqpt_b=1,nntot_q
    1460             : !            CALL close_eig(innerEig_idList(iqpt_b))
    1461             :          ENDDO
    1462             :       ENDIF
    1463             : 
    1464             : 
    1465             : c      if(l_p0) write(*,*)'write uHu file'
    1466             : c      uHu = htr2ev * uHu
    1467             : c      call wann_write_uHu(1,l_p0,fullnkpts,nntot,nntot,wann,
    1468             : c     >          nbnd,bpt,gb,isize,irank,fending,'',uHu,
    1469             : c     >          counts(irank),counts,displs,isize,
    1470             : c     >          wann%l_unformatted,.true.)
    1471           0 :       if(allocated(uHu)) deallocate(uHu)
    1472             : c      if(allocated(uHuold)) deallocate(uHuold)
    1473             : 
    1474           0 :       DEALLOCATE(innerEig_idList)
    1475             : 
    1476           0 : 314   continue ! iqpt, q-points
    1477             : c************************************************c
    1478             : c               END Q LOOP                       c
    1479             : c************************************************c
    1480             : 
    1481             : 
    1482           0 :       deallocate ( kveclo,nv,k1,k2,k3 )
    1483           0 :       deallocate(flo)
    1484           0 :       deallocate ( ff,gg)
    1485           0 :       if (allocated(nv_b))deallocate(kveclo_b,nv_b,k1_b,k2_b,k3_b)
    1486           0 :       if (allocated(nv_b2))deallocate(kveclo_b2,nv_b2,k1_b2,k2_b2,k3_b2)
    1487           0 :       if (wann%l_bzsym)deallocate(irreduc,mapkoper,shiftkpt)
    1488           0 :       if (wann%l_bzsym.AND.l_gwf)deallocate(irreduc_q,mapqoper,shiftqpt)
    1489           0 :       if (allocated(pair_to_do)) deallocate(pair_to_do,maptopair)
    1490           0 :       if (allocated(pair_to_do_q)) deallocate(pair_to_do_q,maptopair_q)
    1491           0 :       if (allocated(kdiff)) deallocate ( kdiff )
    1492           0 :       if (allocated(qdiff)) deallocate(qdiff,zero_qdiff)
    1493           0 :       if(allocated(vpw)) deallocate(vpw)
    1494           0 :       if(allocated(vzxy)) deallocate(vzxy)
    1495           0 :       if(allocated(vz)) deallocate(vz)
    1496           0 :       deallocate(vr)
    1497           0 :       deallocate(tdd,tdu,tud,tuu)
    1498           0 :       deallocate(tdulo,tuulo)
    1499           0 :       deallocate(tulou,tulod)
    1500           0 :       deallocate(tuloulo)
    1501           0 :       deallocate(tdd_soc,tdu_soc,tud_soc,tuu_soc)
    1502           0 :       deallocate(tdulo_soc,tuulo_soc)
    1503           0 :       deallocate(tulou_soc,tulod_soc)
    1504           0 :       deallocate(tuloulo_soc)
    1505           0 :       deallocate(counts,displs)   
    1506             : 
    1507             : 9110  continue
    1508             : 
    1509           0 :       if(l_sgwf .or. l_socgwf) gb_q = gb_q/2
    1510             : 
    1511           0 :       if(l_p0.and.l_gwf) then
    1512             :          call wann_uHu_commat(
    1513             :      >              fullnkpts,nntot,bpt,fullnqpts,nntot_q,bpt_q,
    1514             :      >              gb,gb_q,aux_latt_const,wann%l_unformatted,l_dim,
    1515           0 :      >              nparampts,param_vec/2.0)
    1516             :       endif
    1517             : 
    1518           0 :       if(allocated(gb)) deallocate(gb,bpt)
    1519           0 :       if(allocated(gb_q)) deallocate(gb_q,bpt_q)
    1520             : 
    1521             : #ifdef CPP_MPI
    1522           0 :       call MPI_BARRIER(MPI_COMM_WORLD,ierr)
    1523             : #endif
    1524           0 :       call cpu_time(t1)
    1525           0 :       t_total = t1-t00
    1526           0 :       if(l_p0) then
    1527           0 :        write(*,900)'t_init  =',t_init
    1528           0 :        write(*,900)'t_myTlmplm=',t_myTlmplm
    1529           0 :        write(*,900)'t_eig   =',t_eig   
    1530           0 :        write(*,900)'t_abcof =',t_abcof 
    1531           0 :        write(*,900)'t_int   =',t_int   
    1532           0 :        write(*,900)'t_sph   =',t_sph   
    1533           0 :        write(*,900)'t_vac   =',t_vac   
    1534           0 :        write(*,900)'t_total =',t_total
    1535             :       endif
    1536             : 
    1537             : 900   FORMAT(a,f14.4)
    1538             : 
    1539           0 :       END SUBROUTINE wann_uHu
    1540             : 
    1541             :  
    1542             :       END MODULE m_wann_uHu
    1543             : 
    1544             : 

Generated by: LCOV version 1.13