LCOV - code coverage report
Current view: top level - wannier - wannier.F90 (source / functions) Coverage Total Hit
Test: FLEUR test coverage Lines: 54.1 % 798 432
Test Date: 2025-06-14 04:34:23 Functions: 100.0 % 1 1

            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_wannier
       8              :   USE m_juDFT
       9              : #ifdef CPP_MPI
      10              :    use mpi
      11              : #endif
      12              : CONTAINS
      13            4 :   SUBROUTINE wannier(&
      14              :        fmpi,input,kpts,sym,atoms,stars,vacuum,sphhar ,&
      15              :        wann,noco,nococonv,cell,enpara,banddos,sliceplot,vTot,results,&
      16            4 :        eig_idList,l_real,nkpt)
      17              :     !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      18              :     !     Makes necessary for the construction of the wannier functions
      19              :     !     (W90: Yates, Mostofi, Marzari, Souza, Vanderbilt '06 f90 code)
      20              :     !     ab initio preliminaries: constructs the overlaps of the periodic
      21              :     !     parts of the wavefunctions and the projections of the
      22              :     !     wavefunctions
      23              :     !     onto a set of starting wfs, i.e. atomic-like orbitals.
      24              :     !                                                            YM 06
      25              :     !     Mmn(k,b) = <u_{nk}|u_{m(k+b)}>, u being a periodic part
      26              :     !                        of the wavefunction psi_nk
      27              :     !     A_mn^k = <psi_mk|g_n>, where g_n is a trial orbital
      28              :     !     which are written into the files 'WF1.mmn' and 'WF1.amn'
      29              :     !           Marzari Vanderbilt PRB 56,12847(1997)
      30              :     !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      31              :     !     Parallelization, Optionals, Symmetry, Noco&Soc:
      32              :     !     Frank Freimuth
      33              :     !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      34              :     !     The routine reads the bkpts file, which contains the following
      35              :     !     information:
      36              :     !     1st line: nntot (INT) - number of the nearest neighbors for
      37              :     !                             each k-point in the MP mesh
      38              :     !     2-nkpts*nntot lines containing 5 integers i1,i2,i3,i4,i5:
      39              :     !     i1 - the number of the k-point in the kpts file
      40              :     !     i2 - number of the k-point, which is a periodic image of
      41              :     !          k+b in the 1st BZ
      42              :     !     i3-i5 - coordinates of the G-vector, connecting k-point
      43              :     !             i2 with the actual k+b k-point
      44              :     !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      45              :     !     In general, the number of bands for each k-poin t is
      46              :     !     different: N(k), and also differs from the number of bands
      47              :     !     we are interested in: N (for instance 5 d-bands of Cu among
      48              :     !     the 6 s- and d-bands). While matrices Mmn(k) are square
      49              :     !     for each k-point, matrices Mmn(k,b) can be made so after
      50              :     !     defining the maximum number of bands max(N(k)).
      51              :     !     The matrix Amn is non-diagonal by default (N(k)*N).
      52              :     !ccccc ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      53              :     !     Total number of wannier functions: nwfs
      54              :     !     sliceplot%e1s,sliceplot%e2s: lower and upper boundaries of the energy window:
      55              :     !     Needed for sorting by number and sorting by energy.
      56              :     !     Not needed for sorting by index.
      57              :     !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      58              :     !     Extension to case of higher-dimensional Wannier functions
      59              :     !     according to the formalism in PRB 91, 184413 (2015)
      60              :     !     Jan-Philipp Hanke
      61              :     !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      62              :     USE m_types
      63              :     USE m_constants
      64              :     USE m_wann_mmnk_symm
      65              :     USE m_wann_rw_eig
      66              :     USE m_abcof
      67              :     USE m_radfun
      68              :     USE m_radflo
      69              :     USE m_cdnread
      70              :     USE m_wann_mmk0_vac
      71              :     USE m_wann_mmkb_vac
      72              :     USE m_wann_updown
      73              :     USE m_wann_mmk0_sph
      74              :     USE m_wann_ujugaunt
      75              :     USE m_wann_mmkb_sph
      76              :     USE m_wann_projmethod
      77              :     USE m_wann_amn
      78              :     USE m_wann_abinv
      79              :     USE m_wann_kptsrotate
      80              :     USE m_wann_plot
      81              : !    USE m_wann_read_inp !Call wann_read_inp in fleur_init
      82              :     USE m_wann_plot_symm
      83              :     USE m_wann_mmkb_int
      84              :     USE m_wann_postproc
      85              :     USE m_wann_write_mmnk
      86              :     USE m_wann_write_amn
      87              :     USE m_wann_write_nabla
      88              :     USE m_vsoc
      89              :     USE m_wann_write_matrix4
      90              :     USE m_wann_write_matrix5
      91              :     USE m_wann_orbcomp
      92              :     USE m_wann_anglmom
      93              : #ifdef CPP_TOPO
      94              :     USE m_wann_surfcurr
      95              :     USE m_wann_surfcurr_int2
      96              :     USE m_wann_nabla
      97              :     USE m_wann_nabla_vac
      98              :     USE m_wann_soc_to_mom
      99              : #endif
     100              :     USE m_wann_gwf_tools, ONLY : get_index_kq, gwf_plottemplate
     101              :     USE m_wann_gwf_commat
     102              :     USE m_wann_gwf_anglmom
     103              :     USE m_wann_write_mmnk2
     104              :     USE m_wann_uHu
     105              :     USE m_wann_uHu_dmi
     106              :     USE m_eig66_io
     107              : 
     108              :     IMPLICIT NONE
     109              : 
     110              : #ifdef CPP_MPI
     111              :     INTEGER ierr(3)
     112              :     INTEGER cpu_index
     113              :     INTEGER stt(MPI_STATUS_SIZE)
     114              : #endif
     115              : 
     116              : 
     117              :     TYPE(t_mpi),       INTENT(IN) :: fmpi
     118              :     TYPE(t_input),     INTENT(IN) :: input
     119              :     TYPE(t_kpts),      INTENT(IN) :: kpts
     120              :     TYPE(t_sym),       INTENT(IN) :: sym
     121              :     TYPE(t_atoms),     INTENT(IN) :: atoms
     122              :     TYPE(t_stars),     INTENT(IN) :: stars
     123              :     TYPE(t_vacuum),    INTENT(IN) :: vacuum
     124              :     TYPE(t_sphhar),    INTENT(IN) :: sphhar
     125              :      
     126              :     TYPE(t_noco),      INTENT(IN) :: noco
     127              :     TYPE(t_nococonv),  INTENT(IN) :: nococonv
     128              :     TYPE(t_cell),      INTENT(IN) :: cell
     129              :     TYPE(t_enpara),    INTENT(IN) :: enpara
     130              :     TYPE(t_banddos),   INTENT(IN) :: banddos
     131              :     TYPE(t_sliceplot), INTENT(IN) :: sliceplot
     132              :     TYPE(t_potden),    INTENT(IN) :: vTot
     133              :     TYPE(t_results),   INTENT(IN) :: results
     134              :     TYPE(t_wann),      INTENT(INOUT) :: wann
     135              : 
     136              :     LOGICAL, INTENT (in) :: l_real
     137              :     INTEGER, INTENT (in) :: nkpt
     138              :     INTEGER, INTENT (IN) :: eig_idList(:)!wann%nparampts)
     139              : 
     140              :     !ccccccccccccccccc   local variables   cccccccccccccccccccc
     141              :     INTEGER :: lmd,n,nmat,iter,ikpt,ikpt_b, pc
     142              :     INTEGER :: addnoco,loplod,addnoco2,igvm2,eig_id
     143              :     INTEGER :: noccbd,noccbd_b,nn,nkpts,i,jspin,j,l,i_rec,m,nwf,nwfp
     144              :     INTEGER :: jsp_start,jsp_end,nrec,nrec1,nrec_b,nbands,nbands_b
     145              :     INTEGER :: nodeu,noded,n_size,na,n_rank,nbnd,numbands
     146              :     INTEGER :: i1,i2,i3,in,lda
     147            4 :     INTEGER :: n_bands(0:input%neig),nslibd,nslibd_b
     148              :     CHARACTER(len=8) :: dop,iop,name(10)
     149              :     REAL    :: wronk,phase
     150              :     COMPLEX :: c_phase
     151            4 :     REAL    :: eig(input%neig),eig_b(input%neig)
     152              :     REAL    :: efermi
     153              :     LOGICAL :: l_p0,l_bkpts,l_proj,l_amn,l_mmn
     154              : !!! energy window boundaries
     155            4 :     INTEGER, ALLOCATABLE :: innerEig_idList(:)
     156            4 :     REAL,    ALLOCATABLE :: we(:),we_b(:)
     157              : 
     158            4 :     REAL,    ALLOCATABLE :: eigg(:)
     159            4 :     REAL kpoints(nkpt)
     160              : !!! a and b coeff. constructed for each k-point
     161            4 :     COMPLEX, ALLOCATABLE :: acof(:,:,:),acof_b(:,:,:)
     162            4 :     COMPLEX, ALLOCATABLE :: bcof(:,:,:),bcof_b(:,:,:)
     163            4 :     COMPLEX, ALLOCATABLE :: ccof(:,:,:,:),ccof_b(:,:,:,:)
     164              : !!! the parameters for the number of wfs
     165              :     INTEGER :: nwfs
     166              : !!! the potential in the spheres and the vacuum
     167            4 :     REAL, ALLOCATABLE :: vr(:,:,:),vz(:,:,:)
     168              : !!! auxiliary potentials
     169            4 :     COMPLEX, ALLOCATABLE :: vpw(:,:)
     170              : !!! bkpts data
     171              :     INTEGER nntot,ikpt_help
     172            4 :     INTEGER, ALLOCATABLE :: gb(:,:,:),bpt(:,:)
     173              : !!! radial wavefunctions in the muffin-tins and more ...
     174            4 :     REAL,    ALLOCATABLE :: flo(:,:,:,:,:),vso(:,:,:)
     175            4 :     REAL,    ALLOCATABLE :: ff(:,:,:,:,:),gg(:,:,:,:,:)
     176              : 
     177            4 :     REAL     :: uuilon(atoms%nlod,atoms%ntype)
     178            4 :     REAL     :: duilon(atoms%nlod,atoms%ntype)
     179            4 :     REAL     :: ulouilopn(atoms%nlod,atoms%nlod,atoms%ntype)
     180              : !!! the Mmn matrices
     181            4 :     COMPLEX, ALLOCATABLE :: mmnk(:,:,:,:),mmn(:,:,:)
     182            4 :     COMPLEX, ALLOCATABLE :: amn(:,:,:),nablamat(:,:,:,:)
     183            4 :     COMPLEX, ALLOCATABLE :: soctomom(:,:,:,:)
     184            4 :     COMPLEX, ALLOCATABLE :: surfcurr(:,:,:,:)
     185            4 :     COMPLEX, ALLOCATABLE :: socmmn(:,:,:)
     186              :     COMPLEX, ALLOCATABLE :: a(:)
     187            4 :     COMPLEX, ALLOCATABLE :: psiw(:,:,:)
     188            4 :     COMPLEX, ALLOCATABLE :: anglmom(:,:,:,:)
     189            4 :     COMPLEX, ALLOCATABLE :: orbcomp(:,:,:,:,:)
     190              :     !..wf-hamiltonian in real space (hopping in the same unit cell)
     191            4 :     COMPLEX, ALLOCATABLE :: hwfr(:,:),hwfr2(:,:)
     192              :     !      real, allocatable :: ei(:)
     193              :     COMPLEX, ALLOCATABLE :: work(:)
     194              :     REAL,ALLOCATABLE::centers(:,:,:)
     195              :     LOGICAL :: l_file
     196              :     LOGICAL :: l_amn2, l_conjugate
     197              :     CHARACTER(len=3) :: spin12(2)
     198              :     DATA   spin12/'WF1' , 'WF2'/
     199              :     CHARACTER(len=30)  :: task
     200            4 :     INTEGER,ALLOCATABLE::irreduc(:)
     201            4 :     INTEGER,ALLOCATABLE::mapkoper(:)
     202              :     INTEGER :: fullnkpts,kpt,kptibz,kptibz_b,j1,j2,j3,oper,oper_b,k
     203              :     REAL :: bkrot(3),dirfacs(3)
     204              :     INTEGER :: ios,kplot,kplotoper,plotoper,gfcut
     205              :     COMPLEX :: phasust
     206            4 :     INTEGER,ALLOCATABLE::pair_to_do(:,:)
     207            4 :     INTEGER :: ngopr1(atoms%nat)
     208            4 :     INTEGER,ALLOCATABLE::maptopair(:,:,:)
     209              :     INTEGER :: wannierspin,jspin2,jspin7,jspin2_b
     210              :     REAL, ALLOCATABLE :: rwork(:)
     211            4 :     REAL,ALLOCATABLE::kdiff(:,:)
     212            4 :     INTEGER,ALLOCATABLE :: shiftkpt(:,:)
     213              :     INTEGER :: unigrid(6),gfthick
     214            4 :     COMPLEX,ALLOCATABLE::ujug(:,:,:,:),ujdg(:,:,:,:)
     215            4 :     COMPLEX,ALLOCATABLE::djug(:,:,:,:),djdg(:,:,:,:)
     216            4 :     COMPLEX,ALLOCATABLE::ujulog(:,:,:,:,:)
     217            4 :     COMPLEX,ALLOCATABLE::djulog(:,:,:,:,:)
     218            4 :     COMPLEX,ALLOCATABLE::ulojug(:,:,:,:,:)
     219            4 :     COMPLEX,ALLOCATABLE::ulojdg(:,:,:,:,:)
     220            4 :     COMPLEX,ALLOCATABLE::ulojulog(:,:,:,:,:,:)
     221              :     INTEGER :: n_start,n_end,mlotot,mlolotot,err
     222              :     INTEGER :: mlot_d,mlolot_d,ilo,dir,length
     223              :     CHARACTER(len=2) :: spin012(0:2)
     224              :     DATA spin012/'  ', '.1', '.2'/
     225              :     CHARACTER(len=6) :: filename
     226              :     REAL :: arg,hescale
     227              :     COMPLEX :: nsfactor,nsfactor_b,VALUE
     228              :     REAL :: b1(3),b2(3)
     229              :     REAL,PARAMETER :: bohrtocm=0.529177e-8
     230              :     REAL,PARAMETER :: condquant=7.7480917e-5
     231              :     INTEGER :: npotmatfile,ig3,maxvac,irec,imz,ivac,ipot
     232              :     LOGICAL :: l_orbcompinp
     233              :     INTEGER :: num_angl,nbasfcn,nbasfcn_b,noconbasfcn,noconbasfcn_b
     234            4 :     COMPLEX,ALLOCATABLE :: vxy(:,:,:)
     235              : 
     236              : 
     237              :     !---->gwf
     238              : 
     239              :     ! FURTHER VARIABLES
     240              :     REAL :: qpt_i(3),qptb_i(3)
     241            4 :     REAL :: alph_i(atoms%ntype),alphb_i(atoms%ntype)
     242            4 :     REAL :: beta_i(atoms%ntype),betab_i(atoms%ntype)
     243              :     REAL :: theta_i, thetab_i, phi_i, phib_i
     244              :     REAL :: dalph,db1,db2,coph,siph
     245            4 :     REAL :: zero_taual(3,atoms%nat),bqpt(3)
     246            4 :     REAL :: eig_qb(input%neig)
     247              : 
     248            8 :     REAL,ALLOCATABLE :: qdiff(:,:), we_qb(:)
     249              :     REAL,ALLOCATABLE :: energies(:,:,:)
     250            4 :     REAL,ALLOCATABLE :: zero_qdiff(:,:)
     251              : 
     252              : 
     253            4 :     INTEGER,ALLOCATABLE :: irreduc_q(:),mapqoper(:)
     254            4 :     INTEGER,ALLOCATABLE :: shiftqpt(:,:),pair_to_do_q(:,:)
     255            4 :     INTEGER,ALLOCATABLE :: maptopair_q(:,:,:)
     256            4 :     INTEGER,ALLOCATABLE :: gb_q(:,:,:),bpt_q(:,:)
     257              : 
     258              :     INTEGER :: nntot_q = 1
     259              :     INTEGER :: fullnqpts = 1
     260              :     INTEGER :: funit_start = 5000
     261              :     INTEGER :: qptibz, qptibz_b, oper_q, oper_qb
     262              :     INTEGER :: qpt,iqpt_help, iqpt, iqpt_b
     263              :     INTEGER :: nbands_qb, nmat_qb, nslibd_qb, noccbd_qb
     264              :     INTEGER :: sign_q = 1,band_help
     265              :     INTEGER :: doublespin,jspin_b,jspin3,jspin4,jspin5
     266              :     INTEGER :: doublespin_max,nrec5
     267              :     INTEGER :: count_i,count_j
     268              :     INTEGER :: n1,n2,ii,jj,lmplmd
     269              : 
     270              :     COMPLEX :: interchi,vacchi,amnchi
     271              :     COMPLEX :: phasfac,phasfac2
     272              : 
     273            4 :     COMPLEX,ALLOCATABLE :: chi(:)
     274            4 :     COMPLEX,ALLOCATABLE :: acof_qb(:,:,:)
     275            4 :     COMPLEX,ALLOCATABLE :: bcof_qb(:,:,:)
     276            4 :     COMPLEX,ALLOCATABLE :: ccof_qb(:,:,:,:)
     277            4 :     COMPLEX,ALLOCATABLE :: mmnk_q(:,:,:,:)
     278            4 :     COMPLEX,ALLOCATABLE :: m_int(:,:,:,:)
     279            4 :     COMPLEX,ALLOCATABLE :: m_sph(:,:,:,:)
     280            4 :     COMPLEX,ALLOCATABLE :: m_vac(:,:,:,:)
     281            4 :     COMPLEX,ALLOCATABLE :: ujug_q(:,:,:,:),ujdg_q(:,:,:,:)
     282            4 :     COMPLEX,ALLOCATABLE :: djug_q(:,:,:,:),djdg_q(:,:,:,:)
     283            4 :     COMPLEX,ALLOCATABLE :: ujulog_q(:,:,:,:,:)
     284            4 :     COMPLEX,ALLOCATABLE :: djulog_q(:,:,:,:,:)
     285            4 :     COMPLEX,ALLOCATABLE :: ulojug_q(:,:,:,:,:)
     286            4 :     COMPLEX,ALLOCATABLE :: ulojdg_q(:,:,:,:,:)
     287            4 :     COMPLEX,ALLOCATABLE :: ulojulog_q(:,:,:,:,:,:)
     288              : 
     289              :     CHARACTER(len=30) fname
     290              : 
     291              :     LOGICAL :: l_bqpts,l_gwf,l_nochi
     292              : 
     293            4 :     TYPE(t_usdus) :: usdus
     294            4 :     TYPE(t_mat)   :: zMat, zzMat, zMat_b, zMat_qb
     295            4 :     TYPE(t_lapw)  :: lapw, lapw_b, lapw_qb
     296           68 :     TYPE(t_wann)  :: wannTemp
     297              : 
     298            4 :     eig_id = eig_idList(1)
     299              : 
     300              :     !----<gwf
     301              : 
     302            4 :     lmplmd = (atoms%lmaxd*(atoms%lmaxd+2)* (atoms%lmaxd*(atoms%lmaxd+2)+3))/2
     303              : 
     304              : 
     305              :     !-----initializations
     306           12 :     ngopr1(:)=1
     307           36 :     zero_taual = 0.0
     308              : 
     309            4 :     hescale=sqrt(tpi_const*condquant/bohrtocm/cell%omtil)
     310              : 
     311            4 :     CALL timestart("Wannier total")
     312              : 
     313            4 :     l_p0 = .FALSE.
     314            4 :     IF (fmpi%irank.EQ.0) l_p0 = .TRUE.
     315              : 
     316            4 :     lmd = atoms%lmaxd*(atoms%lmaxd+2)
     317              : 
     318              : !!!   should be changed in case the windows are really used
     319            4 :     nkpts = nkpt
     320              : 
     321              :     ! do we have to construct GWF ?
     322              :     l_gwf = .FALSE.
     323            4 :     l_gwf = wann%l_sgwf.OR.wann%l_socgwf
     324              : 
     325            4 :     l_nochi = .FALSE.
     326            4 :     INQUIRE(file='l_nochi',exist=l_nochi)
     327            4 :     IF(l_gwf.AND.l_p0) WRITE(*,*)'disable chi trafo: ',l_nochi
     328              : 
     329            4 :     IF(l_gwf.AND.l_p0) CALL gwf_plottemplate()
     330            4 :     ALLOCATE( chi(atoms%ntype) )
     331              : 
     332              :     !-----read the input file to determine what to do
     333            4 :     wann%atomlist_num=atoms%nat
     334            4 :     wann%oc_num_orbs=atoms%nat
     335              : 
     336              : !    CALL wann_read_inp(input,l_p0,wann) !Call wann_read_inp in fleur_init
     337              : 
     338              : 
     339              :     !-----input file for orbital decomposition
     340            4 :     IF(wann%l_orbcomp.OR.wann%l_orbcomprs)THEN
     341            0 :        INQUIRE(file='orbcomp_inp',exist=l_orbcompinp)
     342            0 :        IF(l_orbcompinp)THEN
     343            0 :           OPEN(159,file='orbcomp_inp')
     344            0 :           READ(159,*)wann%oc_num_orbs,wann%l_oc_f
     345            0 :           ALLOCATE(wann%oc_orbs(wann%oc_num_orbs))
     346            0 :           DO n=1,wann%oc_num_orbs
     347            0 :              READ(159,*)wann%oc_orbs(n)
     348              :           ENDDO
     349            0 :           CLOSE(159)
     350              :        ELSE !default is all atoms including f
     351            0 :           wann%oc_num_orbs=atoms%nat
     352            0 :           wann%l_oc_f=.TRUE.
     353            0 :           ALLOCATE(wann%oc_orbs(wann%oc_num_orbs))
     354            0 :           DO n=1,wann%oc_num_orbs
     355            4 :              wann%oc_orbs(n)=n
     356              :           ENDDO
     357              :        ENDIF
     358              :     ENDIF
     359              : 
     360            4 :     IF(wann%l_updown)THEN
     361              :        CALL wann_updown(&
     362              :             wann,fmpi,input,kpts,sym,atoms,stars,vacuum,sphhar, &
     363              :              noco,nococonv,cell,vTot,&
     364              :             enpara,eig_idList(1),l_real,&
     365              :             fmpi%mpi_comm,atoms%l_dulo,noco%l_noco,noco%l_ss,&
     366              :             atoms%lmaxd,atoms%ntype,input%neig,atoms%nat,sym%nop,&
     367              :             lapw%dim_nvd(),input%jspins,atoms%llod,&
     368              :             atoms%nlod,atoms%ntype,cell%omtil,atoms%nlo,atoms%llo,&
     369              :             atoms%lapw_l,sym%invtab,sym%mrot,sym%ngopr,atoms%neq,&
     370              :             atoms%lmax,sym%invsat,sym%invsatnr,nkpt,atoms%taual,&
     371              :             atoms%rmt,cell%amat,cell%bmat,cell%bbmat,nococonv%alph,&
     372              :             nococonv%beta,nococonv%qss,&                    ! TODO: adapt if needed&
     373              :             stars%sk2,stars%phi2 ,fmpi%irank,fmpi%isize,&
     374              :             stars%ng3,&
     375              :             vacuum%nmzxyd,vacuum%nmzd,atoms%jmtd,sphhar%nlhd,stars%ng3,&
     376              :             vacuum%nvac,sym%invs,sym%invs2,input%film,sphhar%nlh,&
     377              :             atoms%jri,sphhar%ntypsd,sym%ntypsy,input%jspins,nkpt,&
     378              :             atoms%dx,stars%ng2,atoms%rmsh,sliceplot%e1s,sliceplot%e2s,&
     379              :             atoms%ulo_der,stars%ustep,stars%ig,stars%mx1,&
     380              :             stars%mx2,stars%mx3,stars%rgphs,&
     381              :             sliceplot%slice,sliceplot%kk,sliceplot%nnne,cell%z1,&
     382              :             lapw%dim_nv2d(),vacuum%nmzxy,vacuum%nmz,vacuum%delz,&
     383              :             stars%ig2,cell%area,sym%tau,atoms%zatom,stars%ng2,sym%nop2,&
     384              :             cell%volint,sym%symor,atoms%pos,results%ef,noco%l_soc,&
     385              :             sphhar%memd,atoms%lnonsph,sphhar%clnu,lmplmd,&
     386              :             sphhar%mlh,sphhar%nmem,sphhar%llh,atoms%lo1l,&
     387            0 :             nococonv%theta,nococonv%phi)
     388            0 :        if(wann%l_stopupdown)then
     389            0 :        DO pc = 1, wann%nparampts
     390            0 :           CALL close_eig(eig_idList(pc))
     391              :        END DO
     392              : 
     393            0 :        CALL juDFT_end("updown done",fmpi%irank)
     394              :        endif
     395              :     ENDIF
     396              : 
     397              :     !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     398              :     !  modern theory of orbital magnetization from Wannier functions
     399              :     !  Jan-Philipp Hanke
     400              :     !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     401            4 :     IF(wann%l_matrixuHu)THEN
     402            0 :        wannTemp = wann
     403              :        CALL wann_uHu(&
     404              :             stars,vacuum,atoms,sphhar,input,kpts,sym,fmpi,&
     405              :             banddos ,noco,nococonv,cell,vTot,wannTemp,enpara,eig_idList,&
     406              :             l_real,atoms%l_dulo,noco%l_noco,noco%l_ss,atoms%lmaxd,&
     407              :             atoms%ntype,input%neig,atoms%nat,sym%nop,lapw%dim_nvd(),&
     408              :             input%jspins,atoms%llod,atoms%nlod,&
     409              :             atoms%ntype,cell%omtil,atoms%nlo,atoms%llo,&
     410              :             atoms%lapw_l,sym%invtab,sym%mrot,sym%ngopr,atoms%neq,&
     411              :             atoms%lmax,sym%invsat,sym%invsatnr,nkpt,atoms%taual,&
     412              :             atoms%rmt,cell%amat,cell%bmat,cell%bbmat,nococonv%alph,&
     413              :             nococonv%beta,nococonv%qss,stars%sk2,stars%phi2,&
     414              :             fmpi%irank,&
     415              :             fmpi%isize,stars%ng3,vacuum%nmzxyd,vacuum%nmzd,atoms%jmtd,&
     416              :             sphhar%nlhd,stars%ng3,vacuum%nvac,sym%invs,sym%invs2,&
     417              :             input%film,sphhar%nlh,atoms%jri,sphhar%ntypsd,sym%ntypsy,&
     418              :             input%jspins,nkpt,atoms%dx,stars%ng2,atoms%rmsh,&
     419              :             sliceplot%e1s,sliceplot%e2s,atoms%ulo_der,stars%ustep,&
     420              :             stars%ig,stars%mx1,stars%mx2,stars%mx3,&
     421              :             stars%rgphs,sliceplot%slice,&
     422              :             sliceplot%kk,sliceplot%nnne,cell%z1,lapw%dim_nv2d(),&
     423              :             vacuum%nmzxy,vacuum%nmz,vacuum%delz,stars%ig2,&
     424              :             cell%area,sym%tau,atoms%zatom,stars%ng2,stars%kv2,sym%nop2,&
     425              :             cell%volint,sym%symor,atoms%pos,results%ef,noco%l_soc,&
     426              :             sphhar%memd,atoms%lnonsph,sphhar%clnu,lmplmd,&
     427              :             sphhar%mlh,sphhar%nmem,sphhar%llh,atoms%lo1l,&
     428              :             nococonv%theta,nococonv%phi,&
     429              :             wann%l_ms,wann%l_sgwf,wann%l_socgwf,wann%aux_latt_const,&
     430              :             wann%param_file,wann%param_vec,wann%nparampts,&
     431            0 :             wann%param_alpha,wann%l_dim)
     432              : 
     433              : 
     434            0 :        if(wann%l_stopuhu) then
     435              : 
     436            0 :         DO pc = 1, wann%nparampts
     437            0 :           CALL close_eig(eig_idList(pc))
     438              :         END DO
     439              : 
     440            0 :         CALL juDFT_end("wann_uHu done",fmpi%irank)
     441              :        endif
     442              :     ENDIF
     443              : 
     444              :     !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     445              :     !  modern theory of DMI from higher-dimensional Wannier functions
     446              :     !  Jan-Philipp Hanke
     447              :     !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     448            4 :     IF(wann%l_matrixuHu_dmi)THEN
     449            0 :        wannTemp = wann
     450              :        CALL wann_uHu_dmi(&
     451              :             stars,vacuum,atoms,sphhar,input,kpts,sym,fmpi,&
     452              :             banddos ,noco,nococonv,cell,vTot,wannTemp,enpara,eig_idList,&
     453              :             l_real,atoms%l_dulo,noco%l_noco,noco%l_ss,atoms%lmaxd,&
     454              :             atoms%ntype,input%neig,atoms%nat,sym%nop,lapw%dim_nvd(),&
     455              :             input%jspins,lapw%dim_nbasfcn(),atoms%llod,atoms%nlod,&
     456              :             atoms%ntype,cell%omtil,atoms%nlo,atoms%llo,&
     457              :             atoms%lapw_l,sym%invtab,sym%mrot,sym%ngopr,atoms%neq,&
     458              :             atoms%lmax,sym%invsat,sym%invsatnr,nkpt,atoms%taual,&
     459              :             atoms%rmt,cell%amat,cell%bmat,cell%bbmat,nococonv%alph,&
     460              :             nococonv%beta,nococonv%qss,stars%sk2,stars%phi2,&
     461              :             fmpi%irank,&
     462              :             fmpi%isize,stars%ng3,vacuum%nmzxyd,vacuum%nmzd,atoms%jmtd,&
     463              :             sphhar%nlhd,stars%ng3,vacuum%nvac,sym%invs,sym%invs2,&
     464              :             input%film,sphhar%nlh,atoms%jri,sphhar%ntypsd,sym%ntypsy,&
     465              :             input%jspins,nkpt,atoms%dx,stars%ng2,atoms%rmsh,&
     466              :             sliceplot%e1s,sliceplot%e2s,atoms%ulo_der,stars%ustep,&
     467              :             stars%ig,stars%mx1,stars%mx2,stars%mx3,&
     468              :             stars%rgphs,sliceplot%slice,&
     469              :             sliceplot%kk,sliceplot%nnne,cell%z1,lapw%dim_nv2d(),&
     470              :             vacuum%nmzxy,vacuum%nmz,vacuum%delz,stars%ig2,&
     471              :             cell%area,sym%tau,atoms%zatom,stars%ng2,stars%kv2,sym%nop2,&
     472              :             cell%volint,sym%symor,atoms%pos,results%ef,noco%l_soc,&
     473              :             sphhar%memd,atoms%lnonsph,sphhar%clnu,lmplmd,&
     474              :             sphhar%mlh,sphhar%nmem,sphhar%llh,atoms%lo1l,&
     475              :             nococonv%theta,nococonv%phi,&
     476              :             wann%l_ms,wann%l_sgwf,wann%l_socgwf,wann%aux_latt_const,&
     477              :             wann%param_file,wann%param_vec,wann%nparampts,&
     478            0 :             wann%param_alpha,wann%l_dim,l_nochi)
     479              : 
     480            0 :        DO pc = 1, wann%nparampts
     481            0 :           CALL close_eig(eig_idList(pc))
     482              :        END DO
     483              : 
     484            0 :        CALL juDFT_end("wann_uHu dmi done",fmpi%irank)
     485              :     ENDIF
     486              : 
     487            4 :     IF(wann%l_byenergy.AND.wann%l_byindex) CALL juDFT_error&
     488            0 :          ("byenergy.and.byindex",calledby ="wannier")
     489            4 :     IF(wann%l_byenergy.AND.wann%l_bynumber) CALL juDFT_error&
     490            0 :          ("byenergy.and.bynumber",calledby ="wannier")
     491            4 :     IF(wann%l_bynumber.AND.wann%l_byindex) CALL juDFT_error&
     492            0 :          ("bynumber.and.byindex",calledby ="wannier")
     493            4 :     IF(.NOT.(wann%l_bynumber.OR.wann%l_byindex.OR.wann%l_byenergy))&
     494            0 :          CALL juDFT_error("no rule to sort bands",calledby ="wannier")
     495              : 
     496              : 
     497            4 :     efermi=results%ef
     498            4 :     IF(.NOT.wann%l_fermi)efermi=0.0
     499              : 
     500              : #ifdef CPP_MPI
     501            4 :     CALL MPI_BARRIER(fmpi%mpi_comm,ierr(1))
     502              : #endif
     503              : 
     504              :     !**************************************************************
     505              :     !   for bzsym=.true.: determine mapping between kpts and w90kpts
     506              :     !**************************************************************
     507            4 :     IF (wann%l_bzsym) THEN
     508            0 :        l_file=.FALSE.
     509            0 :        INQUIRE(file='w90kpts',exist=l_file)
     510            0 :        IF(.NOT.l_file)  CALL juDFT_error&
     511            0 :             ("w90kpts not found, needed if bzsym",calledby ="wannier")
     512            0 :        OPEN(412,file='w90kpts',form='formatted')
     513            0 :        READ(412,*)fullnkpts
     514            0 :        CLOSE(412)
     515            0 :        IF(l_p0)PRINT*,"fullnkpts=",fullnkpts
     516            0 :        IF(fullnkpts<nkpts) CALL juDFT_error("fullnkpts.lt.nkpts"&
     517            0 :             ,calledby ="wannier")
     518            0 :        ALLOCATE(irreduc(fullnkpts),mapkoper(fullnkpts))
     519            0 :        ALLOCATE(shiftkpt(3,fullnkpts))
     520            0 :        l_file=.FALSE.
     521            0 :        INQUIRE(file='kptsmap',exist=l_file)
     522            0 :        IF(.NOT.l_file)  CALL juDFT_error&
     523            0 :             ("kptsmap not found, needed if bzsym",calledby ="wannier")
     524            0 :        OPEN(713,file='kptsmap')
     525            0 :        DO i=1,fullnkpts
     526            0 :           READ(713,*)kpt,irreduc(i),mapkoper(i),shiftkpt(:,i)
     527            0 :           IF(kpt/=i) CALL juDFT_error("kpt.ne.i",calledby ="wannier")
     528            0 :           IF(l_p0)PRINT*,i,irreduc(i),mapkoper(i)
     529              :        ENDDO
     530            0 :        CLOSE(713)
     531            0 :        IF(MAXVAL(irreduc(:))/=nkpts) CALL juDFT_error&
     532            0 :             ("max(irreduc(:))/=nkpts",calledby ="wannier")
     533              :     ELSE
     534            4 :        fullnkpts=nkpts
     535            4 :        ALLOCATE(irreduc(fullnkpts),mapkoper(fullnkpts))
     536            4 :        ALLOCATE(shiftkpt(3,fullnkpts))
     537              :     ENDIF
     538              : 
     539              : 
     540            4 :     IF(l_gwf) fullnqpts = wann%nparampts
     541              : 
     542              : 
     543            4 :     nrec = 0
     544            4 :     IF(l_p0)THEN
     545            2 :        WRITE (*,*) 'fermi energy:',efermi
     546            2 :        WRITE (*,*) 'emin,emax=',sliceplot%e1s,sliceplot%e2s
     547              :     ENDIF
     548              : 
     549              :     IF((.NOT.wann%l_matrixmmn).AND.(.NOT.wann%l_wann_plot).AND.&
     550              :          (.NOT.wann%l_matrixamn).AND.(.NOT.wann%l_projmethod).AND.&
     551              :          (.NOT.wann%l_bestproj).AND.(.NOT.wann%l_nabla).AND.&
     552              :          (.NOT.wann%l_mmn0).AND.(.NOT.wann%l_surfcurr).AND.&
     553              :          (.NOT.wann%l_offdiposop).AND.(.NOT.wann%l_anglmom).AND.&
     554            4 :          (.NOT.wann%l_orbcomp).AND.(.NOT.wann%l_perturb) .AND.&
     555              :          (.NOT.wann%l_finishgwf) ) GOTO 1911
     556              : 
     557              :     !**********************************************************
     558              :     !cccccccccccccc   read in the bkpts file  ccccccccccccccccc
     559              :     !**********************************************************
     560            2 :     IF (wann%l_matrixmmn) THEN ! for Omega functional minimization
     561            2 :        l_bkpts = .FALSE.
     562            2 :        INQUIRE (file='bkpts',exist=l_bkpts)
     563            2 :        IF (.NOT.l_bkpts)  CALL juDFT_error("need bkpts for matrixmmn"&
     564            0 :             ,calledby ="wannier")
     565            2 :        OPEN (202,file='bkpts',form='formatted',status='old')
     566            2 :        REWIND (202)
     567            2 :        READ (202,'(i4)') nntot
     568            2 :        IF(l_p0)THEN
     569            1 :           WRITE (*,*) 'nntot=',nntot
     570            1 :           WRITE(*,*) 'fullnkpts=',fullnkpts
     571            1 :           WRITE(*,*) 'nkpts=',nkpts
     572              :        ENDIF
     573            6 :        ALLOCATE ( gb(1:3,1:nntot,1:fullnkpts),bpt(1:nntot,1:fullnkpts))
     574           18 :        DO ikpt=1,fullnkpts
     575          146 :           DO nn=1,nntot
     576              :              READ (202,'(2i6,3x,3i4)')&
     577          128 :                   ikpt_help,bpt(nn,ikpt),(gb(i,nn,ikpt),i=1,3)
     578          128 :              IF (ikpt/=ikpt_help)  CALL juDFT_error("ikpt.ne.ikpt_help"&
     579            0 :                   ,calledby ="wannier")
     580          128 :              IF (bpt(nn,ikpt)>fullnkpts) CALL juDFT_error("bpt.gt.fullnkpts"&
     581           16 :                   ,calledby ="wannier")
     582              :           ENDDO
     583              :        ENDDO
     584            2 :        CLOSE (202)
     585            2 :        ALLOCATE(kdiff(3,nntot))
     586              :     ENDIF
     587              : 
     588              :     !**********************************************************
     589              :     !cccccccccccccc   read in the bqpts file  ccccccccccccccccc
     590              :     !**********************************************************
     591            2 :     IF ((wann%l_matrixmmn).AND.(l_gwf.OR.wann%l_ms)) THEN
     592            0 :        l_bqpts = .FALSE.
     593            0 :        INQUIRE (file='bqpts',exist=l_bqpts)
     594            0 :        IF (.NOT.l_bqpts)  CALL juDFT_error("need bqpts for matrixmmn"&
     595            0 :             ,calledby ="wannier")
     596            0 :        OPEN (202,file='bqpts',form='formatted',status='old')
     597            0 :        REWIND (202)
     598            0 :        READ (202,'(i4)') nntot_q
     599            0 :        IF(l_p0)THEN
     600            0 :           WRITE (*,*) 'nntot_q=',nntot_q
     601            0 :           WRITE(*,*) 'fullnqpts=',fullnqpts
     602              :        ENDIF
     603              :        ALLOCATE ( gb_q(1:3,1:nntot_q,1:fullnqpts),&
     604            0 :             bpt_q(1:nntot_q,1:fullnqpts))
     605            0 :        DO iqpt=1,fullnqpts
     606            0 :           DO nn=1,nntot_q
     607              :              READ (202,'(2i6,3x,3i4)')&
     608            0 :                   iqpt_help,bpt_q(nn,iqpt),(gb_q(i,nn,iqpt),i=1,3)
     609            0 :              IF (iqpt/=iqpt_help)  CALL juDFT_error("iqpt.ne.iqpt_help"&
     610            0 :                   ,calledby ="wannier")
     611            0 :              IF (bpt_q(nn,iqpt)>fullnqpts)&
     612            0 :                   CALL juDFT_error("bpt_q.gt.fullnqpts",calledby ="wannier")
     613              :           ENDDO
     614              :        ENDDO
     615            0 :        CLOSE (202)
     616            0 :        ALLOCATE(qdiff(3,nntot_q))
     617            0 :        ALLOCATE(zero_qdiff(3,nntot_q))
     618            2 :        zero_qdiff=0.0
     619              :     ENDIF
     620              : 
     621              : 
     622              :     ! when treating gen. WF for spin spirals, the Brillouin zone
     623              :     ! of q-points is twice as large compared to k-BZ. Thus,
     624              :     ! the G-vectors connecting neighbors across the boundary
     625              :     ! need to be doubled
     626            2 :     IF(wann%l_sgwf) gb_q = 2*gb_q
     627            2 :     IF(wann%l_socgwf) gb_q = 2*gb_q
     628              : 
     629            2 :     IF(wann%l_finishgwf) GOTO 9110
     630              :     !********************************************************
     631              :     !      find symmetry-related elements in mmkb
     632              :     !********************************************************
     633            2 :     IF(wann%l_matrixmmn)THEN
     634            6 :        ALLOCATE(maptopair(3,fullnkpts,nntot))
     635            2 :        ALLOCATE(pair_to_do(fullnkpts,nntot))
     636              :        CALL wann_mmnk_symm(input,kpts,&
     637              :             fullnkpts,nntot,bpt,gb,wann%l_bzsym,&
     638              :             irreduc,mapkoper,l_p0,input%film,sym%nop,sym%invtab,sym%mrot,&
     639              :             sym%tau,&
     640            2 :             pair_to_do,maptopair,kdiff,.FALSE.,wann%param_file)
     641              :     ENDIF
     642              : 
     643              :     ! do the same for q-points to construct GWFs
     644            2 :     IF(wann%l_matrixmmn.AND.l_gwf)THEN
     645            0 :        ALLOCATE(maptopair_q(3,fullnqpts,nntot_q))
     646            0 :        ALLOCATE(pair_to_do_q(fullnqpts,nntot_q))
     647              :        CALL wann_mmnk_symm(input,kpts,&
     648              :             fullnqpts,nntot_q,bpt_q,gb_q,wann%l_bzsym,&
     649              :             irreduc_q,mapqoper,l_p0,.FALSE.,1,sym%invtab(1),&
     650              :             sym%mrot(:,:,1),sym%tau,&
     651            0 :             pair_to_do_q,maptopair_q,qdiff,.TRUE.,wann%param_file)
     652              :     ENDIF
     653              : 
     654              : 
     655              :     !*********************************************************
     656              :     !ccccccccccccccc   initialize the potential   cccccccccccc
     657              :     !*********************************************************
     658              : 
     659            6 :     ALLOCATE ( vr(atoms%jmtd,atoms%ntype,input%jspins) )
     660            4 :     ALLOCATE ( vso(atoms%jmtd,atoms%nat,2) )
     661            2 :     ALLOCATE ( vz(vacuum%nmzd,2,4) )
     662              :     
     663         4026 :     vz = 0.0
     664              : 
     665            2 :     IF(input%film)THEN 
     666            0 :       vz(:,:,:SIZE(vTot%vac,4)) = REAL(vTot%vac(:,1,:,:))
     667            0 :       IF (SIZE(vTot%vac,4)==3) vz(:,:,4) = AIMAG(vTot%vac(:,1,:,3))
     668              :     END IF
     669              : 
     670            4 :     DO jspin = 1,input%jspins
     671            6 :        DO n = 1, atoms%ntype
     672          946 :           DO j = 1,atoms%jri(n)
     673          944 :              vr(j,n,jspin) = vTot%mt(j,0,n,jspin)
     674              :           ENDDO
     675              :        ENDDO
     676              :     ENDDO
     677              : 
     678            2 :     IF(wann%l_soctomom)THEN
     679            0 :        CALL vsoc(input,atoms,vr,enpara%el0,.TRUE., vso)
     680              :     ENDIF
     681              : 
     682            2 :     IF(noco%l_noco.AND.input%film)THEN
     683            0 :        npotmatfile=25
     684            0 :        ALLOCATE(vpw(stars%ng3,1))
     685            0 :        ALLOCATE( vxy(vacuum%nmzxyd,stars%ng2-1,2) )
     686              : 
     687              :        OPEN (npotmatfile,FILE='potmat',FORM='unformatted',&
     688            0 :             STATUS='old')
     689            0 :        READ (npotmatfile) (vpw(ig3,1),ig3=1,stars%ng3)
     690            0 :        READ (npotmatfile) (vpw(ig3,1),ig3=1,stars%ng3)
     691            0 :        READ (npotmatfile) (vpw(ig3,1),ig3=1,stars%ng3)
     692            0 :        maxvac=2
     693            0 :        DO ivac = 1,maxvac
     694              :           !--->       if the two vacuua are equivalent, the potential file has to
     695              :           !--->       be backspaced, because the potential is the same at both
     696              :           !--->       surfaces of the film
     697            0 :           IF ((ivac.EQ.2) .AND. (vacuum%nvac.EQ.1)) THEN
     698            0 :              DO irec = 1,4
     699            0 :                 BACKSPACE (npotmatfile)
     700              :              ENDDO
     701              :           ENDIF
     702              :           !--->       load the non-warping part of the potential
     703              :           READ (npotmatfile)&
     704            0 :                ((vz(imz,ivac,ipot),imz=1,vacuum%nmzd),ipot=1,4)
     705              : 
     706            0 :              DO ipot = 1,3
     707            0 :                 READ (npotmatfile)((vxy(imz,igvm2,ivac),&
     708            0 :                      imz=1,vacuum%nmzxy),igvm2=1,stars%ng2-1)
     709              :              ENDDO
     710              :        ENDDO
     711            0 :        CLOSE (npotmatfile)
     712              :     ENDIF
     713              : 
     714              :     !ccccccccccccccc   end of the potential part  ccccccccccc
     715            2 :     wannierspin=input%jspins
     716            2 :     IF(noco%l_soc) wannierspin=2
     717              : 
     718              : 
     719            8 :     ALLOCATE ( ff(atoms%ntype,atoms%jmtd,2,0:atoms%lmaxd,2) )
     720            4 :     ALLOCATE ( gg(atoms%ntype,atoms%jmtd,2,0:atoms%lmaxd,2) )
     721            4 :     ALLOCATE ( usdus%us(0:atoms%lmaxd,atoms%ntype,2) )
     722            2 :     ALLOCATE ( usdus%uds(0:atoms%lmaxd,atoms%ntype,2) )
     723            2 :     ALLOCATE ( usdus%dus(0:atoms%lmaxd,atoms%ntype,2) )
     724            2 :     ALLOCATE ( usdus%duds(0:atoms%lmaxd,atoms%ntype,2) )
     725            2 :     ALLOCATE ( usdus%ddn(0:atoms%lmaxd,atoms%ntype,2) )
     726            4 :     ALLOCATE ( usdus%ulos(atoms%nlod,atoms%ntype,2) )
     727            2 :     ALLOCATE ( usdus%dulos(atoms%nlod,atoms%ntype,2) )
     728            2 :     ALLOCATE ( usdus%uulon(atoms%nlod,atoms%ntype,2) )
     729            2 :     ALLOCATE ( usdus%dulon(atoms%nlod,atoms%ntype,2) )
     730            8 :     ALLOCATE ( usdus%uloulopn(atoms%nlod,atoms%nlod,atoms%ntype,2) )
     731              : 
     732            2 :     IF(l_gwf.AND..NOT.(wann%l_wann_plot)) THEN
     733              :        doublespin_max=4!2
     734              :     ELSE
     735            2 :        doublespin_max=wannierspin
     736              :     ENDIF
     737              : 
     738              :     !*****************************************************************c
     739              :     !                         START Q LOOP                            c
     740              :     !   standard functionality of code for fullnqpts = nntot_q = 1    c
     741              :     !        and wann%l_ms = wann%l_sgwf = wann%l_socgwf = F          c
     742              :     !*****************************************************************c
     743            4 :     DO iqpt = 1,fullnqpts  ! loop by q-points starts
     744              : 
     745            2 :        ALLOCATE(innerEig_idList(nntot_q))
     746              : 
     747            2 :        qptibz=iqpt
     748            2 :        IF(wann%l_bzsym .AND. l_gwf) qptibz=irreduc_q(iqpt)
     749              :        IF(wann%l_bzsym .AND. l_gwf) oper_q=mapqoper(iqpt)
     750              : 
     751            8 :        qpt_i = nococonv%qss
     752            4 :        alph_i = nococonv%alph
     753            4 :        beta_i = nococonv%beta
     754            2 :        theta_i = nococonv%theta
     755            2 :        phi_i = nococonv%phi
     756            2 :        IF(wann%l_sgwf.OR.wann%l_ms) THEN
     757            0 :           qpt_i(:) = wann%param_vec(:,qptibz)
     758            0 :           alph_i(:) = wann%param_alpha(:,qptibz)
     759            2 :        ELSEIF(wann%l_socgwf) THEN
     760            0 :           IF(wann%l_dim(2)) phi_i = tpi_const*wann%param_vec(2,qptibz)
     761            0 :           IF(wann%l_dim(3)) theta_i = tpi_const*wann%param_vec(3,qptibz)
     762              :        ENDIF
     763              : 
     764            2 :        IF (l_gwf) THEN
     765            0 :           IF(wann%l_matrixmmn)THEN
     766            0 :              DO iqpt_b=1,nntot_q
     767              : 
     768            0 :                 innerEig_idList(iqpt_b) = eig_idList(bpt_q(iqpt_b,iqpt))
     769              : 
     770              :                 !            WRITE(fending,'("_",i4.4)')bpt_q(iqpt_b,iqpt)
     771              :                 !            innerEig_idList(iqpt_b)=open_eig(fmpi%mpi_comm,
     772              :                 !     +                  lapw%dim_nbasfcn(),input%neig,
     773              :                 !     +                  nkpts,wannierspin,atoms%lmaxd,
     774              :                 !     +                  atoms%nlod,atoms%ntype,atoms%nlotot,
     775              :                 !     +                  noco%l_noco,.FALSE.,l_real,noco%l_soc,.FALSE.,
     776              :                 !     +                  mpi%n_size,filename=trim(fstart)//fending,
     777              :                 !     +                  layers=banddos%layers,nstars=banddos%nstars,
     778              :                 !     +                  ncored=29,nsld=atoms%nat,
     779              :                 !     +                  nat=atoms%nat,l_dos=banddos%dos.OR.input%cdinf,
     780              :                 !     +                  l_mcd=banddos%l_mcd,l_orb=banddos%l_orb)
     781              : 
     782              :              ENDDO
     783              :           ENDIF
     784              : 
     785            0 :           eig_id = eig_idList(qptibz)
     786              : 
     787              :           !        WRITE(fending,'("_",i4.4)')qptibz
     788              :           !        eig_id=open_eig(fmpi%mpi_comm,lapw%dim_nbasfcn(),input%neig,
     789              :           !     +                  nkpts,wannierspin,atoms%lmaxd,
     790              :           !     +                  atoms%nlod,atoms%ntype,atoms%nlotot,
     791              :           !     +                  noco%l_noco,.FALSE.,l_real,noco%l_soc,.FALSE.,
     792              :           !     +                  mpi%n_size,filename=trim(fstart)//fending,
     793              :           !     +                  layers=banddos%layers,nstars=banddos%nstars,
     794              :           !     +                  ncored=29,nsld=atoms%nat,
     795              :           !     +                  nat=atoms%nat,l_dos=banddos%dos.OR.input%cdinf,
     796              :           !     +                  l_mcd=banddos%l_mcd,l_orb=banddos%l_orb)
     797              : 
     798            2 :        ELSEIF(wann%l_ms) THEN
     799              : 
     800            0 :           eig_id = eig_idList(qptibz)
     801              : 
     802              :           !        WRITE(fending,'("_",i4.4)')qptibz
     803              :           !        eig_id=open_eig(fmpi%mpi_comm,lapw%dim_nbasfcn(),input%neig,
     804              :           !     +                  nkpts,wannierspin,atoms%lmaxd,
     805              :           !     +                  atoms%nlod,atoms%ntype,atoms%nlotot,
     806              :           !     +                  noco%l_noco,.FALSE.,l_real,noco%l_soc,.FALSE.,
     807              :           !     +                  mpi%n_size,filename=trim(fstart)//fending,
     808              :           !     +                  layers=banddos%layers,nstars=banddos%nstars,
     809              :           !     +                  ncored=29,nsld=atoms%nat,
     810              :           !     +                  nat=atoms%nat,l_dos=banddos%dos.OR.input%cdinf,
     811              :           !     +                  l_mcd=banddos%l_mcd,l_orb=banddos%l_orb)
     812              : 
     813              :        ENDIF ! l_gwf.or.wann%l_ms
     814            2 :        nrec=0
     815            2 :        nrec_b=0
     816              : 
     817              : 
     818              :        !****************************************************
     819              :        ! cycle by spins starts!
     820              :        !****************************************************
     821            4 :        DO doublespin=1,doublespin_max   ! cycle by spins
     822              : 
     823            2 :           jspin=MOD(doublespin+1,2)+1
     824            2 :           jspin_b=jspin
     825            2 :           IF(doublespin.EQ.3) jspin_b=2
     826            2 :           IF(doublespin.EQ.4) jspin_b=1
     827              : 
     828            2 :           nrec_b = nrec
     829              : 
     830            2 :           IF(.NOT.noco%l_noco) THEN
     831            2 :              nrec = (jspin-1)*nkpts
     832            2 :              nrec_b = (jspin_b-1)*nkpts
     833              :           ENDIF
     834              : 
     835              :           ! spin-dependent sign of the q-dependent phase
     836              :           ! in the generalized Bloch theorem
     837              :           ! -1: spin up, +1: spin down
     838            2 :           sign_q = -sign_q
     839              : 
     840              :           !...read number of bands and wannier functions from file proj
     841              : 
     842              :           !..reading the proj.1 / proj.2 / proj file
     843            2 :           l_proj=.FALSE.
     844            4 :           DO j=jspin,0,-1
     845            4 :              INQUIRE(file=TRIM('proj'//spin012(j)),exist=l_proj)
     846            4 :              IF(l_proj)THEN
     847            2 :                 filename='proj'//spin012(j)
     848            2 :                 EXIT
     849              :              ENDIF
     850              :           ENDDO
     851              : 
     852            2 :           IF(l_proj)THEN
     853            2 :              OPEN (203,file=TRIM(filename),status='old')
     854            2 :              REWIND (203)
     855            2 :              READ (203,*) nwfs,numbands
     856            2 :              REWIND (203)
     857            2 :              CLOSE (203)
     858              :           ELSEIF(wann%l_projmethod.OR.wann%l_bestproj&
     859            0 :                .OR.wann%l_matrixamn)THEN
     860              :              CALL juDFT_error("no proj/proj.1/proj.2"&
     861            0 :                   ,calledby ="wannier")
     862              :           ENDIF
     863              : 
     864              : 
     865            2 :           jspin2=jspin
     866            2 :           IF(noco%l_soc .AND. input%jspins.EQ.1)jspin2=1
     867            2 :           jspin2_b=jspin_b
     868            2 :           IF(noco%l_soc .AND. input%jspins.EQ.1)jspin2_b=1
     869              : 
     870            2 :           jsp_start = jspin ; jsp_end = jspin
     871              : 
     872              :           !ccccccccccc   read in the eigenvalues and vectors   cccccc
     873            2 :           WRITE(*,*)'wannierspin',wannierspin
     874            4 :           DO jspin5=1,wannierspin!1!2
     875              :              !       jspin5=jspin
     876            2 :              jsp_start=jspin5; jsp_end=jspin5
     877            2 :              nrec5=0
     878              :              IF(.NOT.noco%l_noco) nrec5 = (jspin5-1)*nkpts
     879              : 
     880              :              CALL cdn_read0(eig_id,fmpi%irank,fmpi%isize,jspin5,input%jspins, &!wannierspin instead of DIMENSION%jspd?&
     881            4 :                   noco%l_noco, n_bands,n_size)
     882              : 
     883              :           ENDDO
     884              :           !..   now we want to define the maximum number of the bands by all kpts
     885            2 :           nbnd = 0
     886              : 
     887            2 :           i_rec = 0 ; n_rank = 0
     888              :           !*************************************************************
     889              :           !..writing down the eig.1 and/or eig.2 files
     890              : 
     891              :           !..write individual files if multi-spiral mode wann%l_ms=T
     892              :           !*************************************************************
     893            2 :           IF(l_p0)THEN
     894              :              CALL wann_write_eig(&
     895              :                   fmpi,cell,noco,nococonv,input,kpts,sym,atoms,  &
     896              :                   eig_id,l_real,&
     897              :                   atoms%ntype,&
     898              :                   lapw%dim_nvd(),wannierspin,&
     899              :                   fmpi%isize,jspin,&
     900              :                   noco%l_ss,noco%l_noco,nrec,fullnkpts,&
     901              :                   wann%l_bzsym,wann%l_byindex,wann%l_bynumber,&
     902              :                   wann%l_byenergy,&
     903              :                   irreduc ,wann%band_min(jspin),&
     904              :                   wann%band_max(jspin),&
     905              :                   numbands,&
     906              :                   sliceplot%e1s,sliceplot%e2s,efermi,.FALSE.,nkpts,&
     907            1 :                   nbnd,kpoints,l_gwf,iqpt)
     908              : 
     909              :           ENDIF!l_p0
     910              : 
     911              :           ! nbnd is calculated for process zero and is sent here to the others
     912              : #ifdef CPP_MPI
     913            2 :           IF(l_p0)THEN
     914            2 :              DO cpu_index=1,fmpi%isize-1
     915            2 :                 CALL MPI_SEND(nbnd,1,MPI_INTEGER,cpu_index,1,fmpi%mpi_comm,ierr(1))
     916              :              ENDDO
     917              :           ELSE
     918            1 :              CALL MPI_RECV(nbnd,1,MPI_INTEGER,0,1,fmpi%mpi_comm,stt,ierr(1))
     919              :           ENDIF
     920              : #endif
     921              : 
     922            2 :           PRINT*,"process: ",fmpi%irank," nbnd= ",nbnd
     923              :           !##################################################################
     924            2 :           IF(wann%l_mmn0)THEN
     925            6 :              ALLOCATE ( mmn(nbnd,nbnd,fullnkpts) )
     926         1170 :              mmn(:,:,:) = CMPLX(0.,0.)
     927            2 :              IF((noco%l_soc.OR.noco%l_noco) .AND. (doublespin.EQ.1))&
     928            0 :                   ALLOCATE(socmmn(nbnd,nbnd,fullnkpts) )
     929              :           ENDIF
     930            2 :           IF(wann%l_nabla)THEN
     931            0 :              ALLOCATE ( nablamat(3,nbnd,nbnd,fullnkpts) )
     932            0 :              nablamat = CMPLX(0.,0.)
     933              :           ENDIF
     934              : 
     935            2 :           IF(wann%l_soctomom)THEN
     936            0 :              ALLOCATE ( soctomom(3,nbnd,nbnd,fullnkpts) )
     937            0 :              soctomom = CMPLX(0.,0.)
     938              :           ENDIF
     939              : 
     940            2 :           IF(wann%l_surfcurr)THEN
     941            0 :              ALLOCATE ( surfcurr(3,nbnd,nbnd,fullnkpts) )
     942            0 :              surfcurr = CMPLX(0.,0.)
     943              :           ENDIF
     944              : 
     945            2 :           IF(wann%l_anglmom)THEN
     946            0 :              IF(.NOT.ALLOCATED(anglmom))THEN
     947            0 :                 ALLOCATE ( anglmom(3,nbnd,nbnd,fullnkpts) )
     948            0 :                 anglmom=CMPLX(0.,0.)
     949              :              ENDIF
     950              :           ENDIF
     951              : 
     952            2 :           IF(wann%l_orbcomp)THEN
     953            0 :              IF(ALLOCATED(orbcomp))DEALLOCATE(orbcomp)
     954            0 :              IF(wann%l_oc_f)THEN
     955            0 :                 ALLOCATE(orbcomp(16,wann%oc_num_orbs,nbnd,nbnd,fullnkpts))
     956              :              ELSE
     957            0 :                 ALLOCATE(orbcomp(9,wann%oc_num_orbs,nbnd,nbnd,fullnkpts))
     958              :              ENDIF
     959            0 :              orbcomp=CMPLX(0.,0.)
     960              :           ENDIF
     961              : 
     962              :           !write (*,*) 'nwfs=',nwfs
     963            2 :           IF(wann%l_projmethod.OR.wann%l_bestproj.OR.wann%l_matrixamn)THEN
     964            2 :              IF(.NOT.ALLOCATED(amn))THEN
     965            6 :                 ALLOCATE ( amn(nbnd,nwfs,fullnkpts) )
     966         1170 :                 amn(:,:,:) = CMPLX(0.,0.)
     967              :              ENDIF
     968              :           ENDIF
     969              : 
     970            2 :           IF (wann%l_projmethod.OR.wann%l_bestproj) THEN
     971            0 :              ALLOCATE ( psiw(nbnd,nwfs,fullnkpts) )
     972            0 :              psiw(:,:,:) = CMPLX(0.,0.)
     973            0 :              IF(.NOT.ALLOCATED(hwfr))THEN
     974            0 :                 ALLOCATE ( hwfr(nwfs,nwfs) )
     975            0 :                 hwfr(:,:) = CMPLX(0.,0.)
     976              :              ENDIF
     977              :           ENDIF
     978              : 
     979              : 
     980            2 :           IF (wann%l_matrixmmn) THEN
     981            2 :              IF(.NOT.ALLOCATED(mmnk))THEN
     982            8 :                 ALLOCATE ( mmnk(nbnd,nbnd,nntot,fullnkpts) )
     983         9362 :                 mmnk = (0.,0.)
     984              :              ENDIF
     985              :           ENDIF
     986              : 
     987            2 :           IF(wann%l_matrixmmn)THEN
     988            2 :              IF(.NOT.ALLOCATED(mmnk_q).AND.l_gwf)THEN
     989            0 :                 ALLOCATE ( mmnk_q (nbnd,nbnd,nntot_q,fullnkpts) )
     990            0 :                 mmnk_q = (0.,0.)
     991              : 
     992              :                 !             allocate ( m_int(nbnd,nbnd,nntot_q,fullnkpts) )
     993              :                 !             allocate ( m_sph(nbnd,nbnd,nntot_q,fullnkpts) )
     994              :                 !             allocate ( m_vac(nbnd,nbnd,nntot_q,fullnkpts) )
     995              :                 !             m_int = cmplx(0.,0.)
     996              :                 !             m_sph = cmplx(0.,0.)
     997              :                 !             m_vac = cmplx(0.,0.)
     998              :              ENDIF
     999              :           ENDIF
    1000              : 
    1001              : 
    1002            8 :           ALLOCATE ( flo(atoms%ntype,atoms%jmtd,2,atoms%nlod,2) )
    1003              : 
    1004            4 :           DO jspin4=1,wannierspin!2
    1005            2 :              jspin3=jspin4
    1006            2 :              IF(input%jspins.EQ.1) jspin3=1
    1007            2 :              na = 1
    1008            6 :              DO  n = 1,atoms%ntype
    1009           16 :                 DO  l = 0,atoms%lmax(n)
    1010              :                    !...compute the l-dependent, k-independent radial MT- basis functions
    1011              : 
    1012              :                    CALL radfun(&
    1013              :                         l,n,jspin4,enpara%el0(l,n,jspin3),vr(1,n,jspin3),atoms,&
    1014              :                         ff(n,:,:,l,jspin4),gg(n,:,:,l,jspin4),usdus,&
    1015           16 :                         nodeu,noded,wronk)
    1016              : 
    1017              :                 ENDDO
    1018              :                 !...and the local orbital radial functions
    1019            4 :                 DO ilo = 1, atoms%nlo(n)
    1020              : 
    1021              :                    CALL radflo(&
    1022              :                         atoms,n,jspin4,enpara%ello0(:,:,jspin3),vr(1,n,jspin3),&
    1023              :                         ff(n,1:,1:,0:,jspin4),gg(n,1:,1:,0:,jspin4),fmpi,&
    1024            2 :                         usdus,uuilon,duilon,ulouilopn,flo(n,:,:,:,jspin4))
    1025              : 
    1026              :                 ENDDO
    1027              :                 !       na = na + atoms%neq(n)
    1028              :              ENDDO
    1029              :           ENDDO!jspin3
    1030              :           !****************************************************************
    1031              :           !   calculate the k-independent uju*gaunt-matrix needed for
    1032              :           !   mmnmatrix
    1033              :           !****************************************************************
    1034              :           ! TODO: make this more efficient (i.e., compute ujugaunt only once
    1035              :           ! and not for all q-points).
    1036            2 :           IF(wann%l_matrixmmn)THEN
    1037            0 :              ALLOCATE(ujug(0:lmd,0:lmd,&
    1038            8 :                   1:atoms%ntype,1:nntot))
    1039            0 :              ALLOCATE(ujdg(0:lmd,0:lmd,&
    1040            8 :                   1:atoms%ntype,1:nntot))
    1041            0 :              ALLOCATE(djug(0:lmd,0:lmd,&
    1042            6 :                   1:atoms%ntype,1:nntot))
    1043            0 :              ALLOCATE(djdg(0:lmd,0:lmd,&
    1044            6 :                   1:atoms%ntype,1:nntot))
    1045            0 :              ALLOCATE(ujulog(0:lmd,1:atoms%nlod,-atoms%llod:atoms%llod,&
    1046           10 :                   1:atoms%ntype,1:nntot))
    1047            0 :              ALLOCATE(djulog(0:lmd,1:atoms%nlod,-atoms%llod:atoms%llod,&
    1048            8 :                   1:atoms%ntype,1:nntot))
    1049            0 :              ALLOCATE(ulojug(0:lmd,1:atoms%nlod,-atoms%llod:atoms%llod,&
    1050            8 :                   1:atoms%ntype,1:nntot))
    1051            0 :              ALLOCATE(ulojdg(0:lmd,1:atoms%nlod,-atoms%llod:atoms%llod,&
    1052            8 :                   1:atoms%ntype,1:nntot))
    1053            0 :              ALLOCATE(ulojulog(1:atoms%nlod,-atoms%llod:atoms%llod,&
    1054              :                   1:atoms%nlod,-atoms%llod:atoms%llod,&
    1055           12 :                   1:atoms%ntype,1:nntot))
    1056              : 
    1057              :              CALL wann_ujugaunt(&
    1058              :                   atoms%llod,nntot,kdiff,atoms%lmax,atoms%ntype,&
    1059              :                   atoms%ntype,cell%bbmat,cell%bmat,atoms%nlod,atoms%nlo,&
    1060              :                   atoms%llo,flo(:,:,:,:,jspin),&
    1061              :                   flo(:,:,:,:,jspin),&
    1062              :                   ff(:,:,:,:,jspin),&
    1063              :                   ff(:,:,:,:,jspin),&
    1064              :                   gg(:,:,:,:,jspin),&
    1065              :                   gg(:,:,:,:,jspin),atoms%jri,atoms%rmsh,atoms%dx,&
    1066              :                   atoms%jmtd,atoms%lmaxd,lmd,&
    1067              :                   ujug,ujdg,djug,djdg,&
    1068            2 :                   ujulog,djulog,ulojug,ulojdg,ulojulog,.FALSE.,1)
    1069              : 
    1070              :              ! compute integrals of radial solution, according energy derivatives,
    1071              :              ! the spherical Bessel function and the Gaunt coefficients in order
    1072              :              ! to account for the overlap of the lattice periodic parts at
    1073              :              ! neighboring q-points
    1074            2 :              IF(l_gwf)THEN
    1075            0 :                 ALLOCATE(ujug_q(0:lmd,0:lmd,&
    1076            0 :                      1:atoms%ntype,1:nntot_q))
    1077            0 :                 ALLOCATE(ujdg_q(0:lmd,0:lmd,&
    1078            0 :                      1:atoms%ntype,1:nntot_q))
    1079            0 :                 ALLOCATE(djug_q(0:lmd,0:lmd,&
    1080            0 :                      1:atoms%ntype,1:nntot_q))
    1081            0 :                 ALLOCATE(djdg_q(0:lmd,0:lmd,&
    1082            0 :                      1:atoms%ntype,1:nntot_q))
    1083            0 :                 ALLOCATE(ujulog_q(0:lmd,1:atoms%nlod,-atoms%llod:atoms%llod,&
    1084            0 :                      1:atoms%ntype,1:nntot_q))
    1085            0 :                 ALLOCATE(djulog_q(0:lmd,1:atoms%nlod,-atoms%llod:atoms%llod,&
    1086            0 :                      1:atoms%ntype,1:nntot_q))
    1087            0 :                 ALLOCATE(ulojug_q(0:lmd,1:atoms%nlod,-atoms%llod:atoms%llod,&
    1088            0 :                      1:atoms%ntype,1:nntot_q))
    1089            0 :                 ALLOCATE(ulojdg_q(0:lmd,1:atoms%nlod,-atoms%llod:atoms%llod,&
    1090            0 :                      1:atoms%ntype,1:nntot_q))
    1091            0 :                 ALLOCATE(ulojulog_q(1:atoms%nlod,-atoms%llod:atoms%llod,&
    1092              :                      1:atoms%nlod,-atoms%llod:atoms%llod,&
    1093            0 :                      1:atoms%ntype,1:nntot_q))
    1094              : 
    1095              :                 ! we need G(q+b)/2 as argument for the sph. Bessel func.
    1096              :                 ! and additionally a spin-dependent sign (-/+ 1)^{lpp}
    1097            0 :                 IF(wann%l_sgwf) CALL wann_ujugaunt(&
    1098              :                      atoms%llod,nntot_q,qdiff/2.0,atoms%lmax,atoms%ntype,&
    1099              :                      atoms%ntype,cell%bbmat,cell%bmat,atoms%nlod,atoms%nlo,&
    1100              :                      atoms%llo,flo(:,:,:,:,jspin),&
    1101              :                      flo(:,:,:,:,jspin_b),&
    1102              :                      ff(:,:,:,:,jspin),&
    1103              :                      ff(:,:,:,:,jspin_b),&
    1104              :                      gg(:,:,:,:,jspin),&
    1105              :                      gg(:,:,:,:,jspin_b),atoms%jri,atoms%rmsh,atoms%dx,&
    1106              :                      atoms%jmtd,atoms%lmaxd,lmd,&
    1107              :                      ujug_q,ujdg_q,djug_q,djdg_q,&
    1108              :                      ujulog_q,djulog_q,ulojug_q,ulojdg_q,ulojulog_q,.TRUE.,&
    1109            0 :                      sign_q)
    1110              : 
    1111            0 :                 IF(wann%l_socgwf) CALL wann_ujugaunt(&
    1112              :                      atoms%llod,nntot_q,zero_qdiff,atoms%lmax,atoms%ntype,&
    1113              :                      atoms%ntype,cell%bbmat,cell%bmat,atoms%nlod,atoms%nlo,&
    1114              :                      atoms%llo,flo(:,:,:,:,jspin),&
    1115              :                      flo(:,:,:,:,jspin_b),&
    1116              :                      ff(:,:,:,:,jspin),&
    1117              :                      ff(:,:,:,:,jspin_b),&
    1118              :                      gg(:,:,:,:,jspin),&
    1119              :                      gg(:,:,:,:,jspin_b),atoms%jri,atoms%rmsh,atoms%dx,&
    1120              :                      atoms%jmtd,&
    1121              :                      atoms%lmaxd,lmd,ujug_q,ujdg_q,djug_q,djdg_q,&
    1122              :                      ujulog_q,djulog_q,ulojug_q,ulojdg_q,ulojulog_q,&
    1123            0 :                      .FALSE.,1)
    1124              : 
    1125              :              ENDIF ! l_gwf
    1126              : 
    1127              :           ENDIF !l_matrixmmn
    1128              : ! Replace the following by calls to zmat%init further below
    1129              : !          zzMat%l_real = l_real
    1130              : !          zzMat%matsize1 = lapw%dim_nbasfcn()
    1131              : !          zzMat%matsize2 = input%neig
    1132              : !          IF(l_real) THEN
    1133              : !             IF(.NOT.ALLOCATED(zzMat%data_r))&
    1134              : !               ALLOCATE (zzMat%data_r(zzMat%matsize1,zzMat%matsize2))
    1135              : !          ELSE
    1136              : !               IF(.NOT.ALLOCATED(zzMat%data_c))&
    1137              : !               ALLOCATE (zzMat%data_c(zzMat%matsize1,zzMat%matsize2))
    1138              : !          END IF
    1139              : 
    1140              : !          zMat%l_real = zzMat%l_real
    1141              : !          zMat%matsize1 = zzMat%matsize1
    1142              : !          zMat%matsize2 = zzMat%matsize2
    1143              : !          IF (zzMat%l_real) THEN
    1144              : !             IF(.NOT.ALLOCATED(zMat%data_r))&
    1145              : !                  ALLOCATE (zMat%data_r(zMat%matsize1,zMat%matsize2))
    1146              : !             zMat%data_r = 0.0
    1147              : !          ELSE
    1148              : !             IF(.NOT.ALLOCATED(zMat%data_c))&
    1149              : !                  ALLOCATE (zMat%data_c(zMat%matsize1,zMat%matsize2))
    1150              : !             zMat%data_c = CMPLX(0.0,0.0)
    1151              : !          END IF
    1152              : 
    1153              :  !         zMat_b%l_real = zzMat%l_real
    1154              :  !         zMat_b%matsize1 = zzMat%matsize1
    1155              :  !         zMat_b%matsize2 = zzMat%matsize2
    1156              :  !         IF (zzMat%l_real) THEN
    1157              :  !              IF(.NOT.ALLOCATED(zMat_b%data_r))&
    1158              :  !              ALLOCATE (zMat_b%data_r(zMat_b%matsize1,zMat_b%matsize2))
    1159              :  !            zMat_b%data_r = 0.0
    1160              :  !         ELSE
    1161              :  !              IF(.NOT.ALLOCATED(zMat_b%data_c))&
    1162              :  !              ALLOCATE (zMat_b%data_c(zMat_b%matsize1,zMat_b%matsize2))
    1163              :  !            zMat_b%data_c = CMPLX(0.0,0.0)
    1164              :  !         END IF
    1165              : 
    1166            2 :           i_rec = 0 ; n_rank = 0
    1167              : 
    1168              :           !****************************************************************
    1169              :           !.. loop by kpoints starts!      each may be a separate task
    1170              :           !****************************************************************
    1171           18 :           DO ikpt = wann%ikptstart,fullnkpts  ! loop by k-points starts
    1172           16 :              kptibz=ikpt
    1173           16 :              IF(wann%l_bzsym) kptibz=irreduc(ikpt)
    1174           16 :              IF(wann%l_bzsym) oper=mapkoper(ikpt)
    1175              : 
    1176           16 :              i_rec = i_rec + 1
    1177           18 :              IF (MOD(i_rec-1,fmpi%isize).EQ.fmpi%irank) THEN
    1178              : 
    1179            8 :                 ALLOCATE ( we(input%neig),eigg(input%neig) )
    1180              : 
    1181            8 :                 n_start=1
    1182            8 :                 n_end=input%neig
    1183              : 
    1184              : 
    1185              :                 ! read information of diagonalization for fixed q-point iqpt
    1186              :                 ! stored in the eig file on unit 66. the lattice respectively
    1187              :                 ! plane-wave vectors G(k,q) are saved in (k1,k2,k3).
    1188              : 
    1189            8 :                 CALL lapw%init(input,noco,nococonv,kpts,atoms,sym,kptibz,cell,fmpi)
    1190              : 
    1191            8 :           nbasfcn = MERGE(lapw%nv(1)+lapw%nv(2)+2*atoms%nlotot,lapw%nv(1)+atoms%nlotot,noco%l_noco)
    1192            8 :                   CALL zzMat%init(l_real,nbasfcn,input%neig)
    1193            8 :                       CALL zMat%init(l_real,nbasfcn,input%neig)
    1194              : 
    1195              : 
    1196              :                 CALL cdn_read(&
    1197              :                      eig_id,&
    1198              :                      lapw%dim_nvd(),input%jspins,fmpi%irank,fmpi%isize, &!wannierspin instead of DIMENSION%jspd?&
    1199              :                      kptibz,jspin,lapw%dim_nbasfcn(),&
    1200              :                      noco%l_ss,noco%l_noco,input%neig,n_start,n_end,&
    1201            8 :                      nbands,eigg,zzMat)
    1202              : 
    1203              : 
    1204            8 :                 nslibd = 0
    1205              : 
    1206              :                 !...we work only within the energy window
    1207              : 
    1208           72 :                 eig(:) = 0.
    1209              : 
    1210              :                 !      print*,"bands used:"
    1211              : 
    1212           72 :                 DO i = 1,nbands
    1213              :                    IF ((eigg(i).GE.sliceplot%e1s.AND.nslibd.LT.numbands.AND.&
    1214              :                         wann%l_bynumber).OR.&
    1215              :                         (eigg(i).GE.sliceplot%e1s.AND.eigg(i).LE.sliceplot%e2s.AND.&
    1216           64 :                         wann%l_byenergy).OR.(i.GE.wann%band_min(jspin2).AND.&
    1217            8 :                         (i.LE.wann%band_max(jspin2)).AND.wann%l_byindex))THEN
    1218              : 
    1219              :                       !           print*,i
    1220           64 :                       nslibd = nslibd + 1
    1221           64 :                       eig(nslibd) = eigg(i)
    1222           64 :                       we(nslibd) = we(i)
    1223           64 :                       IF(zzMat%l_real) THEN
    1224         6056 :                          zMat%data_r(:,nslibd) = zzMat%data_r(:,i)
    1225              :                       ELSE
    1226            0 :                          zMat%data_c(:,nslibd) = zzMat%data_c(:,i)
    1227              :                       END IF
    1228              :                    ENDIF
    1229              :                 ENDDO
    1230              : 
    1231              :                 !***********************************************************
    1232              :                 !              rotate the wavefunction
    1233              :                 !***********************************************************
    1234            8 :                 IF (wann%l_bzsym.AND.oper.NE.1) THEN  !rotate bkpt
    1235              :                     call wann_kptsrotate(&
    1236              :                            atoms, &
    1237              :                            sym%invsat,&
    1238              :                            noco%l_noco,noco%l_soc,&
    1239              :                            jspin,&
    1240              :                            oper,sym%nop,sym%mrot,lapw%dim_nvd(),&
    1241              :                            shiftkpt(:,ikpt),&
    1242              :                            sym%tau,&
    1243              :                            lapw,&
    1244            0 :                            zMat,nsfactor)
    1245              :                 ELSE
    1246            8 :                    nsfactor=CMPLX(1.0,0.0)
    1247              :                 ENDIF
    1248              : 
    1249              :                 !******************************************************************
    1250              : 
    1251              :                 !...the overlap matrix Mmn which is computed for each k- and b-point
    1252              : 
    1253            8 :                 noccbd = nslibd
    1254              : 
    1255            0 :                 ALLOCATE(acof(noccbd,0:lmd,atoms%nat),&
    1256            0 :                      bcof(noccbd,0:lmd,atoms%nat),&
    1257           56 :                      ccof(-atoms%llod:atoms%llod,noccbd,atoms%nlod,atoms%nat))
    1258              : 
    1259        14152 :                 acof(:,:,:) = CMPLX(0.,0.) ; bcof(:,:,:) = CMPLX(0.,0.)
    1260           24 :                 ccof(:,:,:,:) = CMPLX(0.,0.)
    1261              : 
    1262              :                 !...generation the A,B,C coefficients in the spheres
    1263              :                 !...for the lapws and local orbitals, summed by the basis functions
    1264              : 
    1265              : 
    1266              :                 CALL abcof(input,atoms,sym,cell,lapw,noccbd,usdus,&
    1267            8 :                      noco,nococonv,jspin2 ,acof,bcof,ccof,zMat)
    1268              : 
    1269              : 
    1270            8 :                 CALL wann_abinv(atoms,sym, acof,bcof,ccof)
    1271              : 
    1272              : 
    1273            8 :                 IF((doublespin.EQ.3).OR.(doublespin.EQ.4)) GOTO 9900
    1274              : 
    1275              : 
    1276            8 :                 IF(wann%l_orbcomp)THEN
    1277              :                    CALL wann_orbcomp(atoms,usdus,jspin,acof,bcof,ccof,&
    1278            0 :                         wann%oc_num_orbs,wann%oc_orbs,wann%l_oc_f,orbcomp(:,:,:,:,ikpt))
    1279              :                 ENDIF
    1280              : 
    1281            8 :                 IF(wann%l_anglmom)THEN
    1282            0 :                    CALL wann_anglmom(atoms,usdus,jspin, acof,bcof,ccof,anglmom(:,:,:,ikpt))
    1283              :                 ENDIF
    1284              : 
    1285              : #ifdef CPP_TOPO
    1286              :                 IF(wann%l_surfcurr)THEN
    1287              :                    !         call wann_surfcurr_int(
    1288              :                    !     >        lapw%dim_nv2d(),jspin,stars%ng3,vacuum%nmzxyd,stars%ng2,sphhar%ntypsd,
    1289              :                    !     >        atoms%ntype,atoms%lmaxd,atoms%jmtd,atoms%ntype,atoms%nat,vacuum%nmzd,atoms%neq,stars%ng3,vacuum%nvac,
    1290              :                    !     >        vacuum%nmz,vacuum%nmzxy,stars%ng2,sym%nop,sym%nop2,cell%volint,input%film,sliceplot%slice,sym%symor,
    1291              :                    !     >        sym%invs,sym%invs2,cell%z1,vacuum%delz,sym%ngopr,sym%ntypsy,atoms%jri,atoms%pos,atoms%zatom,
    1292              :                    !     >        atoms%lmax,sym%mrot,sym%tau,atoms%rmsh,sym%invtab,cell%amat,cell%bmat,cell%bbmat,ikpt,sliceplot%nnne,sliceplot%kk,
    1293              :                    !     >        lapw%dim_nvd(),atoms%nlod,atoms%llod,nv(jspin),lmd,lapw%bkpt,cell%omtil,atoms%nlo,atoms%llo,
    1294              :                    !     >        k1(:,jspin),k2(:,jspin),k3(:,jspin),evac(:,jspin),
    1295              :                    !     >        vz(:,:,jspin2),
    1296              :                    !     >        nslibd,lapw%dim_nbasfcn(),input%neig,ff,gg,flo,acof,bcof,ccof,z,
    1297              :                    !     >        surfcurr(:,:,:,ikpt))
    1298              : 
    1299              :                    CALL wann_surfcurr_int2(&
    1300              :                         lapw%dim_nv2d(),jspin ,stars%ng3,&
    1301              :                         vacuum%nmzxyd,&
    1302              :                         stars%ng2,sphhar%ntypsd,atoms%ntype,atoms%lmaxd,&
    1303              :                         atoms%jmtd,atoms%ntype,atoms%nat,vacuum%nmzd,&
    1304              :                         atoms%neq,stars%ng3,vacuum%nvac,vacuum%nmz,&
    1305              :                         vacuum%nmzxy,stars%ng2,sym%nop,sym%nop2,cell%volint,&
    1306              :                         input%film,sliceplot%slice,sym%symor,&
    1307              :                         sym%invs,sym%invs2,cell%z1,vacuum%delz,sym%ngopr,&
    1308              :                         sym%ntypsy,atoms%jri,atoms%pos,atoms%taual,&
    1309              :                         atoms%zatom,atoms%rmt,atoms%lmax,sym%mrot,sym%tau,&
    1310              :                         atoms%rmsh,sym%invtab,cell%amat,cell%bmat,cell%bbmat,&
    1311              :                         ikpt,lapw%dim_nvd(),lapw%nv(jspin),lapw%bkpt,cell%omtil,&
    1312              :                         lapw%k1(:,jspin),lapw%k2(:,jspin),lapw%k3(:,jspin),&
    1313              :                         nslibd,lapw%dim_nbasfcn(),input%neig,z,&
    1314              :                         dirfacs,&
    1315              :                         surfcurr(:,:,:,ikpt))
    1316              : 
    1317              :                    CALL wann_surfcurr(&
    1318              :                         dirfacs,cell%amat,&
    1319              :                         jspin,atoms%ntype,atoms%lmaxd,atoms%lmax,atoms%nat,&
    1320              :                         atoms%neq,noccbd,lmd,atoms%nat,atoms%llod,atoms%nlod,&
    1321              :                         atoms%nlo,atoms%llo, &
    1322              :                         acof,bcof,ccof,&
    1323              :                         us(:,:,jspin),dus(:,:,jspin),duds(:,:,jspin),&
    1324              :                         uds(:,:,jspin),&
    1325              :                         ulos(:,:,jspin),dulos(:,:,jspin),&
    1326              :                         atoms%rmt,atoms%pos, &
    1327              :                         surfcurr(:,:,:,ikpt))
    1328              :                    WRITE(oUnit,*)"dirfacs=",dirfacs
    1329              :                 ENDIF
    1330              : 
    1331              :                 IF(wann%l_soctomom)THEN
    1332              :                    CALL wann_soc_to_mom(&
    1333              :                         jspin,atoms%ntype,atoms%lmaxd,atoms%lmax,atoms%nat,&
    1334              :                         atoms%jmtd,atoms%jri,atoms%rmsh,atoms%dx,atoms%neq,&
    1335              :                         noccbd,lmd,atoms%nat,atoms%llod,atoms%nlod,&
    1336              :                         vso(:,:,1), &
    1337              :                         ff(:,:,:,:,jspin),gg(:,:,:,:,jspin),&
    1338              :                         acof,bcof,ccof,&
    1339              :                         soctomom(:,:,:,ikpt))
    1340              :                 ENDIF
    1341              : 
    1342              :                 IF(wann%l_nabla)THEN
    1343              :                    CALL wann_nabla(&
    1344              :                         atoms%nlo,atoms%llo,&
    1345              :                         jspin,atoms%ntype,atoms%lmaxd,atoms%lmax,atoms%nat,&
    1346              :                         atoms%jmtd,atoms%jri,atoms%rmsh,atoms%dx,atoms%neq,&
    1347              :                         noccbd,lmd,atoms%nat,atoms%llod,atoms%nlod, &
    1348              :                         ff(:,:,:,:,jspin),gg(:,:,:,:,jspin),flo(:,:,:,:,jspin),&
    1349              :                         acof,bcof,ccof,&
    1350              :                         nablamat(:,:,:,ikpt))
    1351              :                    IF(input%film)THEN
    1352              :                       CALL wann_nabla_vac(&
    1353              :                            cell%z1,vacuum%nmzd,lapw%dim_nv2d(),&
    1354              :                            stars%mx1,stars%mx2,stars%mx3,&
    1355              :                            stars%ng3,vacuum%nvac,stars%ig,vacuum%nmz,vacuum%delz,&
    1356              :                            stars%ig2,cell%area,cell%bmat,cell%bbmat,enpara%evac0(:,jspin),&
    1357              :                            lapw%bkpt,vz(:,:,jspin2),nslibd,jspin,lapw%k1,lapw%k2,lapw%k3,&
    1358              :                            wannierspin,lapw%dim_nvd(),lapw%dim_nbasfcn(),input%neig,z,nv,&
    1359              :                            cell%omtil,&
    1360              :                            nablamat(:,:,:,ikpt))
    1361              :                    ENDIF
    1362              :                    addnoco=0
    1363              :                    DO  i = n_rank+1,lapw%nv(jspin),n_size
    1364              :                       b1 = lapw%bkpt+lapw%gvec(:,i,jspin)
    1365              :                       b2(1)=b1(1)*cell%bmat(1,1)+b1(2)*cell%bmat(2,1)+&
    1366              :                            b1(3)*cell%bmat(3,1)
    1367              :                       b2(2)=b1(1)*cell%bmat(1,2)+b1(2)*cell%bmat(2,2)+&
    1368              :                            b1(3)*cell%bmat(3,2)
    1369              :                       b2(3)=b1(1)*cell%bmat(1,3)+b1(2)*cell%bmat(2,3)+&
    1370              :                            b1(3)*cell%bmat(3,3)
    1371              :                       DO  j = n_rank+1,lapw%nv(jspin),n_size
    1372              :                          !-->     determine index and phase factor
    1373              :                          i1 = lapw%k1(j,jspin) - lapw%k1(i,jspin)
    1374              :                          i2 = lapw%k2(j,jspin) - lapw%k2(i,jspin)
    1375              :                          i3 = lapw%k3(j,jspin) - lapw%k3(i,jspin)
    1376              :                          in = stars%ig(i1,i2,i3)
    1377              :                          IF (in.EQ.0) CYCLE
    1378              :                          phase   = stars%rgphs(i1,i2,i3)
    1379              :                          phasust = CMPLX(phase,0.0)*stars%ustep(in)
    1380              : 
    1381              :                          DO m = 1,nslibd
    1382              :                             DO n = 1,nslibd
    1383              :                                DO dir=1,3
    1384              : 
    1385              : #if ( !defined(CPP_INVERSION) || defined(CPP_SOC) )
    1386              :                                   VALUE=phasust*z(i+addnoco,m)*CONJG(z(j+addnoco,n))
    1387              :                                   nablamat(dir,m,n,ikpt) =&
    1388              :                                        nablamat(dir,m,n,ikpt) -&
    1389              :                                        VALUE*b2(dir)
    1390              : #else
    1391              :                                   VALUE=phasust*CMPLX(z(i+addnoco,m)*z(j+addnoco,n),0.0)
    1392              :                                   nablamat(dir,m,n,ikpt) =&
    1393              :                                        nablamat(dir,m,n,ikpt) - &
    1394              :                                        VALUE*b2(dir)
    1395              : #endif
    1396              :                                ENDDO
    1397              :                             ENDDO
    1398              :                          ENDDO
    1399              : 
    1400              :                       ENDDO
    1401              :                    ENDDO
    1402              : 
    1403              :                 ENDIF
    1404              : #endif
    1405              :                 !      goto jump no longer needed?
    1406              :                 !      if ((.not.wann%l_matrixmmn).and.(.not.wann%l_matrixamn).and.
    1407              :                 !          (.not.wann%l_bestproj).and.(.not.wann%l_projmethod).and.
    1408              :                 !          (.not.wann%l_mmn0)) goto 3
    1409              : 
    1410              : 
    1411              :                 !------mmn0-matrix
    1412            8 :                 IF(wann%l_mmn0)THEN
    1413            8 :                    addnoco=0
    1414            8 :                    noconbasfcn=nbasfcn
    1415            8 :                    IF(noco%l_noco.AND.(jspin.EQ.2))THEN
    1416            0 :                       addnoco=lapw%nv(1)+atoms%nlotot
    1417              : !                      noconbasfcn=nbasfcn-addnoco
    1418              :                    ENDIF
    1419              : 
    1420              :                    !-----> interstitial contribution to mmn0-matrix
    1421              : 
    1422              :                    CALL wann_mmkb_int(&
    1423              :                         cmplx_1,addnoco,addnoco,&
    1424              :                         lapw%dim_nvd(),stars%mx1,stars%mx2,stars%mx3,&
    1425              :                         stars%ng3,lapw%k1(:,jspin2),lapw%k2(:,jspin2),lapw%k3(:,jspin2),&
    1426              :                         lapw%nv(jspin2),input%neig,noconbasfcn,noconbasfcn,zMat,nslibd,&
    1427              :                         lapw%k1(:,jspin2),lapw%k2(:,jspin2),lapw%k3(:,jspin2),&
    1428              :                         lapw%nv(jspin2),zMat,nslibd,&
    1429              :                         nbnd,&
    1430              :                         stars%rgphs,stars%ustep,stars%ig,(/ 0,0,0 /),&
    1431            8 :                         mmn(:,:,ikpt))
    1432              : 
    1433              :                    !---> spherical contribution to mmn0-matrix
    1434              : 
    1435              :                    CALL wann_mmk0_sph(&
    1436              :                         atoms%llod,noccbd,atoms%nlod,atoms%nat,atoms%ntype,&
    1437              :                         atoms%lmaxd,atoms%lmax,lmd,atoms%ntype,atoms%neq,&
    1438              :                         atoms%nlo,atoms%llo,acof(1:noccbd,:,:),&
    1439              :                         bcof(1:noccbd,:,:),ccof(:,1:noccbd,:,:),&
    1440              :                         usdus%ddn(:,:,jspin),usdus%uulon(:,:,jspin),&
    1441              :                         usdus%dulon(:,:,jspin),usdus%uloulopn(:,:,:,jspin),&
    1442            8 :                         mmn(:,:,ikpt))
    1443              :                    !---> vacuum contribution to mmn0-matrix
    1444              : 
    1445            8 :                    IF (input%film ) THEN
    1446              : 
    1447              :                       CALL wann_mmk0_vac(&
    1448              :                            noco%l_noco,atoms%nlotot,qpt_i,&
    1449              :                            cell%z1,vacuum%nmzd,lapw%dim_nv2d(),&
    1450              :                            stars%mx1,stars%mx2,stars%mx3,&
    1451              :                            stars%ng3,vacuum%nvac,stars%ig,vacuum%nmz,vacuum%delz,&
    1452              :                            stars%ig2,cell%area,cell%bmat,&
    1453              :                            cell%bbmat,enpara%evac0(:,jspin2),lapw%bkpt,vz(:,:,jspin2),&
    1454              :                            nslibd,jspin2,lapw%k1,lapw%k2,lapw%k3,wannierspin,lapw%dim_nvd(),&
    1455              :                            lapw%dim_nbasfcn(),input%neig,zMat,lapw%nv,cell%omtil,&
    1456            0 :                            mmn(:,:,ikpt))
    1457              :                    ENDIF
    1458              :                 ENDIF !l_mmn0
    1459              : 
    1460              : 
    1461              :                 !---> overlaps with the trial orbitals
    1462            8 :                 IF (wann%l_projmethod.OR.wann%l_bestproj.OR.wann%l_matrixamn) THEN
    1463            8 :                    l_amn2=.FALSE.
    1464            8 :                    amnchi = cmplx_1!cmplx(1.,0.)
    1465              : 
    1466              :                    CALL wann_amn (&
    1467              :                         amnchi,nslibd,nwfs,atoms%ntype,atoms%nlod,atoms%llod,&
    1468              :                         atoms%llo,atoms%nlo,atoms%lmaxd,atoms%jmtd,lmd,&
    1469              :                         atoms%neq,atoms%nat,ikpt,nbnd,&
    1470              :                         atoms%rmsh,atoms%rmt,atoms%jri,atoms%dx,atoms%lmax,&
    1471              :                         usdus%us(:,:,jspin),usdus%dus(:,:,jspin),&
    1472              :                         usdus%uds(:,:,jspin),&
    1473              :                         usdus%duds(:,:,jspin),flo(:,:,:,:,jspin),&
    1474              :                         ff(:,:,:,:,jspin),gg(:,:,:,:,jspin),acof,bcof,ccof,&
    1475              :                         (noco%l_soc.OR.noco%l_noco),jspin,&
    1476            8 :                         l_amn2,amn(:,:,ikpt))
    1477            8 :                    IF(l_amn2)THEN
    1478              :                       CALL wann_amn (&
    1479              :                            amnchi,nslibd,nwfs,atoms%ntype,atoms%nlod,atoms%llod,&
    1480              :                            atoms%llo,atoms%nlo,atoms%lmaxd,atoms%jmtd,lmd,&
    1481              :                            atoms%neq,atoms%nat,ikpt,nbnd,&
    1482              :                            atoms%rmsh,atoms%rmt,atoms%jri,atoms%dx,atoms%lmax,&
    1483              :                            usdus%us(:,:,jspin),usdus%dus(:,:,jspin),&
    1484              :                            usdus%uds(:,:,jspin),&
    1485              :                            usdus%duds(:,:,jspin),flo(:,:,:,:,jspin),&
    1486              :                            ff(:,:,:,:,jspin),gg(:,:,:,:,jspin),acof,bcof,ccof,&
    1487              :                            (noco%l_soc.OR.noco%l_noco),jspin,&
    1488            0 :                            l_amn2,amn(:,:,ikpt),lapw%bkpt)
    1489              :                    ENDIF
    1490              :                    !         amn(ikpt,:,:)=amn(ikpt,:,:)*conjg(nsfactor)
    1491              :                 ENDIF
    1492              : 
    1493              : 
    1494              : 
    1495              :                 !****************************************************************
    1496              :                 !...         vanderbilt mmn matrix
    1497              :                 !***************************************************************
    1498            8 :                 IF (wann%l_matrixmmn .AND.&
    1499              :                      (.NOT.wann%l_skipkov)) THEN   !  vanderbilt procedure Mmn matrix
    1500            8 :                    ALLOCATE ( we_b(input%neig) )
    1501              : 
    1502              : !!! the cycle by the nearest neighbors (nntot) for each kpoint
    1503              : 
    1504           72 :                    DO   ikpt_b = 1,nntot
    1505              : 
    1506           64 :                       IF(pair_to_do(ikpt,ikpt_b).EQ.0)CYCLE !save time by symmetry
    1507           32 :                       kptibz_b=bpt(ikpt_b,ikpt)
    1508           32 :                       IF(wann%l_bzsym) oper_b=mapkoper(kptibz_b)
    1509           32 :                       IF (wann%l_bzsym) kptibz_b=irreduc(kptibz_b)
    1510              : 
    1511              : 
    1512              : 
    1513              :                       ! now we need the wavefunctions for k_b kpoint
    1514              : 
    1515              :                       !        print*,"something to do"
    1516              : 
    1517              : 
    1518           32 :                       n_start=1
    1519           32 :                       n_end=input%neig
    1520           32 :                       call lapw_b%init(input,noco,nococonv,kpts,atoms,sym,kptibz_b,cell,fmpi)
    1521              : 
    1522              : 
    1523           32 :                     nbasfcn_b = MERGE(lapw_b%nv(1)+lapw_b%nv(2)+2*atoms%nlotot,lapw_b%nv(1)+atoms%nlotot,noco%l_noco)
    1524           32 :                               CALL zMat_b%init(l_real,nbasfcn_b,input%neig)
    1525           32 :                               CALL zzMat%init(l_real,nbasfcn_b,input%neig)
    1526              : 
    1527              : 
    1528              :                       CALL cdn_read(&
    1529              :                            eig_id,&
    1530              :                            lapw%dim_nvd(),input%jspins,fmpi%irank,fmpi%isize, &!wannierspin instead of DIMENSION%jspd?&
    1531              :                            kptibz_b,jspin,nbasfcn_b,&
    1532              :                            noco%l_ss,noco%l_noco,input%neig,n_start,n_end,&
    1533           32 :                            nbands_b,eigg,zzMat)
    1534              : 
    1535              : 
    1536           32 :                       nslibd_b = 0
    1537              : 
    1538          288 :                       eig_b(:) = 0.
    1539              : 
    1540          288 :                       DO i = 1,nbands_b
    1541              :                          IF((eigg(i).GE.sliceplot%e1s.AND.nslibd_b.LT.numbands&
    1542              :                               .AND.wann%l_bynumber).OR.&
    1543              :                               (eigg(i).GE.sliceplot%e1s.AND.eigg(i).LE.sliceplot%e2s.AND.&
    1544          256 :                               wann%l_byenergy).OR.(i.GE.wann%band_min(jspin2).AND.&
    1545              :                               (i.LE.wann%band_max(jspin2)).AND.&
    1546           32 :                               wann%l_byindex))THEN
    1547          256 :                             nslibd_b = nslibd_b + 1
    1548          256 :                             eig_b(nslibd_b) = eigg(i)
    1549          256 :                             we_b(nslibd_b) = we_b(i)
    1550          256 :                             IF (zzMat%l_real) THEN
    1551        23552 :                                zMat_b%data_r(:,nslibd_b) = zzMat%data_r(:,i)
    1552              :                             ELSE
    1553            0 :                                zMat_b%data_c(:,nslibd_b) = zzMat%data_c(:,i)
    1554              :                             END IF
    1555              :                          ENDIF
    1556              :                       ENDDO
    1557              : 
    1558              :                       !***********************************************************
    1559              :                       !              Rotate the wavefunction of next neighbor.
    1560              :                       !***********************************************************
    1561           32 :                       IF (wann%l_bzsym .AND. (oper_b.NE.1)  ) THEN
    1562              :                      call wann_kptsrotate(&
    1563              :                            atoms, &
    1564              :                            sym%invsat,&
    1565              :                            noco%l_noco,noco%l_soc,&
    1566              :                            jspin,&
    1567              :                            oper_b,sym%nop,sym%mrot,lapw_b%dim_nvd(),&
    1568              :                            shiftkpt(:,bpt(ikpt_b,ikpt)),&
    1569              :                            sym%tau,&
    1570              :                            lapw_b,&
    1571            0 :                            zMat_b,nsfactor_b)
    1572              :                       ELSE
    1573           32 :                          nsfactor_b=CMPLX(1.0,0.0)
    1574              :                       ENDIF
    1575              :                       !      print*,"kpt2=",bkpt_b
    1576              : 
    1577           32 :                       noccbd_b = nslibd_b
    1578              : 
    1579              :                       !cccc   we start with the Mmn matrix   ccccccccccccc
    1580              : 
    1581              : !!! matrix elements of the interstitial overlap
    1582              : !!! matrix with the weights given by products of
    1583              : !!! the c-coeff. for different G-vectors, bands and k-points
    1584              : !!! Mmn(k,b)(IR) = \sum(G,G')C_G^(k,n)*C_G'^(k+b,m)\theta_(G-G')
    1585              : 
    1586              :                       !ccccccccccc  Spherical Contributions         ccccccccccccc
    1587              : 
    1588            0 :                       ALLOCATE (acof_b(noccbd_b,0:lmd,atoms%nat),&
    1589            0 :                            bcof_b(noccbd_b,0:lmd,atoms%nat),&
    1590            0 :                            ccof_b(-atoms%llod:atoms%llod,noccbd_b,atoms%nlod,&
    1591          224 :                            atoms%nat))
    1592              : 
    1593              : !!! get the band-dependent k-dependent ab coeff.
    1594              : 
    1595              :                       CALL abcof(input,atoms,sym,cell,lapw_b,&
    1596              :                            noccbd_b,usdus,noco,nococonv,jspin2 ,&
    1597           32 :                            acof_b,bcof_b,ccof_b,zMat_b)
    1598              : 
    1599              : 
    1600              :                       CALL wann_abinv(atoms,sym,&
    1601           32 :                            acof_b,bcof_b,ccof_b)
    1602              : 
    1603              :                       !cccccccc  Interstitial  ccccccccccccccccccccccccc
    1604              : !!! matrix elements of the interstitial overlap
    1605              : !!! matrix with the weights given by products of
    1606              : !!! the c-coeff. for different G-vectors, bands and k-points
    1607              : !!! Mmn(k,b)(IR) = \sum(G,G')C_G^(k,n)*C_G'^(k+b,m)\theta_(G-G')
    1608              :                       !... these overlaps are the same as in the mmk0 case, only differ
    1609              :                       !... by the b-dependence of the C-coefficients
    1610              :                       !cccccccccccccccccccccccccccccccccccccccccccccccccccc
    1611              : 
    1612           32 :                       addnoco=0
    1613           32 :                       addnoco2=0
    1614           32 :                       noconbasfcn=nbasfcn
    1615           32 :                       noconbasfcn_b=nbasfcn_b
    1616           32 :                       IF(noco%l_noco.AND.(jspin.EQ.2))THEN
    1617            0 :                          addnoco  = lapw%nv(1)   + atoms%nlotot
    1618            0 :                          addnoco2 = lapw_b%nv(1) + atoms%nlotot
    1619              : !                         noconbasfcn=nbasfcn-addnoco
    1620              : !                         noconbasfcn_b=nbasfcn_b-addnoco2
    1621              :                       ENDIF
    1622              : 
    1623              :                       CALL wann_mmkb_int(&
    1624              :                            cmplx_1,addnoco,addnoco2,&
    1625              :                            lapw%dim_nvd(),stars%mx1,stars%mx2,stars%mx3,&
    1626              :                            stars%ng3,lapw%k1(:,jspin2),lapw%k2(:,jspin2),lapw%k3(:,jspin2),&
    1627              :                            lapw%nv(jspin2),input%neig,noconbasfcn,noconbasfcn_b,zMat,nslibd,&
    1628              :                            lapw_b%k1(:,jspin2),lapw_b%k2(:,jspin2),lapw_b%k3(:,jspin2),&
    1629              :                            lapw_b%nv(jspin2),zMat_b,nslibd_b,&
    1630              :                            nbnd,&
    1631              :                            stars%rgphs,stars%ustep,stars%ig,gb(:,ikpt_b,ikpt),&
    1632           32 :                            mmnk(:,:,ikpt_b,ikpt))
    1633              : 
    1634              : 
    1635              : 
    1636              :                       !ccccccccccc  Spherical Contributions  ccccccccccccc
    1637              : 
    1638           64 :                       chi = cmplx_1
    1639              :                       CALL wann_mmkb_sph(&
    1640              :                            nbnd,atoms%llod,nslibd,nslibd_b,atoms%nlod,atoms%nat,&
    1641              :                            atoms%ntype,lmd,atoms%jmtd,atoms%taual,sym%nop,atoms%lmax,&
    1642              :                            atoms%ntype,atoms%neq,atoms%nlo,atoms%llo,acof,bcof,ccof,&
    1643              :                            lapw_b%bkpt,acof_b,bcof_b,ccof_b,gb(:,ikpt_b,ikpt),lapw%bkpt,&
    1644              :                            ujug,ujdg,&
    1645              :                            djug,djdg,ujulog,djulog,ulojug,ulojdg,ulojulog,kdiff,&
    1646              :                            nntot,chi,&
    1647           32 :                            mmnk(:,:,ikpt_b,ikpt))
    1648              : 
    1649              : 
    1650              : 
    1651              : 
    1652              :                       !...vacuum contributions
    1653           32 :                       IF (input%film ) THEN
    1654              : 
    1655              :                          CALL wann_mmkb_vac(&
    1656              :                               cmplx_1,noco%l_noco,atoms%nlotot,qpt_i,&
    1657              :                               nbnd,cell%z1,vacuum%nmzd,lapw%dim_nv2d(),&
    1658              :                               stars%mx1,stars%mx2,stars%mx3,&
    1659              :                               stars%ng3,vacuum%nvac,stars%ig,vacuum%nmz,&
    1660              :                               vacuum%delz,stars%ig2,cell%area,cell%bmat,&
    1661              :                               cell%bbmat,enpara%evac0(:,jspin2),&
    1662              :                               enpara%evac0(:,jspin2_b),&
    1663              :                               lapw%bkpt,lapw_b%bkpt,vz(:,:,jspin2),vz(:,:,jspin2_b),&
    1664              :                               nslibd,nslibd_b,jspin2,jspin2_b,&
    1665              :                               lapw%k1,lapw%k2,lapw%k3,lapw_b%k1,lapw_b%k2,lapw_b%k3,&
    1666              :                               wannierspin,lapw%dim_nvd(),&
    1667              :                               lapw%dim_nbasfcn(),input%neig,zMat,zMat_b,&
    1668              :                               lapw%nv,lapw_b%nv,cell%omtil,&
    1669              :                               gb(:,ikpt_b,ikpt),&
    1670            0 :                               mmnk(:,:,ikpt_b,ikpt))
    1671              :                       
    1672              :                       ENDIF
    1673              : 
    1674              :                       !        mmnk(:,:,ikpt_b,ikpt)=
    1675              :                       !                mmnk(:,:,ikpt_b,ikpt)*nsfactor*conjg(nsfactor_b)
    1676              : 
    1677              : 
    1678           72 :                       DEALLOCATE ( acof_b,bcof_b,ccof_b )
    1679              : 
    1680              : 
    1681              :                    enddo ! end of loop by the nearest k-neighbors
    1682              : 
    1683            8 :                    DEALLOCATE ( we_b )
    1684              : 
    1685              :                 ENDIF!l_matrixmmn=.true.
    1686              : 
    1687              : 9900            CONTINUE ! jump for doublespin loop
    1688              : 
    1689            8 :                 IF (wann%l_matrixmmn) THEN   !  vanderbilt procedure Mmn matrix
    1690              : 
    1691              :                    !*******************************************c
    1692              :                    !          START Q-NEIGHBOR LOOP            c
    1693              :                    !*******************************************c
    1694            8 :                    ALLOCATE ( we_qb(input%neig) )
    1695              : 
    1696            8 :                    DO iqpt_b=1,nntot_q
    1697            8 :                       IF(.NOT.l_gwf) EXIT              ! old functionality
    1698              : 
    1699            0 :                       qptibz_b = bpt_q(iqpt_b,iqpt)
    1700            0 :                       IF(qptibz_b.EQ.qptibz) CYCLE     ! no need to compute overlaps
    1701              :                       ! with periodic images for now
    1702              : 
    1703            0 :                       qptb_i = nococonv%qss
    1704            0 :                       alphb_i = nococonv%alph
    1705            0 :                       betab_i = nococonv%beta
    1706            0 :                       thetab_i = nococonv%theta
    1707            0 :                       phib_i = nococonv%phi
    1708            0 :                       IF(wann%l_sgwf) THEN
    1709            0 :                          qptb_i(:) = wann%param_vec(:,qptibz_b)
    1710            0 :                          alphb_i(:) = wann%param_alpha(:,qptibz_b)
    1711            0 :                       ELSEIF(wann%l_socgwf) THEN
    1712            0 :                          IF(wann%l_dim(2)) phib_i = tpi_const*wann%param_vec(2,qptibz_b)
    1713            0 :                          IF(wann%l_dim(3)) thetab_i = tpi_const*wann%param_vec(3,qptibz_b)
    1714              :                       ENDIF
    1715              : 
    1716              :                       !if(pair_to_do_q(iqpt,iqpt_b).eq.0)cycle    ! TODO: use symmetry
    1717            0 :                       IF(wann%l_bzsym) oper_qb=mapqoper(qptibz_b)
    1718            0 :                       IF (wann%l_bzsym) qptibz_b=irreduc_q(qptibz_b)
    1719              : 
    1720            0 :                       n_start=1
    1721            0 :                       n_end=input%neig
    1722              : 
    1723              :                       ! read in diagonalization information from corresponding
    1724              :                       ! eig file to q-point iqpt_b at a given k-point ikpt.
    1725              :                       ! as a check verify that bkpt.eq.bqpt (same k).
    1726              :                       ! moreover, the plane-wave vectors G(k,q+b) are stored
    1727              :                       ! in (k1_qb,k2_qb,k3_qb) for later use.
    1728              : 
    1729            0 :                       CALL lapw_qb%init(input,noco,nococonv,kpts,atoms,sym,kptibz,cell,fmpi)
    1730              :                       CALL cdn_read(&
    1731              :                            innerEig_idList(iqpt_b),&
    1732              :                            lapw%dim_nvd(),input%jspins,fmpi%irank,fmpi%isize, &!wannierspin instead of DIMENSION%jspd? !kptibz_b2?&
    1733              :                            kptibz,jspin_b,lapw%dim_nbasfcn(),&
    1734              :                            noco%l_ss,noco%l_noco,input%neig,n_start,&
    1735              :                            n_end,&
    1736              :                            nbands_qb,eigg,   &
    1737            0 :                            zzMat)
    1738              : 
    1739              : 
    1740              :                       ! are we dealing with the same k-point at which
    1741              :                       ! we want to construct A,B,C coefficients etc. ?
    1742            0 :                       IF(ANY(bqpt.NE.lapw%bkpt)) CALL juDFT_error("bqpt.ne.bkpt",&
    1743            0 :                            calledby="wannier")
    1744              : 
    1745            0 :                       zMat_qb%l_real = zzMat%l_real
    1746            0 :                       zMat_qb%matsize1 = zzMat%matsize1
    1747            0 :                       zMat_qb%matsize2 = zzMat%matsize2
    1748            0 :                       IF (zzMat%l_real) THEN
    1749            0 :                          ALLOCATE (zMat_qb%data_r(zMat%matsize1,zMat%matsize2))
    1750            0 :                          zMat_qb%data_r = 0.0
    1751              :                       ELSE
    1752            0 :                          ALLOCATE (zMat_qb%data_c(zMat%matsize1,zMat%matsize2))
    1753            0 :                          zMat_qb%data_c = CMPLX(0.0,0.0)
    1754              :                       END IF
    1755              : 
    1756            0 :                       eig_qb(:) = 0.
    1757              : 
    1758            0 :                       nslibd_qb = 0
    1759            0 :                       DO i = 1,nbands_qb
    1760              :                          IF((eigg(i).GE.sliceplot%e1s.AND.nslibd_qb.LT.numbands&
    1761              :                               .AND.wann%l_bynumber).OR.&
    1762              :                               (eigg(i).GE.sliceplot%e1s.AND.eigg(i).LE.sliceplot%e2s.AND.&
    1763            0 :                               wann%l_byenergy).OR.(i.GE.wann%band_min(jspin).AND.&
    1764            0 :                               (i.LE.wann%band_max(jspin)).AND.wann%l_byindex)) THEN
    1765            0 :                             nslibd_qb = nslibd_qb + 1
    1766            0 :                             eig_qb(nslibd_qb) = eigg(i)
    1767            0 :                             we_qb(nslibd_qb) = we_qb(i)
    1768            0 :                             IF (zzMat%l_real) THEN
    1769            0 :                                zMat_qb%data_r(:,nslibd_qb) = zzMat%data_r(:,i)
    1770              :                             ELSE
    1771            0 :                                zMat_qb%data_c(:,nslibd_qb) = zzMat%data_c(:,i)
    1772              :                             END IF
    1773              :                          ENDIF
    1774              :                       ENDDO
    1775              : 
    1776              :                       ! check that eigenvectors and -values are identical if q=q+b
    1777            0 :                       IF(iqpt.EQ.qptibz_b .AND. jspin.EQ.jspin_b) THEN
    1778            0 :                          IF(zMat%l_real) THEN
    1779            0 :                             IF(ANY(zMat%data_r.NE.zMat_qb%data_r)) &
    1780            0 :                                  WRITE(*,*)'z.ne.z_qb',iqpt,ikpt
    1781              :                          ELSE
    1782            0 :                             IF(ANY(zMat%data_c.NE.zMat_qb%data_c)) &
    1783            0 :                                  WRITE(*,*)'z.ne.z_qb',iqpt,ikpt
    1784              :                          END IF
    1785            0 :                          IF(ANY(eig.NE.eig_qb)) WRITE(*,*)'eig.ne.eiq_qb',iqpt,ikpt
    1786            0 :                          IF(lapw%nv(jspin).NE.lapw_qb%nv(jspin)) WRITE(*,*)'nv!=nv_qb',iqpt,ikpt
    1787              :                       ENDIF
    1788              : 
    1789              :                       ! check that number of bands are the same at (k,q) and (k,q+b)
    1790            0 :                       IF(nslibd.NE.nslibd_qb)&
    1791            0 :                            WRITE(*,*)'nslibd.ne.nslibd_qb',ikpt,iqpt,iqpt_b
    1792              : 
    1793              : 
    1794            0 :                       noccbd_qb = nslibd_qb
    1795            0 :                       nsfactor_b=CMPLX(1.0,0.0)
    1796              : 
    1797              : 
    1798            0 :                       ALLOCATE (acof_qb(noccbd_qb,0:lmd,atoms%nat),&
    1799            0 :                            bcof_qb(noccbd_qb,0:lmd,atoms%nat),&
    1800            0 :                            ccof_qb(-atoms%llod:atoms%llod,noccbd_qb,atoms%nlod,&
    1801            0 :                            atoms%nat))
    1802              : 
    1803            0 :                       acof_qb(:,:,:) = CMPLX(0.,0.) ; bcof_qb(:,:,:) = CMPLX(0.,0.)
    1804            0 :                       ccof_qb(:,:,:,:) = CMPLX(0.,0.)
    1805              : 
    1806              :                       ! construct the A,B,C coefficients of the wave function
    1807              :                       ! at the point (k,q+b) using previously read information
    1808              : 
    1809              :                       CALL abcof(input,atoms,sym,cell,lapw_qb,&
    1810              :                            noccbd_qb,usdus,noco,nococonv,jspin_b ,&
    1811            0 :                            acof_qb,bcof_qb,ccof_qb,zMat_qb)
    1812              : 
    1813              :                       CALL wann_abinv(atoms,sym,&
    1814            0 :                            acof_qb,bcof_qb,ccof_qb)
    1815              : 
    1816              :                       ! check that A,B,C coefficients are the same if q+b = q
    1817            0 :                       IF(l_gwf.AND.(iqpt.EQ.qptibz_b).AND.(jspin.EQ.jspin_b)) THEN
    1818            0 :                          IF(ANY(acof_qb.NE.acof)) WRITE(*,*)'acof',iqpt,ikpt
    1819            0 :                          IF(ANY(bcof_qb.NE.bcof)) WRITE(*,*)'bcof',iqpt,ikpt
    1820            0 :                          IF(ANY(ccof_qb.NE.ccof)) WRITE(*,*)'ccof',iqpt,ikpt
    1821              :                       ENDIF
    1822              : 
    1823            0 :                       addnoco=0
    1824            0 :                       addnoco2=0
    1825            0 :                       IF(noco%l_noco.AND.(jspin.EQ.2))THEN
    1826            0 :                          addnoco  = lapw%nv(1)   + atoms%nlotot
    1827              :                       ENDIF
    1828            0 :                       IF(noco%l_noco.AND.(jspin_b.EQ.2))THEN
    1829            0 :                          addnoco2 = lapw_qb%nv(1) + atoms%nlotot
    1830              :                       ENDIF
    1831              : 
    1832              : 
    1833              :                       ! set up local->global transformation for overlaps
    1834            0 :                       DO n=1,atoms%ntype
    1835            0 :                          IF(wann%l_sgwf) THEN
    1836            0 :                             dalph = alph_i(n)-alphb_i(n)
    1837            0 :                             db1 = beta_i(n)/2.
    1838            0 :                             db2 = betab_i(n)/2.
    1839            0 :                          ELSEIF(wann%l_socgwf) THEN
    1840            0 :                             dalph = phi_i-phib_i
    1841            0 :                             db1 = theta_i/2.
    1842            0 :                             db2 = thetab_i/2.
    1843              :                          ENDIF
    1844              :                          coph = COS(dalph)
    1845              :                          siph = SIN(dalph)
    1846              :                          phasfac = CMPLX(coph,siph)
    1847              :                          phasfac2= CMPLX(coph,-siph)
    1848              : 
    1849            0 :                          IF(l_p0 .AND. dalph.NE.0.0) THEN
    1850            0 :                             WRITE(*,*)'WARNING: include dalph in chi trafo!'
    1851              :                          ENDIF
    1852              : 
    1853            0 :                          IF( (jspin.EQ.1) .AND. (jspin_b.EQ.1) ) THEN ! uu
    1854              :                             !              chi(n) = cos(db1)*cos(db2)*phasfac
    1855              :                             !     >               + sin(db1)*sin(db2)*phasfac2
    1856            0 :                             chi(n) = COS(db2-db1)
    1857            0 :                          ELSEIF( (jspin.EQ.2) .AND. (jspin_b.EQ.2) ) THEN !dd
    1858              :                             !              chi(n) = cos(db1)*cos(db2)*phasfac2
    1859              :                             !     >               + sin(db1)*sin(db2)*phasfac
    1860            0 :                             chi(n) = COS(db2-db1)
    1861            0 :                          ELSEIF( (jspin.EQ.1) .AND. (jspin_b.EQ.2) ) THEN ! ud
    1862              :                             !              chi(n) = sin(db1)*cos(db2)*phasfac2
    1863              :                             !     >               - cos(db1)*sin(db2)*phasfac
    1864            0 :                             chi(n) = -SIN(db2-db1)
    1865            0 :                          ELSEIF( (jspin.EQ.2) .AND. (jspin_b.EQ.1) ) THEN ! du
    1866              :                             !              chi(n) = cos(db1)*sin(db2)*phasfac2
    1867              :                             !     >               - sin(db1)*cos(db2)*phasfac
    1868            0 :                             chi(n) = SIN(db2-db1)
    1869              :                          ELSE
    1870            0 :                             STOP 'problem setting up chi: jspin,jspin_b'
    1871              :                          ENDIF
    1872              :                       ENDDO
    1873            0 :                       chi = CONJG(chi)
    1874              :                       !chi = cmplx_1
    1875              : 
    1876              :                       ! optional: disable chi transformation
    1877              :                       ! instead of computing overlap w.r.t. global frame
    1878              :                       ! only consider wave function overlaps in local frames
    1879            0 :                       IF(l_nochi) THEN
    1880            0 :                          IF(doublespin.LT.3) THEN
    1881            0 :                             chi = CMPLX(1.0,0.0)
    1882              :                          ELSE
    1883            0 :                             chi = CMPLX(0.0,0.0)
    1884              :                          ENDIF
    1885              :                       ENDIF
    1886              : 
    1887            0 :                       IF((iqpt.EQ.1).AND.(ikpt.EQ.1).AND.(iqpt_b.EQ.1)) THEN
    1888            0 :                          WRITE(*,*)'dbs',doublespin,'chi',chi(1)
    1889              :                       ENDIF
    1890              : 
    1891              :                       ! muffin tin contribution to overlap is computed taking into account
    1892              :                       ! the spin-dependent phase exp(+/- tau*b/2) and the ujugaunt integrals
    1893              :                       ! calculated especially for the q-points before. then, q and q+b
    1894              :                       ! take the role of k and k+b and the same for G(q+b) and G(k+b)
    1895            0 :                       IF(wann%l_sgwf) CALL wann_mmkb_sph( &
    1896              :                            nbnd,atoms%llod,nslibd,nslibd_qb,atoms%nlod,atoms%nat,&
    1897              :                            atoms%ntype,lmd,atoms%jmtd,sign_q*atoms%taual/2.0,sym%nop,&
    1898              :                            atoms%lmax,atoms%ntype,atoms%neq,atoms%nlo,atoms%llo,&
    1899              :                            acof,bcof,ccof,qptb_i,&
    1900              :                            acof_qb,bcof_qb,ccof_qb,gb_q(:,iqpt_b,iqpt),qpt_i,&
    1901              :                            ujug_q,ujdg_q,&
    1902              :                            djug_q,djdg_q,ujulog_q,djulog_q,ulojug_q,ulojdg_q,&
    1903              :                            ulojulog_q,qdiff,       &
    1904              :                            nntot_q,chi,&
    1905            0 :                            mmnk_q(:,:,iqpt_b,ikpt))
    1906            0 :                       IF(wann%l_socgwf) CALL wann_mmkb_sph( &
    1907              :                            nbnd,atoms%llod,nslibd,nslibd_qb,atoms%nlod,atoms%nat,&
    1908              :                            atoms%ntype,lmd,atoms%jmtd,&
    1909              :                            zero_taual,sym%nop,atoms%lmax,         &
    1910              :                            atoms%ntype,atoms%neq,atoms%nlo,atoms%llo,acof,bcof,ccof,&
    1911              :                            (/ 0.0, phib_i/tpi_const, thetab_i/tpi_const /),&
    1912              :                            acof_qb,bcof_qb,ccof_qb,gb_q(:,iqpt_b,iqpt),&
    1913              :                            (/ 0.0, phi_i/tpi_const, theta_i/tpi_const /),&
    1914              :                            ujug_q,ujdg_q,&
    1915              :                            djug_q,djdg_q,ujulog_q,djulog_q,ulojug_q,ulojdg_q,&
    1916              :                            ulojulog_q,qdiff,       &
    1917              :                            nntot_q,chi,&
    1918            0 :                            mmnk_q(:,:,iqpt_b,ikpt))
    1919              : 
    1920              : 
    1921              :                       IF(((doublespin.NE.3).AND.(doublespin.NE.4))&
    1922            0 :                            .OR.(.NOT.noco%l_noco)) THEN
    1923              : 
    1924            0 :                          IF(.NOT.noco%l_noco)THEN
    1925            0 :                             interchi=chi(1)
    1926            0 :                             vacchi=chi(1)
    1927              :                          ELSE
    1928            0 :                             interchi=cmplx_1
    1929            0 :                             vacchi=cmplx_1
    1930              :                          ENDIF
    1931              : 
    1932              :                          ! interstitial contribution to overlap is computed using
    1933              :                          ! (-/+ 1)*G(q+b)/2 as G-vector connecting neighbors
    1934              :                          ! and lattice vectors G(k,q) (k1...) and G(k,q+b) (k1_qb...)
    1935            0 :                          IF(wann%l_sgwf) CALL wann_mmkb_int(&
    1936              :                               interchi,addnoco,addnoco2,&
    1937              :                               lapw%dim_nvd(),stars%mx1,stars%mx2,stars%mx3,&
    1938              :                               stars%ng3,lapw%k1(:,jspin),lapw%k2(:,jspin),lapw%k3(:,jspin),&
    1939              :                               lapw%nv(jspin),input%neig,nbasfcn,nbasfcn_b,zMat,nslibd,&
    1940              :                               lapw_qb%k1(:,jspin_b),lapw_qb%k2(:,jspin_b),lapw_qb%k3(:,jspin_b),&
    1941              :                               lapw_qb%nv(jspin_b),zMat_qb,nslibd_qb,&
    1942              :                               nbnd,&
    1943              :                               stars%rgphs,stars%ustep,stars%ig,&
    1944              :                               sign_q*gb_q(:,iqpt_b,iqpt)/2,   &
    1945            0 :                               mmnk_q(:,:,iqpt_b,ikpt))
    1946            0 :                          IF(wann%l_socgwf) CALL wann_mmkb_int(&
    1947              :                               interchi,addnoco,addnoco2,&
    1948              :                               lapw%dim_nvd(),stars%mx1,stars%mx2,stars%mx3,&
    1949              :                               stars%ng3,lapw%k1(:,jspin),lapw%k2(:,jspin),lapw%k3(:,jspin),&
    1950              :                               lapw%nv(jspin),input%neig,nbasfcn,nbasfcn_b,zMat,nslibd,&
    1951              :                               lapw_qb%k1(:,jspin_b),lapw_qb%k2(:,jspin_b),lapw_qb%k3(:,jspin_b),&
    1952              :                               lapw_qb%nv(jspin_b),zMat_qb,nslibd_qb,&
    1953              :                               nbnd,&
    1954              :                               stars%rgphs,stars%ustep,stars%ig,(/ 0, 0, 0 /),  &
    1955            0 :                               mmnk_q(:,:,iqpt_b,ikpt))!m_int(:,:,iqpt_b,ikpt))
    1956              : 
    1957              : 
    1958              :                          ! vacuum contribution in film calculation
    1959            0 :                          IF (input%film ) THEN
    1960            0 :                             IF(wann%l_sgwf) CALL wann_mmkb_vac(&
    1961              :                                  vacchi,noco%l_noco,atoms%nlotot,sign_q*2.*lapw%bkpt,&
    1962              :                                  nbnd,cell%z1,vacuum%nmzd,lapw%dim_nv2d(),&
    1963              :                                  stars%mx1,stars%mx2,stars%mx3,&
    1964              :                                  stars%ng3,vacuum%nvac,stars%ig,vacuum%nmz,&
    1965              :                                  vacuum%delz,stars%ig2,cell%area,cell%bmat,&
    1966              :                                  cell%bbmat,enpara%evac0(:,jspin),enpara%evac0(:,jspin_b),&
    1967              :                                  sign_q*qpt_i/2.,&
    1968              :                                  sign_q*qptb_i/2.,&
    1969              :                                  vz(:,:,jspin2),vz(:,:,jspin2_b),&
    1970              :                                  nslibd,nslibd_qb,jspin,jspin_b,&
    1971              :                                  lapw%k1,lapw%k2,lapw%k3,lapw_qb%k1,lapw_qb%k2,lapw_qb%k3,&
    1972              :                                  wannierspin,lapw%dim_nvd(),&
    1973              :                                  lapw%dim_nbasfcn(),input%neig,zMat,zMat_qb,lapw%nv,&
    1974              :                                  lapw_qb%nv,cell%omtil,&
    1975              :                                  sign_q*gb_q(:,iqpt_b,iqpt)/2,&
    1976            0 :                                  mmnk_q(:,:,iqpt_b,ikpt))
    1977            0 :                             IF(wann%l_socgwf) CALL wann_mmkb_vac(&
    1978              :                                  vacchi,noco%l_noco,atoms%nlotot,qpt_i,&
    1979              :                                  nbnd,cell%z1,vacuum%nmzd,lapw%dim_nv2d(),&
    1980              :                                  stars%mx1,stars%mx2,stars%mx3,&
    1981              :                                  stars%ng3,vacuum%nvac,stars%ig,vacuum%nmz,&
    1982              :                                  vacuum%delz,stars%ig2,cell%area,cell%bmat,&
    1983              :                                  cell%bbmat,enpara%evac0(:,jspin),enpara%evac0(:,jspin_b),&
    1984              :                                  bqpt,bqpt,&
    1985              :                                  vz(:,:,jspin2),vz(:,:,jspin2_b),&
    1986              :                                  nslibd,nslibd_qb,jspin,jspin_b,&
    1987              :                                  lapw%k1,lapw%k2,lapw%k3,lapw_qb%k1,lapw_qb%k2,lapw_qb%k3,&
    1988              :                                  wannierspin,lapw%dim_nvd(),&
    1989              :                                  lapw%dim_nbasfcn(),input%neig,zMat,zMat_qb,lapw%nv,&
    1990              :                                  lapw_qb%nv,cell%omtil,&
    1991              :                                  (/ 0, 0, 0 /),&
    1992            0 :                                  mmnk_q(:,:,iqpt_b,ikpt))
    1993              : 
    1994              :                             ! vacuum contribution in one-dimensional chain calculation where
    1995              :                             ! q-point plays role of k-point with proper prefactor of (-/+ 1/2).
    1996              :                             ! moreover, the k-point ikpt is treated like qss in the subroutine
    1997              :                             ! such that a correction for sign and factor has to appear as well.
    1998              :                             ! lattice vectors G(k,q) (k1...) and G(k,q+b) (k1_qb...) need to
    1999              :                             ! be provided.
    2000              :                          ENDIF!film resp. odi
    2001              : 
    2002              : 
    2003              :                       ENDIF!doublespin
    2004              : 
    2005            8 :                       DEALLOCATE ( acof_qb,bcof_qb,ccof_qb )
    2006              : 
    2007              : 
    2008              :                    ENDDO !iqpt_b, q-neighbors
    2009              :                    !**************************************************c
    2010              :                    !              END Q-NEIGHBOR LOOP                 c
    2011              :                    !**************************************************c
    2012              : 
    2013            8 :                    DEALLOCATE ( we_qb )
    2014              : 
    2015              :                 ENDIF    !   if wann%l_matrixmmn = true
    2016              : 
    2017            8 :                 IF(.NOT.wann%l_bzsym)oper=0
    2018            8 :                 IF(.NOT.wann%l_plot_symm.OR.oper.EQ.1)THEN
    2019            8 :                    IF (wann%l_wann_plot .AND. &
    2020              :                         (doublespin.EQ.1 .OR. doublespin.EQ.2)) THEN
    2021              : 
    2022            0 :                       addnoco=0
    2023            0 :                       IF(noco%l_noco.AND.(jspin.EQ.2))THEN
    2024            0 :                          addnoco=lapw%nv(1)+atoms%nlotot
    2025              :                       ENDIF
    2026              : 
    2027            0 :                       IF (sliceplot%slice) THEN
    2028            0 :                          IF (ikpt.EQ.sliceplot%kk) THEN
    2029              : 
    2030            0 :                             WRITE (oUnit,*) 'nnne=',sliceplot%nnne
    2031            0 :                             WRITE (oUnit,*) 'eig(nnne)=',eig(sliceplot%nnne)
    2032            0 :                             WRITE (oUnit,*) 'we(nnne)=',we(sliceplot%nnne)
    2033              : 
    2034              :                             CALL wann_plot(&
    2035              :                                   vacuum,stars,cell,atoms,&
    2036              :                                  lapw%dim_nv2d(),jspin,stars%ng3,&
    2037              :                                  vacuum%nmzxyd,&
    2038              :                                  stars%ng2,sphhar%ntypsd,atoms%ntype,atoms%lmaxd,&
    2039              :                                  atoms%jmtd,atoms%ntype,atoms%nat,vacuum%nmzd,atoms%neq,&
    2040              :                                  stars%ng3,vacuum%nvac,vacuum%nmz,vacuum%nmzxy,stars%ng2,&
    2041              :                                  sym%nop,sym%nop2,cell%volint,input%film,sliceplot%slice,&
    2042              :                                  sym%symor,sym%invs,sym%invs2,cell%z1,vacuum%delz,&
    2043              :                                  sym%ngopr,sym%ntypsy,atoms%jri,atoms%pos,atoms%zatom,&
    2044              :                                  atoms%lmax,sym%mrot,sym%tau,atoms%rmsh,sym%invtab,&
    2045              :                                  cell%amat,cell%bmat,cell%bbmat,ikpt,sliceplot%nnne,&
    2046              :                                  sliceplot%kk,lapw%dim_nvd(),atoms%nlod,atoms%llod,&
    2047              :                                  lapw%nv(jspin),lmd,lapw%bkpt,cell%omtil,atoms%nlo,atoms%llo,&
    2048              :                                  lapw%k1(:,jspin),lapw%k2(:,jspin),lapw%k3(:,jspin),enpara%evac0(:,jspin),&
    2049              :                                  vz(:,:,jspin2),&
    2050              :                                  nslibd,lapw%dim_nbasfcn(),input%neig,&
    2051              :                                  ff(:,:,:,:,jspin),&
    2052              :                                  gg(:,:,:,:,jspin),flo,acof,bcof,ccof,zMat,&
    2053              :                                  stars%mx1,stars%mx2,stars%mx3,stars%ig,stars%ig2,&
    2054              :                                  stars%sk2,stars%phi2,&
    2055              :                                  noco%l_noco,noco%l_ss,qpt_i,&
    2056            0 :                                  addnoco,get_index_kq(ikpt,iqpt,fullnkpts),wann%l_sgwf)
    2057              : 
    2058              :                          ENDIF
    2059              :                       ELSE ! not sliceplot%slice
    2060              : 
    2061              :                          CALL wann_plot(&
    2062              :                                vacuum,stars,cell,atoms,&
    2063              :                               lapw%dim_nv2d(),jspin ,stars%ng3,&
    2064              :                               vacuum%nmzxyd,&
    2065              :                               stars%ng2,sphhar%ntypsd,atoms%ntype,atoms%lmaxd,&
    2066              :                               atoms%jmtd,atoms%ntype,atoms%nat,vacuum%nmzd,atoms%neq,&
    2067              :                               stars%ng3,vacuum%nvac,vacuum%nmz,vacuum%nmzxy,stars%ng2,&
    2068              :                               sym%nop,sym%nop2,cell%volint,input%film,sliceplot%slice,&
    2069              :                               sym%symor,sym%invs,sym%invs2,cell%z1,vacuum%delz,&
    2070              :                               sym%ngopr,sym%ntypsy,atoms%jri,atoms%pos,atoms%zatom,&
    2071              :                               atoms%lmax,sym%mrot,sym%tau,atoms%rmsh,sym%invtab,&
    2072              :                               cell%amat,cell%bmat,cell%bbmat,ikpt,sliceplot%nnne,&
    2073              :                               sliceplot%kk,lapw%dim_nvd(),atoms%nlod,atoms%llod,&
    2074              :                               lapw%nv(jspin),lmd,lapw%bkpt,cell%omtil,atoms%nlo,atoms%llo,&
    2075              :                               lapw%k1(:,jspin),lapw%k2(:,jspin),lapw%k3(:,jspin),enpara%evac0(:,jspin),&
    2076              :                               vz(:,:,jspin2),&
    2077              :                               nslibd,lapw%dim_nbasfcn(),input%neig,&
    2078              :                               ff(:,:,:,:,jspin),&
    2079              :                               gg(:,:,:,:,jspin),flo,acof,bcof,ccof,zMat,&
    2080              :                               stars%mx1,stars%mx2,stars%mx3,stars%ig,stars%ig2,&
    2081              :                               stars%sk2,stars%phi2,&
    2082              :                               noco%l_noco,noco%l_ss,qpt_i,&
    2083            0 :                               addnoco,get_index_kq(ikpt,iqpt,fullnkpts),wann%l_sgwf)
    2084              : 
    2085            0 :                          IF(wann%l_plot_symm.AND.wann%l_bzsym)THEN
    2086            0 :                             DO kplot=1,fullnkpts
    2087            0 :                                IF(irreduc(kplot).EQ.kptibz)THEN
    2088            0 :                                   plotoper=mapkoper(kplot)
    2089            0 :                                   IF(plotoper.LT.0)THEN
    2090            0 :                                      plotoper=-plotoper
    2091            0 :                                      l_conjugate=.TRUE.
    2092              :                                   ELSE
    2093            0 :                                      l_conjugate=.FALSE.
    2094              :                                   ENDIF
    2095            0 :                                   kplotoper=sym%invtab(plotoper)
    2096              :                                   CALL wann_plot_symm(jspin,sym%mrot(:,:,kplotoper),ikpt,&
    2097            0 :                                        kplot,l_conjugate)
    2098              :                                ENDIF
    2099              :                             ENDDO
    2100              :                          ENDIF
    2101              : 
    2102              :                       ENDIF
    2103              :                    ENDIF
    2104              :                 ENDIF !wann%l_plot_symm
    2105            8 :                 DEALLOCATE ( acof,bcof,ccof,we,eigg )
    2106              : 
    2107            8 :                 IF(wann%l_projmethod.OR.wann%l_bestproj)THEN
    2108              :                    CALL wann_projmethod(&
    2109              :                         fullnkpts,&
    2110              :                         wann%l_projmethod,wann%l_bestproj,&
    2111              :                         ikpt,nwfs,nslibd,amn,eig,&
    2112            0 :                         psiw,hwfr)
    2113              :                 ENDIF ! projmethod
    2114              : 
    2115              :              ENDIF   ! loop by processors
    2116              : 
    2117              :           ENDDO ! end of cycle by the k-points
    2118              : 
    2119              : 
    2120            2 :           IF(wann%l_matrixmmn)&
    2121            0 :                DEALLOCATE(ujug,ujdg,djug,djdg,&
    2122            2 :                ujulog,djulog,ulojulog,ulojug,ulojdg)
    2123            2 :           IF(wann%l_matrixmmn.AND.l_gwf)&
    2124            0 :                DEALLOCATE(ujug_q,ujdg_q,djug_q,&
    2125            0 :                djdg_q,ujulog_q,djulog_q,ulojulog_q,ulojug_q,&
    2126            0 :                ulojdg_q)
    2127              : 
    2128              : #ifdef CPP_MPI
    2129            2 :           CALL MPI_BARRIER(fmpi%mpi_comm,ierr(1))
    2130              : #endif
    2131              : 
    2132              : 
    2133              :           !******************************************************
    2134              :           !     Write down the projections.
    2135              :           !******************************************************
    2136              : 5         CONTINUE
    2137              : 
    2138            2 :           IF(doublespin.EQ.3 .OR. doublespin.EQ.4) GOTO 912
    2139              : 
    2140            2 :           IF(wann%l_nabla)THEN
    2141              : 
    2142            0 :              nablamat=nablamat*hescale
    2143              :              CALL wann_write_nabla(&
    2144              :                   fmpi%mpi_comm,l_p0,spin12(jspin)//'.nabl',&
    2145              :                   'Matrix elements of nabla operator',&
    2146              :                   nbnd,fullnkpts,nbnd,&
    2147              :                   fmpi%irank,fmpi%isize,wann%l_unformatted,&
    2148            0 :                   nablamat)
    2149              :           ENDIF
    2150              : 
    2151            2 :           IF(wann%l_soctomom)THEN
    2152            0 :              soctomom=soctomom*hescale
    2153              :              CALL wann_write_nabla(&
    2154              :                   fmpi%mpi_comm,l_p0,spin12(jspin)//'.stm',&
    2155              :                   'Matrix elements of stm operator',&
    2156              :                   nbnd,fullnkpts,nbnd,&
    2157              :                   fmpi%irank,fmpi%isize,wann%l_unformatted,&
    2158            0 :                   soctomom)
    2159              :           ENDIF
    2160              : 
    2161            2 :           IF(wann%l_surfcurr)THEN
    2162            0 :              surfcurr=surfcurr*hescale
    2163              :              CALL wann_write_nabla(&
    2164              :                   fmpi%mpi_comm,l_p0,spin12(jspin)//'.surfcurr',&
    2165              :                   'Surface currents',&
    2166              :                   nbnd,fullnkpts,nbnd,&
    2167              :                   fmpi%irank,fmpi%isize,wann%l_unformatted,&
    2168            0 :                   surfcurr)
    2169              :           ENDIF
    2170              : 
    2171            2 :           IF((noco%l_soc.OR.noco%l_noco).AND.wann%l_mmn0)THEN
    2172              :              CALL wann_write_amn(&
    2173              :                   fmpi%mpi_comm,&
    2174              :                   l_p0,spin12(jspin)//'.socmmn0',&
    2175              :                   'Overlaps of the wavefunct. at the same kpoint',&
    2176              :                   nbnd,fullnkpts,nbnd,&
    2177              :                   fmpi%irank,fmpi%isize,.FALSE.,.TRUE.,&
    2178            0 :                   mmn,wann%mmn0fmt==2)
    2179              :           ENDIF !noco%l_soc and l_mmn0
    2180              : 
    2181            2 :           IF(wann%l_orbcomp)THEN
    2182            0 :              num_angl=9
    2183            0 :              IF(wann%l_oc_f)num_angl=16
    2184              :              CALL wann_write_matrix5(&
    2185              :                   fmpi%mpi_comm,l_p0,spin12(jspin)//'.orbcomp',&
    2186              :                   'angular components',&
    2187              :                   nbnd,nbnd,&
    2188              :                   num_angl,wann%oc_num_orbs,fullnkpts,&
    2189              :                   fmpi%irank,fmpi%isize,wann%l_unformatted,&
    2190            0 :                   orbcomp)
    2191              :           ENDIF
    2192              : 
    2193            2 :           IF((noco%l_soc.OR.noco%l_noco) .AND. (doublespin.EQ.1)) THEN
    2194            0 :              IF(wann%l_mmn0) socmmn(:,:,:)=mmn(:,:,:)
    2195              :              GOTO 912
    2196              :           ENDIF
    2197              : 
    2198            2 :           IF(noco%l_soc.OR.noco%l_noco) THEN
    2199            0 :              jspin2=1
    2200            0 :              IF(wann%l_mmn0)      mmn(:,:,:)=socmmn(:,:,:)+mmn(:,:,:)
    2201            0 :              IF(wann%l_mmn0)      DEALLOCATE(socmmn)
    2202              :           ENDIF
    2203              : 
    2204            2 :           IF (wann%l_matrixamn)THEN
    2205              :              CALL wann_write_amn(&
    2206              :                   fmpi%mpi_comm,&
    2207              :                   l_p0,spin12(jspin2)//'.amn',&
    2208              :                   'Overlaps of the wavefunct. with the trial orbitals',&
    2209              :                   nbnd,fullnkpts,nwfs,&
    2210              :                   fmpi%irank,fmpi%isize,.FALSE.,.FALSE.,&
    2211            2 :                   amn(:,:,:),wann%matrixamnfmt==2)
    2212              :           ENDIF !wann%l_matrixamn
    2213              : 
    2214            2 :           IF(wann%l_anglmom)THEN
    2215              :              CALL wann_write_matrix4(&
    2216              :                   fmpi%mpi_comm,&
    2217              :                   l_p0,spin12(jspin2)//'.anglmom',&
    2218              :                   'Matrix elements of angular momentum',&
    2219              :                   nbnd,nbnd,3,fullnkpts,&
    2220              :                   fmpi%irank,fmpi%isize,&
    2221            0 :                   anglmom)
    2222              :           ENDIF
    2223              : 
    2224            2 :           IF (l_proj) THEN
    2225              :              !**************************************************************
    2226              :              !            for projmethod: write down WF1.umn
    2227              :              !*************************************************************
    2228            2 :              IF((wann%l_projmethod.OR.wann%l_bestproj))THEN
    2229              :                 CALL wann_write_amn(&
    2230              :                      fmpi%mpi_comm,l_p0,spin12(jspin2)//'.umn',&
    2231              :                      'transformation to first guess Wannier functions',&
    2232              :                      nbnd,fullnkpts,nwfs,&
    2233              :                      fmpi%irank,fmpi%isize,.FALSE.,.TRUE.,&
    2234            0 :                      psiw,.FALSE.)
    2235              : #ifdef CPP_MPI
    2236            0 :                 ALLOCATE( hwfr2(SIZE(hwfr,1),SIZE(hwfr,2)) )
    2237            0 :                 length=nwfs*nwfs
    2238              :                 CALL MPI_REDUCE(&
    2239              :                      hwfr,hwfr2,length,&
    2240              :                      MPI_DOUBLE_COMPLEX,MPI_SUM,0,&
    2241            0 :                      fmpi%mpi_comm,ierr(1))
    2242            0 :                 hwfr=hwfr2
    2243            0 :                 DEALLOCATE(hwfr2)
    2244              : #endif
    2245              :                 !********************************************************
    2246              :                 !        projmethod: hamiltonian matrix in real space
    2247              :                 !********************************************************
    2248            0 :                 IF(l_p0)THEN
    2249            0 :                    WRITE (oUnit,*) 'the hamiltonian matrix in real space:'
    2250            0 :                    DO i = 1,nwfs
    2251            0 :                       DO j = 1,nwfs
    2252            0 :                          WRITE (oUnit,*) '   WFs:',i,'and',j
    2253            0 :                          WRITE (oUnit,*) '     matrix element:',hwfr(i,j)
    2254              :                       ENDDO
    2255              :                    ENDDO
    2256              :                 ENDIF !l_p0
    2257            0 :                 DEALLOCATE(hwfr)
    2258              :              ENDIF !wann%l_projmethod or wann%l_bestproj
    2259              :           ENDIF !l_proj
    2260              : 
    2261              :           !*********************************************************
    2262              :           !.....write down the mmn0 matrix
    2263              :           !*********************************************************
    2264              : 
    2265            2 :           IF(wann%l_mmn0)THEN
    2266              :              CALL wann_write_amn(&
    2267              :                   fmpi%mpi_comm,&
    2268              :                   l_p0,spin12(jspin2)//'.mmn0',&
    2269              :                   'Overlaps of the wavefunct. at the same kpoint',&
    2270              :                   nbnd,fullnkpts,nbnd,&
    2271              :                   fmpi%irank,fmpi%isize,.FALSE.,.TRUE.,&
    2272            2 :                   mmn,wann%mmn0fmt==2)
    2273              :           ENDIF !wann%l_mmn0
    2274              : 
    2275              : 
    2276              :           !*****************************************************
    2277              :           !.....write down the matrix M^{k,b}_{mn}
    2278              :           !*****************************************************
    2279              : 
    2280              :           ! 912   continue
    2281            2 :           IF(wann%l_matrixmmn.AND.(.NOT.wann%l_skipkov))THEN
    2282              :              CALL wann_write_mmnk(&
    2283              :                   fmpi%mpi_comm,jspin2,l_p0,fullnkpts,nntot,wann,&
    2284              :                   maptopair,pair_to_do,nbnd,bpt,gb,&
    2285              :                   fmpi%isize,fmpi%irank,"            ",&
    2286            2 :                   mmnk,wann%matrixmmnfmt==2) !    wann%l_unformatted)
    2287              :           ENDIF !wann%l_matrixmmn
    2288              : 
    2289              : 912       CONTINUE
    2290              : 
    2291            2 :           IF(l_gwf .AND. wann%l_matrixmmn) THEN
    2292              :              !      mmnk_q = mmnk_q + m_sph+m_int+m_vac
    2293            0 :              IF(doublespin.EQ.doublespin_max) THEN
    2294            0 :                 WRITE(fname,'("param_",i4.4,".mmn")')iqpt
    2295              :                 CALL wann_write_mmnk2(l_p0,fullnkpts,nntot_q,wann,&
    2296              :                      nbnd,bpt_q(:,iqpt),gb_q(:,:,iqpt)/2,&
    2297              :                      fmpi%isize,fmpi%irank,fname,mmnk_q,&
    2298            0 :                      wann%l_unformatted)
    2299              :              ENDIF
    2300              : 
    2301              :              IF(.FALSE.) THEN
    2302              :                 WRITE(fname,'("param_",i4.4,"_",i1,".mmn")')iqpt,doublespin
    2303              :                 CALL wann_write_mmnk2(l_p0,fullnkpts,nntot_q,wann,&
    2304              :                      nbnd,bpt_q(:,iqpt),gb_q(:,:,iqpt)/2,&
    2305              :                      fmpi%isize,fmpi%irank,fname,&
    2306              :                      m_sph+m_int+m_vac,wann%l_unformatted)
    2307              : 
    2308              :                 WRITE(fname,'("param_",i4.4,"_",i1,"_int.mmn")')iqpt,doublespin
    2309              :                 CALL wann_write_mmnk2(l_p0,fullnkpts,nntot_q,wann,&
    2310              :                      nbnd,bpt_q(:,iqpt),gb_q(:,:,iqpt)/2,&
    2311              :                      fmpi%isize,fmpi%irank,fname,m_int,&
    2312              :                      wann%l_unformatted)
    2313              :                 WRITE(fname,'("param_",i4.4,"_",i1,"_sph.mmn")')iqpt,doublespin
    2314              :                 CALL wann_write_mmnk2(l_p0,fullnkpts,nntot_q,wann,&
    2315              :                      nbnd,bpt_q(:,iqpt),gb_q(:,:,iqpt)/2,&
    2316              :                      fmpi%isize,fmpi%irank,fname,m_sph,&
    2317              :                      wann%l_unformatted)
    2318              :                 WRITE(fname,'("param_",i4.4,"_",i1,"_vac.mmn")')iqpt,doublespin
    2319              :                 CALL wann_write_mmnk2(l_p0,fullnkpts,nntot_q,wann,&
    2320              :                      nbnd,bpt_q(:,iqpt),gb_q(:,:,iqpt)/2,&
    2321              :                      fmpi%isize,fmpi%irank,fname,m_vac,&
    2322              :                      wann%l_unformatted)
    2323              :              ENDIF
    2324              :              !      m_int = cmplx(0.,0.)
    2325              :              !      m_sph = cmplx(0.,0.)
    2326              :              !      m_vac = cmplx(0.,0.)
    2327              : 
    2328              :           ENDIF
    2329              : 
    2330              : 
    2331            2 :           IF( ALLOCATED (nablamat) ) DEALLOCATE( nablamat )
    2332            2 :           IF( ALLOCATED (soctomom) ) DEALLOCATE( soctomom )
    2333              : 
    2334            2 :           IF( ALLOCATED (surfcurr) ) DEALLOCATE( surfcurr )
    2335            2 :           IF( ALLOCATED( mmn ) ) DEALLOCATE(mmn)
    2336            2 :           IF( ALLOCATED( amn ) ) THEN
    2337            2 :              IF(.NOT.(noco%l_soc.OR.noco%l_noco))DEALLOCATE(amn)
    2338              :           ENDIF
    2339            2 :           IF ( ALLOCATED (psiw) ) DEALLOCATE ( psiw )
    2340            2 :           IF (wann%l_matrixmmn) THEN
    2341            2 :              IF(.NOT.(noco%l_soc.OR.noco%l_noco))DEALLOCATE (mmnk)
    2342              :           ENDIF
    2343            2 :           IF (wann%l_anglmom .AND. ALLOCATED(anglmom))THEN
    2344            0 :              IF(.NOT.(noco%l_soc.OR.noco%l_noco))DEALLOCATE (anglmom)
    2345              :           ENDIF
    2346            2 :           DEALLOCATE(flo)
    2347              : 
    2348            4 :           IF(.NOT.noco%l_noco)nrec=nrec+nkpts
    2349              : 
    2350              :           !#ifdef CPP_MPI
    2351              :           !      call MPI_BARRIER(fmpi%mpi_comm,ierr(1))
    2352              :           !#endif
    2353              : 
    2354              :        ENDDO ! end of cycle by spins
    2355              : 
    2356              : 
    2357              : #ifdef CPP_MPI
    2358            2 :        CALL MPI_BARRIER(fmpi%mpi_comm,ierr(1))
    2359              : #endif
    2360              : 
    2361              :        ! close eig files
    2362            2 :        IF (l_gwf) THEN
    2363              :           !         CALL close_eig(eig_id)
    2364            0 :           IF(wann%l_matrixmmn)THEN
    2365            0 :              DO iqpt_b=1,nntot_q
    2366              :                 !               CALL close_eig(innerEig_idList(iqpt_b))
    2367              :              ENDDO
    2368              :           ENDIF
    2369              :        ENDIF
    2370              : 
    2371            2 :        IF (wann%l_matrixmmn.AND.ALLOCATED(mmnk))THEN
    2372            0 :           DEALLOCATE ( mmnk )
    2373              :        ENDIF
    2374              : 
    2375            2 :        IF(ALLOCATED(mmnk_q)) DEALLOCATE(mmnk_q)
    2376              :        IF(ALLOCATED(m_int)) DEALLOCATE(m_int)
    2377              :        IF(ALLOCATED(m_sph)) DEALLOCATE(m_sph)
    2378              :        IF(ALLOCATED(m_vac)) DEALLOCATE(m_vac)
    2379              : 
    2380              :        IF ((wann%l_projmethod.OR.wann%l_bestproj.OR.wann%l_matrixamn)&
    2381            2 :             .AND.ALLOCATED(amn))THEN
    2382            0 :           DEALLOCATE ( amn )
    2383              :        ENDIF
    2384              : 
    2385            2 :        IF(wann%l_anglmom.AND.ALLOCATED(anglmom))THEN
    2386            0 :           DEALLOCATE ( anglmom )
    2387              :        ENDIF
    2388              : 
    2389              :        IF ((wann%l_projmethod.OR.wann%l_bestproj)&
    2390            2 :             .AND.ALLOCATED(hwfr)) THEN
    2391            0 :           DEALLOCATE ( hwfr )
    2392              :        ENDIF
    2393              : 
    2394            4 :        DEALLOCATE(innerEig_idList)
    2395              : 
    2396              :     ENDDO ! iqpt, q-points
    2397              :        !************************************************c
    2398              :        !               END Q LOOP                       c
    2399              :        !************************************************c
    2400              : 
    2401            2 :        IF(ALLOCATED(pair_to_do))DEALLOCATE(pair_to_do,maptopair)
    2402              : 
    2403            2 :        DEALLOCATE ( vr,vz)
    2404            2 :        DEALLOCATE ( ff,gg )
    2405            2 :        IF (wann%l_bzsym)  DEALLOCATE(irreduc,mapkoper)
    2406            2 :        IF (wann%l_bzsym.AND.l_gwf)  DEALLOCATE(irreduc_q,mapqoper)
    2407            2 :        IF(ALLOCATED(pair_to_do_q))&
    2408            0 :             DEALLOCATE(pair_to_do_q,maptopair_q)
    2409              : 
    2410            2 :        IF (ALLOCATED(kdiff)) DEALLOCATE ( kdiff )
    2411            2 :        IF (ALLOCATED(qdiff)) DEALLOCATE(qdiff,zero_qdiff)
    2412              : 
    2413              : 
    2414              : 9110   CONTINUE  ! jump for l_finishgwf=T
    2415              : 
    2416              :        ! correct for previously introduced factor of 2 in the
    2417              :        ! G-vectors connecting neighbors across the BZ boundary
    2418            2 :        IF(wann%l_sgwf) gb_q = gb_q/2
    2419            2 :        IF(wann%l_socgwf) gb_q = gb_q/2
    2420              : 
    2421              :        ! set up input files for wannier90 --> HDWFs
    2422              :        IF(l_p0.AND.l_gwf.AND.(wann%l_matrixmmn.OR.wann%l_matrixamn)&
    2423            2 :             .AND.(.NOT.wann%l_skipkov)) THEN
    2424              :           CALL wann_gwf_commat(fullnkpts,nntot,bpt,fullnqpts,&
    2425              :                nntot_q,bpt_q,gb,gb_q,wann%aux_latt_const,&
    2426              :                wann%l_unformatted,wann%l_matrixamn,&
    2427              :                wann%l_matrixmmn,wann%l_dim,&
    2428            0 :                wann%nparampts,wann%param_vec/2.0)
    2429              :        ENDIF
    2430              : 
    2431            2 :        IF(l_p0.AND.l_gwf.AND.wann%l_anglmom) THEN
    2432            0 :           CALL wann_gwf_anglmom(fullnkpts,fullnqpts,wann%l_unformatted)
    2433              :        ENDIF
    2434              : 
    2435            2 :        IF (wann%l_matrixmmn) THEN
    2436            2 :           DEALLOCATE (gb,bpt)
    2437            2 :           IF(l_gwf) DEALLOCATE (gb_q,bpt_q)
    2438              :        ENDIF
    2439              : 
    2440              : 
    2441              : 1911   CONTINUE
    2442              : 
    2443            4 :        IF(ALLOCATED(chi)) DEALLOCATE(chi)
    2444              : 
    2445            4 :        CALL timeStop("Wannier total")
    2446              :        CALL wann_postproc(&
    2447              :             stars,vacuum,atoms,sphhar,input,kpts,sym,fmpi,&
    2448              :             lapw ,noco,nococonv,cell,vTot,enpara,sliceplot,eig_id,l_real,&                !eig_id is used here after closing the files?!&
    2449            4 :             wann,fullnkpts,l_proj,results%ef,wann%l_sgwf,fullnqpts)
    2450              : 
    2451              : #ifdef CPP_MPI
    2452            4 :        CALL MPI_BARRIER(fmpi%mpi_comm,ierr(1))
    2453              : #endif
    2454              : 
    2455            4 :        IF(.NOT.wann%l_ldauwan) THEN
    2456            4 :           DO pc = 1, wann%nparampts
    2457            4 :              CALL close_eig(eig_idList(pc))
    2458              :           END DO
    2459              : 
    2460            4 :           CALL juDFT_end("wannier good",fmpi%irank)
    2461              :        END IF
    2462              : 
    2463            0 :      END SUBROUTINE wannier
    2464              : 
    2465              : 
    2466              :    END MODULE m_wannier
        

Generated by: LCOV version 2.0-1