TABLE OF CONTENTS
- ABINIT/m_qparticles
- m_qparticles/rdgw
- m_qparticles/rdqps
- m_qparticles/show_QP
- m_qparticles/updt_m_ks_to_qp
- m_qparticles/wrqps
ABINIT/m_qparticles [ Modules ]
NAME
m_qparticles
FUNCTION
This module contains tools for the IO of the QP file and other procedures related to the calculation of the quasiparticle amplitudes represented in terms of KS states.
COPYRIGHT
Copyright (C) 2008-2024 ABINIT group (FB, MG) 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
18 #if defined HAVE_CONFIG_H 19 #include "config.h" 20 #endif 21 22 #include "abi_common.h" 23 24 MODULE m_qparticles 25 26 use defs_basis 27 use m_abicore 28 use m_hdr 29 use m_errors 30 use m_nctk 31 use m_distribfft 32 33 34 use defs_datatypes, only : pseudopotential_type, ebands_t 35 use defs_abitypes, only : MPI_type 36 use m_io_tools, only : open_file, file_exists, isncfile 37 use m_fstrings, only : int2char10, itoa, sjoin 38 use m_numeric_tools, only : linfit, c2r, set2unit, interpol3d_0d, rhophi 39 use m_gwdefs, only : sigparams_t 40 use m_crystal, only : crystal_t 41 use m_bz_mesh, only : kmesh_t 42 use m_ebands, only : ebands_get_valence_idx 43 use m_sigma, only : sigma_t 44 use m_pawtab, only : pawtab_type 45 use m_pawrhoij, only : pawrhoij_type, pawrhoij_alloc, pawrhoij_io, pawrhoij_inquire_dim 46 use m_fourier_interpol,only : fourier_interpol 47 48 implicit none 49 50 private 51 52 public :: wrqps ! Write a QPS file. 53 public :: rdqps ! Read a QPS file. 54 public :: show_QP ! Report the components of a QP amplitude in terms of KS eigenstates. 55 public :: rdgw ! Read GW corrections from an external file. 56 public :: updt_m_ks_to_qp ! Updates the matrix of unitary transformation from lda to qp states. 57 58 CONTAINS !=======================================================================================
m_qparticles/rdgw [ Functions ]
[ Top ] [ m_qparticles ] [ Functions ]
NAME
rdgw
FUNCTION
This subroutine reads the GW corrections from a _GW file.
INPUTS
[extrapolate]= if .TRUE., the routine extrapolates the
GW corrections for the states that have not been explicitly evaluated (default).
If .FALSE., only the GW states that have been calculated will be used to replace
the input eigenvalues stored in Bst%eig
Bst<ebands_t>=type describing the Band structure.
%nbnds=number of bands.
%nkpt=number of irred k-points.
%nsppol=number of spin
%kptns(3,nkpt)=irreducible k-points
SIDE EFFECTS
Bst%eig(%mband,%nkpt,%nsppol)=Overwritten with GW energies according to extrapolate flag.
OUTPUT
igwene(Bst%mband,Bst%nkpt,Bst%nsppol)= The imaginary part of the QP energies.
SOURCE
683 subroutine rdgw(Bst,fname,igwene,extrapolate) 684 685 !Arguments ------------------------------------ 686 !scalars 687 character(len=*),intent(in) :: fname 688 logical,optional,intent(in) :: extrapolate 689 type(ebands_t),intent(inout) :: Bst 690 !arrays 691 real(dp),intent(out) :: igwene(Bst%mband,Bst%nkpt,Bst%nsppol) 692 693 !Local variables ------------------------------ 694 !scalars 695 integer :: ib,ibr,ik,ikibz,ikr,is,nn,nbandR,nkibzR,nsppolR,unt,nbv 696 real(dp) :: alpha,beta,degw,egw_r,egw_i,smrt 697 logical :: do_extrapolate 698 character(len=500) :: msg 699 !arrays 700 integer,allocatable :: vbik(:,:),seen(:) 701 real(dp) :: kread(3) 702 real(dp),allocatable :: gwcorr(:,:,:) 703 704 !************************************************************************ 705 706 call wrtout(std_out,'Reading GW corrections from file: '//TRIM(fname)) 707 ABI_CHECK(ALL(Bst%nband==Bst%mband),"nband must be constant") 708 709 if (open_file(fname,msg,newunit=unt,status='old') /=0) then 710 ABI_ERROR(msg) 711 end if 712 713 read(unt,*)nkibzR,nsppolR 714 715 ABI_CHECK(nsppolR==Bst%nsppol,"mismatch in nsppol") 716 if (nkibzR/=Bst%nkpt) then 717 write(msg,'(a,i4,a,i4,2a)')& 718 & 'Found less k-points than that required ',nkibzR,'/',Bst%nkpt,ch10,& 719 & 'Some k-points will be skipped. Continuing anyway ' 720 ABI_WARNING(msg) 721 end if 722 723 ABI_MALLOC(gwcorr,(Bst%mband,Bst%nkpt,Bst%nsppol)) 724 ABI_MALLOC(seen,(Bst%nkpt)) 725 gwcorr=zero 726 igwene=zero 727 728 do is=1,Bst%nsppol 729 seen=0 730 731 do ikr=1,nkibzR 732 read(unt,*)kread(:) 733 read(unt,*)nbandR 734 ikibz=0 735 do ik=1,Bst%nkpt 736 if (ALL(ABS(kread(:)-Bst%kptns(:,ik))<0.0001)) then 737 ikibz=ik 738 seen(ik) = seen(ik) + 1 739 end if 740 end do 741 do ib=1,nbandR 742 read(unt,*)ibr,egw_r,degw,egw_i 743 if (ibr<=Bst%mband .and. ikibz/=0) then 744 gwcorr(ibr,ikibz,is)=degw/Ha_eV 745 igwene(ibr,ikibz,is)=egw_i/Ha_eV 746 end if 747 end do 748 end do 749 750 if (ANY(seen/=1)) then 751 do ik=1,Bst%nkpt 752 if (seen(ik)/=1) then 753 write(msg,'(a,3f8.3,a)')" k-point: ",Bst%kptns(:,ik)," not found in the GW file!" 754 ABI_WARNING(msg) 755 end if 756 end do 757 end if 758 759 end do 760 761 ABI_FREE(seen) 762 close(unt) 763 764 do_extrapolate=.TRUE.; if (PRESENT(extrapolate)) do_extrapolate=extrapolate 765 766 if (.not. do_extrapolate) then ! Only the bands calculated are updated. 767 Bst%eig = Bst%eig + gwcorr 768 769 else 770 771 if (ANY(ABS(igwene)>tol6)) then 772 write(msg,'(4a)')ch10,& 773 & "The GW file contains QP energies with non-zero imaginary part",ch10,& 774 & "Extrapolation not coded, change the source! " 775 ABI_ERROR(msg) 776 end if 777 778 ABI_MALLOC(vbik,(BSt%nkpt,BSt%nsppol)) 779 vbik(:,:) = ebands_get_valence_idx(BSt) 780 781 do is=1,Bst%nsppol 782 do ik=1,Bst%nkpt 783 784 nbv=vbik(ik,is) ! Index of the (valence band| Fermi band) for each spin 785 nn=Bst%mband-nbv 786 787 do ib=nbv+1,Bst%mband 788 if ( ABS(gwcorr(ib,ik,is)) < tol16) then 789 nn=ib-1-nbv 790 if (nn>1) then 791 call wrtout(std_out,"Linear extrapolating (conduction) GW corrections beyond the read values") 792 smrt=linfit(nn,Bst%eig(nbv+1:nbv+nn,ik,is),gwcorr(nbv+1:nbv+nn,ik,is),alpha,beta) 793 else 794 call wrtout(std_out,"Assuming constant (conduction) GW corrections beyond the read values") 795 alpha=zero 796 beta =gwcorr(nbv+nn,ik,is) 797 end if 798 EXIT !ib loop 799 end if 800 end do !ib 801 802 do ib=nbv+nn+1,Bst%mband 803 gwcorr(ib,ik,is)= alpha*Bst%eig(ib,ik,is) + beta 804 end do 805 806 nn=nbv 807 do ib=nbv,1,-1 808 if ( ABS(gwcorr(ib,ik,is)) < tol16) then 809 nn=nbv-ib 810 if (nn>1) then 811 call wrtout(std_out,"Linear extrapolating (valence) GW corrections beyond the read values") 812 smrt=linfit(nn,Bst%eig(nbv-nn+1:nbv,ik,is),gwcorr(nbv-nn+1:nbv,ik,is),alpha,beta) 813 else 814 call wrtout(std_out,"Assuming constant (valence) GW corrections beyond the read values") 815 alpha=zero 816 beta =gwcorr(nbv,ik,is) 817 end if 818 EXIT !ib 819 end if 820 end do !ib 821 822 do ib=1,nbv-nn 823 gwcorr(ib,ik,is)=alpha*Bst%eig(ib,ik,is) + beta 824 end do 825 826 end do !ik 827 end do !is 828 829 call wrtout(std_out,' k s GW corrections [eV] ') 830 do is=1,Bst%nsppol 831 do ik=1,Bst%nkpt 832 write(msg,'(i3,1x,i3,10f7.2/50(10x,10f7.2/))')ik,is,(Ha_eV*gwcorr(ib,ik,is),ib=1,Bst%mband) 833 call wrtout(std_out,msg) 834 end do 835 end do 836 Bst%eig = Bst%eig + gwcorr 837 ABI_FREE(vbik) 838 end if 839 840 call wrtout(std_out,' k s GW eigenvalues [eV]') 841 do is=1,Bst%nsppol 842 do ik=1,Bst%nkpt 843 write(std_out,'(2(i3,1x),7x,10f7.2/50(15x,10f7.2/))')ik,is,(Ha_eV*Bst%eig(ib,ik,is),ib=1,Bst%mband) 844 end do 845 end do 846 847 ABI_FREE(gwcorr) 848 849 end subroutine rdgw
m_qparticles/rdqps [ Functions ]
[ Top ] [ m_qparticles ] [ Functions ]
NAME
rdqps
FUNCTION
Read a _QPS file containing the QP energies of the previous iteration, the coefficients defining the QP amplitudes in terms of the KS basis set and the QP density for mixing.
INPUTS
nfftot=Total number of FFT points for density ngfftf(18)=Info on the FFT mesh for the density. nspden=Number of SPin-DENsity components. usepaw=1 if we are using PAW. fname=Name of the file dimrho=1 if density has to be read, 0 otherwise BSt<ebands_t>=Structure containing the initial band structure. ucvol=Volume of the unit cell
OUTPUT
nbsc=number of bands used to describe the QP amplitudes
nscf=number of iterations that have been performed (==0 if we start from a KS calculation)
m_ks_to_qp(mband,mband,nibz,nsppol)=matrix giving the decomposition of the QP
wavefunction in the mainfold generated by the KS wavefunctions
(i.e. $ m_ks_to_qp(ib,jb,k,s) := <\psi_{ib,k,s}^{KS}|\psi_{jb,k,s}^{QP}>$
rhor_out(nfftot,nspden)=quasiparticle density
SIDE EFFECTS
BSt<ebands_t>=Structure containing the initial band structure. %en_qp(mband,nkpt,nsppol)=QP energies at iteration nscf
TODO
The value of nspden is not reported in the QPS file thus we have a possible undetected error.
SOURCE
259 subroutine rdqps(BSt,fname,usepaw,nspden,dimrho,nscf,& 260 & nfftot,ngfftf,ucvol,Cryst,Pawtab,MPI_enreg,nbsc,m_ks_to_qp,rhor_out,Pawrhoij) 261 262 !Arguments ------------------------------------ 263 !scalars 264 integer,intent(in) :: nfftot,nspden,usepaw,dimrho 265 integer,intent(out) :: nbsc,nscf 266 real(dp),intent(in) :: ucvol 267 character(len=*),intent(in) :: fname 268 type(crystal_t),intent(in) :: Cryst 269 type(ebands_t),intent(inout) :: BSt 270 type(MPI_type),intent(inout) :: MPI_enreg 271 !arrays 272 integer,intent(in) :: ngfftf(18) 273 real(dp),intent(out) :: rhor_out(nfftot,nspden*dimrho) 274 complex(dpc),intent(out) :: m_ks_to_qp(BSt%mband,BSt%mband,BSt%nkpt,BSt%nsppol) 275 type(Pawtab_type),intent(in) :: Pawtab(Cryst%ntypat*usepaw) 276 type(Pawrhoij_type),intent(inout) :: Pawrhoij(Cryst%natom*usepaw) 277 278 !Local variables------------------------------- 279 !scalars 280 integer,parameter :: master=0 281 integer :: ib,ii,ik,isppol,nbandR,nkibzR,nsppolR,unqps,my_rank,ispden 282 integer :: ifft,n1,n2,n3,ir1,ir2,ir3,ios 283 integer :: cplex_fft,optin,optout,nfft_found 284 integer :: iatom,natomR,nspdenR,ntypatR,itypat 285 real(dp) :: uerr,nelect_qps,ratio 286 logical,parameter :: use_FFT_interpolation=.TRUE. 287 logical :: ltest 288 character(len=500) :: msg 289 !arrays 290 integer :: ngfft_found(18) 291 integer,allocatable :: nlmn_type(:),typatR(:) 292 real(dp) :: kibz(3),rr(3),rhogdum(1,1) 293 real(dp),allocatable :: en_tmp(:) 294 real(dp),allocatable :: rhor_tmp(:,:) 295 complex(dpc),allocatable :: mtmp(:,:),utest(:,:) 296 297 ! ************************************************************************* 298 299 DBG_ENTER("COLL") 300 301 ABI_CHECK(ALL(BSt%nband==BSt%nband(1)),"No. of bands must be constant") 302 ABI_CHECK(dimrho==0.or.dimrho==1,'dimrho must be 0 or 1') 303 304 ! This does not work in parallel !!? 305 !% my_rank = xmpi_comm_rank(MPI_enreg%spaceComm) 306 my_rank = MPI_enreg%me_kpt 307 308 ! Check whether file exists or not. 309 write(msg,'(5a)')ch10,& 310 ' rdqps: reading QP wavefunctions of the previous step ',ch10,& 311 ' looking for file ',TRIM(fname) 312 call wrtout([std_out, ab_out], msg) 313 314 if (.not.file_exists(fname)) then 315 write(msg,'(2a)')' file not found, 1st iteration initialized with KS eigenelements ',ch10 316 call wrtout([std_out, ab_out], msg) 317 nscf=0; RETURN 318 end if 319 320 if (.not.isncfile(fname)) then 321 if (open_file(fname,msg,newunit=unqps,form='formatted',status='unknown') /= 0) then 322 ABI_ERROR(msg) 323 end if 324 325 ! TODO the _QPS file should contain additional information 326 read(unqps,*)nscf 327 write(msg,'(a,i4,a)')' Number of iteration(s) already performed: ',nscf,ch10 328 call wrtout([std_out, ab_out], msg) 329 330 read(unqps,*)nkibzR 331 if (nkibzR/=BSt%nkpt) then 332 write(msg,'(2(a,i0))')'Wrong number of k-points; Expected: ',BSt%nkpt,', Found: ',nkibzR 333 ABI_ERROR(msg) 334 end if 335 336 read(unqps,*)nbandR 337 nbsc=MIN(nbandR,BSt%mband) 338 339 if (nbsc/=BSt%mband) then 340 write(msg,'(3a,i4,a,i4)')& 341 'QPS file contains less bands than that used in the present calculation ',ch10,& 342 'Required: ',BSt%mband,', Found: ',nbandR 343 ABI_WARNING(msg) 344 end if 345 346 if (nbsc/=nbandR) then 347 write(msg,'(3a,i4,a)')& 348 'The QPS file contains more bands than that used in the present calculation ',ch10,& 349 'only the first ',nbandR,' bands will be read' 350 ABI_COMMENT(msg) 351 end if 352 353 ABI_MALLOC(mtmp,(nbandR,nbandR)) 354 ABI_MALLOC(en_tmp,(nbandR)) 355 read(unqps,*)nsppolR 356 357 ABI_CHECK(nsppolR==BSt%nsppol,"QPS generated with different nsppol") 358 359 ! Read energies and transformation for each k-point and spin. 360 ! TODO: The format of the QPS file must be standardized ! 361 ! For example we might add the occupation numbers. 362 do isppol=1,BSt%nsppol 363 do ik=1,BSt%nkpt 364 read(unqps,*)kibz(:) 365 write(msg,'(a,i5,a,3(f6.3,1x),4x,a,i2)')' Reading ik ',ik,') k = ',kibz(:),' is = ',isppol 366 call wrtout(std_out,msg) 367 ltest=(ALL(ABS(kibz(:)-BSt%kptns(:,ik))<0.001)) 368 ABI_CHECK(ltest,'Wrong k-point read') 369 do ib=1,nbandR 370 read(unqps,*)en_tmp(ib) 371 read(unqps,*)mtmp(:,ib) 372 end do 373 374 ! Store transformation and update energies. 375 m_ks_to_qp(1:nbsc,1:nbsc,ik,isppol)=mtmp(1:nbsc,1:nbsc) 376 BSt%eig(1:nbsc,ik,isppol)=en_tmp(1:nbsc) 377 378 ! Chech if matrix is unitary. 379 ABI_MALLOC(utest,(nbsc,nbsc)) 380 utest(:,:) = TRANSPOSE(mtmp(1:nbsc,1:nbsc)) !this is just for the buggy gfortran 381 utest(:,:) = MATMUL(CONJG(utest),mtmp(1:nbsc,1:nbsc)) 382 do ii=1,nbsc 383 utest(ii,ii)=utest(ii,ii)-one 384 end do 385 uerr=MAXVAL(ABS(utest)) 386 if (uerr>tol6) then 387 write(msg,'(a,es16.8)')' KS -> QP matrix is not unitary, MAX error = ',uerr 388 ABI_WARNING(msg) 389 end if 390 ABI_FREE(utest) 391 end do !ik 392 end do !isppol 393 394 ABI_FREE(mtmp) 395 ABI_FREE(en_tmp) 396 397 ! Read the QP density. 398 ! The two FFT grids might differ. In case perform an FFT interpolation to have rhor on the input mesh. 399 if (dimrho==1) then 400 read(unqps,*)n1,n2,n3 401 402 if (all(ngfftf(1:3)== [n1, n2, n3]) ) then 403 read(unqps,*)rhor_out(:,:) 404 else 405 write(msg,'(2a,a,5(i3,a),i3)')& 406 'FFT meshes differ. Performing Fourier interpolation. ',ch10,& 407 'Found: ',n1,' x',n2,' x',n3,'; Expected: ',ngfftf(1),' x',ngfftf(2),' x',ngfftf(3) 408 ABI_COMMENT(msg) 409 410 ABI_MALLOC(rhor_tmp,(n1*n2*n3,nspden)) 411 read(unqps,*)rhor_tmp(:,:) 412 413 if (use_FFT_interpolation) then 414 ngfft_found(1:3)=(/n1,n2,n3/) 415 ngfft_found(4)=2*(ngfft_found(1)/2)+1 ! 4:18 are not used, anyway! 416 ngfft_found(5)=2*(ngfft_found(2)/2)+1 417 ngfft_found(6)=ngfft_found(3) 418 ngfft_found(7:18)=ngfftf(7:18) 419 nfft_found=PRODUCT(ngfft_found(1:3)) !no FFT para 420 421 cplex_fft =1 ! Real quantities. 422 optin =0 ! Input is taken from rhor. 423 optout =0 ! Output is only in real space. 424 call destroy_distribfft(MPI_enreg%distribfft) 425 call init_distribfft(MPI_enreg%distribfft,'c',MPI_enreg%nproc_fft,ngfftf(2),ngfftf(3)) 426 call init_distribfft(MPI_enreg%distribfft,'f',MPI_enreg%nproc_fft,ngfft_found(2),ngfft_found(3)) 427 428 call fourier_interpol(cplex_fft,nspden,optin,optout,nfft_found,ngfft_found,nfftot,ngfftf,& 429 MPI_enreg,rhor_tmp,rhor_out,rhogdum,rhogdum) 430 431 else 432 ! Linear interpolation. 433 do ispden=1,nspden 434 do ir3=0,ngfftf(3)-1 435 rr(3)=DBLE(ir3)/n3 436 do ir2=0,ngfftf(2)-1 437 rr(2)=DBLE(ir2)/n2 438 do ir1=0,ngfftf(1)-1 439 rr(1)=DBLE(ir1)/n1 440 ifft = 1 +ir1 +ir2*ngfftf(1) +ir3*ngfftf(1)*ngfftf(2) 441 rhor_out(ifft,ispden) = interpol3d_0d(rr,n1,n2,n3,rhor_tmp(:,ispden)) 442 end do 443 end do 444 end do 445 end do 446 end if 447 448 ABI_FREE(rhor_tmp) 449 end if 450 451 ! Test the normalization of the QPS density. 452 ! There might be errors due to the interpolation or the truncation of the G basis set 453 ! Density will be renormalized in the caller since for PAW we still have to add the onsite contribution. 454 if (usepaw==0) then 455 nelect_qps=SUM(rhor_out(:,1))*ucvol/nfftot; ratio=BSt%nelect/nelect_qps 456 write(msg,'(3(a,f9.4))')& 457 ' Number of electrons calculated using the QPS density = ',nelect_qps,' Expected = ',BSt%nelect,' ratio = ',ratio 458 call wrtout(std_out, msg) 459 !!rhor_out(:,:)=ratio*rhor_out(:,:) 460 end if 461 462 if (usepaw==1) then 463 ! Write QP_rhoij for on-site density mixing. 464 read(unqps,*,iostat=ios)natomR,ntypatR 465 if (ios/=0) then 466 msg="Old version of QPS file found. DO NOT USE rhoqpmix for this run." 467 ABI_WARNING(msg) 468 call wrtout(ab_out,msg) 469 ! Init dummy rhoij just to avoid problems in sigma when rhoij is freed. 470 call pawrhoij_inquire_dim(nspden_rhoij=nspdenR, nspden=nspden) 471 call pawrhoij_alloc(Pawrhoij,1,nspdenR,BSt%nspinor,BSt%nsppol,Cryst%typat,pawtab=Pawtab) 472 close(unqps) 473 RETURN 474 end if 475 476 ABI_CHECK(natomR ==Cryst%natom, "mismatch in natom") 477 ABI_CHECK(ntypatR==Cryst%ntypat,"mismatch in ntypat") 478 ABI_MALLOC(nlmn_type,(ntypatR)) 479 ABI_MALLOC(typatR,(ntypatR)) 480 481 read(unqps,*)(typatR(iatom), iatom=1,natomR) 482 ABI_CHECK(ALL(Cryst%typat==typatR),"mismatch in typat") 483 484 read(unqps,*)(nlmn_type(itypat), itypat=1,ntypatR) 485 do itypat =1,Cryst%ntypat 486 if (nlmn_type(itypat)/=Pawtab(itypat)%lmn_size) then 487 ABI_ERROR("mismatch in nlmn_type, check QPS file") 488 end if 489 end do 490 491 read(unqps,*) nsppolR,nspdenR 492 ABI_CHECK(nsppolR==BSt%nsppol,"mismatch in nsppol") 493 ABI_CHECK(nspdenR==nspden ,"mismatch in nspden") 494 495 call pawrhoij_io(pawrhoij,unqps,BSt%nsppol,BSt%nspinor,nspden,nlmn_type,Cryst%typat,& 496 HDR_LATEST_HEADFORM,"Read",form="formatted") 497 !% call pawrhoij_io(pawrhoij,std_out,BSt%nsppol,BSt%nspinor,nspden,nlmn_type,Cryst%typat,HDR_LATEST_HEADFORM,"Echo") 498 499 ABI_FREE(nlmn_type) 500 ABI_FREE(typatR) 501 end if ! usepaw 502 503 end if !dimrho=1 504 505 close(unqps) 506 507 else 508 ABI_ERROR("netdf format not implemented") 509 end if 510 511 DBG_EXIT("COLL") 512 513 end subroutine rdqps
m_qparticles/show_QP [ Functions ]
[ Top ] [ m_qparticles ] [ Functions ]
NAME
show_QP
FUNCTION
Print in a nice format (?) the expansion coefficients of the quasiparticle amplitudes in terms of the KS eigenvectors
INPUTS
Bst<ebands_t>=Description of the band structure. %nsppol=1 for unpolarized, 2 for spin-polarized. %mband=Max number of bands (in GW doesn"t depend on k an spin) %nkpt=number of irreducible k-points. %eig(mband,nkpt,nsppol)= QP energies for each k-point, band and spin. m_ks_to_qp(nbnds,nbnds,nkibz,nsppol)=matrix giving the decomposition of the QP amplitued in the mainfold generated by the KS wavefunctions (i.e $ m_ks_to_qp(ib,jb,k,s) := \langle \psi_{ib,k,s}^{KS}| \psi_{jb,k,s}^{QP}\rangle $ fromb,tob=initial and final band index for QP, only states in this range are printed prtvol=Verbosity level (not used) unit=Unit number of the output file tolmat[Optional]=Only components whose coefficient has modulus larger than tolmat are shown (default is 0.01)
OUTPUT
Only printing
NOTES
Only master node should call this routine.
SOURCE
548 subroutine show_QP(Bst,m_ks_to_qp,fromb,tob,unit,prtvol,tolmat,kmask) 549 550 !Arguments ------------------------------------ 551 !scalars 552 integer,optional,intent(in) :: fromb,tob 553 integer,optional,intent(in) :: prtvol,unit 554 real(dp),optional,intent(in) :: tolmat 555 type(ebands_t),intent(in) :: Bst 556 !arrays 557 logical,optional,intent(in) :: kmask(Bst%nkpt) 558 complex(dpc),intent(in) :: m_ks_to_qp(Bst%mband,Bst%mband,Bst%nkpt,Bst%nsppol) 559 560 !Local variables------------------------------- 561 !scalars 562 integer,parameter :: NBRA=5 563 logical,parameter :: use_rhophi=.True. 564 integer :: ib_start,ib_stop,my_prtvol,counter,ib_KS,ib_QP,ikibz,isp,nspace,my_unt,nband_k 565 real(dp) :: my_tolmat,rho,phi 566 character(len=10) :: bks,bqp,k_tag,spin_tag 567 character(len=500) :: KS_row,KS_ket,tmpstr,QP_ket 568 !arrays 569 real(dp) :: cx(2) 570 571 ! ********************************************************************* 572 573 my_unt =std_out ; if (PRESENT(unit )) my_unt =unit 574 my_prtvol=0 ; if (PRESENT(prtvol)) my_prtvol=prtvol 575 ib_start =1 ; if (PRESENT(fromb )) ib_start =fromb 576 ib_stop =Bst%mband; if (PRESENT(tob )) ib_stop =tob 577 my_tolmat=0.001 ; if (PRESENT(tolmat)) my_tolmat=ABS(tolmat) 578 579 ! I suppose nband_k is constant thus the check is done here. 580 if (ib_start<=0 ) ib_start=1 581 if (ib_start>Bst%mband) ib_start=Bst%mband 582 if (ib_stop<=0 ) ib_stop=1 583 if (ib_stop>Bst%mband ) ib_stop=Bst%mband 584 585 ! Have to follow rules 7.f. 586 write(my_unt,'(/,a,/,a,/,a,f6.3,a,/,a)')& 587 ' '//REPEAT('*',76),& 588 & ' ***** QP amplitudes expressed as linear combination of KS eigenstates. *****',& 589 & ' ***** Only KS components whose modulus is larger than ',my_tolmat,' are shown ***** ',& 590 & ' '//REPEAT('*',76) 591 if (use_rhophi) then 592 write(my_unt,"(a)")"Complex coefficients given in (rho, phi) polar representation." 593 else 594 write(my_unt,"(a)")"Complex coefficients given in (Re, Im) representation." 595 end if 596 597 if (PRESENT(kmask)) then 598 if (.not.ALL(kmask)) write(my_unt,'(/,a,i3,a)')' Only ',COUNT(kmask),' k-points are reported ' 599 end if 600 601 do isp=1,Bst%nsppol 602 call int2char10(isp,spin_tag) 603 write(my_unt,'(/,a,i2,a,/)')' >>>>> Begin block for spin ',isp,' <<<<< ' 604 605 do ikibz=1,Bst%nkpt 606 if (PRESENT(kmask)) then 607 if (.not.kmask(ikibz)) CYCLE 608 end if 609 call int2char10(ikibz,k_tag) 610 nband_k=Bst%nband(ikibz+(isp-1)*Bst%nkpt) 611 write(my_unt,'(a,i4,a,3es16.8,a,f6.3,/)')' k-point: ',ikibz,') ',Bst%kptns(:,ikibz),'; wtk= ',Bst%wtk(ikibz) 612 613 do ib_QP=ib_start,ib_stop 614 call int2char10(ib_QP,bqp) 615 QP_ket=' |QP: b='//TRIM(bqp)//'; s='//TRIM(spin_tag)//'> = ' 616 write(my_unt,'(a)')TRIM(QP_ket) 617 nspace=LEN(TRIM(QP_ket)) 618 619 counter=0 ; KS_row=REPEAT('',nspace+2) 620 do ib_KS=1,Bst%mband 621 if (ABS(m_ks_to_qp(ib_KS,ib_QP,ikibz,isp))<my_tolmat) CYCLE 622 counter=counter+1 623 call int2char10(ib_KS,bks) 624 write(tmpstr,'(3a)')' |',TRIM(bks),'>' 625 626 if (use_rhophi) then 627 ! coefficient as (rho, phi) 628 cx(1) = real(m_ks_to_qp(ib_KS,ib_QP,ikibz,isp)) 629 cx(2) = aimag(m_ks_to_qp(ib_KS,ib_QP,ikibz,isp)) 630 call rhophi(cx, phi, rho) 631 write(KS_ket,'(1x,2f7.3,a,1x)')rho, phi, TRIM(tmpstr) 632 else 633 ! coefficient as (Re, Im) 634 write(KS_ket,'(1x,2f7.3,a,1x)')m_ks_to_qp(ib_KS,ib_QP,ikibz,isp),TRIM(tmpstr) 635 end if 636 KS_row=TRIM(KS_row)//TRIM(KS_ket) 637 if (MOD(counter,NBRA)==0) then ! nbra KS kets per row 638 write(my_unt,'(a)')TRIM(KS_row) 639 KS_row=REPEAT('',nspace+2) 640 end if 641 end do 642 643 if (MOD(counter,NBRA)/=0) write(my_unt,'(a)')TRIM(KS_row) ! Last row, if any 644 write(my_unt,'(a)')'' 645 end do !ib_QP 646 647 end do !ikibz 648 end do !isp 649 650 write(my_unt,'(a,/)')' '//REPEAT('*',76) 651 652 end subroutine show_QP
m_qparticles/updt_m_ks_to_qp [ Functions ]
[ Top ] [ m_qparticles ] [ Functions ]
NAME
updt_m_ks_to_qp
FUNCTION
Updates the matrix containing the unitary transformation from the lda states to the quasiparticle states.
INPUTS
Sigp<sigparams_t>=Parameters characterizing the self-energy calculation. %nsppol=1 for unpolarized, 2 for spin-polarized %nbnds=number of bands used for sigma Sr<sigma_t>=Structure containing the results of the sigma run. %en_qp_diago(nbnds,nibz,nsppol)= NEW quasi-particle energies %eigvec_qp(nbnds,nbnds,nibz,nsppol)= NEW QP amplitudes in the KS basis set obtained by diagonalizing H0 + Herm(Sigma). Kmesh<kmesh_t>=information on the k-point sampling. %nibz=number of irreducible k-points %ibz(3,kibz)=reduced coordinates of the irreducible k-points nscf=Number of self consistent cycles performed
OUTPUT
(see side effects)
SIDE EFFECTS
m_ks_to_qp(nbnds,nbnds,nibz,nsppol)= overwritten with the new QP amplitudes
in terms of KS wavefunctions
NOTES
Only master node should call this routine.
SOURCE
887 subroutine updt_m_ks_to_qp(Sigp,Kmesh,nscf,Sr,m_ks_to_qp) 888 889 !Arguments ------------------------------------ 890 !scalars 891 integer,intent(in) :: nscf 892 type(kmesh_t),intent(in) :: Kmesh 893 type(sigparams_t),intent(in) :: Sigp 894 type(sigma_t),intent(in) :: Sr 895 !arrays 896 complex(dpc),intent(inout) :: m_ks_to_qp(Sigp%nbnds,Sigp%nbnds,Kmesh%nibz,Sigp%nsppol) 897 898 !Local variables------------------------------- 899 !scalars 900 integer :: ik,is 901 !arrays 902 complex(dpc),allocatable :: mtmp(:,:) 903 904 ! ************************************************************************* 905 906 if (nscf>=0) then 907 ! Calculate the new m_ks_to_qp 908 ABI_MALLOC(mtmp,(Sigp%nbnds,Sigp%nbnds)) 909 do is=1,Sigp%nsppol 910 do ik=1,Kmesh%nibz 911 mtmp(:,:)=m_ks_to_qp(:,:,ik,is) 912 m_ks_to_qp(:,:,ik,is)=MATMUL(mtmp(:,:),Sr%eigvec_qp(:,:,ik,is)) 913 end do 914 end do 915 ABI_FREE(mtmp) 916 end if 917 918 end subroutine updt_m_ks_to_qp 919 920 !---------------------------------------------------------------------- 921 922 END MODULE m_qparticles
m_qparticles/wrqps [ Functions ]
[ Top ] [ m_qparticles ] [ Functions ]
NAME
wrqps
FUNCTION
Write the _QPS file containing information on the quasi-particles energies and wavefunctions.
INPUTS
fname=The name of the file Sigp<sigparams_t>=Parameters characterizing the self-energy calculation. %nsppol=1 for unpolarized, 2 for spin-polarized %nbnds=number of bands used for sigma Sr<sigma_t>=Structure containing the results of the sigma run. %en_qp_diago(nbnds,nibz,nsppol)= NEW quasi-particle energies %eigvec_qp(nbnds,nbnds,nibz,nsppol)= NEW QP amplitudes in the KS basis set obtained by diagonalizing H0 + Herm(Sigma). m_ks_to_qp(nbnds,nbnds,nibz,nsppol)= expansion of the OLD QP amplitudes in terms of KS wavefunctions Kmesh<kmesh_t>=information on the k-point sampling. %nibz=number of irreducible k-points %ibz(3,kibz)=reduced coordinates of the irreducible k-points nfftot=Total number of FFT points for density ngfftf(18)=Info on the FFT mesh for the density. nscf=Number of self consistent cycles performed nspden=number of spin-density components Cryst<crystal_t>=Structure defining the crystal structure. Psps<type(pseudopotential_type)>=variables related to pseudopotentials. Pawrhoij(Cryst%natom*Psps%usepaw)<type(pawrhoij_type)>= rhoij datastructure. BSt<ebands_t>=Structure containing the band structure energies (only used is nscf==-1)
OUTPUT
Only writing
NOTES
Old QPS fileformat:
|
| No. of QPSCF cycles already performed.
| No. of k-points in the IBZ.
| Total number of bands used to construct the Green's function (nbnds)
| nsppol
| For each spin and k-point in the IBZ:
| Reduced coordinates of the k-point.
| for each band:
| QP energies obtained by diagonalizing the QPSCGW Hamiltonian.
| <\psi_{ib,k,s}^{KS}|\psi_{jb,k,s}^{QP}>$, ib=1,nbnds
| FFT dimensions of the fine grid
| QP density in real space.
SOURCE
110 subroutine wrqps(fname,Sigp,Cryst,Kmesh,Psps,Pawtab,Pawrhoij,nspden,nscf,nfftot,ngfftf,Sr,Bst,m_ks_to_qp,rho_qp) 111 112 !Arguments ------------------------------------ 113 !scalars 114 integer,intent(in) :: nfftot,nscf,nspden 115 character(len=*),intent(in) :: fname 116 type(kmesh_t),intent(in) :: Kmesh 117 type(ebands_t),intent(in) :: BSt 118 type(sigparams_t),intent(in) :: Sigp 119 type(sigma_t),intent(in) :: Sr 120 type(crystal_t),intent(in) :: Cryst 121 type(Pseudopotential_type),intent(in) :: Psps 122 !arrays 123 integer,intent(in) :: ngfftf(18) 124 real(dp),intent(in) :: rho_qp(nfftot,nspden) 125 complex(dpc),intent(in) :: m_ks_to_qp(Sigp%nbnds,Sigp%nbnds,Kmesh%nibz,Sigp%nsppol) 126 type(Pawrhoij_type),intent(inout) :: Pawrhoij(Cryst%natom*Psps%usepaw) 127 type(Pawtab_type),intent(in) :: Pawtab(Psps%ntypat*Psps%usepaw) 128 129 !Local variables------------------------------- 130 !scalars 131 integer :: ib,ik,is,unqps,iatom,itypat 132 character(len=500) :: msg 133 !arrays 134 integer,allocatable :: nlmn_type(:) 135 complex(dpc),allocatable :: mtmp(:,:) 136 137 ! ************************************************************************* 138 139 DBG_ENTER("COLL") 140 141 if (nscf >= 0) then 142 write(msg,'(3a)')ch10,' writing QP data on file : ',TRIM(fname) 143 call wrtout([std_out, ab_out], msg) 144 end if 145 146 if (open_file(fname,msg,newunit=unqps,form='formatted',status='unknown') /= 0) then 147 ABI_ERROR(msg) 148 end if 149 150 write(unqps,*)nscf+1 151 write(unqps,*)Kmesh%nibz 152 write(unqps,*)Sigp%nbnds 153 write(unqps,*)Sigp%nsppol 154 155 ABI_MALLOC(mtmp,(Sigp%nbnds,Sigp%nbnds)) 156 157 if (nscf>=0) then 158 ! Write the new m_ks_to_qp on file. 159 do is=1,Sigp%nsppol 160 do ik=1,Kmesh%nibz 161 write(unqps,*)Kmesh%ibz(:,ik) 162 do ib=1,Sigp%nbnds 163 write(unqps,*)Sr%en_qp_diago(ib,ik,is) 164 write(unqps,*)m_ks_to_qp(:,ib,ik,is) 165 end do 166 end do 167 end do 168 else if (nscf==-1) then 169 ! Write fake QPS file with KS band structure (Mainly used for G0W) 170 call set2unit(mtmp) 171 do is=1,Sigp%nsppol 172 do ik=1,Kmesh%nibz 173 write(unqps,*)Kmesh%ibz(:,ik) 174 do ib=1,Sigp%nbnds 175 write(unqps,*)BSt%eig(ib,ik,is) 176 write(unqps,*)mtmp(:,ib) 177 end do 178 end do 179 end do 180 else 181 ABI_ERROR(sjoin("Wrong nscf ",itoa(nscf))) 182 end if 183 184 ABI_FREE(mtmp) 185 186 write(msg,'(a,f9.4)')' (wrqps) planewave contribution to nelect: ',SUM(rho_qp(:,1))*Cryst%ucvol/nfftot 187 call wrtout(std_out,msg) 188 if (nspden == 4) then 189 write(msg,'(a,3f9.4)')' mx, my, mz: ',& 190 SUM(rho_qp(:,2))*Cryst%ucvol/nfftot,SUM(rho_qp(:,3))*Cryst%ucvol/nfftot,SUM(rho_qp(:,4))*Cryst%ucvol/nfftot 191 call wrtout(std_out,msg) 192 end if 193 194 ! Write FFT dimensions and QP density 195 write(unqps,*)ngfftf(1:3) 196 write(unqps,*)rho_qp(:,:) 197 198 if (Psps%usepaw==1) then 199 ! Write QP rhoij to be used for on-site density mixing. 200 ABI_MALLOC(nlmn_type,(Cryst%ntypat)) 201 do itypat=1,Cryst%ntypat 202 nlmn_type(itypat)=Pawtab(itypat)%lmn_size 203 end do 204 205 write(unqps,*) Cryst%natom, Cryst%ntypat 206 write(unqps,*) (Cryst%typat(iatom), iatom=1,Cryst%natom) 207 write(unqps,*) (nlmn_type(itypat), itypat=1,Cryst%ntypat) 208 write(unqps,*) Pawrhoij(1)%nsppol, Pawrhoij(1)%nspden 209 210 call pawrhoij_io(pawrhoij,unqps,Sigp%nsppol,Sigp%nspinor,nspden,nlmn_type,Cryst%typat,& 211 & HDR_LATEST_HEADFORM,"Write",form="formatted") 212 ABI_FREE(nlmn_type) 213 end if 214 215 close(unqps) 216 217 DBG_EXIT("COLL") 218 219 end subroutine wrqps