TABLE OF CONTENTS


ABINIT/m_exc_diago [ Modules ]

[ Top ] [ Modules ]

NAME

 m_exc_diago

FUNCTION

COPYRIGHT

 Copyright (C) 2009-2018 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 .

PARENTS

SOURCE

18 #if defined HAVE_CONFIG_H
19 #include "config.h"
20 #endif
21 
22 #include "abi_common.h"
23 
24 MODULE m_exc_diago
25 
26  use defs_basis
27  use defs_datatypes
28  use m_slk
29  use m_bs_defs
30  use m_abicore
31  use m_errors
32  use m_xmpi
33 #if defined HAVE_MPI2
34  use mpi
35 #endif
36  use m_hdr
37  use m_sort
38 
39  use defs_abitypes,     only : hdr_type
40  use m_io_tools,        only : open_file
41  use m_fstrings,        only : int2char4
42  use m_numeric_tools,   only : print_arr, hermitianize
43  use m_crystal,         only : crystal_t
44  use m_kpts,            only : listkk
45  use m_bz_mesh,         only : kmesh_t
46  use m_ebands,          only : ebands_report_gap
47  use m_eprenorms,       only : eprenorms_t
48  use m_wfd,             only : wfd_t
49  use m_paw_hr,          only : pawhur_t
50  use m_pawtab,          only : pawtab_type
51  use m_exc_itdiago,     only : exc_iterative_diago
52  use m_hide_lapack,     only : xheev, xheevx, xgeev, xhegvx, xginv, xhdp_invert, xhegv
53  use m_hide_blas,       only : xdotc, xgemm
54  use m_bse_io,          only : exc_fullh_from_blocks, offset_in_file, rrs_of_glob, ccs_of_glob, &
55 &                              exc_read_bshdr, exc_skip_bshdr_mpio, exc_read_rblock_fio
56  use m_exc_spectra,     only : build_spectra
57 
58  implicit none
59 
60  private
61 
62 #if defined HAVE_MPI1
63  include 'mpif.h'
64 #endif
65 
66 !#define DEV_MG_DEBUG_THIS
67 
68  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.

PARENTS

      m_exc_diago

CHILDREN

      destruction_matrix_scalapack,end_scalapack,exc_fullh_from_blocks
      exc_read_bshdr,exc_skip_bshdr_mpio,hermitianize,idx_glob
      init_matrix_scalapack,init_scalapack,mpi_file_close,mpi_file_open
      mpi_file_read_all,mpi_file_set_view,mpi_type_free,slk_pzgemm
      slk_pzhegvx,slk_single_fview_read_mask,slk_write,slk_zinvert,wrtout
      xgemm,xhdp_invert,xhegv,xhegvx,xmpi_barrier,xmpio_read_frm

SOURCE

 900 subroutine exc_diago_coupling(Bsp,BS_files,Hdr_bse,prtvol,comm)
 901 
 902 
 903 !This section has been created automatically by the script Abilint (TD).
 904 !Do not modify the following lines by hand.
 905 #undef ABI_FUNC
 906 #define ABI_FUNC 'exc_diago_coupling'
 907 !End of the abilint section
 908 
 909  implicit none
 910 
 911 !Arguments ------------------------------------
 912 !scalars
 913  integer,intent(in) :: comm,prtvol
 914  type(excfiles),intent(in) :: BS_files
 915  type(excparam),intent(in) :: BSp
 916  type(Hdr_type),intent(in) :: Hdr_bse
 917 
 918 !Local variables ------------------------------
 919 !scalars
 920  integer,parameter :: master=0,ldvl=1
 921  integer(i8b) :: bsize_ham
 922  integer :: ii,exc_size,hreso_unt,hcoup_unt,eig_unt,nsppol,nstates
 923  integer :: bsz,block,bs1,bs2,jj
 924  integer :: fform,row_sign
 925  integer :: mi,it,nprocs,my_rank !itp
 926  integer :: nene_printed,ierr
 927  real(dp) :: exc_gap,exc_maxene,temp
 928  logical :: diago_is_real,do_full_diago
 929  logical :: do_ep_lifetime
 930  character(len=500) :: msg
 931  character(len=fnlen) :: hreso_fname,hcoup_fname,bseig_fname
 932 !arrays
 933  complex(dpc),allocatable :: exc_ham(:,:),exc_rvect(:,:),exc_ene(:),ovlp(:,:)
 934  complex(dpc),allocatable :: cbuff(:,:)
 935  complex(dpc) :: vl_dpc(ldvl,1)
 936 
 937 !************************************************************************
 938 
 939  my_rank = xmpi_comm_rank(comm); nprocs  = xmpi_comm_size(comm)
 940 
 941  nsppol = Hdr_bse%nsppol
 942  if (nsppol==2) then
 943    MSG_WARNING("nsppol==2 with coupling is still under development")
 944  end if
 945 
 946  if (nprocs > 1) then
 947    MSG_WARNING("Scalapack does not provide ZGEEV, diagonalization is done in sequential!")
 948  end if
 949 
 950  exc_size = 2*SUM(BSp%nreh)
 951  nstates  = BSp%nstates
 952  do_full_diago = (exc_size==nstates)
 953  ABI_CHECK(do_full_diago,"Partial diago not coded yet")
 954 
 955  bseig_fname = BS_files%out_eig
 956  if (BS_files%in_eig /= BSE_NOFILE) then
 957    MSG_ERROR("BS_files%in_eig is defined!")
 958  end if
 959  !
 960  ! Only master performs the diagonalization since ScaLAPACK does not provide the parallel version of ZGEEV.
 961  if (my_rank/=master) GOTO 10
 962 
 963  write(msg,'(a,i0)')' Direct diagonalization of the full excitonic Hamiltonian, Matrix size= ',exc_size
 964  call wrtout(std_out,msg,"COLL")
 965  call wrtout(ab_out,msg,"COLL")
 966 
 967  bsize_ham = 2*dpc*exc_size**2
 968  write(msg,'(a,f9.2,a)')' Allocating full excitonic Hamiltonian. Memory requested: ',bsize_ham*b2Gb,' Gb. '
 969  call wrtout(std_out,msg,"COLL")
 970 
 971  ABI_STAT_MALLOC(exc_ham,(exc_size,exc_size), ierr)
 972  ABI_CHECK(ierr==0, 'out of memory: full excitonic hamiltonian')
 973 
 974  write(msg,'(3a,f8.1,3a,f8.1,a)')&
 975 &  ' Allocating excitonic eigenvalues and eigenvectors. ',ch10,&
 976 &  ' Memory-space requested: ',2*dpc*exc_size*b2Gb,' Gb. ',ch10,&
 977 &  ' Memory-space requested: ',bsize_ham*b2Gb,' Gb. '
 978  call wrtout(std_out,msg,"COLL")
 979 
 980  ABI_STAT_MALLOC(exc_ene,(exc_size), ierr)
 981  ABI_CHECK(ierr==0, 'out of memory: exc_ene')
 982 
 983  if (BS_files%in_hreso /= BSE_NOFILE) then
 984    hreso_fname = BS_files%in_hreso
 985  else
 986    hreso_fname = BS_files%out_hreso
 987  end if
 988 
 989  call wrtout(std_out,' Reading resonant excitonic Hamiltonian from '//TRIM(hreso_fname),"COLL")
 990 
 991  if (open_file(hreso_fname,msg,newunit=hreso_unt,form="unformatted",status="old",action="read") /= 0) then
 992    MSG_ERROR(msg)
 993  end if
 994  !
 995  ! Read the header and perform consistency checks.
 996  call exc_read_bshdr(hreso_unt,Bsp,fform,ierr)
 997  ABI_CHECK(ierr==0,"Wrong header")
 998  !
 999  ! Construct resonant and anti-resonant part of the excitonic Hamiltonian using Hermiticity. File is always in double precision.
1000  ! Fill exc_ham with ( R  0 )
1001  !                   ( 0 -R*)
1002 !BEGINDEBUG
1003  exc_ham = HUGE(one)
1004 !ENDDEBUG
1005 
1006  row_sign=-1; diago_is_real=(.not.BSp%have_complex_ene)
1007  call exc_fullh_from_blocks(hreso_unt,"Resonant",nsppol,row_sign,diago_is_real,BSp%nreh,exc_size,exc_ham)
1008  close(hreso_unt)
1009 
1010  if (BS_files%in_hcoup /= BSE_NOFILE) then
1011    hcoup_fname =  BS_files%in_hcoup
1012  else
1013    hcoup_fname =  BS_files%out_hcoup
1014  end if
1015 
1016  call wrtout(std_out,' Reading coupling excitonic Hamiltonian from '//TRIM(hcoup_fname),"COLL")
1017  if (open_file(hcoup_fname,msg,newunit=hcoup_unt,form="unformatted",status="old",action="read") /= 0) then
1018    MSG_ERROR(msg)
1019  end if
1020  !
1021  ! Read the header and perform consistency checks.
1022  call exc_read_bshdr(hcoup_unt,Bsp,fform,ierr)
1023  ABI_CHECK(ierr==0,"Wrong header")
1024  !
1025  ! Fill exc_ham with ( 0  C) to have ( R   C )
1026  !                   (-C* 0)         (-C* -R*)
1027  row_sign=-1; diago_is_real=(.not.BSp%have_complex_ene) ! not used here
1028  call exc_fullh_from_blocks(hcoup_unt,"Coupling",nsppol,row_sign,diago_is_real,BSp%nreh,exc_size,exc_ham)
1029 
1030 !BEGINDEBUG
1031  if (ANY(exc_ham==HUGE(one))) then
1032    write(msg,'(a,2(1x,i0))')"There is a bug in exc_fullh_from_blocks",COUNT(exc_ham==HUGE(one)),exc_size**2
1033    MSG_WARNING(msg)
1034    bsz = Bsp%nreh(1)
1035    ABI_MALLOC(cbuff,(bsz,bsz))
1036    block=0
1037    do jj=1,2*nsppol
1038      do ii=1,2*nsppol
1039        block=block+1
1040        bs1 = (ii-1)*bsz+1
1041        bs2 = (jj-1)*bsz+1
1042        cbuff = exc_ham(bs1:bs1+bsz-1,bs2:bs2+bsz-1)
1043        if (ANY(cbuff==HUGE(one))) then
1044          write(std_out,*)" for block ",ii,jj," found ",COUNT(cbuff==HUGE(one))," wrong entries"
1045        end if
1046      end do
1047    end do
1048 
1049    ABI_FREE(cbuff)
1050    MSG_ERROR("Cannot continue")
1051  end if
1052 !ENDDEBUG
1053 
1054  close(hcoup_unt)
1055  !
1056  ! ======================================================
1057  ! ==== Calculate right eigenvectors and eigenvalues ====
1058  ! ======================================================
1059  ABI_STAT_MALLOC(exc_rvect,(exc_size,exc_size), ierr)
1060  ABI_CHECK(ierr==0, "out of memory: excitonic eigenvectors")
1061 
1062  if (do_full_diago) then
1063    call wrtout(std_out,"Complete direct diagonalization with xgeev...","COLL")
1064    call xgeev("No_left_eigen","Vectors",exc_size,exc_ham,exc_size,exc_ene,vl_dpc,ldvl,exc_rvect,exc_size)
1065  else
1066    MSG_ERROR("Not implemented error")
1067  end if
1068 
1069  ABI_FREE(exc_ham)
1070 
1071  exc_gap    = MINVAL(ABS(DBLE (exc_ene(1:nstates))))
1072  exc_maxene = MAXVAL(ABS(DBLE (exc_ene(1:nstates))))
1073  temp       = MAXVAL(ABS(AIMAG(exc_ene(1:nstates))))
1074 
1075  write(msg,'(2(a,f7.2,2a),a,es9.2,2a)')&
1076 &  " First excitonic eigenvalue: ",exc_gap*Ha_eV,   " [eV].",ch10,&
1077 &  " Last  excitonic eigenvalue: ",exc_maxene*Ha_eV," [eV].",ch10,&
1078 &  " Largest imaginary part:     ",temp*Ha_eV,      " [eV] ",ch10
1079  call wrtout(std_out,msg,"COLL")
1080  call wrtout(ab_out,msg,"COLL")
1081 
1082  nene_printed = MIN(32*nsppol,nstates); if (prtvol>10) nene_printed = nstates
1083 
1084  ! This is not portable as the the eigenvalues calculated by ZGEEV are not sorted.
1085  ! Even two subsequent calculations with the same input on the same machine
1086  ! might produce different orderings. Might sort the eigenvalues though, just for printing.
1087 
1088  write(msg,'(a,i0)')' Complex excitonic eigenvalues in eV up to n= ',nene_printed
1089  call wrtout(std_out,msg,"PERS")
1090 
1091  do it=0,(nene_printed-1)/4
1092    write(msg,'(8f10.5)') ( exc_ene(ii)*Ha_eV, ii=1+it*4,MIN(it*4+4,nene_printed) )
1093    call wrtout(std_out,msg,"PERS")
1094  end do
1095 
1096  call wrtout(std_out,ch10//" Writing eigenvalues and eigenvectors on file "//TRIM(bseig_fname),"COLL")
1097 
1098  if (open_file(bseig_fname,msg,newunit=eig_unt,form="unformatted",action="write") /= 0) then
1099    MSG_ERROR(msg)
1100  end if
1101 
1102 !YG : new version with lifetime
1103  do_ep_lifetime = .FALSE.
1104  write(eig_unt) do_ep_lifetime
1105 
1106  write(eig_unt)exc_size,nstates
1107  write(eig_unt)CMPLX(exc_ene(1:nstates),kind=dpc)
1108  do mi=1,nstates
1109    write(eig_unt) exc_rvect(:,mi)
1110  end do
1111 
1112  ABI_FREE(exc_ene)
1113 
1114  ABI_STAT_MALLOC(ovlp,(nstates,nstates), ierr)
1115  ABI_CHECK(ierr==0, 'out of memory in ovlp matrix')
1116 
1117  call wrtout(std_out,' Calculating overlap matrix... ',"COLL")
1118 
1119  !do itp=1,nstates
1120  !  do it=1,nstates
1121  !    ovlp(it,itp) = xdotc(exc_size,exc_rvect(:,it),1,exc_rvect(:,itp),1)
1122  !  end do
1123  !end do
1124  call xgemm("C","N",exc_size,nstates,nstates,cone,exc_rvect,exc_size,exc_rvect,exc_size,czero,ovlp,nstates)
1125  ABI_FREE(exc_rvect)
1126 
1127  call wrtout(std_out," Inverting overlap matrix... ","COLL")
1128 
1129  ! Version for generic complex matrix.
1130  !call xginv(ovlp,exc_size)
1131 
1132  ! The overlap is Hermitian definite positive.
1133  call xhdp_invert("Upper",ovlp,nstates)
1134  call hermitianize(ovlp,"Upper")
1135 
1136  call wrtout(std_out,' Writing overlap matrix S^-1 on file: '//TRIM(bseig_fname),"COLL")
1137 
1138  do it=1,nstates
1139    write(eig_unt) CMPLX(ovlp(:,it),kind=dpc)
1140  end do
1141 
1142  ABI_FREE(ovlp)
1143 
1144  close(eig_unt)
1145 
1146 10 call xmpi_barrier(comm)
1147 
1148 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 on file BS_files%out_eig.

PARENTS

      m_exc_diago

CHILDREN

      destruction_matrix_scalapack,end_scalapack,exc_fullh_from_blocks
      exc_read_bshdr,exc_skip_bshdr_mpio,hermitianize,idx_glob
      init_matrix_scalapack,init_scalapack,mpi_file_close,mpi_file_open
      mpi_file_read_all,mpi_file_set_view,mpi_type_free,slk_pzgemm
      slk_pzhegvx,slk_single_fview_read_mask,slk_write,slk_zinvert,wrtout
      xgemm,xhdp_invert,xhegv,xhegvx,xmpi_barrier,xmpio_read_frm

SOURCE

1184 subroutine exc_diago_coupling_hegv(Bsp,BS_files,Hdr_bse,prtvol,comm)
1185 
1186 
1187 !This section has been created automatically by the script Abilint (TD).
1188 !Do not modify the following lines by hand.
1189 #undef ABI_FUNC
1190 #define ABI_FUNC 'exc_diago_coupling_hegv'
1191 !End of the abilint section
1192 
1193  implicit none
1194 
1195 !Arguments ------------------------------------
1196 !scalars
1197  integer,intent(in) :: comm,prtvol
1198  type(excparam),intent(in) :: BSp
1199  type(excfiles),intent(in) :: BS_files
1200  type(Hdr_type),intent(in) :: Hdr_bse
1201 
1202 !Local variables ------------------------------
1203 !scalars
1204  integer,parameter :: master=0
1205  integer(i8b) :: bsize_ham
1206  integer :: itype,il,iu,spin,row1,row2,pad_r1,pad_r2,neh1,neh2
1207  integer :: ii,exc_size,hreso_unt,hcoup_unt,eig_unt
1208  integer :: fform,neig_found,nstates
1209  integer :: mi,it,nprocs,my_rank
1210  integer :: nene_printed,nsppol,row_sign,ierr
1211  real(dp) :: exc_gap,exc_maxene,abstol,vl,vu
1212  character(len=500) :: msg
1213  character(len=fnlen) :: reso_fname,coup_fname,bseig_fname
1214  logical :: use_scalapack,do_full_diago,diago_is_real
1215  logical :: do_ep_lifetime
1216 !arrays
1217  real(dp),allocatable :: exc_ene(:) !,test_ene(:)
1218  complex(dpc),allocatable :: exc_ham(:,:),exc_rvect(:,:),fmat(:,:),ovlp(:,:)
1219 #if defined HAVE_LINALG_SCALAPACK && defined HAVE_MPI_IO
1220  integer,parameter :: istwfk1=1
1221  integer :: amode,mpi_fh,tbloc,tmp_unt,mene_found,mpi_err,my_nel,nsblocks
1222  integer :: iloc,jj,jloc,iglob,jglob,etype,slk_mask_type,offset_err,el,rrs_kind,ccs_kind
1223  integer :: max_r,max_c
1224  integer(XMPI_OFFSET_KIND) :: ehdr_offset,fmarker,my_offset
1225  integer :: gsub(2,2)
1226  logical,parameter :: is_fortran_file=.TRUE.
1227  complex(dpc) :: ctmp
1228  integer,allocatable :: sub_block(:,:,:)
1229  integer,pointer :: myel2loc(:,:)
1230  complex(dpc),allocatable :: tmp_cbuffer(:)
1231  character(50) :: uplo
1232  real(dp),external :: PDLAMCH
1233  type(matrix_scalapack)    :: Slk_F,Slk_Hbar,Slk_vec,Slk_ovlp,Slk_tmp
1234  type(processor_scalapack) :: Slk_processor
1235 #endif
1236 
1237 !************************************************************************
1238 
1239  my_rank = xmpi_comm_rank(comm); nprocs  = xmpi_comm_size(comm)
1240 
1241  nsppol = Hdr_bse%nsppol
1242  if (nsppol==2) then
1243    MSG_WARNING("nsppol==2 is still under development!")
1244  end if
1245 
1246  neh1 = BSp%nreh(1); neh2=neh1
1247  if (nsppol==2) neh2 = BSp%nreh(2)
1248 
1249  exc_size = 2*SUM(Bsp%nreh)
1250  nstates  = Bsp%nstates
1251  do_full_diago=(nstates==exc_size)
1252 
1253  write(msg,'(a,i0)')'. Direct diagonalization of the full excitonic Hamiltonian, Matrix size= ',exc_size
1254  call wrtout(std_out,msg,"COLL")
1255  call wrtout(ab_out,msg,"COLL")
1256 
1257  bseig_fname = BS_files%out_eig
1258  if (BS_files%in_eig /= BSE_NOFILE) then
1259    MSG_ERROR("BS_files%in_eig is defined!")
1260  end if
1261 
1262  if (BS_files%in_hreso /= BSE_NOFILE) then
1263    reso_fname = BS_files%in_hreso
1264  else
1265    reso_fname = BS_files%out_hreso
1266  end if
1267  call wrtout(std_out,' Reading resonant excitonic Hamiltonian from '//TRIM(reso_fname),"COLL")
1268 
1269  if (BS_files%in_hcoup /= BSE_NOFILE) then
1270    coup_fname =  BS_files%in_hcoup
1271  else
1272    coup_fname =  BS_files%out_hcoup
1273  end if
1274  call wrtout(std_out,' Reading coupling excitonic Hamiltonian from '//TRIM(coup_fname),"COLL")
1275 
1276  use_scalapack = .FALSE.
1277 #ifdef HAVE_LINALG_SCALAPACK
1278  use_scalapack = (nprocs > 1)
1279 #endif
1280  !use_scalapack = .FALSE.
1281  !use_scalapack = .TRUE.
1282 
1283  if (.not.use_scalapack .and. my_rank/=master) GOTO 10
1284 
1285  ABI_STAT_MALLOC(exc_ene,(exc_size), ierr)
1286  ABI_CHECK(ierr==0, 'out of memory: exc_ene')
1287 
1288  SELECT CASE (use_scalapack)
1289 
1290  CASE (.FALSE.)
1291    write(msg,'(a)')". Using LAPACK sequential version to solve FHv = ev with H positive definite. "
1292    call wrtout(std_out,msg,"PERS")
1293    call wrtout(ab_out,msg,"COLL")
1294 
1295    bsize_ham = 2*dpc*exc_size**2
1296    write(msg,'(a,f9.2,a)')' Allocating full excitonic Hamiltonian. Memory requested: ',2*bsize_ham*b2Gb,' Gb. '
1297    call wrtout(std_out,msg,"COLL")
1298 
1299    ABI_STAT_MALLOC(exc_ham,(exc_size,exc_size), ierr)
1300    ABI_CHECK(ierr==0, 'out of memory: full excitonic hamiltonian')
1301 
1302    ABI_STAT_MALLOC(fmat,(exc_size,exc_size), ierr)
1303    ABI_CHECK(ierr==0, 'out of memory: fmat')
1304 
1305    write(msg,'(3a,f8.1,3a,f8.1,a)')&
1306 &    ' Allocating excitonic eigenvalues and eigenvectors. ',ch10,&
1307 &    ' Memory-space requested: ',2*dpc*exc_size*b2Gb,' Gb. ',ch10,&
1308 &    ' Memory-space requested: ',bsize_ham*b2Gb,' Gb. '
1309    call wrtout(std_out,msg,"COLL")
1310 
1311    if (open_file(reso_fname,msg,newunit=hreso_unt,form="unformatted",status="old",action="read") /= 0) then
1312      MSG_ERROR(msg)
1313    end if
1314    !
1315    ! Read the header and perform consistency checks.
1316    call exc_read_bshdr(hreso_unt,Bsp,fform,ierr)
1317    ABI_CHECK(ierr==0,"Wrong header")
1318    !
1319    ! Construct Hbar = ( R   C )
1320    !                  ( C*  R*)
1321    !
1322    row_sign=+1; diago_is_real=(.not.BSp%have_complex_ene)
1323    call exc_fullh_from_blocks(hreso_unt,"Resonant",nsppol,row_sign,diago_is_real,Bsp%nreh,exc_size,exc_ham)
1324    close(hreso_unt)
1325 
1326    if (open_file(coup_fname,msg,newunit=hcoup_unt,form="unformatted",status="old",action="read") /= 0) then
1327      MSG_ERROR(msg)
1328    end if
1329    !
1330    ! Read the header and perform consistency checks.
1331    call exc_read_bshdr(hcoup_unt,Bsp,fform,ierr)
1332    ABI_CHECK(ierr==0,"Wrong header")
1333 
1334    row_sign=+1; diago_is_real=(.not.BSp%have_complex_ene) ! not used here.
1335    call exc_fullh_from_blocks(hcoup_unt,"Coupling",nsppol,row_sign,diago_is_real,Bsp%nreh,exc_size,exc_ham)
1336    close(hcoup_unt)
1337 
1338 #ifdef DEV_MG_DEBUG_THIS
1339 write(666)exc_ham
1340 #endif
1341    !
1342    ! Fill fmat = (1  0)
1343    !             (0 -1)
1344    fmat = czero
1345    do spin=1,nsppol
1346      pad_r1 = (spin-1)*Bsp%nreh(1)
1347      pad_r2 = SUM(Bsp%nreh)
1348      if (spin==2) pad_r2 = pad_r2 + Bsp%nreh(1)
1349      do it=1,Bsp%nreh(spin)
1350        row1 = it + pad_r1
1351        row2 = it + pad_r2
1352        fmat(row1,row1) =  cone
1353        fmat(row2,row2) = -cone
1354      end do
1355    end do
1356    !
1357    ! ==================================================
1358    ! ==== Solve generalized EV problem F H u = e u ====
1359    ! ==================================================
1360    ! 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.
1361    !
1362    itype=2
1363    if (do_full_diago) then
1364      call wrtout(std_out," Full diagonalization with XHEGV... ","COLL")
1365      call xhegv(itype,"Vectors","Upper",exc_size,fmat,exc_ham,exc_ene)
1366    else
1367      call wrtout(std_out," Partial diagonalization with XHEGVX... ","COLL")
1368      ABI_STAT_MALLOC(exc_rvect,(exc_size,nstates), ierr)
1369      ABI_CHECK(ierr==0, "out of memory: excitonic eigenvectors")
1370      il=1; iu=1; abstol=zero
1371      call xhegvx(itype,"Vectors","All","Upper",exc_size,fmat,exc_ham,vl,vu,il,iu,abstol,neig_found,exc_ene,exc_rvect,exc_size)
1372    end if
1373 
1374    ABI_FREE(exc_ham)
1375 
1376    if (do_full_diago) then
1377      ABI_MALLOC(exc_rvect,(exc_size,nstates))
1378      exc_rvect = fmat(:,1:nstates)
1379    end if
1380 
1381    ABI_FREE(fmat)
1382 
1383    call wrtout(std_out," Writing eigenvalues and eigenvectors on file: "//TRIM(bseig_fname),"COLL")
1384 
1385    if (open_file(bseig_fname,msg,newunit=eig_unt,form="unformatted",action="write") /= 0) then
1386      MSG_ERROR(msg)
1387    end if
1388 
1389    do_ep_lifetime = .FALSE.
1390    write(eig_unt) do_ep_lifetime
1391    write(eig_unt) exc_size, nstates
1392    write(eig_unt) CMPLX(exc_ene(1:nstates),kind=dpc)
1393    do mi=1,nstates
1394      write(eig_unt) CMPLX(exc_rvect(:,mi),kind=dpc)
1395    end do
1396 
1397 #ifdef DEV_MG_DEBUG_THIS
1398    write(888)exc_rvect
1399    write(888)exc_ene
1400 #endif
1401 
1402    ABI_STAT_MALLOC(ovlp,(nstates,nstates), ierr)
1403    ABI_CHECK(ierr==0, 'out of memory in ovlp matrix')
1404 
1405    call wrtout(std_out,' Calculating overlap matrix...',"COLL")
1406 
1407    call xgemm("C","N",exc_size,nstates,nstates,cone,exc_rvect,exc_size,exc_rvect,exc_size,czero,ovlp,nstates)
1408    ABI_FREE(exc_rvect)
1409 
1410 #ifdef DEV_MG_DEBUG_THIS
1411 write(667)ovlp
1412 #endif
1413 
1414    call wrtout(std_out," Inverting overlap matrix... ","COLL")
1415    !
1416    ! The overlap is Hermitian definite positive.
1417    call xhdp_invert("Upper",ovlp,nstates)
1418    call hermitianize(ovlp,"Upper")
1419 
1420    ! Version for generic complex matrix.
1421    !call xginv(ovlp,nstates)
1422 
1423 #ifdef DEV_MG_DEBUG_THIS
1424 write(668,*)ovlp
1425 #endif
1426 
1427    call wrtout(std_out,' Writing overlap matrix O^-1 on file: '//TRIM(bseig_fname),"COLL")
1428    do it=1,nstates
1429      write(eig_unt) ovlp(:,it)
1430    end do
1431 
1432    ABI_FREE(ovlp)
1433    close(eig_unt)
1434 
1435  CASE (.TRUE.)
1436 
1437 #if defined HAVE_LINALG_SCALAPACK && defined HAVE_MPI_IO
1438    !
1439    ! Init scaLAPACK matrix Hbar = ( R   C )
1440    !                              ( C*  R*)
1441    ! Battle plan:
1442    !   Here the reading is complicated by the fact that R and C are stored on two different files
1443    !   and moreover the matrices are in packed storage mode.
1444    !   For initializing the local part of the resonant and anti-resonant block we have to allocate
1445    !   a temporary buffer. Then we read the buffer from file and the corresponding elements of the
1446    !   scaLAPACK matrix are initialized taking into account the symmetries of R (Hermitian)
1447    !   The same procedure is used to read the coupling and the anti-coupling part (Symmetric).
1448    !
1449    tbloc=50
1450    write(msg,'(2(a,i0))')". Using MPI-IO + scaLAPACK version with nprocs= ",nprocs,"; block size= ",tbloc
1451    call wrtout(std_out,msg,"PERS",do_flush=.True.)
1452    call wrtout(ab_out,msg,"COLL",do_flush=.True.)
1453    !
1454    ! Init scaLAPACK environment.
1455    call init_scalapack(Slk_processor,comm)
1456    !
1457    ! Open the Resonant file with MPI-IO and skip the record.
1458    amode=MPI_MODE_RDONLY
1459 
1460    call MPI_FILE_OPEN(comm, reso_fname, amode, MPI_INFO_NULL, mpi_fh, mpi_err)
1461    msg = " MPI_IO error opening file: "//TRIM(reso_fname)
1462    ABI_CHECK_MPI(mpi_err,msg)
1463    !
1464    ! Skip the header and find the offset for reading the matrix.
1465    call exc_skip_bshdr_mpio(mpi_fh,xmpio_collective,ehdr_offset)
1466    !
1467    ! Read  = ( R  - )
1468    !         ( -  R*)
1469    call init_matrix_scalapack(Slk_Hbar,exc_size,exc_size,Slk_processor,istwfk1,tbloc=tbloc)
1470 
1471    nullify(myel2loc)
1472    nsblocks=nsppol
1473    ABI_MALLOC(sub_block,(2,2,nsblocks))
1474    ABI_CHECK(nsppol==1,"nsppol==2 not coded yet")
1475 
1476    call slk_single_fview_read_mask(Slk_Hbar,rrs_of_glob,offset_in_file,nsblocks,sub_block,my_nel,myel2loc,etype,slk_mask_type,&
1477 &    offset_err,is_fortran_file)
1478 
1479    if (offset_err/=0) then
1480      write(msg,"(3a)")&
1481 &      " Global position index cannot be stored in a standard Fortran integer ",ch10,&
1482 &      " Excitonic matrix cannot be read with a single MPI-IO call."
1483      MSG_ERROR(msg)
1484    end if
1485 
1486    ! Shift the offset because the view starts at the fist matrix element!
1487    ! TODO should rationalize the treatment of the offset
1488    my_offset = ehdr_offset + xmpio_bsize_frm
1489    call MPI_FILE_SET_VIEW(mpi_fh, my_offset, etype, slk_mask_type, 'native', MPI_INFO_NULL, mpi_err)
1490    ABI_CHECK_MPI(mpi_err,"SET_VIEW")
1491 
1492    call MPI_TYPE_FREE(slk_mask_type,mpi_err)
1493    ABI_CHECK_MPI(mpi_err,"MPI_TYPE_FREE")
1494    !
1495    ! Read my portion of the R,-R* sublocks and store the values in a temporary buffer.
1496    ABI_STAT_MALLOC(tmp_cbuffer,(my_nel), ierr)
1497    ABI_CHECK(ierr==0, " out of memory tmp_cbuffer")
1498 
1499    call xmpi_barrier(comm)
1500 
1501    call MPI_FILE_READ_ALL(mpi_fh, tmp_cbuffer, my_nel, MPI_DOUBLE_COMPLEX, MPI_STATUS_IGNORE, mpi_err)
1502    ABI_CHECK_MPI(mpi_err,"READ_ALL")
1503    !
1504    ! Symmetrize my Resonant part.
1505    do el=1,my_nel
1506      iloc = myel2loc(1,el)
1507      jloc = myel2loc(2,el)
1508      call idx_glob(Slk_Hbar,iloc,jloc,iglob,jglob)
1509      ctmp = tmp_cbuffer(el)
1510      if (iglob==jglob.and..not.Bsp%have_complex_ene) ctmp = DBLE(ctmp) ! Force the diagonal to be real.
1511      rrs_kind = rrs_of_glob(iglob,jglob,Slk_Hbar%sizeb_global)
1512      if (rrs_kind==1.and.jglob<iglob) then ! Lower resonant
1513        ctmp = DCONJG(ctmp)
1514      else if (rrs_kind==-1.and.jglob>=iglob) then  ! Lower Anti-resonant (Diagonal is included).
1515        ctmp = DCONJG(ctmp)
1516      end if
1517      Slk_Hbar%buffer_cplx(iloc,jloc) = ctmp
1518    end do
1519 
1520    ABI_FREE(tmp_cbuffer)
1521    ABI_FREE(myel2loc)
1522 
1523    call MPI_FILE_CLOSE(mpi_fh, mpi_err)
1524    ABI_CHECK_MPI(mpi_err,"FILE_CLOSE")
1525    !
1526    ! Read  = ( -  C)
1527    !         (-C* -)
1528    !
1529    call MPI_FILE_OPEN(comm, coup_fname, amode, MPI_INFO_NULL, mpi_fh, mpi_err)
1530    msg = " MPI_IO error opening file: "//TRIM(coup_fname)
1531    ABI_CHECK_MPI(mpi_err,msg)
1532    !
1533    ! Skip the header and find the offset for reading the matrix.
1534    call exc_skip_bshdr_mpio(mpi_fh,xmpio_collective,ehdr_offset)
1535 
1536    nullify(myel2loc)
1537    call slk_single_fview_read_mask(Slk_Hbar,ccs_of_glob,offset_in_file,nsblocks,sub_block,my_nel,myel2loc,etype,slk_mask_type,&
1538 &    offset_err,is_fortran_file)
1539 
1540    ABI_FREE(sub_block)
1541 
1542    if (offset_err/=0) then
1543      write(msg,"(3a)")&
1544 &      " Global position index cannot be stored in a standard Fortran integer ",ch10,&
1545 &      " Excitonic matrix cannot be read with a single MPI-IO call."
1546      MSG_ERROR(msg)
1547    end if
1548    !
1549    ! Shift the offset because the view starts at the fist matrix element!
1550    ! TODO should rationalize the treatment of the offset so that the client code
1551    ! will automatically receive my_offset.
1552    my_offset = ehdr_offset + xmpio_bsize_frm
1553    call MPI_FILE_SET_VIEW(mpi_fh, my_offset, etype, slk_mask_type, 'native', MPI_INFO_NULL, mpi_err)
1554    ABI_CHECK_MPI(mpi_err,"SET_VIEW")
1555 
1556    call MPI_TYPE_FREE(slk_mask_type,mpi_err)
1557    ABI_CHECK_MPI(mpi_err,"MPI_TYPE_FREE")
1558    !
1559    ! Read my portion of the C-C* blocks and store the values in a temporary buffer.
1560    ABI_STAT_MALLOC(tmp_cbuffer,(my_nel), ierr)
1561    ABI_CHECK(ierr==0, " out of memory tmp_cbuffer")
1562 
1563    call MPI_FILE_READ_ALL(mpi_fh, tmp_cbuffer, my_nel, MPI_DOUBLE_COMPLEX, MPI_STATUS_IGNORE, mpi_err)
1564    ABI_CHECK_MPI(mpi_err,"READ_ALL")
1565    !
1566    ! Symmetrize my coupling part.
1567    ! Coupling block is symmetric => No symmetrization of the lower triangle.
1568    do el=1,my_nel
1569      iloc = myel2loc(1,el)
1570      jloc = myel2loc(2,el)
1571      call idx_glob(Slk_Hbar,iloc,jloc,iglob,jglob)
1572      ccs_kind = ccs_of_glob(iglob,jglob,Slk_Hbar%sizeb_global)
1573      ctmp = tmp_cbuffer(el)
1574      if (ccs_kind==-1) ctmp = DCONJG(ctmp) ! Anti-coupling (Diagonal is included).
1575      Slk_Hbar%buffer_cplx(iloc,jloc) = ctmp
1576    end do
1577 
1578    ABI_FREE(tmp_cbuffer)
1579    ABI_FREE(myel2loc)
1580 
1581    !max_r=20; max_c=10
1582    !call print_arr(Slk_Hbar%buffer_cplx,max_r=max_r,max_c=max_c,unit=std_out)
1583 
1584 #ifdef DEV_MG_DEBUG_THIS
1585    ABI_MALLOC(exc_ham,(exc_size,exc_size))
1586    read(666)exc_ham
1587 
1588    write(std_out,*)"Error Hbar: ",MAXVAL(ABS(exc_ham-Slk_Hbar%buffer_cplx))
1589    ABI_FREE(exc_ham)
1590 #endif
1591 
1592    call MPI_FILE_CLOSE(mpi_fh, mpi_err)
1593    ABI_CHECK_MPI(mpi_err,"FILE_CLOSE")
1594    !
1595    ! Init scaLAPACK matrix F
1596    call init_matrix_scalapack(Slk_F,exc_size,exc_size,Slk_processor,istwfk1,tbloc=tbloc)
1597    !
1598    ! Global F = (1  0)
1599    !            (0 -1)
1600    do jloc=1,Slk_F%sizeb_local(2)
1601      do iloc=1,Slk_F%sizeb_local(1)
1602        call idx_glob(Slk_F,iloc,jloc,iglob,jglob)
1603        if (iglob==jglob) then
1604          if (iglob<=SUM(Bsp%nreh)) then
1605            Slk_F%buffer_cplx(iloc,jloc) =  cone
1606          else
1607            Slk_F%buffer_cplx(iloc,jloc) = -cone
1608          end if
1609        else
1610          Slk_F%buffer_cplx(iloc,jloc) =  czero
1611        end if
1612      end do
1613    end do
1614    !
1615    ! ===========================================================
1616    ! ==== Solve generalized EV problem H u = F Hbar u = e u ====
1617    ! ===========================================================
1618    call init_matrix_scalapack(Slk_vec,exc_size,exc_size,Slk_processor,istwfk1,tbloc=tbloc)
1619    !
1620    itype=2; vl=1; vu=1; il=1; iu=nstates
1621    abstol=zero !ABSTOL = PDLAMCH(comm,'U')
1622 
1623 #if 1
1624    if (do_full_diago) then
1625      call slk_pzhegvx(itype,"Vectors","All","Upper",Slk_F,Slk_Hbar,vl,vu,il,iu,abstol,Slk_vec,mene_found,exc_ene)
1626    else
1627      MSG_WARNING("Partial diago is still under testing")
1628      call slk_pzhegvx(itype,"Vectors","Index","Upper",Slk_F,Slk_Hbar,vl,vu,il,iu,abstol,Slk_vec,mene_found,exc_ene)
1629    end if
1630 #else
1631    call xhegv(itype,"Vectors","Upper",exc_size,Slk_F%buffer_cplx,Slk_Hbar%buffer_cplx,exc_ene)
1632    Slk_vec%buffer_cplx = Slk_F%buffer_cplx
1633 #endif
1634 
1635 #ifdef DEV_MG_DEBUG_THIS
1636    if (PRODUCT(Slk_Hbar%sizeb_local) /= exc_size**2) then
1637      MSG_ERROR("Wrong size")
1638    end if
1639 
1640    ABI_MALLOC(exc_ham,(exc_size,exc_size))
1641    read(888)exc_ham
1642 
1643    write(std_out,*)"Error rvec: ",MAXVAL(ABS(exc_ham-Slk_vec%buffer_cplx))
1644    ABI_FREE(exc_ham)
1645 
1646    ABI_MALLOC(test_ene,(exc_size))
1647    read(888)test_ene
1648    write(std_out,*)"Error ene: ",MAXVAL(ABS(exc_ene-test_ene))
1649    ABI_FREE(test_ene)
1650 #endif
1651 
1652    call destruction_matrix_scalapack(Slk_F)
1653    call destruction_matrix_scalapack(Slk_Hbar)
1654 
1655    call wrtout(std_out,ch10//" Writing eigenvalues and eigenvectors on file: "//TRIM(bseig_fname),"COLL")
1656    !
1657    ! Open the file with Fortran-IO to write the Header.
1658    if (my_rank==master) then
1659      if (open_file(bseig_fname,msg,newunit=eig_unt,form="unformatted",action="write") /= 0) then
1660        MSG_ERROR(msg)
1661      end if
1662 
1663      write(eig_unt) exc_size,nstates
1664      write(eig_unt) CMPLX(exc_ene(1:nstates),kind=dpc)
1665      !do mi=1,exc_size
1666      !  write(eig_unt) CMPLX(exc_rvect(:,mi),kind=dpc)
1667      !end do
1668      close(eig_unt)
1669    end if
1670    !
1671    ! Open the file with MPI-IO and write the distributed eigevectors.
1672    call xmpi_barrier(comm)
1673    amode=MPI_MODE_RDWR
1674    call MPI_FILE_OPEN(comm, bseig_fname, amode, MPI_INFO_NULL, mpi_fh, mpi_err)
1675    ABI_CHECK_MPI(mpi_err,"FILE_OPEN: "//TRIM(bseig_fname))
1676    !
1677    ! Skip the header and find the offset for writing the matrix.
1678    ehdr_offset=0
1679    !call hdr_mpio_skip(mpi_fh,fform,ehdr_offset)
1680 
1681    call xmpio_read_frm(mpi_fh,ehdr_offset,xmpio_collective,fmarker,mpi_err)
1682    write(std_out,*)" fmarker1 = ",fmarker
1683    call xmpio_read_frm(mpi_fh,ehdr_offset,xmpio_collective,fmarker,mpi_err)
1684    write(std_out,*)" fmarker2 = ",fmarker
1685 
1686    write(std_out,*)" Writing nstates ",nstates
1687    gsub(:,1) = (/1,1/)
1688    gsub(:,2) = (/exc_size,nstates/)
1689    call slk_write(Slk_vec,"All",is_fortran_file,mpi_fh=mpi_fh,offset=ehdr_offset,glob_subarray=gsub)
1690 
1691    call wrtout(std_out,' Calculating overlap matrix... ',"COLL")
1692    if (.not.do_full_diago) then
1693      MSG_ERROR(" Init of Slk_ovlp is wrong")
1694    end if
1695 
1696    call init_matrix_scalapack(Slk_ovlp,exc_size,exc_size,Slk_processor,istwfk1,tbloc=tbloc)
1697 
1698    ! Calculate the overlap matrix.
1699    ! FIXME
1700    ! The ESLL manual says that "matrices matrix1 and matrix2 must have no common elements; otherwise, results are unpredictable."
1701    ! However the official scaLAPACK documentation does not report this (severe) limitation.
1702 
1703    !call init_matrix_scalapack(Slk_tmp,exc_size,exc_size,Slk_processor,istwfk1,tbloc=tbloc)
1704    !Slk_tmp%buffer_cplx = Slk_vec%buffer_cplx
1705    !call slk_pzgemm("C","N",Slk_tmp,cone,Slk_vec,czero,Slk_ovlp)
1706    !call destruction_matrix_scalapack(Slk_tmp)
1707 
1708    call slk_pzgemm("C","N",Slk_vec,cone,Slk_vec,czero,Slk_ovlp)
1709 
1710 #ifdef DEV_MG_DEBUG_THIS
1711    ABI_MALLOC(exc_ham,(exc_size,exc_size))
1712    read(667)exc_ham
1713 
1714    write(std_out,*)"Error Ovlp: ",MAXVAL(ABS(exc_ham-Slk_ovlp%buffer_cplx))
1715    !Slk_ovlp%buffer_cplx = exc_ham
1716 #endif
1717 
1718    !max_r=20; max_c=10
1719    !call print_arr(Slk_ovlp%buffer_cplx,max_r=max_r,max_c=max_c,unit=std_out)
1720 
1721    call destruction_matrix_scalapack(Slk_vec)
1722 
1723    call wrtout(std_out," Inverting overlap matrix... ","COLL")
1724    uplo="Upper"
1725 
1726 #if 0
1727 !DEBUG
1728    call xhdp_invert(uplo,Slk_ovlp%buffer_cplx,exc_size)
1729 
1730    !call slk_symmetrize(Slk_ovlp,uplo,"Hermitian")
1731    call hermitianize(Slk_ovlp%buffer_cplx,uplo)
1732 
1733    exc_ham = MATMUL(exc_ham,Slk_ovlp%buffer_cplx)
1734    do it=1,exc_size
1735      exc_ham(it,it) = exc_ham(it,it) - cone
1736    end do
1737 
1738    write(std_out,*)"Error Inversion: ",MAXVAL(ABS(exc_ham))
1739    ABI_FREE(exc_ham)
1740 !END DEBUG
1741 
1742 #else
1743    ! call slk_zdhp_invert(Slk_ovlp,uplo)
1744    ! call hermitianize(Slk_ovlp%buffer_cplx,uplo)
1745    ! !call slk_symmetrize(Slk_ovlp,uplo,"Hermitian")
1746 
1747    call slk_zinvert(Slk_ovlp)  ! Version for generic complex matrix.
1748 #endif
1749 
1750    if (allocated(exc_ham))  then
1751      ABI_FREE(exc_ham)
1752    end if
1753 
1754 #ifdef DEV_MG_DEBUG_THIS
1755    ABI_MALLOC(exc_ham,(exc_size,exc_size))
1756    read(668)exc_ham
1757    write(std_out,*)"Error in Inv Ovlp: ",MAXVAL(ABS(exc_ham-Slk_ovlp%buffer_cplx))
1758 
1759    !exc_ham = exc_ham-Slk_ovlp%buffer_cplx
1760    !do it=1,exc_size
1761    !  if ( MAXVAL(ABS(exc_ham(:,it))) > 0.1 ) write(std_out,*)"it: ",it,exc_ham(:,it)
1762    !end do
1763 
1764    !Slk_ovlp%buffer_cplx = exc_ham
1765    ABI_FREE(exc_ham)
1766 
1767    !write(std_out,*)"MAX ERR",MAXVAL(ABS(Slk_ovlp%buffer_cplx - TRANSPOSE(DCONJG(Slk_ovlp%buffer_cplx))))
1768 #endif
1769 
1770    call wrtout(std_out,' Writing overlap matrix S^-1 on file: '//TRIM(bseig_fname),"COLL")
1771 
1772    call slk_write(Slk_ovlp,"All",is_fortran_file,mpi_fh=mpi_fh,offset=ehdr_offset)
1773 
1774    call MPI_FILE_CLOSE(mpi_fh, mpi_err)
1775    ABI_CHECK_MPI(mpi_err,"FILE_CLOSE")
1776 
1777    call destruction_matrix_scalapack(Slk_ovlp)
1778 
1779    call end_scalapack(Slk_processor)
1780 #else
1781    MSG_BUG("You should not be here!")
1782 #endif
1783 
1784  END SELECT
1785 
1786  exc_gap    = MINVAL(ABS(exc_ene(1:nstates)))
1787  exc_maxene = MAXVAL(ABS(exc_ene(1:nstates)))
1788 
1789  write(msg,'(2(a,f7.2,2a))')&
1790 &  " First excitonic eigenvalue: ",exc_gap*Ha_eV,   " [eV].",ch10,&
1791 &  " Last  excitonic eigenvalue: ",exc_maxene*Ha_eV," [eV].",ch10
1792  call wrtout(std_out,msg,"COLL")
1793  call wrtout(ab_out,msg,"COLL")
1794 
1795  nene_printed = MIN(32*nsppol,nstates); if (prtvol>10) nene_printed = nstates
1796  write(msg,'(a,i0)')' Complex excitonic eigenvalues in eV up to n= ',nene_printed
1797  call wrtout(std_out,msg,"PERS")
1798 
1799  do it=0,(nene_printed-1)/4
1800    write(msg,'(4f10.5)') ( exc_ene(ii)*Ha_eV, ii=1+it*4,MIN(it*4+4,nene_printed) )
1801    call wrtout(std_out,msg,"COLL")
1802  end do
1803 
1804  ABI_FREE(exc_ene)
1805 
1806 10 call xmpi_barrier(comm)
1807 
1808 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.

PARENTS

      bethe_salpeter

CHILDREN

      destruction_matrix_scalapack,end_scalapack,exc_fullh_from_blocks
      exc_read_bshdr,exc_skip_bshdr_mpio,hermitianize,idx_glob
      init_matrix_scalapack,init_scalapack,mpi_file_close,mpi_file_open
      mpi_file_read_all,mpi_file_set_view,mpi_type_free,slk_pzgemm
      slk_pzhegvx,slk_single_fview_read_mask,slk_write,slk_zinvert,wrtout
      xgemm,xhdp_invert,xhegv,xhegvx,xmpi_barrier,xmpio_read_frm

SOURCE

101 subroutine exc_diago_driver(Wfd,Bsp,BS_files,KS_BSt,QP_BSt,Cryst,Kmesh,Psps,&
102 &  Pawtab,Hur,Hdr_bse,drude_plsmf,Epren)
103 
104 
105 !This section has been created automatically by the script Abilint (TD).
106 !Do not modify the following lines by hand.
107 #undef ABI_FUNC
108 #define ABI_FUNC 'exc_diago_driver'
109 !End of the abilint section
110 
111  implicit none
112 
113 !Arguments ------------------------------------
114 !scalars
115  real(dp),intent(in) :: drude_plsmf
116  type(excparam),intent(in) :: BSp
117  type(excfiles),intent(in) ::  BS_files
118  type(Hdr_type),intent(in) :: Hdr_bse
119  type(crystal_t),intent(in) :: Cryst
120  type(pseudopotential_type),intent(in) :: Psps
121  type(kmesh_t),intent(in) :: Kmesh
122  type(ebands_t),intent(in) :: KS_BSt,QP_BSt
123  type(wfd_t),intent(inout) :: Wfd
124  type(eprenorms_t),intent(in) :: Epren
125 !arrays
126  type(pawtab_type),intent(in) :: Pawtab(Cryst%ntypat*Wfd%usepaw)
127  type(pawhur_t),intent(in) :: Hur(Cryst%natom*Wfd%usepaw)
128 
129 !Local variables ------------------------------
130 !scalars
131  integer :: my_rank,master,comm,prtvol
132  complex(dpc) :: exc_gap,gw_gap
133  logical :: eval_eigenstates
134  character(len=500) :: msg
135 !arrays
136  real(dp) :: gaps(3,QP_BSt%nsppol)
137 
138 !************************************************************************
139 
140  DBG_ENTER("COLL")
141 
142  comm    = Wfd%comm
143  my_rank = Wfd%my_rank
144  master  = Wfd%master
145  prtvol  = Wfd%prtvol
146 
147  if (BSp%have_complex_ene) then
148    MSG_ERROR("Complex energies are not supported yet")
149  end if
150  !
151  ! This trick is needed to restart a CG run, use DDIAGO to calculate the spectra reusing an old BSEIG file.
152  eval_eigenstates = (BS_files%in_eig == BSE_NOFILE) .or. (Bsp%algorithm == BSE_ALGO_CG)
153 
154  if (eval_eigenstates) then
155    !
156    select case (BSp%algorithm)
157    case (BSE_ALGO_DDIAGO)
158      if (BSp%use_coupling==0) then
159        call exc_diago_resonant(BSp,BS_files,Hdr_bse,prtvol,comm)
160        if(Bsp%do_ep_renorm) then
161          call exc_diago_resonant(BSp,BS_files,Hdr_bse,prtvol,comm,Epren=Epren,Kmesh=Kmesh,Cryst=Cryst,elph_lifetime=.TRUE.)
162        end if
163      else
164        if (Bsp%have_complex_ene) then
165          ! Solve Hv = ev with generic complex matrix.
166          call exc_diago_coupling(BSp,BS_files,Hdr_bse,prtvol,comm)
167        else
168          ! Solve generalized eigenvalue problem F Hbar with Hbar Hermitian definitive positive matrix.
169          call exc_diago_coupling_hegv(BSp,BS_files,Hdr_bse,prtvol,comm)
170        end if
171      end if
172 
173    case (BSE_ALGO_CG)
174      if (BSp%use_coupling==0) then
175        call exc_iterative_diago(Bsp,BS_files,Hdr_bse,prtvol,comm)
176      else
177        MSG_ERROR("CG + coupling not coded")
178      end if
179 
180    case default
181      write(msg,'(a,i0)')" Wrong value for Bsp%algorithm: ",Bsp%algorithm
182      MSG_ERROR(msg)
183    end select
184    !
185    if (my_rank==master) then
186      call ebands_report_gap(QP_BSt,header="QP bands",unit=std_out,gaps=gaps)
187      gw_gap = MINVAL(gaps(2,:))
188      call exc_print_eig(BSp,BS_files%out_eig,gw_gap,exc_gap)
189    end if
190    call xmpi_barrier(comm)
191    !
192  end if
193 
194  call build_spectra(BSp,BS_files,Cryst,Kmesh,KS_BSt,QP_BSt,Psps,Pawtab,Wfd,Hur,drude_plsmf,comm)
195 
196  ! Electron-phonon renormalization !
197  if (BSp%algorithm == BSE_ALGO_DDIAGO .and. BSp%use_coupling == 0 .and. BSp%do_ep_renorm) then
198    call build_spectra(BSp,BS_files,Cryst,Kmesh,KS_BSt,QP_BSt,Psps,Pawtab,Wfd,Hur,drude_plsmf,comm,Epren=Epren)
199  end if
200 
201 
202  DBG_EXIT("COLL")
203 
204 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

PARENTS

      m_exc_diago

CHILDREN

      destruction_matrix_scalapack,end_scalapack,exc_fullh_from_blocks
      exc_read_bshdr,exc_skip_bshdr_mpio,hermitianize,idx_glob
      init_matrix_scalapack,init_scalapack,mpi_file_close,mpi_file_open
      mpi_file_read_all,mpi_file_set_view,mpi_type_free,slk_pzgemm
      slk_pzhegvx,slk_single_fview_read_mask,slk_write,slk_zinvert,wrtout
      xgemm,xhdp_invert,xhegv,xhegvx,xmpi_barrier,xmpio_read_frm

SOURCE

239 subroutine exc_diago_resonant(Bsp,BS_files,Hdr_bse,prtvol,comm,Epren,Kmesh,Cryst,elph_lifetime)
240 
241 
242 !This section has been created automatically by the script Abilint (TD).
243 !Do not modify the following lines by hand.
244 #undef ABI_FUNC
245 #define ABI_FUNC 'exc_diago_resonant'
246 !End of the abilint section
247 
248  implicit none
249 
250 !Arguments ------------------------------------
251 !scalars
252  integer,intent(in) :: comm,prtvol
253  logical,optional,intent(in) :: elph_lifetime
254  type(excparam),intent(in) :: BSp
255  type(excfiles),intent(in) :: BS_files
256  type(Hdr_type),intent(in) :: Hdr_bse
257  type(eprenorms_t),optional,intent(in) :: Epren
258  type(kmesh_t),optional,intent(in) :: Kmesh
259  type(crystal_t),optional,intent(in) :: Cryst
260 
261 !Local variables ------------------------------
262 !scalars
263  integer,parameter :: master=0
264  integer :: ii,it,mi,hreso_unt,eig_unt,exc_size,neh1,neh2,j
265  integer :: nsppol,il,iu,mene_found,nstates
266  integer :: nprocs,my_rank,fform,nene_printed,ierr
267  real(dp) :: exc_gap,exc_maxene,abstol
268  real(dp) :: vl,vu
269  logical :: use_scalapack,do_full_diago,diagonal_is_real
270  character(len=500) :: msg
271  character(len=fnlen) :: hreso_fname,bseig_fname
272 !arrays
273  real(dp),allocatable :: exc_ene(:)
274  complex(dpc),allocatable :: exc_mat(:,:),exc_vec(:,:)
275 #if defined HAVE_LINALG_SCALAPACK && defined HAVE_MPI_IO
276  integer :: amode,mpi_fh,istwf_k,tbloc,tmp_unt
277  integer :: itloc,jj,jtloc,itglob,jtglob
278  integer(XMPI_OFFSET_KIND) :: ehdr_offset,fmarker
279  integer :: block_sizes(2,3),array_of_sizes(2),gsub(2,2)
280  logical,parameter :: is_fortran_file=.TRUE.
281  real(dp),external :: PDLAMCH
282  type(matrix_scalapack)    :: Slk_mat,Slk_vec
283  type(processor_scalapack) :: Slk_processor
284 #endif
285 
286  integer :: ik, ic, iv, isppol, ireh, ep_ik, itemp
287  complex(dpc) :: en
288 
289  real(dp) :: dksqmax
290  integer,allocatable :: bs2eph(:,:)
291  integer :: sppoldbl, timrev
292  logical :: do_ep_renorm, do_ep_lifetime
293  integer :: ntemp
294  character(len=4) :: ts
295  complex(dpc),allocatable :: exc_vl(:,:),exc_ene_c(:)
296  complex(dpc) :: ctemp
297 !! complex(dpc),allocatable :: ovlp(:,:)
298 
299 !************************************************************************
300 
301  DBG_ENTER("PERS")
302 
303  nprocs  = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
304 
305  if (BSp%have_complex_ene) then ! QP lifetimes are not included
306    MSG_ERROR("complex energies not coded yet")
307  end if
308 
309  if (ANY(Bsp%nreh/=Bsp%nreh(1))) then
310    write(std_out,*)" Bsp%nreh: ",Bsp%nreh
311    write(msg,'(a)')" BSE code does not support different number of transitions for the two spin channels"
312    MSG_WARNING(msg)
313  end if
314 
315  nsppol   = Hdr_bse%nsppol
316  exc_size = SUM(BSp%nreh)
317  nstates  = BSp%nstates; do_full_diago=(Bsp%nstates==exc_size)
318 
319  neh1 = Bsp%nreh(1); neh2 = neh1
320  if (Hdr_bse%nsppol==2) neh2 = Bsp%nreh(2)
321 
322  use_scalapack = .FALSE.
323 #if defined HAVE_LINALG_SCALAPACK
324  use_scalapack = (nprocs > 1)
325 #endif
326  !use_scalapack = .FALSE.
327  !use_scalapack = .TRUE.
328  if (use_scalapack .and. nsppol == 2) then
329    use_scalapack = .False.
330    msg = "Scalapack with nsppol==2 not yet available. Using sequential version"
331    MSG_WARNING(msg)
332  end if
333 
334  if (.not.use_scalapack .and. my_rank/=master) GOTO 10 ! Inversion is done by master only.
335 
336  nene_printed = MIN(32*nsppol,nstates); if (prtvol>10) nene_printed = nstates
337 
338  if (BS_files%in_hreso /= BSE_NOFILE) then
339    hreso_fname = BS_files%in_hreso
340  else
341    hreso_fname = BS_files%out_hreso
342  end if
343 
344  bseig_fname = BS_files%out_eig
345  if (BS_files%in_eig /= BSE_NOFILE) then
346    MSG_ERROR("BS_files%in_eig is defined!")
347  end if
348 
349  write(msg,'(a,i0)')' Direct diagonalization of the resonant excitonic Hamiltonian, Matrix size= ',exc_size
350  call wrtout(std_out,msg,"COLL")
351  call wrtout(ab_out,msg,"COLL")
352 
353  ABI_STAT_MALLOC(exc_ene,(exc_size), ierr)
354  ABI_CHECK(ierr==0, 'out of memory: excitonic eigenvalues')
355 
356  ABI_STAT_MALLOC(exc_ene_c,(exc_size), ierr)
357  ABI_CHECK(ierr==0,'out of memory: excitonic complex eigenvalues')
358 
359  do_ep_renorm = .FALSE.
360  ntemp = 1
361  do_ep_lifetime = .FALSE.
362  if(BSp%do_ep_renorm .and. present(Epren)) then
363    do_ep_renorm = .TRUE.
364    ntemp = Epren%ntemp
365    if(present(elph_lifetime)) then
366      do_ep_lifetime = elph_lifetime
367    end if
368  end if
369 
370  if (do_ep_renorm) then
371    ABI_CHECK(nsppol == 1, "Nsppol == 2 not supported with elphon renormalizations")
372  end if
373 
374  SELECT CASE (use_scalapack)
375  CASE (.FALSE.)
376 
377    write(msg,'(a)')". Using LAPACK sequential version. "
378    call wrtout(std_out,msg,"PERS")
379    call wrtout(ab_out,msg,"COLL")
380 
381    write(msg,'(a,f8.1,a)')' Allocating excitonic eigenvalues. Memory required: ', exc_size*dp*b2Mb,' Mb. '
382    call wrtout(std_out,msg,"COLL")
383 
384    write(msg,'(a,f8.1,a)')' Allocating excitonic hamiltonian.  Memory required: ',exc_size**2*dpc*b2Mb,' Mb.'
385    call wrtout(std_out,msg,"COLL",do_flush=.True.)
386 
387    ABI_STAT_MALLOC(exc_mat,(exc_size,exc_size), ierr)
388    ABI_CHECK(ierr==0, 'out of memory: excitonic hamiltonian')
389 
390    if (do_ep_renorm) then
391      ABI_STAT_MALLOC(exc_vl,(exc_size,exc_size),ierr)
392      ABI_CHECK(ierr==0,'out of memory for left eigenvectors !')
393    end if
394    !exc_mat = HUGE(zero)
395    !
396    ! Read data from file.
397    if (open_file(hreso_fname,msg,newunit=hreso_unt,form="unformatted",status="old",action="read") /= 0) then
398      MSG_ERROR(msg)
399    end if
400    !
401    ! Read the header and perform consistency checks.
402    call exc_read_bshdr(hreso_unt,Bsp,fform,ierr)
403    ABI_CHECK(ierr==0,"Fatal error, cannot continue")
404    !
405    ! Construct full resonant block using Hermiticity.
406    diagonal_is_real = .not.Bsp%have_complex_ene
407    call exc_read_rblock_fio(hreso_unt,diagonal_is_real,nsppol,Bsp%nreh,exc_size,exc_mat,ierr)
408    ABI_CHECK(ierr==0,"Fatal error, cannot continue")
409 
410    close(hreso_unt)
411 
412 
413    if (do_ep_renorm) then
414      write(std_out,'(a)') "Mapping kpts from bse to eph"
415      sppoldbl = 1 !; if (any(Cryst%symafm == -1) .and. Epren%nsppol == 1) nsppoldbl=2
416      ABI_MALLOC(bs2eph, (Kmesh%nbz*sppoldbl, 6))
417      timrev = 1
418      call listkk(dksqmax, Cryst%gmet, bs2eph, Epren%kpts, Kmesh%bz, Epren%nkpt, Kmesh%nbz, Cryst%nsym, &
419 &       sppoldbl, Cryst%symafm, Cryst%symrel, timrev, use_symrec=.False.)
420    end if
421 
422    do itemp = 1, ntemp
423 
424      !TODO should find a way not to read again and again !
425      !  but without storing it twice !!!
426      !exc_mat = HUGE(zero)
427      !
428      ! Read data from file.
429      if (open_file(hreso_fname,msg,newunit=hreso_unt,form="unformatted",status="old",action="read") /= 0) then
430        MSG_ERROR(msg)
431      end if
432      !
433      ! Read the header and perform consistency checks.
434      call exc_read_bshdr(hreso_unt,Bsp,fform,ierr)
435      ABI_CHECK(ierr==0,"Fatal error, cannot continue")
436      !
437      ! Construct full resonant block using Hermiticity.
438      diagonal_is_real = .not.Bsp%have_complex_ene
439      call exc_read_rblock_fio(hreso_unt,diagonal_is_real,nsppol,Bsp%nreh,exc_size,exc_mat,ierr)
440      ABI_CHECK(ierr==0,"Fatal error, cannot continue")
441 
442      close(hreso_unt)
443 
444      bseig_fname = BS_files%out_eig
445 
446      if (do_ep_renorm) then
447        write(std_out,'(a,i4)') "Will perform elphon renormalization for itemp = ",itemp
448 
449        call int2char4(itemp,ts)
450 
451        bseig_fname = TRIM(BS_files%out_eig) // TRIM("_T") // ts
452 
453        ! Should patch the diagonal of exc_mat
454 
455        do isppol = 1, BSp%nsppol
456          do ireh = 1, BSp%nreh(isppol)
457            ic = BSp%Trans(ireh,isppol)%c
458            iv = BSp%Trans(ireh,isppol)%v
459            ik = BSp%Trans(ireh,isppol)%k ! In the full bz
460            en = BSp%Trans(ireh,isppol)%en
461 
462            ep_ik = bs2eph(ik,1)
463 
464            !TODO support multiple spins !
465            if(ABS(en - (Epren%eigens(ic,ep_ik,isppol)-Epren%eigens(iv,ep_ik,isppol)+BSp%mbpt_sciss)) > tol3) then
466              MSG_ERROR("Eigen from the transition does not correspond to the EP file !")
467            end if
468            exc_mat(ireh,ireh) = exc_mat(ireh,ireh) + (Epren%renorms(1,ic,ik,isppol,itemp) - Epren%renorms(1,iv,ik,isppol,itemp))
469 
470            ! Add lifetime
471            if(do_ep_lifetime) then
472              exc_mat(ireh,ireh) = exc_mat(ireh,ireh) - j_dpc*(Epren%lifetimes(1,ic,ik,isppol,itemp) &
473 &                + Epren%lifetimes(1,iv,ik,isppol,itemp))
474            end if
475 
476          end do
477        end do
478 
479      end if
480 
481      if (do_full_diago) then
482        if(do_ep_renorm) then
483          call wrtout(std_out," Full diagonalization with XGEEV... ","COLL")
484          ABI_MALLOC(exc_vec,(exc_size,exc_size))
485          call xgeev('V','V',exc_size,exc_mat,exc_size,exc_ene_c,exc_vl,exc_size,exc_vec,exc_size)
486          exc_mat(:,1:nstates) = exc_vec
487          ABI_FREE(exc_vec)
488        else
489          call wrtout(std_out," Full diagonalization with XHEEV... ","COLL")
490          call xheev("Vectors","Upper",exc_size,exc_mat,exc_ene)
491          exc_ene_c(:) = exc_ene(:)
492        end if
493      else
494        call wrtout(std_out," Partial diagonalization with XHEEVX... ","COLL")
495        abstol=zero; il=1; iu=nstates
496        ABI_STAT_MALLOC(exc_vec,(exc_size,nstates),ierr)
497        ABI_CHECK(ierr==0,"out of memory in exc_vec")
498        call xheevx("Vectors","Index","Upper",exc_size,exc_mat,vl,vu,il,iu,abstol,mene_found,exc_ene,exc_vec,exc_size)
499        exc_mat(:,1:nstates) = exc_vec
500        exc_ene_c(:) = exc_ene(:)
501        ABI_FREE(exc_vec)
502      end if
503      !
504      ! ==============================================
505      ! === Now exc_mat contains the eigenvectors ====
506      ! ==============================================
507 
508      ! * Write the final results.
509      call wrtout(std_out,' Writing eigenvalues and eigenvectors on file: '//TRIM(bseig_fname),"COLL")
510 
511      if (open_file(bseig_fname,msg,newunit=eig_unt,form="unformatted",action="write") /= 0) then
512        MSG_ERROR(msg)
513      end if
514 
515      !!! !DBYG
516      !!! !Compute overlap matrix
517      !!! ABI_MALLOC(ovlp,(exc_size,exc_size))
518      !!! do mi=1,nstates
519      !!!   do ireh=1,nstates
520      !!!     ovlp(mi,ireh) = xdotc(exc_size,exc_vl(:,mi),1,exc_mat(:,ireh),1)
521      !!!     if(mi==ireh) then
522      !!!       !if(ABS(ovlp(mi,ireh)) < 0.999) then
523      !!!       !  write(*,*) "it,itp = ",mi,ireh,"ovlp = ",ovlp(mi,ireh)
524      !!!       !end if
525      !!!     else
526      !!!       if(ABS(ovlp(mi,ireh)) > 0.001) then
527      !!!         write(*,*) "it,itp = ",mi,ireh,"ovlp = ",ovlp(mi,ireh)
528      !!!       end if
529      !!!     end if
530      !!!   end do
531      !!! end do
532      !!! !call xgemm("C","N",exc_size,nstates,nstates,cone,exc_vl,exc_size,exc_mat,exc_size,czero,ovlp,nstates)
533 
534      !!! write(777,*) ovlp
535      !!! ABI_FREE(ovlp)
536      !!! !ENDDBYG
537 
538      !% fform = 1002 ! FIXME
539      !% call hdr_io_int(fform,Hdr_bse,2,eig_unt)
540 
541      write(eig_unt) do_ep_lifetime
542      write(eig_unt) exc_size, nstates
543      write(eig_unt) exc_ene_c(1:nstates)
544      do mi=1,nstates
545        write(eig_unt) exc_mat(1:exc_size,mi)
546        if(do_ep_lifetime) then
547          write(eig_unt) exc_vl(1:exc_size,mi)
548        end if
549      end do
550 
551      close(eig_unt)
552 
553    end do ! itemp
554 
555    ABI_FREE(exc_mat)
556    if (do_ep_renorm) then
557      ABI_FREE(exc_vl)
558      ABI_FREE(bs2eph)
559    end if
560 
561  CASE (.TRUE.)
562 
563 #if defined HAVE_LINALG_SCALAPACK && defined HAVE_MPI_IO
564    if (nsppol==2) then
565      MSG_WARNING("nsppol==2 + scalapack not coded yet")
566    end if
567 
568    istwf_k=1; tbloc=50
569    write(msg,'(2(a,i0))')". Using scaLAPACK version with nprocs= ",nprocs,"; block size= ",tbloc
570    call wrtout(std_out,msg,"PERS")
571    call wrtout(ab_out,msg,"COLL")
572 
573    write(msg,'(a,f8.1,a)')' Allocating excitonic eigenvalues. Memory required: ',exc_size*dp*b2Mb,' Mb. '
574    call wrtout(std_out,msg,"PERS")
575    !
576    ! Init scaLAPACK environment.
577    call init_scalapack(Slk_processor,comm)
578    !
579    ! Init scaLAPACK matrices
580    call init_matrix_scalapack(Slk_mat,exc_size,exc_size,Slk_processor,istwf_k,tbloc=tbloc)
581 
582    call init_matrix_scalapack(Slk_vec,exc_size,exc_size,Slk_processor,istwf_k,tbloc=tbloc)
583    !
584    ! Open the file with MPI-IO and skip the record.
585    amode=MPI_MODE_RDONLY
586 
587    call MPI_FILE_OPEN(comm, hreso_fname, amode, MPI_INFO_NULL, mpi_fh, ierr)
588    ABI_CHECK_MPI(ierr,"MPI_IO error opening file: "//TRIM(hreso_fname))
589 
590    ! Skip the header and find the offset for reading the matrix.
591    call exc_skip_bshdr_mpio(mpi_fh,xmpio_collective,ehdr_offset)
592    !
593    ! Read scaLAPACK matrix from the file.
594    if (nsppol==1) then
595      call slk_read(Slk_mat,"Upper","Hermitian",is_fortran_file,mpi_fh=mpi_fh,offset=ehdr_offset)
596    else
597      array_of_sizes = (/exc_size,exc_size/)
598      block_sizes(:,1) = (/neh1,neh1/)
599      block_sizes(:,2) = (/neh2,neh2/)
600      block_sizes(:,3) = (/neh1,neh2/)
601      MSG_ERROR("Not tested")
602      !call slk_read_from_blocks(Slk_mat,array_of_sizes,block_sizes,is_fortran_file,mpi_fh=mpi_fh,offset=ehdr_offset)
603    end if
604 
605    call MPI_FILE_CLOSE(mpi_fh, ierr)
606    ABI_CHECK_MPI(ierr,"FILE_CLOSE")
607 
608    if (do_full_diago) then
609      call wrtout(std_out," Performing full diagonalization with scaLAPACK...","COLL")
610 
611      call slk_pzheev("Vectors","Upper",Slk_mat,Slk_vec,exc_ene)
612    else
613      call wrtout(std_out," Performing partial diagonalization with scaLAPACK...","COLL")
614      il=1; iu=nstates; abstol=zero !ABSTOL = PDLAMCH(comm,'U')
615      call slk_pzheevx("Vectors","Index","Upper",Slk_mat,vl,vu,il,iu,abstol,Slk_vec,mene_found,exc_ene)
616    end if
617 
618    exc_ene_c(:) = exc_ene(:)
619 
620    call destruction_matrix_scalapack(Slk_mat)
621 
622    call wrtout(std_out,' Writing eigenvalues/vectors on file: '//TRIM(bseig_fname),"COLL", do_flush=.True.)
623 
624    ! Write distributed matrix on file bseig_fname with Fortran records.
625    if (my_rank==master) then ! Write exc eigenvalues. Vectors will be appended in slk_write.
626      if (open_file(bseig_fname,msg,newunit=eig_unt,form="unformatted",action="write") /= 0) then
627        MSG_ERROR(msg)
628      end if
629      write(eig_unt) exc_size, nstates
630      write(eig_unt) exc_ene_c(1:nstates)
631      close(eig_unt)
632    end if
633 
634    call xmpi_barrier(comm)
635    !
636    ! Open the file with MPI-IO and skip the record.
637    amode=MPI_MODE_RDWR
638 
639    call MPI_FILE_OPEN(comm, bseig_fname, amode, MPI_INFO_NULL, mpi_fh, ierr)
640    ABI_CHECK_MPI(ierr,"MPI_IO error opening file: "//TRIM(hreso_fname))
641 
642    !call MPI_FILE_SYNC(mpi_fh,ierr)
643 
644    ehdr_offset = 0
645    call xmpio_read_frm(mpi_fh,ehdr_offset,xmpio_collective,fmarker,ierr)
646    call xmpio_read_frm(mpi_fh,ehdr_offset,xmpio_collective,fmarker,ierr)
647 
648    write(std_out,*)"Writing nstates ",nstates
649    gsub(:,1) = (/1,1/)
650    gsub(:,2) = (/exc_size,nstates/)
651    call slk_write(Slk_vec,"All",is_fortran_file,mpi_fh=mpi_fh,offset=ehdr_offset,glob_subarray=gsub)
652 
653    call MPI_FILE_CLOSE(mpi_fh, ierr)
654    ABI_CHECK_MPI(ierr,"FILE_CLOSE")
655 
656    call destruction_matrix_scalapack(Slk_vec)
657    call end_scalapack(Slk_processor)
658    call xmpi_barrier(comm)
659 #else
660    MSG_BUG("You should not be here!")
661 #endif
662 
663  END SELECT
664 
665  ! Order the eigenvalues
666  do ii=nstates,2,-1
667    do j=1,ii-1
668      if (DBLE(exc_ene_c(j)) > DBLE(exc_ene_c(j+1))) then
669        ctemp = exc_ene_c(j)
670        exc_ene_c(j) = exc_ene_c(j+1)
671        exc_ene_c(j+1) = ctemp
672      end if
673    end do
674  end do
675 
676 
677  write(msg,'(a,i4)')' Excitonic eigenvalues in eV up to n= ',nene_printed
678  call wrtout(std_out,msg,"PERS")
679  call wrtout(ab_out,msg,"COLL")
680 
681  do it=0,(nene_printed-1)/8
682    write(msg,'(8f10.5)') ( DBLE(exc_ene_c(ii))*Ha_eV, ii=1+it*8,MIN(it*8+8,nene_printed) )
683    call wrtout(std_out,msg,"PERS")
684    call wrtout(ab_out,msg,"COLL")
685  end do
686 
687  exc_gap    = MINVAL(DBLE(exc_ene_c(1:nstates)))
688  exc_maxene = MAXVAL(DBLE(exc_ene_c(1:nstates)))
689 
690  write(msg,'(a,2(a,f7.2,2a),a)')ch10,&
691 &  " First excitonic eigenvalue= ",exc_gap*Ha_eV,   " [eV]",ch10,&
692 &  " Last  excitonic eigenvalue= ",exc_maxene*Ha_eV," [eV]",ch10,ch10
693  call wrtout(std_out,msg,"COLL",do_flush=.True.)
694  call wrtout(ab_out,msg,"COLL",do_flush=.True.)
695 
696  ABI_FREE(exc_ene_c)
697  ABI_FREE(exc_ene)
698 
699 10 call xmpi_barrier(comm)
700 
701  DBG_EXIT("PERS")
702 
703 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.

PARENTS

      m_exc_diago

CHILDREN

      destruction_matrix_scalapack,end_scalapack,exc_fullh_from_blocks
      exc_read_bshdr,exc_skip_bshdr_mpio,hermitianize,idx_glob
      init_matrix_scalapack,init_scalapack,mpi_file_close,mpi_file_open
      mpi_file_read_all,mpi_file_set_view,mpi_type_free,slk_pzgemm
      slk_pzhegvx,slk_single_fview_read_mask,slk_write,slk_zinvert,wrtout
      xgemm,xhdp_invert,xhegv,xhegvx,xmpi_barrier,xmpio_read_frm

SOURCE

736 subroutine exc_print_eig(BSp,bseig_fname,gw_gap,exc_gap)
737 
738 
739 !This section has been created automatically by the script Abilint (TD).
740 !Do not modify the following lines by hand.
741 #undef ABI_FUNC
742 #define ABI_FUNC 'exc_print_eig'
743 !End of the abilint section
744 
745  implicit none
746 
747 !Arguments ------------------------------------
748 !scalars
749  complex(dpc),intent(in) :: gw_gap
750  complex(dpc),intent(out) :: exc_gap
751  character(len=*),intent(in) :: bseig_fname
752  type(excparam),intent(in) :: BSp
753 
754 !Local variables ------------------------------
755 !scalars
756  integer :: nstates_read,ii,j,k,eig_unt,ieig,hsize_exp
757  integer :: hsize_read !,nstates
758  complex(dpc) :: bind_energy,ctemp
759  character(len=500) :: msg
760  !type(Hdr_type) :: tmp_Hdr
761 !arrays
762  integer,allocatable :: iperm(:)
763  real(dp),allocatable :: exc_rene(:)
764  complex(dpc),allocatable :: exc_cene(:)
765 
766 !************************************************************************
767 
768  exc_gap = czero
769 
770  if (open_file(bseig_fname,msg,newunit=eig_unt,form="unformatted",status="old",action="read") /= 0) then
771    MSG_ERROR(msg)
772  end if
773 
774  read(eig_unt) ! do_ep_lifetime
775  read(eig_unt) hsize_read, nstates_read
776 
777  if (BSp%use_coupling==0) hsize_exp =   SUM(Bsp%nreh)
778  if (BSp%use_coupling>0)  hsize_exp = 2*SUM(Bsp%nreh)
779 
780  if (hsize_exp /= hsize_read) then
781    write(msg,'(2(a,i0))')" Wrong dimension: read: ",hsize_read," expected= ",hsize_exp
782    MSG_ERROR(msg)
783  end if
784 
785  ABI_MALLOC(exc_cene,(nstates_read))
786  read(eig_unt) exc_cene(:)
787 
788  ABI_MALLOC(exc_rene,(nstates_read))
789  exc_rene = DBLE(exc_cene)
790 
791  ABI_MALLOC(iperm,(nstates_read))
792  iperm = (/(ii, ii=1,nstates_read)/)
793 
794  call sort_dp(nstates_read,exc_rene,iperm,tol6)
795 
796  ABI_FREE(exc_rene)
797  ABI_FREE(iperm)
798 
799  ! put in ascending order
800  do ii=nstates_read,2,-1
801    do j=1,ii-1
802      if (DBLE(exc_cene(j)) > DBLE(exc_cene(j+1))) then
803        ctemp = exc_cene(j)
804        exc_cene(j) = exc_cene(j+1)
805        exc_cene(j+1) = ctemp
806      end if
807    end do
808  end do
809 
810  exc_gap = DCMPLX(ABS(DBLE(exc_cene(1))),AIMAG(exc_cene(1)))
811 
812  do ii=1,nstates_read
813    if (ABS(DBLE(exc_cene(ii))) < DBLE(exc_gap)) then
814      exc_gap = DCMPLX(ABS(DBLE(exc_cene(ii))),AIMAG(exc_cene(ii)))
815    end if
816  end do
817 
818  bind_energy = gw_gap - exc_gap
819 
820  write(msg,"(3(a,2f6.2,2a))")&
821 &  " GW  direct gap     ",gw_gap*Ha_eV,     " [eV] ",ch10,&
822 &  " EXC direct gap     ",exc_gap*Ha_eV,    " [eV] ",ch10,&
823 &  " EXC binding energy ",bind_energy*Ha_eV," [eV] ",ch10
824  call wrtout(std_out,msg,"COLL")
825  call wrtout(ab_out,msg,"COLL")
826 
827  msg=' Excitonic eigenvalues up to the GW energy gap [eV]'
828  call wrtout(std_out,msg,"COLL")
829  call wrtout(ab_out,msg,"COLL")
830 
831  do ii=1,nstates_read
832    if (DBLE(exc_cene(ii)) > zero) EXIT
833  end do
834 
835  do j=ii,nstates_read
836    if (DBLE(exc_cene(j)) > DBLE(gw_gap)) EXIT
837  end do
838  j=j-1
839 
840  do ieig=ii,j
841    write(msg,'(i3,a,2f6.2,a)')ieig," (",exc_cene(ieig)*Ha_eV,")"
842    call wrtout(std_out,msg,"COLL")
843    call wrtout(ab_out,msg,"COLL")
844  end do
845 
846  ii=ii-1
847  do j=ii,1,-1
848    if (ABS(DBLE(exc_cene(j))) > DBLE(gw_gap)) EXIT
849  end do
850  j=j+1
851 
852  ! This coding is not portable, write to ab_out has been disabled.
853  if (ii>0) then
854    do k=ii,j,-1
855      write(msg,'(i3,a,2f6.2,a)')k," (",exc_cene(k)*Ha_eV,")"
856      call wrtout(std_out,msg,"COLL")
857    end do
858  end if
859 
860  ABI_FREE(exc_cene)
861 
862  close(eig_unt)
863 
864 end subroutine exc_print_eig