LCOV - code coverage report
Current view: top level - wannier - wann_updown.F (source / functions) Coverage Total Hit
Test: FLEUR test coverage Lines: 0.0 % 488 0
Test Date: 2025-06-14 04:34:23 Functions: 0.0 % 1 0

            Line data    Source code
       1              : !--------------------------------------------------------------------------------
       2              : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
       3              : ! This file is part of FLEUR and available as free software under the conditions
       4              : ! of the MIT license as expressed in the LICENSE file in more detail.
       5              : !--------------------------------------------------------------------------------
       6              : 
       7              :       MODULE m_wann_updown
       8              :       use m_juDFT
       9              : #ifdef CPP_MPI
      10              :       use mpi 
      11              : #endif
      12              :       CONTAINS
      13            0 :       SUBROUTINE wann_updown(
      14              :      >      wann,fmpi,input,kpts,sym,atoms,stars,vacuum,
      15              :      >      sphhar ,noco,nococonv,cell,vTot,
      16              :      >      enpara,
      17              :      >      eig_id,l_real,mpi_communicatior,l_dulo,l_noco,l_ss,lmaxd,
      18              :      >      ntypd,
      19              :      >      neigd,natd,nop,nvd,jspd,llod,nlod,ntype,
      20            0 :      >      omtil,nlo,llo,lapw_l,invtab,mrot,ngopr,neq,lmax,
      21            0 :      >      invsat,invsatnr,nkpt,taual,rmt,amat,bmat,bbmat,alph,
      22            0 :      >      beta,qss,sk2,phi2,irank,isize,n3d,nmzxyd,nmzd,
      23            0 :      >      jmtd,nlhd,nq3,nvac,invs,invs2,film,nlh,jri,ntypsd,
      24            0 :      >      ntypsy,jspins,nkptd,dx,n2d,rmsh,e1s,e2s,ulo_der,
      25            0 :      >      ustep,ig,k1d,k2d,k3d,rgphs,slice,kk,nnne,
      26            0 :      >      z1,nv2d,nmzxy,nmz,delz,ig2,area,tau,zatom,nq2,nop2,
      27            0 :      >      volint,symor,pos,ef,l_soc,
      28            0 :      >      memd,lnonsph,clnu,lmplmd,mlh,nmem,llh,lo1l,
      29              :      >      theta,phi)
      30              : c****************************************************************************
      31              : c     Calculate spin-up-spin-down matrix elements.
      32              : c     Useful in the case of spin-orbit coupling and noncollinear magnetism.
      33              : c     Implemented matrix elements:
      34              : c         (i)   Overlaps => Pauli matrix elements.
      35              : c         (ii)  Momentum operator => Spin current.
      36              : c         (iii) H_{SOC} => Include SOC in the basis set of Wannier functions.
      37              : c         (iv)  Spin-current into the MT-spheres => Torque.
      38              : c         (v)   Matrix elements of exchange field => Torque.
      39              : c     Frank Freimuth
      40              : c****************************************************************************
      41              :       use m_types
      42              :       use m_juDFT
      43              :       use m_wann_rw_eig
      44              :       use m_abcof
      45              :       use m_radfun
      46              :       use m_radflo
      47              :       use m_cdnread, only : cdn_read0, cdn_read
      48              :       use m_constants
      49              : !      use m_wann_mmk0_od_vac
      50              :       use m_wann_mmk0_vac
      51              :       use m_wann_mmk0_updown_sph
      52              :       use m_wann_mmk0_updown_sph_at
      53              :       use m_wann_abinv
      54              :       use m_wann_kptsrotate
      55              :       use m_wann_read_inp
      56              :       use m_wann_radovlp_integrals
      57              :       use m_wann_socmat
      58              :       use m_wann_socmat_vec
      59              :       use m_wann_write_matrix5
      60              :       use m_wann_write_nabla
      61              :       use m_wann_write_matrix6
      62              :       use m_wann_write_amn
      63              :       use m_wann_mmkb_int
      64              :       use m_wann_perpmag
      65              :       use m_wann_gabeffgagaunt
      66              :       use m_wann_torque
      67              : #ifdef CPP_TOPO
      68              :       use m_wann_nabla_updown
      69              :       use m_wann_surfcurr_updown
      70              :       use m_wann_surfcurr_int_updown
      71              : #endif
      72              :       IMPLICIT NONE
      73              : 
      74              : #ifdef CPP_MPI
      75              :       integer ierr
      76              :       integer cpu_index
      77              :       integer stt(MPI_STATUS_SIZE)
      78              : 
      79              : #endif
      80              : 
      81              : 
      82              :       TYPE(t_wann),      INTENT(INOUT) :: wann
      83              :       TYPE(t_mpi),INTENT(IN)    :: fmpi
      84              :       TYPE(t_input),INTENT(IN)  :: input
      85              :       TYPE(t_kpts),      INTENT(IN) :: kpts
      86              :       TYPE(t_sym),INTENT(IN)    :: sym
      87              :       TYPE(t_atoms),INTENT(IN)  :: atoms
      88              :       TYPE(t_stars),INTENT(IN)  :: stars
      89              :       TYPE(t_vacuum),INTENT(IN) :: vacuum
      90              :       TYPE(t_sphhar),INTENT(IN) :: sphhar
      91              :        
      92              :       TYPE(t_noco),INTENT(IN)   :: noco
      93              :       TYPE(t_nococonv),INTENT(IN):: nococonv
      94              :       TYPE(t_cell),INTENT(IN)   :: cell
      95              :       TYPE(t_potden),INTENT(IN) :: vTot
      96              :       TYPE(t_enpara),INTENT(IN) :: enpara
      97              : 
      98              :       logical, intent (in) :: invs,invs2,film,slice,symor,l_real
      99              :       integer, intent (in) :: lmaxd,ntypd,neigd,nkptd,kk,nnne
     100              :       integer, intent (in) :: mpi_communicatior
     101              :       integer, intent (in) :: natd,nop,nvd,jspd,nq2,nop2,eig_id
     102              :       integer, intent (in) :: llod,nlod,ntype,n3d,n2d
     103              :       integer, intent (in) :: nmzxyd,nmzd,jmtd,nlhd,nq3,nvac
     104              :       integer, intent (in) :: ntypsd,jspins,k1d,k2d,k3d
     105              :       real,    intent (in) :: omtil,e1s,e2s,delz,area,z1,volint
     106              :       integer, intent (in) :: ig(-k1d:k1d,-k2d:k2d,-k3d:k3d)
     107              :       complex, intent (in) :: rgphs(-k1d:k1d,-k2d:k2d,-k3d:k3d)
     108              :       integer, intent (in) :: nlh(ntypsd),jri(ntypd),ntypsy(natd)
     109              :       integer, intent (in) :: nlo(ntypd),llo(nlod,ntypd),lapw_l(ntypd)
     110              :       integer, intent (in) :: invtab(nop),mrot(3,3,nop),ngopr(natd)
     111              :       integer, intent (in) :: neq(ntypd),lmax(ntypd)
     112              :       integer, intent (in) :: invsat(natd),invsatnr(natd),nkpt
     113              :       integer, intent (in) :: irank,isize,nv2d,nmzxy,nmz
     114              :       integer, intent (in) :: ulo_der(nlod,ntypd),ig2(n3d)
     115              :       real,    intent (in) :: taual(3,natd),rmt(ntypd),dx(ntypd)
     116              :       real,    intent (in) :: amat(3,3),bmat(3,3),bbmat(3,3)
     117              :       real,    intent (in) :: rmsh(jmtd,ntypd),tau(3,nop),zatom(ntype)
     118              :       real,    intent (in) :: alph(ntypd),beta(ntypd),qss(3)
     119              :       real,    intent (in) :: pos(3,natd)
     120              :       real, intent(in)     :: ef
     121              :       complex, intent (in) :: ustep(n3d)
     122              :       logical, intent (in) :: l_dulo(nlod,ntypd),l_noco,l_ss
     123              :       logical, intent (in) :: l_soc
     124              :       integer, intent (in) :: memd
     125              :       integer, intent (in) :: lnonsph(ntypd)
     126              :       complex, intent (in) :: clnu(memd,0:nlhd,ntypsd)
     127              :       integer, intent (in) :: lmplmd
     128              :       integer, intent (in) :: mlh(memd,0:nlhd,ntypsd)
     129              :       integer, intent (in) :: nmem(0:nlhd,ntypsd)
     130              :       integer, intent (in) :: llh(0:nlhd,ntypsd)
     131              :       integer, intent (in) :: lo1l(0:llod,ntypd)
     132              :       real,    intent (in) :: sk2(n2d),phi2(n2d)
     133              :       real, intent(in) :: theta,phi
     134              : 
     135              : cccccccccccccccccc   local variables   cccccccccccccccccccc
     136              : !      type(t_wann) :: wann  !now input variable
     137              :       integer :: lmd,nlotot,n,nmat,iter,ikpt
     138              :       integer :: addnoco,addnoco2,funbas,loplod,at
     139              :       integer :: noccbd,noccbd_b,nn,nkpts,i,jspin,j,l,i_rec,m,nwf,nwfp
     140              :       integer :: jsp_start,jsp_end,nrec,nrec1,nbands,nbands_b
     141              :       integer :: nodeu,noded,n_size,na,n_rank,nbnd,numbands
     142              :       integer :: i1,i2,i3,in,lda
     143            0 :       integer :: n_bands(0:neigd),nslibd,nslibd_b
     144              :       character(len=8) :: dop,iop,name(10)
     145              :       real    :: sfp,tpi,wronk,wk,wk_b,phase
     146              :       complex :: c_phase, zzConjgTemp
     147            0 :       real    :: eig(neigd),cp_time(9),efermi
     148              :       logical :: l_p0,l_bkpts,l_proj,l_amn,l_mmn
     149              : !!! energy window boundaries
     150              : !      integer, allocatable :: nv(:)
     151              : !      integer, allocatable :: nv_b(:)
     152              : !      integer, allocatable :: k1(:,:),k2(:,:),k3(:,:)
     153              : !      integer, allocatable :: k1_b(:,:),k2_b(:,:),k3_b(:,:)
     154              : 
     155            0 :       real,    allocatable :: eigg(:,:)
     156            0 :       real kpoints(nkptd)
     157              : !!! a and b coeff. constructed for each k-point
     158            0 :       complex, allocatable :: acof(:,:,:,:)
     159            0 :       complex, allocatable :: bcof(:,:,:,:)
     160            0 :       complex, allocatable :: ccof(:,:,:,:,:)
     161              : !!! the parameters for the number of wfs
     162              :       integer :: nwfs
     163              : !!! the potential in the spheres and the vacuum
     164            0 :       real, allocatable :: vr(:,:,:),vz(:,:,:)
     165              : !!! bkpts data
     166              :       integer nntot,ikpt_help
     167              :       integer, allocatable :: gb(:,:,:),bpt(:,:)
     168              : !!! radial wavefunctions in the muffin-tins and more ...
     169            0 :       real,    allocatable :: flo(:,:,:,:,:)
     170            0 :       real,    allocatable :: ff(:,:,:,:,:),gg(:,:,:,:,:)
     171              : 
     172            0 :       real     :: uuilon(nlod,ntypd),duilon(nlod,ntypd)
     173            0 :       real     :: ulouilopn(nlod,nlod,ntypd)
     174              : !!! energy parameters
     175              :       real    :: ello(nlod,ntypd,max(2,jspd)),evac(2,max(2,jspd))
     176              :       real    :: epar(0:lmaxd,ntypd,max(2,jspd)),evdu(2,max(jspd,2))
     177              : !!! the Mmn matrices
     178            0 :       complex, allocatable :: surfcurr(:,:,:,:)
     179            0 :       complex, allocatable :: hsomtx(:,:,:,:,:),hsomtx_l(:,:,:,:)
     180            0 :       complex, allocatable :: hsomtxvec(:,:,:,:,:,:)
     181            0 :       complex, allocatable :: mmnk(:,:,:,:),mmn(:,:,:)
     182            0 :       complex, allocatable :: mmn0at(:,:,:,:)
     183            0 :       complex, allocatable :: perpmag(:,:,:)
     184            0 :       complex, allocatable :: amn(:,:,:),nablamat(:,:,:,:)
     185              :       complex, allocatable :: socamn(:,:,:)
     186            0 :       complex, allocatable :: socmmn(:,:,:)
     187            0 :       complex, allocatable :: torque(:,:,:,:,:)
     188              :       complex, allocatable :: socmmnk(:,:,:,:)
     189              :       complex, allocatable :: omn(:,:,:),vec(:,:),a(:)
     190              :       complex, allocatable :: smn(:,:,:)
     191              :       complex, allocatable :: psiw(:,:,:)
     192              : c..wf-hamiltonian in reciprocal space Hwf(k)
     193              :       complex, allocatable :: hwfk(:,:,:)
     194              : c..egenvalues of the wf-hamiltonian Hwf(k)
     195              :       real, allocatable :: eigw(:,:)
     196              : c..wf-hamiltonian in real space (hopping in the same unit cell)
     197              :       complex, allocatable :: hwfr(:,:)
     198              :       real, allocatable :: ei(:)
     199              :       complex, allocatable :: work(:)
     200              :       real,allocatable::centers(:,:,:)
     201              :       integer, allocatable :: iwork(:)
     202              :       integer :: info,lwork,lrwork,liwork,lee
     203              :       character :: jobz, uplo
     204              :       character(len=10)  :: torquename
     205              :       character(len=14)  :: torquenamekpts
     206              :       logical :: l_file
     207              :       logical :: l_amn2, l_conjugate
     208              :       character(len=3) :: spin12(2)
     209              :       data   spin12/'WF1' , 'WF2'/
     210              :       character(len=30)  :: task
     211            0 :       integer,allocatable::irreduc(:)
     212            0 :       integer,allocatable::mapkoper(:)
     213              :       integer :: fullnkpts,kpt,kptibz,kptibz_b,j1,j2,j3,oper,oper_b,k
     214              :       real :: bkrot(3)
     215              :       integer :: ios,kplot,kplotoper,plotoper,gfcut
     216              :       integer :: funbas_start,funbas_stop,jspin5
     217              :       complex :: phasust
     218              :       integer,allocatable::pair_to_do(:,:)
     219            0 :       integer :: ngopr1(natd),nbnd1p2
     220              :       integer,allocatable::maptopair(:,:,:)
     221              :       integer :: ldu,ldv,wannierspin,jspin2
     222              :       complex, allocatable :: amn_copy(:,:)
     223              :       complex, allocatable :: sunit(:,:)
     224              :       real,    allocatable :: sdiag(:)
     225              :       complex, allocatable :: umn(:,:),vmn(:,:)
     226              :       real, allocatable :: rwork(:)
     227              :       character :: jobu*1,jobv*1
     228              :       real,allocatable::kdiff(:,:)
     229            0 :       integer,allocatable :: shiftkpt(:,:)
     230              :       integer :: unigrid(6),gfthick
     231              :       complex,allocatable::ujug(:,:,:,:),ujdg(:,:,:,:)
     232              :       complex,allocatable::djug(:,:,:,:),djdg(:,:,:,:)
     233              :       complex,allocatable::ujulog(:,:,:,:)
     234              :       complex,allocatable::djulog(:,:,:,:)
     235              :       complex,allocatable::ulojug(:,:,:,:)
     236              :       complex,allocatable::ulojdg(:,:,:,:)
     237              :       complex,allocatable::ulojulog(:,:,:,:)
     238              :       real delta,delta1,time_interstitial,time_mmn
     239              :       real time_total,delta2,delta3
     240              :       real time_lapw_expand,time_rw,time_symm,time_film
     241              :       real time_lapw_plot,time_ujugaunt,time_abcof
     242              :       integer :: n_start,n_end,mlotot,mlolotot,err
     243              :       integer :: mlot_d,mlolot_d,ilo,dir,l1,l2,ii,jj
     244              :       character(len=2)    :: spin012(0:2)
     245              :       data spin012/'  ', '.1', '.2'/
     246              :       character(len=6)    :: filename
     247              :       real                :: arg,hescale
     248              :       complex             :: nsfactor,nsfactor_b,value
     249              :       real                :: b1(3),b2(3),b3(3)
     250            0 :       real,allocatable    :: radial1_ff(:,:,:,:,:)
     251            0 :       real,allocatable    :: radial1_gg(:,:,:,:,:)
     252            0 :       real,allocatable    :: radial1_fg(:,:,:,:,:)
     253            0 :       real,allocatable    :: radial1_gf(:,:,:,:,:)
     254            0 :       real,allocatable    :: radial1_flo(:,:,:,:,:)
     255            0 :       real,allocatable    :: radial1_glo(:,:,:,:,:)
     256            0 :       real,allocatable    :: radial1_lof(:,:,:,:,:)
     257            0 :       real,allocatable    :: radial1_log(:,:,:,:,:)
     258            0 :       real,allocatable    :: radial1_lolo(:,:,:,:,:)
     259              : 
     260              :       real,parameter      :: bohrtocm=0.529177e-8
     261              :       real,parameter      :: condquant=7.7480917e-5
     262              :       real                :: dirfacs(3)
     263            0 :       complex,allocatable :: ubeffug(:,:,:)
     264            0 :       complex,allocatable :: ubeffdg(:,:,:)
     265            0 :       complex,allocatable :: dbeffug(:,:,:)
     266            0 :       complex,allocatable :: dbeffdg(:,:,:)
     267            0 :       complex,allocatable :: ubeffulog(:,:,:,:)
     268            0 :       complex,allocatable :: dbeffulog(:,:,:,:)
     269            0 :       complex,allocatable :: ulobeffug(:,:,:,:)
     270            0 :       complex,allocatable :: ulobeffdg(:,:,:,:)
     271            0 :       complex,allocatable :: ulobeffulog(:,:,:,:,:)
     272              :       integer :: spinloop1,spinloop2,ilop
     273              :       integer :: alerr,nbasfcn
     274              : 
     275            0 :       TYPE(t_usdus) :: usdus
     276            0 :       TYPE(t_mat)   :: zMat(2),zzMat(2)
     277            0 :       TYPE(t_lapw)  :: lapw
     278              : 
     279              : !      IF(l_ss) CALL juDFT_error("spin-spiral not yet in this version"
     280              : !     +     ,calledby ="wann_updown")
     281              : 
     282              : c-----initializations
     283              :       
     284            0 :       call timestart("wann_updown")
     285            0 :       time_interstitial=0.0
     286            0 :       time_mmn=0.0
     287            0 :       time_total=0.0
     288            0 :       time_ujugaunt=0.0
     289            0 :       time_abcof=0.0
     290            0 :       time_rw=0.0
     291            0 :       time_symm=0.0
     292            0 :       time_film=0.0
     293            0 :       ngopr1(:)=1
     294              : 
     295            0 :       call cpu_time(time_total)
     296              : 
     297            0 :       l_p0 = .false.
     298            0 :       if (irank.eq.0) l_p0 = .true.
     299              : 
     300            0 :       sfp = 2* sqrt( pimach() )
     301            0 :       tpi = 2* pimach()
     302            0 :       lmd = lmaxd*(lmaxd+2)
     303              : 
     304              : !!!   should be changed in case the windows are really used
     305            0 :       nkpts = nkpt
     306              : 
     307              : c-----read the input file to determine what to do
     308              : !      call wann_read_inp(input,l_p0,wann)  !call wann_read_inp now in fleur_init
     309              : 
     310            0 :       IF(wann%l_byenergy.AND.wann%l_byindex) CALL juDFT_error
     311            0 :      +     ("byenergy.and.byindex",calledby ="wann_updown")
     312            0 :       IF(wann%l_byenergy.AND.wann%l_bynumber) CALL juDFT_error
     313            0 :      +     ("byenergy.and.bynumber",calledby ="wann_updown")
     314            0 :       IF(wann%l_bynumber.AND.wann%l_byindex) CALL juDFT_error
     315            0 :      +     ("bynumber.and.byindex",calledby ="wann_updown")
     316            0 :       IF(.NOT.(wann%l_bynumber.OR.wann%l_byindex.OR.wann%l_byenergy))
     317              :      &     CALL juDFT_error("no rule to sort bands",calledby
     318            0 :      +     ="wann_updown")
     319              : 
     320            0 :       efermi=ef
     321            0 :       if(.not.wann%l_fermi)efermi=0.0
     322              : 
     323              : #ifdef CPP_MPI
     324            0 :          call MPI_BARRIER(mpi_communicatior,ierr)
     325              : #endif
     326              : 
     327              : c**************************************************************
     328              : c   for bzsym=.true.: determine mapping between kpts and w90kpts
     329              : c**************************************************************
     330            0 :       if (wann%l_bzsym) then
     331            0 :          l_file=.false.
     332            0 :          inquire(file='w90kpts',exist=l_file)
     333            0 :          IF(.NOT.l_file)  CALL juDFT_error
     334              :      +        ("w90kpts not found, needed if bzsym",calledby
     335            0 :      +        ="wann_updown")
     336            0 :          open(412,file='w90kpts',form='formatted')
     337            0 :          read(412,*)fullnkpts
     338            0 :          close(412)
     339            0 :          if(l_p0)print*,"fullnkpts=",fullnkpts
     340            0 :          IF(fullnkpts<nkpts) CALL juDFT_error("fullnkpts.lt.nkpts"
     341            0 :      +        ,calledby ="wann_updown")
     342            0 :          allocate(irreduc(fullnkpts),stat=alerr)
     343            0 :          if (alerr /= 0) call juDFT_error('alerr: 1',
     344            0 :      &         calledby='wann_updown')
     345            0 :          allocate(mapkoper(fullnkpts),stat=alerr)
     346            0 :          if (alerr /= 0) call juDFT_error('alerr: 2',
     347            0 :      &         calledby='wann_updown')
     348            0 :          allocate(shiftkpt(3,fullnkpts),stat=alerr)
     349            0 :          if (alerr /= 0) call juDFT_error('alerr: 3',
     350            0 :      &         calledby='wann_updown')
     351            0 :          l_file=.false.
     352            0 :          inquire(file='kptsmap',exist=l_file)
     353            0 :          IF(.NOT.l_file)  CALL juDFT_error
     354              :      +        ("kptsmap not found, needed if bzsym",calledby
     355            0 :      +        ="wann_updown")
     356            0 :          open(713,file='kptsmap')
     357            0 :          do i=1,fullnkpts
     358            0 :             read(713,*)kpt,irreduc(i),mapkoper(i),shiftkpt(:,i)
     359            0 :             IF(kpt/=i) CALL juDFT_error("kpt.ne.i",
     360            0 :      +                                  calledby ="wann_updown")
     361            0 :             if(l_p0)print*,i,irreduc(i),mapkoper(i)
     362              :          enddo
     363            0 :          close(713)
     364            0 :          IF(MAXVAL(irreduc(:))/=nkpts) CALL juDFT_error
     365            0 :      +        ("max(irreduc(:))/=nkpts",calledby ="wann_updown")
     366              :       else
     367            0 :          fullnkpts=nkpts
     368              :       endif
     369              : 
     370            0 :       nrec = 0
     371            0 :       if(l_p0)then
     372            0 :       write (*,*) 'fermi energy:',efermi
     373            0 :       write (*,*) 'emin,emax=',e1s,e2s
     374              :       endif
     375            0 :       nlotot = 0
     376            0 :       do n = 1, ntype
     377            0 :         do l = 1,nlo(n)
     378            0 :           nlotot = nlotot + neq(n) * ( 2*llo(l,n) + 1 )
     379              :         enddo
     380              :       enddo
     381              : 
     382              : c*********************************************************
     383              : cccccccccccccccc   initialize the potential   cccccccccccc
     384              : c*********************************************************
     385              : 
     386            0 :       allocate (vz(nmzd,2,jspd),stat=alerr)
     387            0 :          if (alerr /= 0) call juDFT_error('alerr: 4',
     388            0 :      &         calledby='wann_updown')
     389            0 :       allocate (vr(jmtd,ntypd,jspd),stat=alerr)
     390            0 :          if (alerr /= 0) call juDFT_error('alerr: 5',
     391            0 :      &         calledby='wann_updown')
     392              : 
     393            0 :       vz = REAL(vTot%vac(:,1,:,:))
     394              : 
     395            0 :       do jspin = 1,jspins
     396            0 :         do n = 1, ntype
     397            0 :           do j = 1,jri(n)
     398            0 :             vr(j,n,jspin) = vTot%mt(j,0,n,jspin)
     399              :           enddo
     400              :         enddo
     401              :       enddo
     402              : 
     403              : cccccccccccccccc   end of the potential part  ccccccccccc
     404            0 :       wannierspin=jspd
     405            0 :       if(l_soc) wannierspin=2
     406              : 
     407              : !      allocate ( nv(wannierspin) )
     408              : !      allocate ( k1(nvd,wannierspin),k2(nvd,wannierspin),
     409              : !     &           k3(nvd,wannierspin) )
     410              : 
     411              : 
     412            0 :       allocate ( ff(ntypd,jmtd,2,0:lmaxd,2),stat=alerr )
     413            0 :          if (alerr /= 0) call juDFT_error('alerr: 6',
     414            0 :      &         calledby='wann_updown')
     415            0 :       allocate ( gg(ntypd,jmtd,2,0:lmaxd,2),stat=alerr )
     416            0 :          if (alerr /= 0) call juDFT_error('alerr: 7',
     417            0 :      &         calledby='wann_updown')
     418            0 :       allocate ( usdus%us(0:lmaxd,ntypd,2),stat=alerr )
     419            0 :          if (alerr /= 0) call juDFT_error('alerr: 8',
     420            0 :      &         calledby='wann_updown')
     421            0 :       allocate ( usdus%uds(0:lmaxd,ntypd,2),stat=alerr )
     422            0 :          if (alerr /= 0) call juDFT_error('alerr: 9',
     423            0 :      &         calledby='wann_updown')
     424            0 :       allocate ( usdus%dus(0:lmaxd,ntypd,2),stat=alerr )
     425            0 :          if (alerr /= 0) call juDFT_error('alerr: 10',
     426            0 :      &         calledby='wann_updown')
     427            0 :       allocate ( usdus%duds(0:lmaxd,ntypd,2),stat=alerr )
     428            0 :          if (alerr /= 0) call juDFT_error('alerr: 11',
     429            0 :      &         calledby='wann_updown')
     430            0 :       allocate ( usdus%ddn(0:lmaxd,ntypd,2),stat=alerr )
     431            0 :          if (alerr /= 0) call juDFT_error('alerr: 12',
     432            0 :      &         calledby='wann_updown')
     433            0 :       allocate ( usdus%ulos(nlod,ntypd,2),stat=alerr )
     434            0 :          if (alerr /= 0) call juDFT_error('alerr: 13',
     435            0 :      &         calledby='wann_updown')
     436            0 :       allocate ( usdus%dulos(nlod,ntypd,2),stat=alerr )
     437            0 :          if (alerr /= 0) call juDFT_error('alerr: 14',
     438            0 :      &         calledby='wann_updown')
     439            0 :       allocate ( usdus%uulon(nlod,ntypd,2),stat=alerr )
     440            0 :          if (alerr /= 0) call juDFT_error('alerr: 15',
     441            0 :      &         calledby='wann_updown')
     442            0 :       allocate ( usdus%dulon(nlod,ntypd,2),stat=alerr )
     443            0 :          if (alerr /= 0) call juDFT_error('alerr: 16',
     444            0 :      &         calledby='wann_updown')
     445            0 :       allocate ( usdus%uloulopn(nlod,nlod,ntypd,2),stat=alerr )
     446            0 :          if (alerr /= 0) call juDFT_error('alerr: 17',
     447            0 :      &         calledby='wann_updown')
     448              : c****************************************************
     449              : c     jspin=1
     450              : c****************************************************
     451            0 :       jspin=1
     452              : 
     453              : c...read number of bands and wannier functions from file proj
     454              : 
     455              : c..reading the proj.1 / proj.2 / proj file
     456            0 :        l_proj=.false.
     457            0 :        do j=jspin,0,-1
     458            0 :          inquire(file=trim('proj'//spin012(j)),exist=l_proj)
     459            0 :          if(l_proj)then
     460            0 :             filename='proj'//spin012(j)
     461            0 :             exit
     462              :          endif
     463              :        enddo
     464              : 
     465            0 :        if(l_proj)then
     466            0 :          open (203,file=trim(filename),status='old')
     467            0 :          rewind (203)
     468            0 :          read (203,*) nwfs,numbands
     469            0 :          rewind (203)
     470            0 :          close (203)
     471              :        else
     472              :           CALL juDFT_error("no proj/proj.1/proj.2",
     473            0 :      +                     calledby ="wann_updown")
     474              :        endif
     475              : 
     476              : 
     477              :        jspin2=jspin
     478            0 :        if(l_soc .and. jspins.eq.1)jspin2=1
     479            0 :        jsp_start = jspin ; jsp_end = jspin
     480              : 
     481              : cccccccccccc   read in the eigenvalues and vectors   cccccc
     482              : 
     483              : 
     484            0 :        do jspin5=1,jspins
     485              : !          CALL judft_error("TODO:wann_updown")
     486              : !         call cdn_read0(eig_id,irank,isize,jspin,jspd,l_noco,
     487              :  !    <                  ello,evac,epar,wk,n_bands,n_size)
     488              : 
     489              :         CALL cdn_read0(eig_id,fmpi%irank,fmpi%isize,jspin5,
     490              :      >   input%jspins, !wannierspin instead of DIMENSION%jspd?&
     491            0 :      >  noco%l_noco, n_bands,n_size)
     492              : 
     493              : 
     494              :        enddo
     495              : 
     496              :        jspin=1
     497              : 
     498              : c.. now we want to define the maximum number of the bands by all kpts
     499              : 
     500            0 :       nbnd = 0
     501              : 
     502            0 :       i_rec = 0 ; n_rank = 0
     503              : c*************************************************************
     504              : c..writing down the eig.1 and/or eig.2 files
     505              : c*************************************************************
     506            0 :       if(l_p0)then
     507              :          call wann_write_eig(
     508              :      >           fmpi,cell,noco,nococonv,input,kpts,sym,atoms,           
     509              :      >            eig_id,l_real,
     510              :      >            ntypd,nvd,wannierspin,
     511              :      >            isize,jspin,
     512              :      >            l_ss,l_noco,nrec,fullnkpts,
     513              :      >            wann%l_bzsym,wann%l_byindex,wann%l_bynumber,
     514              :      >            wann%l_byenergy,
     515              :      >            irreduc,wann%band_min(jspin),
     516              :      >            wann%band_max(jspin),
     517              :      >            numbands,
     518              :      >            e1s,e2s,efermi,wann%l_pauli,nkpts,
     519            0 :      <            nbnd,kpoints,.false.,-1)
     520              :        
     521              :       endif!l_p0
     522              : 
     523              : ! nbnd is calculated for process zero and is sent here to the others
     524              : #ifdef CPP_MPI
     525            0 :       if(l_p0)then
     526            0 :          do cpu_index=1,isize-1
     527              :       call MPI_SEND(nbnd,1,MPI_INTEGER,cpu_index,1,mpi_communicatior,
     528            0 :      & ierr)
     529              :          enddo
     530              :       else
     531              :        call MPI_RECV(nbnd,1,MPI_INTEGER,0,1,mpi_communicatior,stt,
     532            0 :      & ierr)
     533              :       endif
     534              : #endif
     535              : 
     536            0 :        print*,"process: ",irank," nbnd= ",nbnd
     537              : c##################################################################
     538              : 
     539            0 :       allocate ( mmn(nbnd,nbnd,fullnkpts),stat=alerr  )
     540            0 :          if (alerr /= 0) call juDFT_error('alerr: 18',
     541            0 :      &         calledby='wann_updown')
     542            0 :       mmn(:,:,:) = cmplx(0.,0.)
     543            0 :       if(wann%l_mmn0at)then
     544            0 :          allocate( mmn0at(nbnd,nbnd,wann%atomlist_num,fullnkpts))
     545            0 :                             mmn0at = cmplx(0.,0.)
     546              :       endif ! wann%l_mmn0at
     547              : 
     548            0 :       if((l_soc.or.l_noco) .and. (jspin.eq.1))
     549            0 :      &   allocate( socmmn(nbnd,nbnd,fullnkpts),stat=alerr  )
     550            0 :          if (alerr /= 0) call juDFT_error('alerr: 19',
     551            0 :      &         calledby='wann_updown')
     552            0 :       if(wann%l_nabla)then
     553            0 :          allocate ( nablamat(3,nbnd,nbnd,fullnkpts),stat=alerr  )
     554            0 :          if (alerr /= 0) call juDFT_error('alerr: 20',
     555            0 :      &         calledby='wann_updown')
     556            0 :          nablamat = cmplx(0.,0.)
     557              :       endif
     558              : 
     559            0 :       if(wann%l_surfcurr)then
     560            0 :          allocate ( surfcurr(3,nbnd,nbnd,fullnkpts),stat=alerr  )
     561            0 :          if (alerr /= 0) call juDFT_error('alerr: 21',
     562            0 :      &         calledby='wann_updown')
     563            0 :          surfcurr = cmplx(0.,0.)
     564              :       endif
     565              : 
     566            0 :       if(wann%l_socmat)then
     567            0 :          allocate ( hsomtx(2,2,nbnd,nbnd,fullnkpts),stat=alerr  )
     568            0 :          if (alerr /= 0) call juDFT_error('alerr: 22',
     569            0 :      &         calledby='wann_updown')
     570            0 :          allocate ( hsomtx_l(nbnd,nbnd,2,2),stat=alerr  )
     571            0 :          if (alerr /= 0) call juDFT_error('alerr: 23',
     572            0 :      &         calledby='wann_updown')
     573              :       endif
     574              : 
     575            0 :       if(wann%l_socmatvec)then
     576            0 :          allocate ( hsomtxvec(2,2,nbnd,nbnd,3,fullnkpts),stat=alerr  )
     577            0 :          if (alerr /= 0) call juDFT_error('alerr: 24',
     578            0 :      &         calledby='wann_updown')
     579            0 :          hsomtxvec=0.0
     580              :       endif
     581              : 
     582            0 :       write (*,*) 'nwfs=',nwfs
     583              : 
     584            0 :       allocate ( flo(ntypd,jmtd,2,nlod,2),stat=alerr )
     585            0 :          if (alerr /= 0) call juDFT_error('alerr: 25',
     586            0 :      &         calledby='wann_updown')
     587              : 
     588              : c***********************************************************
     589              : c     Calculate the radial wave functions for the two spins.
     590              : c***********************************************************
     591            0 :       do jspin=1,wannierspin
     592            0 :         jspin2=jspin
     593            0 :         if(l_soc .and. jspins.eq.1)jspin2=1
     594            0 :         do n = 1,ntype
     595            0 :           do l = 0,lmax(n)
     596              :              call radfun(
     597              :      > l,n,jspin,enpara%el0(l,n,jspin2),vr(1,n,jspin2),atoms,
     598              :      < ff(n,:,:,l,jspin),gg(n,:,:,l,jspin),usdus,
     599            0 :      < nodeu,noded,wronk)
     600              :           enddo !l
     601            0 :           do ilo = 1, nlo(n)
     602              : 
     603              :              call radflo(
     604              :      > atoms,n,jspin,enpara%ello0(:,:,jspin2),vr(1,n,jspin2),
     605              :      >  ff(n,1:,1:,0:,jspin),gg(n,1:,1:,0:,jspin),fmpi,
     606            0 :      <  usdus,uuilon,duilon,ulouilopn,flo(n,:,:,:,jspin))
     607              : 
     608              :           enddo !ilo
     609              :         enddo !n
     610              :       enddo !jspin
     611              : 
     612              : c****************************************************************
     613              : c     Calculate the radial overlap integrals.
     614              : c****************************************************************
     615            0 :       if(wann%l_mmn0.or.wann%l_mmn0at)then
     616            0 :        allocate( radial1_ff(2,2,0:lmaxd,0:lmaxd,ntype),stat=alerr  )
     617            0 :          if (alerr /= 0) call juDFT_error('alerr: 26',
     618            0 :      &         calledby='wann_updown')
     619            0 :        allocate( radial1_fg(2,2,0:lmaxd,0:lmaxd,ntype),stat=alerr  )
     620            0 :          if (alerr /= 0) call juDFT_error('alerr: 27',
     621            0 :      &         calledby='wann_updown')
     622            0 :        allocate( radial1_gf(2,2,0:lmaxd,0:lmaxd,ntype),stat=alerr  )
     623            0 :          if (alerr /= 0) call juDFT_error('alerr: 28',
     624            0 :      &         calledby='wann_updown')
     625            0 :        allocate( radial1_gg(2,2,0:lmaxd,0:lmaxd,ntype),stat=alerr  )
     626            0 :          if (alerr /= 0) call juDFT_error('alerr: 29',
     627            0 :      &         calledby='wann_updown')
     628            0 :        allocate( radial1_flo(2,2,0:lmaxd,1:nlod,ntype),stat=alerr  )
     629            0 :          if (alerr /= 0) call juDFT_error('alerr: 30',
     630            0 :      &         calledby='wann_updown')
     631            0 :        allocate( radial1_glo(2,2,0:lmaxd,1:nlod,ntype),stat=alerr  )
     632            0 :          if (alerr /= 0) call juDFT_error('alerr: 31',
     633            0 :      &         calledby='wann_updown')
     634            0 :        allocate( radial1_lof(2,2,1:nlod,0:lmaxd,ntype),stat=alerr  )
     635            0 :          if (alerr /= 0) call juDFT_error('alerr: 32',
     636            0 :      &         calledby='wann_updown')
     637            0 :        allocate( radial1_log(2,2,1:nlod,0:lmaxd,ntype),stat=alerr  )
     638            0 :          if (alerr /= 0) call juDFT_error('alerr: 33',
     639            0 :      &         calledby='wann_updown')
     640            0 :        allocate( radial1_lolo(2,2,1:nlod,1:nlod,ntype),stat=alerr  )
     641            0 :          if (alerr /= 0) call juDFT_error('alerr: 34',
     642            0 :      &         calledby='wann_updown')
     643              : 
     644            0 :        do n=1,ntype
     645            0 :         do l1=0,lmaxd
     646            0 :          do l2=0,lmaxd
     647            0 :           do spinloop1=1,2
     648            0 :            do spinloop2=1,2
     649              :             call wann_radovlp_integrals(
     650              :      >         rmsh(:,n),jri(n),dx(n),
     651              :      >         ff(n,:,:,l1,spinloop2),
     652              :      >         ff(n,:,:,l2,spinloop1),
     653            0 :      <         radial1_ff(spinloop2,spinloop1,l1,l2,n))
     654              :             call wann_radovlp_integrals(
     655              :      >         rmsh(:,n),jri(n),dx(n),
     656              :      >         ff(n,:,:,l1,spinloop2),
     657              :      >         gg(n,:,:,l2,spinloop1),
     658            0 :      <         radial1_fg(spinloop2,spinloop1,l1,l2,n))
     659              :             call wann_radovlp_integrals(
     660              :      >         rmsh(:,n),jri(n),dx(n),
     661              :      >         gg(n,:,:,l1,spinloop2),
     662              :      >         ff(n,:,:,l2,spinloop1),
     663            0 :      <         radial1_gf(spinloop2,spinloop1,l1,l2,n))
     664              :             call wann_radovlp_integrals(
     665              :      >         rmsh(:,n),jri(n),dx(n),
     666              :      >         gg(n,:,:,l1,spinloop2),
     667              :      >         gg(n,:,:,l2,spinloop1),
     668            0 :      <         radial1_gg(spinloop2,spinloop1,l1,l2,n))
     669              :            enddo !spinloop2
     670              :           enddo !spinloop1
     671              :          enddo !l2
     672              :         enddo !l1
     673              :        enddo !n
     674              : 
     675            0 :        do n=1,ntype
     676            0 :         do ilo=1,nlo(n)
     677            0 :          do ilop=1,nlo(n)
     678            0 :           do spinloop1=1,2
     679            0 :            do spinloop2=1,2
     680              :              call wann_radovlp_integrals(
     681              :      > rmsh(:,n),jri(n),dx(n),
     682              :      > flo(n,:,:,ilo,spinloop2),
     683              :      > flo(n,:,:,ilop,spinloop1),
     684            0 :      < radial1_lolo(spinloop2,spinloop1,ilo,ilop,n))
     685              :            enddo !spinloop2
     686              :           enddo !spinloop1
     687              :          enddo!ilop
     688              :         enddo !ilo
     689              :        enddo !n
     690              : 
     691            0 :        do n=1,ntype
     692            0 :         do ilo=1,nlo(n)
     693            0 :          do l2=0,lmaxd
     694            0 :           do spinloop1=1,2
     695            0 :            do spinloop2=1,2
     696              :              call wann_radovlp_integrals(
     697              :      > rmsh(:,n),jri(n),dx(n),
     698              :      > flo(n,:,:,ilo,spinloop2),
     699              :      > ff(n,:,:,l2,spinloop1),
     700            0 :      < radial1_lof(spinloop2,spinloop1,ilo,l2,n))
     701              :            enddo !spinloop2
     702              :           enddo !spinloop1
     703              :          enddo!l2
     704              :         enddo !ilo
     705              :        enddo !n
     706              : 
     707            0 :        do n=1,ntype
     708            0 :         do ilo=1,nlo(n)
     709            0 :          do l2=0,lmaxd
     710            0 :           do spinloop1=1,2
     711            0 :            do spinloop2=1,2
     712              :              call wann_radovlp_integrals(
     713              :      > rmsh(:,n),jri(n),dx(n),
     714              :      > flo(n,:,:,ilo,spinloop2),
     715              :      > gg(n,:,:,l2,spinloop1),
     716            0 :      < radial1_log(spinloop2,spinloop1,ilo,l2,n))
     717              :            enddo !spinloop2
     718              :           enddo !spinloop1
     719              :          enddo!l2
     720              :         enddo !ilo
     721              :        enddo !n
     722              : 
     723            0 :        do n=1,ntype
     724            0 :         do ilo=1,nlo(n)
     725            0 :          do l2=0,lmaxd
     726            0 :           do spinloop1=1,2
     727            0 :            do spinloop2=1,2
     728              :              call wann_radovlp_integrals(
     729              :      > rmsh(:,n),jri(n),dx(n),
     730              :      > ff(n,:,:,l2,spinloop2),
     731              :      > flo(n,:,:,ilo,spinloop1),
     732            0 :      < radial1_flo(spinloop2,spinloop1,l2,ilo,n))
     733              :            enddo !spinloop2
     734              :           enddo !spinloop1
     735              :          enddo!l2
     736              :         enddo !ilo
     737              :        enddo !n
     738              : 
     739            0 :        do n=1,ntype
     740            0 :         do ilo=1,nlo(n)
     741            0 :          do l2=0,lmaxd
     742            0 :           do spinloop1=1,2
     743            0 :            do spinloop2=1,2
     744              :              call wann_radovlp_integrals(
     745              :      > rmsh(:,n),jri(n),dx(n),
     746              :      > gg(n,:,:,l2,spinloop2),
     747              :      > flo(n,:,:,ilo,spinloop1),
     748            0 :      < radial1_glo(spinloop2,spinloop1,l2,ilo,n))
     749              :            enddo !spinloop2
     750              :           enddo !spinloop1
     751              :          enddo!l2
     752              :         enddo !ilo
     753              :        enddo !n
     754              : 
     755              :       endif
     756              : 
     757              : 
     758            0 :       jspin=1
     759              : 
     760              : c****************************************************************
     761              : c     compute matrix elements of exchange field
     762              : c****************************************************************
     763              : 
     764            0 :       if(wann%l_perpmag)then
     765            0 :          allocate( perpmag(nbnd,nbnd,fullnkpts),stat=alerr  )
     766            0 :          if (alerr /= 0) call juDFT_error('alerr: 35',
     767            0 :      &         calledby='wann_updown')
     768            0 :          perpmag=0.0
     769            0 :          allocate( ubeffug(0:lmd,0:lmd,1:ntype),stat=alerr  )
     770            0 :          if (alerr /= 0) call juDFT_error('alerr: 36',
     771            0 :      &         calledby='wann_updown')
     772            0 :          allocate( ubeffdg(0:lmd,0:lmd,1:ntype),stat=alerr  )
     773            0 :          if (alerr /= 0) call juDFT_error('alerr: 37',
     774            0 :      &         calledby='wann_updown')
     775            0 :          allocate( dbeffug(0:lmd,0:lmd,1:ntype),stat=alerr  )
     776            0 :          if (alerr /= 0) call juDFT_error('alerr: 38',
     777            0 :      &         calledby='wann_updown')
     778            0 :          allocate( dbeffdg(0:lmd,0:lmd,1:ntype),stat=alerr  )
     779            0 :          if (alerr /= 0) call juDFT_error('alerr: 39',
     780            0 :      &         calledby='wann_updown')
     781              :          allocate( ubeffulog(0:lmd,1:nlod,-llod:llod,1:ntype),
     782            0 :      & stat=alerr  )
     783            0 :          if (alerr /= 0) call juDFT_error('alerr: 40',
     784            0 :      &         calledby='wann_updown')
     785              :          allocate( dbeffulog(0:lmd,1:nlod,-llod:llod,1:ntype),
     786            0 :      &         stat=alerr  )
     787            0 :          if (alerr /= 0) call juDFT_error('alerr: 41',
     788            0 :      &         calledby='wann_updown')
     789              :          allocate( ulobeffug(0:lmd,1:nlod,-llod:llod,1:ntype),
     790            0 :      &          stat=alerr  )
     791            0 :          if (alerr /= 0) call juDFT_error('alerr: 42',
     792            0 :      &         calledby='wann_updown')
     793              :          allocate( ulobeffdg(0:lmd,1:nlod,-llod:llod,1:ntype),
     794            0 :      &         stat=alerr  )
     795              :          allocate( ulobeffulog(1:nlod,-llod:llod,
     796            0 :      &                     1:nlod,-llod:llod,1:ntype),stat=alerr  )
     797            0 :          if (alerr /= 0) call juDFT_error('alerr: 43',
     798            0 :      &         calledby='wann_updown')
     799              : 
     800              : 
     801              :          call wann_gabeffgagaunt(
     802              :      >            vTot,
     803              :      >            wann%l_perpmagatlres,wann%perpmagl,
     804              :      >            neq,mlh,nlhd,nlh,ntypsy,llh,lmax,
     805              :      >            nmem,ntype,ntypd,bbmat,bmat,
     806              :      >            nlod,llod,nlo,llo,flo,ff,gg,jri,rmsh,dx,jmtd,
     807              :      >            lmaxd,lmd,clnu,
     808              :      <            ubeffug,ubeffdg,dbeffug,dbeffdg,ubeffulog,
     809              :      <            dbeffulog,
     810            0 :      <            ulobeffug,ulobeffdg,ulobeffulog)
     811              :       endif
     812              : 
     813              : 
     814            0 :       if(wann%l_torque)then
     815            0 :          if(l_soc)then
     816              :             allocate( torque(3,nbnd,nbnd,
     817            0 :      &             wann%atomlist_num,fullnkpts) )
     818              :          else
     819            0 :             nbnd1p2=nbnd+nbnd
     820              :             allocate( torque(3,nbnd1p2,nbnd1p2,
     821            0 :      &             wann%atomlist_num,fullnkpts) )
     822              :          endif
     823            0 :          torque=cmplx(0.0,0.0)
     824              :       endif
     825              : 
     826              : 
     827            0 :       i_rec = 0 ; n_rank = 0
     828              : 
     829              : c****************************************************************
     830              : c..   loop by kpoints starts!      each may be a separate task
     831              : c****************************************************************
     832            0 :       do 10 ikpt = wann%ikptstart,fullnkpts  ! loop by k-points starts
     833            0 :         kptibz=ikpt
     834            0 :         if(wann%l_bzsym) kptibz=irreduc(ikpt)
     835            0 :         if(wann%l_bzsym) oper=mapkoper(ikpt)
     836              : 
     837            0 :         i_rec = i_rec + 1
     838            0 :       if (mod(i_rec-1,isize).eq.irank) then
     839              : 
     840              : c****************************************************************
     841              : c     Obtain eigenvectors for the two spin channels.
     842              : c****************************************************************
     843              : 
     844            0 :       allocate ( eigg(neigd,2),stat=alerr  )
     845            0 :          if (alerr /= 0) call juDFT_error('alerr: 44',
     846            0 :      &         calledby='wann_updown')
     847              : 
     848              : 
     849            0 :       call cpu_time(delta)
     850            0 :       do jspin=1,wannierspin
     851              :          if(jspin.eq.1)nrec=0
     852              :          if(jspin.eq.2)nrec=nkpts
     853              :          if(l_noco)nrec=0
     854            0 :          n_start=1
     855            0 :          n_end=neigd
     856              : 
     857              : 
     858              : 
     859              :         CALL lapw%init(input,noco,nococonv,kpts,atoms,sym,
     860            0 :      &  kptibz,cell,fmpi)
     861              : 
     862              : 
     863              :         nbasfcn = MERGE(lapw%nv(1)+lapw%nv(2)+2*atoms%nlotot,
     864            0 :      &  lapw%nv(1)+atoms%nlotot,noco%l_noco)
     865              :         
     866            0 :         CALL zzMat(jspin)%init(l_real,nbasfcn,input%neig)
     867            0 :         CALL zMat(jspin)%init(l_real,nbasfcn,input%neig)
     868              : 
     869              : 
     870              :         CALL cdn_read(
     871              :      & eig_id,
     872              :      & lapw%dim_nvd(),input%jspins,fmpi%irank,fmpi%isize, !wannierspin instead of DIMENSION%jspd?
     873              :      & kptibz,jspin,lapw%dim_nbasfcn(),
     874              :      & noco%l_ss,noco%l_noco,input%neig,n_start,n_end,
     875            0 :      & nbands,eigg(:,jspin),zzMat(jspin))
     876              : 
     877              : 
     878              : !         CALL cdn_read(
     879              : !     >              eig_id,
     880              : !     >              nvd,jspd,irank,isize,kptibz,jspin,nbasfcn,
     881              : !     >              l_ss,l_noco,neigd,n_start,n_end,
     882              : !     <              nmat,nv,ello,evdu,epar,
     883              : !     <              k1,k2,k3,bkpt,wk,nbands,eigg(:,jspin),zzMat(jspin))
     884              : 
     885              : 
     886              : 
     887              : 
     888              : 
     889              : !      call judft_error("cdn_read in wann_updown")
     890            0 :          nrec=0
     891              :       enddo !jspin
     892            0 :       call cpu_time(delta1)
     893            0 :       time_rw=time_rw+delta1-delta
     894            0 :       nslibd = 0
     895              : 
     896            0 :       jspin=1
     897              : 
     898              : c...we work only within the energy window
     899              : 
     900              :  
     901              : 
     902            0 :       eig(:) = 0.
     903              : 
     904            0 :       do jspin=1,wannierspin
     905            0 :        jspin2=min(input%jspins,jspin)
     906            0 :        nslibd=0
     907              :        !print*,"bands used:"
     908            0 :        do i = 1,nbands
     909              :         if ((eigg(i,1).ge.e1s.and.nslibd.lt.numbands.
     910              :      &       and.wann%l_bynumber)
     911              :      &    .or.(eigg(i,1).ge.e1s.and.eigg(i,1).le.e2s.
     912              :      &       and.wann%l_byenergy)
     913            0 :      &    .or.(i.ge.wann%band_min(jspin2).and.
     914              :      &        (i.le.wann%band_max(jspin2))
     915            0 :      &       .and.wann%l_byindex))then
     916              : 
     917              :            !print*,i
     918            0 :            nslibd = nslibd + 1
     919            0 :            eig(nslibd) = eigg(i,1)
     920            0 :            if(l_noco)then
     921            0 :              funbas=lapw%nv(1)+atoms%nlotot
     922            0 :              funbas=funbas+lapw%nv(2)+atoms%nlotot
     923              :            else
     924            0 :              funbas=lapw%nv(jspin2)+atoms%nlotot
     925              :            endif
     926            0 :            IF (zzMat(jspin)%l_real) THEN
     927            0 :               do j = 1,funbas
     928            0 :                  zMat(jspin)%data_r(j,nslibd) = zzMat(jspin)%data_r(j,i)
     929              :               enddo
     930              :            ELSE
     931            0 :               do j = 1,funbas
     932            0 :                  zMat(jspin)%data_c(j,nslibd) = zzMat(jspin)%data_c(j,i)
     933              :               enddo
     934              :            END IF
     935              :         endif
     936              :        enddo
     937              :       enddo
     938            0 :       jspin=1
     939              : 
     940              : c***********************************************************
     941              : c              rotate the wavefunction
     942              : c***********************************************************
     943            0 :       do jspin=1,wannierspin
     944            0 :         if (wann%l_bzsym.and.oper.ne.1) then  !rotate bkpt
     945              :            call wann_kptsrotate(
     946              :      >            atoms,
     947              :      >            invsat,
     948              :      >            l_noco,l_soc,
     949              :      >            jspin,
     950              :      >            oper,nop,mrot,nvd,
     951              :      >            shiftkpt(:,ikpt),
     952              :      >          tau,
     953              :      >          lapw,    
     954              : !     x            lapw%bkpt,lapw%k1(:,:),lapw%k2(:,:),lapw%k3(:,:),
     955            0 :      x            zMat(jspin),nsfactor)
     956              :         else
     957            0 :           nsfactor=cmplx(1.0,0.0)
     958              :         endif
     959              :       enddo !jspin
     960              : c      print*,"bkpt1=",bkpt
     961              : 
     962              : c******************************************************************
     963              : 
     964              : c...the overlap matrix Mmn which is computed for each k- and b-point
     965              : 
     966            0 :       noccbd = nslibd
     967              : 
     968              :       allocate ( acof(noccbd,0:lmd,natd,wannierspin),
     969              :      & bcof(noccbd,0:lmd,natd,wannierspin),
     970            0 :      & ccof(-llod:llod,noccbd,nlod,natd,wannierspin),stat=alerr )
     971            0 :          if (alerr /= 0) call juDFT_error('alerr: 49',
     972            0 :      &         calledby='wann_updown')
     973              : 
     974              : 
     975              : 
     976            0 :       acof(:,:,:,:) = cmplx(0.,0.) ; bcof(:,:,:,:) = cmplx(0.,0.)
     977            0 :       ccof(:,:,:,:,:) = cmplx(0.,0.)
     978              : 
     979              : c...generation the A,B,C coefficients in the spheres
     980              : c...for the lapws and local orbitals, summed by the basis functions
     981              : 
     982              : !      ALLOCATE(lapw%k1(SIZE(k1,1),SIZE(k1,2)))
     983              : !      ALLOCATE(lapw%k2(SIZE(k1,1),SIZE(k1,2)))
     984              : !      ALLOCATE(lapw%k3(SIZE(k1,1),SIZE(k1,2)))
     985              : !      lapw%nv = nv
     986              : !      lapw%nmat = nmat
     987              : !      lapw%k1 = k1
     988              : !      lapw%k2 = k2
     989              : !      lapw%k3 = k3
     990              :       ! I think the other variables of lapw are not needed here.
     991              : 
     992            0 :       call cpu_time(delta)
     993            0 :       do jspin=1,wannierspin
     994              : 
     995              :          CALL abcof(input,atoms,sym,cell,lapw,noccbd,usdus,
     996              :      +      noco,nococonv,min(jspin,input%jspins) ,
     997              :      +      acof(1:,0:,1:,jspin),
     998              :      +              bcof(1:,0:,1:,jspin),ccof(-llod:,1:,1:,1:,jspin),
     999            0 :      +              zMat(jspin))
    1000              : 
    1001              :          call wann_abinv(atoms,sym,
    1002              :      X        acof(:,0:,:,jspin),bcof(:,0:,:,jspin),
    1003            0 :      <       ccof(-llod:,:,:,:,jspin))
    1004              :       enddo !jspin
    1005              : 
    1006              : !      DEALLOCATE(lapw%k1,lapw%k2,lapw%k3)
    1007              : 
    1008            0 :       jspin=1
    1009            0 :       call cpu_time(delta1)
    1010            0 :       time_abcof=time_abcof+delta1-delta
    1011              : 
    1012            0 :       if(wann%l_socmat)then
    1013              :          call wann_socmat(
    1014              :      >       fmpi,enpara,input,noco,nococonv,atoms,
    1015              :      >       lmaxd,natd,noccbd,
    1016              :      >       llod,jspd,
    1017              :      >       theta,phi,jspins,irank,
    1018              :      >       vTot%mt,
    1019              :      >       acof,bcof,ccof,
    1020            0 :      <       hsomtx_l)
    1021            0 :         do i1 = 1,2
    1022            0 :          do i2 = 1,2
    1023            0 :           do i = 1,nbnd
    1024            0 :            do j = 1,nbnd
    1025            0 :              hsomtx(i1,i2,i,j,ikpt) = hsomtx_l(i,j,i1,i2)
    1026              :            enddo
    1027              :           enddo
    1028              :          enddo
    1029              :         enddo
    1030              :       endif
    1031              : 
    1032            0 :       if(wann%l_socmatvec)then
    1033              :         call wann_socmat_vec(
    1034              :      >       fmpi,enpara,input,noco,nococonv,atoms,
    1035              :      >       lmaxd,natd,noccbd,
    1036              :      >       llod,jspd,
    1037              :      >       jspins,irank,
    1038              :      >       vTot%mt,! vr,
    1039              :      >       acof,bcof,ccof,
    1040            0 :      <       hsomtxvec(:,:,:,:,:,ikpt))
    1041              :       endif
    1042              : 
    1043            0 :       if(wann%l_perpmag)then
    1044              :          call wann_perpmag(
    1045              :      >            nslibd,neq,mlh,nlhd,nlh,ntypsy,llh,lmax,
    1046              :      >            nmem,ntype,ntypd,bbmat,bmat,
    1047              :      >            nlod,llod,nlo,llo,flo,ff,gg,jri,rmsh,dx,jmtd,
    1048              :      >            lmaxd,lmd,clnu,
    1049              :      >            ubeffug,ubeffdg,dbeffug,dbeffdg,ubeffulog,
    1050              :      >            dbeffulog,
    1051              :      >            ulobeffug,ulobeffdg,ulobeffulog,
    1052              :      >            acof,bcof,ccof,
    1053            0 :      <            perpmag(:,:,ikpt))
    1054              :       endif
    1055              :       
    1056            0 :       if(wann%l_torque)then
    1057              :          call wann_torque(
    1058              :      >       amat,
    1059              :      >       ntype,lmaxd,lmax,natd,
    1060              :      >       neq,noccbd,lmd,natd,llod,nlod,
    1061              :      >       nlo,llo, 
    1062              :      >       acof,bcof,ccof,
    1063              :      >       usdus%us,usdus%dus,usdus%duds,usdus%uds,
    1064              :      >       usdus%ulos,usdus%dulos,
    1065              :      >       rmt,pos,wann%atomlist_num,wann%atomlist, 
    1066              :      >       (l_soc.or.l_noco),nbnd,
    1067            0 :      &       torque(:,:,:,:,ikpt))
    1068              :       endif
    1069              :       
    1070              : 
    1071              : #ifdef CPP_TOPO
    1072              : 
    1073              :       if(wann%l_surfcurr)then
    1074              :          call wann_surfcurr_int_updown(
    1075              :      >       nv2d,jspin,n3d,nmzxyd,n2d,ntypsd,
    1076              :      >       ntype,lmaxd,jmtd,ntypd,natd,nmzd,neq,nq3,nvac,
    1077              :      >       nmz,nmzxy,nq2,nop,nop2,volint,film,slice,symor,
    1078              :      >       invs,invs2,z1,delz,ngopr,ntypsy,jri,
    1079              :      >       pos,taual,zatom,rmt,
    1080              :      >       lmax,mrot,tau,rmsh,invtab,amat,bmat,bbmat,ikpt,
    1081              :      >       nvd,nv(jspin),lapw%bkpt,omtil,
    1082              :      >       lapw%k1(:,jspin),lapw%k2(:,jspin),lapw%k3(:,jspin),
    1083              :      >       nslibd,nbasfcn,neigd,zMat,
    1084              :      >       dirfacs,
    1085              :      >       surfcurr(:,:,:,ikpt))
    1086              :          call wann_surfcurr_updown(
    1087              :      >       dirfacs,amat,
    1088              :      >       ntype,lmaxd,lmax,natd,
    1089              :      >       neq,noccbd,lmd,natd,llod,nlod,
    1090              :      >       nlo,llo,
    1091              :      >       acof,bcof,ccof,
    1092              :      >       us,dus,duds,uds,
    1093              :      >       ulos,dulos,
    1094              :      >       rmt,pos,
    1095              :      &       surfcurr(:,:,:,ikpt))
    1096              :         write(oUnit,*)"dirfacs=",dirfacs
    1097              :       endif
    1098              : 
    1099              :       if(wann%l_nabla)then
    1100              :          call wann_nabla_updown(
    1101              :      >       jspin,ntype,lmaxd,lmax,natd,
    1102              :      >       jmtd,jri,rmsh,dx,neq,
    1103              :      >       noccbd,lmd,natd,llod,nlod,
    1104              :      >       ff(:,:,:,:,1),gg(:,:,:,:,1),
    1105              :      >       ff(:,:,:,:,2),gg(:,:,:,:,2),
    1106              :      >       acof(:,:,:,1),bcof(:,:,:,1),
    1107              :      >       ccof(:,:,:,:,1),
    1108              :      >       acof(:,:,:,2),bcof(:,:,:,2),
    1109              :      >       ccof(:,:,:,:,2),
    1110              :      &       nablamat(:,:,:,ikpt))
    1111              : 
    1112              :          addnoco=0
    1113              :          do 41 i = n_rank+1,nv(jspin),n_size
    1114              :            b1(1)=lapw%bkpt(1)+lapw%k1(i,jspin)
    1115              :            b1(2)=lapw%bkpt(2)+lapw%k2(i,jspin)
    1116              :            b1(3)=lapw%bkpt(3)+lapw%k3(i,jspin)
    1117              :            b2(1)=b1(1)*bmat(1,1)+b1(2)*bmat(2,1)+b1(3)*bmat(3,1)
    1118              :            b2(2)=b1(1)*bmat(1,2)+b1(2)*bmat(2,2)+b1(3)*bmat(3,2)
    1119              :            b2(3)=b1(1)*bmat(1,3)+b1(2)*bmat(2,3)+b1(3)*bmat(3,3)
    1120              :            do 42 j = n_rank+1,nv(jspin),n_size
    1121              : c            b1(1)=bkpt(1)+k1(j,jspin)
    1122              : c            b1(2)=bkpt(2)+k2(j,jspin)
    1123              : c            b1(3)=bkpt(3)+k3(j,jspin)
    1124              : c            b3(1)=b1(1)*bmat(1,1)+b1(2)*bmat(2,1)+b1(3)*bmat(3,1)
    1125              : c            b3(2)=b1(1)*bmat(1,2)+b1(2)*bmat(2,2)+b1(3)*bmat(3,2)
    1126              : c            b3(3)=b1(1)*bmat(1,3)+b1(2)*bmat(2,3)+b1(3)*bmat(3,3)
    1127              : c-->     determine index and phase factor
    1128              :             i1 = lapw%k1(j,jspin) - lapw%k1(i,jspin)
    1129              :             i2 = lapw%k2(j,jspin) - lapw%k2(i,jspin)
    1130              :             i3 = lapw%k3(j,jspin) - lapw%k3(i,jspin)
    1131              :             in = ig(i1,i2,i3)
    1132              :             if (in.eq.0) goto 42
    1133              :             phase   = rgphs(i1,i2,i3)
    1134              :             phasust = cmplx(phase,0.0)*ustep(in)
    1135              : 
    1136              :             do m = 1,nslibd
    1137              :              do n = 1,nslibd
    1138              :               do dir=1,3
    1139              :                IF (zMat%l_real) THEN
    1140              :                   zzConjgTemp = cmplx(zMat(1)%data_r(i+addnoco,m) *
    1141              :      +                                zMat(2)%data_r(j+addnoco,n), 0.0)
    1142              :                ELSE
    1143              :                   zzConjgTemp =       zMat(1)%data_c(i+addnoco,m) *
    1144              :      +                          CONJG(zMat(2)%data_c(j+addnoco,n))
    1145              :                END IF
    1146              :                value=phasust*zzConjgTemp
    1147              :                nablamat(dir,m,n,ikpt) = nablamat(dir,m,n,ikpt) -
    1148              :      +                                  value*b2(dir)
    1149              :               enddo
    1150              :              enddo
    1151              :             enddo
    1152              : 
    1153              :  42        continue
    1154              :  41      continue
    1155              :       endif
    1156              : #endif
    1157              : 
    1158              : c***************************************************
    1159              : c     Calculate matrix elements of the Pauli matrix.
    1160              : c***************************************************
    1161            0 :       if(wann%l_mmn0at)then
    1162              :         call wann_mmk0_updown_sph_at(
    1163              :      >      l_noco,alph,beta,
    1164              :      >      llod,noccbd,nlod,natd,ntypd,lmaxd,lmax,lmd,
    1165              :      >      ntype,neq,nlo,llo,
    1166              :      >      radial1_ff,radial1_gg,radial1_fg,radial1_gf,
    1167              :      >      acof(1:noccbd,:,:,:),
    1168              :      >      bcof(1:noccbd,:,:,:),ccof(:,1:noccbd,:,:,:),
    1169              :      >      usdus%ddn,usdus%uulon,usdus%dulon,usdus%uloulopn,
    1170              :      >      wann%atomlist_num,wann%atomlist(:),
    1171            0 :      =      mmn0at(:,:,:,ikpt))
    1172              :       endif
    1173              : 
    1174              : 
    1175              : 
    1176              : 
    1177            0 :       if(wann%l_mmn0)then
    1178              : 
    1179              : c------> interstitial contribution to updown.mmn0-matrix
    1180            0 :        addnoco=0
    1181            0 :        if(l_noco)then
    1182            0 :           addnoco=lapw%nv(1)+nlotot
    1183              :        endif
    1184              : 
    1185              : c$$$       do i = n_rank+1,nv(jspin),n_size
    1186              : c$$$         do 22 j = n_rank+1,nv(jspin),n_size
    1187              : c$$$
    1188              : c$$$c-->     determine index and phase factor
    1189              : c$$$            i1 = k1(j,jspin) - k1(i,jspin)
    1190              : c$$$            i2 = k2(j,jspin) - k2(i,jspin)
    1191              : c$$$            i3 = k3(j,jspin) - k3(i,jspin)
    1192              : c$$$            in = ig(i1,i2,i3)
    1193              : c$$$            if (in.eq.0) goto 22
    1194              : c$$$            phase   = rgphs(i1,i2,i3)
    1195              : c$$$            phasust = cmplx(phase,0.0)*ustep(in)
    1196              : c$$$            do m = 1,nslibd
    1197              : c$$$             do n = 1,nslibd
    1198              : c$$$#if ( !defined(CPP_INVERSION) || defined(CPP_SOC) )
    1199              : c$$$              mmn(m,n,ikpt) =
    1200              : c$$$     =        mmn(m,n,ikpt) +
    1201              : c$$$     +          phasust*z(i,m,1)*conjg(z(j+addnoco,n,2))
    1202              : c$$$#else
    1203              : c$$$              mmn(m,n,ikpt) =
    1204              : c$$$     =        mmn(m,n,ikpt) +
    1205              : c$$$     +          phasust*cmplx(z(i,m,1)*z(j+addnoco,n,2),0.0)
    1206              : c$$$#endif
    1207              : c$$$             enddo !n
    1208              : c$$$            enddo !m
    1209              : c$$$
    1210              : c$$$  22     continue
    1211              : c$$$       enddo !i
    1212              : 
    1213            0 :       addnoco2 = addnoco      ! TODO: correct for addnoco2
    1214              :       call wann_mmkb_int(
    1215              :      >         cmplx(1.,0.),addnoco,addnoco2,
    1216              :      >         nvd,stars%mx1,stars%mx2,stars%mx3,
    1217              :      > stars%ng3,lapw%k1(:,jspin),lapw%k2(:,jspin),lapw%k3(:,jspin),
    1218              :      > lapw%nv(jspin),neigd,nbasfcn,nbasfcn,zMat(1),nslibd,
    1219              :      > lapw%k1(:,jspin),lapw%k2(:,jspin),lapw%k3(:,jspin),
    1220              :      >         lapw%nv(jspin),zMat(2),nslibd,
    1221              :      >         nbnd,
    1222              :      >         rgphs,ustep,ig,(/ 0,0,0 /),
    1223            0 :      <         mmn(:,:,ikpt))
    1224              : 
    1225              : 
    1226              : c---> spherical contribution to updown.mmn0-matrix
    1227              :        call wann_mmk0_updown_sph(
    1228              :      >          l_noco,alph,beta,
    1229              :      >          llod,noccbd,nlod,natd,ntypd,lmaxd,lmax,lmd,
    1230              :      >          ntype,neq,nlo,llo,
    1231              :      >          radial1_ff,radial1_gg,radial1_fg,radial1_gf,
    1232              :      >          radial1_flo,radial1_glo,
    1233              :      >          radial1_lof,radial1_log,
    1234              :      >          radial1_lolo,
    1235              :      >          acof(1:noccbd,:,:,:),
    1236              :      >          bcof(1:noccbd,:,:,:),ccof(:,1:noccbd,:,:,:),
    1237              :      >          usdus%ddn,usdus%uulon,usdus%dulon,usdus%uloulopn,
    1238            0 :      =          mmn(1:noccbd,1:noccbd,ikpt))
    1239              : 
    1240              : c---> vacuum contribution to updown.mmn0-matrix
    1241              : 
    1242              :        if (film ) then
    1243              : 
    1244              : c         call wann_mmk0_vac(
    1245              : c     >           z1,nmzd,nv2d,k1d,k2d,k3d,n3d,nvac,
    1246              : c     >           ig,nmz,delz,ig2,area,bmat,
    1247              : c     >           bbmat,evac(:,jspin),bkpt,vz(:,:,jspin2),
    1248              : c     >           nslibd,jspin,lapw%k1(:,jspin),lapw%k2(:,jspin),
    1249              : c     >           lapw%k3(:,jspin),wannierspin,nvd,
    1250              : c     >           nbasfcn,neigd,z,nv,omtil,
    1251              : c     <           mmn(1:noccbd,1:noccbd,ikpt))
    1252              : 
    1253              :  
    1254              :        endif
    1255              :       endif !l_mmn0
    1256              : 
    1257            0 :       deallocate ( acof,bcof,ccof,eigg )
    1258              : 
    1259            0 :       write (*,*) 'nslibd=',nslibd
    1260              : 
    1261              : 
    1262              :       endif   ! loop by processors
    1263              : 
    1264            0 : 10    continue ! end of cycle by the k-points
    1265              : 
    1266              : #ifdef CPP_MPI
    1267            0 :       call MPI_BARRIER(mpi_communicatior,ierr)
    1268              : #endif
    1269              : c******************************************************
    1270              : c     Write down the projections.
    1271              : c******************************************************
    1272              :   5   continue
    1273              : 
    1274            0 :       if(wann%l_nabla)then
    1275            0 :         hescale=tpi*condquant/bohrtocm/omtil
    1276            0 :         hescale=sqrt(hescale)
    1277            0 :         nablamat=nablamat*hescale
    1278              :         call wann_write_nabla(
    1279              :      >      mpi_communicatior,l_p0,'nabl.updown',
    1280              :      >      'Matrix elements of nabla operator',
    1281              :      >      nbnd,fullnkpts,nbnd,
    1282              :      >      irank,isize,wann%l_unformatted,
    1283            0 :      <      nablamat)
    1284              :       endif
    1285              : 
    1286            0 :       if(wann%l_socmat)then
    1287              : !        call wann_write_hsomtx(
    1288              : !     >      mpi_communicatior,l_p0,'WF1.hsomtx',
    1289              : !     >      'Matrix elements of spin-orbit coupling',
    1290              : !     >      2,2,nbnd,nbnd,fullnkpts,
    1291              : !     >      irank,isize,
    1292              : !     <      hsomtx)
    1293              : 
    1294              :         call wann_write_matrix5(
    1295              :      >      mpi_communicatior,l_p0,'WF1.hsomtx',
    1296              :      >      'Matrix elements of spin-orbit coupling',
    1297              :      >      nbnd,nbnd,2,2,fullnkpts,
    1298              :      >      irank,isize,wann%socmatfmt==2,
    1299            0 :      <      hsomtx)
    1300              : 
    1301              : 
    1302              :       endif
    1303              : 
    1304            0 :       if(wann%l_socmatvec)then
    1305              :         call wann_write_matrix6(
    1306              :      >      mpi_communicatior,l_p0,'WF1.hsomtxvec',
    1307              :      >      'Vectors for spin-orbit coupling',
    1308              :      >      2,2,nbnd,nbnd,3,fullnkpts,
    1309              :      >      irank,isize,wann%socmatvecfmt==2,
    1310            0 :      <      hsomtxvec)
    1311              :       endif
    1312              : 
    1313            0 :       if(wann%l_mmn0)then
    1314              :         call wann_write_amn(
    1315              :      >      mpi_communicatior,l_p0,'updown.mmn0',
    1316              :      >      'Overlaps of the wavefunct. at the same kpoint',
    1317              :      >      nbnd,fullnkpts,nbnd,
    1318              :      >      irank,isize,.false.,.true.,
    1319            0 :      <      mmn,wann%mmn0fmt==2)
    1320              :       endif !l_soc and l_mmn0
    1321              : 
    1322            0 :       if(wann%l_perpmag)then
    1323              :         call wann_write_amn(
    1324              :      >      mpi_communicatior,l_p0,'updown.perpmag',
    1325              :      >      'Overlaps of the wavefunct. with Beff',
    1326              :      >      nbnd,fullnkpts,nbnd,
    1327              :      >      irank,isize,.true.,.true.,
    1328            0 :      <      perpmag,wann%perpmagfmt==2)
    1329              :       endif !l_soc and l_mmn0
    1330              : 
    1331            0 :       if(wann%l_mmn0at)then
    1332            0 :         do at=1,wann%atomlist_num
    1333            0 :          write(torquename,fmt='("mmn0up_",i3.3)')wann%atomlist(at)
    1334              :          call wann_write_amn(
    1335              :      >      mpi_communicatior,l_p0,torquename,
    1336              :      >      'Atom resolved overlap',
    1337              :      >      nbnd,fullnkpts,nbnd,
    1338              :      >      irank,isize,.true.,wann%l_unformatted,
    1339            0 :      <      mmn0at(:,:,at,:),wann%l_unformatted)
    1340              :         enddo
    1341              :       endif
    1342              : 
    1343              : 
    1344            0 :       if(wann%l_surfcurr)then
    1345            0 :         surfcurr=conjg(surfcurr/cmplx(0.0,2.0))
    1346            0 :         hescale=tpi*condquant/bohrtocm/omtil
    1347            0 :         hescale=sqrt(hescale)
    1348            0 :         surfcurr=surfcurr*hescale
    1349              :         call wann_write_nabla(
    1350              :      >      mpi_communicatior,l_p0,'surfcurr.updown',
    1351              :      >      'Surface currents',
    1352              :      >      nbnd,fullnkpts,nbnd,
    1353              :      >      irank,isize,wann%l_unformatted,
    1354            0 :      <      surfcurr)
    1355              :       endif
    1356              : 
    1357            0 :       if(wann%l_torque)then
    1358            0 :        if(l_soc)then
    1359            0 :           nbnd1p2=nbnd
    1360              :        else
    1361            0 :           nbnd1p2=nbnd+nbnd
    1362              :        endif
    1363            0 :        do at=1,wann%atomlist_num
    1364            0 :          write(torquename,fmt='("torque_",i3.3)')wann%atomlist(at)
    1365              :          call wann_write_nabla(
    1366              :      >      mpi_communicatior,l_p0,torquename,
    1367              :      >      'Spin currents into MT-sphere',
    1368              :      >      nbnd1p2,fullnkpts,nbnd1p2,
    1369              :      >      irank,isize,wann%l_unformatted,
    1370            0 :      <      torque(:,:,:,at,:)) 
    1371              :        enddo
    1372              :       endif
    1373              : 
    1374            0 :       if( allocated (nablamat) ) deallocate( nablamat )
    1375            0 :       deallocate ( mmn )
    1376              : 
    1377            0 :       deallocate(flo)
    1378              : 
    1379            0 :       deallocate ( vr,vz )
    1380            0 :       deallocate ( ff,gg )
    1381            0 :       if (wann%l_bzsym)  deallocate(irreduc,mapkoper)
    1382              : 
    1383            0 :       call cpu_time(delta)
    1384            0 :       time_total=delta-time_total
    1385            0 :       write(oUnit,*)"time_total=",time_total
    1386            0 :       write(oUnit,*)"time_mmn=",time_mmn
    1387            0 :       write(oUnit,*)"time_interstitial=",time_interstitial
    1388            0 :       write(oUnit,*)"time_ujugaunt=",time_ujugaunt
    1389            0 :       write(oUnit,*)"time_abcof=",time_abcof
    1390            0 :       write(oUnit,*)"time_symm=",time_symm
    1391            0 :       write(oUnit,*)"time_rw=",time_rw
    1392            0 :       write(oUnit,*)"time_film=",time_film
    1393              : 
    1394              : #ifdef CPP_MPI
    1395            0 :       call MPI_BARRIER(mpi_communicatior,ierr)
    1396              : #endif
    1397              : 
    1398              : #ifdef CPP_MPI
    1399            0 :       call MPI_BARRIER(mpi_communicatior,ierr)
    1400              : #endif
    1401            0 :       call timestop("wann_updown")
    1402              : 
    1403            0 :       END SUBROUTINE wann_updown
    1404              :       END MODULE m_wann_updown
        

Generated by: LCOV version 2.0-1