TABLE OF CONTENTS


ABINIT/m_qparticles [ Modules ]

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

PARENTS

CHILDREN

SOURCE

23 #if defined HAVE_CONFIG_H
24 #include "config.h"
25 #endif
26 
27 #include "abi_common.h"
28 
29 MODULE m_qparticles
30 
31  use defs_basis
32  use defs_datatypes
33  use defs_abitypes
34  use m_abicore
35  use m_hdr
36  use m_errors
37  use m_nctk
38 
39  use m_io_tools,       only : open_file, file_exists, isncfile
40  use m_fstrings,       only : int2char10, itoa, sjoin
41  use m_numeric_tools,  only : linfit, c2r, set2unit, interpol3d, rhophi
42  use m_gwdefs,         only : sigparams_t
43  use m_crystal,        only : crystal_t
44  use m_crystal_io,     only : crystal_ncwrite
45  use m_bz_mesh,        only : kmesh_t
46  use m_ebands,         only : get_valence_idx
47  use m_sigma,          only : sigma_t
48  use m_pawtab,         only : pawtab_type
49  use m_pawrhoij,       only : pawrhoij_type, pawrhoij_alloc, pawrhoij_io
50  use m_fourier_interpol,only : fourier_interpol
51 
52  implicit none
53 
54  private
55 
56  public :: wrqps             ! Write a QPS file.
57  public :: rdqps             ! Read a QPS file.
58  public :: show_QP           ! Report the components of a QP amplitude in terms of KS eigenstates.
59  public :: rdgw              ! Read GW corrections from an external file.
60  public :: updt_m_lda_to_qp  ! Updates the matrix of unitary transformation from lda to qp states.
61 
62 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.

PARENTS

      mlwfovlp_qp,screening,setup_bse,sigma

CHILDREN

SOURCE

737 subroutine rdgw(Bst,fname,igwene,extrapolate)
738 
739 
740 !This section has been created automatically by the script Abilint (TD).
741 !Do not modify the following lines by hand.
742 #undef ABI_FUNC
743 #define ABI_FUNC 'rdgw'
744 !End of the abilint section
745 
746  implicit none
747 
748 !Arguments ------------------------------------
749 !scalars
750  character(len=*),intent(in) :: fname
751  logical,optional,intent(in) :: extrapolate
752  type(ebands_t),intent(inout) :: Bst
753 !arrays
754  real(dp),intent(out) :: igwene(Bst%mband,Bst%nkpt,Bst%nsppol)
755 
756 !Local variables ------------------------------
757 !scalars
758  integer :: ib,ibr,ik,ikibz,ikr,is,nn,nbandR,nkibzR,nsppolR,unt,nbv
759  real(dp) :: alpha,beta,degw,egw_r,egw_i,smrt
760  logical :: do_extrapolate
761  character(len=500) :: msg
762 !arrays
763  integer,allocatable :: vbik(:,:),seen(:)
764  real(dp) :: kread(3)
765  real(dp),allocatable :: gwcorr(:,:,:)
766 
767 !************************************************************************
768 
769  call wrtout(std_out,'Reading GW corrections from file: '//TRIM(fname),'COLL')
770  ABI_CHECK(ALL(Bst%nband==Bst%mband),"nband must be constant")
771 
772  if (open_file(fname,msg,newunit=unt,status='old') /=0) then
773    MSG_ERROR(msg)
774  end if
775 
776  read(unt,*)nkibzR,nsppolR
777 
778  ABI_CHECK(nsppolR==Bst%nsppol,"mismatch in nsppol")
779  if (nkibzR/=Bst%nkpt) then
780    write(msg,'(a,i4,a,i4,2a)')&
781 &   'Found less k-points than that required ',nkibzR,'/',Bst%nkpt,ch10,&
782 &   'Some k-points will be skipped. Continuing anyway '
783    MSG_WARNING(msg)
784  end if
785 
786  ABI_MALLOC(gwcorr,(Bst%mband,Bst%nkpt,Bst%nsppol))
787  ABI_MALLOC(seen,(Bst%nkpt))
788  gwcorr=zero
789  igwene=zero
790 
791  do is=1,Bst%nsppol
792    seen=0
793 
794    do ikr=1,nkibzR
795      read(unt,*)kread(:)
796      read(unt,*)nbandR
797      ikibz=0
798      do ik=1,Bst%nkpt
799        if (ALL(ABS(kread(:)-Bst%kptns(:,ik))<0.0001)) then
800          ikibz=ik
801          seen(ik) = seen(ik) + 1
802        end if
803      end do
804      do ib=1,nbandR
805        read(unt,*)ibr,egw_r,degw,egw_i
806        if (ibr<=Bst%mband .and. ikibz/=0) then
807          gwcorr(ibr,ikibz,is)=degw/Ha_eV
808          igwene(ibr,ikibz,is)=egw_i/Ha_eV
809        end if
810      end do
811    end do
812 
813    if (ANY(seen/=1)) then
814      do ik=1,Bst%nkpt
815        if (seen(ik)/=1) then
816          write(msg,'(a,3f8.3,a)')" k-point: ",Bst%kptns(:,ik)," not found in the GW file!"
817          MSG_WARNING(msg)
818        end if
819      end do
820    end if
821 
822  end do
823 
824  ABI_FREE(seen)
825  close(unt)
826 
827  do_extrapolate=.TRUE.; if (PRESENT(extrapolate)) do_extrapolate=extrapolate
828 
829  if (.not. do_extrapolate) then ! Only the bands calculated are updated.
830    Bst%eig = Bst%eig + gwcorr
831 
832  else
833 
834    if (ANY(ABS(igwene)>tol6)) then
835      write(msg,'(4a)')ch10,&
836 &      "The GW file contains QP energies with non-zero imaginary part",ch10,&
837 &      "Extrapolation not coded, change the source! "
838      MSG_ERROR(msg)
839    end if
840 
841    ABI_MALLOC(vbik,(BSt%nkpt,BSt%nsppol))
842    vbik(:,:) = get_valence_idx(BSt)
843 
844    do is=1,Bst%nsppol
845      do ik=1,Bst%nkpt
846 
847       nbv=vbik(ik,is) ! Index of the (valence band| Fermi band) for each spin
848       nn=Bst%mband-nbv
849 
850       do ib=nbv+1,Bst%mband
851         if ( ABS(gwcorr(ib,ik,is)) < tol16) then
852           nn=ib-1-nbv
853           if (nn>1) then
854             call wrtout(std_out,"Linear extrapolating (conduction) GW corrections beyond the read values","COLL")
855             smrt=linfit(nn,Bst%eig(nbv+1:nbv+nn,ik,is),gwcorr(nbv+1:nbv+nn,ik,is),alpha,beta)
856           else
857             call wrtout(std_out,"Assuming constant (conduction) GW corrections beyond the read values",'COLL')
858             alpha=zero
859             beta =gwcorr(nbv+nn,ik,is)
860           end if
861           EXIT !ib loop
862         end if
863       end do !ib
864 
865       do ib=nbv+nn+1,Bst%mband
866         gwcorr(ib,ik,is)= alpha*Bst%eig(ib,ik,is) + beta
867       end do
868 
869       nn=nbv
870       do ib=nbv,1,-1
871         if ( ABS(gwcorr(ib,ik,is)) < tol16) then
872          nn=nbv-ib
873          if (nn>1) then
874            call wrtout(std_out,"Linear extrapolating (valence) GW corrections beyond the read values","COLL")
875            smrt=linfit(nn,Bst%eig(nbv-nn+1:nbv,ik,is),gwcorr(nbv-nn+1:nbv,ik,is),alpha,beta)
876          else
877            call wrtout(std_out,"Assuming constant (valence) GW corrections beyond the read values","COLL")
878            alpha=zero
879            beta =gwcorr(nbv,ik,is)
880          end if
881          EXIT !ib
882         end if
883       end do !ib
884 
885       do ib=1,nbv-nn
886         gwcorr(ib,ik,is)=alpha*Bst%eig(ib,ik,is) + beta
887       end do
888 
889      end do !ik
890    end do !is
891 
892    call wrtout(std_out,' k  s     GW corrections [eV] ','COLL')
893    do is=1,Bst%nsppol
894      do ik=1,Bst%nkpt
895        write(msg,'(i3,1x,i3,10f7.2/50(10x,10f7.2/))')ik,is,(Ha_eV*gwcorr(ib,ik,is),ib=1,Bst%mband)
896        call wrtout(std_out,msg,"COLL")
897      end do
898    end do
899    Bst%eig = Bst%eig + gwcorr
900    ABI_FREE(vbik)
901  end if
902 
903  call wrtout(std_out,' k   s    GW eigenvalues [eV]',"COLL")
904  do is=1,Bst%nsppol
905    do ik=1,Bst%nkpt
906      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)
907    end do
908  end do
909 
910  ABI_FREE(gwcorr)
911 
912 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_lda_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_lda_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.

PARENTS

      bethe_salpeter,mlwfovlp_qp,screening,sigma

CHILDREN

SOURCE

283 subroutine rdqps(BSt,fname,usepaw,nspden,dimrho,nscf,&
284 & nfftot,ngfftf,ucvol,paral_kgb,Cryst,Pawtab,MPI_enreg,nbsc,m_lda_to_qp,rhor_out,Pawrhoij)
285 
286 
287 !This section has been created automatically by the script Abilint (TD).
288 !Do not modify the following lines by hand.
289 #undef ABI_FUNC
290 #define ABI_FUNC 'rdqps'
291 !End of the abilint section
292 
293  implicit none
294 
295 !Arguments ------------------------------------
296 !scalars
297  integer,intent(in) :: nfftot,nspden,usepaw,paral_kgb,dimrho
298  integer,intent(out) :: nbsc,nscf
299  real(dp),intent(in) :: ucvol
300  character(len=*),intent(in) :: fname
301  type(crystal_t),intent(in) :: Cryst
302  type(ebands_t),intent(inout) :: BSt
303  type(MPI_type),intent(inout) :: MPI_enreg
304 !arrays
305  integer,intent(in) :: ngfftf(18)
306  real(dp),intent(out) :: rhor_out(nfftot,nspden*dimrho)
307  complex(dpc),intent(out) :: m_lda_to_qp(BSt%mband,BSt%mband,BSt%nkpt,BSt%nsppol)
308  type(Pawtab_type),intent(in) :: Pawtab(Cryst%ntypat*usepaw)
309  type(Pawrhoij_type),intent(inout) :: Pawrhoij(Cryst%natom*usepaw)
310 
311 !Local variables-------------------------------
312 !scalars
313  integer,parameter :: master=0
314  integer :: ib,ii,ik,isppol,nbandR,nkibzR,nsppolR,unqps,my_rank,ispden
315  integer :: ifft,n1,n2,n3,ir1,ir2,ir3,ios
316  integer :: cplex_fft,optin,optout,nfft_found
317  integer :: iatom,natomR,nspdenR,ntypatR,itypat
318  real(dp) :: uerr,nelect_qps,ratio
319  logical,parameter :: use_FFT_interpolation=.TRUE.
320  logical :: ltest
321  character(len=500) :: msg
322 !arrays
323  integer :: ngfft_found(18)
324  integer,allocatable :: nlmn_type(:),typatR(:)
325  real(dp) :: kibz(3),rr(3),rhogdum(1,1)
326  real(dp),allocatable :: en_tmp(:)
327  real(dp),allocatable :: rhor_tmp(:,:)
328  complex(dpc),allocatable :: mtmp(:,:),utest(:,:)
329 
330 ! *************************************************************************
331 
332  DBG_ENTER("COLL")
333 
334  ABI_CHECK(ALL(BSt%nband==BSt%nband(1)),"No. of bands must be constant")
335  ABI_CHECK(dimrho==0.or.dimrho==1,'dimrho must be 0 or 1')
336 
337  ! This does not work in parallel !!?
338  !% my_rank = xmpi_comm_rank(MPI_enreg%spaceComm)
339  my_rank = MPI_enreg%me_kpt
340 
341  ! Check whether file exists or not.
342  write(msg,'(5a)')ch10,&
343 &  ' rdqps: reading QP wavefunctions of the previous step ',ch10,&
344 &  '        looking for file ',TRIM(fname)
345  call wrtout(std_out,msg,'COLL')
346  call wrtout(ab_out,msg,'COLL')
347 
348  if (.not.file_exists(fname)) then
349    write(msg,'(2a)')' file not found, 1st iteration initialized with KS eigenelements ',ch10
350    call wrtout(std_out,msg,'COLL')
351    call wrtout(ab_out, msg,'COLL')
352    nscf=0; RETURN
353  end if
354 
355  if (.not.isncfile(fname)) then
356    if (open_file(fname,msg,newunit=unqps,form='formatted',status='unknown') /= 0) then
357      MSG_ERROR(msg)
358    end if
359 
360    ! TODO the _QPS file should contain additional information
361    read(unqps,*)nscf
362    write(msg,'(a,i4,a)')' Number of iteration(s) already performed: ',nscf,ch10
363    call wrtout(std_out,msg,'COLL')
364    call wrtout(ab_out,msg,'COLL')
365 
366    read(unqps,*)nkibzR
367    if (nkibzR/=BSt%nkpt) then
368      write(msg,'(2(a,i0))')'Wrong number of k-points; Expected: ',BSt%nkpt,', Found: ',nkibzR
369      MSG_ERROR(msg)
370    end if
371 
372    read(unqps,*)nbandR
373    nbsc=MIN(nbandR,BSt%mband)
374 
375    if (nbsc/=BSt%mband) then
376      write(msg,'(3a,i4,a,i4)')&
377 &      'QPS file contains less bands than that used in the present calculation ',ch10,&
378 &      'Required: ',BSt%mband,', Found: ',nbandR
379      MSG_WARNING(msg)
380    end if
381 
382    if (nbsc/=nbandR) then
383      write(msg,'(3a,i4,a)')&
384 &      'The QPS file contains more bands than that used in the present calculation ',ch10,&
385 &      'only the first ',nbandR,' bands will be read'
386      MSG_COMMENT(msg)
387    end if
388 
389    ABI_MALLOC(mtmp,(nbandR,nbandR))
390    ABI_MALLOC(en_tmp,(nbandR))
391    read(unqps,*)nsppolR
392 
393    ABI_CHECK(nsppolR==BSt%nsppol,"QPS generated with different nsppol")
394 
395    ! Read energies and transformation for each k-point and spin.
396    ! TODO: The format of the QPS file must be standardized !
397    ! For example we might add the occupation numbers.
398    do isppol=1,BSt%nsppol
399      do ik=1,BSt%nkpt
400        read(unqps,*)kibz(:)
401        write(msg,'(a,i5,a,3(f6.3,1x),4x,a,i2)')' Reading ik ',ik,')  k = ',kibz(:),' is = ',isppol
402        call wrtout(std_out,msg,'COLL')
403        ltest=(ALL(ABS(kibz(:)-BSt%kptns(:,ik))<0.001))
404        ABI_CHECK(ltest,'Wrong k-point read')
405        do ib=1,nbandR
406          read(unqps,*)en_tmp(ib)
407          read(unqps,*)mtmp(:,ib)
408        end do
409 
410        ! Store transformation and update energies.
411        m_lda_to_qp(1:nbsc,1:nbsc,ik,isppol)=mtmp(1:nbsc,1:nbsc)
412        BSt%eig(1:nbsc,ik,isppol)=en_tmp(1:nbsc)
413 
414        ! Chech if matrix is unitary.
415        ABI_MALLOC(utest,(nbsc,nbsc))
416        utest(:,:) = TRANSPOSE(mtmp(1:nbsc,1:nbsc)) !this is just for the buggy gfortran
417        utest(:,:) = MATMUL(CONJG(utest),mtmp(1:nbsc,1:nbsc))
418        do ii=1,nbsc
419          utest(ii,ii)=utest(ii,ii)-one
420        end do
421        uerr=MAXVAL(ABS(utest))
422        if (uerr>tol6) then
423          write(msg,'(a,es16.8)')' KS -> QP matrix is not unitary, MAX error = ',uerr
424          MSG_WARNING(msg)
425        end if
426        ABI_FREE(utest)
427      end do !ik
428    end do !isppol
429 
430    ABI_FREE(mtmp)
431    ABI_FREE(en_tmp)
432 
433    ! Read the QP density.
434    ! The two FFT grids might differ. In case perform an FFT interpolation to have rhor on the input mesh.
435    if (dimrho==1) then
436      read(unqps,*)n1,n2,n3
437 
438      if (all(ngfftf(1:3)== [n1, n2, n3]) ) then
439        read(unqps,*)rhor_out(:,:)
440      else
441        write(msg,'(2a,a,5(i3,a),i3)')&
442 &        'FFT meshes differ. Performing Fourier interpolation. ',ch10,&
443 &        'Found: ',n1,' x',n2,' x',n3,'; Expected: ',ngfftf(1),' x',ngfftf(2),' x',ngfftf(3)
444        MSG_COMMENT(msg)
445 
446        ABI_MALLOC(rhor_tmp,(n1*n2*n3,nspden))
447        read(unqps,*)rhor_tmp(:,:)
448 
449        if (use_FFT_interpolation) then
450          ngfft_found(1:3)=(/n1,n2,n3/)
451          ngfft_found(4)=2*(ngfft_found(1)/2)+1 ! 4:18 are not used, anyway!
452          ngfft_found(5)=2*(ngfft_found(2)/2)+1
453          ngfft_found(6)=ngfft_found(3)
454          ngfft_found(7:18)=ngfftf(7:18)
455          nfft_found=PRODUCT(ngfft_found(1:3)) !no FFT para
456 
457          cplex_fft =1 ! Real quantities.
458          optin     =0 ! Input is taken from rhor.
459          optout    =0 ! Output is only in real space.
460          call destroy_distribfft(MPI_enreg%distribfft)
461          call init_distribfft(MPI_enreg%distribfft,'c',MPI_enreg%nproc_fft,ngfftf(2),ngfftf(3))
462          call init_distribfft(MPI_enreg%distribfft,'f',MPI_enreg%nproc_fft,ngfft_found(2),ngfft_found(3))
463 
464          call fourier_interpol(cplex_fft,nspden,optin,optout,nfft_found,ngfft_found,nfftot,ngfftf,&
465 &          paral_kgb,MPI_enreg,rhor_tmp,rhor_out,rhogdum,rhogdum)
466 
467        else
468          ! Linear interpolation.
469          do ispden=1,nspden
470            do ir3=0,ngfftf(3)-1
471              rr(3)=DBLE(ir3)/n3
472              do ir2=0,ngfftf(2)-1
473                rr(2)=DBLE(ir2)/n2
474                do ir1=0,ngfftf(1)-1
475                  rr(1)=DBLE(ir1)/n1
476                  ifft = 1 +ir1 +ir2*ngfftf(1) +ir3*ngfftf(1)*ngfftf(2)
477                  rhor_out(ifft,ispden) = interpol3d(rr,n1,n2,n3,rhor_tmp(:,ispden))
478                end do
479              end do
480            end do
481          end do
482        end if
483 
484        ABI_FREE(rhor_tmp)
485      end if
486 
487      ! Test the normalization of the QPS density.
488      ! There might be errors due to the interpolation or the truncation of the G basis set
489      ! Density will be renormalized in the caller since for PAW we still have to add the onsite contribution.
490      if (usepaw==0) then
491        nelect_qps=SUM(rhor_out(:,1))*ucvol/nfftot; ratio=BSt%nelect/nelect_qps
492        write(msg,'(3(a,f9.4))')&
493 &        ' Number of electrons calculated using the QPS density = ',nelect_qps,' Expected = ',BSt%nelect,' ratio = ',ratio
494        call wrtout(std_out,msg,'COLL')
495        !!rhor_out(:,:)=ratio*rhor_out(:,:)
496      end if
497 
498      if (usepaw==1) then
499        ! Write QP_rhoij for on-site density mixing.
500        read(unqps,*,iostat=ios)natomR,ntypatR
501        if (ios/=0) then
502          msg="Old version of QPS file found. DO NOT USE rhoqpmix for this run."
503          MSG_WARNING(msg)
504          call wrtout(ab_out,msg,"COLL")
505          ! Init dummy rhoij just to avoid problems in sigma when rhoij is freed.
506          call pawrhoij_alloc(Pawrhoij,1,nspden,BSt%nspinor,BSt%nsppol,Cryst%typat,pawtab=Pawtab)
507          close(unqps)
508          RETURN
509        end if
510 
511        ABI_CHECK(natomR ==Cryst%natom, "mismatch in natom")
512        ABI_CHECK(ntypatR==Cryst%ntypat,"mismatch in ntypat")
513        ABI_MALLOC(nlmn_type,(ntypatR))
514        ABI_MALLOC(typatR,(ntypatR))
515 
516        read(unqps,*)(typatR(iatom), iatom=1,natomR)
517        ABI_CHECK(ALL(Cryst%typat==typatR),"mismatch in typat")
518 
519        read(unqps,*)(nlmn_type(itypat), itypat=1,ntypatR)
520        do itypat =1,Cryst%ntypat
521          if (nlmn_type(itypat)/=Pawtab(itypat)%lmn_size) then
522            MSG_ERROR("mismatch in nlmn_type, check QPS file")
523          end if
524        end do
525 
526        read(unqps,*) nsppolR,nspdenR
527        ABI_CHECK(nsppolR==BSt%nsppol,"mismatch in nsppol")
528        ABI_CHECK(nspdenR==nspden    ,"mismatch in nspden")
529 
530        call pawrhoij_io(pawrhoij,unqps,BSt%nsppol,BSt%nspinor,nspden,nlmn_type,Cryst%typat,&
531 &                    HDR_LATEST_HEADFORM,"Read",form="formatted")
532        !% call pawrhoij_io(pawrhoij,std_out,BSt%nsppol,BSt%nspinor,nspden,nlmn_type,Cryst%typat,HDR_LATEST_HEADFORM,"Echo")
533 
534        ABI_FREE(nlmn_type)
535        ABI_FREE(typatR)
536      end if ! usepaw
537 
538    end if !dimrho=1
539 
540    close(unqps)
541 
542  else
543    MSG_ERROR("netdf format not implemented")
544  end if
545 
546  DBG_EXIT("COLL")
547 
548 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_lda_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_lda_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.

PARENTS

      sigma

CHILDREN

SOURCE

588 subroutine show_QP(Bst,m_lda_to_qp,fromb,tob,unit,prtvol,tolmat,kmask)
589 
590 
591 !This section has been created automatically by the script Abilint (TD).
592 !Do not modify the following lines by hand.
593 #undef ABI_FUNC
594 #define ABI_FUNC 'show_QP'
595 !End of the abilint section
596 
597  implicit none
598 
599 !Arguments ------------------------------------
600 !scalars
601  integer,optional,intent(in) :: fromb,tob
602  integer,optional,intent(in) :: prtvol,unit
603  real(dp),optional,intent(in) :: tolmat
604  type(ebands_t),intent(in) :: Bst
605 !arrays
606  logical,optional,intent(in) :: kmask(Bst%nkpt)
607  complex(dpc),intent(in) :: m_lda_to_qp(Bst%mband,Bst%mband,Bst%nkpt,Bst%nsppol)
608 
609 !Local variables-------------------------------
610 !scalars
611  integer,parameter :: NBRA=5
612  logical,parameter :: use_rhophi=.True.
613  integer :: ib_start,ib_stop,my_prtvol,counter,ib_KS,ib_QP,ikibz,isp,nspace,my_unt,nband_k
614  real(dp) :: my_tolmat,rho,phi
615  character(len=10) :: bks,bqp,k_tag,spin_tag
616  character(len=500) :: KS_row,KS_ket,tmpstr,QP_ket
617 !arrays
618  real(dp) :: cx(2)
619 
620 ! *********************************************************************
621 
622  my_unt   =std_out  ; if (PRESENT(unit  )) my_unt   =unit
623  my_prtvol=0        ; if (PRESENT(prtvol)) my_prtvol=prtvol
624  ib_start =1        ; if (PRESENT(fromb )) ib_start =fromb
625  ib_stop  =Bst%mband; if (PRESENT(tob   )) ib_stop  =tob
626  my_tolmat=0.001    ; if (PRESENT(tolmat)) my_tolmat=ABS(tolmat)
627 
628  ! I suppose nband_k is constant thus the check is done here.
629  if (ib_start<=0       ) ib_start=1
630  if (ib_start>Bst%mband) ib_start=Bst%mband
631  if (ib_stop<=0        ) ib_stop=1
632  if (ib_stop>Bst%mband ) ib_stop=Bst%mband
633 
634  ! Have to follow rules 7.f.
635  write(my_unt,'(/,a,/,a,/,a,f6.3,a,/,a)')&
636    ' '//REPEAT('*',76),&
637 &  ' ***** QP amplitudes expressed as linear combination of KS eigenstates. *****',&
638 &  ' ***** Only KS components whose modulus is larger than ',my_tolmat,' are shown  ***** ',&
639 &  ' '//REPEAT('*',76)
640  if (use_rhophi) then
641    write(my_unt,"(a)")"Complex coefficients given in (rho, phi) polar representation."
642  else
643    write(my_unt,"(a)")"Complex coefficients given in (Re, Im) representation."
644  end if
645 
646  if (PRESENT(kmask)) then
647    if (.not.ALL(kmask)) write(my_unt,'(/,a,i3,a)')' Only ',COUNT(kmask),' k-points are reported '
648  end if
649 
650  do isp=1,Bst%nsppol
651    call int2char10(isp,spin_tag)
652    write(my_unt,'(/,a,i2,a,/)')' >>>>> Begin block for spin ',isp,' <<<<< '
653 
654    do ikibz=1,Bst%nkpt
655      if (PRESENT(kmask)) then
656        if (.not.kmask(ikibz)) CYCLE
657      end if
658      call int2char10(ikibz,k_tag)
659      nband_k=Bst%nband(ikibz+(isp-1)*Bst%nkpt)
660      write(my_unt,'(a,i4,a,3es16.8,a,f6.3,/)')' k-point: ',ikibz,') ',Bst%kptns(:,ikibz),'; wtk= ',Bst%wtk(ikibz)
661 
662      do ib_QP=ib_start,ib_stop
663        call int2char10(ib_QP,bqp)
664        QP_ket=' |QP: b='//TRIM(bqp)//'; s='//TRIM(spin_tag)//'> = '
665        write(my_unt,'(a)')TRIM(QP_ket)
666        nspace=LEN(TRIM(QP_ket))
667 
668        counter=0 ; KS_row=REPEAT('',nspace+2)
669        do ib_KS=1,Bst%mband
670          if (ABS(m_lda_to_qp(ib_KS,ib_QP,ikibz,isp))<my_tolmat) CYCLE
671          counter=counter+1
672          call int2char10(ib_KS,bks)
673          write(tmpstr,'(3a)')' |',TRIM(bks),'>'
674 
675          if (use_rhophi) then
676            ! coefficient as (rho, phi)
677            cx(1) = real(m_lda_to_qp(ib_KS,ib_QP,ikibz,isp))
678            cx(2) = aimag(m_lda_to_qp(ib_KS,ib_QP,ikibz,isp))
679            call rhophi(cx, phi, rho)
680            write(KS_ket,'(1x,2f7.3,a,1x)')rho, phi, TRIM(tmpstr)
681          else
682            ! coefficient as (Re, Im)
683            write(KS_ket,'(1x,2f7.3,a,1x)')m_lda_to_qp(ib_KS,ib_QP,ikibz,isp),TRIM(tmpstr)
684          end if
685          KS_row=TRIM(KS_row)//TRIM(KS_ket)
686          if (MOD(counter,NBRA)==0) then  ! nbra KS kets per row
687            write(my_unt,'(a)')TRIM(KS_row)
688            KS_row=REPEAT('',nspace+2)
689          end if
690        end do
691 
692        if (MOD(counter,NBRA)/=0) write(my_unt,'(a)')TRIM(KS_row) ! Last row, if any
693        write(my_unt,'(a)')''
694      end do !ib_QP
695 
696    end do !ikibz
697  end do !isp
698 
699  write(my_unt,'(a,/)')' '//REPEAT('*',76)
700 
701 end subroutine show_QP

m_qparticles/updt_m_lda_to_qp [ Functions ]

[ Top ] [ m_qparticles ] [ Functions ]

NAME

 updt_m_lda_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_lda_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.

PARENTS

      sigma

CHILDREN

SOURCE

955 subroutine updt_m_lda_to_qp(Sigp,Kmesh,nscf,Sr,m_lda_to_qp)
956 
957 
958 !This section has been created automatically by the script Abilint (TD).
959 !Do not modify the following lines by hand.
960 #undef ABI_FUNC
961 #define ABI_FUNC 'updt_m_lda_to_qp'
962 !End of the abilint section
963 
964  implicit none
965 
966 !Arguments ------------------------------------
967 !scalars
968  integer,intent(in) :: nscf
969  type(kmesh_t),intent(in) :: Kmesh
970  type(sigparams_t),intent(in) :: Sigp
971  type(sigma_t),intent(in) :: Sr
972 !arrays
973  complex(dpc),intent(inout) :: m_lda_to_qp(Sigp%nbnds,Sigp%nbnds,Kmesh%nibz,Sigp%nsppol)
974 
975 !Local variables-------------------------------
976 !scalars
977  integer :: ik,is
978 !arrays
979  complex(dpc),allocatable :: mtmp(:,:)
980 
981 ! *************************************************************************
982 
983  if (nscf>=0) then
984    ! Calculate the new m_lda_to_qp
985    ABI_MALLOC(mtmp,(Sigp%nbnds,Sigp%nbnds))
986    do is=1,Sigp%nsppol
987      do ik=1,Kmesh%nibz
988        mtmp(:,:)=m_lda_to_qp(:,:,ik,is)
989        m_lda_to_qp(:,:,ik,is)=MATMUL(mtmp(:,:),Sr%eigvec_qp(:,:,ik,is))
990      end do
991    end do
992    ABI_FREE(mtmp)
993  end if
994 
995 end subroutine updt_m_lda_to_qp
996 
997 !----------------------------------------------------------------------
998 
999 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_lda_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.

PARENTS

      sigma

CHILDREN

SOURCE

119 subroutine wrqps(fname,Sigp,Cryst,Kmesh,Psps,Pawtab,Pawrhoij,nspden,nscf,nfftot,ngfftf,Sr,Bst,m_lda_to_qp,rho_qp)
120 
121 
122 !This section has been created automatically by the script Abilint (TD).
123 !Do not modify the following lines by hand.
124 #undef ABI_FUNC
125 #define ABI_FUNC 'wrqps'
126 !End of the abilint section
127 
128  implicit none
129 
130 !Arguments ------------------------------------
131 !scalars
132  integer,intent(in) :: nfftot,nscf,nspden
133  character(len=*),intent(in) :: fname
134  type(kmesh_t),intent(in) :: Kmesh
135  type(ebands_t),intent(in) :: BSt
136  type(sigparams_t),intent(in) :: Sigp
137  type(sigma_t),intent(in) :: Sr
138  type(crystal_t),intent(in) :: Cryst
139  type(Pseudopotential_type),intent(in) :: Psps
140 !arrays
141  integer,intent(in) :: ngfftf(18)
142  real(dp),intent(in) :: rho_qp(nfftot,nspden)
143  complex(dpc),intent(in) :: m_lda_to_qp(Sigp%nbnds,Sigp%nbnds,Kmesh%nibz,Sigp%nsppol)
144  type(Pawrhoij_type),intent(inout) :: Pawrhoij(Cryst%natom*Psps%usepaw)
145  type(Pawtab_type),intent(in) :: Pawtab(Psps%ntypat*Psps%usepaw)
146 
147 !Local variables-------------------------------
148 !scalars
149  integer :: ib,ik,is,unqps,iatom,itypat
150  character(len=500) :: msg
151 !arrays
152  integer,allocatable :: nlmn_type(:)
153  complex(dpc),allocatable :: mtmp(:,:)
154 
155 ! *************************************************************************
156 
157  DBG_ENTER("COLL")
158 
159  if (nscf >= 0) then
160    write(msg,'(3a)')ch10,' writing QP data on file : ',TRIM(fname)
161    call wrtout(std_out,msg,'COLL')
162    call wrtout(ab_out,msg,'COLL')
163  end if
164 
165  if (open_file(fname,msg,newunit=unqps,form='formatted',status='unknown') /= 0) then
166    MSG_ERROR(msg)
167  end if
168 
169  write(unqps,*)nscf+1
170  write(unqps,*)Kmesh%nibz
171  write(unqps,*)Sigp%nbnds
172  write(unqps,*)Sigp%nsppol
173 
174  ABI_MALLOC(mtmp,(Sigp%nbnds,Sigp%nbnds))
175 
176  if (nscf>=0) then
177    ! Write the new m_lda_to_qp on file.
178    do is=1,Sigp%nsppol
179      do ik=1,Kmesh%nibz
180        write(unqps,*)Kmesh%ibz(:,ik)
181        do ib=1,Sigp%nbnds
182          write(unqps,*)Sr%en_qp_diago(ib,ik,is)
183          write(unqps,*)m_lda_to_qp(:,ib,ik,is)
184        end do
185      end do
186    end do
187  else if (nscf==-1) then
188    ! Write fake QPS file with KS band structure (Mainly used for G0W)
189    call set2unit(mtmp)
190    do is=1,Sigp%nsppol
191      do ik=1,Kmesh%nibz
192        write(unqps,*)Kmesh%ibz(:,ik)
193        do ib=1,Sigp%nbnds
194          write(unqps,*)BSt%eig(ib,ik,is)
195          write(unqps,*)mtmp(:,ib)
196        end do
197      end do
198    end do
199  else
200    MSG_ERROR(sjoin("Wrong nscf ",itoa(nscf)))
201  end if
202 
203  ABI_FREE(mtmp)
204 
205  write(msg,'(a,f9.4)')' (wrqps) planewave contribution to nelect: ',SUM(rho_qp(:,1))*Cryst%ucvol/nfftot
206  call wrtout(std_out,msg,'COLL')
207  if (nspden == 4) then
208    write(msg,'(a,3f9.4)')' mx, my, mz: ',&
209      SUM(rho_qp(:,2))*Cryst%ucvol/nfftot,SUM(rho_qp(:,3))*Cryst%ucvol/nfftot,SUM(rho_qp(:,4))*Cryst%ucvol/nfftot
210    call wrtout(std_out,msg,'COLL')
211  end if
212 
213  ! Write FFT dimensions and QP density
214  write(unqps,*)ngfftf(1:3)
215  write(unqps,*)rho_qp(:,:)
216 
217  if (Psps%usepaw==1) then
218    ! Write QP rhoij to be used for on-site density mixing.
219    ABI_MALLOC(nlmn_type,(Cryst%ntypat))
220    do itypat=1,Cryst%ntypat
221      nlmn_type(itypat)=Pawtab(itypat)%lmn_size
222    end do
223 
224    write(unqps,*) Cryst%natom, Cryst%ntypat
225    write(unqps,*) (Cryst%typat(iatom), iatom=1,Cryst%natom)
226    write(unqps,*) (nlmn_type(itypat), itypat=1,Cryst%ntypat)
227    write(unqps,*) Pawrhoij(1)%nsppol, Pawrhoij(1)%nspden
228 
229    call pawrhoij_io(pawrhoij,unqps,Sigp%nsppol,Sigp%nspinor,nspden,nlmn_type,Cryst%typat,&
230 &                HDR_LATEST_HEADFORM,"Write",form="formatted")
231    ABI_FREE(nlmn_type)
232  end if
233 
234  close(unqps)
235 
236  DBG_EXIT("COLL")
237 
238 end subroutine wrqps