LCOV - code coverage report
Current view: top level - wannier - wann_updown.F (source / functions) Hit Total Coverage
Test: FLEUR test coverage Lines: 0 488 0.0 %
Date: 2024-03-29 04:21:46 Functions: 0 1 0.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------
       2             : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
       3             : ! This file is part of FLEUR and available as free software under the conditions
       4             : ! of the MIT license as expressed in the LICENSE file in more detail.
       5             : !--------------------------------------------------------------------------------
       6             : 
       7             :       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 1.14