TABLE OF CONTENTS


ABINIT/m_mlwfovlp_qp [ Modules ]

[ Top ] [ Modules ]

NAME

  m_mlwfovlp_qp

FUNCTION

  Interpolate GW corrections with Wannier functions

COPYRIGHT

  Copyright (C) 2008-2022 ABINIT group (DRH)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

SOURCE

17 #if defined HAVE_CONFIG_H
18 #include "config.h"
19 #endif
20 
21 #include "abi_common.h"
22 
23 module m_mlwfovlp_qp
24 
25  use defs_basis
26  use defs_wannier90
27  use m_errors
28  use m_abicore
29  use m_xmpi
30  use m_hdr
31  use m_dtset
32  use m_dtfil
33 
34  use defs_datatypes,   only : ebands_t
35  use defs_abitypes,    only : MPI_type
36  use m_mpinfo,         only : destroy_mpi_enreg, initmpi_seq
37  use m_pawtab,         only : pawtab_type
38  use m_pawcprj,        only : pawcprj_type, paw_overlap, pawcprj_getdim, pawcprj_alloc, pawcprj_free
39  use m_pawrhoij,       only : pawrhoij_type
40  use m_numeric_tools,  only : isordered
41  use m_geometry,       only : metric
42  use m_crystal,        only : crystal_t
43  use m_kpts,           only : listkk
44  use m_bz_mesh,        only : kmesh_t
45  use m_ebands,         only : ebands_init, ebands_free
46  use m_qparticles,     only : rdqps, rdgw
47  use m_sort,           only : sort_dp
48 
49  implicit none
50 
51  private

ABINIT/update_cprj [ Functions ]

[ Top ] [ Functions ]

NAME

 update_cprj

FUNCTION

  Update the matrix elements of the PAW projectors in case of self-consistent GW.

INPUTS

  dimlmn(natom)=number of (l,m,n) components for each atom (only for PAW)
  nkibz=number of k-points
  nsppol=number of spin
  nbnds=number of bands in the present GW calculation
  m_ks_to_qp(nbnds,nbnds,nkibz,nsppol)= expansion of the QP amplitudes in terms of KS wavefunctions
  natom=number of atomd in unit cell

OUTPUT

  Cprj_ibz(natom,nspinor*nkibz*nbnds*nsppol) <type(pawcprj_type)>=projected wave functions
   <Proj_i|Cnk> with all NL projectors. On exit, it contains the projections onto the
   QP amplitudes.

TODO

 To be moved to cprj_utils, although here we use complex variables.

SOURCE

547 subroutine update_cprj(natom,nkibz,nbnds,nsppol,nspinor,m_ks_to_qp,dimlmn,Cprj_ibz)
548 
549 !Arguments ------------------------------------
550 !scalars
551  integer,intent(in) :: natom,nbnds,nkibz,nsppol,nspinor
552 !arrays
553  integer,intent(in) :: dimlmn(natom)
554  complex(dpc),intent(in) :: m_ks_to_qp(nbnds,nbnds,nkibz,nsppol)
555  type(pawcprj_type),intent(inout) :: Cprj_ibz(natom,nspinor*nbnds*nkibz*nsppol)
556 
557 !Local variables-------------------------------
558 !scalars
559  integer :: iat,ib,ik,is,shift,indx_kibz,ilmn,nlmn,ispinor,ibsp,spad,ibdx
560 !arrays
561  real(dp),allocatable :: re_p(:),im_p(:),vect(:,:),umat(:,:,:)
562  type(pawcprj_type),allocatable :: Cprj_ks(:,:)
563 
564 !************************************************************************
565 
566  DBG_ENTER("COLL")
567 
568  ABI_MALLOC(Cprj_ks,(natom,nspinor*nbnds))
569  call pawcprj_alloc(Cprj_ks,0,dimlmn)
570 
571  ABI_MALLOC(re_p,(nbnds))
572  ABI_MALLOC(im_p,(nbnds))
573  ABI_MALLOC(vect,(2,nbnds))
574  ABI_MALLOC(umat,(2,nbnds,nbnds))
575  !
576  ! $ \Psi^{QP}_{r,b} = \sum_n \Psi^{KS}_{r,n} M_{n,b} $
577  !
578  ! therefore the updated PAW projections are given by:
579  !
580  ! $ \<\tprj_j|\Psi^{QP}_a\> = sum_b M_{b,a} <\tprj_j|\Psi^{KS}_b\> $.
581  !
582  do is=1,nsppol
583    do ik=1,nkibz
584 
585     shift=nspinor*nbnds*nkibz*(is-1)
586     indx_kibz=nspinor*nbnds*(ik-1)+shift
587     ibsp=0
588     do ib=1,nbnds
589       do ispinor=1,nspinor
590         ibsp=ibsp+1
591         do iat=1,natom
592           Cprj_ks(iat,ibsp)%cp(:,:)=Cprj_ibz(iat,indx_kibz+ibsp)%cp(:,:)
593         end do
594       end do
595     end do
596 
597     umat(1,:,:)=TRANSPOSE( REAL (m_ks_to_qp(:,:,ik,is)) )
598     umat(2,:,:)=TRANSPOSE( AIMAG(m_ks_to_qp(:,:,ik,is)) )
599 
600     do iat=1,natom
601       nlmn=dimlmn(iat)
602       do ilmn=1,nlmn
603 
604         do ispinor=1,nspinor
605            ! * Retrieve projections for this spinor component, at fixed atom and ilmn.
606            spad=(ispinor-1)
607            ibdx=0
608            do ib=1,nbnds*nspinor,nspinor
609             ibdx=ibdx+1
610             vect(1,ibdx)=Cprj_ks(iat,ib+spad)%cp(1,ilmn)
611             vect(2,ibdx)=Cprj_ks(iat,ib+spad)%cp(2,ilmn)
612            end do
613 
614            re_p(:)= &
615 &            MATMUL(umat(1,:,:),vect(1,:)) &
616 &           -MATMUL(umat(2,:,:),vect(2,:))
617 
618            im_p(:)= &
619 &            MATMUL(umat(1,:,:),vect(2,:)) &
620 &           +MATMUL(umat(2,:,:),vect(1,:))
621 
622            ! === Save values ===
623            ibdx=0
624            do ib=1,nbnds*nspinor,nspinor
625             ibdx=ibdx+1
626             Cprj_ibz(iat,indx_kibz+spad+ib)%cp(1,ilmn)=re_p(ibdx)
627             Cprj_ibz(iat,indx_kibz+spad+ib)%cp(2,ilmn)=im_p(ibdx)
628            end do
629         end do !ispinor
630 
631       end do !ilmn
632     end do !iat
633 
634    end do !ik
635  end do !is
636 
637  ABI_FREE(re_p)
638  ABI_FREE(im_p)
639  ABI_FREE(vect)
640  ABI_FREE(umat)
641 
642  call pawcprj_free(Cprj_ks)
643  ABI_FREE(Cprj_ks)
644 
645  DBG_EXIT("COLL")
646 
647 end subroutine update_cprj

m_mlwfovlp_qp/mlwfovlp_qp [ Functions ]

[ Top ] [ m_mlwfovlp_qp ] [ Functions ]

NAME

 mlwfovlp_qp

FUNCTION

 Routine which computes replaces DFT wave functions and eigenvalues with
 GW quasiparticle ones using previously computed qp wave functions in
 DFT bloch function representation for Wannier code (www.wannier.org f90 version).

INPUTS

  dtset <type(dataset_type)>=all input variables for this dataset
  dtfil <type(datafiles_type)>=variables related to files
  mband=maximum number of bands
  mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
  mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol
  mkmem =number of k points treated by this node.
  mpw=maximum dimensioned size of npw.
  natom=number of atoms in cell.
  nkpt=number of k points.
  npwarr(nkpt)=number of planewaves in basis at this k point
  nspden=number of spin-density components
  nsppol=1 for unpolarized, 2 for spin-polarized
  rprimd(3,3)=dimensional primitive translations for real space (bohr)
  Hdr<Hdr_type>=The m_mlwfovlp_qp header.
  MPI_enreg=information about MPI parallelization
  Cprj_BZ(natom,nspinor*mband*mkmem*nsppol)= <p_lmn|Cnk> coefficients for each WF |Cnk> and each |p_lmn> non-local projector

OUTPUT

SIDE EFFECTS

  cg(2,mcg)=planewave coefficients of wavefunctions
   replaced by quasiparticle wavefunctions
  eigen(mband*nkpt*nsppol)=array for holding eigenvalues replaced by qp eigenvalues(hartree)

NOTES

  Number of bands for wannier calculation must be identical to number used
   for gw calculation.  Bands not wanted for wannier calculation must be
   excluded in exclude_band statement in wannier90.win file.
  Full plane-wave basis for DFT wavefunctions must be used in GW calculation,
   or inaccuracies may result.
  This is at best a beta version of this code, with little consistency
   checking, so the user must be very careful or the results may be invalid.

SOURCE

105 subroutine mlwfovlp_qp(cg,Cprj_BZ,dtset,dtfil,eigen,mband,mcg,mcprj,mkmem,mpw,natom,&
106 & nkpt,npwarr,nspden,nsppol,ntypat,Hdr,Pawtab,rprimd,MPI_enreg)
107 
108 !Arguments ------------------------------------
109 !scalars
110  integer,intent(in) :: mband,mcg,mcprj,mkmem,mpw,nkpt,nspden,natom,ntypat
111  integer,intent(in) :: nsppol
112  type(dataset_type),intent(in) :: dtset
113  type(datafiles_type),intent(in) :: dtfil
114  type(Hdr_type),intent(in) :: Hdr
115  type(MPI_type),intent(in) :: MPI_enreg
116  type(pawcprj_type),target,intent(inout) :: Cprj_BZ(natom,mcprj)
117  type(Pawtab_type),intent(in) :: Pawtab(ntypat*Dtset%usepaw)
118 !arrays
119  integer,intent(in) :: npwarr(nkpt)
120  real(dp),intent(inout) :: cg(2,mcg)
121  real(dp),intent(inout) :: eigen(mband*nkpt*nsppol)
122  real(dp),intent(in) :: rprimd(3,3)
123 
124 !Local variables-------------------------------
125 !scalars
126  integer,parameter :: from_QPS_FILE=1,from_GW_FILE=2
127  integer :: sppoldbl,timrev,bantot_ibz,ikibz,ikbz,dimrho
128  integer :: iband,icg,icg_shift,ii,ipw,isppol,my_nspinor,nband_k,ord_iband
129  integer :: nfftot,ikpt,irzkpt,npw_k,ikg
130  integer :: nscf,nbsc,itimrev,band_index,nkibz,nkbz
131  integer :: input !,jb_idx,ib_idx,ijpack, jband,
132  integer :: nprocs,ios
133  real(dp) :: TOL_SORT=tol12
134  real(dp) :: dksqmax,ucvol !ortho_err,
135  logical :: ltest,qpenek_is_ordered,g0w0_exists
136  character(len=500) :: msg
137  character(len=fnlen) :: gw_fname
138  type(ebands_t) :: QP_bst
139  type(crystal_t)  :: Cryst
140  type(kmesh_t) :: Kibz_mesh
141  type(MPI_type) :: MPI_enreg_seq
142 !arrays
143  integer :: indkk(nkpt,6),my_ngfft(18)
144  integer,allocatable :: npwarr_ibz(:),nband_ibz(:),ibz2bz(:,:),istwfk_ibz(:)
145  integer,allocatable :: dimlmn(:),iord(:),nattyp_dum(:)
146  real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3) !,paw_ovlp(2)
147  real(dp),allocatable :: qp_rhor(:,:),sorted_qpene(:)
148  real(dp),allocatable :: kibz(:,:),wtk_ibz(:)
149  real(dp),allocatable :: doccde_ibz(:),occfact_ibz(:),eigen_ibz(:)
150  real(dp),allocatable ::  igwene(:,:,:)
151  complex(dpc),allocatable :: m_ks_to_qp(:,:,:,:),m_ks_to_qp_BZ(:,:,:,:) !,ortho(:)
152  complex(dpc),allocatable :: m_tmp(:,:),cg_k(:,:),cg_qpk(:,:)
153  type(Pawrhoij_type),allocatable :: prev_Pawrhoij(:)
154  !type(pawcprj_type),pointer :: Cp1(:,:),Cp2(:,:)
155 
156 !************************************************************************
157 
158  ABI_UNUSED(mkmem)
159 
160  DBG_ENTER("COLL")
161 
162  write(msg,'(17a)')ch10,&
163   ' mlwfovlp_qp: WARNING',ch10,&
164   '  The input *_WFK file of DFT wavefunctions to be  converted',ch10,&
165   '  to GW quasiparticle wavefunctions MUST have been written in',ch10,&
166   '  the run that produced the GW *_KSS file using kssform 3,',ch10,&
167   '  the ONLY value of kssform permitted for GW Wannier functions.',ch10,&
168   '  Otherwise, the *_QPS file needed here will be inconsistent,',ch10,&
169   '  and the output quasiparticle wavefunctions will be garbage.',ch10,&
170   '  No internal check that can verify this is presently possible.',ch10
171  call wrtout(std_out,msg,'COLL')
172 
173  ! === Some features are not implemented yet ===
174  ABI_CHECK(Dtset%nspinor==1,'nspinor==2 not implemented')
175  ABI_CHECK(Dtset%nsppol==1,'nsppol==2 not implemented, check wannier90')
176  ltest=ALL(Dtset%nband(1:Dtset%nkpt*Dtset%nsppol)==Dtset%nband(1))
177  ABI_CHECK(ltest,'nband(:) should be constant')
178  !
179  ! MPI initialization
180  nprocs=MPI_enreg%nproc_cell
181 
182  if (nprocs/=1) then
183    ABI_ERROR("mlwfovlp_qp not programmed for parallel execution")
184  end if
185 
186  ! Compute reciprocal space metric gmet for unit cell of disk wf
187  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
188 
189  ! Compute k points from gw irreducible set equivalent to full-zone wannier set
190  sppoldbl=1 ; timrev=1 ; my_nspinor=max(1,Dtset%nspinor/MPI_enreg%nproc_spinor)
191  call listkk(dksqmax,gmet,indkk,dtset%kptgw,dtset%kpt,dtset%nkptgw,nkpt,&
192 &  dtset%nsym,sppoldbl,dtset%symafm,dtset%symrel,timrev,xmpi_comm_self)
193 
194  if (dksqmax>tol8) then
195    write(msg,'(5a)')&
196 &    'Set of GW irreducible-zone kptgw in input file is inconsistent',ch10,&
197 &    'with full-zone set being used for wannier90 setup.',ch10,&
198 &    'Action: correct input data'
199    ABI_ERROR(msg)
200  end if
201  !
202  ! === Initialize object defining the Band strucuture ===
203  ! * Initialize with KS results using IBZ indexing.
204  ! * After rdqps, QP_bst will contain the QP amplitudes.
205  nkibz=Dtset%nkptgw
206  ABI_MALLOC(kibz,(3,nkibz))
207  ABI_MALLOC(wtk_ibz,(nkibz))
208  kibz=Dtset%kptgw(:,1:Dtset%nkptgw)
209 
210  ! MG: This part is needed to get the IBZ weight that will be reported
211  ! on ab_out thus we should be consistent. Ideally Cryst should be
212  ! one of the basic abinit objects and it should be passed to this routine.
213 
214  !different conventions are used in GW and abinit!!
215  cryst = hdr%get_crystal(gw_timrev=timrev+1)
216  call Kibz_mesh%init(Cryst,nkibz,kibz,Dtset%kptopt)
217  wtk_ibz=Kibz_mesh%wt
218  call cryst%free()
219  call Kibz_mesh%free()
220 
221  ABI_MALLOC(ibz2bz,(nkibz,6))
222  call listkk(dksqmax,gmet,ibz2bz,dtset%kpt,dtset%kptgw,nkpt,dtset%nkptgw,&
223 &  dtset%nsym,sppoldbl,dtset%symafm,dtset%symrel,timrev,xmpi_comm_self)
224 
225  ltest=ALL(ibz2bz(:,2)==1)
226  ABI_CHECK(ltest,'Not able to found irreducible points in the BZ set!')
227 
228  if (dksqmax>tol8) then
229     write(msg,'(5a)')&
230      'Set of GW irreducible-zone kptgw in input file is inconsistent',ch10,&
231      'with full-zone set being used for wannier90 setup.',ch10,&
232      'Action: correct input data'
233     ABI_ERROR(msg)
234  end if
235 
236  ABI_MALLOC(npwarr_ibz,(nkibz))
237  ABI_MALLOC(istwfk_ibz,(nkibz))
238  ABI_MALLOC(nband_ibz,(nkibz*nsppol))
239 
240  do isppol=1,nsppol
241    do ikibz=1,nkibz
242      ikbz=ibz2bz(ikibz+(sppoldbl-1)*(isppol-1)*nkibz,1)
243      npwarr_ibz(ikibz)=      npwarr(ikbz)
244      istwfk_ibz(ikibz)=Dtset%istwfk(ikbz)
245      nband_ibz(ikibz+(isppol-1)*nkibz)=Dtset%nband(ikbz+(isppol-1)*nkpt)
246    end do
247  end do
248 
249  bantot_ibz=SUM(nband_ibz)
250  ABI_MALLOC(doccde_ibz,(bantot_ibz))
251  ABI_MALLOC(eigen_ibz,(bantot_ibz))
252  ABI_MALLOC(occfact_ibz,(bantot_ibz))
253  doccde_ibz(:)=zero ; eigen_ibz(:)=zero ; occfact_ibz(:)=zero
254 
255  band_index=0
256  do isppol=1,nsppol
257    do ikibz=1,nkibz
258      ikbz=ibz2bz(ikibz+(sppoldbl-1)*(isppol-1)*nkibz,1)
259      nband_k=nband_ibz(ikibz+(isppol-1)*nkibz)
260      ii=SUM(Dtset%nband(1:ikbz+(isppol-1)*nkpt))-nband_k
261      eigen_ibz(band_index+1:band_index+nband_k)=eigen(ii+1:ii+nband_k)
262      band_index=band_index+nband_k
263    end do
264  end do
265 
266  ! CP modified
267 ! call ebands_init(bantot_ibz,QP_bst,Dtset%nelect,doccde_ibz,eigen_ibz,istwfk_ibz,kibz,nband_ibz,&
268 !  nkibz,npwarr_ibz,nsppol,Dtset%nspinor,Dtset%tphysel,Dtset%tsmear,Dtset%occopt,occfact_ibz,wtk_ibz,&
269 !  dtset%cellcharge(1),dtset%kptopt,dtset%kptrlatt_orig,dtset%nshiftk_orig,dtset%shiftk_orig,&
270 !  dtset%kptrlatt,dtset%nshiftk,dtset%shiftk)
271  call ebands_init(bantot_ibz,QP_bst,Dtset%nelect,Dtset%ne_qFD,Dtset%nh_qFD,Dtset%ivalence,&
272   doccde_ibz,eigen_ibz,istwfk_ibz,kibz,nband_ibz,&
273   nkibz,npwarr_ibz,nsppol,Dtset%nspinor,Dtset%tphysel,Dtset%tsmear,Dtset%occopt,occfact_ibz,wtk_ibz,&
274   dtset%cellcharge(1),dtset%kptopt,dtset%kptrlatt_orig,dtset%nshiftk_orig,dtset%shiftk_orig,&
275   dtset%kptrlatt,dtset%nshiftk,dtset%shiftk)
276  ! End CP modified
277 
278  ABI_FREE(kibz)
279  ABI_FREE(wtk_ibz)
280  ABI_FREE(ibz2bz)
281  ABI_FREE(npwarr_ibz)
282  ABI_FREE(istwfk_ibz)
283  ABI_FREE(nband_ibz)
284  ABI_FREE(doccde_ibz)
285  ABI_FREE(eigen_ibz)
286  ABI_FREE(occfact_ibz)
287 
288  ! === Read in quasiparticle information ===
289  ! * Initialize QP amplitudes with KS, QP_bst% presently contains KS energies.
290  ! * If file not found return, everything has been already initialized with KS values
291  !   Here qp_rhor is not needed thus dimrho=0
292  ABI_MALLOC(m_ks_to_qp,(mband,mband,dtset%nkptgw,nsppol))
293  m_ks_to_qp=czero
294  do iband=1,mband
295    m_ks_to_qp(iband,iband,:,:)=cone
296  end do
297 
298  ! Fake MPI_type for rdqps
299  call initmpi_seq(MPI_enreg_seq)
300 
301  my_ngfft=Dtset%ngfft; if (Dtset%usepaw==1.and.ALL(Dtset%ngfftdg(1:3)/=0)) my_ngfft=Dtset%ngfftdg
302  nfftot=PRODUCT(my_ngfft(1:3)); dimrho=0
303 
304  ! Change gw_fname to read a GW file instead of the QPS file.
305  ! TODO not so sure that wannier90 can handle G0W0 eigenvalues that are not ordered, though!
306  gw_fname = "g0w0"
307  g0w0_exists = .FALSE.
308  inquire(file=gw_fname,iostat=ios,exist=g0w0_exists)
309  if (ios/=0) then
310    ABI_ERROR('File g0w0 exists but iostat returns nonzero value!')
311  end if
312 
313  if (.not.g0w0_exists) then ! read QPS file (default behavior).
314    input = from_QPS_FILE
315    ABI_MALLOC(prev_Pawrhoij,(Cryst%natom*Dtset%usepaw))
316    ABI_MALLOC(qp_rhor,(nfftot,nspden*dimrho))
317 
318    call rdqps(QP_bst,Dtfil%fnameabi_qps,Dtset%usepaw,Dtset%nspden,dimrho,nscf,&
319     nfftot,my_ngfft,ucvol,Cryst,Pawtab,MPI_enreg_seq,nbsc,m_ks_to_qp,qp_rhor,prev_Pawrhoij)
320 
321    ABI_FREE(qp_rhor)
322    ABI_FREE(prev_Pawrhoij)
323 
324  else
325    ! Read GW file (m_ks_to_qp has been already set to 1, no extrapolation is performed)
326    ABI_WARNING(' READING GW CORRECTIONS FROM FILE g0w0 !')
327    input = from_GW_FILE
328    ABI_MALLOC(igwene,(QP_bst%mband,QP_bst%nkpt,QP_bst%nsppol))
329    call rdgw(QP_bst,gw_fname,igwene,extrapolate=.FALSE.)
330    ABI_FREE(igwene)
331  end if
332 
333  ! === Begin big loop over full-zone k points and spin (not implemented) ===
334  ! * Wannier90 treats only a single spin, changes in wannier90 are needed
335  ABI_MALLOC(cg_k,(mpw,mband))
336  ABI_MALLOC(cg_qpk,(mpw,mband))
337  ABI_MALLOC(m_tmp,(mband,mband))
338 
339  band_index=0 ; icg=0 ; ikg=0
340  do isppol=1,nsppol
341    do ikpt=1,nkpt
342 
343     irzkpt =indkk(ikpt+(sppoldbl-1)*(isppol-1)*nkpt,1)
344     itimrev=indkk(ikpt+(sppoldbl-1)*(isppol-1)*nkpt,6)
345     npw_k=npwarr(ikpt)
346     nband_k=dtset%nband(ikpt+(isppol-1)*nkpt)
347 
348     if (nband_k/=mband) then
349       write(msg,'(a,i0,7a)')&
350        'Number of bands for k point ',ikpt,' is inconsistent with number',ch10,&
351        'specified for wannier90 calculation',ch10,&
352        'Action: correct input so all band numbers are equal for GW',ch10,&
353        'and wannier90 datasets.'
354       ABI_ERROR(msg)
355     end if
356 
357     ! Load KS states for this kbz and spin
358     do iband=1,nband_k
359       icg_shift=npw_k*my_nspinor*(iband-1)+icg
360       do ipw=1,npw_k
361         cg_k(ipw,iband)=DCMPLX(cg(1,ipw+icg_shift),cg(2,ipw+icg_shift))
362       end do
363     end do
364 
365     ! If time reversal is used for relating ikpt to irzkpt, then multiply by
366     ! the complex conjugate of the lda-to-qp transformation matrix
367     if (itimrev==0) then
368       m_tmp(:,:)=m_ks_to_qp(:,:,irzkpt,isppol)
369     else if (itimrev==1) then
370       m_tmp(:,:)=conjg(m_ks_to_qp(:,:,irzkpt,isppol))
371     else
372       write(msg,'(2(a,i0))')'Invalid indkk(ikpt,6) ',itimrev,'from routine listkk for k-point ',ikpt
373       ABI_BUG(msg)
374     end if
375 
376     call ZGEMM('N','N',npw_k,mband,mband,cone,cg_k,mpw,m_tmp,mband,czero,cg_qpk,mpw)
377 
378     ! === Orthonormality test ===
379     ! * nband >= maxval(bndgw) for this to pass, but may be less than nband used in GW.
380     ! * Unfortunately, does not test WFK and QPS consistency.
381     !allocate(ortho(nband_k*(nband_k+1)/2))
382     !ortho=czero; ijpack=0
383     !do jband=1,nband_k
384     !  jb_idx=band_index+jband
385     !  if (dtset%usepaw==1) Cp2 => Cprj_BZ(:,jband:jband+(my_nspinor-1))
386     !  do iband=1,jband
387     !    ib_idx=band_index+iband
388     !    ijpack=ijpack+1
389     !    ortho(ijpack)=sum(conjg(cg_qpk(1:npw_k,iband))*cg_qpk(1:npw_k,jband))
390     !    if (dtset%usepaw==1) then
391     !      Cp1 => Cprj_BZ(:,iband:iband+(my_nspinor-1))
392     !      paw_ovlp = paw_overlap(Cp2,Cp1,Cryst%typat,Pawtab)
393     !      ortho(ijpack) = ortho(ijpack) + CMPLX(paw_ovlp(1),paw_ovlp(2))
394     !    end if
395     !    if (jband==iband) ortho(ijpack)=ortho(ijpack)-cone
396     !  end do
397     !end do
398     !ortho_err=maxval(abs(ortho))
399 
400     !write(std_out,*)' drh - mlwfovlp_qp: ikpt,ortho_err',ikpt,ortho_err
401     !if (ortho_err>tol6) then
402     !  write(msg, '(3a,i4,a,i6,a,1p,e8.1,3a)' )&
403     !&    '  orthonormality error for quasiparticle wave functions.',ch10,&
404     !&    '  spin=',isppol,'  k point=',ikpt,'  ortho_err=',ortho_err,' >1E-6',ch10,&
405     !&    '  Action: Be sure input nband>=maxval(bndgw)'
406     !  ABI_ERROR(msg)
407     !end if
408     !deallocate(ortho)
409 
410     ! Replace lda wave functions and eigenvalues with quasiparticle ones.
411     qpenek_is_ordered = isordered(nband_k,QP_bst%eig(:,irzkpt,isppol),">",TOL_SORT)
412 
413     if (input==from_QPS_FILE .and. .not.qpenek_is_ordered) then
414       write(msg,'(3a)')&
415       " QP energies read from QPS file are not ordered, likely nband_k>nbdgw. ",ch10,&
416       " Change nband in the input file so that it equals the number of GW states calculated"
417       ABI_WARNING(msg)
418     end if
419 
420     if ( .TRUE. ) then
421       do iband=1,nband_k
422         icg_shift=npw_k*my_nspinor*(iband-1)+icg
423         eigen(iband+band_index)=QP_bst%eig(iband,irzkpt,isppol)
424         do ipw=1,npw_k
425           cg(1,ipw+icg_shift)= real(cg_qpk(ipw,iband))
426           cg(2,ipw+icg_shift)=aimag(cg_qpk(ipw,iband))
427         end do
428       end do
429     else
430       ! FIXME There's a problem in twannier90 since nband_k > nbdgw and therefore we also read KS states from the QPS file!
431       ! Automatic test has to be changed!
432       write(msg,'(2a,3f8.4,3a)')ch10,&
433         "QP energies at k-point ",QP_bst%kptns(:,irzkpt)," are not sorted in ascending numerical order!",ch10,&
434         "Performing reordering of energies and wavefunctions to be written on the final WKF file."
435       ABI_ERROR(msg)
436       !write(std_out,*)"eig",(QP_bst%eig(ii,irzkpt,isppol),ii=1,nband_k)
437       ABI_MALLOC(sorted_qpene,(nband_k))
438       ABI_MALLOC(iord,(nband_k))
439       sorted_qpene = QP_bst%eig(1:nband_k,irzkpt,isppol)
440       iord = (/(ii, ii=1,nband_k)/)
441 
442       call sort_dp(nband_k,sorted_qpene,iord,TOL_SORT)
443       do ii=1,nband_k
444         write(std_out,*)"%eig, sorted_qpene, iord",QP_bst%eig(ii,irzkpt,isppol)*Ha_eV,sorted_qpene(ii)*Ha_eV,iord(ii)
445       end do
446 
447       do iband=1,nband_k
448         ord_iband = iord(iband)
449         icg_shift=npw_k*my_nspinor*(iband-1)+icg
450         !eigen(iband+band_index)=QP_bst%eig(iband,irzkpt,isppol)
451         eigen(iband+band_index)=QP_bst%eig(ord_iband,irzkpt,isppol)
452         do ipw=1,npw_k
453           !cg(1,ipw+icg_shift)= real(cg_qpk(ipw,iband))
454           !cg(2,ipw+icg_shift)=aimag(cg_qpk(ipw,iband))
455           cg(1,ipw+icg_shift)= real(cg_qpk(ipw,ord_iband))
456           cg(2,ipw+icg_shift)=aimag(cg_qpk(ipw,ord_iband))
457         end do
458       end do
459       ABI_FREE(sorted_qpene)
460       ABI_FREE(iord)
461     end if
462 
463     band_index=band_index+nband_k
464     icg=icg+npw_k*my_nspinor*nband_k
465     ikg=ikg+npw_k
466    end do !ikpt
467  end do !isppol
468 
469  ABI_FREE(cg_k)
470  ABI_FREE(cg_qpk)
471  ABI_FREE(m_tmp)
472 
473  ! === If PAW, update projections in BZ ===
474  ! * Since I am lazy and here I do not care about memory, I just reconstruct m_ks_to_qp in the BZ.
475  ! * update_cprj will take care of updating the PAW projections to get <p_lmn|QP_{nks]>
476  !   This allows some CPU saving, no need to call ctocprj.
477  ! FIXME this part should be tested, automatic test to be provided
478  if (Dtset%usepaw==1) then
479    ABI_MALLOC(dimlmn,(natom))
480    call pawcprj_getdim(dimlmn,dtset%natom,nattyp_dum,ntypat,Dtset%typat,pawtab,'R')
481 
482    nkbz=nkpt
483    ABI_MALLOC(m_ks_to_qp_BZ,(mband,mband,nkbz,nsppol))
484    do isppol=1,nsppol
485      do ikbz=1,nkbz
486        ikibz  =indkk(ikibz+(sppoldbl-1)*(isppol-1)*nkbz,1)
487        itimrev=indkk(ikibz+(sppoldbl-1)*(isppol-1)*nkbz,6)
488        select case (itimrev)
489        case (0)
490          m_ks_to_qp_BZ(:,:,ikbz,isppol)=m_ks_to_qp(:,:,ikibz,isppol)
491        case (1)
492          m_ks_to_qp_BZ(:,:,ikbz,isppol)=CONJG(m_ks_to_qp(:,:,ikibz,isppol))
493        case default
494          write(msg,'(a,i3)')"Wrong itimrev= ",itimrev
495          ABI_BUG(msg)
496        end select
497      end do
498    end do
499 
500    call update_cprj(natom,nkbz,mband,nsppol,my_nspinor,m_ks_to_qp_BZ,dimlmn,Cprj_BZ)
501    ABI_FREE(dimlmn)
502    ABI_FREE(m_ks_to_qp_BZ)
503  end if !PAW
504 
505  write(msg,'(6a)')ch10,&
506   ' mlwfovlp_qp: Input KS wavefuctions have been converted',ch10,&
507   '  to GW quasiparticle wavefunctions for maximally localized wannier',ch10,&
508   '  function construction by wannier90.'
509  call wrtout(ab_out,msg,'COLL')
510  call wrtout(std_out,msg,'COLL')
511 
512  ABI_FREE(m_ks_to_qp)
513  call ebands_free(QP_bst)
514  call destroy_mpi_enreg(MPI_enreg_seq)
515 
516  DBG_EXIT("COLL")
517 
518 end subroutine mlwfovlp_qp