TABLE OF CONTENTS
ABINIT/m_mlwfovlp_qp [ Modules ]
NAME
m_mlwfovlp_qp
FUNCTION
Interpolate GW corrections with Wannier functions
COPYRIGHT
Copyright (C) 2008-2024 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
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_mlwfovlp_qp 23 24 use defs_basis 25 use defs_wannier90 26 use m_errors 27 use m_abicore 28 use m_xmpi 29 use m_hdr 30 use m_dtset 31 use m_dtfil 32 33 use defs_datatypes, only : ebands_t 34 use defs_abitypes, only : MPI_type 35 use m_mpinfo, only : destroy_mpi_enreg, initmpi_seq 36 use m_pawtab, only : pawtab_type 37 use m_pawcprj, only : pawcprj_type, paw_overlap, pawcprj_getdim, pawcprj_alloc, pawcprj_free 38 use m_pawrhoij, only : pawrhoij_type 39 use m_numeric_tools, only : isordered 40 use m_geometry, only : metric 41 use m_crystal, only : crystal_t 42 use m_kpts, only : listkk 43 use m_bz_mesh, only : kmesh_t 44 use m_ebands, only : ebands_init, ebands_free 45 use m_qparticles, only : rdqps, rdgw 46 use m_sort, only : sort_dp 47 48 implicit none 49 50 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
540 subroutine update_cprj(natom,nkibz,nbnds,nsppol,nspinor,m_ks_to_qp,dimlmn,Cprj_ibz) 541 542 !Arguments ------------------------------------ 543 !scalars 544 integer,intent(in) :: natom,nbnds,nkibz,nsppol,nspinor 545 !arrays 546 integer,intent(in) :: dimlmn(natom) 547 complex(dpc),intent(in) :: m_ks_to_qp(nbnds,nbnds,nkibz,nsppol) 548 type(pawcprj_type),intent(inout) :: Cprj_ibz(natom,nspinor*nbnds*nkibz*nsppol) 549 550 !Local variables------------------------------- 551 !scalars 552 integer :: iat,ib,ik,is,shift,indx_kibz,ilmn,nlmn,ispinor,ibsp,spad,ibdx 553 !arrays 554 real(dp),allocatable :: re_p(:),im_p(:),vect(:,:),umat(:,:,:) 555 type(pawcprj_type),allocatable :: Cprj_ks(:,:) 556 557 !************************************************************************ 558 559 DBG_ENTER("COLL") 560 561 ABI_MALLOC(Cprj_ks,(natom,nspinor*nbnds)) 562 call pawcprj_alloc(Cprj_ks,0,dimlmn) 563 564 ABI_MALLOC(re_p,(nbnds)) 565 ABI_MALLOC(im_p,(nbnds)) 566 ABI_MALLOC(vect,(2,nbnds)) 567 ABI_MALLOC(umat,(2,nbnds,nbnds)) 568 ! 569 ! $ \Psi^{QP}_{r,b} = \sum_n \Psi^{KS}_{r,n} M_{n,b} $ 570 ! 571 ! therefore the updated PAW projections are given by: 572 ! 573 ! $ \<\tprj_j|\Psi^{QP}_a\> = sum_b M_{b,a} <\tprj_j|\Psi^{KS}_b\> $. 574 ! 575 do is=1,nsppol 576 do ik=1,nkibz 577 578 shift=nspinor*nbnds*nkibz*(is-1) 579 indx_kibz=nspinor*nbnds*(ik-1)+shift 580 ibsp=0 581 do ib=1,nbnds 582 do ispinor=1,nspinor 583 ibsp=ibsp+1 584 do iat=1,natom 585 Cprj_ks(iat,ibsp)%cp(:,:)=Cprj_ibz(iat,indx_kibz+ibsp)%cp(:,:) 586 end do 587 end do 588 end do 589 590 umat(1,:,:)=TRANSPOSE( REAL (m_ks_to_qp(:,:,ik,is)) ) 591 umat(2,:,:)=TRANSPOSE( AIMAG(m_ks_to_qp(:,:,ik,is)) ) 592 593 do iat=1,natom 594 nlmn=dimlmn(iat) 595 do ilmn=1,nlmn 596 597 do ispinor=1,nspinor 598 ! * Retrieve projections for this spinor component, at fixed atom and ilmn. 599 spad=(ispinor-1) 600 ibdx=0 601 do ib=1,nbnds*nspinor,nspinor 602 ibdx=ibdx+1 603 vect(1,ibdx)=Cprj_ks(iat,ib+spad)%cp(1,ilmn) 604 vect(2,ibdx)=Cprj_ks(iat,ib+spad)%cp(2,ilmn) 605 end do 606 607 re_p(:)= & 608 & MATMUL(umat(1,:,:),vect(1,:)) & 609 & -MATMUL(umat(2,:,:),vect(2,:)) 610 611 im_p(:)= & 612 & MATMUL(umat(1,:,:),vect(2,:)) & 613 & +MATMUL(umat(2,:,:),vect(1,:)) 614 615 ! === Save values === 616 ibdx=0 617 do ib=1,nbnds*nspinor,nspinor 618 ibdx=ibdx+1 619 Cprj_ibz(iat,indx_kibz+spad+ib)%cp(1,ilmn)=re_p(ibdx) 620 Cprj_ibz(iat,indx_kibz+spad+ib)%cp(2,ilmn)=im_p(ibdx) 621 end do 622 end do !ispinor 623 624 end do !ilmn 625 end do !iat 626 627 end do !ik 628 end do !is 629 630 ABI_FREE(re_p) 631 ABI_FREE(im_p) 632 ABI_FREE(vect) 633 ABI_FREE(umat) 634 635 call pawcprj_free(Cprj_ks) 636 ABI_FREE(Cprj_ks) 637 638 DBG_EXIT("COLL") 639 640 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
104 subroutine mlwfovlp_qp(cg,Cprj_BZ,dtset,dtfil,eigen,mband,mcg,mcprj,mkmem,mpw,natom,& 105 & nkpt,npwarr,nspden,nsppol,ntypat,Hdr,Pawtab,rprimd,MPI_enreg) 106 107 !Arguments ------------------------------------ 108 !scalars 109 integer,intent(in) :: mband,mcg,mcprj,mkmem,mpw,nkpt,nspden,natom,ntypat 110 integer,intent(in) :: nsppol 111 type(dataset_type),intent(in) :: dtset 112 type(datafiles_type),intent(in) :: dtfil 113 type(Hdr_type),intent(in) :: Hdr 114 type(MPI_type),intent(in) :: MPI_enreg 115 type(pawcprj_type),target,intent(inout) :: Cprj_BZ(natom,mcprj) 116 type(Pawtab_type),intent(in) :: Pawtab(ntypat*Dtset%usepaw) 117 !arrays 118 integer,intent(in) :: npwarr(nkpt) 119 real(dp),intent(inout) :: cg(2,mcg) 120 real(dp),intent(inout) :: eigen(mband*nkpt*nsppol) 121 real(dp),intent(in) :: rprimd(3,3) 122 123 !Local variables------------------------------- 124 !scalars 125 integer,parameter :: from_QPS_FILE=1,from_GW_FILE=2 126 integer :: sppoldbl,timrev,bantot_ibz,ikibz,ikbz,dimrho 127 integer :: iband,icg,icg_shift,ii,ipw,isppol,my_nspinor,nband_k,ord_iband 128 integer :: nfftot,ikpt,irzkpt,npw_k,ikg 129 integer :: nscf,nbsc,itimrev,band_index,nkibz,nkbz 130 integer :: input !,jb_idx,ib_idx,ijpack, jband, 131 integer :: nprocs,ios 132 real(dp) :: TOL_SORT=tol12 133 real(dp) :: dksqmax,ucvol !ortho_err, 134 logical :: ltest,qpenek_is_ordered,g0w0_exists 135 character(len=500) :: msg 136 character(len=fnlen) :: gw_fname 137 type(ebands_t) :: QP_bst 138 type(crystal_t) :: Cryst 139 type(kmesh_t) :: Kibz_mesh 140 type(MPI_type) :: MPI_enreg_seq 141 !arrays 142 integer :: indkk(nkpt,6),my_ngfft(18) 143 integer,allocatable :: npwarr_ibz(:),nband_ibz(:),ibz2bz(:,:),istwfk_ibz(:) 144 integer,allocatable :: dimlmn(:),iord(:),nattyp_dum(:) 145 real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3) !,paw_ovlp(2) 146 real(dp),allocatable :: qp_rhor(:,:),sorted_qpene(:) 147 real(dp),allocatable :: kibz(:,:),wtk_ibz(:) 148 real(dp),allocatable :: doccde_ibz(:),occfact_ibz(:),eigen_ibz(:) 149 real(dp),allocatable :: igwene(:,:,:) 150 complex(dpc),allocatable :: m_ks_to_qp(:,:,:,:),m_ks_to_qp_BZ(:,:,:,:) !,ortho(:) 151 complex(dpc),allocatable :: m_tmp(:,:),cg_k(:,:),cg_qpk(:,:) 152 type(Pawrhoij_type),allocatable :: prev_Pawrhoij(:) 153 !type(pawcprj_type),pointer :: Cp1(:,:),Cp2(:,:) 154 155 !************************************************************************ 156 157 ABI_UNUSED(mkmem) 158 159 DBG_ENTER("COLL") 160 161 write(msg,'(17a)')ch10,& 162 ' mlwfovlp_qp: WARNING',ch10,& 163 ' The input *_WFK file of DFT wavefunctions to be converted',ch10,& 164 ' to GW quasiparticle wavefunctions MUST have been written in',ch10,& 165 ' the run that produced the GW *_KSS file using kssform 3,',ch10,& 166 ' the ONLY value of kssform permitted for GW Wannier functions.',ch10,& 167 ' Otherwise, the *_QPS file needed here will be inconsistent,',ch10,& 168 ' and the output quasiparticle wavefunctions will be garbage.',ch10,& 169 ' No internal check that can verify this is presently possible.',ch10 170 call wrtout(std_out,msg,'COLL') 171 172 ! === Some features are not implemented yet === 173 ABI_CHECK(Dtset%nspinor==1,'nspinor==2 not implemented') 174 ABI_CHECK(Dtset%nsppol==1,'nsppol==2 not implemented, check wannier90') 175 ltest=ALL(Dtset%nband(1:Dtset%nkpt*Dtset%nsppol)==Dtset%nband(1)) 176 ABI_CHECK(ltest,'nband(:) should be constant') 177 ! 178 ! MPI initialization 179 nprocs=MPI_enreg%nproc_cell 180 181 if (nprocs/=1) then 182 ABI_ERROR("mlwfovlp_qp not programmed for parallel execution") 183 end if 184 185 ! Compute reciprocal space metric gmet for unit cell of disk wf 186 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 187 188 ! Compute k points from gw irreducible set equivalent to full-zone wannier set 189 sppoldbl=1 ; timrev=1 ; my_nspinor=max(1,Dtset%nspinor/MPI_enreg%nproc_spinor) 190 call listkk(dksqmax,gmet,indkk,dtset%kptgw,dtset%kpt,dtset%nkptgw,nkpt,& 191 & dtset%nsym,sppoldbl,dtset%symafm,dtset%symrel,timrev,xmpi_comm_self) 192 193 if (dksqmax>tol8) then 194 write(msg,'(5a)')& 195 & 'Set of GW irreducible-zone kptgw in input file is inconsistent',ch10,& 196 & 'with full-zone set being used for wannier90 setup.',ch10,& 197 & 'Action: correct input data' 198 ABI_ERROR(msg) 199 end if 200 ! 201 ! === Initialize object defining the Band strucuture === 202 ! * Initialize with KS results using IBZ indexing. 203 ! * After rdqps, QP_bst will contain the QP amplitudes. 204 nkibz=Dtset%nkptgw 205 ABI_MALLOC(kibz,(3,nkibz)) 206 ABI_MALLOC(wtk_ibz,(nkibz)) 207 kibz=Dtset%kptgw(:,1:Dtset%nkptgw) 208 209 ! MG: This part is needed to get the IBZ weight that will be reported 210 ! on ab_out thus we should be consistent. Ideally Cryst should be 211 ! one of the basic abinit objects and it should be passed to this routine. 212 213 !different conventions are used in GW and abinit!! 214 cryst = hdr%get_crystal(gw_timrev=timrev+1) 215 call Kibz_mesh%init(Cryst,nkibz,kibz,Dtset%kptopt) 216 wtk_ibz=Kibz_mesh%wt 217 call cryst%free() 218 call Kibz_mesh%free() 219 220 ABI_MALLOC(ibz2bz,(nkibz,6)) 221 call listkk(dksqmax,gmet,ibz2bz,dtset%kpt,dtset%kptgw,nkpt,dtset%nkptgw,& 222 & dtset%nsym,sppoldbl,dtset%symafm,dtset%symrel,timrev,xmpi_comm_self) 223 224 ltest=ALL(ibz2bz(:,2)==1) 225 ABI_CHECK(ltest,'Not able to found irreducible points in the BZ set!') 226 227 if (dksqmax>tol8) then 228 write(msg,'(5a)')& 229 'Set of GW irreducible-zone kptgw in input file is inconsistent',ch10,& 230 'with full-zone set being used for wannier90 setup.',ch10,& 231 'Action: correct input data' 232 ABI_ERROR(msg) 233 end if 234 235 ABI_MALLOC(npwarr_ibz,(nkibz)) 236 ABI_MALLOC(istwfk_ibz,(nkibz)) 237 ABI_MALLOC(nband_ibz,(nkibz*nsppol)) 238 239 do isppol=1,nsppol 240 do ikibz=1,nkibz 241 ikbz=ibz2bz(ikibz+(sppoldbl-1)*(isppol-1)*nkibz,1) 242 npwarr_ibz(ikibz)= npwarr(ikbz) 243 istwfk_ibz(ikibz)=Dtset%istwfk(ikbz) 244 nband_ibz(ikibz+(isppol-1)*nkibz)=Dtset%nband(ikbz+(isppol-1)*nkpt) 245 end do 246 end do 247 248 bantot_ibz=SUM(nband_ibz) 249 ABI_MALLOC(doccde_ibz,(bantot_ibz)) 250 ABI_MALLOC(eigen_ibz,(bantot_ibz)) 251 ABI_MALLOC(occfact_ibz,(bantot_ibz)) 252 doccde_ibz(:)=zero ; eigen_ibz(:)=zero ; occfact_ibz(:)=zero 253 254 band_index=0 255 do isppol=1,nsppol 256 do ikibz=1,nkibz 257 ikbz=ibz2bz(ikibz+(sppoldbl-1)*(isppol-1)*nkibz,1) 258 nband_k=nband_ibz(ikibz+(isppol-1)*nkibz) 259 ii=SUM(Dtset%nband(1:ikbz+(isppol-1)*nkpt))-nband_k 260 eigen_ibz(band_index+1:band_index+nband_k)=eigen(ii+1:ii+nband_k) 261 band_index=band_index+nband_k 262 end do 263 end do 264 265 call ebands_init(bantot_ibz,QP_bst,Dtset%nelect,Dtset%ne_qFD,Dtset%nh_qFD,Dtset%ivalence,& 266 doccde_ibz,eigen_ibz,istwfk_ibz,kibz,nband_ibz,& 267 nkibz,npwarr_ibz,nsppol,Dtset%nspinor,Dtset%tphysel,Dtset%tsmear,Dtset%occopt,occfact_ibz,wtk_ibz,& 268 dtset%cellcharge(1),dtset%kptopt,dtset%kptrlatt_orig,dtset%nshiftk_orig,dtset%shiftk_orig,& 269 dtset%kptrlatt,dtset%nshiftk,dtset%shiftk) 270 271 ABI_FREE(kibz) 272 ABI_FREE(wtk_ibz) 273 ABI_FREE(ibz2bz) 274 ABI_FREE(npwarr_ibz) 275 ABI_FREE(istwfk_ibz) 276 ABI_FREE(nband_ibz) 277 ABI_FREE(doccde_ibz) 278 ABI_FREE(eigen_ibz) 279 ABI_FREE(occfact_ibz) 280 281 ! === Read in quasiparticle information === 282 ! * Initialize QP amplitudes with KS, QP_bst% presently contains KS energies. 283 ! * If file not found return, everything has been already initialized with KS values 284 ! Here qp_rhor is not needed thus dimrho=0 285 ABI_MALLOC(m_ks_to_qp,(mband,mband,dtset%nkptgw,nsppol)) 286 m_ks_to_qp=czero 287 do iband=1,mband 288 m_ks_to_qp(iband,iband,:,:)=cone 289 end do 290 291 ! Fake MPI_type for rdqps 292 call initmpi_seq(MPI_enreg_seq) 293 294 my_ngfft=Dtset%ngfft; if (Dtset%usepaw==1.and.ALL(Dtset%ngfftdg(1:3)/=0)) my_ngfft=Dtset%ngfftdg 295 nfftot=PRODUCT(my_ngfft(1:3)); dimrho=0 296 297 ! Change gw_fname to read a GW file instead of the QPS file. 298 ! TODO not so sure that wannier90 can handle G0W0 eigenvalues that are not ordered, though! 299 gw_fname = "g0w0" 300 g0w0_exists = .FALSE. 301 inquire(file=gw_fname,iostat=ios,exist=g0w0_exists) 302 if (ios/=0) then 303 ABI_ERROR('File g0w0 exists but iostat returns nonzero value!') 304 end if 305 306 if (.not.g0w0_exists) then ! read QPS file (default behavior). 307 input = from_QPS_FILE 308 ABI_MALLOC(prev_Pawrhoij,(Cryst%natom*Dtset%usepaw)) 309 ABI_MALLOC(qp_rhor,(nfftot,nspden*dimrho)) 310 311 call rdqps(QP_bst,Dtfil%fnameabi_qps,Dtset%usepaw,Dtset%nspden,dimrho,nscf,& 312 nfftot,my_ngfft,ucvol,Cryst,Pawtab,MPI_enreg_seq,nbsc,m_ks_to_qp,qp_rhor,prev_Pawrhoij) 313 314 ABI_FREE(qp_rhor) 315 ABI_FREE(prev_Pawrhoij) 316 317 else 318 ! Read GW file (m_ks_to_qp has been already set to 1, no extrapolation is performed) 319 ABI_WARNING(' READING GW CORRECTIONS FROM FILE g0w0 !') 320 input = from_GW_FILE 321 ABI_MALLOC(igwene,(QP_bst%mband,QP_bst%nkpt,QP_bst%nsppol)) 322 call rdgw(QP_bst,gw_fname,igwene,extrapolate=.FALSE.) 323 ABI_FREE(igwene) 324 end if 325 326 ! === Begin big loop over full-zone k points and spin (not implemented) === 327 ! * Wannier90 treats only a single spin, changes in wannier90 are needed 328 ABI_MALLOC(cg_k,(mpw,mband)) 329 ABI_MALLOC(cg_qpk,(mpw,mband)) 330 ABI_MALLOC(m_tmp,(mband,mband)) 331 332 band_index=0 ; icg=0 ; ikg=0 333 do isppol=1,nsppol 334 do ikpt=1,nkpt 335 336 irzkpt =indkk(ikpt+(sppoldbl-1)*(isppol-1)*nkpt,1) 337 itimrev=indkk(ikpt+(sppoldbl-1)*(isppol-1)*nkpt,6) 338 npw_k=npwarr(ikpt) 339 nband_k=dtset%nband(ikpt+(isppol-1)*nkpt) 340 341 if (nband_k/=mband) then 342 write(msg,'(a,i0,7a)')& 343 'Number of bands for k point ',ikpt,' is inconsistent with number',ch10,& 344 'specified for wannier90 calculation',ch10,& 345 'Action: correct input so all band numbers are equal for GW',ch10,& 346 'and wannier90 datasets.' 347 ABI_ERROR(msg) 348 end if 349 350 ! Load KS states for this kbz and spin 351 do iband=1,nband_k 352 icg_shift=npw_k*my_nspinor*(iband-1)+icg 353 do ipw=1,npw_k 354 cg_k(ipw,iband)=DCMPLX(cg(1,ipw+icg_shift),cg(2,ipw+icg_shift)) 355 end do 356 end do 357 358 ! If time reversal is used for relating ikpt to irzkpt, then multiply by 359 ! the complex conjugate of the lda-to-qp transformation matrix 360 if (itimrev==0) then 361 m_tmp(:,:)=m_ks_to_qp(:,:,irzkpt,isppol) 362 else if (itimrev==1) then 363 m_tmp(:,:)=conjg(m_ks_to_qp(:,:,irzkpt,isppol)) 364 else 365 write(msg,'(2(a,i0))')'Invalid indkk(ikpt,6) ',itimrev,'from routine listkk for k-point ',ikpt 366 ABI_BUG(msg) 367 end if 368 369 call ZGEMM('N','N',npw_k,mband,mband,cone,cg_k,mpw,m_tmp,mband,czero,cg_qpk,mpw) 370 371 ! === Orthonormality test === 372 ! * nband >= maxval(bndgw) for this to pass, but may be less than nband used in GW. 373 ! * Unfortunately, does not test WFK and QPS consistency. 374 !allocate(ortho(nband_k*(nband_k+1)/2)) 375 !ortho=czero; ijpack=0 376 !do jband=1,nband_k 377 ! jb_idx=band_index+jband 378 ! if (dtset%usepaw==1) Cp2 => Cprj_BZ(:,jband:jband+(my_nspinor-1)) 379 ! do iband=1,jband 380 ! ib_idx=band_index+iband 381 ! ijpack=ijpack+1 382 ! ortho(ijpack)=sum(conjg(cg_qpk(1:npw_k,iband))*cg_qpk(1:npw_k,jband)) 383 ! if (dtset%usepaw==1) then 384 ! Cp1 => Cprj_BZ(:,iband:iband+(my_nspinor-1)) 385 ! paw_ovlp = paw_overlap(Cp2,Cp1,Cryst%typat,Pawtab) 386 ! ortho(ijpack) = ortho(ijpack) + CMPLX(paw_ovlp(1),paw_ovlp(2)) 387 ! end if 388 ! if (jband==iband) ortho(ijpack)=ortho(ijpack)-cone 389 ! end do 390 !end do 391 !ortho_err=maxval(abs(ortho)) 392 393 !write(std_out,*)' drh - mlwfovlp_qp: ikpt,ortho_err',ikpt,ortho_err 394 !if (ortho_err>tol6) then 395 ! write(msg, '(3a,i4,a,i6,a,1p,e8.1,3a)' )& 396 !& ' orthonormality error for quasiparticle wave functions.',ch10,& 397 !& ' spin=',isppol,' k point=',ikpt,' ortho_err=',ortho_err,' >1E-6',ch10,& 398 !& ' Action: Be sure input nband>=maxval(bndgw)' 399 ! ABI_ERROR(msg) 400 !end if 401 !deallocate(ortho) 402 403 ! Replace lda wave functions and eigenvalues with quasiparticle ones. 404 qpenek_is_ordered = isordered(nband_k,QP_bst%eig(:,irzkpt,isppol),">",TOL_SORT) 405 406 if (input==from_QPS_FILE .and. .not.qpenek_is_ordered) then 407 write(msg,'(3a)')& 408 " QP energies read from QPS file are not ordered, likely nband_k>nbdgw. ",ch10,& 409 " Change nband in the input file so that it equals the number of GW states calculated" 410 ABI_WARNING(msg) 411 end if 412 413 if ( .TRUE. ) then 414 do iband=1,nband_k 415 icg_shift=npw_k*my_nspinor*(iband-1)+icg 416 eigen(iband+band_index)=QP_bst%eig(iband,irzkpt,isppol) 417 do ipw=1,npw_k 418 cg(1,ipw+icg_shift)= real(cg_qpk(ipw,iband)) 419 cg(2,ipw+icg_shift)=aimag(cg_qpk(ipw,iband)) 420 end do 421 end do 422 else 423 ! FIXME There's a problem in twannier90 since nband_k > nbdgw and therefore we also read KS states from the QPS file! 424 ! Automatic test has to be changed! 425 write(msg,'(2a,3f8.4,3a)')ch10,& 426 "QP energies at k-point ",QP_bst%kptns(:,irzkpt)," are not sorted in ascending numerical order!",ch10,& 427 "Performing reordering of energies and wavefunctions to be written on the final WKF file." 428 ABI_ERROR(msg) 429 !write(std_out,*)"eig",(QP_bst%eig(ii,irzkpt,isppol),ii=1,nband_k) 430 ABI_MALLOC(sorted_qpene,(nband_k)) 431 ABI_MALLOC(iord,(nband_k)) 432 sorted_qpene = QP_bst%eig(1:nband_k,irzkpt,isppol) 433 iord = (/(ii, ii=1,nband_k)/) 434 435 call sort_dp(nband_k,sorted_qpene,iord,TOL_SORT) 436 do ii=1,nband_k 437 write(std_out,*)"%eig, sorted_qpene, iord",QP_bst%eig(ii,irzkpt,isppol)*Ha_eV,sorted_qpene(ii)*Ha_eV,iord(ii) 438 end do 439 440 do iband=1,nband_k 441 ord_iband = iord(iband) 442 icg_shift=npw_k*my_nspinor*(iband-1)+icg 443 !eigen(iband+band_index)=QP_bst%eig(iband,irzkpt,isppol) 444 eigen(iband+band_index)=QP_bst%eig(ord_iband,irzkpt,isppol) 445 do ipw=1,npw_k 446 !cg(1,ipw+icg_shift)= real(cg_qpk(ipw,iband)) 447 !cg(2,ipw+icg_shift)=aimag(cg_qpk(ipw,iband)) 448 cg(1,ipw+icg_shift)= real(cg_qpk(ipw,ord_iband)) 449 cg(2,ipw+icg_shift)=aimag(cg_qpk(ipw,ord_iband)) 450 end do 451 end do 452 ABI_FREE(sorted_qpene) 453 ABI_FREE(iord) 454 end if 455 456 band_index=band_index+nband_k 457 icg=icg+npw_k*my_nspinor*nband_k 458 ikg=ikg+npw_k 459 end do !ikpt 460 end do !isppol 461 462 ABI_FREE(cg_k) 463 ABI_FREE(cg_qpk) 464 ABI_FREE(m_tmp) 465 466 ! === If PAW, update projections in BZ === 467 ! * Since I am lazy and here I do not care about memory, I just reconstruct m_ks_to_qp in the BZ. 468 ! * update_cprj will take care of updating the PAW projections to get <p_lmn|QP_{nks]> 469 ! This allows some CPU saving, no need to call ctocprj. 470 ! FIXME this part should be tested, automatic test to be provided 471 if (Dtset%usepaw==1) then 472 ABI_MALLOC(dimlmn,(natom)) 473 call pawcprj_getdim(dimlmn,dtset%natom,nattyp_dum,ntypat,Dtset%typat,pawtab,'R') 474 475 nkbz=nkpt 476 ABI_MALLOC(m_ks_to_qp_BZ,(mband,mband,nkbz,nsppol)) 477 do isppol=1,nsppol 478 do ikbz=1,nkbz 479 ikibz =indkk(ikibz+(sppoldbl-1)*(isppol-1)*nkbz,1) 480 itimrev=indkk(ikibz+(sppoldbl-1)*(isppol-1)*nkbz,6) 481 select case (itimrev) 482 case (0) 483 m_ks_to_qp_BZ(:,:,ikbz,isppol)=m_ks_to_qp(:,:,ikibz,isppol) 484 case (1) 485 m_ks_to_qp_BZ(:,:,ikbz,isppol)=CONJG(m_ks_to_qp(:,:,ikibz,isppol)) 486 case default 487 write(msg,'(a,i3)')"Wrong itimrev= ",itimrev 488 ABI_BUG(msg) 489 end select 490 end do 491 end do 492 493 call update_cprj(natom,nkbz,mband,nsppol,my_nspinor,m_ks_to_qp_BZ,dimlmn,Cprj_BZ) 494 ABI_FREE(dimlmn) 495 ABI_FREE(m_ks_to_qp_BZ) 496 end if !PAW 497 498 write(msg,'(6a)')ch10,& 499 ' mlwfovlp_qp: Input KS wavefuctions have been converted',ch10,& 500 ' to GW quasiparticle wavefunctions for maximally localized wannier',ch10,& 501 ' function construction by wannier90.' 502 call wrtout(ab_out,msg,'COLL') 503 call wrtout(std_out,msg,'COLL') 504 505 ABI_FREE(m_ks_to_qp) 506 call ebands_free(QP_bst) 507 call destroy_mpi_enreg(MPI_enreg_seq) 508 509 DBG_EXIT("COLL") 510 511 end subroutine mlwfovlp_qp