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

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

Generated by: LCOV version 1.13