TABLE OF CONTENTS


ABINIT/mlwfovlp_qp [ Functions ]

[ Top ] [ 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).

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 .

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 abinit 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

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