TABLE OF CONTENTS
- ABINIT/m_sigma
- ABINIT/mels_get_haene
- ABINIT/mels_get_kiene
- m_sigma/find_wpoles_for_cd
- m_sigma/gw_spectral_function
- m_sigma/print_Sigma_perturbative
- m_sigma/print_Sigma_QPSC
- m_sigma/sigma_distribute_bks
- m_sigma/sigma_free
- m_sigma/sigma_get_excene
- m_sigma/sigma_get_exene
- m_sigma/sigma_init
- m_sigma/sigma_ncwrite
- m_sigma/sigma_t
- m_sigma/write_sigma_header
- m_sigma/write_sigma_results
ABINIT/m_sigma [ Modules ]
NAME
m_sigma
FUNCTION
This module provides the definition of the sigma_t data type used to store results of the GW calculation as well as as methods bound to the object.
COPYRIGHT
Copyright (C) 2008-2024 ABINIT group (MG, FB, GMR, VO, LR, RWG) 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_sigma 25 26 use defs_basis 27 use m_xmpi 28 use m_abicore 29 use m_errors 30 use, intrinsic :: iso_c_binding 31 use m_nctk 32 use m_yaml 33 use m_melemts 34 #ifdef HAVE_NETCDF 35 use netcdf 36 #endif 37 use m_wfd 38 39 use defs_datatypes, only : ebands_t 40 use defs_abitypes, only : MPI_type 41 use m_numeric_tools, only : c2r 42 use m_gwdefs, only : unt_gw, unt_sig, unt_sgr, unt_sgm, unt_gwdiag, sigparams_t, sigma_needs_w, unt_sigc ! MRM 43 use m_crystal, only : crystal_t 44 use m_bz_mesh, only : kmesh_t, littlegroup_t, findqg0 45 use m_screening, only : epsilonm1_results 46 use m_stream_string, only : stream_string 47 48 implicit none 49 50 private
ABINIT/mels_get_haene [ Functions ]
NAME
mels_get_haene
FUNCTION
Compute the Hartree energy
INPUTS
Kmesh <kmesh_t>=Structure describing the k-point sampling. Sr=sigma_t (see the definition of this structured datatype) bands=<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 Mels %vhartr=matrix elements of $v_H$.
OUTPUT
Compute the Hartree energy on eh_energy
SOURCE
1195 real(dp) pure function mels_get_haene(sigma,Mels,kmesh,bands) result(eh_energy) 1196 1197 !Arguments ------------------------------------ 1198 !scalars 1199 class(sigma_t),intent(in) :: sigma 1200 type(kmesh_t),intent(in) :: kmesh 1201 type(ebands_t),intent(in) :: bands 1202 type(melements_t),intent(in) :: Mels 1203 1204 !Local variables------------------------------- 1205 !scalars 1206 integer :: ik,ib,spin 1207 real(dp) :: wtk,occ_bks 1208 1209 ! ************************************************************************* 1210 1211 eh_energy=zero 1212 1213 do spin=1,sigma%nsppol 1214 do ik=1,sigma%nkibz 1215 wtk = kmesh%wt(ik) 1216 do ib=sigma%b1gw,sigma%b2gw 1217 occ_bks = bands%occ(ib,ik,spin) 1218 if (sigma%nsig_ab==1) then ! Only closed-shell restricted is programed 1219 eh_energy=eh_energy+occ_bks*wtk*Mels%vhartree(ib,ib,ik,spin) 1220 end if 1221 end do 1222 end do 1223 end do 1224 1225 eh_energy=half*eh_energy 1226 1227 end function mels_get_haene
ABINIT/mels_get_kiene [ Functions ]
NAME
mels_get_kiene
FUNCTION
Compute the kinetic energy
INPUTS
Kmesh <kmesh_t>=Structure describing the k-point sampling. Sr=sigma_t (see the definition of this structured datatype) bands=<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 Mels %kinetic=matrix elements of $T$.
OUTPUT
Compute the kinetic energy on ek_energy
SOURCE
1253 real(dp) pure function mels_get_kiene(sigma,Mels,kmesh,bands) result(ek_energy) 1254 1255 !Arguments ------------------------------------ 1256 !scalars 1257 type(sigma_t),intent(in) :: sigma 1258 type(kmesh_t),intent(in) :: kmesh 1259 type(ebands_t),intent(in) :: bands 1260 type(melements_t),intent(in) :: Mels 1261 1262 !Local variables------------------------------- 1263 !scalars 1264 integer :: ik,ib,spin 1265 real(dp) :: wtk,occ_bks 1266 1267 ! ************************************************************************* 1268 1269 ek_energy=zero 1270 1271 do spin=1,sigma%nsppol 1272 do ik=1,sigma%nkibz 1273 wtk = kmesh%wt(ik) 1274 do ib=sigma%b1gw,sigma%b2gw 1275 occ_bks = bands%occ(ib,ik,spin) 1276 if (sigma%nsig_ab==1) then ! Only closed-shell restricted is programed 1277 ek_energy=ek_energy+occ_bks*wtk*Mels%kinetic(ib,ib,ik,spin) 1278 end if 1279 end do 1280 end do 1281 end do 1282 1283 ek_energy=ek_energy 1284 1285 end function mels_get_kiene
m_sigma/find_wpoles_for_cd [ Functions ]
[ Top ] [ m_sigma ] [ Functions ]
NAME
find_wpoles_for_cd
FUNCTION
Find the max frequency needed to account for all the poles of GW used in the contour deformation technique.
INPUTS
Sigp=sigparams_t
OUTPUT
omega_max
NOTES
SOURCE
1308 subroutine find_wpoles_for_cd(Sigp,Sr,Kmesh,BSt,omega_max) 1309 1310 !Arguments ------------------------------------ 1311 !scalars 1312 type(sigparams_t),intent(in) :: Sigp 1313 type(sigma_t),intent(in) :: Sr 1314 type(ebands_t),intent(in) :: Bst 1315 type(kmesh_t),intent(in) :: Kmesh 1316 real(dp),intent(out) :: omega_max 1317 1318 !Local variables------------------------------- 1319 !scalars 1320 integer :: spin,ik_ibz,band_gr,bgw_start,bgw_stop,io,ioe0j 1321 integer :: ikgw,ikgw_ibz,ikgw_bz,band_gw,nomega_tot 1322 real(dp) :: e_green,e_screen,theta_mu_minus_e0i,e_qp 1323 real(dp) :: fact_sp 1324 !character(len=500) :: msg 1325 !arrays 1326 real(dp),allocatable :: omegame0i(:) 1327 1328 ! ************************************************************************* 1329 1330 omega_max = smallest_real 1331 ! 1332 ! === Normalization of theta_mu_minus_e0i === 1333 ! * If nsppol==2, qp_occ $\in [0,1]$ 1334 fact_sp=one 1335 if (Bst%nsppol==1) then 1336 fact_sp=half; if (Bst%nspinor==2) fact_sp=one 1337 end if 1338 ! 1339 ! Total number of frequencies for sigma (Spectral function + mesh for the derivative). 1340 nomega_tot=Sr%nomega_r+Sr%nomega4sd 1341 ABI_MALLOC(omegame0i,(nomega_tot)) 1342 1343 ioe0j=Sr%nomega4sd/2+1 1344 ! 1345 ! Loop over bands used to construct the Green function. 1346 do spin=1,Bst%nsppol 1347 do ik_ibz=1,Bst%nkpt 1348 do band_gr=1,Bst%nband(ik_ibz+(spin-1)*Bst%nkpt) 1349 e_green = Bst%eig(band_gr,ik_ibz,spin) 1350 theta_mu_minus_e0i= Bst%occ(band_gr,ik_ibz,spin)*fact_sp 1351 ! 1352 ! Loop over GW states. 1353 do ikgw=1,Sigp%nkptgw 1354 bgw_start=Sigp%minbnd(ikgw,spin) 1355 bgw_stop =Sigp%minbnd(ikgw,spin) 1356 ikgw_bz =Sigp%kptgw2bz(ikgw_bz) 1357 ikgw_ibz =Kmesh%tab(ikgw_bz) 1358 1359 do band_gw=bgw_start,bgw_stop 1360 e_qp = Bst%eig(band_gw,ikgw_ibz,spin) 1361 ! 1362 ! Get frequencies $\omega$-\epsilon_in$ to evaluate $d\Sigma/dE$, note the spin 1363 ! subtract e_KS since we have stored e_KS+ Delta \omega in Sr%omega4sd, not required for AC 1364 if (Sr%nomega_r>0) omegame0i(1:Sr%nomega_r)=DBLE(Sigp%omega_r(1:Sr%nomega_r))-e_green 1365 do io=Sr%nomega_r+1,nomega_tot 1366 !omegame0i(io)=DBLE(Sr%omega4sd(band_gw,ikgw_ibz,io-Sr%nomega_r,spin)) - e_green 1367 !Sr%omega4sd(jb,ik_ibz,io,spin)=Sr%egw(jb,ik_ibz,spin)+Sigp%deltae*(io-ioe0j) 1368 omegame0i(io) = e_qp + Sigp%deltae*(io-ioe0j) - e_green 1369 end do 1370 1371 do io=1,nomega_tot 1372 e_screen = ABS(omegame0i(io)) 1373 if (omegame0i(io)>tol12) then 1374 !ket(spadc+ig,ios)=ket(spadc+ig,ios)+ct*(one-theta_mu_minus_e0i) 1375 if ( (one-theta_mu_minus_e0i) > tol12 ) omega_max = MAX(omega_max, e_screen) 1376 end if 1377 if (omegame0i(io)<-tol12) then 1378 !ket(spadc+ig,ios)=ket(spadc+ig,ios)-ct*theta_mu_minus_e0i 1379 if ( theta_mu_minus_e0i > tol12) omega_max = MAX(omega_max, e_screen) 1380 end if 1381 end do 1382 1383 end do 1384 end do 1385 ! 1386 end do 1387 end do 1388 end do 1389 1390 ABI_FREE(omegame0i) 1391 1392 end subroutine find_wpoles_for_cd
m_sigma/gw_spectral_function [ Functions ]
[ Top ] [ m_sigma ] [ Functions ]
NAME
gw_spectral_function
FUNCTION
Compute the spectral function
INPUTS
io,ib,ikibz,spin=Frequency, band, k-point, spin index Sr=sigma results datatype
OUTPUT
SOURCE
594 real(dp) pure function gw_spectral_function(Sr, io, ib, ikibz, spin) result(aw) 595 596 !Arguments ------------------------------------ 597 integer,intent(in) :: io,ib,ikibz,spin 598 type(sigma_t),intent(in) :: Sr 599 600 ! ********************************************************************* 601 602 aw = one / pi * abs(aimag(Sr%sigcme(ib,ikibz,io,spin))) & 603 /( (real(Sr%omega_r(io) - Sr%hhartree(ib,ib,ikibz,spin) - Sr%sigxcme(ib,ikibz,io,spin)))**2 & 604 +(aimag(Sr%sigcme(ib,ikibz,io,spin))) ** 2) / Ha_eV 605 606 end function gw_spectral_function
m_sigma/print_Sigma_perturbative [ Functions ]
[ Top ] [ m_sigma ] [ Functions ]
NAME
print_Sigma_perturbative
FUNCTION
Write the results of the GW calculation done with the perturbative approach
INPUTS
OUTPUT
SOURCE
624 subroutine print_Sigma_perturbative(Sr,ik_ibz,iband,isp,unit,prtvol,mode_paral,witheader,ydoc) 625 626 !Arguments ------------------------------------ 627 !scalars 628 integer,intent(in) :: iband,ik_ibz,isp 629 integer,optional,intent(in) :: prtvol,unit 630 character(len=*),optional,intent(in) :: mode_paral 631 logical,optional,intent(in) :: witheader 632 type(sigma_t),intent(in) :: Sr 633 type(yamldoc_t),intent(inout),optional :: ydoc 634 635 !Local variables------------------------------- 636 !scalars 637 integer :: my_unt,verbose 638 character(len=4) :: my_mode 639 character(len=500) :: msg 640 641 ! ********************************************************************* 642 643 my_unt =std_out; if (PRESENT(unit )) my_unt =unit 644 verbose=0 ; if (PRESENT(prtvol )) verbose=prtvol 645 my_mode='COLL' ; if (PRESENT(mode_paral)) my_mode=mode_paral 646 647 if (PRESENT(witheader)) then 648 if (witheader) then 649 call wrtout(my_unt,' Band E0 <VxcDFT> SigX SigC(E0) Z dSigC/dE Sig(E) E-E0 E ',my_mode) 650 end if 651 end if 652 653 if (Sr%usepawu==0) then 654 655 if (Sr%nsig_ab/=1) then 656 write(msg,'(i5,9f8.3)') & 657 iband, & 658 Sr%e0 (iband,ik_ibz,1)*Ha_eV, & 659 SUM(Sr%vxcme (iband,ik_ibz,:))*Ha_eV, & 660 SUM(Sr%sigxme (iband,ik_ibz,:))*Ha_eV, & 661 REAL(SUM(Sr%sigcmee0(iband,ik_ibz,:)))*Ha_eV,& 662 REAL(Sr%ze0 (iband,ik_ibz,1)), & 663 REAL(SUM(Sr%dsigmee0(iband,ik_ibz,:))), & 664 REAL(SUM(Sr%sigmee (iband,ik_ibz,:)))*Ha_eV,& 665 REAL(Sr%degw (iband,ik_ibz,1))*Ha_eV, & 666 REAL(Sr%egw (iband,ik_ibz,1))*Ha_eV 667 call wrtout(my_unt,msg,my_mode) 668 if (present(ydoc)) call ydoc%add_tabular_line(msg) 669 if (verbose/=0) then 670 write(msg,'(i5,9f8.3)') & 671 iband, & 672 zero, & 673 zero, & 674 zero, & 675 AIMAG(SUM(Sr%sigcmee0(iband,ik_ibz,:)))*Ha_eV,& 676 AIMAG(Sr%ze0 (iband,ik_ibz,1)), & 677 AIMAG(SUM(Sr%dsigmee0(iband,ik_ibz,:))), & 678 AIMAG(SUM(Sr%sigmee (iband,ik_ibz,:)))*Ha_eV,& 679 AIMAG(Sr%degw (iband,ik_ibz,1))*Ha_eV, & 680 AIMAG(Sr%egw (iband,ik_ibz,1))*Ha_eV 681 call wrtout(my_unt,msg,my_mode) 682 if(present(ydoc)) call ydoc%add_tabular_line(msg) 683 end if 684 else 685 write(msg,'(i5,9f8.3)') & 686 iband, & 687 Sr%e0 (iband,ik_ibz,isp)*Ha_eV, & 688 Sr%vxcme (iband,ik_ibz,isp)*Ha_eV, & 689 Sr%sigxme (iband,ik_ibz,isp)*Ha_eV, & 690 REAL(Sr%sigcmee0(iband,ik_ibz,isp))*Ha_eV,& 691 REAL(Sr%ze0 (iband,ik_ibz,isp)), & 692 REAL(Sr%dsigmee0(iband,ik_ibz,isp)), & 693 REAL(Sr%sigmee (iband,ik_ibz,isp))*Ha_eV,& 694 REAL(Sr%degw (iband,ik_ibz,isp))*Ha_eV,& 695 REAL(Sr%egw (iband,ik_ibz,isp))*Ha_eV 696 call wrtout(my_unt,msg,my_mode) 697 if (present(ydoc)) call ydoc%add_tabular_line(msg) 698 699 if (verbose/=0) then 700 write(msg,'(i5,9f8.3)') & 701 iband, & 702 zero, & 703 zero, & 704 zero, & 705 AIMAG(Sr%sigcmee0(iband,ik_ibz,isp))*Ha_eV,& 706 AIMAG(Sr%ze0 (iband,ik_ibz,isp)), & 707 AIMAG(Sr%dsigmee0(iband,ik_ibz,isp)), & 708 AIMAG(Sr%sigmee (iband,ik_ibz,isp))*Ha_eV,& 709 AIMAG(Sr%degw (iband,ik_ibz,isp))*Ha_eV,& 710 AIMAG(Sr%egw (iband,ik_ibz,isp))*Ha_eV 711 call wrtout(my_unt,msg,my_mode) 712 if (present(ydoc)) call ydoc%add_tabular_line(msg) 713 end if 714 end if 715 716 else 717 ! PAW+U+GW calculation. 718 ABI_CHECK(Sr%nsig_ab==1,'DFT+U with spinor not implemented') 719 write(msg,'(i5,10f8.3)') & 720 iband, & 721 Sr%e0 (iband,ik_ibz,isp)*Ha_eV, & 722 Sr%vxcme (iband,ik_ibz,isp)*Ha_eV, & 723 Sr%vUme (iband,ik_ibz,isp)*Ha_eV, & 724 Sr%sigxme (iband,ik_ibz,isp)*Ha_eV, & 725 REAL(Sr%sigcmee0(iband,ik_ibz,isp))*Ha_eV,& 726 REAL(Sr%ze0 (iband,ik_ibz,isp)), & 727 REAL(Sr%dsigmee0(iband,ik_ibz,isp)), & 728 REAL(Sr%sigmee (iband,ik_ibz,isp))*Ha_eV,& 729 REAL(Sr%degw (iband,ik_ibz,isp))*Ha_eV,& 730 REAL(Sr%egw (iband,ik_ibz,isp))*Ha_eV 731 call wrtout(my_unt,msg,my_mode) 732 if(present(ydoc)) call ydoc%add_tabular_line(msg) 733 734 if (verbose/=0) then 735 write(msg,'(i5,10f8.3)') & 736 iband, & 737 zero, & 738 zero, & 739 zero, & 740 zero, & 741 AIMAG(Sr%sigcmee0(iband,ik_ibz,isp))*Ha_eV,& 742 AIMAG(Sr%ze0 (iband,ik_ibz,isp)), & 743 AIMAG(Sr%dsigmee0(iband,ik_ibz,isp)), & 744 AIMAG(Sr%sigmee (iband,ik_ibz,isp))*Ha_eV,& 745 AIMAG(Sr%degw (iband,ik_ibz,isp))*Ha_eV,& 746 AIMAG(Sr%egw (iband,ik_ibz,isp))*Ha_eV 747 call wrtout(my_unt,msg,my_mode) 748 if(present(ydoc)) call ydoc%add_tabular_line(msg) 749 end if 750 end if 751 752 end subroutine print_Sigma_perturbative
m_sigma/print_Sigma_QPSC [ Functions ]
[ Top ] [ m_sigma ] [ Functions ]
NAME
print_Sigma_QPSC
FUNCTION
Write the results of the GW calculation in case of self-consistency
INPUTS
OUTPUT
SOURCE
770 subroutine print_Sigma_QPSC(Sr,ik_ibz,iband,isp,KS_BSt,unit,prtvol,mode_paral,ydoc) 771 772 !Arguments ------------------------------------ 773 !scalars 774 integer,intent(in) :: iband,ik_ibz,isp 775 integer,intent(in),optional :: prtvol,unit 776 character(len=*),intent(in),optional :: mode_paral 777 type(sigma_t),intent(in) :: Sr 778 type(ebands_t),intent(in) :: KS_BSt 779 type(yamldoc_t),intent(inout),optional :: ydoc 780 781 !Local variables------------------------------- 782 !scalars 783 integer :: my_unt,verbose 784 character(len=4) :: my_mode 785 character(len=500) :: msg 786 787 ! ********************************************************************* 788 789 my_unt =std_out; if (PRESENT(unit )) my_unt =unit 790 verbose=0 ; if (PRESENT(prtvol )) verbose=prtvol 791 my_mode='COLL' ; if (PRESENT(mode_paral)) my_mode=mode_paral 792 793 ! write(msg,'(a)')& 794 !& ' Band E_DFT <VxcDFT> E(N-1) <Hhartree> SigX SigC[E(N-1)]',& 795 !& ' Z dSigC/dE Sig[E(N)] DeltaE E(N)_pert E(N)_diago' 796 797 if (Sr%usepawu==0 .or. .TRUE.) then 798 if (Sr%nsig_ab/=1) then 799 write(msg,'(i5,12(2x,f8.3))') & 800 iband, & 801 KS_BSt%eig (iband,ik_ibz,1)*Ha_eV, & 802 SUM(Sr%vxcme (iband,ik_ibz,:))*Ha_eV, & 803 Sr%e0 (iband,ik_ibz,1)*Ha_eV, & 804 REAL(SUM(Sr%hhartree(iband,iband,ik_ibz,:)))*Ha_eV,& 805 SUM(Sr%sigxme (iband,ik_ibz,:))*Ha_eV, & 806 REAL(SUM(Sr%sigcmee0(iband,ik_ibz,:)))*Ha_eV, & 807 REAL(Sr%ze0 (iband,ik_ibz,1)), & 808 REAL(SUM(Sr%dsigmee0(iband,ik_ibz,:))), & 809 REAL(SUM(Sr%sigmee (iband,ik_ibz,:)))*Ha_eV, & 810 REAL(Sr%degw (iband,ik_ibz,1))*Ha_eV, & 811 REAL(Sr%egw (iband,ik_ibz,1))*Ha_eV, & 812 Sr%en_qp_diago (iband,ik_ibz,1)*Ha_eV 813 call wrtout(my_unt,msg,my_mode) 814 if (present(ydoc)) call ydoc%add_tabular_line(msg) 815 816 write(msg,'(i5,12(2x,f8.3))') & 817 iband, & 818 zero, & 819 zero, & 820 zero, & 821 AIMAG(SUM(Sr%hhartree(iband,iband,ik_ibz,:)))*Ha_eV,& 822 zero, & 823 AIMAG(SUM(Sr%sigcmee0(iband,ik_ibz,:)))*Ha_eV, & 824 AIMAG(Sr%ze0 (iband,ik_ibz,1)), & 825 AIMAG(SUM(Sr%dsigmee0(iband,ik_ibz,:))), & 826 AIMAG(SUM(Sr%sigmee (iband,ik_ibz,:)))*Ha_eV, & 827 AIMAG(Sr%degw (iband,ik_ibz,1))*Ha_eV, & 828 AIMAG(Sr%egw (iband,ik_ibz,1))*Ha_eV, & 829 zero 830 if (verbose/=0) then 831 call wrtout(my_unt,msg,my_mode) 832 if( present(ydoc)) call ydoc%add_tabular_line(msg) 833 end if 834 else 835 write(msg,'(i5,12(2x,f8.3))') & 836 iband, & 837 KS_BSt%eig (iband,ik_ibz,isp)*Ha_eV, & 838 Sr%vxcme (iband,ik_ibz,isp)*Ha_eV, & 839 Sr%e0 (iband,ik_ibz,isp)*Ha_eV, & 840 REAL(Sr%hhartree (iband,iband,ik_ibz,isp))*Ha_eV,& 841 Sr%sigxme (iband,ik_ibz,isp)*Ha_eV, & 842 REAL(Sr%sigcmee0 (iband,ik_ibz,isp))*Ha_eV, & 843 REAL(Sr%ze0 (iband,ik_ibz,isp)), & 844 REAL(Sr%dsigmee0 (iband,ik_ibz,isp)), & 845 REAL(Sr%sigmee (iband,ik_ibz,isp))*Ha_eV, & 846 REAL(Sr%degw (iband,ik_ibz,isp))*Ha_eV, & 847 REAL(Sr%egw (iband,ik_ibz,isp))*Ha_eV, & 848 Sr%en_qp_diago(iband,ik_ibz,isp)*Ha_eV 849 call wrtout(my_unt,msg,my_mode) 850 if (present(ydoc)) call ydoc%add_tabular_line(msg) 851 852 write(msg,'(i5,12(2x,f8.3))') & 853 iband, & 854 zero, & 855 zero, & 856 zero, & 857 AIMAG(Sr%hhartree (iband,iband,ik_ibz,isp))*Ha_eV,& 858 zero, & 859 AIMAG(Sr%sigcmee0 (iband,ik_ibz,isp))*Ha_eV, & 860 AIMAG(Sr%ze0 (iband,ik_ibz,isp)), & 861 AIMAG(Sr%dsigmee0 (iband,ik_ibz,isp)), & 862 AIMAG(Sr%sigmee (iband,ik_ibz,isp))*Ha_eV, & 863 AIMAG(Sr%degw (iband,ik_ibz,isp))*Ha_eV, & 864 AIMAG(Sr%egw (iband,ik_ibz,isp))*Ha_eV, & 865 zero 866 if (verbose/=0) then 867 call wrtout(my_unt,msg,my_mode) 868 if (present(ydoc)) call ydoc%add_tabular_line(msg) 869 end if 870 end if 871 872 else 873 ! PAW+U+GW calculation. 874 ABI_ERROR("PAW+U+GW not yet implemented") 875 end if 876 877 end subroutine print_Sigma_QPSC
m_sigma/sigma_distribute_bks [ Functions ]
[ Top ] [ m_sigma ] [ Functions ]
NAME
sigma_distribute_bks
FUNCTION
Distribute the loop over (b,k,s) used to calculate the self-energy matrix elements taking into account the MPI distribution of the wavefunctions and the use of symmetries to reduce the BZ sum to an appropriate irreducible wedge.
INPUTS
nsppol can_symmetrize(nsppol)=.TRUE if symmetries can be used to reduce the number of k-points to be summed. Kmesh<kmesh_t> Qmesh<kmesh_t> Ltg_kgw<littlegroup_t> Wfd(wfdgw_t) mg0(3) kptgw(3) [bks_mask(Wfd%mband,Kmesh%nbz,nsppol)] [got(Wfd%nproc)]=The number of tasks already assigned to the nodes. [global]=If true, an MPI global communication is performed such that each node will have the same table. Useful if for implementing algorithms in which each node needs to know the global distribution of the tasks, not only the task it has to complete. Defaults to .FALSE.
OUTPUT
my_nbks proc_distrb(Wfd%mband,Kmesh%nbz,nsppol)
SIDE EFFECTS
Wfd%bks_tab
SOURCE
1693 subroutine sigma_distribute_bks(Wfd,Kmesh,Ltg_kgw,Qmesh,nsppol,can_symmetrize,kptgw,mg0,my_nbks,proc_distrb,got,bks_mask,global) 1694 1695 !Arguments ------------------------------------ 1696 !scalars 1697 integer,intent(in) :: nsppol 1698 integer,intent(out) :: my_nbks 1699 logical,optional,intent(in) :: global 1700 type(kmesh_t),intent(in) :: Kmesh,Qmesh 1701 type(littlegroup_t),intent(in) :: Ltg_kgw 1702 type(wfdgw_t),intent(inout) :: Wfd 1703 !arrays 1704 integer,intent(in) :: mg0(3) 1705 integer,optional,intent(inout) :: got(Wfd%nproc) 1706 integer,intent(out) :: proc_distrb(Wfd%mband,Kmesh%nbz,nsppol) 1707 real(dp),intent(in) :: kptgw(3) 1708 logical,intent(in) :: can_symmetrize(Wfd%nsppol) 1709 logical,optional,intent(in) :: bks_mask(Wfd%mband,Kmesh%nbz,nsppol) 1710 1711 !Local variables------------------------------- 1712 !scalars 1713 integer :: ierr,ik_bz,ik_ibz,spin,iq_bz,my_nband 1714 !character(len=500) :: msg 1715 !arrays 1716 integer :: g0(3) 1717 real(dp) :: kgwmk(3) 1718 integer :: get_more(Wfd%nproc),my_band_list(Wfd%mband) 1719 !integer :: test(Wfd%mband,Kmesh%nbz,nsppol) 1720 logical :: bmask(Wfd%mband) 1721 1722 !************************************************************************ 1723 1724 call wfd%update_bkstab() 1725 1726 get_more=0; if (PRESENT(got)) get_more=got 1727 1728 ! Different distribution of tasks depending whether symmetries can be used or not. 1729 proc_distrb= xmpi_undefined_rank 1730 1731 do spin=1,Wfd%nsppol 1732 1733 if (can_symmetrize(spin)) then 1734 do ik_bz=1,Kmesh%nbz 1735 ik_ibz = Kmesh%tab(ik_bz) 1736 kgwmk= kptgw-Kmesh%bz(:,ik_bz) ! kptgw must be inside the BZ 1737 call findqg0(iq_bz,g0,kgwmk,Qmesh%nbz,Qmesh%bz,mG0) ! <- (mg0=mG0) Identify q_bz and G0 where q_bz+G0=k_gw-k_bz 1738 if (Ltg_kgw%ibzq(iq_bz)==1) then 1739 bmask=.FALSE.; bmask(1:Wfd%nband(ik_ibz,spin))=.TRUE. 1740 if (PRESENT(bks_mask)) bmask = bks_mask(:,ik_bz,spin) 1741 call wfd%distribute_bands(ik_ibz,spin,my_nband,my_band_list,got=get_more,bmask=bmask) 1742 if (my_nband>0) proc_distrb(my_band_list(1:my_nband),ik_bz,spin)=Wfd%my_rank 1743 end if 1744 end do 1745 1746 else 1747 ! No symmetries for this spin. Divide the full BZ among procs. 1748 do ik_bz=1,Kmesh%nbz 1749 ik_ibz = Kmesh%tab(ik_bz) 1750 bmask=.FALSE.; bmask(1:Wfd%nband(ik_ibz,spin))=.TRUE. 1751 if (PRESENT(bks_mask)) bmask = bks_mask(:,ik_bz,spin) 1752 call wfd%distribute_bands(ik_ibz,spin,my_nband,my_band_list,got=get_more,bmask=bmask) 1753 if (my_nband>0) proc_distrb(my_band_list(1:my_nband),ik_bz,spin)=Wfd%my_rank 1754 end do 1755 end if 1756 end do ! spin 1757 1758 if (PRESENT(global)) then 1759 if (global) then ! Each node will have the same table so that it will know how the tasks are distributed. 1760 proc_distrb = proc_distrb + 1 1761 where (proc_distrb == xmpi_undefined_rank+1) 1762 proc_distrb = 0 1763 end where 1764 call xmpi_sum(proc_distrb,Wfd%comm,ierr) 1765 where (proc_distrb == 0) 1766 proc_distrb = xmpi_undefined_rank 1767 elsewhere 1768 proc_distrb = proc_distrb -1 1769 end where 1770 !where (proc_distrb /= xmpi_undefined_rank) 1771 ! ltest = (ANY(proc_distrb == (/(ii,ii=0,Wfd%nproc-1)/))) 1772 !end where 1773 !if (.not.ltest) then 1774 ! write(std_out,*)proc_distrb 1775 ! ABI_BUG("Bug in the generation of proc_distrb table") 1776 !end if 1777 end if 1778 end if 1779 1780 my_nbks = COUNT(proc_distrb==Wfd%my_rank) 1781 if (PRESENT(got)) got=get_more 1782 1783 end subroutine sigma_distribute_bks
m_sigma/sigma_free [ Functions ]
[ Top ] [ m_sigma ] [ Functions ]
NAME
sigma_free
FUNCTION
Deallocate all associated pointers defined in the sigma_t data type.
INPUTS
OUTPUT
SOURCE
1022 subroutine sigma_free(Sr) 1023 1024 !Arguments ------------------------------------ 1025 class(sigma_t),intent(inout) :: Sr 1026 1027 ! ************************************************************************* 1028 1029 !@sigma_t 1030 !integer 1031 ABI_SFREE(Sr%maxbnd) 1032 ABI_SFREE(Sr%minbnd) 1033 1034 !real 1035 ABI_SFREE(Sr%degwgap) 1036 ABI_SFREE(Sr%egwgap) 1037 ABI_SFREE(Sr%en_qp_diago) 1038 ABI_SFREE(Sr%e0) 1039 ABI_SFREE(Sr%e0gap) 1040 ABI_SFREE(Sr%omega_r) 1041 ABI_SFREE(Sr%kptgw) 1042 ABI_SFREE(Sr%sigxme) 1043 ABI_SFREE(Sr%sigxcnofme) 1044 ABI_SFREE(Sr%x_mat) 1045 ABI_SFREE(Sr%vxcme) 1046 ABI_SFREE(Sr%vUme) 1047 1048 !complex 1049 ABI_SFREE(Sr%degw) 1050 ABI_SFREE(Sr%dsigmee0) 1051 ABI_SFREE(Sr%egw) 1052 ABI_SFREE(Sr%eigvec_qp) 1053 ABI_SFREE(Sr%m_ks_to_qp) 1054 ABI_SFREE(Sr%hhartree) 1055 ABI_SFREE(Sr%sigcme) 1056 ABI_SFREE(Sr%sigmee) 1057 ABI_SFREE(Sr%sigcmee0) 1058 ABI_SFREE(Sr%sigcmesi) 1059 ABI_SFREE(Sr%sigcme4sd) 1060 ABI_SFREE(Sr%sigxcme) 1061 ABI_SFREE(Sr%sigxcmesi) 1062 ABI_SFREE(Sr%sigxcme4sd) 1063 ABI_SFREE(Sr%ze0) 1064 ABI_SFREE(Sr%omega_i) 1065 ABI_SFREE(Sr%omega4sd) 1066 1067 end subroutine sigma_free
m_sigma/sigma_get_excene [ Functions ]
[ Top ] [ m_sigma ] [ Functions ]
NAME
sigma_get_excene
FUNCTION
Compute exchange correlation energy using MBB (nat. orb. functional approx.).
INPUTS
sigma<sigma_t>=Sigma results kmesh<kmesh_t>=BZ sampling. bands<band_t>=Bands with occupation factors
SOURCE
1135 real(dp) pure function sigma_get_excene(sigma, kmesh, bands) result(exc_energy) 1136 1137 !Arguments ------------------------------------ 1138 class(sigma_t),intent(in) :: sigma 1139 type(kmesh_t),intent(in) :: kmesh 1140 type(ebands_t),intent(in) :: bands 1141 1142 !Local variables------------------------------- 1143 !scalars 1144 integer :: ik,ib,spin 1145 real(dp) :: wtk,occ_bks 1146 1147 ! ************************************************************************* 1148 1149 exc_energy = zero 1150 1151 do spin=1,sigma%nsppol 1152 do ik=1,sigma%nkibz 1153 wtk = kmesh%wt(ik) 1154 do ib=sigma%b1gw,sigma%b2gw 1155 occ_bks = bands%occ(ib,ik,spin) 1156 if (sigma%nsig_ab==1) then 1157 if (sigma%nsppol==1) then 1158 exc_energy = exc_energy + sqrt( abs( half * occ_bks ) ) * wtk * sigma%sigxcnofme(ib,ik,spin) ! 2*sqrt(occ_i), occ in [0,2] -> [0,1]. 1159 else 1160 exc_energy = exc_energy + half * sqrt( abs( occ_bks ) ) * wtk * sigma%sigxcnofme(ib,ik,spin) ! 2*sqrt(occ_i), occ in [0,1] -> [0,1]. 1161 end if 1162 else 1163 exc_energy = exc_energy + half * sqrt( abs( occ_bks ) ) * wtk * SUM(sigma%sigxcnofme(ib,ik,:)) ! 2*sqrt(occ_i), occ in [0,1]. 1164 end if 1165 end do 1166 end do 1167 end do 1168 1169 end function sigma_get_excene
m_sigma/sigma_get_exene [ Functions ]
[ Top ] [ m_sigma ] [ Functions ]
NAME
sigma_get_exene
FUNCTION
Compute exchange energy.
INPUTS
sigma<sigma_t>=Sigma results kmesh<kmesh_t>=BZ sampling. bands<band_t>=Bands with occupation factors
SOURCE
1086 real(dp) pure function sigma_get_exene(sigma, kmesh, bands) result(ex_energy) 1087 1088 !Arguments ------------------------------------ 1089 class(sigma_t),intent(in) :: sigma 1090 type(kmesh_t),intent(in) :: kmesh 1091 type(ebands_t),intent(in) :: bands 1092 1093 !Local variables------------------------------- 1094 !scalars 1095 integer :: ik,ib,spin 1096 real(dp) :: wtk,occ_bks 1097 1098 ! ************************************************************************* 1099 1100 ex_energy = zero 1101 1102 do spin=1,sigma%nsppol 1103 do ik=1,sigma%nkibz 1104 wtk = kmesh%wt(ik) 1105 do ib=sigma%b1gw,sigma%b2gw 1106 occ_bks = bands%occ(ib,ik,spin) 1107 if (sigma%nsig_ab == 1) then 1108 ex_energy = ex_energy + half * occ_bks * wtk * sigma%sigxme(ib,ik,spin) 1109 else 1110 ex_energy = ex_energy + half * occ_bks * wtk * SUM(sigma%sigxme(ib,ik,:)) 1111 end if 1112 end do 1113 end do 1114 end do 1115 1116 end function sigma_get_exene
m_sigma/sigma_init [ Functions ]
[ Top ] [ m_sigma ] [ Functions ]
NAME
sigma_init
FUNCTION
Main creation method for the sigma_t data type.
INPUTS
usepawu= /=0 if we used DFT+U as starting point (only for PAW)
OUTPUT
TODO
Write documentation.
SOURCE
899 subroutine sigma_init(Sigp,nkibz,usepawu,Sr) 900 901 !Arguments ------------------------------------ 902 integer,intent(in) :: nkibz,usepawu 903 !scalars 904 type(sigparams_t),intent(in) :: Sigp 905 type(sigma_t),intent(inout) :: Sr 906 907 !Local variables------------------------------- 908 !scalars 909 integer :: b1gw,b2gw,mod10 910 911 ! ************************************************************************* 912 913 !@sigma_t 914 mod10=MOD(Sigp%gwcalctyp,10) 915 916 ! Copy important dimensions 917 Sr%nkptgw =Sigp%nkptgw 918 Sr%gwcalctyp =Sigp%gwcalctyp 919 Sr%deltae =Sigp%deltae 920 Sr%maxomega4sd=Sigp%maxomega4sd 921 Sr%maxomega_r =Sigp%maxomega_r 922 Sr%scissor_ene=Sigp%mbpt_sciss 923 924 !FIXME this should be done in sigma_allocate 925 ABI_MALLOC(Sr%minbnd,(Sr%nkptgw,Sigp%nsppol)) 926 ABI_MALLOC(Sr%maxbnd,(Sr%nkptgw,Sigp%nsppol)) 927 Sr%minbnd=Sigp%minbnd; Sr%maxbnd=Sigp%maxbnd 928 ABI_MALLOC(Sr%kptgw,(3,Sr%nkptgw)) 929 Sr%kptgw=Sigp%kptgw 930 931 Sr%b1gw =Sigp%minbdgw ! * min and Max GW band index over k and spin. 932 Sr%b2gw =Sigp%maxbdgw ! Used to dimension arrays. 933 Sr%nbnds =Sigp%nbnds 934 Sr%nkibz =nkibz 935 Sr%nsppol =Sigp%nsppol 936 Sr%nsig_ab =Sigp%nsig_ab 937 Sr%nomega_r =Sigp%nomegasr !FIXME change name 938 Sr%nomega_i =Sigp%nomegasi 939 Sr%nomega4sd=Sigp%nomegasrd 940 Sr%usepawu =usepawu 941 942 !====================================================== 943 ! === Allocate arrays in the sigma_t datatype === 944 !====================================================== 945 b1gw=Sr%b1gw 946 b2gw=Sr%b2gw 947 948 !TODO write routine to allocate all this stuff 949 950 ! hhartree(b1,b2,k,s)= <b1,k,s|T+v_{loc}+v_{nl}+v_{H}|b2,k,s> 951 ABI_CALLOC(Sr%hhartree,(b1gw:b2gw,b1gw:b2gw,Sr%nkibz,Sr%nsppol*Sr%nsig_ab)) 952 953 ! QP amplitudes and energies 954 ABI_CALLOC(Sr%en_qp_diago,(Sr%nbnds,Sr%nkibz,Sr%nsppol)) 955 ABI_CALLOC(Sr%eigvec_qp,(Sr%nbnds,Sr%nbnds,Sr%nkibz,Sr%nsppol)) 956 957 ! Dont know if it is better to do this here or in the sigma 958 ! * Initialize with KS wavefunctions and energies 959 !do ib=1,Sr%nbnds 960 ! Sr%en_qp_diago(ib,:,:)=en(:,ib,:) 961 ! Sr%eigvec_qp(ib,ib,:,:)=cone 962 !end do 963 964 ABI_CALLOC(Sr%vxcme, (b1gw:b2gw,Sr%nkibz,Sr%nsppol*Sr%nsig_ab)) 965 ABI_CALLOC(Sr%vUme, (b1gw:b2gw,Sr%nkibz,Sr%nsppol*Sr%nsig_ab)) 966 ABI_CALLOC(Sr%sigxme, (b1gw:b2gw,Sr%nkibz,Sr%nsppol*Sr%nsig_ab)) 967 ABI_CALLOC(Sr%sigxcnofme, (b1gw:b2gw,Sr%nkibz,Sr%nsppol*Sr%nsig_ab)) 968 ABI_CALLOC(Sr%x_mat, (b1gw:b2gw,b1gw:b2gw,Sr%nkibz,Sr%nsppol*Sr%nsig_ab)) 969 ABI_CALLOC(Sr%sigcme, (b1gw:b2gw,Sr%nkibz,Sr%nomega_r,Sr%nsppol*Sr%nsig_ab)) 970 ABI_CALLOC(Sr%sigxcme, (b1gw:b2gw,Sr%nkibz,Sr%nomega_r,Sr%nsppol*Sr%nsig_ab)) 971 ABI_CALLOC(Sr%sigcmee0, (b1gw:b2gw,Sr%nkibz,Sr%nsppol*Sr%nsig_ab)) 972 ABI_CALLOC(Sr%ze0, (b1gw:b2gw,Sr%nkibz,Sr%nsppol)) 973 ABI_CALLOC(Sr%dsigmee0, (b1gw:b2gw,Sr%nkibz,Sr%nsppol*Sr%nsig_ab)) 974 ABI_CALLOC(Sr%sigmee, (b1gw:b2gw,Sr%nkibz,Sr%nsppol*Sr%nsig_ab)) 975 ABI_CALLOC(Sr%degw, (b1gw:b2gw,Sr%nkibz,Sr%nsppol)) 976 ABI_CALLOC(Sr%e0, (Sr%nbnds,Sr%nkibz,Sr%nsppol)) 977 ABI_CALLOC(Sr%egw, (Sr%nbnds,Sr%nkibz,Sr%nsppol)) 978 ABI_CALLOC(Sr%e0gap, (Sr%nkibz,Sr%nsppol)) 979 ABI_CALLOC(Sr%degwgap, (Sr%nkibz,Sr%nsppol)) 980 ABI_CALLOC(Sr%egwgap, (Sr%nkibz,Sr%nsppol)) 981 982 ! These quantities are used to evaluate $\Sigma(E)$ around the KS\QP eigenvalue 983 ABI_CALLOC(Sr%omega4sd,(b1gw:b2gw,Sr%nkibz,Sr%nomega4sd,Sr%nsppol)) 984 ABI_CALLOC(Sr%sigcme4sd,(b1gw:b2gw,Sr%nkibz,Sr%nomega4sd,Sr%nsppol*Sr%nsig_ab)) 985 ABI_CALLOC(Sr%sigxcme4sd,(b1gw:b2gw,Sr%nkibz,Sr%nomega4sd,Sr%nsppol*Sr%nsig_ab)) 986 987 !TODO Find better treatment 988 ! Mesh along the real axis. 989 if (Sr%nomega_r>0) then 990 ABI_MALLOC(Sr%omega_r,(Sr%nomega_r)) 991 Sr%omega_r(:)=Sigp%omega_r(:) 992 end if 993 994 ! Analytical Continuation 995 ! FIXME omegasi should not be in Sigp% here we should construct the mesh 996 if (mod10==1) then 997 ABI_MALLOC(Sr%omega_i,(Sr%nomega_i)) 998 Sr%omega_i=Sigp%omegasi 999 ABI_MALLOC(Sr%sigcmesi ,(b1gw:b2gw,Sr%nkibz,Sr%nomega_i,Sr%nsppol*Sr%nsig_ab)) 1000 ABI_MALLOC(Sr%sigxcmesi,(b1gw:b2gw,Sr%nkibz,Sr%nomega_i,Sr%nsppol*Sr%nsig_ab)) 1001 Sr%sigcmesi=czero; Sr%sigxcmesi=czero 1002 end if 1003 1004 end subroutine sigma_init
m_sigma/sigma_ncwrite [ Functions ]
[ Top ] [ m_sigma ] [ Functions ]
NAME
sigma_ncwrite
FUNCTION
Save the data stored in the sigma_t data type on a NETCDF file.
INPUTS
filename
OUTPUT
SOURCE
1411 integer function sigma_ncwrite(Sigp,Er,Sr,ncid) result (ncerr) 1412 1413 !Arguments ------------------------------------ 1414 !scalars 1415 integer,intent(in) :: ncid 1416 type(sigparams_t),target,intent(in) :: Sigp 1417 type(Epsilonm1_results),target,intent(in) :: Er 1418 type(sigma_t),target,intent(in) :: Sr 1419 1420 !Local variables --------------------------------------- 1421 #ifdef HAVE_NETCDF 1422 !scalars 1423 integer :: nbgw,ndim_sig,b1gw,b2gw,cplex 1424 !character(len=500) :: msg 1425 !arrays 1426 real(dp),allocatable :: rdata2(:,:),rdata4(:,:,:,:),rdata5(:,:,:,:,:) 1427 1428 ! ************************************************************************* 1429 1430 !@sigma_t 1431 cplex=2; b1gw=Sr%b1gw; b2gw=Sr%b2gw; nbgw=b2gw-b1gw+1 1432 ndim_sig=Sr%nsppol*Sr%nsig_ab 1433 1434 ncerr = nctk_def_dims(ncid, [& 1435 nctkdim_t("cplex", cplex), nctkdim_t("b1gw", sr%b1gw), nctkdim_t("b2gw", sr%b2gw),& 1436 nctkdim_t("nbgw", nbgw), nctkdim_t("nkptgw", sr%nkptgw), nctkdim_t("ndim_sig", ndim_sig), & 1437 nctkdim_t("nomega4sd", sr%nomega4sd), nctkdim_t("nsig_ab", sr%nsig_ab)], defmode=.True.) 1438 NCF_CHECK(ncerr) 1439 1440 ! No. of real frequencies, might be zero. 1441 if (Sr%nomega_r > 0) then 1442 NCF_CHECK(nctk_def_dims(ncid, nctkdim_t("nomega_r", Sr%nomega_r))) 1443 end if 1444 1445 ! No. of imaginary frequencies, might be zero. 1446 if (Sr%nomega_i > 0) then 1447 NCF_CHECK(nctk_def_dims(ncid, nctkdim_t("nomega_i", Sr%nomega_i))) 1448 end if 1449 1450 ! ======================= 1451 ! == Define variables === 1452 ! ======================= 1453 ! parameters of the calculation. 1454 1455 ncerr = nctk_def_iscalars(ncid, [character(len=nctk_slen) :: 'sigma_nband', 'scr_nband', 'gwcalctyp', 'usepawu']) 1456 NCF_CHECK(ncerr) 1457 1458 ncerr = nctk_def_dpscalars(ncid, [character(len=nctk_slen) :: & 1459 & 'ecutwfn', 'ecuteps', 'ecutsigx', 'omegasrdmax', 'deltae', 'omegasrmax', 'scissor_ene']) 1460 NCF_CHECK(ncerr) 1461 1462 ! TODO: Decrease size of file: Remove arrays whose size scale as mband ** 2 1463 ! especially those that are not commonly used e.g. hhartree. 1464 ncerr = nctk_def_arrays(ncid, [& 1465 nctkarr_t("kptgw", "dp", "number_of_reduced_dimensions, nkptgw"),& 1466 nctkarr_t("minbnd", "i", "nkptgw, number_of_spins"),& 1467 nctkarr_t("maxbnd", "i", "nkptgw, number_of_spins"), & 1468 nctkarr_t('degwgap', "dp", 'number_of_kpoints, number_of_spins'),& 1469 nctkarr_t('egwgap', "dp", 'number_of_kpoints, number_of_spins'),& 1470 nctkarr_t('en_qp_diago', "dp",'max_number_of_states, number_of_kpoints, number_of_spins'),& 1471 nctkarr_t('e0', "dp", 'max_number_of_states, number_of_kpoints, number_of_spins'),& 1472 nctkarr_t('e0gap', "dp", 'number_of_kpoints, number_of_spins'),& 1473 nctkarr_t('sigxme', "dp", 'nbgw, number_of_kpoints, ndim_sig'),& 1474 nctkarr_t('vxcme', "dp", 'nbgw, number_of_kpoints, ndim_sig'),& 1475 nctkarr_t('degw', "dp", 'cplex, nbgw, number_of_kpoints, number_of_spins'),& 1476 nctkarr_t('dsigmee0', "dp", 'cplex, nbgw, number_of_kpoints, ndim_sig'),& 1477 nctkarr_t('egw', "dp",'cplex, max_number_of_states, number_of_kpoints, number_of_spins'),& 1478 nctkarr_t('eigvec_qp', "dp",'cplex, max_number_of_states, max_number_of_states, number_of_kpoints, number_of_spins'),& 1479 nctkarr_t('hhartree', "dp",'cplex, nbgw, nbgw, number_of_kpoints, ndim_sig'),& 1480 nctkarr_t('sigmee', "dp", 'cplex, nbgw, number_of_kpoints, ndim_sig'),& 1481 nctkarr_t('sigcmee0', "dp",'cplex, nbgw, number_of_kpoints, ndim_sig'),& 1482 nctkarr_t('sigcme4sd', "dp",'cplex, nbgw, number_of_kpoints, nomega4sd, ndim_sig'),& 1483 nctkarr_t('sigxcme4sd', "dp", 'cplex, nbgw, number_of_kpoints, nomega4sd, ndim_sig'),& 1484 nctkarr_t('ze0',"dp", 'cplex, nbgw, number_of_kpoints, number_of_spins'),& 1485 nctkarr_t('omega4sd', "dp", 'cplex, nbgw, number_of_kpoints, nomega4sd, number_of_spins')]) 1486 NCF_CHECK(ncerr) 1487 1488 if (Sr%usepawu==0) then 1489 ncerr = nctk_def_arrays(ncid, nctkarr_t("vUme", "dp", 'nbgw, number_of_kpoints, ndim_sig')) 1490 NCF_CHECK(ncerr) 1491 end if 1492 1493 if (Sr%nomega_r > 0) then 1494 ncerr = nctk_def_arrays(ncid, [& 1495 nctkarr_t('omega_r', "dp", "nomega_r"),& 1496 nctkarr_t('sigcme', "dp", 'cplex, nbgw, number_of_kpoints, nomega_r, ndim_sig'),& 1497 nctkarr_t('sigxcme', "dp", 'cplex, nbgw, number_of_kpoints, nomega_r, ndim_sig')]) 1498 NCF_CHECK(ncerr) 1499 end if 1500 1501 if (Sr%nomega_i > 0) then 1502 ncerr = nctk_def_arrays(ncid, [& 1503 nctkarr_t('sigxcmesi', "dp", 'cplex, nbgw, number_of_kpoints, nomega_i, ndim_sig'),& 1504 nctkarr_t('sigcmesi', "dp",'cplex, nbgw, number_of_kpoints, nomega_i, ndim_sig'),& 1505 nctkarr_t('omega_i', "dp", 'cplex, nomega_i')]) 1506 NCF_CHECK(ncerr) 1507 end if 1508 1509 if (allocated(sr%m_ks_to_qp)) then 1510 ncerr = nctk_def_arrays(ncid, [nctkarr_t('m_ks_to_qp', "dp", & 1511 "cplex, max_number_of_states, max_number_of_states, number_of_kpoints, number_of_spins")]) 1512 NCF_CHECK(ncerr) 1513 end if 1514 1515 ! ===================== 1516 ! === Start writing === 1517 ! ===================== 1518 NCF_CHECK(nctk_set_datamode(ncid)) 1519 NCF_CHECK(nf90_put_var(ncid, vid('ecutwfn'), Sigp%ecutwfn)) 1520 NCF_CHECK(nf90_put_var(ncid, vid('ecuteps'), Sigp%ecuteps)) 1521 NCF_CHECK(nf90_put_var(ncid, vid('ecutsigx'), Sigp%ecutsigx)) 1522 NCF_CHECK(nf90_put_var(ncid, vid('sigma_nband'), Sigp%nbnds)) 1523 NCF_CHECK(nf90_put_var(ncid, vid('scr_nband'), Er%Hscr%nbnds_used)) 1524 NCF_CHECK(nf90_put_var(ncid, vid('gwcalctyp'), Sr%gwcalctyp)) 1525 NCF_CHECK(nf90_put_var(ncid, vid('usepawu'), Sr%usepawu)) 1526 NCF_CHECK(nf90_put_var(ncid, vid('kptgw'), Sr%kptgw)) 1527 NCF_CHECK(nf90_put_var(ncid, vid('minbnd'), Sr%minbnd)) 1528 NCF_CHECK(nf90_put_var(ncid, vid('maxbnd'),Sr%maxbnd)) 1529 NCF_CHECK(nf90_put_var(ncid, vid('omegasrdmax'), Sr%maxomega4sd*Ha_eV)) 1530 NCF_CHECK(nf90_put_var(ncid, vid('deltae'), Sr%deltae*Ha_eV)) 1531 NCF_CHECK(nf90_put_var(ncid, vid('omegasrmax'), Sr%maxomega_r*Ha_eV)) 1532 NCF_CHECK(nf90_put_var(ncid, vid('scissor_ene'), Sr%scissor_ene*Ha_eV)) 1533 NCF_CHECK(nf90_put_var(ncid, vid('degwgap'), Sr%degwgap*Ha_eV)) 1534 NCF_CHECK(nf90_put_var(ncid, vid('egwgap'), Sr%egwgap*Ha_eV)) 1535 NCF_CHECK(nf90_put_var(ncid, vid('en_qp_diago'), Sr%en_qp_diago*Ha_eV)) 1536 NCF_CHECK(nf90_put_var(ncid, vid('e0'), Sr%e0*Ha_eV)) 1537 NCF_CHECK(nf90_put_var(ncid, vid('e0gap'), Sr%e0gap*Ha_eV)) 1538 1539 if (Sr%nomega_r>0) then 1540 NCF_CHECK(nf90_put_var(ncid, vid('omega_r'), Sr%omega_r*Ha_eV)) 1541 end if 1542 1543 NCF_CHECK(nf90_put_var(ncid, vid('sigxme'), Sr%sigxme*Ha_eV)) 1544 NCF_CHECK(nf90_put_var(ncid, vid('vxcme'), Sr%vxcme*Ha_eV)) 1545 NCF_CHECK(nf90_put_var(ncid, vid('vUme'), Sr%vUme*Ha_eV)) 1546 1547 ! Have to transfer complex arrays 1548 ABI_MALLOC(rdata4,(cplex,b1gw:b2gw,Sr%nkibz,Sr%nsppol)) 1549 rdata4=c2r(Sr%degw) 1550 NCF_CHECK(nf90_put_var(ncid, vid('degw'), rdata4*Ha_eV)) 1551 ABI_FREE(rdata4) 1552 1553 ABI_MALLOC(rdata4,(cplex,b1gw:b2gw,Sr%nkibz,Sr%nsppol*Sr%nsig_ab)) 1554 rdata4=c2r(Sr%dsigmee0) 1555 NCF_CHECK(nf90_put_var(ncid, vid('dsigmee0'), rdata4)) 1556 ABI_FREE(rdata4) 1557 1558 ABI_MALLOC(rdata4,(cplex,Sr%nbnds,Sr%nkibz,Sr%nsppol)) 1559 rdata4=c2r(Sr%egw) 1560 NCF_CHECK(nf90_put_var(ncid, vid('egw'), rdata4*Ha_eV)) 1561 ABI_FREE(rdata4) 1562 1563 ABI_MALLOC(rdata5,(cplex,Sr%nbnds,Sr%nbnds,Sr%nkibz,Sr%nsppol)) 1564 rdata5=c2r(Sr%eigvec_qp) 1565 NCF_CHECK(nf90_put_var(ncid, vid('eigvec_qp'), rdata5)) 1566 ABI_FREE(rdata5) 1567 1568 ABI_MALLOC(rdata5,(cplex,nbgw,nbgw,Sr%nkibz,Sr%nsppol*Sr%nsig_ab)) 1569 rdata5=c2r(Sr%hhartree) 1570 NCF_CHECK(nf90_put_var(ncid, vid('hhartree'), rdata5*Ha_eV)) 1571 ABI_FREE(rdata5) 1572 1573 if (Sr%nomega_r>0) then 1574 ABI_MALLOC(rdata5,(cplex,nbgw,Sr%nkibz,Sr%nomega_r,Sr%nsppol*Sr%nsig_ab)) 1575 rdata5=c2r(Sr%sigcme) 1576 NCF_CHECK(nf90_put_var(ncid, vid('sigcme'), rdata5*Ha_eV)) 1577 ABI_FREE(rdata5) 1578 end if 1579 1580 ABI_MALLOC(rdata4,(cplex,nbgw,Sr%nkibz,Sr%nsppol*Sr%nsig_ab)) 1581 rdata4=c2r(Sr%sigmee) 1582 NCF_CHECK(nf90_put_var(ncid, vid('sigmee'), rdata4*Ha_eV)) 1583 ABI_FREE(rdata4) 1584 1585 ABI_MALLOC(rdata4,(cplex,nbgw,Sr%nkibz,Sr%nsppol*Sr%nsig_ab)) 1586 rdata4=c2r(Sr%sigcmee0) 1587 NCF_CHECK(nf90_put_var(ncid, vid('sigcmee0'), rdata4*Ha_eV)) 1588 ABI_FREE(rdata4) 1589 1590 if (Sr%nomega_i>0) then 1591 ABI_MALLOC(rdata5,(cplex,nbgw,Sr%nkibz,Sr%nomega_i,Sr%nsppol*Sr%nsig_ab)) 1592 rdata5=c2r(Sr%sigcmesi) 1593 NCF_CHECK(nf90_put_var(ncid, vid('sigcmesi'), rdata5*Ha_eV)) 1594 ABI_FREE(rdata5) 1595 end if 1596 1597 ABI_MALLOC(rdata5,(cplex,nbgw,Sr%nkibz,Sr%nomega4sd,Sr%nsppol*Sr%nsig_ab)) 1598 rdata5=c2r(Sr%sigcme4sd) 1599 NCF_CHECK(nf90_put_var(ncid, vid('sigcme4sd'), rdata5*Ha_eV)) 1600 ABI_FREE(rdata5) 1601 1602 if (Sr%nomega_r>0) then 1603 ABI_MALLOC(rdata5,(cplex,nbgw,Sr%nkibz,Sr%nomega_r,Sr%nsppol*Sr%nsig_ab)) 1604 rdata5=c2r(Sr%sigxcme) 1605 NCF_CHECK(nf90_put_var(ncid, vid('sigxcme'), rdata5*Ha_eV)) 1606 ABI_FREE(rdata5) 1607 end if 1608 1609 if (Sr%nomega_i>0) then 1610 ABI_MALLOC(rdata5,(cplex,nbgw,Sr%nkibz,Sr%nomega_i,Sr%nsppol*Sr%nsig_ab)) 1611 rdata5=c2r(Sr%sigxcmesi) 1612 NCF_CHECK(nf90_put_var(ncid, vid('sigxcmesi'), rdata5*Ha_eV)) 1613 ABI_FREE(rdata5) 1614 end if 1615 1616 if (allocated(sr%m_ks_to_qp)) then 1617 ABI_MALLOC(rdata5,(cplex,Sr%nbnds,Sr%nbnds,Sr%nkibz,Sr%nsppol)) 1618 rdata5=c2r(Sr%m_ks_to_qp) 1619 NCF_CHECK(nf90_put_var(ncid, vid('m_ks_to_qp'), rdata5)) 1620 ABI_FREE(rdata5) 1621 end if 1622 1623 ABI_MALLOC(rdata5,(cplex,nbgw,Sr%nkibz,Sr%nomega4sd,Sr%nsppol*Sr%nsig_ab)) 1624 rdata5=c2r(Sr%sigxcme4sd) 1625 NCF_CHECK(nf90_put_var(ncid, vid('sigxcme4sd'), rdata5*Ha_eV)) 1626 ABI_FREE(rdata5) 1627 1628 ABI_MALLOC(rdata4,(cplex,nbgw,Sr%nkibz,Sr%nsppol)) 1629 rdata4=c2r(Sr%ze0) 1630 NCF_CHECK(nf90_put_var(ncid, vid('ze0'), rdata4)) 1631 ABI_FREE(rdata4) 1632 1633 if (Sr%nomega_i > 0) then 1634 ABI_MALLOC(rdata2,(cplex,Sr%nomega_i)) 1635 rdata2=c2r(Sr%omega_i) 1636 NCF_CHECK(nf90_put_var(ncid, vid('omega_i'), rdata2*Ha_eV)) 1637 ABI_FREE(rdata2) 1638 end if 1639 1640 ABI_MALLOC(rdata5,(cplex,nbgw,Sr%nkibz,Sr%nomega4sd,Sr%nsppol)) 1641 rdata5=c2r(Sr%omega4sd) 1642 NCF_CHECK(nf90_put_var(ncid, vid('omega4sd'), rdata5*Ha_eV)) 1643 ABI_FREE(rdata5) 1644 1645 #else 1646 ABI_ERROR('netcdf support is not activated.') 1647 #endif 1648 1649 contains 1650 integer function vid(vname) 1651 character(len=*),intent(in) :: vname 1652 vid = nctk_idname(ncid, vname) 1653 end function vid 1654 1655 end function sigma_ncwrite
m_sigma/sigma_t [ Types ]
NAME
sigma_t
FUNCTION
For the GW part of ABINIT, the sigma_t structured datatype gather the results of a sigma calculation.
TODO
ragged arrays (nk,nsppol) --> values ?
SOURCE
68 type,public :: sigma_t 69 70 integer :: b1gw, b2gw ! min and Max gw band indices over spin and k-points (used to dimension arrays) 71 integer :: gwcalctyp ! Flag defining the calculation type. 72 integer :: nkptgw ! No. of points calculated 73 integer :: nkibz ! No. of irreducible k-points. 74 integer :: nbnds ! Total number of bands 75 integer :: nomega_r ! No. of real frequencies for the spectral function. 76 integer :: nomega_i ! No. of frequencies along the imaginary axis. 77 integer :: nomega4sd ! No. of real frequencies to evaluate the derivative of $\Sigma(E)$. 78 integer :: nsig_ab ! 1 if nspinor=1,4 for noncollinear case. 79 integer :: nsppol ! No. of spin polarizations. 80 integer :: usepawu ! 1 if we are using DFT+U as starting point (only for PAW) 81 82 real(dp) :: deltae ! Frequency step for the calculation of d\Sigma/dE 83 real(dp) :: maxomega4sd ! Max frequency around E_ks for d\Sigma/dE. 84 real(dp) :: maxomega_r ! Max frequency for spectral function. 85 real(dp) :: scissor_ene ! Scissor energy value. zero for None. 86 87 integer,allocatable :: maxbnd(:,:) 88 ! (nkptgw,nsppol) 89 ! Max band index considered in GW for this k-point. 90 91 integer,allocatable :: minbnd(:,:) 92 ! (nkptgw,nsppol) 93 ! Min band index considered in GW for this k-point. 94 95 !real(dp),allocatable :: ame(:,:,:) 96 ! (nbnds,nkibz,nomega)) 97 ! Diagonal matrix elements of the spectral function. 98 ! Commented out, it can be calculated from the other quantities 99 100 real(dp),allocatable :: degwgap(:,:) 101 ! (nkibz, nsppol) 102 ! Difference btw the QP and the KS direct gap. 103 104 real(dp),allocatable :: egwgap(:,:) 105 ! (nkibz,nsppol)) 106 ! QP direct gap at each k-point and spin. 107 108 real(dp),allocatable :: en_qp_diago(:,:,:) 109 ! (nbnds,nkibz, nsppol)) 110 ! QP energies obtained from the diagonalization of the Hermitian approximation to Sigma (QPSCGW) 111 112 real(dp),allocatable :: e0(:,:,:) 113 ! (nbnds,nkibz,nsppol) 114 ! KS eigenvalues for each band, k-point and spin. In case of self-consistent? 115 116 real(dp),allocatable :: e0gap(:,:) 117 ! (nkibz,nsppol), 118 ! KS gap at each k-point, for each spin. 119 120 real(dp),allocatable :: omega_r(:) 121 ! (nomega_r) 122 ! real frequencies used for the self energy. 123 124 real(dp),allocatable :: kptgw(:,:) 125 ! (3,nkptgw) 126 ! ! TODO there is a similar array in sigparams_t 127 ! List of calculated k-points. 128 129 real(dp),allocatable :: sigxme(:,:,:) 130 ! (b1gw:b2gw,nkibz,nsppol*nsig_ab)) 131 ! Diagonal matrix elements $\<nks|\Sigma_x|nks\>$ 132 133 real(dp),allocatable :: sigxcnofme(:,:,:) 134 ! (b1gw:b2gw,nkibz,nsppol*nsig_ab)) 135 ! Diagonal matrix elements $\<nks|\Sigma_xc|nks\>$ taking sqrt(occs) in \Sigma_x, occs in [0,1] 136 137 complex(dp),allocatable :: x_mat(:,:,:,:) 138 ! (b1gw:b2gw, b1gw:b2gw, nkibz, nsppol*nsig_ab) 139 ! Matrix elements of $\<nks|\Sigma_x|mks\>$ 140 141 real(dp),allocatable :: vxcme(:,:,:) 142 ! (b1gw:b2gw,nkibz,nsppol*nsig_ab)) 143 ! $\<nks|v_{xc}[n_val]|nks\>$ matrix elements of vxc 144 ! NB: valence-only contribution i.e. computed without model core charge 145 146 real(dp),allocatable :: vUme(:,:,:) 147 ! (b1gw:b2gw,nkibz,nsppol*nsig_ab)) 148 ! $\<nks|v_{U}|nks\>$ for DFT+U. 149 150 complex(dpc),allocatable :: degw(:,:,:) 151 ! (b1gw:b2gw,nkibz,nsppol)) 152 ! Difference between the QP and the KS energies. 153 154 complex(dpc),allocatable :: dsigmee0(:,:,:) 155 ! (b1gw:b2gw,nkibz,nsppol*nsig_ab)) 156 ! Derivative of $\Sigma_c(E)$ calculated at the KS eigenvalue. 157 158 complex(dpc),allocatable :: egw(:,:,:) 159 ! (nbnds,nkibz,nsppol)) 160 ! QP energies, $\epsilon_{nks}^{QP}$. 161 162 complex(dpc),allocatable :: eigvec_qp(:,:,:,:) 163 ! (nbnds,nbnds,nkibz,nsppol)) 164 ! Expansion of the QP amplitudes in the QP basis set of the previous iteration. 165 166 complex(dpc),allocatable :: m_ks_to_qp(:,:,:,:) 167 ! (nbnds,nbnds,nkibz,nsppol)) 168 ! m_ks_to_qp(ib,jb,k,s) := <\psi_{ib,k,s}^{KS}|\psi_{jb,k,s}^{QP}> 169 170 complex(dpc),allocatable :: hhartree(:,:,:,:) 171 ! (b1gw:b2gw,b1gw:b2gw,nkibz,nsppol*nsig_ab) 172 ! $\<nks|T+v_H+v_{loc}+v_{nl}|mks\>$ 173 ! Note that v_{loc} does not include the contribution to vxc(r) given by the model core charge. 174 175 complex(dpc),allocatable :: sigcme(:,:,:,:) 176 ! (b1gw:b2gw,nkibz,nomega_r,nsppol*nsig_ab)) 177 ! $\<nks|\Sigma_{c}(E)|nks\>$ at each nomega_r frequency 178 179 complex(dpc),allocatable :: sigmee(:,:,:) 180 ! (b1gw:b2gw,nkibz,nsppol*nsig_ab)) 181 ! $\Sigma_{xc}E_{KS} + (E_{QP}- E_{KS})*dSigma/dE_KS 182 183 complex(dpc),allocatable :: sigcmee0(:,:,:) 184 ! (b1gw:b2gw,nkibz,nsppol*nsig_ab)) 185 ! Diagonal matrix elements of $\Sigma_c(E)$ calculated at the KS energy $E_{KS}$ 186 187 complex(dpc),allocatable :: sigcmesi(:,:,:,:) 188 ! (b1gw:b2gw,nkibz,nomega_i,nsppol*nsig_ab)) 189 ! Matrix elements of $\Sigma_c$ along the imaginary axis. 190 ! Only used in case of analytical continuation. 191 192 complex(dpc),allocatable :: sigcme4sd(:,:,:,:) 193 ! (b1gw:b2gw,nkibz,nomega4sd,nsppol*nsig_ab)) 194 ! Diagonal matrix elements of \Sigma_c around the zeroth order eigenvalue (usually KS). 195 196 complex(dpc),allocatable :: sigxcme(:,:,:,:) 197 ! (b1gw:b2gw,nkibz,nomega_r,nsppol*nsig_ab)) 198 ! $\<nks|\Sigma_{xc}(E)|nks\>$ at each real frequency frequency. 199 200 complex(dpc),allocatable :: sigxcmesi(:,:,:,:) 201 ! (b1gw:b2gw,nkibz,nomega_i,nsppol*nsig_ab)) 202 ! Matrix elements of $\Sigma_{xc}$ along the imaginary axis. 203 ! Only used in case of analytical continuation. 204 205 complex(dpc),allocatable :: sigxcme4sd(:,:,:,:) 206 ! (b1gw:b2gw,nkibz,nomega4sd,nsppol*nsig_ab)) 207 ! Diagonal matrix elements of \Sigma_xc for frequencies around the zeroth order eigenvalues. 208 209 complex(dpc),allocatable :: ze0(:,:,:) 210 ! (b1gw:b2gw,nkibz,nsppol)) 211 ! renormalization factor. $(1-\dfrac{\partial\Sigma_c} {\partial E_{KS}})^{-1}$ 212 213 complex(dpc),allocatable :: omega_i(:) 214 ! (nomega_i) 215 ! Frequencies along the imaginary axis used for the analytical continuation. 216 217 complex(dpc),allocatable :: omega4sd(:,:,:,:) 218 ! (b1gw:b2gw,nkibz,nomega4sd,nsppol). 219 ! Frequencies used to evaluate the Derivative of Sigma. 220 221 end type sigma_t 222 223 public :: sigma_init ! Initialize the object 224 public :: sigma_free ! Deallocate memory 225 public :: sigma_get_exene ! Compute exchange energy. 226 public :: sigma_get_excene ! Compute exchange-correlation MBB (Nat. Orb. Funct. Approx.) energy. 227 public :: mels_get_haene ! Compute hartree energy. 228 public :: mels_get_kiene ! Compute kinetic energy. 229 public :: sigma_ncwrite ! Write data in netcdf format. 230 public :: write_sigma_header 231 public :: write_sigma_results 232 public :: print_Sigma_perturbative 233 public :: print_Sigma_QPSC 234 public :: sigma_distribute_bks
m_sigma/write_sigma_header [ Functions ]
[ Top ] [ m_sigma ] [ Functions ]
NAME
write_sigma_header
FUNCTION
Write basic info and dimensions used during the calculation of the QP correctoions (optdriver==4).
INPUTS
Sigp=sigparams_t Cryst<crystal_t>= Info on the Crystal structure Kmesh<kmesh_t>= Description of the BZ sampling.
OUTPUT
(for writing routines, no output) otherwise, should be described
NOTES
SOURCE
261 subroutine write_sigma_header(Sigp,Er,Cryst,Kmesh,Qmesh) 262 263 !Arguments ------------------------------------ 264 !scalars 265 type(kmesh_t),intent(in) :: Kmesh,Qmesh 266 type(crystal_t),intent(in) :: Cryst 267 type(Epsilonm1_results),intent(in) :: Er 268 type(sigparams_t),intent(in) :: Sigp 269 270 !Local variables------------------------------- 271 !scalars 272 integer :: gwcalctyp,mod10 273 character(len=500) :: msg 274 275 ! ************************************************************************* 276 277 call wrtout([std_out, ab_out], ' SIGMA fundamental parameters:') 278 279 gwcalctyp=Sigp%gwcalctyp 280 mod10=MOD(Sigp%gwcalctyp,10) 281 282 SELECT CASE (mod10) 283 CASE (0) 284 write(msg,'(a,i2)')' PLASMON POLE MODEL ',Sigp%ppmodel 285 CASE (1) 286 write(msg,'(a)')' ANALYTIC CONTINUATION' 287 CASE (2) 288 write(msg,'(a)')' CONTOUR DEFORMATION' 289 CASE (5) 290 write(msg,'(a)')' Hartree-Fock' 291 CASE (6) 292 write(msg,'(a)')' Screened Exchange' 293 CASE (7) 294 write(msg,'(a)')' COHSEX' 295 CASE (8) 296 write(msg,'(a,i2)')' MODEL GW with PLASMON POLE MODEL ',Sigp%ppmodel 297 CASE (9) 298 write(msg,'(a)')' MODEL GW without PLASMON POLE MODEL' 299 CASE DEFAULT 300 write(msg,'(a,i3)')' Wrong value for Sigp%gwcalctyp = ',Sigp%gwcalctyp 301 ABI_BUG(msg) 302 END SELECT 303 call wrtout([std_out, ab_out], msg) 304 305 write(msg,'(a,i12)')' number of plane-waves for SigmaX ',Sigp%npwx 306 call wrtout([std_out, ab_out], msg) 307 write(msg,'(a,i12)')' number of plane-waves for SigmaC and W ',Sigp%npwc 308 call wrtout([std_out, ab_out], msg) 309 write(msg,'(a,i12)')' number of plane-waves for wavefunctions ',Sigp%npwwfn 310 call wrtout([std_out, ab_out], msg) 311 write(msg,'(a,i12)')' number of bands ',Sigp%nbnds 312 call wrtout([std_out, ab_out], msg) 313 write(msg,'(a,i12)')' number of independent spin polarizations ',Sigp%nsppol 314 call wrtout([std_out, ab_out], msg) 315 write(msg,'(a,i12)')' number of spinorial components ',Sigp%nspinor 316 call wrtout([std_out, ab_out], msg) 317 write(msg,'(a,i12)')' number of k-points in IBZ ',Kmesh%nibz 318 call wrtout([std_out, ab_out], msg) 319 write(msg,'(a,i12)')' number of q-points in IBZ ',Qmesh%nibz 320 call wrtout([std_out, ab_out], msg) 321 write(msg,'(a,i12)')' number of symmetry operations ',Cryst%nsym 322 call wrtout([std_out, ab_out], msg) 323 write(msg,'(a,i12)')' number of k-points in BZ ',Kmesh%nbz 324 call wrtout([std_out, ab_out], msg) 325 write(msg,'(a,i12)')' number of q-points in BZ ',Qmesh%nbz 326 call wrtout([std_out, ab_out], msg) 327 write(msg,'(a,i12)')' number of frequencies for dSigma/dE ',Sigp%nomegasrd 328 call wrtout([std_out, ab_out], msg) 329 write(msg,'(a,f12.2)')' frequency step for dSigma/dE [eV] ',Sigp%deltae*Ha_eV 330 call wrtout([std_out, ab_out], msg) 331 write(msg,'(a,i12)')' number of omega for Sigma on real axis ',Sigp%nomegasr 332 call wrtout([std_out, ab_out], msg) 333 write(msg,'(a,f12.2)')' max omega for Sigma on real axis [eV] ',Sigp%maxomega_r*Ha_eV 334 call wrtout([std_out, ab_out], msg) 335 write(msg,'(a,f12.2)')' zcut for avoiding poles [eV] ',Sigp%zcut*Ha_eV 336 call wrtout([std_out, ab_out], msg) 337 338 if (Sigp%mbpt_sciss>0.1d-4) then 339 write(msg,'(a,f12.2)')' scissor energy [eV] ',Sigp%mbpt_sciss*Ha_eV 340 call wrtout([std_out, ab_out], msg) 341 end if 342 343 if (mod10==1) then 344 write(msg,'(a,i12)')' number of imaginary frequencies for Sigma',Sigp%nomegasi 345 call wrtout([std_out, ab_out], msg) 346 ! MRM not needed for GW 1RDM 347 if(gwcalctyp/=21) then 348 write(msg,'(a,f12.2)')' max omega for Sigma on imag axis [eV] ',Sigp%omegasimax*Ha_eV 349 call wrtout([std_out, ab_out], msg) 350 endif 351 end if 352 353 if (sigma_needs_w(Sigp)) then 354 write(msg,'(2a)')ch10,' EPSILON^-1 parameters (SCR file):' 355 call wrtout([std_out, ab_out], msg) 356 write(msg,'(a,i12)')' dimension of the eps^-1 matrix on file ',Er%Hscr%npwe 357 call wrtout([std_out, ab_out], msg) 358 write(msg,'(a,i12)')' dimension of the eps^-1 matrix used ',Er%npwe 359 call wrtout([std_out, ab_out], msg) 360 write(msg,'(a,i12)')' number of plane-waves for wavefunctions ',Er%Hscr%npwwfn_used 361 call wrtout([std_out, ab_out], msg) 362 write(msg,'(a,i12)')' number of bands ',Er%Hscr%nbnds_used 363 call wrtout([std_out, ab_out], msg) 364 write(msg,'(a,i12)')' number of q-points in IBZ ',Qmesh%nibz 365 call wrtout([std_out, ab_out], msg) 366 write(msg,'(a,i12)')' number of frequencies ',Er%nomega 367 call wrtout([std_out, ab_out], msg) 368 write(msg,'(a,i12)')' number of real frequencies ',Er%nomega_r 369 call wrtout([std_out, ab_out], msg) 370 write(msg,'(a,i12)')' number of imag frequencies ',Er%nomega_i 371 call wrtout([std_out, ab_out], msg) 372 end if 373 374 ! MRM not needed for GW 1RDM 375 if(gwcalctyp/=21) then 376 write(msg,'(3a)')ch10,' matrix elements of self-energy operator (all in [eV])',ch10 377 call wrtout([std_out, ab_out], msg) 378 endif 379 380 if (gwcalctyp<10) then 381 write(msg,'(a)')' Perturbative Calculation' 382 else if (gwcalctyp<20) then 383 write(msg,'(a)')' Self-Consistent on Energies only' 384 else 385 write(msg,'(a)')' Self-Consistent on Energies and Wavefunctions' 386 end if 387 call wrtout([std_out, ab_out], msg) 388 389 end subroutine write_sigma_header
m_sigma/write_sigma_results [ Functions ]
[ Top ] [ m_sigma ] [ Functions ]
NAME
write_sigma_results
FUNCTION
Write the final results of the GW calculation.
INPUTS
KS_BSt<ebands_t>=Info on the KS band structure energies. %eig(mband,nkibz,nsppol)= KS energies ikibz= index of the k-point in the array kibz, where GW corrections are calculated ikcalc= index of the k-point in the array Sigp%kptgw2bz Sigp=sigparams_t datatype sr=sigma results datatype
OUTPUT
(for writing routines, no output) otherwise, should be described
SOURCE
415 subroutine write_sigma_results(ikcalc,ikibz,Sigp,Sr,KS_BSt) 416 417 !Arguments ------------------------------------ 418 !scalars 419 integer,intent(in) :: ikcalc,ikibz 420 type(ebands_t),intent(in) :: KS_BSt 421 type(sigparams_t),intent(in) :: Sigp 422 type(sigma_t),intent(in) :: Sr 423 424 !Local variables------------------------------- 425 !scalars 426 integer :: ib,io,is 427 integer :: gwcalctyp,mod10 428 character(len=500) :: msg 429 type(yamldoc_t) :: ydoc 430 !arrays 431 character(len=12) :: tag_spin(2) 432 433 ! ************************************************************************* 434 435 gwcalctyp=Sigp%gwcalctyp 436 mod10=MOD(Sigp%gwcalctyp,10) 437 438 ! unt_gw: File with GW corrections. 439 ! unt_sig: Self-energy as a function of frequency. 440 ! unt_sgr: Derivative wrt omega of the Self-energy. 441 ! unt_sigc: Sigma_c(eik) MRM 442 ! unt_sgm: Sigma on the Matsubara axis (imag axis) 443 444 tag_spin=(/' ',' '/); if (Sr%nsppol==2) tag_spin=(/', SPIN UP ',', SPIN DOWN'/) 445 446 do is=1,Sr%nsppol 447 write(msg,'(2a,3f8.3,a)')ch10,' k = ',Sigp%kptgw(:,ikcalc),tag_spin(is) 448 call wrtout(std_out,msg) 449 !call wrtout(ab_out,msg) 450 451 msg = ' Band E0 <VxcDFT> SigX SigC(E0) Z dSigC/dE Sig(E) E-E0 E' 452 if (Sr%usepawu/=0) then 453 msg = ' Band E0 <VxcDFT> <H_U> SigX SigC(E0) Z dSigC/dE Sig(E) E-E0 E' 454 end if 455 456 if (gwcalctyp>=10) then 457 write(msg,'(2a)')& 458 ' Band E_DFT <VxcDFT> E(N-1) <Hhartree> SigX SigC[E(N-1)]',& 459 ' Z dSigC/dE Sig[E(N)] DeltaE E(N)_pert E(N)_diago' 460 end if 461 call wrtout(std_out,msg) 462 463 ydoc = yamldoc_open('SelfEnergy_ee', width=11, real_fmt='(3f8.3)') 464 call ydoc%add_real1d('kpoint', Sigp%kptgw(:,ikcalc)) 465 call ydoc%add_int('spin', is, int_fmt="(i1)") 466 call ydoc%add_real('KS_gap', Sr%e0gap(ikibz,is)*Ha_eV) 467 call ydoc%add_real('QP_gap', Sr%egwgap(ikibz,is)*Ha_eV) 468 call ydoc%add_real('Delta_QP_KS', Sr%degwgap(ikibz,is)*Ha_eV) 469 call ydoc%open_tabular('data', tag='SigmaeeData') 470 call ydoc%add_tabular_line(msg) 471 472 write(unt_gw,'(3f10.6)')Sigp%kptgw(:,ikcalc) 473 write(unt_gw,'(i4)')Sigp%maxbnd(ikcalc,is)-Sigp%minbnd(ikcalc,is)+1 474 475 write(unt_gwdiag,'(3f10.6)')Sigp%kptgw(:,ikcalc) 476 write(unt_gwdiag,'(i4)')Sigp%maxbnd(ikcalc,is)-Sigp%minbnd(ikcalc,is)+1 477 478 write(unt_sig,'("# k = ",3f10.6)')Sigp%kptgw(:,ikcalc) 479 write(unt_sig,'("# b = ",2i10)')Sigp%minbnd(ikcalc,is),Sigp%maxbnd(ikcalc,is) 480 481 write(unt_sgr,'("# k = ",3f10.6)')Sigp%kptgw(:,ikcalc) 482 write(unt_sgr,'("# b = ",2i10)')Sigp%minbnd(ikcalc,is),Sigp%maxbnd(ikcalc,is) 483 484 write(unt_sigc,'("# k = ",3f10.6)')Sigp%kptgw(:,ikcalc) 485 write(unt_sigc,'("# b = ",2i10)')Sigp%minbnd(ikcalc,is),Sigp%maxbnd(ikcalc,is) 486 487 do ib=Sigp%minbnd(ikcalc,is),Sigp%maxbnd(ikcalc,is) 488 if (gwcalctyp >= 10) then 489 call print_Sigma_QPSC(Sr,ikibz,ib,is,KS_BSt,unit=dev_null, ydoc=ydoc) 490 call print_Sigma_QPSC(Sr,ikibz,ib,is,KS_BSt,unit=std_out,prtvol=1) 491 492 write(unt_gwdiag,'(i6,3f9.4)') & 493 ib, & 494 Sr%en_qp_diago(ib,ikibz,is)*Ha_eV, & 495 (Sr%en_qp_diago(ib,ikibz,is) - KS_BSt%eig(ib,ikibz,is))*Ha_eV,& 496 zero 497 498 else 499 ! If not ppmodel, write out also the imaginary part in ab_out 500 SELECT CASE(mod10) 501 CASE(1,2) 502 call print_Sigma_perturbative(Sr,ikibz,ib,is,unit=dev_null,ydoc=ydoc,prtvol=1) 503 CASE DEFAULT 504 call print_Sigma_perturbative(Sr,ikibz,ib,is,unit=dev_null,ydoc=ydoc) 505 END SELECT 506 call print_Sigma_perturbative(Sr,ikibz,ib,is,unit=std_out,prtvol=1) 507 end if 508 509 write(unt_gw,'(i6,3f9.4)') & 510 ib, & 511 REAL (Sr%egw (ib,ikibz,is))*Ha_eV,& 512 REAL (Sr%degw(ib,ikibz,is))*Ha_eV,& 513 AIMAG(Sr%egw (ib,ikibz,is))*Ha_eV 514 end do !ib 515 516 if (Sr%e0gap(ikibz,is)**2+Sr%egwgap(ikibz,is)**2+Sr%degwgap(ikibz,is)**2 > tol10) then 517 ! Output the direct gap for each spin 518 ! If all the gaps are zero, this means that it could not be computed in the calling routine 519 write(msg,'(2a,f8.3)')ch10,' E^0_gap ',Sr%e0gap(ikibz,is)*Ha_eV 520 call wrtout(std_out,msg) 521 write(msg,'(a,f8.3)') ' E^GW_gap ',Sr%egwgap(ikibz,is)*Ha_eV 522 call wrtout(std_out,msg) 523 write(msg,'(a,f8.3,a)') ' DeltaE^GW_gap ',Sr%degwgap(ikibz,is)*Ha_eV,ch10 524 call wrtout(std_out,msg) 525 end if 526 527 call ydoc%write_and_free(ab_out) 528 529 ! Output of the spectral function 530 do io=1,Sr%nomega_r 531 write(unt_sig,'(100(e12.5,2x))')& 532 REAL(Sr%omega_r(io))*Ha_eV,& 533 (REAL(Sr%sigxcme(ib,ikibz,io,is))*Ha_eV,& 534 AIMAG(Sr%sigxcme(ib,ikibz,io,is))*Ha_eV,& 535 gw_spectral_function(Sr,io,ib,ikibz,is),& 536 ib=Sigp%minbnd(ikcalc,is),Sigp%maxbnd(ikcalc,is)) 537 end do 538 ! 539 do ib=Sigp%minbnd(ikcalc,is),Sigp%maxbnd(ikcalc,is) 540 write(unt_sgr,'("# ik, ib",2i5)')ikibz,ib 541 do io=1,Sr%nomega4sd 542 write(unt_sgr,'(100(e12.5,2x))') & 543 REAL (Sr%omega4sd (ib,ikibz,io,is)) *Ha_eV,& 544 REAL (Sr%sigxcme4sd(ib,ikibz,io,is)) *Ha_eV,& 545 AIMAG(Sr%sigxcme4sd(ib,ikibz,io,is)) *Ha_eV 546 end do 547 end do 548 !MRM 549 do ib=Sigp%minbnd(ikcalc,is),Sigp%maxbnd(ikcalc,is) 550 write(unt_sigc,'("# ik, ib",2i5)')ikibz,ib 551 do io=1,Sr%nomega4sd 552 write(unt_sigc,'(100(e12.5,2x))') & 553 REAL (Sr%omega4sd (ib,ikibz,io,is)) *Ha_eV,& 554 REAL (Sr%sigcme4sd(ib,ikibz,io,is)) *Ha_eV,& 555 AIMAG(Sr%sigcme4sd(ib,ikibz,io,is)) *Ha_eV 556 end do 557 end do 558 ! 559 if (mod10 == 1) then 560 ! For AC, write sigma matrix elements along the imaginary axis 561 do ib=Sigp%minbnd(ikcalc,is),Sigp%maxbnd(ikcalc,is) 562 write(unt_sgm,'("# ik, ib",2i5)')ikibz,ib 563 do io=1,Sr%nomega_i 564 write(unt_sgm,'(3(e12.5,2x))') & 565 AIMAG(Sr%omega_i(io)) *Ha_eV,& 566 REAL (Sr%sigxcmesi(ib,ikibz,io,is))*Ha_eV,& 567 AIMAG(Sr%sigxcmesi(ib,ikibz,io,is))*Ha_eV 568 end do 569 end do 570 end if 571 572 end do !is 573 574 end subroutine write_sigma_results