TABLE OF CONTENTS
ABINIT/m_mlwfovlp_qp [ 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 ]
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