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

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

Generated by: LCOV version 1.13