TABLE OF CONTENTS
- ABINIT/Calc_Ec_GM_k
- ABINIT/calc_rdmc
- ABINIT/calc_rdmx
- ABINIT/change_matrix
- ABINIT/get_chkprdm
- ABINIT/ks2no
- ABINIT/m_gwrdm
- ABINIT/natoccs
- ABINIT/no2ks
- ABINIT/print_band_energies
- ABINIT/print_chkprdm
- ABINIT/print_tot_occ
- ABINIT/print_total_energy
- ABINIT/printrdm_k
- ABINIT/quadrature_sigma_cw
- ABINIT/rotate_ks_no
- ABINIT/update_hdr_bst
ABINIT/Calc_Ec_GM_k [ Functions ]
NAME
calc_Ec_GM_k
FUNCTION
Calculate Galitskii-Migdal corr. energy integrated in the Imaginary axis Ec = 1/pi \sum_i \int Gii(iv)*Sigma_c,ii(iv) + cc. dv
INPUTS
ib1=min band for given k ib2=max band for given k. ik_ibz= the label of k-point in the IBZ whose Galitskii-Migdal contribution is accounted. weights=array containing the weights used in the quadrature. sigcme_k=array containing Sigma(iw) as Sigma(iw,ib1:ib2,ib1:ib2,nspin) rdm_k=density matrix, matrix (i,j), where i and j belong to the k-point k (see m_sigma_driver for more details). ebands=<ebands_t>=Datatype gathering info on the QP energies (KS if one shot) eig(Sigp%nbnds,Kmesh%nibz,Wfd%nsppol)=KS or QP energies for k-points, bands and spin occ(Sigp%nbnds,Kmesh%nibz,Wfd%nsppol)=occupation numbers, for each k point in IBZ, each band and spin Sr=sigma_t (see the definition of this structured datatype)
OUTPUT
Compute the Galitskii-Migdal corr energy contribution of this k-point: Ec ^k = 1/(4*pi) * fact_spin * int _{ -Inf }^{ +Inf } dv Sigma_c ^k (iv) * G0(iv) = 1/(4*pi) * fact_spin * int _{ 0 }^{ +Inf } dv 2 * Re{ Sigma_c ^k (iv) * G0(iv) }
SOURCE
144 function calc_Ec_GM_k(ib1,ib2,ik_ibz,Sr,weights,sigcme_k,ebands) result(Ec_GM_k) 145 146 !Arguments ------------------------------------ 147 !scalars 148 real(dp) :: Ec_GM_k 149 integer,intent(in) :: ib1,ib2,ik_ibz 150 type(ebands_t),target,intent(in) :: ebands 151 type(sigma_t),intent(in) :: Sr 152 !arrays 153 real(dp),intent(in) :: weights(:) 154 complex(dpc),intent(in) :: sigcme_k(:,:,:,:) 155 !Local variables ------------------------------ 156 !scalars 157 integer :: ibdm!,unitt 158 real(dp) :: ec_integrated,spin_fact,fact 159 character(len=500) :: msg 160 !arrays 161 !************************************************************************ 162 163 ec_integrated=zero 164 spin_fact=two 165 fact=spin_fact*(one/(two_pi*two)) 166 167 if (ib1/=1) then 168 msg="Unable to compute the Galitskii-Migdal correlation energy because the first band was " // & 169 & "not included in bdgw interval. Restart the calculation starting bdgw from 1." 170 ABI_WARNING(msg) 171 else 172 ! WARNING: Sigma_c(iv) produced from a previous integration at the screening stage, is numerically not much stable and introduces bumps. 173 ! Unfortunately, the Green's function times Sigma_c(iv) does not decay fast enough with iv to overcome the bumps. These bumps are 174 ! not pronouced for the linearized density matrix update, as two Green's functions are multiplied making the decay much faster with iv. 175 ! If a better way to produce more stable Sigma_c(iv) values is found, this subroutine can be use to evaluate GM Ecorr in the future. TODO 176 do ibdm=1,ib2 177 ! Sigma_pp(iv)/[(iv - e_ibdm,k)] + [Sigma_pp(iv)/[(iv - e_ibdm,k)]]^* = 2 Re [Sigma_pp(iv)/(iv - e_ibdm,k)] 178 ec_integrated=ec_integrated+two*real( sum(weights(:)*sigcme_k(:,ibdm,ibdm,1)/(Sr%omega_i(:)-ebands%eig(ibdm,ik_ibz,1)) ) ) 179 end do 180 endif 181 182 Ec_GM_k=fact*ec_integrated 183 184 end function calc_Ec_GM_k
ABINIT/calc_rdmc [ Functions ]
NAME
calc_rdmc
FUNCTION
Calculate density matrix corrections for G = Go + int Go(iw) Sigma_c(iw) Go(iw) dw
INPUTS
ib1=min band for given k ib2=max band for given k. ik_ibz= the label of k-point in the IBZ. omega_i=Frequencies along the imaginary axis. weights=array containing the weights used in the quadrature. sigcme_k=array containing Sigma(iw) as Sigma(iw,ib1:ib2,ib1:ib2,nspin) rdm_k=density matrix, matrix (i,j), where i and j belong to the k-point k (see m_sigma_driver for more details). ebands=<ebands_t>=Datatype gathering info on the QP energies (KS if one shot)
OUTPUT
Updated rdm_k matrix array with int Go(iw) Sigma_c(iw) Go(iw) dw
SOURCE
272 subroutine calc_rdmc(ib1,ib2,ik_ibz,omega_i,weights,sigcme_k,ebands,rdm_k) 273 274 !Arguments ------------------------------------ 275 !scalars 276 integer,intent(in) :: ib1,ib2,ik_ibz 277 type(ebands_t),target,intent(in) :: ebands 278 complex(dpc),intent(in) :: omega_i(:) 279 !arrays 280 real(dp),intent(in) :: weights(:) 281 complex(dpc),intent(inout) :: rdm_k(:,:) 282 complex(dpc),intent(in) :: sigcme_k(:,:,:,:) 283 284 !Local variables ------------------------------ 285 !scalars 286 real(dp) :: spin_fact,fact 287 integer :: ib1dm, ib2dm, units(2) 288 character(len=500) :: msg 289 !************************************************************************ 290 291 spin_fact = two 292 fact = spin_fact * (one/two_pi) 293 units = [std_out, ab_out] 294 295 write(msg,'(a58,3f10.5)')' Computing the 1-RDM correction for Sc(iw) and k-point: ',ebands%kptns(1:,ik_ibz) 296 call wrtout(units , msg) 297 write(msg,'(a11,i5,a8,i5)')'from band ',ib1,' to band',ib2 298 call wrtout(units , msg) 299 300 rdm_k(:,:)=czero 301 do ib1dm=ib1,ib2 302 do ib2dm=ib1dm,ib2 303 ! Sigma_pq/[(denominator)] + [Sigma_qp/[(denominator)]]^* 304 rdm_k(1+(ib1dm-ib1),1+(ib2dm-ib1))=fact*sum(weights(:)*( sigcme_k(:,1+(ib1dm-ib1),1+(ib2dm-ib1),1)/& 305 &( (omega_i(:)-ebands%eig(ib1dm,ik_ibz,1))*(omega_i(:)-ebands%eig(ib2dm,ik_ibz,1)) )& 306 +conjg( sigcme_k(:,1+(ib2dm-ib1),1+(ib1dm-ib1),1)/& 307 &( (omega_i(:)-ebands%eig(ib1dm,ik_ibz,1))*(omega_i(:)-ebands%eig(ib2dm,ik_ibz,1)) ) ) ) ) 308 ! Dji = Dij^* 309 rdm_k(1+(ib2dm-ib1),1+(ib1dm-ib1))=conjg(rdm_k(1+(ib1dm-ib1),1+(ib2dm-ib1))) 310 end do 311 end do 312 313 end subroutine calc_rdmc
ABINIT/calc_rdmx [ Functions ]
NAME
calc_rdmx
FUNCTION
Calculate density matrix corrections for G = Go + Go (Sigma_x - alpha*Sigma_x - Vxc) Go
INPUTS
ib1=min band for given k ib2=max band for given k. ik_ibz= the label of k-point in the IBZ. rdm_k=density matrix, matrix (i,j), where i and j belong to the k-point k (see m_sigma_driver for more details). pot=Self-energy-Potential difference, matrix size (i,j), where i and j belong to k. ebands=<ebands_t>=Datatype gathering info on the QP energies (KS if one shot) eig(Sigp%nbnds,Kmesh%nibz,Wfd%nsppol)=KS or QP energies for k-points, bands and spin occ(Sigp%nbnds,Kmesh%nibz,Wfd%nsppol)=occupation numbers, for each k point in IBZ, each band and spin
OUTPUT
Updated rdm_k matrix array with Go (Sigma_x - alpha*Sigma_x - Vxc) Go
SOURCE
209 subroutine calc_rdmx(ib1,ib2,ik_ibz,pot,rdm_k,ebands) 210 211 !Arguments ------------------------------------ 212 !scalars 213 integer,intent(in) :: ib1,ib2,ik_ibz 214 type(ebands_t),target,intent(in) :: ebands 215 !arrays 216 complex(dpc),intent(in) :: pot(:,:) 217 complex(dpc),intent(inout) :: rdm_k(:,:) 218 219 !Local variables ------------------------------ 220 !scalars 221 character(len=500) :: msg 222 integer :: ib1dm,ib2dm, units(2) 223 real(dp) :: spin_fact,tol8 224 !************************************************************************ 225 226 tol8=1.0e-8 227 spin_fact=two 228 units = [std_out, ab_out] 229 230 write(msg,'(a58,3f10.5)')' Computing the 1-RDM correction for Sx-Vxc and k-point: ',ebands%kptns(:,ik_ibz) 231 call wrtout(units, msg) 232 write(msg,'(a11,i5,a8,i5)')'from band ',ib1,' to band',ib2 233 call wrtout(units, msg) 234 235 rdm_k(:,:)=czero 236 do ib1dm=ib1,ib2-1 237 do ib2dm=ib1dm+1,ib2 238 if ((ebands%occ(ib1dm,ik_ibz,1)>tol8) .and. (ebands%occ(ib2dm,ik_ibz,1)<tol8)) then 239 rdm_k(1+(ib1dm-ib1),1+(ib2dm-ib1))=spin_fact& 240 &*pot(1+(ib1dm-ib1),1+(ib2dm-ib1))/(ebands%eig(ib1dm,ik_ibz,1)-ebands%eig(ib2dm,ik_ibz,1)+tol8) 241 ! Dji = Dij^* 242 rdm_k(1+(ib2dm-ib1),1+(ib1dm-ib1))=conjg(rdm_k(1+(ib1dm-ib1),1+(ib2dm-ib1))) 243 end if 244 end do 245 end do 246 247 end subroutine calc_rdmx
ABINIT/change_matrix [ Functions ]
NAME
change_matrix
FUNCTION
Transform integrals from KS -> NO and NO -> KS orbitals
Transform <NO_i|K[NO]|NO_j> -> <KS_i|K[NO]|KS_j>,
<KS_i|J[NO]|KS_j> -> <NO_i|J[NO]|NO_j>,
and <KS_i|T|KS_j> -> <NO_i|T|NO_j>
INPUTS
Kmesh <kmesh_t>=Structure describing the k-point sampling. Sigp<sigparams_t>=Parameters governing the self-energy calculation. nateigv = natural orbital eigenvectors nateigv(Wfd%mband,Wfd%mband,Wfd%nkibz,Sigp%nsppol))
OUTPUT
Mels %kinetic=matrix elements of $t$. %vhartr =matrix elements of $v_H$. Sr=sigma_t (see the definition of this structured datatype)
SOURCE
838 subroutine change_matrix(Sigp,Sr,Mels,Kmesh,nateigv) 839 840 !Arguments ------------------------------------ 841 !scalars 842 type(kmesh_t),intent(in) :: Kmesh 843 type(sigparams_t),intent(in) :: Sigp 844 type(sigma_t),intent(inout) :: Sr 845 type(melements_t),intent(inout) :: Mels 846 !arrays 847 complex(dpc),intent(in) :: nateigv(:,:,:,:) 848 !Local variables------------------------------- 849 !scalars 850 integer :: ikcalc,ik_ibz,ib1,ib2,ib1dm,ib2dm 851 !arrays 852 complex(dpc),allocatable :: mat2rot(:,:),Umat(:,:) 853 ! ************************************************************************* 854 855 do ikcalc=1,Sigp%nkptgw 856 ik_ibz=Kmesh%tab(Sigp%kptgw2bz(ikcalc)) ! Index of the irreducible k-point for GW 857 ib1=MINVAL(Sigp%minbnd(ikcalc,:)) ! min and max band indices for GW corrections (for this k-point) 858 ib2=MAXVAL(Sigp%maxbnd(ikcalc,:)) 859 ABI_MALLOC(mat2rot,(ib2-ib1+1,ib2-ib1+1)) 860 ABI_MALLOC(Umat,(ib2-ib1+1,ib2-ib1+1)) 861 ! <NO_i|K[NO]|NO_j> -> <KS_i|K[NO]|KS_j> 862 do ib1dm=1,ib2-ib1+1 863 do ib2dm=1,ib2-ib1+1 864 Umat(ib1dm,ib2dm)=nateigv(ib1+(ib1dm-1),ib1+(ib2dm-1),ik_ibz,1) 865 mat2rot(ib1dm,ib2dm)=Sr%x_mat(ib1+(ib1dm-1),ib1+(ib2dm-1),ik_ibz,1) 866 end do 867 end do 868 call rotate_ks_no(ib1,ib2,mat2rot,Umat,0) 869 do ib1dm=1,ib2-ib1+1 870 do ib2dm=1,ib2-ib1+1 871 Sr%x_mat(ib1+(ib1dm-1),ib1+(ib2dm-1),ik_ibz,1)=mat2rot(ib1dm,ib2dm) 872 end do 873 end do 874 ! <KS_i|J[NO]|KS_j> -> <NO_i|J[NO]|NO_j> 875 do ib1dm=1,ib2-ib1+1 876 do ib2dm=1,ib2-ib1+1 877 mat2rot(ib1dm,ib2dm)=Mels%vhartree(ib1+(ib1dm-1),ib1+(ib2dm-1),ik_ibz,1) 878 end do 879 end do 880 call rotate_ks_no(ib1,ib2,mat2rot,Umat,1) 881 do ib1dm=1,ib2-ib1+1 882 do ib2dm=1,ib2-ib1+1 883 Mels%vhartree(ib1+(ib1dm-1),ib1+(ib2dm-1),ik_ibz,1)=mat2rot(ib1dm,ib2dm) 884 end do 885 end do 886 ! <KS_i|T|KS_j> -> <NO_i|T|NO_j> 887 do ib1dm=1,ib2-ib1+1 888 do ib2dm=1,ib2-ib1+1 889 mat2rot(ib1dm,ib2dm)=Mels%kinetic(ib1+(ib1dm-1),ib1+(ib2dm-1),ik_ibz,1) 890 end do 891 end do 892 call rotate_ks_no(ib1,ib2,mat2rot,Umat,1) 893 do ib1dm=1,ib2-ib1+1 894 do ib2dm=1,ib2-ib1+1 895 Mels%kinetic(ib1+(ib1dm-1),ib1+(ib2dm-1),ik_ibz,1)=mat2rot(ib1dm,ib2dm) 896 end do 897 end do 898 ABI_FREE(Umat) 899 ABI_FREE(mat2rot) 900 end do 901 end subroutine change_matrix
ABINIT/get_chkprdm [ Functions ]
NAME
get_chkprdm
FUNCTION
Read all checkpoint files built on previous runs
INPUTS
Wfd<wfd_t>=Wave function descriptor see file 69_wfd/m_wfd.F90 Kmesh <kmesh_t>=Structure describing the k-point sampling. Sigp<sigparams_t>=Parameters governing the self-energy calculation. ebands=<ebands_t>=Datatype gathering info on the QP energies (KS if one shot) eig(Sigp%nbnds,Kmesh%nibz,Wfd%nsppol)=KS or QP energies for k-points, bands and spin occ(Sigp%nbnds,Kmesh%nibz,Wfd%nsppol)=occupation numbers, for each k point in IBZ, each band and spin occs = occ. numbers array occs(Wfd%mband,Wfd%nkibz) nateigv = natural orbital eigenvectors nateigv(Wfd%mband,Wfd%mband,Wfd%nkibz,Sigp%nsppol)) sigmak_todo = integer array initialized to 1 and its components are set to 0 if the kpoint is read from the checkpoint sigmak_todo(Wfd%nkibz) my_rank = rank of the mpi process.
OUTPUT
occ are updated if they are read from any checkpoint file nateigv are stored if they are read from any checkpoint file sigmak_todo components set to 1 if the kpoint is read from any checkpoint file
SOURCE
618 subroutine get_chkprdm(Wfd,Kmesh,Sigp,ebands,occs,nateigv,sigmak_todo,my_rank,gw1rdm_fname_in) 619 !Arguments ------------------------------------ 620 !scalars 621 integer,intent(in) :: my_rank 622 class(wfd_t),intent(in) :: Wfd 623 type(kmesh_t),intent(in) :: Kmesh 624 type(sigparams_t),intent(in) :: Sigp 625 type(ebands_t),intent(in) :: ebands 626 character(len=fnlen),intent(in) :: gw1rdm_fname_in 627 !arrays 628 integer,intent(inout) :: sigmak_todo(:) 629 real(dp),intent(inout) :: occs(:,:) 630 complex(dpc),intent(inout) :: nateigv(:,:,:,:) 631 !Local variables------------------------------- 632 !scalars 633 integer,parameter :: master=0,iunit=666314 634 integer :: ierr,ib1,ib2,ib3,ikcalc,istat,ik_ibz,ik_ibz_read,iread,iread_eigv 635 real(dp) :: auxl_read 636 character(len=fnlen) :: gw1rdm_fname 637 character(len=500) :: msg 638 !arrays 639 real(dp),allocatable :: occ_tmp(:),eigvect_tmp(:) 640 ! ************************************************************************* 641 642 if (my_rank==master) then 643 iread_eigv=Wfd%mband 644 iread_eigv=iread_eigv*(2*iread_eigv) 645 ABI_MALLOC(occ_tmp,(Wfd%mband)) 646 ABI_MALLOC(eigvect_tmp,(iread_eigv)) 647 648 do ikcalc=1,Sigp%nkptgw 649 ik_ibz=Kmesh%tab(Sigp%kptgw2bz(ikcalc)) ! Irred k-point for GW 650 if(ik_ibz<10) then 651 write(gw1rdm_fname,"(a,i1)") trim(gw1rdm_fname_in),ik_ibz 652 else if(ik_ibz<100 .and. ik_ibz>=10) then 653 write(gw1rdm_fname,"(a,i2)") trim(gw1rdm_fname_in),ik_ibz 654 else if(ik_ibz<1000 .and. ik_ibz>=100) then 655 write(gw1rdm_fname,"(a,i3)") trim(gw1rdm_fname_in),ik_ibz 656 else 657 ABI_ERROR("The maximum k-point label for the checkpoint file to read is 999.") 658 end if 659 write(msg,'(a1)')' ' 660 call wrtout(std_out,msg) 661 write(msg,'(a25,a)')' Reading checkpoint file ',gw1rdm_fname 662 call wrtout(std_out,msg) 663 write(msg,'(a1)')' ' 664 call wrtout(std_out,msg) 665 occ_tmp(:)=zero;eigvect_tmp(:)=zero; 666 open(unit=iunit,form='unformatted',file=gw1rdm_fname,iostat=istat,status='old') 667 iread=0;ik_ibz_read=0; 668 if (istat==0) then 669 do 670 if (iread<Wfd%mband) then 671 iread=iread+1 672 read(iunit,iostat=istat) auxl_read 673 if (istat==0) then 674 occ_tmp(iread)=auxl_read 675 end if 676 else if (iread<(iread_eigv+Wfd%mband)) then 677 iread=iread+1 678 read(iunit,iostat=istat) auxl_read 679 if (istat==0) then 680 eigvect_tmp(iread-Wfd%mband)=auxl_read 681 end if 682 else 683 read(iunit,iostat=istat) ik_ibz_read 684 if (istat==0 .and. ik_ibz_read/=0) then 685 iread=0 686 sigmak_todo(ik_ibz_read)=0 687 ib3=1 688 do ib1=1,Wfd%mband 689 occs(ib1,ik_ibz_read)=occ_tmp(ib1) 690 do ib2=1,Wfd%mband 691 nateigv(ib1,ib2,ik_ibz_read,1)=cmplx(eigvect_tmp(ib3),eigvect_tmp(ib3+1)) 692 ib3=ib3+2 693 end do 694 end do 695 ik_ibz_read=0 696 occ_tmp=zero;eigvect_tmp=zero; 697 end if 698 end if 699 if(istat/=0) then 700 exit 701 end if 702 end do 703 end if 704 close(iunit) 705 end do 706 write(msg,'(a1)')' ' 707 call wrtout(std_out,msg) 708 write(msg,'(a49)')' List of k-points read from all checkpoint files ' 709 call wrtout(std_out,msg) 710 do ikcalc=1,Sigp%nkptgw 711 ik_ibz=Kmesh%tab(Sigp%kptgw2bz(ikcalc)) ! Irred k-point for GW 712 if (sigmak_todo(ik_ibz)==0) then 713 write(msg,'(3f10.5)') ebands%kptns(1:,ik_ibz) 714 call wrtout(std_out,msg) 715 end if 716 enddo 717 write(msg,'(a1)')' ' 718 call wrtout(std_out,msg) 719 ABI_FREE(occ_tmp) 720 ABI_FREE(eigvect_tmp) 721 end if 722 723 ! Broadcast from master the information stored in occs and nateigv to all processes. 724 call xmpi_barrier(Wfd%comm) 725 ierr=0 726 call xmpi_bcast(sigmak_todo(:),master,Wfd%comm,ierr) 727 if(ierr/=0) then 728 ABI_ERROR("Error distributing the sigmak_todo table.") 729 endif 730 call xmpi_bcast(occs(:,:),master,Wfd%comm,ierr) 731 if(ierr/=0) then 732 ABI_ERROR("Error distributing the occs read from checkpoint file(s).") 733 endif 734 call xmpi_bcast(nateigv(:,:,:,:),master,Wfd%comm,ierr) 735 if(ierr/=0) then 736 ABI_ERROR("Error distributing the natural orbital eigenvectors read from checkpoint file(s).") 737 endif 738 739 end subroutine get_chkprdm
ABINIT/ks2no [ Functions ]
NAME
ks2no
FUNCTION
Transform the matrix mat from KS to NO basis
INPUTS
dim=dimension of the matrices mat=array in the KS basis rot=unitary matrix containg the eigenvectors in NO basis
OUTPUT
mat=array in the NO basis
SOURCE
1125 subroutine ks2no(ndim,mat,rot) 1126 !Arguments ------------------------------------ 1127 !scalars 1128 integer,intent(in) :: ndim 1129 !arrays 1130 complex(dpc),dimension(:,:),intent(in) :: rot 1131 complex(dpc),dimension(:,:),intent(inout) :: mat 1132 !Local variables ------------------------------ 1133 !scalars 1134 !arrays 1135 complex(dpc),allocatable :: res(:,:) 1136 !************************************************************************ 1137 1138 ABI_MALLOC(res,(ndim,ndim)) 1139 res=czero 1140 1141 ! <NO|Op|NO> = (U^t)* <KS|Op|KS> U 1142 res=matmul(conjg(transpose(rot)),mat) 1143 mat=matmul(res,rot) 1144 1145 ABI_FREE(res) 1146 1147 end subroutine ks2no
ABINIT/m_gwrdm [ Modules ]
NAME
m_gwrdm
FUNCTION
Compute density matrix correction Galitskii-Migdal Ecorr, G = Go + Go Sigma Go (imaginary freqs. are used in Sigma_c) and associated quantities (natural orbitals, matrix elements, etc.).
SOURCE
11 #if defined HAVE_CONFIG_H 12 #include "config.h" 13 #endif 14 15 #include "abi_common.h" 16 17 module m_gwrdm 18 19 use defs_basis 20 use m_gwdefs 21 use m_abicore 22 use m_xmpi 23 use m_errors 24 use m_hide_blas 25 use m_time 26 use m_wfd 27 use m_hdr 28 use m_dtset 29 30 use m_fstrings, only : sjoin, itoa 31 use m_melemts, only : melements_t 32 use m_bz_mesh, only : kmesh_t 33 use defs_datatypes, only : ebands_t 34 use m_sigma, only : sigma_t 35 use m_xctk, only : xcden 36 use m_gaussian_quadrature, only: cgqf 37 38 implicit none 39 40 private :: no2ks,ks2no,printrdm_k,rotate_ks_no
ABINIT/natoccs [ Functions ]
NAME
natoccs
FUNCTION
Calculate natural orbitals and occ. numbers for a given k-point
INPUTS
ib1=min band for given k ib2=max band for given k. ik_ibz= the label of k-point in the IBZ. iinfo=use Sigma_x or Sigma_c phaser weights=array containing the weights used in the quadrature. nateigv=array containing the natural eigenvectors in columns (nbands,nband,k-point,nspin) rdm_k=density matrix, matrix (i,j), where i and j belong to the k-point k (see m_sigma_driver for more details). occs = array containing the occ numbers for a given k-point occs(nband,k-point). ebands=<ebands_t>=Datatype gathering info on the QP energies (KS if one shot) eig(Sigp%nbnds,Kmesh%nibz,Wfd%nsppol)=KS or QP energies for k-points, bands and spin occ(Sigp%nbnds,Kmesh%nibz,Wfd%nsppol)=occupation numbers, for each k point in IBZ, each band and spin checksij=check the orthonormality of the nat. orbitals
OUTPUT
Compute the nat. orbitals and occ. numbers from the rdm_k matrix (for exchange and correlations)
SOURCE
342 subroutine natoccs(ib1,ib2,rdm_k,nateigv,occs,ebands,ik_ibz,iinfo,checksij) 343 !Arguments ------------------------------------ 344 !scalars 345 integer,intent(in) :: ib1,ib2,ik_ibz,iinfo 346 integer,intent(in),optional :: checksij 347 type(ebands_t),target,intent(in) :: ebands 348 !arrays 349 real(dp),intent(inout) :: occs(:,:) 350 complex(dpc),intent(inout) :: rdm_k(:,:),nateigv(:,:,:,:) 351 !Local variables ------------------------------ 352 !scalars 353 integer:: ndim,ib1dm,ib2dm,ib3dm,lwork,info 354 logical:: check_Sijmat 355 character(len=500) :: msg 356 real(dp) :: toccs_k,tol10 357 complex(dp) :: Sib1k_ib2k 358 !arrays 359 integer :: units(2) 360 real(dp),allocatable :: occs_tmp(:),occs_tmp2(:),rwork(:) 361 complex(dpc),allocatable :: work(:),tmp_mat(:,:),eigenvect(:,:) 362 !************************************************************************ 363 364 check_Sijmat=.false.; if (present(checksij)) check_Sijmat=.true. 365 units = [std_out, ab_out] 366 tol10=1.0e-10 367 368 ndim=ib2-ib1+1 369 lwork=2*ndim-1 370 ABI_MALLOC(occs_tmp,(ndim)) 371 ABI_MALLOC(occs_tmp2,(ndim)) 372 ABI_MALLOC(work,(lwork)) 373 ABI_MALLOC(tmp_mat,(ndim,ndim)) 374 ABI_MALLOC(eigenvect,(ndim,ndim)) 375 ABI_MALLOC(rwork,(3*ndim-2)) 376 377 tmp_mat=zero 378 do ib2dm=1,ndim 379 do ib1dm=ib2dm,ndim 380 tmp_mat(ib1dm,ib2dm)=rdm_k(ib1dm,ib2dm) 381 ! Dji = Dij^* 382 tmp_mat(ib2dm,ib1dm)=conjg(tmp_mat(ib1dm,ib2dm)) 383 end do 384 end do 385 386 work=zero 387 occs_tmp=zero 388 info=0 389 call zheev('v','u',ndim,tmp_mat,ndim,occs_tmp,work,lwork,rwork,info) 390 ABI_CHECK(info == 0, sjoin("Failed the diagonalization of the updated GW 1-RDM with info:", itoa(info))) 391 392 ! Sort in descending order 393 do ib1dm=1,ndim 394 occs_tmp2(ib1dm)=occs_tmp(ndim-(ib1dm-1)) 395 do ib2dm=1,ndim 396 eigenvect(ib2dm,ib1dm)=tmp_mat(ib2dm,(ndim-(ib1dm-1))) 397 end do 398 if (abs(occs_tmp2(ib1dm))<tol10) then 399 occs_tmp2(ib1dm)=zero 400 end if 401 end do 402 403 ! Check orthonormality? 404 if (check_Sijmat) then 405 do ib1dm=1,ndim 406 do ib2dm=1,ib1dm 407 Sib1k_ib2k=czero 408 do ib3dm=1,ndim 409 Sib1k_ib2k=Sib1k_ib2k+conjg(eigenvect(ib3dm,ib1dm))*eigenvect(ib3dm,ib2dm) 410 end do 411 if (ib1dm==ib2dm) then 412 if(abs(Sib1k_ib2k-cmplx(one,zero))>tol10) then 413 write(msg,'(a45,i5,a1,i5,f10.5)') 'Large deviation from identity for bands ',ib1dm,' ',ib2dm,real(Sib1k_ib2k) 414 call wrtout(std_out,msg) 415 endif 416 else 417 if (abs(Sib1k_ib2k)>tol10) then 418 write(msg,'(a45,i5,a1,i5,f10.5)') 'Large deviation from identity for bands ',ib1dm,' ',ib2dm,real(Sib1k_ib2k) 419 call wrtout(std_out,msg) 420 end if 421 end if 422 end do 423 end do 424 end if 425 426 ! Print results 427 if (info==0) then 428 if (iinfo==0) then 429 write(msg,'(a51,3f10.5)') 'Occs. after updating with Sx-Vxc corr. at k-point:',ebands%kptns(1:,ik_ibz) 430 else 431 write(msg,'(a51,3f10.5)') 'Occs. after updating with S_c correct. at k-point:',ebands%kptns(1:,ik_ibz) 432 endif 433 call wrtout(units, msg) 434 ib1dm=ndim-(ndim/10)*10 435 do ib2dm=1,(ndim/10)*10,10 436 write(msg,'(f11.5,9f10.5)') occs_tmp2(ib2dm:ib2dm+9) 437 call wrtout(units, msg) 438 end do 439 ib1dm=(ndim/10)*10+1 440 write(msg,'(f11.5,*(f10.5))') occs_tmp2(ib1dm:) 441 call wrtout(units, msg) 442 else 443 write(msg,'(a36,3f10.5)') 'Error computing occs. for k-point: ',ebands%kptns(1:,ik_ibz) 444 call wrtout(units, msg) 445 end if 446 447 ! Store natural orbital eigenvectors matrix and occs. Also compute total number of electrons for this k-point 448 toccs_k=zero 449 do ib1dm=1,ndim 450 do ib2dm=1,ndim 451 nateigv(ib1+(ib1dm-1),ib1+(ib2dm-1),ik_ibz,1)=eigenvect(ib1dm,ib2dm) 452 end do 453 occs(ib1+(ib1dm-1),ik_ibz)=occs_tmp2(ib1dm) ! Overwrite the initial KS-DFT occs from ib1 to ib2 454 toccs_k=toccs_k+occs_tmp2(ib1dm) 455 end do 456 457 write(msg,'(a22,i5,a3,i5,a21,f10.5)') ' Total occ. from band ',ib1,' to', ib2,' at current k-point: ',toccs_k 458 call wrtout(units, msg) 459 write(msg,'(a5)') ' ' 460 call wrtout(units, msg) 461 462 ABI_FREE(rwork) 463 ABI_FREE(work) 464 ABI_FREE(tmp_mat) 465 ABI_FREE(eigenvect) 466 ABI_FREE(occs_tmp) 467 ABI_FREE(occs_tmp2) 468 469 end subroutine natoccs
ABINIT/no2ks [ Functions ]
NAME
no2ks
FUNCTION
Transform the matrix mat from NO to KS basis
INPUTS
dim=dimension of the matrices mat=array in the KS basis rot=unitary matrix containg the eigenvectors in NO basis
OUTPUT
mat=array in the KS basis
SOURCE
1167 subroutine no2ks(ndim,mat,rot) 1168 !Arguments ------------------------------------ 1169 !scalars 1170 integer,intent(in) :: ndim 1171 !arrays 1172 complex(dpc),dimension(:,:),intent(in) :: rot 1173 complex(dpc),dimension(:,:),intent(inout) :: mat 1174 !Local variables ------------------------------ 1175 !scalars 1176 !arrays 1177 complex(dpc),allocatable :: res(:,:) 1178 !************************************************************************ 1179 1180 ABI_MALLOC(res,(ndim,ndim)) 1181 res=czero 1182 1183 ! <KS|Op|KS> = U <NO|Op|NO> (U^t)* 1184 res=matmul(rot,mat) 1185 mat=matmul(res,conjg(transpose(rot))) 1186 1187 ABI_FREE(res) 1188 1189 end subroutine no2ks
ABINIT/print_band_energies [ Functions ]
NAME
print_band_energies
FUNCTION
Print updated band energies
INPUTS
Kmesh <kmesh_t>=Structure describing the k-point sampling. Sigp<sigparams_t>=Parameters governing the self-energy calculation. Mels %kinetic=matrix elements of $t$. %vhartr =matrix elements of $v_H$. Sr=sigma_t (see the definition of this structured datatype) ebands=<ebands_t>=Datatype gathering info on the QP energies (KS if one shot) eig(Sigp%nbnds,Kmesh%nibz,Wfd%nsppol)=KS or QP energies for k-points, bands and spin occ(Sigp%nbnds,Kmesh%nibz,Wfd%nsppol)=occupation numbers, for each k point in IBZ, each band and spin
OUTPUT
SOURCE
1002 subroutine print_band_energies(b1gw,b2gw,Sr,Sigp,Mels,Kmesh,ebands,new_hartr,old_purex) 1003 !Arguments ------------------------------------ 1004 !scalars 1005 type(kmesh_t),intent(in) :: Kmesh 1006 type(sigparams_t),intent(in) :: Sigp 1007 type(sigma_t),intent(in) :: Sr 1008 type(ebands_t),intent(in) :: ebands 1009 type(melements_t),intent(in) :: Mels 1010 integer,intent(in) :: b1gw,b2gw 1011 !arrays 1012 complex(dpc),intent(in) :: old_purex(:,:),new_hartr(:,:) 1013 !Local variables------------------------------- 1014 !scalars 1015 integer :: ib,ikcalc,ik_ibz, units(2) 1016 real(dp) :: eik_new 1017 complex(dpc) :: delta_band_ibik 1018 character(len=500) :: msg 1019 !************************************************************************ 1020 1021 units = [std_out, ab_out] 1022 1023 write(msg,'(a1)') ' ' 1024 call wrtout(units, msg) 1025 write(msg,'(a42)') ' Computing band corrections Delta eik (eV)' 1026 call wrtout(units, msg) 1027 write(msg,'(a42)') ' -----------------------------------------' 1028 call wrtout(units, msg) 1029 write(msg,'(a1)') ' ' 1030 call wrtout(units, msg) 1031 write(msg,'(a1)') ' ' 1032 call wrtout(units, msg) 1033 write(msg,'(a110)') ' Band corrections Delta eik = <KS_i|K[NO]-a*K[KS]+vH[NO]& 1034 &-vH[KS]-Vxc[KS]|KS_i> and eik^new = eik^GS + Delta eik' 1035 call wrtout(units, msg) 1036 write(msg,'(a1)') ' ' 1037 call wrtout(units, msg) 1038 do ikcalc=1,Sigp%nkptgw 1039 ik_ibz=Kmesh%tab(Sigp%kptgw2bz(ikcalc)) ! Index of the irreducible k-point for GW 1040 write(msg,'(a127)')'---------------------------------------------------------& 1041 &--------------------------------------------------------------------' 1042 call wrtout(units, msg) 1043 write(msg,'(a)')' k-point band eik^GS eik^new Delta eik & 1044 & K[NO] a*K[KS] Vxc[KS] vH[NO] vH[KS]' 1045 call wrtout(units, msg) 1046 do ib=b1gw,b2gw 1047 delta_band_ibik=(new_hartr(ib,ikcalc)-Mels%vhartree(ib,ib,ik_ibz,1))& 1048 &+Sr%x_mat(ib,ib,ik_ibz,1)-Mels%vxcval(ib,ib,ik_ibz,1)-old_purex(ib,ikcalc) 1049 eik_new=real(ebands%eig(ib,ik_ibz,1))+real(delta_band_ibik) 1050 write(msg,'(i5,4x,i5,8(4x,f10.3))') & 1051 & ik_ibz,ib,real(ebands%eig(ib,ik_ibz,1))*Ha_eV,eik_new*Ha_eV,real(delta_band_ibik)*Ha_eV,& 1052 & real(Sr%x_mat(ib,ib,ik_ibz,1))*Ha_eV,real(old_purex(ib,ikcalc))*Ha_eV,& 1053 & real(Mels%vxcval(ib,ib,ik_ibz,1))*Ha_eV,& 1054 & real(new_hartr(ib,ikcalc))*Ha_eV,real(Mels%vhartree(ib,ib,ik_ibz,1))*Ha_eV 1055 call wrtout(units, msg) 1056 enddo 1057 enddo 1058 write(msg,'(a127)')'---------------------------------------------------------& 1059 &--------------------------------------------------------------------' 1060 call wrtout(units, msg) 1061 1062 end subroutine print_band_energies
ABINIT/print_chkprdm [ Functions ]
NAME
print_chkprdm
FUNCTION
Write the checkpoint file for a given k-point
INPUTS
Wfd<wfd_t>=Wave function descriptor see file 69_wfd/m_wfd.F90 occs = occ. numbers array occs(Wfd%mband,Wfd%nkibz) nateigv = natural orbital eigenvectors nateigv(Wfd%mband,Wfd%mband,Wfd%nkibz,Sigp%nsppol)) my_rank = rank of the mpi process. gw1rdm_fname_out = name of the gw1rdm checkpoint out without k-point extension
OUTPUT
SOURCE
760 subroutine print_chkprdm(Wfd,occs,nateigv,ik_ibz,my_rank,gw1rdm_fname_out) 761 !Arguments ------------------------------------ 762 !scalars 763 integer,intent(in) :: ik_ibz,my_rank 764 class(wfd_t),intent(in) :: Wfd 765 character(len=fnlen),intent(in) :: gw1rdm_fname_out 766 !arrays 767 real(dp),intent(in) :: occs(:,:) 768 complex(dpc),intent(in) :: nateigv(:,:,:,:) 769 !Local variables------------------------------- 770 !scalars 771 integer,parameter :: master=0,iunit=666314 772 integer :: iwrite,iwrite2 773 character(len=fnlen) :: gw1rdm_fname 774 character(len=500) :: msg 775 !arrays 776 ! ************************************************************************* 777 778 if (my_rank==master) then 779 if(ik_ibz<10) then 780 write(gw1rdm_fname,"(a,i1)") trim(gw1rdm_fname_out),ik_ibz 781 else if(ik_ibz<100 .and. ik_ibz>=10) then 782 write(gw1rdm_fname,"(a,i2)") trim(gw1rdm_fname_out),ik_ibz 783 else if(ik_ibz<1000 .and. ik_ibz>=100) then 784 write(gw1rdm_fname,"(a,i3)") trim(gw1rdm_fname_out),ik_ibz 785 else 786 ABI_ERROR("The maximum k-point label for the checkpoint file to write is 999.") 787 end if 788 write(msg,'(a1)')' ' 789 call wrtout(std_out,msg) 790 write(msg,'(a25,a)')' Writing checkpoint file ',gw1rdm_fname 791 call wrtout(std_out,msg) 792 write(msg,'(a1)')' ' 793 call wrtout(std_out,msg) 794 open(unit=iunit,form='unformatted',file=gw1rdm_fname) 795 do iwrite=1,Wfd%mband 796 write(iunit) occs(iwrite,ik_ibz) 797 end do 798 do iwrite=1,Wfd%mband 799 do iwrite2=1,Wfd%mband 800 write(iunit) real(nateigv(iwrite,iwrite2,ik_ibz,1)) 801 write(iunit) aimag(nateigv(iwrite,iwrite2,ik_ibz,1)) 802 end do 803 end do 804 write(iunit) ik_ibz 805 close(iunit) 806 end if 807 808 call xmpi_barrier(Wfd%comm) 809 810 end subroutine print_chkprdm
ABINIT/print_tot_occ [ Functions ]
NAME
print_tot_occ
FUNCTION
Compute and print the total (averaged) occ. from all k-points
INPUTS
ebands=<ebands_t>=Datatype gathering info on the QP energies (KS if one shot) eig(Sigp%nbnds,%nibz,Wfd%nsppol)=KS or QP energies for k-points, bands and spin occ(Sigp%nbnds,%nibz,Wfd%nsppol)=occupation numbers, for each k point in IBZ, each band and spin
OUTPUT
Print the total (averaged) occ. = sum_k weight_k * Nelec_k
SOURCE
555 subroutine print_tot_occ(ebands) 556 557 !Arguments ------------------------------------ 558 type(ebands_t),intent(in) :: ebands 559 560 !Local variables------------------------------- 561 !scalars 562 character(len=500) :: msg 563 integer :: ik,spin, units(2) 564 real(dp) :: wtk,occ_bks,tot_occ 565 ! ************************************************************************* 566 567 units = [std_out, ab_out] 568 569 tot_occ=zero 570 571 do spin=1,ebands%nsppol 572 do ik=1,ebands%nkpt 573 wtk = ebands%wtk(ik) 574 occ_bks = sum(ebands%occ(:,ik,spin)) 575 !if (sigma%nsig_ab==1) then ! Only closed-shell restricted is programed 576 tot_occ=tot_occ+occ_bks*wtk 577 !end if 578 end do 579 end do 580 581 write(msg,'(a1)') ' ' 582 call wrtout(units, msg) 583 write(msg,'(a39,f10.5)') ' Total averaged occ. from all k-points: ',tot_occ 584 call wrtout(units, msg) 585 write(msg,'(a1)') ' ' 586 call wrtout(units, msg) 587 588 end subroutine print_tot_occ
ABINIT/print_total_energy [ Functions ]
NAME
print_total_energy
FUNCTION
Print total energy and energy components
INPUTS
all energy terms are self-explanatory
OUTPUT
SOURCE
918 subroutine print_total_energy(ekin_energy,evext_energy,evextnl_energy,e_corepsp,eh_energy,ex_energy,& 919 exc_mbb_energy,e_ewald,etot,etot2,den_int) 920 !Arguments ------------------------------------ 921 !scalars 922 real(dp),intent(in) :: ekin_energy,evext_energy,evextnl_energy,e_corepsp,eh_energy,ex_energy 923 real(dp),intent(in) :: exc_mbb_energy,e_ewald,etot,etot2,den_int 924 925 !Local variables------------------------------- 926 character(len=500) :: msg 927 integer :: units(2) 928 929 !************************************************************************ 930 931 units = [std_out, ab_out] 932 933 write(msg,'(a1)')' ' 934 call wrtout(units, msg) 935 write(msg,'(a98)')'---------------------------------------------------------------& 936 &----------------------------------' 937 call wrtout(units, msg) 938 write(msg,'(a,f10.5,a,f10.3,a)')' Ekinetic = : ',ekin_energy,' Ha ,',ekin_energy*Ha_eV,' eV' 939 call wrtout(units, msg) 940 write(msg,'(a,f10.5,a,f10.3,a)')' Evext_l = : ',evext_energy,' Ha ,',evext_energy*Ha_eV,' eV' 941 call wrtout(units, msg) 942 write(msg,'(a,f10.5,a,f10.3,a)')' Evext_nl = : ',evextnl_energy,' Ha ,',evextnl_energy*Ha_eV,' eV' 943 call wrtout(units, msg) 944 write(msg,'(a,f10.5,a,f10.3,a)')' Epsp_core = : ',e_corepsp,' Ha ,',e_corepsp*Ha_eV,' eV' 945 call wrtout(units, msg) 946 write(msg,'(a,f10.5,a,f10.3,a)')' Ehartree = : ',eh_energy,' Ha ,',eh_energy*Ha_eV,' eV' 947 call wrtout(units, msg) 948 write(msg,'(a,f10.5,a,f10.3,a)')' Ex[SD] = : ',ex_energy,' Ha ,',ex_energy*Ha_eV,' eV' 949 call wrtout(units, msg) 950 write(msg,'(a,f10.5,a,f10.3,a)')' Exc[MBB] = : ',exc_mbb_energy,' Ha ,',exc_mbb_energy*Ha_eV,' eV' 951 call wrtout(units, msg) 952 write(msg,'(a,f10.5,a,f10.3,a)')' Enn = : ',e_ewald,' Ha ,',e_ewald*Ha_eV,' eV' 953 call wrtout(units, msg) 954 write(msg,'(a98)')'-----------------------------------------------------------------& 955 &--------------------------------' 956 call wrtout(units, msg) 957 write(msg,'(a,f10.5,a,f10.3,a)')' Etot[SD] = : ',etot,' Ha ,',etot*Ha_eV,' eV' 958 call wrtout(units, msg) 959 write(msg,'(a,f10.5,a,f10.3,a)')' Etot[MBB] = : ',etot2,' Ha ,',etot2*Ha_eV,' eV' 960 call wrtout(units, msg) 961 write(msg,'(a,f10.5,a,f10.3,a)')' Vee[SD] = : ',(ex_energy+eh_energy),' Ha ,',(ex_energy+eh_energy)*Ha_eV,' eV' 962 call wrtout(units, msg) 963 write(msg,'(a,f10.5,a,f10.3,a)')' Vee[MBB] = : ',(exc_mbb_energy+eh_energy),' Ha ,',& 964 &(exc_mbb_energy+eh_energy)*Ha_eV,' eV' 965 call wrtout(units, msg) 966 write(msg,'(a,f10.5)') ' Density = : ',den_int 967 call wrtout(units, msg) 968 write(msg,'(a)')' Vee[SD] (= Ehartree + Ex[SD]) energy obtained using GW 1-RDM:' 969 call wrtout(units, msg) 970 write(msg,'(a)')' Vee[MBB] (= Ehartree + Exc[MBB]) energy obtained using GW 1-RDM:' 971 call wrtout(units, msg) 972 write(msg,'(a98)')'-------------------------------------------------------------------& 973 &------------------------------' 974 call wrtout(units, msg) 975 976 end subroutine print_total_energy
ABINIT/printrdm_k [ Functions ]
NAME
printrdm_k
FUNCTION
Print the DM1 matrix
INPUTS
ib1=min band. ib2=max band. rdm_k=array containing the 1-RDM matrix
OUTPUT
Print the 1-RDM matrix
SOURCE
1208 subroutine printrdm_k(ib1,ib2,rdm_k) ! Only used for debug on this file, do not use it with large arrays! 1209 !Arguments ------------------------------------ 1210 !scalars 1211 integer,intent(in) :: ib1,ib2 1212 !arrays 1213 complex(dpc),intent(in) :: rdm_k(:,:) 1214 !Local variables ------------------------------ 1215 !scalars 1216 integer::ib1dm 1217 character(len=500) :: msg 1218 !arrays 1219 !************************************************************************ 1220 1221 do ib1dm=ib1,ib2 1222 write(msg,'(*(f12.5))') real(rdm_k(ib1dm,ib1:ib2)) 1223 call wrtout(std_out, msg) 1224 end do 1225 1226 end subroutine printrdm_k
ABINIT/quadrature_sigma_cw [ Functions ]
NAME
quadrature_sigma_cw
FUNCTION
Quadrature frequencies used for Sigma_c(iw) integration
INPUTS
Sigp<sigparams_t>=Parameters governing the self-energy calculation. Sr=sigma_t (see the definition of this structured datatype) weights=real quadrature weights.
OUTPUT
Update Sigp and Sr imaginary frequencies with iw, and weights with the quadrature weights
SOURCE
67 subroutine quadrature_sigma_cw(Sigp,Sr,weights) 68 !Arguments ------------------------------------ 69 !scalars 70 type(sigparams_t),intent(inout) :: Sigp 71 type(sigma_t),intent(inout) :: Sr 72 !arrays 73 real(dp),intent(inout) :: weights(:) 74 75 !Local variables ------------------------------ 76 !scalars 77 integer :: ifreqs,order_int,gaussian_kind,units(2) 78 real(dp) :: gwalpha,gwbeta,wmin,wmax 79 character(len=500) :: msg 80 !arrays 81 real(dp),allocatable :: freqs(:) 82 !************************************************************************ 83 84 units = [std_out, ab_out] 85 86 order_int=Sigp%nomegasi 87 write(msg,'(a45,i9)')' number of imaginary frequencies for Sigma_c ',order_int 88 call wrtout(units, msg) 89 write(msg,'(a1)')' ' 90 call wrtout(units, msg) 91 order_int=Sigp%nomegasi 92 ABI_MALLOC(freqs,(order_int)) 93 gaussian_kind=1 94 gwalpha=zero 95 gwbeta=zero 96 wmin=zero 97 wmax=one 98 call cgqf(order_int,gaussian_kind,gwalpha,gwbeta,wmin,wmax,freqs,weights) 99 ! From 0 to 1 -> 0 to infinity 100 weights(:)=weights(:)/(one-freqs(:))**two 101 freqs(:)=freqs(:)/(one-freqs(:)) 102 ! Form complex frequencies from 0 to iInf and print them in the log file 103 write(msg,'(a52)')' Re(iw) Im(iw) Weight ' 104 call wrtout(std_out,msg) 105 write(msg,'(a52)')' -------- -------- -------- ' 106 call wrtout(std_out,msg) 107 do ifreqs=1,order_int 108 Sigp%omegasi(ifreqs)=cmplx(zero,freqs(ifreqs)) 109 Sr%omega_i(ifreqs)=Sigp%omegasi(ifreqs) 110 write(msg,'(3f17.5)') Sr%omega_i(ifreqs),weights(ifreqs) 111 call wrtout(std_out,msg) 112 enddo 113 ABI_FREE(freqs) 114 115 end subroutine quadrature_sigma_cw
ABINIT/rotate_ks_no [ Functions ]
NAME
rotate_ks_no
FUNCTION
Rotate a matrix from KS to NO basis and vicerversa.
INPUTS
ib1=min band for given k ib2=max band for given k. Umat=array containing the eigenvectors in Columns (a unitary matrix) Mat=initially an array containing the matrix elements in KS or NO basis option=0 rotate from NO -> KS | 1 rotate from KS -> NO
OUTPUT
Rotate a matrix from KS to NO basis and vicerversa and save the new matrix on Mat. Mat=at the end an array containing the matrix elements in NO or KS basis
SOURCE
1085 subroutine rotate_ks_no(ib1,ib2,Mat,Umat,option) 1086 !Arguments ------------------------------------ 1087 !scalars 1088 integer,intent(in) :: ib1,ib2,option 1089 !arrays 1090 complex(dpc),intent(in) :: Umat(:,:) 1091 complex(dpc),intent(inout) :: Mat(:,:) 1092 !Local variables ------------------------------ 1093 !scalars 1094 integer:: ndim 1095 !arrays 1096 !************************************************************************ 1097 1098 ndim=ib2-ib1+1 1099 if (option==0) then 1100 call no2ks(ndim,Mat,Umat) 1101 else 1102 call ks2no(ndim,Mat,Umat) 1103 end if 1104 1105 end subroutine rotate_ks_no
ABINIT/update_hdr_bst [ Functions ]
NAME
update_hdr_bst
FUNCTION
Update the Hdr for the WFK and DEN files and the occ. numbers in the ebands file for a given k-point
INPUTS
Wfd<wfd_t>=Datatype gathering data on QP amplitudes. ngfft_in(18)=information on the fine FFT grid used for densities and potentials. b1gw=min band for given k in the interval where we update. b2gw=max band for given k in the interval where we update. occs= array containing the occ numbers for a given k-point occs_ks(nband,k-point). ebands=<ebands_t>=Datatype gathering info on the QP energies (KS if one shot) eig(Sigp%nbnds,Kmesh%nibz,Wfd%nsppol)=KS or QP energies for k-points, bands and spin occ(Sigp%nbnds,Kmesh%nibz,Wfd%nsppol)=occupation numbers, for each k point in IBZ, each band and spin
OUTPUT
Updated Hdr and ebands information
SOURCE
494 subroutine update_hdr_bst(Wfd,occs,b1gw,b2gw,ebands,Hdr,ngfft_in) 495 496 !Arguments ------------------------------------ 497 !scalars 498 integer,intent(in) :: b1gw,b2gw 499 integer,intent(in),dimension(3) :: ngfft_in 500 type(ebands_t),target,intent(inout) :: ebands 501 type(Hdr_type),intent(inout) :: Hdr 502 class(wfd_t),intent(in) :: Wfd 503 !arrays 504 real(dp),intent(in) :: occs(:,:) 505 !Local variables ------------------------------ 506 !scalars 507 integer :: ib1dm,ib2dm,dim_bands,ikpoint 508 !arrays 509 !************************************************************************ 510 511 ! ebands occ (QP_ebands ones) are changed and never recoverd 512 do ikpoint=1,ebands%nkpt 513 ebands%occ(b1gw:b2gw,ikpoint,1) = occs(b1gw:b2gw,ikpoint) ! Spins summed, occ in [0:2] 514 enddo 515 ABI_COMMENT("QP_ebands: occupancies were updated with nat. orb. ones") 516 if ((size(Hdr%occ(:))/ebands%nkpt) < (b2gw-b1gw+1)) then 517 !Actually, we should never reach this point because the code should stop during Wfd initialization in m_sigma_driver 518 ABI_ERROR("Impossible to use the existing read WFK to build a new one!") 519 end if 520 521 ! Update occ in Hdr before printing 522 ib1dm=1 523 do ikpoint=1,ebands%nkpt 524 dim_bands=size(ebands%occ(:,ikpoint,1)) 525 do ib2dm=1,dim_bands 526 Hdr%occ(ib1dm)=ebands%occ(ib2dm,ikpoint,1) ! Because Hdr%occ is a 1-D array 527 ib1dm=ib1dm+1 528 end do 529 end do 530 531 Hdr%npwarr(:)=Wfd%npwarr(:) ! Use the npw and ngfft = ones used in GW calc 532 Hdr%ngfft(1:3)=ngfft_in(1:3) 533 ABI_COMMENT("Hdr_sigma: occupancies, npw, and ngfft were updated") 534 535 end subroutine update_hdr_bst