TABLE OF CONTENTS
- ABINIT/m_exc_diago
- m_exc_diago/exc_diago_coupling
- m_exc_diago/exc_diago_coupling_hegv
- m_exc_diago/exc_diago_driver
- m_exc_diago/exc_diago_resonant
- m_exc_diago/exc_print_eig
ABINIT/m_exc_diago [ Modules ]
NAME
m_exc_diago
FUNCTION
COPYRIGHT
Copyright (C) 2009-2024 ABINIT and EXC groups (L.Reining, V.Olevano, F.Sottile, S.Albrecht, G.Onida, M.Giantomassi) 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
15 #if defined HAVE_CONFIG_H 16 #include "config.h" 17 #endif 18 19 #include "abi_common.h" 20 21 MODULE m_exc_diago 22 23 use defs_basis 24 use m_bs_defs 25 use m_abicore 26 use m_errors 27 use m_xmpi 28 #if defined HAVE_MPI2 29 use mpi 30 #endif 31 use m_hdr 32 use m_sort 33 use m_slk 34 35 use defs_datatypes, only : pseudopotential_type, ebands_t 36 use m_io_tools, only : open_file 37 use m_fstrings, only : int2char4 38 use m_numeric_tools, only : print_arr, hermitianize 39 !use m_slk, only : matrix_scalapack, processor_scalapack 40 use m_crystal, only : crystal_t 41 use m_kpts, only : listkk 42 use m_bz_mesh, only : kmesh_t 43 use m_ebands, only : ebands_report_gap 44 use m_eprenorms, only : eprenorms_t 45 use m_wfd, only : wfdgw_t 46 use m_paw_hr, only : pawhur_t 47 use m_pawtab, only : pawtab_type 48 use m_exc_itdiago, only : exc_iterative_diago 49 use m_hide_lapack, only : xheev, xheevx, xgeev, xhegvx, xginv, xhdp_invert, xhegv 50 use m_hide_blas, only : xdotc, xgemm 51 use m_bse_io, only : exc_fullh_from_blocks, offset_in_file, rrs_of_glob, ccs_of_glob, & 52 & exc_read_bshdr, exc_skip_bshdr_mpio, exc_read_rblock_fio 53 use m_exc_spectra, only : build_spectra 54 55 implicit none 56 57 private 58 59 #if defined HAVE_MPI1 60 include 'mpif.h' 61 #endif 62 63 !#define DEV_MG_DEBUG_THIS 64 65 public :: exc_diago_driver ! Driver routine for the direct diagonalization of the BSE Hamiltonian (main entry point)
m_exc_diago/exc_diago_coupling [ Functions ]
[ Top ] [ m_exc_diago ] [ Functions ]
NAME
exc_diago_coupling
FUNCTION
Calculate excitonic eigenvalues and eigenvectors by performing a direct diagonalization. of the non Hermitian excitonic Hamiltoninan (resonant + coupling).
INPUTS
bseig_fname=The name of the output file.
Bsp
neh=Rank of the resonant block of the Hamiltoninan (equal to the rank of the coupling part)
comm=MPI communicator.
BS_files<excfiles>=Datatype storing names and files used in the Bethe-Salpeter code.
OUTPUT
Excitonic eigenvectors and eigenvalues are written on file BS_files%out_eig.
SOURCE
807 subroutine exc_diago_coupling(Bsp,BS_files,Hdr_bse,prtvol,comm) 808 809 !Arguments ------------------------------------ 810 !scalars 811 integer,intent(in) :: comm,prtvol 812 type(excfiles),intent(in) :: BS_files 813 type(excparam),intent(in) :: BSp 814 type(Hdr_type),intent(in) :: Hdr_bse 815 816 !Local variables ------------------------------ 817 !scalars 818 integer,parameter :: master=0,ldvl=1 819 integer(i8b) :: bsize_ham 820 integer :: ii,exc_size,hreso_unt,hcoup_unt,eig_unt,nsppol,nstates 821 integer :: bsz,block,bs1,bs2,jj 822 integer :: fform,row_sign 823 integer :: mi,it,nprocs,my_rank !itp 824 integer :: nene_printed,ierr 825 real(dp) :: exc_gap,exc_maxene,temp 826 logical :: diago_is_real,do_full_diago 827 logical :: do_ep_lifetime 828 character(len=500) :: msg 829 character(len=fnlen) :: hreso_fname,hcoup_fname,bseig_fname 830 !arrays 831 complex(dpc),allocatable :: exc_ham(:,:),exc_rvect(:,:),exc_ene(:),ovlp(:,:) 832 complex(dpc),allocatable :: cbuff(:,:) 833 complex(dpc) :: vl_dpc(ldvl,1) 834 835 !************************************************************************ 836 837 my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm) 838 839 nsppol = Hdr_bse%nsppol 840 if (nsppol==2) then 841 ABI_WARNING("nsppol==2 with coupling is still under development") 842 end if 843 844 if (nprocs > 1) then 845 ABI_WARNING("Scalapack does not provide ZGEEV, diagonalization is done in sequential!") 846 end if 847 848 exc_size = 2*SUM(BSp%nreh) 849 nstates = BSp%nstates 850 do_full_diago = (exc_size==nstates) 851 ABI_CHECK(do_full_diago,"Partial diago not coded yet") 852 853 bseig_fname = BS_files%out_eig 854 if (BS_files%in_eig /= BSE_NOFILE) then 855 ABI_ERROR("BS_files%in_eig is defined!") 856 end if 857 ! 858 ! Only master performs the diagonalization since ScaLAPACK does not provide the parallel version of ZGEEV. 859 if (my_rank/=master) GOTO 10 860 861 write(msg,'(a,i0)')' Direct diagonalization of the full excitonic Hamiltonian, Matrix size= ',exc_size 862 call wrtout([std_out, ab_out], msg) 863 864 bsize_ham = 2*dpc*exc_size**2 865 write(msg,'(a,f9.2,a)')' Allocating full excitonic Hamiltonian. Memory requested: ',bsize_ham*b2Gb,' Gb. ' 866 call wrtout(std_out, msg) 867 868 ABI_MALLOC_OR_DIE(exc_ham,(exc_size,exc_size), ierr) 869 870 write(msg,'(3a,f8.1,3a,f8.1,a)')& 871 ' Allocating excitonic eigenvalues and eigenvectors. ',ch10,& 872 ' Memory-space requested: ',2*dpc*exc_size*b2Gb,' Gb. ',ch10,& 873 ' Memory-space requested: ',bsize_ham*b2Gb,' Gb. ' 874 call wrtout(std_out, msg) 875 876 ABI_MALLOC_OR_DIE(exc_ene,(exc_size), ierr) 877 878 if (BS_files%in_hreso /= BSE_NOFILE) then 879 hreso_fname = BS_files%in_hreso 880 else 881 hreso_fname = BS_files%out_hreso 882 end if 883 884 call wrtout(std_out,' Reading resonant excitonic Hamiltonian from '//TRIM(hreso_fname)) 885 886 if (open_file(hreso_fname,msg,newunit=hreso_unt,form="unformatted",status="old",action="read") /= 0) then 887 ABI_ERROR(msg) 888 end if 889 ! 890 ! Read the header and perform consistency checks. 891 call exc_read_bshdr(hreso_unt,Bsp,fform,ierr) 892 ABI_CHECK(ierr==0,"Wrong header") 893 ! 894 ! Construct resonant and anti-resonant part of the excitonic Hamiltonian using Hermiticity. File is always in double precision. 895 ! Fill exc_ham with ( R 0 ) 896 ! ( 0 -R*) 897 !BEGINDEBUG 898 exc_ham = HUGE(one) 899 !ENDDEBUG 900 901 row_sign=-1; diago_is_real=(.not.BSp%have_complex_ene) 902 call exc_fullh_from_blocks(hreso_unt,"Resonant",nsppol,row_sign,diago_is_real,BSp%nreh,exc_size,exc_ham) 903 close(hreso_unt) 904 905 if (BS_files%in_hcoup /= BSE_NOFILE) then 906 hcoup_fname = BS_files%in_hcoup 907 else 908 hcoup_fname = BS_files%out_hcoup 909 end if 910 911 call wrtout(std_out,' Reading coupling excitonic Hamiltonian from '//TRIM(hcoup_fname)) 912 if (open_file(hcoup_fname,msg,newunit=hcoup_unt,form="unformatted",status="old",action="read") /= 0) then 913 ABI_ERROR(msg) 914 end if 915 ! 916 ! Read the header and perform consistency checks. 917 call exc_read_bshdr(hcoup_unt,Bsp,fform,ierr) 918 ABI_CHECK(ierr==0,"Wrong header") 919 ! 920 ! Fill exc_ham with ( 0 C) to have ( R C ) 921 ! (-C* 0) (-C* -R*) 922 row_sign=-1; diago_is_real=(.not.BSp%have_complex_ene) ! not used here 923 call exc_fullh_from_blocks(hcoup_unt,"Coupling",nsppol,row_sign,diago_is_real,BSp%nreh,exc_size,exc_ham) 924 925 !BEGINDEBUG 926 if (ANY(exc_ham==HUGE(one))) then 927 write(msg,'(a,2(1x,i0))')"There is a bug in exc_fullh_from_blocks",COUNT(exc_ham==HUGE(one)),exc_size**2 928 ABI_WARNING(msg) 929 bsz = Bsp%nreh(1) 930 ABI_MALLOC(cbuff,(bsz,bsz)) 931 block=0 932 do jj=1,2*nsppol 933 do ii=1,2*nsppol 934 block=block+1 935 bs1 = (ii-1)*bsz+1 936 bs2 = (jj-1)*bsz+1 937 cbuff = exc_ham(bs1:bs1+bsz-1,bs2:bs2+bsz-1) 938 if (ANY(cbuff==HUGE(one))) then 939 write(std_out,*)" for block ",ii,jj," found ",COUNT(cbuff==HUGE(one))," wrong entries" 940 end if 941 end do 942 end do 943 944 ABI_FREE(cbuff) 945 ABI_ERROR("Cannot continue") 946 end if 947 !ENDDEBUG 948 949 close(hcoup_unt) 950 ! 951 ! ====================================================== 952 ! ==== Calculate right eigenvectors and eigenvalues ==== 953 ! ====================================================== 954 ABI_MALLOC_OR_DIE(exc_rvect,(exc_size,exc_size), ierr) 955 956 if (do_full_diago) then 957 call wrtout(std_out,"Complete direct diagonalization with xgeev...") 958 call xgeev("No_left_eigen","Vectors",exc_size,exc_ham,exc_size,exc_ene,vl_dpc,ldvl,exc_rvect,exc_size) 959 else 960 ABI_ERROR("Not implemented error") 961 end if 962 963 ABI_FREE(exc_ham) 964 965 exc_gap = MINVAL(ABS(DBLE (exc_ene(1:nstates)))) 966 exc_maxene = MAXVAL(ABS(DBLE (exc_ene(1:nstates)))) 967 temp = MAXVAL(ABS(AIMAG(exc_ene(1:nstates)))) 968 969 write(msg,'(2(a,f7.2,2a),a,es9.2,2a)')& 970 " First excitonic eigenvalue: ",exc_gap*Ha_eV, " [eV].",ch10,& 971 " Last excitonic eigenvalue: ",exc_maxene*Ha_eV," [eV].",ch10,& 972 " Largest imaginary part: ",temp*Ha_eV, " [eV] ",ch10 973 call wrtout([std_out, ab_out], msg) 974 975 nene_printed = MIN(32*nsppol,nstates); if (prtvol>10) nene_printed = nstates 976 977 ! This is not portable as the the eigenvalues calculated by ZGEEV are not sorted. 978 ! Even two subsequent calculations with the same input on the same machine 979 ! might produce different orderings. Might sort the eigenvalues though, just for printing. 980 981 write(msg,'(a,i0)')' Complex excitonic eigenvalues in eV up to n= ',nene_printed 982 call wrtout(std_out,msg) 983 984 do it=0,(nene_printed-1)/4 985 write(msg,'(8f10.5)') ( exc_ene(ii)*Ha_eV, ii=1+it*4,MIN(it*4+4,nene_printed) ) 986 call wrtout(std_out,msg) 987 end do 988 989 call wrtout(std_out,ch10//" Writing eigenvalues and eigenvectors on file "//TRIM(bseig_fname)) 990 991 if (open_file(bseig_fname,msg,newunit=eig_unt,form="unformatted",action="write") /= 0) then 992 ABI_ERROR(msg) 993 end if 994 995 !YG : new version with lifetime 996 do_ep_lifetime = .FALSE. 997 write(eig_unt) do_ep_lifetime 998 999 write(eig_unt)exc_size,nstates 1000 write(eig_unt)CMPLX(exc_ene(1:nstates),kind=dpc) 1001 do mi=1,nstates 1002 write(eig_unt) exc_rvect(:,mi) 1003 end do 1004 1005 ABI_FREE(exc_ene) 1006 1007 ABI_MALLOC_OR_DIE(ovlp,(nstates,nstates), ierr) 1008 1009 call wrtout(std_out,' Calculating overlap matrix... ') 1010 1011 !do itp=1,nstates 1012 ! do it=1,nstates 1013 ! ovlp(it,itp) = xdotc(exc_size,exc_rvect(:,it),1,exc_rvect(:,itp),1) 1014 ! end do 1015 !end do 1016 call xgemm("C","N",exc_size,nstates,nstates,cone,exc_rvect,exc_size,exc_rvect,exc_size,czero,ovlp,nstates) 1017 ABI_FREE(exc_rvect) 1018 1019 call wrtout(std_out," Inverting overlap matrix... ") 1020 1021 ! Version for generic complex matrix. 1022 !call xginv(ovlp,exc_size) 1023 1024 ! The overlap is Hermitian definite positive. 1025 call xhdp_invert("Upper",ovlp,nstates) 1026 call hermitianize(ovlp,"Upper") 1027 1028 call wrtout(std_out,' Writing overlap matrix S^-1 on file: '//TRIM(bseig_fname)) 1029 1030 do it=1,nstates 1031 write(eig_unt) CMPLX(ovlp(:,it),kind=dpc) 1032 end do 1033 1034 ABI_FREE(ovlp) 1035 1036 close(eig_unt) 1037 1038 10 call xmpi_barrier(comm) 1039 1040 end subroutine exc_diago_coupling
m_exc_diago/exc_diago_coupling_hegv [ Functions ]
[ Top ] [ m_exc_diago ] [ Functions ]
NAME
exc_diago_coupling_hegv
FUNCTION
Calculate excitonic eigenvalues and eigenvectors by performing a direct diagonalization. of the non Hermitian excitonic Hamiltonian (resonant + coupling).
INPUTS
bseig_fname=The name of the output file.
Bsp
neh=Rank of the resonant block of the Hamiltoninan (equal to the rank of the coupling part)
comm=MPI communicator.
BS_files<excfiles>=Datatype storing names and files used in the Bethe-Salpeter code.
OUTPUT
Excitonic eigenvectors and eigenvalues are written to file BS_files%out_eig.
SOURCE
1065 subroutine exc_diago_coupling_hegv(Bsp,BS_files,Hdr_bse,prtvol,comm) 1066 1067 !Arguments ------------------------------------ 1068 !scalars 1069 integer,intent(in) :: comm,prtvol 1070 type(excparam),intent(in) :: BSp 1071 type(excfiles),intent(in) :: BS_files 1072 type(Hdr_type),intent(in) :: Hdr_bse 1073 1074 !Local variables ------------------------------ 1075 !scalars 1076 integer,parameter :: master=0 1077 integer(i8b) :: bsize_ham 1078 integer :: itype,il,iu,spin,row1,row2,pad_r1,pad_r2,neh1,neh2 1079 integer :: ii,exc_size,hreso_unt,hcoup_unt,eig_unt 1080 integer :: fform,neig_found,nstates 1081 integer :: mi,it,nprocs,my_rank 1082 integer :: nene_printed,nsppol,row_sign,ierr 1083 real(dp) :: exc_gap,exc_maxene,abstol,vl,vu 1084 character(len=500) :: msg 1085 character(len=fnlen) :: reso_fname,coup_fname,bseig_fname 1086 logical :: use_scalapack,do_full_diago,diago_is_real 1087 logical :: do_ep_lifetime 1088 !arrays 1089 real(dp),allocatable :: exc_ene(:) !,test_ene(:) 1090 complex(dpc),allocatable :: exc_ham(:,:),exc_rvect(:,:),fmat(:,:),ovlp(:,:) 1091 #if defined HAVE_LINALG_SCALAPACK && defined HAVE_MPI_IO 1092 integer,parameter :: istwfk1=1 1093 integer :: amode,mpi_fh,tbloc,mene_found,mpi_err,my_nel,nsblocks 1094 integer :: iloc,jloc,iglob,jglob,etype,slk_mask_type,offset_err,el,rrs_kind,ccs_kind 1095 !integer :: max_r,max_c 1096 integer(XMPI_OFFSET_KIND) :: ehdr_offset,fmarker,my_offset 1097 integer :: gsub(2,2) 1098 logical,parameter :: is_fortran_file=.TRUE. 1099 complex(dpc) :: ctmp 1100 integer,allocatable :: sub_block(:,:,:) 1101 integer,pointer :: myel2loc(:,:) 1102 complex(dpc),allocatable :: tmp_cbuffer(:) 1103 character(50) :: uplo 1104 real(dp),external :: PDLAMCH 1105 type(matrix_scalapack) :: Slk_F,Slk_Hbar,Slk_vec,Slk_ovlp !,Slk_tmp 1106 type(processor_scalapack) :: Slk_processor 1107 #endif 1108 1109 !************************************************************************ 1110 1111 my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm) 1112 1113 nsppol = Hdr_bse%nsppol 1114 if (nsppol==2) then 1115 ABI_WARNING("nsppol==2 is still under development!") 1116 end if 1117 1118 neh1 = BSp%nreh(1); neh2=neh1 1119 if (nsppol==2) neh2 = BSp%nreh(2) 1120 1121 exc_size = 2*SUM(Bsp%nreh) 1122 nstates = Bsp%nstates 1123 do_full_diago=(nstates==exc_size) 1124 1125 write(msg,'(a,i0)')'. Direct diagonalization of the full excitonic Hamiltonian, Matrix size= ',exc_size 1126 call wrtout([std_out, ab_out], msg) 1127 1128 bseig_fname = BS_files%out_eig 1129 if (BS_files%in_eig /= BSE_NOFILE) then 1130 ABI_ERROR("BS_files%in_eig is defined!") 1131 end if 1132 1133 if (BS_files%in_hreso /= BSE_NOFILE) then 1134 reso_fname = BS_files%in_hreso 1135 else 1136 reso_fname = BS_files%out_hreso 1137 end if 1138 call wrtout(std_out,' Reading resonant excitonic Hamiltonian from '//TRIM(reso_fname)) 1139 1140 if (BS_files%in_hcoup /= BSE_NOFILE) then 1141 coup_fname = BS_files%in_hcoup 1142 else 1143 coup_fname = BS_files%out_hcoup 1144 end if 1145 call wrtout(std_out,' Reading coupling excitonic Hamiltonian from '//TRIM(coup_fname)) 1146 1147 ! TODO: Reintegrate SCALAPACK: use new format 1148 ! --- !ERROR 1149 ! src_file: m_exc_diago.F90 1150 ! src_line: 1398 1151 ! mpi_rank: 0 1152 ! message: | 1153 ! SET_VIEW 1154 ! Other I/O error , error stack: 1155 ! ADIO_Set_view(48): **iobadoverlap displacements of filetype must be in a monotonically nondecreasing order 1156 ! ... 1157 1158 use_scalapack = .FALSE. 1159 #ifdef HAVE_LINALG_SCALAPACK 1160 ! This is alway false. I use this trick so that the second case below is always compiled 1161 ! to avoid regressions. 1162 use_scalapack = nprocs > 1 .and. nsppol > 5 1163 #endif 1164 !use_scalapack = .FALSE. 1165 !use_scalapack = .TRUE. 1166 1167 if (.not.use_scalapack .and. my_rank/=master) GOTO 10 1168 1169 ABI_MALLOC_OR_DIE(exc_ene,(exc_size), ierr) 1170 1171 SELECT CASE (use_scalapack) 1172 1173 CASE (.FALSE.) 1174 write(msg,'(a)')". Using LAPACK sequential version to solve FHv = ev with H positive definite. " 1175 call wrtout([std_out, ab_out], msg) 1176 1177 bsize_ham = 2*dpc*exc_size**2 1178 write(msg,'(a,f9.2,a)')' Allocating full excitonic Hamiltonian. Memory requested: ',2*bsize_ham*b2Gb,' Gb. ' 1179 call wrtout(std_out, msg) 1180 1181 ABI_MALLOC_OR_DIE(exc_ham,(exc_size,exc_size), ierr) 1182 ABI_MALLOC_OR_DIE(fmat,(exc_size,exc_size), ierr) 1183 1184 write(msg,'(3a,f8.1,3a,f8.1,a)')& 1185 ' Allocating excitonic eigenvalues and eigenvectors. ',ch10,& 1186 ' Memory-space requested: ',2*dpc*exc_size*b2Gb,' Gb. ',ch10,& 1187 ' Memory-space requested: ',bsize_ham*b2Gb,' Gb. ' 1188 call wrtout(std_out, msg) 1189 1190 if (open_file(reso_fname,msg,newunit=hreso_unt,form="unformatted",status="old",action="read") /= 0) then 1191 ABI_ERROR(msg) 1192 end if 1193 ! 1194 ! Read the header and perform consistency checks. 1195 call exc_read_bshdr(hreso_unt,Bsp,fform,ierr) 1196 ABI_CHECK(ierr==0,"Wrong header") 1197 ! 1198 ! Construct Hbar = ( R C ) 1199 ! ( C* R*) 1200 ! 1201 row_sign=+1; diago_is_real=(.not.BSp%have_complex_ene) 1202 call exc_fullh_from_blocks(hreso_unt,"Resonant",nsppol,row_sign,diago_is_real,Bsp%nreh,exc_size,exc_ham) 1203 close(hreso_unt) 1204 1205 if (open_file(coup_fname,msg,newunit=hcoup_unt,form="unformatted",status="old",action="read") /= 0) then 1206 ABI_ERROR(msg) 1207 end if 1208 ! 1209 ! Read the header and perform consistency checks. 1210 call exc_read_bshdr(hcoup_unt,Bsp,fform,ierr) 1211 ABI_CHECK(ierr==0,"Wrong header") 1212 1213 row_sign=+1; diago_is_real=(.not.BSp%have_complex_ene) ! not used here. 1214 call exc_fullh_from_blocks(hcoup_unt,"Coupling",nsppol,row_sign,diago_is_real,Bsp%nreh,exc_size,exc_ham) 1215 close(hcoup_unt) 1216 1217 !#ifdef DEV_MG_DEBUG_THIS 1218 !write(666)exc_ham 1219 !#endif 1220 ! 1221 ! Fill fmat = (1 0) 1222 ! (0 -1) 1223 fmat = czero 1224 do spin=1,nsppol 1225 pad_r1 = (spin-1)*Bsp%nreh(1) 1226 pad_r2 = SUM(Bsp%nreh) 1227 if (spin==2) pad_r2 = pad_r2 + Bsp%nreh(1) 1228 do it=1,Bsp%nreh(spin) 1229 row1 = it + pad_r1 1230 row2 = it + pad_r2 1231 fmat(row1,row1) = cone 1232 fmat(row2,row2) = -cone 1233 end do 1234 end do 1235 ! 1236 ! ================================================== 1237 ! ==== Solve generalized EV problem F H u = e u ==== 1238 ! ================================================== 1239 ! The eigenvectors Z are normalized as follows: if ITYPE = 1 or 2, Z**T*B*Z = I; if ITYPE = 3, Z**T*inv(B)*Z = I. 1240 ! 1241 itype=2 1242 if (do_full_diago) then 1243 call wrtout(std_out," Full diagonalization with XHEGV... ") 1244 call xhegv(itype,"Vectors","Upper",exc_size,fmat,exc_ham,exc_ene) 1245 else 1246 call wrtout(std_out," Partial diagonalization with XHEGVX... ") 1247 ABI_MALLOC_OR_DIE(exc_rvect,(exc_size,nstates), ierr) 1248 il=1; iu=1; abstol=zero 1249 call xhegvx(itype,"Vectors","All","Upper",exc_size,fmat,exc_ham,vl,vu,il,iu,abstol,neig_found,exc_ene,exc_rvect,exc_size) 1250 end if 1251 1252 ABI_FREE(exc_ham) 1253 1254 if (do_full_diago) then 1255 ABI_MALLOC(exc_rvect,(exc_size,nstates)) 1256 exc_rvect = fmat(:,1:nstates) 1257 end if 1258 1259 ABI_FREE(fmat) 1260 1261 call wrtout(std_out," Writing eigenvalues and eigenvectors on file: "//TRIM(bseig_fname)) 1262 1263 if (open_file(bseig_fname,msg,newunit=eig_unt,form="unformatted",action="write") /= 0) then 1264 ABI_ERROR(msg) 1265 end if 1266 1267 do_ep_lifetime = .FALSE. 1268 write(eig_unt) do_ep_lifetime 1269 write(eig_unt) exc_size, nstates 1270 write(eig_unt) CMPLX(exc_ene(1:nstates),kind=dpc) 1271 do mi=1,nstates 1272 write(eig_unt) CMPLX(exc_rvect(:,mi),kind=dpc) 1273 end do 1274 1275 !#ifdef DEV_MG_DEBUG_THIS 1276 ! write(888)exc_rvect 1277 ! write(888)exc_ene 1278 !#endif 1279 1280 ABI_MALLOC_OR_DIE(ovlp,(nstates,nstates), ierr) 1281 1282 call wrtout(std_out,' Calculating overlap matrix...') 1283 1284 call xgemm("C","N",exc_size,nstates,nstates,cone,exc_rvect,exc_size,exc_rvect,exc_size,czero,ovlp,nstates) 1285 ABI_FREE(exc_rvect) 1286 1287 !#ifdef DEV_MG_DEBUG_THIS 1288 !write(667)ovlp 1289 !#endif 1290 1291 call wrtout(std_out," Inverting overlap matrix... ") 1292 ! 1293 ! The overlap is Hermitian definite positive. 1294 call xhdp_invert("Upper",ovlp,nstates) 1295 call hermitianize(ovlp,"Upper") 1296 1297 ! Version for generic complex matrix. 1298 !call xginv(ovlp,nstates) 1299 1300 !#ifdef DEV_MG_DEBUG_THIS 1301 !write(668,*)ovlp 1302 !#endif 1303 1304 call wrtout(std_out,' Writing overlap matrix O^-1 on file: '//TRIM(bseig_fname)) 1305 do it=1,nstates 1306 write(eig_unt) ovlp(:,it) 1307 end do 1308 1309 ABI_FREE(ovlp) 1310 close(eig_unt) 1311 1312 CASE (.TRUE.) 1313 1314 #if defined HAVE_LINALG_SCALAPACK && defined HAVE_MPI_IO 1315 ! 1316 ! Init scaLAPACK matrix Hbar = ( R C ) 1317 ! ( C* R*) 1318 ! Battle plan: 1319 ! Here the reading is complicated by the fact that R and C are stored on two different files 1320 ! and moreover the matrices are in packed storage mode. 1321 ! For initializing the local part of the resonant and anti-resonant block we have to allocate 1322 ! a temporary buffer. Then we read the buffer from file and the corresponding elements of the 1323 ! scaLAPACK matrix are initialized taking into account the symmetries of R (Hermitian) 1324 ! The same procedure is used to read the coupling and the anti-coupling part (Symmetric). 1325 ! 1326 tbloc=50 1327 write(msg,'(2(a,i0))')". Using MPI-IO + scaLAPACK version with nprocs= ",nprocs,"; block size= ",tbloc 1328 call wrtout([std_out, ab_out], msg, do_flush=.True.) 1329 ! 1330 ! Init scaLAPACK environment. 1331 call Slk_processor%init(comm) 1332 ! 1333 ! Open the Resonant file with MPI-IO and skip the record. 1334 amode=MPI_MODE_RDONLY 1335 1336 call MPI_FILE_OPEN(comm, reso_fname, amode, MPI_INFO_NULL, mpi_fh, mpi_err) 1337 msg = " MPI_IO error opening file: "//TRIM(reso_fname) 1338 ABI_CHECK_MPI(mpi_err,msg) 1339 ! 1340 ! Skip the header and find the offset for reading the matrix. 1341 call exc_skip_bshdr_mpio(mpi_fh,xmpio_collective,ehdr_offset) 1342 ! 1343 ! Read = ( R - ) 1344 ! ( - R*) 1345 call Slk_Hbar%init(exc_size,exc_size,Slk_processor,istwfk1) 1346 1347 nullify(myel2loc) 1348 nsblocks=nsppol 1349 ABI_MALLOC(sub_block,(2,2,nsblocks)) 1350 ABI_CHECK(nsppol==1,"nsppol==2 not coded yet") 1351 1352 call slk_single_fview_read_mask(Slk_Hbar,rrs_of_glob,offset_in_file,nsblocks,sub_block,my_nel,myel2loc,etype,slk_mask_type,& 1353 & offset_err,is_fortran_file) 1354 1355 if (offset_err/=0) then 1356 write(msg,"(3a)")& 1357 & " Global position index cannot be stored in a standard Fortran integer ",ch10,& 1358 & " Excitonic matrix cannot be read with a single MPI-IO call." 1359 ABI_ERROR(msg) 1360 end if 1361 1362 ! Shift the offset because the view starts at the fist matrix element! 1363 ! TODO should rationalize the treatment of the offset 1364 my_offset = ehdr_offset + xmpio_bsize_frm 1365 call MPI_FILE_SET_VIEW(mpi_fh, my_offset, etype, slk_mask_type, 'native', MPI_INFO_NULL, mpi_err) 1366 ABI_CHECK_MPI(mpi_err,"SET_VIEW") 1367 1368 call MPI_TYPE_FREE(slk_mask_type,mpi_err) 1369 ABI_CHECK_MPI(mpi_err,"MPI_TYPE_FREE") 1370 ! 1371 ! Read my portion of the R,-R* sublocks and store the values in a temporary buffer. 1372 ABI_MALLOC_OR_DIE(tmp_cbuffer,(my_nel), ierr) 1373 1374 call xmpi_barrier(comm) 1375 1376 call MPI_FILE_READ_ALL(mpi_fh, tmp_cbuffer, my_nel, MPI_DOUBLE_COMPLEX, MPI_STATUS_IGNORE, mpi_err) 1377 ABI_CHECK_MPI(mpi_err,"READ_ALL") 1378 ! 1379 ! Symmetrize my Resonant part. 1380 do el=1,my_nel 1381 iloc = myel2loc(1,el) 1382 jloc = myel2loc(2,el) 1383 call Slk_Hbar%loc2glob(iloc,jloc,iglob,jglob) 1384 ctmp = tmp_cbuffer(el) 1385 if (iglob==jglob.and..not.Bsp%have_complex_ene) ctmp = DBLE(ctmp) ! Force the diagonal to be real. 1386 rrs_kind = rrs_of_glob(iglob,jglob,Slk_Hbar%sizeb_global) 1387 if (rrs_kind==1.and.jglob<iglob) then ! Lower resonant 1388 ctmp = DCONJG(ctmp) 1389 else if (rrs_kind==-1.and.jglob>=iglob) then ! Lower Anti-resonant (Diagonal is included). 1390 ctmp = DCONJG(ctmp) 1391 end if 1392 Slk_Hbar%buffer_cplx(iloc,jloc) = ctmp 1393 end do 1394 1395 ABI_FREE(tmp_cbuffer) 1396 ABI_FREE(myel2loc) 1397 1398 call MPI_FILE_CLOSE(mpi_fh, mpi_err) 1399 ABI_CHECK_MPI(mpi_err,"FILE_CLOSE") 1400 ! 1401 ! Read = ( - C) 1402 ! (-C* -) 1403 ! 1404 call MPI_FILE_OPEN(comm, coup_fname, amode, MPI_INFO_NULL, mpi_fh, mpi_err) 1405 msg = " MPI_IO error opening file: "//TRIM(coup_fname) 1406 ABI_CHECK_MPI(mpi_err,msg) 1407 ! 1408 ! Skip the header and find the offset for reading the matrix. 1409 call exc_skip_bshdr_mpio(mpi_fh,xmpio_collective,ehdr_offset) 1410 1411 nullify(myel2loc) 1412 call slk_single_fview_read_mask(Slk_Hbar,ccs_of_glob,offset_in_file,nsblocks,sub_block,my_nel,myel2loc,etype,slk_mask_type,& 1413 & offset_err,is_fortran_file) 1414 1415 ABI_FREE(sub_block) 1416 1417 if (offset_err/=0) then 1418 write(msg,"(3a)")& 1419 & " Global position index cannot be stored in a standard Fortran integer ",ch10,& 1420 & " Excitonic matrix cannot be read with a single MPI-IO call." 1421 ABI_ERROR(msg) 1422 end if 1423 ! 1424 ! Shift the offset because the view starts at the fist matrix element! 1425 ! TODO should rationalize the treatment of the offset so that the client code 1426 ! will automatically receive my_offset. 1427 my_offset = ehdr_offset + xmpio_bsize_frm 1428 call MPI_FILE_SET_VIEW(mpi_fh, my_offset, etype, slk_mask_type, 'native', MPI_INFO_NULL, mpi_err) 1429 ABI_CHECK_MPI(mpi_err,"SET_VIEW") 1430 1431 call MPI_TYPE_FREE(slk_mask_type,mpi_err) 1432 ABI_CHECK_MPI(mpi_err,"MPI_TYPE_FREE") 1433 ! 1434 ! Read my portion of the C-C* blocks and store the values in a temporary buffer. 1435 ABI_MALLOC_OR_DIE(tmp_cbuffer,(my_nel), ierr) 1436 1437 call MPI_FILE_READ_ALL(mpi_fh, tmp_cbuffer, my_nel, MPI_DOUBLE_COMPLEX, MPI_STATUS_IGNORE, mpi_err) 1438 ABI_CHECK_MPI(mpi_err,"READ_ALL") 1439 ! 1440 ! Symmetrize my coupling part. 1441 ! Coupling block is symmetric => No symmetrization of the lower triangle. 1442 do el=1,my_nel 1443 iloc = myel2loc(1,el) 1444 jloc = myel2loc(2,el) 1445 call Slk_Hbar%loc2glob(iloc, jloc, iglob, jglob) 1446 ccs_kind = ccs_of_glob(iglob,jglob,Slk_Hbar%sizeb_global) 1447 ctmp = tmp_cbuffer(el) 1448 if (ccs_kind==-1) ctmp = DCONJG(ctmp) ! Anti-coupling (Diagonal is included). 1449 Slk_Hbar%buffer_cplx(iloc,jloc) = ctmp 1450 end do 1451 1452 ABI_FREE(tmp_cbuffer) 1453 ABI_FREE(myel2loc) 1454 1455 !max_r=20; max_c=10 1456 !call print_arr(Slk_Hbar%buffer_cplx,max_r=max_r,max_c=max_c,unit=std_out) 1457 1458 !#ifdef DEV_MG_DEBUG_THIS 1459 ! ABI_MALLOC(exc_ham,(exc_size,exc_size)) 1460 ! read(666)exc_ham 1461 ! 1462 ! write(std_out,*)"Error Hbar: ",MAXVAL(ABS(exc_ham-Slk_Hbar%buffer_cplx)) 1463 ! ABI_FREE(exc_ham) 1464 !#endif 1465 1466 call MPI_FILE_CLOSE(mpi_fh, mpi_err) 1467 ABI_CHECK_MPI(mpi_err,"FILE_CLOSE") 1468 ! 1469 ! Init scaLAPACK matrix F 1470 call Slk_F%init(exc_size,exc_size,Slk_processor,istwfk1) 1471 ! 1472 ! Global F = (1 0) 1473 ! (0 -1) 1474 do jloc=1,Slk_F%sizeb_local(2) 1475 do iloc=1,Slk_F%sizeb_local(1) 1476 call Slk_F%loc2glob(iloc, jloc, iglob, jglob) 1477 if (iglob==jglob) then 1478 if (iglob<=SUM(Bsp%nreh)) then 1479 Slk_F%buffer_cplx(iloc,jloc) = cone 1480 else 1481 Slk_F%buffer_cplx(iloc,jloc) = -cone 1482 end if 1483 else 1484 Slk_F%buffer_cplx(iloc,jloc) = czero 1485 end if 1486 end do 1487 end do 1488 ! 1489 ! =========================================================== 1490 ! ==== Solve generalized EV problem H u = F Hbar u = e u ==== 1491 ! =========================================================== 1492 call Slk_vec%init(exc_size,exc_size,Slk_processor,istwfk1) 1493 ! 1494 itype=2; vl=1; vu=1; il=1; iu=nstates 1495 abstol=zero !ABSTOL = PDLAMCH(comm,'U') 1496 1497 !#if 1 1498 if (do_full_diago) then 1499 call slk_F%pzhegvx(itype,"Vectors","All","Upper",Slk_Hbar,vl,vu,il,iu,abstol,Slk_vec,mene_found,exc_ene) 1500 else 1501 ABI_WARNING("Partial diago is still under testing") 1502 call slk_F%pzhegvx(itype,"Vectors","Index","Upper",Slk_Hbar,vl,vu,il,iu,abstol,Slk_vec,mene_found,exc_ene) 1503 end if 1504 !#else 1505 ! call xhegv(itype,"Vectors","Upper",exc_size,Slk_F%buffer_cplx,Slk_Hbar%buffer_cplx,exc_ene) 1506 ! Slk_vec%buffer_cplx = Slk_F%buffer_cplx 1507 !#endif 1508 1509 !#ifdef DEV_MG_DEBUG_THIS 1510 ! if (PRODUCT(Slk_Hbar%sizeb_local) /= exc_size**2) then 1511 ! ABI_ERROR("Wrong size") 1512 ! end if 1513 ! 1514 ! ABI_MALLOC(exc_ham,(exc_size,exc_size)) 1515 ! read(888)exc_ham 1516 ! 1517 ! write(std_out,*)"Error rvec: ",MAXVAL(ABS(exc_ham-Slk_vec%buffer_cplx)) 1518 ! ABI_FREE(exc_ham) 1519 ! 1520 ! ABI_MALLOC(test_ene,(exc_size)) 1521 ! read(888)test_ene 1522 ! write(std_out,*)"Error ene: ",MAXVAL(ABS(exc_ene-test_ene)) 1523 ! ABI_FREE(test_ene) 1524 !#endif 1525 1526 call Slk_F%free() 1527 call Slk_Hbar%free() 1528 1529 call wrtout(std_out,ch10//" Writing eigenvalues and eigenvectors on file: "//TRIM(bseig_fname)) 1530 ! 1531 ! Open the file with Fortran-IO to write the Header. 1532 if (my_rank==master) then 1533 if (open_file(bseig_fname,msg,newunit=eig_unt,form="unformatted",action="write") /= 0) then 1534 ABI_ERROR(msg) 1535 end if 1536 1537 write(eig_unt) exc_size,nstates 1538 write(eig_unt) CMPLX(exc_ene(1:nstates),kind=dpc) 1539 !do mi=1,exc_size 1540 ! write(eig_unt) CMPLX(exc_rvect(:,mi),kind=dpc) 1541 !end do 1542 close(eig_unt) 1543 end if 1544 ! 1545 ! Open the file with MPI-IO and write the distributed eigevectors. 1546 call xmpi_barrier(comm) 1547 amode=MPI_MODE_RDWR 1548 call MPI_FILE_OPEN(comm, bseig_fname, amode, MPI_INFO_NULL, mpi_fh, mpi_err) 1549 ABI_CHECK_MPI(mpi_err,"FILE_OPEN: "//TRIM(bseig_fname)) 1550 ! 1551 ! Skip the header and find the offset for writing the matrix. 1552 ehdr_offset=0 1553 !call hdr_mpio_skip(mpi_fh,fform,ehdr_offset) 1554 1555 call xmpio_read_frm(mpi_fh,ehdr_offset,xmpio_collective,fmarker,mpi_err) 1556 write(std_out,*)" fmarker1 = ",fmarker 1557 call xmpio_read_frm(mpi_fh,ehdr_offset,xmpio_collective,fmarker,mpi_err) 1558 write(std_out,*)" fmarker2 = ",fmarker 1559 1560 write(std_out,*)" Writing nstates ",nstates 1561 gsub(:,1) = (/1,1/) 1562 gsub(:,2) = (/exc_size,nstates/) 1563 call slk_write(Slk_vec,"All",is_fortran_file,mpi_fh=mpi_fh,offset=ehdr_offset,glob_subarray=gsub) 1564 1565 call wrtout(std_out,' Calculating overlap matrix... ') 1566 if (.not.do_full_diago) then 1567 ABI_ERROR(" Init of Slk_ovlp is wrong") 1568 end if 1569 1570 call Slk_ovlp%init(exc_size,exc_size,Slk_processor,istwfk1) 1571 1572 ! Calculate the overlap matrix. 1573 ! FIXME 1574 ! The ESLL manual says that "matrices matrix1 and matrix2 must have no common elements; 1575 ! otherwise, results are unpredictable." 1576 ! However the official scaLAPACK documentation does not report this (severe) limitation. 1577 1578 !call Slk_tmp%init(exc_size,exc_size,Slk_processor,istwfk1) 1579 !Slk_tmp%buffer_cplx = Slk_vec%buffer_cplx 1580 !call slk_pgemm("C","N",Slk_tmp,cone,Slk_vec,czero,Slk_ovlp) 1581 !call Slk_tmp%free() 1582 1583 call slk_pgemm("C","N",Slk_vec,cone,Slk_vec,czero,Slk_ovlp) 1584 1585 !#ifdef DEV_MG_DEBUG_THIS 1586 ! ABI_MALLOC(exc_ham,(exc_size,exc_size)) 1587 ! read(667)exc_ham 1588 ! 1589 ! write(std_out,*)"Error Ovlp: ",MAXVAL(ABS(exc_ham-Slk_ovlp%buffer_cplx)) 1590 ! !Slk_ovlp%buffer_cplx = exc_ham 1591 !#endif 1592 1593 !max_r=20; max_c=10 1594 !call print_arr(Slk_ovlp%buffer_cplx,max_r=max_r,max_c=max_c,unit=std_out) 1595 1596 call Slk_vec%free() 1597 1598 call wrtout(std_out," Inverting overlap matrix... ") 1599 uplo="Upper" 1600 1601 !#if 0 1602 !!DEBUG 1603 ! call xhdp_invert(uplo,Slk_ovlp%buffer_cplx,exc_size) 1604 ! 1605 ! !call slk_symmetrize(Slk_ovlp,uplo,"Hermitian") 1606 ! call hermitianize(Slk_ovlp%buffer_cplx,uplo) 1607 ! 1608 ! exc_ham = MATMUL(exc_ham,Slk_ovlp%buffer_cplx) 1609 ! do it=1,exc_size 1610 ! exc_ham(it,it) = exc_ham(it,it) - cone 1611 ! end do 1612 ! 1613 ! write(std_out,*)"Error Inversion: ",MAXVAL(ABS(exc_ham)) 1614 ! ABI_FREE(exc_ham) 1615 !!END DEBUG 1616 ! 1617 !#else 1618 ! call Slk_ovlp%hpd_invert(uplo) 1619 ! call hermitianize(Slk_ovlp%buffer_cplx,uplo) 1620 ! !call slk_symmetrize(Slk_ovlp,uplo,"Hermitian") 1621 1622 call Slk_ovlp%invert() ! Version for generic complex matrix. 1623 !#endif 1624 1625 !#ifdef DEV_MG_DEBUG_THIS 1626 ! ABI_MALLOC(exc_ham,(exc_size,exc_size)) 1627 ! read(668)exc_ham 1628 ! write(std_out,*)"Error in Inv Ovlp: ",MAXVAL(ABS(exc_ham-Slk_ovlp%buffer_cplx)) 1629 ! 1630 ! !exc_ham = exc_ham-Slk_ovlp%buffer_cplx 1631 ! !do it=1,exc_size 1632 ! ! if ( MAXVAL(ABS(exc_ham(:,it))) > 0.1 ) write(std_out,*)"it: ",it,exc_ham(:,it) 1633 ! !end do 1634 ! 1635 ! !Slk_ovlp%buffer_cplx = exc_ham 1636 ! ABI_FREE(exc_ham) 1637 ! 1638 ! !write(std_out,*)"MAX ERR",MAXVAL(ABS(Slk_ovlp%buffer_cplx - TRANSPOSE(DCONJG(Slk_ovlp%buffer_cplx)))) 1639 !#endif 1640 1641 call wrtout(std_out,' Writing overlap matrix S^-1 on file: '//TRIM(bseig_fname)) 1642 1643 call slk_write(Slk_ovlp,"All",is_fortran_file,mpi_fh=mpi_fh,offset=ehdr_offset) 1644 1645 call MPI_FILE_CLOSE(mpi_fh, mpi_err) 1646 ABI_CHECK_MPI(mpi_err,"FILE_CLOSE") 1647 1648 call Slk_ovlp%free() 1649 call Slk_processor%free() 1650 #else 1651 ABI_BUG("You should not be here!") 1652 #endif 1653 1654 END SELECT 1655 1656 exc_gap = MINVAL(ABS(exc_ene(1:nstates))) 1657 exc_maxene = MAXVAL(ABS(exc_ene(1:nstates))) 1658 1659 write(msg,'(2(a,f7.2,2a))')& 1660 " First excitonic eigenvalue: ",exc_gap*Ha_eV, " [eV].",ch10,& 1661 " Last excitonic eigenvalue: ",exc_maxene*Ha_eV," [eV].",ch10 1662 call wrtout([std_out, ab_out], msg) 1663 1664 nene_printed = MIN(32*nsppol,nstates); if (prtvol>10) nene_printed = nstates 1665 write(msg,'(a,i0)')' Complex excitonic eigenvalues in eV up to n= ',nene_printed 1666 call wrtout(std_out, msg) 1667 1668 do it=0,(nene_printed-1)/4 1669 write(msg,'(4f10.5)') ( exc_ene(ii)*Ha_eV, ii=1+it*4,MIN(it*4+4,nene_printed) ) 1670 call wrtout(std_out, msg) 1671 end do 1672 1673 ABI_FREE(exc_ene) 1674 1675 10 call xmpi_barrier(comm) 1676 1677 end subroutine exc_diago_coupling_hegv
m_exc_diago/exc_diago_driver [ Functions ]
[ Top ] [ m_exc_diago ] [ Functions ]
NAME
exc_diago_driver
FUNCTION
Driver routine for the direct diagonalization of the Hermitian excitonic Hamiltonian.
INPUTS
neh=Rank of the resonant block of the Hamiltonian. BS_files<excfiles>=Datatype storing names and files used in the Bethe-Salpeter code. %exh=Name of the file storing the excitonic resonant part.
OUTPUT
Eigenvalues and eigenvectors are written on file.
SOURCE
87 subroutine exc_diago_driver(Wfd,Bsp,BS_files,KS_BSt,QP_BSt,Cryst,Kmesh,Psps,& 88 & Pawtab,Hur,Hdr_bse,drude_plsmf,Epren) 89 90 !Arguments ------------------------------------ 91 !scalars 92 real(dp),intent(in) :: drude_plsmf 93 type(excparam),intent(in) :: BSp 94 type(excfiles),intent(in) :: BS_files 95 type(Hdr_type),intent(in) :: Hdr_bse 96 type(crystal_t),intent(in) :: Cryst 97 type(pseudopotential_type),intent(in) :: Psps 98 type(kmesh_t),intent(in) :: Kmesh 99 type(ebands_t),intent(in) :: KS_BSt,QP_BSt 100 type(wfdgw_t),intent(inout) :: Wfd 101 type(eprenorms_t),intent(in) :: Epren 102 !arrays 103 type(pawtab_type),intent(in) :: Pawtab(Cryst%ntypat*Wfd%usepaw) 104 type(pawhur_t),intent(in) :: Hur(Cryst%natom*Wfd%usepaw) 105 106 !Local variables ------------------------------ 107 !scalars 108 integer :: my_rank,master,comm,prtvol 109 complex(dpc) :: exc_gap,gw_gap 110 logical :: eval_eigenstates 111 character(len=500) :: msg 112 !arrays 113 real(dp) :: gaps(3,QP_BSt%nsppol) 114 115 !************************************************************************ 116 117 DBG_ENTER("COLL") 118 119 comm = Wfd%comm 120 my_rank = Wfd%my_rank 121 master = Wfd%master 122 prtvol = Wfd%prtvol 123 124 if (BSp%have_complex_ene) then 125 ABI_ERROR("Complex energies are not supported yet") 126 end if 127 ! 128 ! This trick is needed to restart a CG run, use DDIAGO to calculate the spectra reusing an old BSEIG file. 129 eval_eigenstates = (BS_files%in_eig == BSE_NOFILE) .or. (Bsp%algorithm == BSE_ALGO_CG) 130 131 if (eval_eigenstates) then 132 ! 133 select case (BSp%algorithm) 134 case (BSE_ALGO_DDIAGO) 135 if (BSp%use_coupling==0) then 136 call exc_diago_resonant(BSp,BS_files,Hdr_bse,prtvol,comm) 137 if(Bsp%do_ep_renorm) then 138 call exc_diago_resonant(BSp,BS_files,Hdr_bse,prtvol,comm,Epren=Epren,Kmesh=Kmesh,Cryst=Cryst,elph_lifetime=.TRUE.) 139 end if 140 else 141 if (Bsp%have_complex_ene) then 142 ! Solve Hv = ev with generic complex matrix. 143 call exc_diago_coupling(BSp,BS_files,Hdr_bse,prtvol,comm) 144 else 145 ! Solve generalized eigenvalue problem F Hbar with Hbar Hermitian definitive positive matrix. 146 call exc_diago_coupling_hegv(BSp,BS_files,Hdr_bse,prtvol,comm) 147 end if 148 end if 149 150 case (BSE_ALGO_CG) 151 if (BSp%use_coupling==0) then 152 call exc_iterative_diago(Bsp,BS_files,Hdr_bse,prtvol,comm) 153 else 154 ABI_ERROR("CG + coupling not coded") 155 end if 156 157 case default 158 write(msg,'(a,i0)')" Wrong value for Bsp%algorithm: ",Bsp%algorithm 159 ABI_ERROR(msg) 160 end select 161 ! 162 if (my_rank==master) then 163 call ebands_report_gap(QP_BSt,header="QP bands",unit=std_out,gaps=gaps) 164 gw_gap = MINVAL(gaps(2,:)) 165 call exc_print_eig(BSp,BS_files%out_eig,gw_gap,exc_gap) 166 end if 167 call xmpi_barrier(comm) 168 ! 169 end if 170 171 call build_spectra(BSp,BS_files,Cryst,Kmesh,KS_BSt,QP_BSt,Psps,Pawtab,Wfd,Hur,drude_plsmf,comm) 172 173 ! Electron-phonon renormalization ! 174 if (BSp%algorithm == BSE_ALGO_DDIAGO .and. BSp%use_coupling == 0 .and. BSp%do_ep_renorm) then 175 call build_spectra(BSp,BS_files,Cryst,Kmesh,KS_BSt,QP_BSt,Psps,Pawtab,Wfd,Hur,drude_plsmf,comm,Epren=Epren) 176 end if 177 178 DBG_EXIT("COLL") 179 180 end subroutine exc_diago_driver
m_exc_diago/exc_diago_resonant [ Functions ]
[ Top ] [ m_exc_diago ] [ Functions ]
NAME
exc_diago_resonant
FUNCTION
Calculates eigenvalues and eigenvectors of the Hermitian excitonic Hamiltonian (coupling is neglected).
INPUTS
Bsp BS_files<excfiles>=Datatype storing names and files used in the Bethe-Salpeter code. comm=MPI communicator. bseig_fname=The name of the output file prtvol=Verbosity level.
OUTPUT
Eigenvalues and eigenvectors are written on file bseig_fname
SOURCE
204 subroutine exc_diago_resonant(Bsp,BS_files,Hdr_bse,prtvol,comm,Epren,Kmesh,Cryst,elph_lifetime) 205 206 !Arguments ------------------------------------ 207 !scalars 208 integer,intent(in) :: comm,prtvol 209 logical,optional,intent(in) :: elph_lifetime 210 type(excparam),intent(in) :: BSp 211 type(excfiles),intent(in) :: BS_files 212 type(Hdr_type),intent(in) :: Hdr_bse 213 type(eprenorms_t),optional,intent(in) :: Epren 214 type(kmesh_t),optional,intent(in) :: Kmesh 215 type(crystal_t),optional,intent(in) :: Cryst 216 217 !Local variables ------------------------------ 218 !scalars 219 integer,parameter :: master=0 220 integer :: ii,it,mi,hreso_unt,eig_unt,exc_size,neh1,neh2,j 221 integer :: nsppol,il,iu,mene_found,nstates 222 integer :: nprocs,my_rank,fform,nene_printed,ierr 223 real(dp) :: exc_gap,exc_maxene,abstol 224 real(dp) :: vl,vu 225 logical :: use_scalapack,do_full_diago,diagonal_is_real 226 character(len=500) :: msg 227 character(len=fnlen) :: hreso_fname,bseig_fname 228 !arrays 229 real(dp),allocatable :: exc_ene(:) 230 complex(dpc),allocatable :: exc_mat(:,:),exc_vec(:,:) 231 #if defined HAVE_LINALG_SCALAPACK && defined HAVE_MPI_IO 232 integer :: amode,mpi_fh,istwf_k,tbloc 233 integer(XMPI_OFFSET_KIND) :: ehdr_offset,fmarker 234 integer :: block_sizes(2,3),array_of_sizes(2),gsub(2,2) 235 logical,parameter :: is_fortran_file=.TRUE. 236 real(dp),external :: PDLAMCH 237 type(matrix_scalapack) :: Slk_mat,Slk_vec 238 type(processor_scalapack) :: Slk_processor 239 #endif 240 241 integer :: ik, ic, iv, isppol, ireh, ep_ik, itemp 242 complex(dpc) :: en 243 244 real(dp) :: dksqmax 245 integer,allocatable :: bs2eph(:,:) 246 integer :: sppoldbl, timrev 247 logical :: do_ep_renorm, do_ep_lifetime 248 integer :: ntemp 249 character(len=4) :: ts 250 complex(dpc),allocatable :: exc_vl(:,:),exc_ene_c(:) 251 complex(dpc) :: ctemp 252 !! complex(dpc),allocatable :: ovlp(:,:) 253 254 !************************************************************************ 255 256 DBG_ENTER("PERS") 257 258 nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm) 259 260 if (BSp%have_complex_ene) then ! QP lifetimes are not included 261 ABI_ERROR("complex energies not coded yet") 262 end if 263 264 if (ANY(Bsp%nreh/=Bsp%nreh(1))) then 265 write(std_out,*)" Bsp%nreh: ",Bsp%nreh 266 write(msg,'(a)')" BSE code does not support different number of transitions for the two spin channels" 267 ABI_WARNING(msg) 268 end if 269 270 nsppol = Hdr_bse%nsppol 271 exc_size = SUM(BSp%nreh) 272 nstates = BSp%nstates; do_full_diago=(Bsp%nstates==exc_size) 273 274 neh1 = Bsp%nreh(1); neh2 = neh1 275 if (Hdr_bse%nsppol==2) neh2 = Bsp%nreh(2) 276 277 ! Scalapack is disabled due to portability issues in slk_read 278 ! This part should be rewritten with hdf5 + mpi-io 279 280 use_scalapack = .FALSE. 281 !#if defined HAVE_LINALG_SCALAPACK 282 ! use_scalapack = (nprocs > 1) 283 !#endif 284 if (use_scalapack .and. nsppol == 2) then 285 use_scalapack = .False. 286 msg = "Scalapack with nsppol==2 not yet available. Using sequential version" 287 ABI_WARNING(msg) 288 end if 289 290 if (.not.use_scalapack .and. my_rank/=master) GOTO 10 ! Inversion is done by master only. 291 292 nene_printed = MIN(32*nsppol,nstates); if (prtvol>10) nene_printed = nstates 293 294 if (BS_files%in_hreso /= BSE_NOFILE) then 295 hreso_fname = BS_files%in_hreso 296 else 297 hreso_fname = BS_files%out_hreso 298 end if 299 300 bseig_fname = BS_files%out_eig 301 if (BS_files%in_eig /= BSE_NOFILE) then 302 ABI_ERROR("BS_files%in_eig is defined!") 303 end if 304 305 write(msg,'(a,i0)')' Direct diagonalization of the resonant excitonic Hamiltonian, Matrix size= ',exc_size 306 call wrtout([std_out, ab_out], msg) 307 308 ABI_MALLOC_OR_DIE(exc_ene,(exc_size), ierr) 309 310 ABI_MALLOC_OR_DIE(exc_ene_c,(exc_size), ierr) 311 312 do_ep_renorm = .FALSE. 313 ntemp = 1 314 do_ep_lifetime = .FALSE. 315 if(BSp%do_ep_renorm .and. present(Epren)) then 316 do_ep_renorm = .TRUE. 317 ntemp = Epren%ntemp 318 if(present(elph_lifetime)) then 319 do_ep_lifetime = elph_lifetime 320 end if 321 end if 322 323 if (do_ep_renorm) then 324 ABI_CHECK(nsppol == 1, "Nsppol == 2 not supported with elphon renormalizations") 325 end if 326 327 SELECT CASE (use_scalapack) 328 CASE (.FALSE.) 329 330 write(msg,'(a)')". Using LAPACK sequential version. " 331 call wrtout([std_out, ab_out], msg) 332 write(msg,'(a,f8.1,a)')' Allocating excitonic eigenvalues. Memory required: ', exc_size*dp*b2Mb,' Mb. ' 333 call wrtout(std_out, msg) 334 335 write(msg,'(a,f8.1,a)')' Allocating excitonic hamiltonian. Memory required: ',exc_size**2*dpc*b2Mb,' Mb.' 336 call wrtout(std_out, msg, do_flush=.True.) 337 338 ABI_MALLOC_OR_DIE(exc_mat,(exc_size,exc_size), ierr) 339 340 if (do_ep_renorm) then 341 ABI_MALLOC_OR_DIE(exc_vl,(exc_size,exc_size),ierr) 342 end if 343 !exc_mat = HUGE(zero) 344 ! 345 ! Read data from file. 346 if (open_file(hreso_fname,msg,newunit=hreso_unt,form="unformatted",status="old",action="read") /= 0) then 347 ABI_ERROR(msg) 348 end if 349 ! 350 ! Read the header and perform consistency checks. 351 call exc_read_bshdr(hreso_unt,Bsp,fform,ierr) 352 ABI_CHECK(ierr==0,"Fatal error, cannot continue") 353 ! 354 ! Construct full resonant block using Hermiticity. 355 diagonal_is_real = .not.Bsp%have_complex_ene 356 call exc_read_rblock_fio(hreso_unt,diagonal_is_real,nsppol,Bsp%nreh,exc_size,exc_mat,ierr) 357 ABI_CHECK(ierr==0,"Fatal error, cannot continue") 358 359 close(hreso_unt) 360 361 if (do_ep_renorm) then 362 write(std_out,'(a)') "Mapping kpts from bse to eph" 363 sppoldbl = 1 !; if (any(Cryst%symafm == -1) .and. Epren%nsppol == 1) nsppoldbl=2 364 ABI_MALLOC(bs2eph, (Kmesh%nbz*sppoldbl, 6)) 365 timrev = 1 366 call listkk(dksqmax, Cryst%gmet, bs2eph, Epren%kpts, Kmesh%bz, Epren%nkpt, Kmesh%nbz, Cryst%nsym, & 367 sppoldbl, Cryst%symafm, Cryst%symrel, timrev, xmpi_comm_self, use_symrec=.False.) 368 end if 369 370 do itemp = 1, ntemp 371 372 !TODO should find a way not to read again and again ! 373 ! but without storing it twice !!! 374 !exc_mat = HUGE(zero) 375 ! 376 ! Read data from file. 377 if (open_file(hreso_fname,msg,newunit=hreso_unt,form="unformatted",status="old",action="read") /= 0) then 378 ABI_ERROR(msg) 379 end if 380 ! 381 ! Read the header and perform consistency checks. 382 call exc_read_bshdr(hreso_unt,Bsp,fform,ierr) 383 ABI_CHECK(ierr==0,"Fatal error, cannot continue") 384 ! 385 ! Construct full resonant block using Hermiticity. 386 diagonal_is_real = .not.Bsp%have_complex_ene 387 call exc_read_rblock_fio(hreso_unt,diagonal_is_real,nsppol,Bsp%nreh,exc_size,exc_mat,ierr) 388 ABI_CHECK(ierr==0,"Fatal error, cannot continue") 389 390 close(hreso_unt) 391 392 bseig_fname = BS_files%out_eig 393 394 if (do_ep_renorm) then 395 write(std_out,'(a,i4)') "Will perform elphon renormalization for itemp = ",itemp 396 397 call int2char4(itemp,ts) 398 399 bseig_fname = TRIM(BS_files%out_eig) // TRIM("_T") // ts 400 401 ! Should patch the diagonal of exc_mat 402 403 do isppol = 1, BSp%nsppol 404 do ireh = 1, BSp%nreh(isppol) 405 ic = BSp%Trans(ireh,isppol)%c 406 iv = BSp%Trans(ireh,isppol)%v 407 ik = BSp%Trans(ireh,isppol)%k ! In the full bz 408 en = BSp%Trans(ireh,isppol)%en 409 410 ep_ik = bs2eph(ik,1) 411 412 !TODO support multiple spins ! 413 if(ABS(en - (Epren%eigens(ic,ep_ik,isppol)-Epren%eigens(iv,ep_ik,isppol)+BSp%mbpt_sciss)) > tol3) then 414 ABI_ERROR("Eigen from the transition does not correspond to the EP file !") 415 end if 416 exc_mat(ireh,ireh) = exc_mat(ireh,ireh) + (Epren%renorms(1,ic,ik,isppol,itemp) - Epren%renorms(1,iv,ik,isppol,itemp)) 417 418 ! Add lifetime 419 if(do_ep_lifetime) then 420 exc_mat(ireh,ireh) = exc_mat(ireh,ireh) - j_dpc*(Epren%linewidth(1,ic,ik,isppol,itemp) & 421 & + Epren%linewidth(1,iv,ik,isppol,itemp)) 422 end if 423 424 end do 425 end do 426 427 end if 428 429 if (do_full_diago) then 430 if(do_ep_renorm) then 431 call wrtout(std_out," Full diagonalization with XGEEV... ") 432 ABI_MALLOC(exc_vec,(exc_size,exc_size)) 433 call xgeev('V','V',exc_size,exc_mat,exc_size,exc_ene_c,exc_vl,exc_size,exc_vec,exc_size) 434 exc_mat(:,1:nstates) = exc_vec 435 ABI_FREE(exc_vec) 436 else 437 call wrtout(std_out," Full diagonalization with XHEEV... ") 438 call xheev("Vectors","Upper",exc_size,exc_mat,exc_ene) 439 exc_ene_c(:) = exc_ene(:) 440 end if 441 else 442 call wrtout(std_out," Partial diagonalization with XHEEVX... ") 443 abstol=zero; il=1; iu=nstates 444 ABI_MALLOC_OR_DIE(exc_vec,(exc_size,nstates),ierr) 445 call xheevx("Vectors","Index","Upper",exc_size,exc_mat,vl,vu,il,iu,abstol,mene_found,exc_ene,exc_vec,exc_size) 446 exc_mat(:,1:nstates) = exc_vec 447 exc_ene_c(:) = exc_ene(:) 448 ABI_FREE(exc_vec) 449 end if 450 ! 451 ! ============================================== 452 ! === Now exc_mat contains the eigenvectors ==== 453 ! ============================================== 454 455 ! * Write the final results. 456 call wrtout(std_out,' Writing eigenvalues and eigenvectors to file: '//TRIM(bseig_fname)) 457 458 if (open_file(bseig_fname,msg,newunit=eig_unt,form="unformatted",action="write") /= 0) then 459 ABI_ERROR(msg) 460 end if 461 462 !!! !DBYG 463 !!! !Compute overlap matrix 464 !!! ABI_MALLOC(ovlp,(exc_size,exc_size)) 465 !!! do mi=1,nstates 466 !!! do ireh=1,nstates 467 !!! ovlp(mi,ireh) = xdotc(exc_size,exc_vl(:,mi),1,exc_mat(:,ireh),1) 468 !!! if(mi==ireh) then 469 !!! !if(ABS(ovlp(mi,ireh)) < 0.999) then 470 !!! ! write(*,*) "it,itp = ",mi,ireh,"ovlp = ",ovlp(mi,ireh) 471 !!! !end if 472 !!! else 473 !!! if(ABS(ovlp(mi,ireh)) > 0.001) then 474 !!! write(*,*) "it,itp = ",mi,ireh,"ovlp = ",ovlp(mi,ireh) 475 !!! end if 476 !!! end if 477 !!! end do 478 !!! end do 479 !!! !call xgemm("C","N",exc_size,nstates,nstates,cone,exc_vl,exc_size,exc_mat,exc_size,czero,ovlp,nstates) 480 481 !!! write(777,*) ovlp 482 !!! ABI_FREE(ovlp) 483 !!! !ENDDBYG 484 485 !% fform = 1002 ! FIXME 486 !% call hdr_io_int(fform,Hdr_bse,2,eig_unt) 487 488 write(eig_unt) do_ep_lifetime 489 write(eig_unt) exc_size, nstates 490 write(eig_unt) exc_ene_c(1:nstates) 491 do mi=1,nstates 492 write(eig_unt) exc_mat(1:exc_size,mi) 493 if(do_ep_lifetime) then 494 write(eig_unt) exc_vl(1:exc_size,mi) 495 end if 496 end do 497 498 close(eig_unt) 499 500 end do ! itemp 501 502 ABI_FREE(exc_mat) 503 if (do_ep_renorm) then 504 ABI_FREE(exc_vl) 505 ABI_FREE(bs2eph) 506 end if 507 508 CASE (.TRUE.) 509 510 #if defined HAVE_LINALG_SCALAPACK && defined HAVE_MPI_IO 511 if (nsppol==2) then 512 ABI_WARNING("nsppol==2 + scalapack not coded yet") 513 end if 514 515 istwf_k=1; tbloc=50 516 write(msg,'(2(a,i0))')". Using scaLAPACK version with nprocs= ",nprocs,"; block size= ",tbloc 517 call wrtout([std_out, ab_out], msg) 518 519 write(msg,'(a,f8.1,a)')' Allocating excitonic eigenvalues. Memory required: ',exc_size*dp*b2Mb,' Mb. ' 520 call wrtout(std_out, msg) 521 ! 522 ! Init scaLAPACK environment. 523 call Slk_processor%init(comm) 524 ! 525 ! Init scaLAPACK matrices 526 call Slk_mat%init(exc_size,exc_size,Slk_processor,istwf_k) 527 call Slk_vec%init(exc_size,exc_size,Slk_processor,istwf_k) 528 ! 529 ! Open the file with MPI-IO and skip the record. 530 amode=MPI_MODE_RDONLY 531 532 call MPI_FILE_OPEN(comm, hreso_fname, amode, MPI_INFO_NULL, mpi_fh, ierr) 533 ABI_CHECK_MPI(ierr,"MPI_IO error opening file: "//TRIM(hreso_fname)) 534 535 ! Skip the header and find the offset for reading the matrix. 536 call exc_skip_bshdr_mpio(mpi_fh,xmpio_collective,ehdr_offset) 537 ! 538 ! Read scaLAPACK matrix from the file. 539 if (nsppol==1) then 540 call slk_read(Slk_mat,"Upper","Hermitian",is_fortran_file,mpi_fh=mpi_fh,offset=ehdr_offset) 541 else 542 array_of_sizes = (/exc_size,exc_size/) 543 block_sizes(:,1) = (/neh1,neh1/) 544 block_sizes(:,2) = (/neh2,neh2/) 545 block_sizes(:,3) = (/neh1,neh2/) 546 ABI_ERROR("Not tested") 547 !call slk_read_from_blocks(Slk_mat,array_of_sizes,block_sizes,is_fortran_file,mpi_fh=mpi_fh,offset=ehdr_offset) 548 end if 549 550 call MPI_FILE_CLOSE(mpi_fh, ierr) 551 ABI_CHECK_MPI(ierr,"FILE_CLOSE") 552 553 if (do_full_diago) then 554 call wrtout(std_out," Performing full diagonalization with scaLAPACK...") 555 call slk_mat%heev("Vectors","Upper",Slk_vec,exc_ene) 556 else 557 call wrtout(std_out," Performing partial diagonalization with scaLAPACK...") 558 il=1; iu=nstates; abstol=zero !ABSTOL = PDLAMCH(comm,'U') 559 call slk_mat%pzheevx("Vectors","Index","Upper",vl,vu,il,iu,abstol,Slk_vec,mene_found,exc_ene) 560 end if 561 562 exc_ene_c(:) = exc_ene(:) 563 564 call Slk_mat%free() 565 566 call wrtout(std_out,' Writing eigenvalues/vectors to file: '//TRIM(bseig_fname), do_flush=.True.) 567 568 ! Write distributed matrix on file bseig_fname with Fortran records. 569 if (my_rank==master) then ! Write exc eigenvalues. Vectors will be appended in slk_write. 570 if (open_file(bseig_fname,msg,newunit=eig_unt,form="unformatted",action="write") /= 0) then 571 ABI_ERROR(msg) 572 end if 573 write(eig_unt) exc_size, nstates 574 write(eig_unt) exc_ene_c(1:nstates) 575 close(eig_unt) 576 end if 577 578 call xmpi_barrier(comm) 579 ! 580 ! Open the file with MPI-IO and skip the record. 581 amode=MPI_MODE_RDWR 582 583 call MPI_FILE_OPEN(comm, bseig_fname, amode, MPI_INFO_NULL, mpi_fh, ierr) 584 ABI_CHECK_MPI(ierr,"MPI_IO error opening file: "//TRIM(hreso_fname)) 585 586 !call MPI_FILE_SYNC(mpi_fh,ierr) 587 588 ehdr_offset = 0 589 call xmpio_read_frm(mpi_fh,ehdr_offset,xmpio_collective,fmarker,ierr) 590 call xmpio_read_frm(mpi_fh,ehdr_offset,xmpio_collective,fmarker,ierr) 591 592 write(std_out,*)"Writing nstates ",nstates 593 gsub(:,1) = (/1,1/) 594 gsub(:,2) = (/exc_size,nstates/) 595 call slk_write(Slk_vec,"All",is_fortran_file,mpi_fh=mpi_fh,offset=ehdr_offset,glob_subarray=gsub) 596 597 call MPI_FILE_CLOSE(mpi_fh, ierr) 598 ABI_CHECK_MPI(ierr,"FILE_CLOSE") 599 600 call Slk_vec%free() 601 call Slk_processor%free() 602 call xmpi_barrier(comm) 603 #else 604 ABI_BUG("You should not be here!") 605 #endif 606 607 END SELECT 608 609 ! Order the eigenvalues 610 do ii=nstates,2,-1 611 do j=1,ii-1 612 if (DBLE(exc_ene_c(j)) > DBLE(exc_ene_c(j+1))) then 613 ctemp = exc_ene_c(j) 614 exc_ene_c(j) = exc_ene_c(j+1) 615 exc_ene_c(j+1) = ctemp 616 end if 617 end do 618 end do 619 620 621 write(msg,'(a,i4)')' Excitonic eigenvalues in eV up to n= ',nene_printed 622 call wrtout([std_out, ab_out], msg) 623 624 do it=0,(nene_printed-1)/8 625 write(msg,'(8f10.5)') ( DBLE(exc_ene_c(ii))*Ha_eV, ii=1+it*8,MIN(it*8+8,nene_printed) ) 626 call wrtout([std_out, ab_out], msg) 627 end do 628 629 exc_gap = MINVAL(DBLE(exc_ene_c(1:nstates))) 630 exc_maxene = MAXVAL(DBLE(exc_ene_c(1:nstates))) 631 632 write(msg,'(a,2(a,f7.2,2a),a)')ch10,& 633 " First excitonic eigenvalue= ",exc_gap*Ha_eV, " [eV]",ch10,& 634 " Last excitonic eigenvalue= ",exc_maxene*Ha_eV," [eV]",ch10,ch10 635 call wrtout([std_out, ab_out], msg, do_flush=.True.) 636 637 ABI_FREE(exc_ene_c) 638 ABI_FREE(exc_ene) 639 640 10 call xmpi_barrier(comm) 641 642 DBG_EXIT("PERS") 643 644 end subroutine exc_diago_resonant
m_exc_diago/exc_print_eig [ Functions ]
[ Top ] [ m_exc_diago ] [ Functions ]
NAME
exc_print_eig
FUNCTION
Print excitonic eigenvalues on std_out and ab_out.
INPUTS
gw_gap=GW direct gap. bseig_fname=The name of file containing eigenvalues and eigenvectors
OUTPUT
exc_gap=Excitonic direct gap. Additional info on the Excitonic spectrum are reported on standard output.
SOURCE
666 subroutine exc_print_eig(BSp,bseig_fname,gw_gap,exc_gap) 667 668 !Arguments ------------------------------------ 669 !scalars 670 complex(dpc),intent(in) :: gw_gap 671 complex(dpc),intent(out) :: exc_gap 672 character(len=*),intent(in) :: bseig_fname 673 type(excparam),intent(in) :: BSp 674 675 !Local variables ------------------------------ 676 !scalars 677 integer :: nstates_read,ii,j,k,eig_unt,ieig,hsize_exp 678 integer :: hsize_read !,nstates 679 complex(dpc) :: bind_energy,ctemp 680 character(len=500) :: msg 681 !type(Hdr_type) :: tmp_Hdr 682 !arrays 683 integer,allocatable :: iperm(:) 684 real(dp),allocatable :: exc_rene(:) 685 complex(dpc),allocatable :: exc_cene(:) 686 687 !************************************************************************ 688 689 exc_gap = czero 690 691 if (open_file(bseig_fname,msg,newunit=eig_unt,form="unformatted",status="old",action="read") /= 0) then 692 ABI_ERROR(msg) 693 end if 694 695 read(eig_unt) ! do_ep_lifetime 696 read(eig_unt) hsize_read, nstates_read 697 698 if (BSp%use_coupling==0) hsize_exp = SUM(Bsp%nreh) 699 if (BSp%use_coupling>0) hsize_exp = 2*SUM(Bsp%nreh) 700 701 if (hsize_exp /= hsize_read) then 702 write(msg,'(2(a,i0))')" Wrong dimension: read: ",hsize_read," expected= ",hsize_exp 703 ABI_ERROR(msg) 704 end if 705 706 ABI_MALLOC(exc_cene,(nstates_read)) 707 read(eig_unt) exc_cene(:) 708 709 ABI_MALLOC(exc_rene,(nstates_read)) 710 exc_rene = DBLE(exc_cene) 711 712 ABI_MALLOC(iperm,(nstates_read)) 713 iperm = (/(ii, ii=1,nstates_read)/) 714 715 call sort_dp(nstates_read,exc_rene,iperm,tol6) 716 717 ABI_FREE(exc_rene) 718 ABI_FREE(iperm) 719 720 ! put in ascending order 721 do ii=nstates_read,2,-1 722 do j=1,ii-1 723 if (DBLE(exc_cene(j)) > DBLE(exc_cene(j+1))) then 724 ctemp = exc_cene(j) 725 exc_cene(j) = exc_cene(j+1) 726 exc_cene(j+1) = ctemp 727 end if 728 end do 729 end do 730 731 exc_gap = DCMPLX(ABS(DBLE(exc_cene(1))),AIMAG(exc_cene(1))) 732 733 do ii=1,nstates_read 734 if (ABS(DBLE(exc_cene(ii))) < DBLE(exc_gap)) then 735 exc_gap = DCMPLX(ABS(DBLE(exc_cene(ii))),AIMAG(exc_cene(ii))) 736 end if 737 end do 738 739 bind_energy = gw_gap - exc_gap 740 741 write(msg,"(3(a,2f6.2,2a))")& 742 & " GW direct gap ",gw_gap*Ha_eV, " [eV] ",ch10,& 743 & " EXC direct gap ",exc_gap*Ha_eV, " [eV] ",ch10,& 744 & " EXC binding energy ",bind_energy*Ha_eV," [eV] ",ch10 745 call wrtout([std_out, ab_out], msg) 746 747 msg=' Excitonic eigenvalues up to the GW energy gap [eV]' 748 call wrtout([std_out, ab_out], msg) 749 750 do ii=1,nstates_read 751 if (DBLE(exc_cene(ii)) > zero) EXIT 752 end do 753 754 do j=ii,nstates_read 755 if (DBLE(exc_cene(j)) > DBLE(gw_gap)) EXIT 756 end do 757 j=j-1 758 759 do ieig=ii,j 760 write(msg,'(i3,a,2f6.2,a)')ieig," (",exc_cene(ieig)*Ha_eV,")" 761 call wrtout([std_out, ab_out], msg) 762 end do 763 764 ii=ii-1 765 do j=ii,1,-1 766 if (ABS(DBLE(exc_cene(j))) > DBLE(gw_gap)) EXIT 767 end do 768 j=j+1 769 770 ! This coding is not portable, write to ab_out has been disabled. 771 if (ii>0) then 772 do k=ii,j,-1 773 write(msg,'(i3,a,2f6.2,a)')k," (",exc_cene(k)*Ha_eV,")" 774 call wrtout(std_out, msg) 775 end do 776 end if 777 778 ABI_FREE(exc_cene) 779 780 close(eig_unt) 781 782 end subroutine exc_print_eig