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

Generated by: LCOV version 1.13