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-2018 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 .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 module m_mlwfovlp_qp
28 
29  use defs_basis
30  use defs_datatypes
31  use defs_abitypes
32  use defs_wannier90
33  use m_errors
34  use m_abicore
35 
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_numeric_tools,  only : isordered
40  use m_geometry,       only : metric
41  use m_crystal,        only : crystal_t, crystal_free
42  use m_crystal_io,     only : crystal_from_hdr
43  use m_kpts,           only : listkk
44  use m_bz_mesh,        only : kmesh_t, kmesh_init, kmesh_free
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 
50  implicit none
51 
52  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_lda_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.

PARENTS

      mlwfovlp_qp

CHILDREN

      pawcprj_alloc,pawcprj_free

SOURCE

566 subroutine update_cprj(natom,nkibz,nbnds,nsppol,nspinor,m_lda_to_qp,dimlmn,Cprj_ibz)
567 
568 
569 !This section has been created automatically by the script Abilint (TD).
570 !Do not modify the following lines by hand.
571 #undef ABI_FUNC
572 #define ABI_FUNC 'update_cprj'
573 !End of the abilint section
574 
575  implicit none
576 
577 !Arguments ------------------------------------
578 !scalars
579  integer,intent(in) :: natom,nbnds,nkibz,nsppol,nspinor
580 !arrays
581  integer,intent(in) :: dimlmn(natom)
582  complex(dpc),intent(in) :: m_lda_to_qp(nbnds,nbnds,nkibz,nsppol)
583  type(pawcprj_type),intent(inout) :: Cprj_ibz(natom,nspinor*nbnds*nkibz*nsppol)
584 
585 !Local variables-------------------------------
586 !scalars
587  integer :: iat,ib,ik,is,shift,indx_kibz,ilmn,nlmn,ispinor,ibsp,spad,ibdx
588 !arrays
589  real(dp),allocatable :: re_p(:),im_p(:),vect(:,:),umat(:,:,:)
590  type(pawcprj_type),allocatable :: Cprj_ks(:,:)
591 
592 !************************************************************************
593 
594  DBG_ENTER("COLL")
595 
596  ABI_DATATYPE_ALLOCATE(Cprj_ks,(natom,nspinor*nbnds))
597  call pawcprj_alloc(Cprj_ks,0,dimlmn)
598 
599  ABI_ALLOCATE(re_p,(nbnds))
600  ABI_ALLOCATE(im_p,(nbnds))
601  ABI_ALLOCATE(vect,(2,nbnds))
602  ABI_ALLOCATE(umat,(2,nbnds,nbnds))
603  !
604  ! $ \Psi^{QP}_{r,b} = \sum_n \Psi^{KS}_{r,n} M_{n,b} $
605  !
606  ! therefore the updated PAW projections are given by:
607  !
608  ! $ \<\tprj_j|\Psi^{QP}_a\> = sum_b M_{b,a} <\tprj_j|\Psi^{KS}_b\> $.
609  !
610  do is=1,nsppol
611    do ik=1,nkibz
612 
613     shift=nspinor*nbnds*nkibz*(is-1)
614     indx_kibz=nspinor*nbnds*(ik-1)+shift
615     ibsp=0
616     do ib=1,nbnds
617       do ispinor=1,nspinor
618         ibsp=ibsp+1
619         do iat=1,natom
620           Cprj_ks(iat,ibsp)%cp(:,:)=Cprj_ibz(iat,indx_kibz+ibsp)%cp(:,:)
621         end do
622       end do
623     end do
624 
625     umat(1,:,:)=TRANSPOSE( REAL (m_lda_to_qp(:,:,ik,is)) )
626     umat(2,:,:)=TRANSPOSE( AIMAG(m_lda_to_qp(:,:,ik,is)) )
627 
628     do iat=1,natom
629       nlmn=dimlmn(iat)
630       do ilmn=1,nlmn
631 
632         do ispinor=1,nspinor
633            ! * Retrieve projections for this spinor component, at fixed atom and ilmn.
634            spad=(ispinor-1)
635            ibdx=0
636            do ib=1,nbnds*nspinor,nspinor
637             ibdx=ibdx+1
638             vect(1,ibdx)=Cprj_ks(iat,ib+spad)%cp(1,ilmn)
639             vect(2,ibdx)=Cprj_ks(iat,ib+spad)%cp(2,ilmn)
640            end do
641 
642            re_p(:)= &
643 &            MATMUL(umat(1,:,:),vect(1,:)) &
644 &           -MATMUL(umat(2,:,:),vect(2,:))
645 
646            im_p(:)= &
647 &            MATMUL(umat(1,:,:),vect(2,:)) &
648 &           +MATMUL(umat(2,:,:),vect(1,:))
649 
650            ! === Save values ===
651            ibdx=0
652            do ib=1,nbnds*nspinor,nspinor
653             ibdx=ibdx+1
654             Cprj_ibz(iat,indx_kibz+spad+ib)%cp(1,ilmn)=re_p(ibdx)
655             Cprj_ibz(iat,indx_kibz+spad+ib)%cp(2,ilmn)=im_p(ibdx)
656            end do
657         end do !ispinor
658 
659       end do !ilmn
660     end do !iat
661 
662    end do !ik
663  end do !is
664 
665  ABI_DEALLOCATE(re_p)
666  ABI_DEALLOCATE(im_p)
667  ABI_DEALLOCATE(vect)
668  ABI_DEALLOCATE(umat)
669 
670  call pawcprj_free(Cprj_ks)
671  ABI_DATATYPE_DEALLOCATE(Cprj_ks)
672 
673  DBG_EXIT("COLL")
674 
675 end subroutine update_cprj

m_mlwfovlp_qp/mlwfovlp_qp [ Functions ]

[ Top ] [ m_mlwfovlp_qp ] [ Functions ]

NAME

 mlwfovlp_qp

FUNCTION

 Routine which computes replaces LDA wave functions and eigenvalues with
 GW quasiparticle ones using previously computed qp wave functions in
 LDA 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 LDA 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.

PARENTS

      outscfcv

CHILDREN

      crystal_free,crystal_from_hdr,destroy_mpi_enreg,ebands_free,ebands_init
      initmpi_seq,kmesh_free,kmesh_init,listkk,metric,pawcprj_getdim,rdgw
      rdqps,sort_dp,update_cprj,wrtout,zgemm

SOURCE

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