TABLE OF CONTENTS


ABINIT/m_bse_io [ Modules ]

[ Top ] [ Modules ]

NAME

 m_bse_io

FUNCTION

  This module provides routines to read the Bethe-Salpeter Hamiltonian from file

COPYRIGHT

  Copyright (C) 2008-2018 ABINIT group (MG)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 MODULE m_bse_io
28 
29  use defs_basis
30  use m_xmpi
31  use m_errors
32  use m_abicore
33 #if defined HAVE_MPI2
34  use mpi
35 #endif
36 #if defined HAVE_NETCDF
37  use netcdf
38 #endif
39  use m_nctk
40  use iso_c_binding
41  use m_hdr
42 
43  use defs_abitypes,    only : Hdr_type
44  use m_time,           only : cwtime
45  use m_fstrings,       only : toupper
46  use m_io_tools,       only : open_file
47  use m_numeric_tools,  only : arth
48  use m_special_funcs,  only : dirac_delta
49  use m_bs_defs,        only : excparam
50  use m_bz_mesh,        only : kmesh_t
51 
52  implicit none
53 
54 #if defined HAVE_MPI1
55  include 'mpif.h'
56 #endif
57 
58  private
59 
60  public  :: exc_read_rblock_fio      ! Reads the entire resonant sub-block from file using Fortran IO.
61  public  :: exc_read_rcblock         ! Reads a distributed sub-block of the excitonic Hamiltonian from file.
62  public  :: exc_fullh_from_blocks    ! Initialize the specified sub-blocks of the *full* matrix (reso+anti-reso) from file.
63  public  :: rrs_of_glob              ! [+1,-1,0] if (row_glob,col_glob) belongs to the [ resonant, anti-resonant, (anti)coupling block ]
64  public  :: ccs_of_glob              ! [+1,-1,0] if (row_glob,col_glob) belongs to the [ coupling, anti-coupling, (anti)resonant block ]
65  public  :: offset_in_file           ! Function used to describe the way the Hamiltonian is stored on disk.
66  public  :: exc_write_bshdr          ! Writes the Header of the (BSR|BSC) files storing the excitonic Hamiltonian.
67  public  :: exc_read_bshdr           ! Reads the Header of the (BSR|BSC) files.
68  public  :: exc_skip_bshdr           ! Skip the Header of the (BSR|BSC) files. Fortran version.
69  public  :: exc_skip_bshdr_mpio      ! Skip the Header of the (BSR|BSC) files. MPI-IO  version.
70  public  :: exc_read_eigen           ! Read selected energies and eigenvectors from the BSEIG file.
71  public  :: exc_amplitude            ! Calculate the amplitude function F(w) = \sum_t |<t|exc_vec>|^2 \delta(ww- ene_t) where t is the eh transition.
72  public  :: exc_write_optme          ! Writes the OME file storing the optical matrix elements
73  public  :: exc_ham_ncwrite          ! Writes the hamiltonian in NETCDF format

m_bse_io/ccs_of_glob [ Functions ]

[ Top ] [ m_bse_io ] [ Functions ]

NAME

FUNCTION

  [+1,-1,0] if (row_glob,col_glob) belongs to the [ coupling, anti-coupling, (anti)resonant block ]

INPUTS

OUTPUT

PARENTS

SOURCE

1142 pure function ccs_of_glob(row_glob,col_glob,size_glob)
1143 
1144 
1145 !This section has been created automatically by the script Abilint (TD).
1146 !Do not modify the following lines by hand.
1147 #undef ABI_FUNC
1148 #define ABI_FUNC 'ccs_of_glob'
1149 !End of the abilint section
1150 
1151  implicit none
1152 
1153 !Arguments ------------------------------------
1154  integer :: ccs_of_glob
1155  integer,intent(in) :: row_glob,col_glob
1156  integer,intent(in) :: size_glob(2)
1157 
1158 !Local variables ------------------------------
1159 !scalars
1160  integer :: nreh1,nreh2
1161 
1162 ! *************************************************************************
1163 
1164  nreh1=size_glob(1)/2 ! Matrix is square and nreh1==nreh2 but oh well.
1165  nreh2=size_glob(2)/2
1166 
1167  if (row_glob<=nreh1 .and. col_glob >nreh2) then      ! Coupling.
1168    ccs_of_glob = +1
1169  else if (row_glob >nreh1 .and. col_glob<=nreh2) then ! anti-Coupling
1170    ccs_of_glob = -1
1171  else
1172    ccs_of_glob = 0
1173  end if
1174 
1175 end function ccs_of_glob

m_bse_io/exc_amplitude [ Functions ]

[ Top ] [ m_bse_io ] [ Functions ]

NAME

  exc_amplitude

FUNCTION

  Calculate the amplitude function of the excitonic eigenstate |exc_vec\>
    F(w) = \sum_t |<t|exc_vec>|^2 \delta(ww- ene_t) where the sum over t is done
  of the full set of transitions used to construct the BS Hamiltoniam.

INPUTS

  Bsp<excparam>=Structure storing the parameters of the run.
  eig_fname=The name of the file storing the excitonic eigenvectors.
  nvec=Number of excitonic states to analyze.
  vec_idx(nvec)=List with the indeces of the excitonic states sorted in ascending order.
  out_fname=The name of the file where the results are written.

OUTPUT

  Only writing.

PARENTS

CHILDREN

      c_f_pointer

SOURCE

1401 subroutine exc_amplitude(Bsp,eig_fname,nvec,vec_idx,out_fname)
1402 
1403 
1404 !This section has been created automatically by the script Abilint (TD).
1405 !Do not modify the following lines by hand.
1406 #undef ABI_FUNC
1407 #define ABI_FUNC 'exc_amplitude'
1408 !End of the abilint section
1409 
1410  implicit none
1411 
1412 !Arguments ------------------------------------
1413 !scalars
1414  integer,intent(in) :: nvec
1415  character(len=*),intent(in) :: eig_fname,out_fname
1416  type(excparam),intent(in) :: BSp
1417 ! arrays
1418  integer,intent(in) :: vec_idx(nvec)
1419 
1420 !Local variables ------------------------------
1421 !scalars
1422  integer :: vec,art_idx,ierr
1423  integer :: spin,iw,it,nw,pos_w,neg_w,out_unt,rt_idx,hsize
1424  real(dp) :: ene_rt,ampl_eh,ampl_he,xx,stdev,w_max,w_min,step
1425  character(len=500) :: msg
1426 !arrays
1427  real(dp),allocatable :: wmesh(:),amplitude(:),ene_list(:)
1428  complex(dpc),allocatable :: vec_list(:,:)
1429 
1430  ! *************************************************************************
1431 
1432  ! Setup of the frequency mesh for F(w).
1433  w_min=greatest_real; w_max=smallest_real
1434  do spin=1,BSp%nsppol
1435    do it=1,BSp%nreh(spin)
1436      ene_rt  = Bsp%Trans(it,spin)%en
1437      w_min = MIN(w_min,ene_rt)
1438      w_max = MAX(w_max,ene_rt)
1439    end do
1440  end do
1441 
1442  step = Bsp%domega
1443  if (Bsp%use_coupling==0) then
1444    nw = (w_max - w_min)/step + 1
1445    ABI_MALLOC(wmesh,(nw))
1446    wmesh = arth(w_min,step,nw)
1447  else
1448    ! Both positive and negative frequencies are needed.
1449    pos_w = (w_max - w_min)/step + 1
1450    neg_w = pos_w; if (ABS(w_min) < tol6) neg_w=neg_w-1 ! zero should not included twice.
1451    nw = neg_w + pos_w
1452    ABI_MALLOC(wmesh,(nw))
1453    wmesh(1:neg_w)  = arth(-w_max,step,neg_w)
1454    wmesh(neg_w+1:) = arth( w_min,step,pos_w)
1455  end if
1456  !
1457  ! Read selected eigenvectors.
1458  hsize = SUM(Bsp%nreh); if (Bsp%use_coupling>0) hsize=2*hsize
1459 
1460  ABI_MALLOC(ene_list,(nvec))
1461  ABI_STAT_MALLOC(vec_list,(hsize,nvec), ierr)
1462  ABI_CHECK(ierr==0, "out of memory in vec_list")
1463 
1464  call exc_read_eigen(eig_fname,hsize,nvec,vec_idx,vec_list,ene_list,Bsp=Bsp)
1465 
1466  ABI_FREE(ene_list)
1467 
1468  if (open_file(out_fname,msg,newunit=out_unt,form="formatted",action="write") /= 0) then
1469    MSG_ERROR(msg)
1470  end if
1471 
1472  write(out_unt,*)"# Amplitude functions F(w) = \sum_t |<t|exc_vec>|^2 \delta(ww- ene_t), w is given in eV. "
1473  write(out_unt,*)"# Number of excitonic eigenvectors analyzed: ",nvec
1474 
1475  ABI_MALLOC(amplitude,(nw))
1476  stdev = BSp%broad ! Broadening for the gaussian.
1477 
1478  do vec=1,nvec
1479    !
1480    ! amplitude(ww) = \sum_t |<t|exc_vec>|^2 \delta(ww- ene_t)
1481    amplitude = zero
1482    do spin=1,BSp%nsppol
1483      do it=1,BSp%nreh(spin)
1484       ene_rt  = Bsp%Trans(it,spin)%en
1485       rt_idx  = it + (spin-1)*Bsp%nreh(1)
1486       ampl_eh = (ABS(vec_list(rt_idx,vec)))**2
1487       if (Bsp%use_coupling>0) then ! Need the index and the amplitude of the anti-resonant component.
1488         art_idx  = it + (spin-1)*Bsp%nreh(1) + SUM(Bsp%nreh)
1489         ampl_he = (ABS(vec_list(art_idx,vec)))**2
1490       end if
1491       !
1492       do iw=1,nw ! Accumulate
1493         xx = wmesh(iw) - ene_rt
1494         amplitude(iw) = amplitude(iw) + ampl_eh * dirac_delta(xx,stdev)
1495         if (Bsp%use_coupling>0) then
1496           xx = wmesh(iw) + ene_rt
1497           amplitude(iw) = amplitude(iw) + ampl_he * dirac_delta(xx,stdev)
1498         end if
1499       end do
1500       !
1501      end do
1502    end do
1503    !
1504    ! Write results
1505    write(out_unt,*)"# Amplitude function F(w) for exc_vector index ",vec_idx(vec)
1506    do iw=1,nw
1507      write(out_unt,*)wmesh(iw)*Ha_eV,amplitude(iw)
1508    end do
1509    write(out_unt,*)"#"
1510  end do
1511 
1512  close(out_unt)
1513 
1514  ABI_FREE(amplitude)
1515  ABI_FREE(wmesh)
1516  ABI_FREE(vec_list)
1517 
1518 end subroutine exc_amplitude

m_bse_io/exc_fullh_from_blocks [ Functions ]

[ Top ] [ m_bse_io ] [ Functions ]

NAME

  exc_fullh_from_blocks

FUNCTION

   Construct the matrix F H

INPUTS

  funt
  nsppol
  block_type
    "Resonant"
    "Coupling"
  row_sign
   -1 to read ( R   C )
              (-C* -R*)

   +1 to read ( R   C )
              ( C*  R*)
  diago_is_real=Used when block_type=resonat to specify wheter the diagonal matrix elements are
  real or complex (when QP linewidth are included)
  nreh(nsppol)
  exc_size=Size of the full excitonic Hamiltonian.

SIDE EFFECTS

  exc_ham(exc_size,exc_size)

NOTES

PARENTS

      m_exc_diago

CHILDREN

      c_f_pointer

SOURCE

 829 subroutine exc_fullh_from_blocks(funt,block_type,nsppol,row_sign,diago_is_real,nreh,exc_size,exc_ham)
 830 
 831 
 832 !This section has been created automatically by the script Abilint (TD).
 833 !Do not modify the following lines by hand.
 834 #undef ABI_FUNC
 835 #define ABI_FUNC 'exc_fullh_from_blocks'
 836 !End of the abilint section
 837 
 838  implicit none
 839 
 840 !Arguments ------------------------------------
 841 !scalars
 842  integer,intent(in) :: funt,exc_size,nsppol,row_sign
 843  logical,intent(in) :: diago_is_real
 844  character(len=*),intent(in) :: block_type
 845 !arrays
 846  integer,intent(in) :: nreh(nsppol)
 847  complex(dpc),intent(inout) :: exc_ham(exc_size,exc_size)
 848 
 849 !Local variables-------------------------------
 850 !scalars
 851  integer :: it,itp,szbuf,neh,pad_c1,pad_r1,spin_dim,spad_r,spad_c
 852  integer :: block,spad,row1,col1,row2,col2,spin_stride,ierr
 853  complex(dpc) :: cttp
 854  character(len=500) :: errmsg
 855 !arrays
 856  complex(dpc),allocatable :: cbuff_dpc(:)
 857 
 858 ! *********************************************************************
 859 
 860  szbuf=exc_size ! FIXME oversized!
 861  neh = nreh(1)
 862 
 863  if (nsppol==2) then
 864    MSG_WARNING("nsppol==2 is very experimental")
 865  end if
 866  if (ANY(nreh(1)/=nreh)) then
 867    MSG_ERROR(" different nreh are not supported")
 868  end if
 869 
 870  ABI_STAT_MALLOC(cbuff_dpc,(exc_size), ierr)
 871  ABI_CHECK(ierr==0, "out of memory in cbuff_dpc")
 872  !
 873  ! The two cases nsppol==1,2 can be merged but we keep them
 874  ! separated to keep to code readable.
 875 
 876  SELECT CASE (toupper(block_type))
 877  CASE ("RESONANT")
 878 
 879    if (nsppol==1) then
 880      !
 881      ! Construct resonant block from the upper triangle stored on file.
 882      neh = nreh(1)
 883      do itp=1,neh
 884        read(funt, err=10, iomsg=errmsg) cbuff_dpc(1:itp)
 885        do it=1,itp-1 ! The diagonal is treated below.
 886          cttp = cbuff_dpc(it)
 887          exc_ham(it     ,itp)     =                cttp   ! R
 888          exc_ham(itp    ,it )     =          CONJG(cttp)  ! row_sign R*
 889          exc_ham(neh+it ,neh+itp) = row_sign*CONJG(cttp)
 890          exc_ham(neh+itp,neh+it ) = row_sign*cttp
 891        end do
 892        if (diago_is_real) then
 893          exc_ham(itp    ,itp)     =          DBLE(cbuff_dpc(itp))
 894          exc_ham(neh+itp,neh+itp) = row_sign*DBLE(cbuff_dpc(itp))
 895        else
 896          exc_ham(itp,itp)         =           cbuff_dpc(itp)
 897          exc_ham(neh+itp,neh+itp) = row_sign*(CONJG(cbuff_dpc(itp)))
 898        end if
 899      end do
 900      !
 901      !
 902    else
 903      ! FIXME this part wont work if we have a different number of e-h pairs
 904      ABI_CHECK(ALL(nreh==nreh(1)),"Different number of transitions")
 905      ! The file contains
 906      ! A) The up-up and the down-down block in packed form
 907      ! (only the upper triangle is stored since these blocks are Hermitian)
 908      ! B) The entire up-down exchange block (no symmetry here)
 909      !
 910      ! The resonant block is given by:
 911      !     |  (v'c' up)    | (v'c' dwn) |
 912      !     ------------------------------           where v_{-+} = v_{+-}^H when the momentum of the photon is neglected.
 913      !     | [T-W+v]++     |      v+-   | (vc up)   Note that v_{+-} is not Hermitian due to the presence of different spins.
 914      ! R = ------------------------------           Actually it reduces to a Hermitian matrix when the system is not spin polarized.
 915      !     |     v-+       | [T-W+v]--  | (vc dwn)  [T-W+v] is Hermitian provided the the QP energies are purely real.
 916      !     ------------------------------
 917      !
 918      ! *) Fill the diagonal blocks.
 919      !    only the upper triangle is stored on file.
 920      !    row1,col1 refer to the resonant block.
 921      !    row2,col2 refer to the anti-resonant block.
 922      do block=1,2
 923        if (block==1) then
 924          spad=0
 925          spin_stride=SUM(nreh)
 926        end if
 927        if (block==2) then
 928          spad=nreh(1)
 929          spin_stride=2*nreh(1) + nreh(2)
 930        end if
 931        do itp=1,nreh(block)
 932          read(funt, err=10, iomsg=errmsg) cbuff_dpc(1:itp)
 933          col1 = itp+spad
 934          col2 = itp+spin_stride
 935          do it=1,itp-1
 936            cttp = cbuff_dpc(it)
 937            row1 = it + spad
 938            row2 = it + spin_stride
 939            exc_ham(row1,col1) =                cttp    ! [T-W+v]
 940            exc_ham(col1,row1) =          CONJG(cttp)
 941            exc_ham(row2,col2) = row_sign*CONJG(cttp)   ! row_sign [T-W+v]*
 942            exc_ham(col2,row2) = row_sign*cttp
 943          end do
 944          if (diago_is_real) then
 945            exc_ham(col1,col1) =          DBLE(cbuff_dpc(itp))
 946            exc_ham(col2,col2) = row_sign*DBLE(cbuff_dpc(itp))
 947          else
 948            exc_ham(col1,col1) = cbuff_dpc(itp)
 949            exc_ham(col2,col2) = row_sign*CONJG(cbuff_dpc(itp))
 950          end if
 951        end do
 952      end do
 953      !
 954      ! Read v+- and reconstruct resonant and anti-resonat blocks.
 955      ! TODO recheck this
 956      pad_r1=SUM(nreh)
 957      pad_c1=2*nreh(1) + nreh(2)
 958 
 959      spin_dim=nreh(1)
 960      do itp=1,spin_dim
 961        read(funt, err=10, iomsg=errmsg) cbuff_dpc(1:spin_dim)
 962        exc_ham(1:spin_dim,nreh(1)+itp) = cbuff_dpc(1:spin_dim)                              ! upper reso
 963        exc_ham(1+pad_r1:pad_r1+spin_dim,pad_c1+itp) = row_sign*CONJG(cbuff_dpc(1:spin_dim)) ! upper anti-reso
 964        col1 = itp+nreh(1)
 965        col2 = itp+(2*nreh(1) + nreh(2))
 966        do it=1,spin_dim
 967          cttp = cbuff_dpc(it)
 968          row1 = it
 969          exc_ham(col1,row1) =    CONJG(cttp)  ! lower reso.
 970          row2 = it + SUM(nreh)
 971          exc_ham(col2,row2) = row_sign*cttp   ! lower anti-reso.
 972        end do
 973      end do
 974    end if
 975 
 976  CASE ("COUPLING")
 977    !
 978    if (nsppol==1) then
 979      !
 980      do itp=1,neh
 981        read(funt, err=10, iomsg=errmsg) cbuff_dpc(1:itp)
 982        do it=1,itp
 983          cttp = cbuff_dpc(it)
 984          exc_ham(it    ,neh+itp) = cttp
 985          exc_ham(itp   ,neh+it ) = cttp
 986          exc_ham(neh+it ,itp   ) = row_sign*CONJG(cttp)
 987          exc_ham(neh+itp,it    ) = row_sign*CONJG(cttp)
 988        end do
 989      end do
 990      !
 991    else
 992      !  The coupling block is given by:
 993      !      |  (c'v' up)   |    (c'v dwn)     |
 994      !      -----------------------------------           where v_{-+} = v_{+-}^t when the momentum of the photon is neglected.
 995      !      | [-W+v]++     |      v+-         | (vc up)   The entire matrix v_{+-} is stored on file.
 996      !  C = -----------------------------------
 997      !      |     v-+      |    [-W+v]--      | (vc dwn)
 998      !      -----------------------------------
 999      !
1000      ! *) Fill blocks that are diagonal in spin coordinates.
1001      !    row1,col1 refer to the resonat block.
1002      !    row2,col2 refer to the anti-resonant block.
1003      do block=1,2
1004        if (block==1) then
1005          spad_r=0
1006          spad_c=SUM(nreh)
1007        end if
1008        if (block==2) then
1009          spad_r=nreh(1)
1010          spad_c=2*nreh(1)+nreh(2)
1011        end if
1012        do itp=1,nreh(block)
1013          read(funt, err=10, iomsg=errmsg) cbuff_dpc(1:itp)
1014          col1 = itp+spad_c
1015          row2 = itp+spad_r
1016          do it=1,itp-1
1017            cttp = cbuff_dpc(it)
1018            row1 = it + spad_r
1019            col2 = it + spad_c
1020            exc_ham(row1,col1) =                cttp  ! upper coupling
1021            exc_ham(row2,col2) =                cttp  ! lower coupling (symmetric)
1022            exc_ham(col1,row1) = row_sign*CONJG(cttp) ! lower anti-coupling
1023            exc_ham(col2,row2) = row_sign*CONJG(cttp) ! upper anti-couling
1024          end do
1025          ! TODO recheck this
1026          row1 = itp+spad_r
1027          exc_ham(row1,col1) = cbuff_dpc(itp)                  ! Diagonals of the block.
1028          col2 = itp+spad_c
1029          exc_ham(col2,row2) = row_sign*CONJG(cbuff_dpc(itp))
1030        end do
1031      end do
1032      !
1033      ! Read Full v+- and reconstruct resonant and anti-resonat blocks.
1034      ! TODO recheck this
1035      spad=2*nreh(1) + nreh(2)
1036      pad_r1=SUM(nreh)
1037      pad_c1=2*nreh(1) + nreh(2)
1038 
1039      spin_dim=nreh(1)
1040      do itp=1,spin_dim
1041        read(funt, err=10, iomsg=errmsg) cbuff_dpc(1:spin_dim)
1042        exc_ham(1:spin_dim,spad+itp) = cbuff_dpc(1:spin_dim)                                  ! upper block reso
1043        exc_ham(1+pad_r1:pad_r1+spin_dim,nreh(1)+itp) = row_sign*CONJG(cbuff_dpc(1:spin_dim)) ! upper block anti-reso
1044        col1 = itp+spad
1045        row2 = itp+nreh(1)
1046        do it=1,spin_dim
1047          cttp = cbuff_dpc(it)
1048          row1 = it
1049          exc_ham(col1,row1) =  row_sign*CONJG(cttp)  ! lower anti-reso.
1050          col2 = it + SUM(nreh)
1051          exc_ham(row2,col2) = cttp                   ! lower reso.
1052        end do
1053      end do
1054      !
1055    end if
1056 
1057  CASE DEFAULT
1058    MSG_ERROR("Unknown block_type: "//TRIM(block_type))
1059  END SELECT
1060 
1061  ABI_FREE(cbuff_dpc)
1062 
1063  return
1064 
1065 ! Handle IO Error
1066 10 continue
1067  MSG_ERROR(errmsg)
1068 
1069 end subroutine exc_fullh_from_blocks

m_bse_io/exc_read_bshdr [ Functions ]

[ Top ] [ m_bse_io ] [ Functions ]

NAME

  exc_read_bshdr

FUNCTION

  Reads the header of the (BSR|BSC) files storing the excitonic Hamiltonian.
  and performs basilar consistency checks.

INPUTS

  funt=Unit number.
  Bsp<excparam>=Structure storing the parameters of the run.
  Hdr<hdr_type>=The abinit header.

OUTPUT

  fform=Integer defining the file format.
  ierr=Status error.

PARENTS

      m_bse_io,m_exc_diago

CHILDREN

      c_f_pointer

SOURCE

167 subroutine exc_read_bshdr(funt,Bsp,fform,ierr)
168 
169 
170 !This section has been created automatically by the script Abilint (TD).
171 !Do not modify the following lines by hand.
172 #undef ABI_FUNC
173 #define ABI_FUNC 'exc_read_bshdr'
174 !End of the abilint section
175 
176  implicit none
177 
178  !Arguments ------------------------------------
179  integer,intent(in) :: funt
180  integer,intent(out) :: fform,ierr
181  type(excparam),intent(in) :: BSp
182 
183 !Local variables ------------------------------
184 !scalars
185  integer :: nkbz_read
186  character(len=500) :: errmsg
187  type(hdr_type) :: Hdr
188 !arrays
189  integer :: nreh_read(SIZE(BSp%nreh))
190 
191  ! *************************************************************************
192 
193  ierr=0
194 
195  ! Read the header and perform consistency checks.
196  call hdr_fort_read(hdr, funt, fform, rewind=.True.)
197  ABI_CHECK(fform /= 0, "hdr_fort_read returned fform == 0")
198 
199  read(funt, err=10, iomsg=errmsg) nreh_read, nkbz_read
200 
201  call hdr_free(Hdr)
202 
203  if (ANY(nreh_read/=BSp%nreh)) then
204    call wrtout(std_out,"Wrong number of e-h transitions","COLL")
205    ierr = ierr + 1
206  end if
207 
208  return
209 
210 10 ierr = 1
211  MSG_WARNING(errmsg)
212 
213 end subroutine exc_read_bshdr

m_bse_io/exc_read_eigen [ Functions ]

[ Top ] [ m_bse_io ] [ Functions ]

NAME

  exc_read_eigen

FUNCTION

  Read selected energies and eigenvectors from the BSEIG file.

INPUTS

  eig_fname=The name of the file storing the excitonic eigenvectors.
  hsize=Size of the Hamiltonian.
  nvec=Number of excitonic states to analyze.
  vec_idx(nvec)=List with the indeces of the excitonic states sorted in ascending order.
  [Bsp]<excparam>=Structure storing the parameters of the run. If present the
    routine will perform additional consistency checks to make sure that
    the content of the file is consistent with the present run.

OUTPUT

  [ene_list(nvec)]=Excitonic energies
  vec_list(hsize,nvec)=Excitonic eigenvectors.

PARENTS

      exc_plot,m_bse_io

CHILDREN

      c_f_pointer

SOURCE

365 subroutine exc_read_eigen(eig_fname,hsize,nvec,vec_idx,vec_list,ene_list,Bsp)
366 
367 
368 !This section has been created automatically by the script Abilint (TD).
369 !Do not modify the following lines by hand.
370 #undef ABI_FUNC
371 #define ABI_FUNC 'exc_read_eigen'
372 !End of the abilint section
373 
374  implicit none
375 
376 !Arguments ------------------------------------
377 !scalars
378  integer,intent(in) :: nvec,hsize
379  character(len=*),intent(in) :: eig_fname
380  type(excparam),optional,intent(in) :: BSp
381 ! arrays
382  integer,intent(in) :: vec_idx(nvec)
383  real(dp),optional,intent(out) :: ene_list(nvec)
384  complex(dpc),intent(out) :: vec_list(hsize,nvec)
385 
386 !Local variables ------------------------------
387 !scalars
388  integer :: eig_unt,hsize_read,neig_read,ll,vec
389  character(len=500) :: msg,errmsg
390 !arrays
391  !real(dp),allocatable :: exc_ene(:)
392  complex(dpc),allocatable :: exc_ene_cplx(:)
393 
394  ! *************************************************************************
395 
396  ABI_UNUSED(BSp%nline)
397 
398  if (open_file(eig_fname,msg,newunit=eig_unt,form="unformatted",status="old",action="read") /= 0) then
399    MSG_ERROR(msg)
400  end if
401 
402  read(eig_unt, err=10, iomsg=errmsg)hsize_read,neig_read
403  !write(std_out,*)hsize_read, neig_read
404 
405  if (hsize_read/=hsize) then
406    write(msg,'(a,2(1x,i0))')" hsize_read/=hsize: ",hsize_read,hsize
407    MSG_ERROR(msg)
408  end if
409 
410  ! Read eigenvalues, ignore possibly small imaginary part.
411  ABI_MALLOC(exc_ene_cplx,(neig_read))
412  read(eig_unt, err=10, iomsg=errmsg) exc_ene_cplx
413 
414  if (PRESENT(ene_list)) then
415    do vec=1,nvec
416      ll = vec_idx(vec)
417      ene_list(vec) = DBLE(exc_ene_cplx(ll))
418    end do
419  end if
420  ABI_FREE(exc_ene_cplx)
421 
422  vec=1
423  do ll=1,neig_read ! Read the selected excitons.
424    if (ll==vec_idx(vec))  then
425      read(eig_unt, err=10, iomsg=errmsg) vec_list(:,vec)
426      if (vec==nvec) EXIT
427      vec=vec+1
428    else
429      read(eig_unt, err=10, iomsg=errmsg)
430    end if
431  end do
432 
433  close(eig_unt, err=10, iomsg=errmsg)
434 
435  if (vec/=nvec) then
436    write(msg,'(a,2(1x,i0))')" vec_idx is wrong, vec/=nvec ",vec,nvec+1
437    MSG_ERROR(msg)
438  end if
439 
440  return
441 
442  ! Handle IO-error
443 10 continue
444  MSG_ERROR(errmsg)
445 
446 end subroutine exc_read_eigen

m_bse_io/exc_read_rblock_fio [ Functions ]

[ Top ] [ m_bse_io ] [ Functions ]

NAME

  exc_read_rblock_fio

FUNCTION

  Reads the resonant block from file using Fortran IO.

INPUTS

  funt=Fortran unit number.
  nsppol=Number of spins
  exc_size=Size of the resonant bock.
  diago_is_real=.TRUE. if diagonal elements are real.
  nreh(nsppol)=Number of resonant transitions for each spin.

OUTPUT

  ierr=Status error
  exc_mat(exc_size,exc_size)=The resonant block.

PARENTS

      m_exc_diago

CHILDREN

      c_f_pointer

SOURCE

1271 subroutine exc_read_rblock_fio(funt,diago_is_real,nsppol,nreh,exc_size,exc_mat,ierr)
1272 
1273 
1274 !This section has been created automatically by the script Abilint (TD).
1275 !Do not modify the following lines by hand.
1276 #undef ABI_FUNC
1277 #define ABI_FUNC 'exc_read_rblock_fio'
1278 !End of the abilint section
1279 
1280  implicit none
1281 
1282 !Arguments ------------------------------------
1283 !scalars
1284  integer,intent(in) :: funt,nsppol,exc_size
1285  logical,intent(in) :: diago_is_real
1286  integer,intent(out) :: ierr
1287 !arrays
1288  integer,intent(in) :: nreh(nsppol)
1289  complex(dpc),intent(out) :: exc_mat(exc_size,exc_size)
1290 
1291 !Local variables ------------------------------
1292 !scalars
1293  integer :: itp,it,block,col,row,spad
1294  complex(dpc) :: ctemp
1295  character(len=500) :: errmsg
1296 !arrays
1297  complex(dpc),allocatable :: cbuff_dpc(:)
1298 
1299  ! *************************************************************************
1300 
1301  ierr=0
1302 
1303  ! Construct full resonant block using Hermiticity. File is always in double precision.
1304  ABI_MALLOC(cbuff_dpc,(exc_size))
1305 
1306  if (nsppol==1) then ! Construct resonant block from the upper triangle stored on file.
1307    !
1308    do itp=1,nreh(1)
1309      read(funt, err=10, iomsg=errmsg) cbuff_dpc(1:itp)
1310      do it=1,itp-1 ! Diagonal is treated below.
1311        ctemp = cbuff_dpc(it)
1312        exc_mat(it,itp) = ctemp
1313        exc_mat(itp,it) = CONJG(ctemp)
1314      end do
1315      if (diago_is_real) then
1316        exc_mat(itp,itp) = DBLE(cbuff_dpc(itp))
1317      else
1318        exc_mat(itp,itp) = cbuff_dpc(itp)
1319      end if
1320    end do
1321    !
1322  else
1323    ! The file contains
1324    ! A) The up-up and the down-down block in packed form
1325    ! (only the upper triangle is stored since these blocks are Hermitian)
1326    ! B) The entire up-down exchange block (no symmetry here)
1327    !
1328    ! A) Construct resonant blocks from the upper triangles stored on file.
1329    ! FIXME this part wont work if we have a different number of e-h pairs
1330    !ABI_CHECK(ALL(nreh==nreh(1)),"Different number of transitions")
1331 
1332    do block=1,2
1333      if (block==1) spad=0
1334      if (block==2) spad=nreh(1)
1335      do itp=1,nreh(block)
1336        read(funt, err=10, iomsg=errmsg) cbuff_dpc(1:itp)
1337        col = itp+spad
1338        do it=1,itp-1 ! diagonal is treated below.
1339          ctemp = cbuff_dpc(it)
1340          row = it + spad
1341          exc_mat(row,col) = ctemp
1342          exc_mat(col,row) = CONJG(ctemp)
1343        end do
1344        if (diago_is_real) then
1345          exc_mat(col,col) = DBLE(cbuff_dpc(itp))
1346        else
1347          exc_mat(col,col) = cbuff_dpc(itp)
1348        end if
1349      end do
1350    end do
1351    !
1352    ! read v+- that is a matrix with shape nreh(1) X nreh(2)
1353    spad=nreh(1)
1354    do itp=1,nreh(2)
1355      read(funt, err=10, iomsg=errmsg) cbuff_dpc(1:nreh(1))
1356      exc_mat(1:nreh(1),spad+itp) = cbuff_dpc(1:nreh(1))
1357    end do
1358 
1359  end if
1360 
1361  ABI_FREE(cbuff_dpc)
1362 
1363  return
1364 
1365  ! Raise the error.
1366 10 continue
1367  ierr = 1
1368  MSG_WARNING(errmsg)
1369 
1370 end subroutine exc_read_rblock_fio

m_bse_io/exc_read_rcblock [ Functions ]

[ Top ] [ m_bse_io ] [ Functions ]

NAME

 exc_read_rcblock

FUNCTION

  Reads the excitonic Hamiltonian from file

INPUTS

  fname=File name.
  diago_is_real=.TRUE. if diagonal elements are real (used only if is_resonant==.TRUE.)
  nreh(nsppol)=Number of resonant transition for the two spins.
  is_resonant=Set to .TRUE. if the block is resonant.
  hsize=Dimension of the block.
  nsppol=2 for spin polarized systems. 1 otherwise.
  my_t1,my_t2=The first and the last colums of the matrix treated by this node.
  use_mpio=.TRUE. is MPI-IO routines are used.
  comm=MPI communicator.

OUTPUT

  hmat(hsize,my_t1:my_t2)=The block read from file fname.

TODO

 Remove Bsp

PARENTS

      m_exc_itdiago,m_hexc

CHILDREN

      c_f_pointer

SOURCE

483 subroutine exc_read_rcblock(fname,Bsp,is_resonant,diago_is_real,nsppol,nreh,hsize,my_t1,my_t2,hmat,use_mpio,comm)
484 
485 
486 !This section has been created automatically by the script Abilint (TD).
487 !Do not modify the following lines by hand.
488 #undef ABI_FUNC
489 #define ABI_FUNC 'exc_read_rcblock'
490 !End of the abilint section
491 
492  implicit none
493 
494 !Arguments ------------------------------------
495 !scalars
496  integer,intent(in) :: comm,hsize,my_t1,my_t2,nsppol
497  logical,intent(in) :: is_resonant,use_mpio,diago_is_real
498  character(len=*),intent(in) :: fname
499  type(excparam),intent(in) :: Bsp
500 !arrays
501  integer,intent(in) :: nreh(nsppol)
502  complex(dpc),intent(out) :: hmat(hsize,my_t1:my_t2)
503 
504 !Local variables ------------------------------
505 !scalars
506  integer,parameter :: master=0
507  integer :: it,itp,funit,nproc,my_rank,neh1,neh2
508  integer :: fform,my_nt
509  integer :: row,col,block,spad,spin_dim,ierr,size_exp
510  real(dp) :: cputime,walltime,gflops
511  character(len=500) :: msg,errmsg
512  !type(Hdr_type) :: bse_Hdr
513 !arrays
514  complex(dpc),allocatable :: buffer_dpc(:)
515  logical :: have_row,have_col
516 #ifdef HAVE_MPI_IO
517  integer :: mpierr,mpifh,ham_type,my_nel,old_type,etype,offset_err,amode
518  integer :: irec,nrec !,ncount
519  integer(XMPI_OFFSET_KIND) :: ehdr_offset,my_offset,my_offpad,fsize
520  integer(XMPI_OFFSET_KIND),allocatable :: bsize_frecord(:)
521  integer :: glob_sizes(2),my_cols(2)
522  integer :: block_sizes(2,3)
523  integer :: status(MPI_STATUS_SIZE)
524 #endif
525 
526 !************************************************************************
527 
528  nproc = xmpi_comm_size(comm); my_rank= xmpi_comm_rank(comm)
529 
530  ! Compute the (Expected) size of the hamiltonian.
531  neh1 = nreh(1)
532  neh2 = neh1; if (nsppol==2) neh2=nreh(2)
533 
534  size_exp=neh1; if (nsppol==2) size_exp=SUM(nreh)
535 
536  ABI_CHECK(hsize==size_exp,"Wrong hsize")
537  if (neh1/=neh2) then
538    msg = "BSE code does not support different number of transitions for the two spin channels"
539    MSG_ERROR(msg)
540  end if
541 
542  my_nt = my_t2-my_t1+1
543 
544 !BEGINDEBUG
545 ! hmat = HUGE(zero)
546 !ENDDEBUG
547 
548  if (.not.use_mpio) then
549 
550    if (is_resonant) then
551      call wrtout(std_out,". Reading resonant block from file: "//TRIM(fname)//" using Fortran-IO","COLL")
552    else
553      call wrtout(std_out,". Reading coupling block from file: "//TRIM(fname)//" using Fortran-IO","COLL")
554    end if
555 
556    if (open_file(fname,msg,newunit=funit,form="unformatted",status="old",action="read") /= 0) then
557      MSG_ERROR(msg)
558    end if
559    !
560    ! Read the header and perform consistency checks.
561    call exc_read_bshdr(funit,Bsp,fform,ierr)
562    ABI_CHECK(ierr==0,"Wrong BSE header")
563    !
564    ! Construct full excitonic Hamiltonian using symmetries.
565    if (nsppol==1) then
566      call cwtime(cputime,walltime,gflops,"start")
567      !
568      ABI_MALLOC(buffer_dpc,(neh1))
569      do itp=1,hsize
570        read(funit, err=10, iomsg=errmsg) buffer_dpc(1:itp)
571        !
572        ! Fill the upper triangle if I have this column.
573        if (itp>=my_t1 .and. itp<=my_t2) then
574          do it=1,itp
575            hmat(it,itp) = buffer_dpc(it)
576          end do
577          ! Force the diagonal to be real.
578          if (is_resonant .and.diago_is_real) hmat(itp,itp) = DBLE(hmat(itp,itp))
579        end if
580        !
581        ! Reconstruct the rows below the diagonal (diagonal part is not touched here).
582        if (is_resonant) then ! Use Hermiticity.
583          do it=1,itp-1
584            if (it>=my_t1 .and. it<=my_t2) hmat(itp,it) = CONJG(buffer_dpc(it))
585          end do
586        else  ! Coupling is symmetric.
587          do it=1,itp-1
588            if (it>=my_t1 .and. it<=my_t2) hmat(itp,it) = buffer_dpc(it)
589          end do
590        end if
591        !
592      end do ! itp
593      ABI_FREE(buffer_dpc)
594 
595      call cwtime(cputime,walltime,gflops,"stop")
596      write(msg,'(2(a,f9.1),a)')" Fortran-IO completed. cpu_time ",cputime,"[s], walltime ",walltime," [s]"
597      call wrtout(std_out,msg,"COLL", do_flush=.True.)
598    else
599      ! Spin polarized case.
600      !
601      ! The file contains
602      ! A) The up-up and the down-down block in packed form
603      ! (only the upper triangle is stored since the blocks are Hermitian)
604      ! B) The entire up-down exchange block (no symmetry here)
605      !
606      ! A) Construct resonant blocks from the upper triangles stored on file.
607      ! FIXME this part wont work if we have a different number of e-h pairs
608      if (.not.is_resonant) then
609        MSG_ERROR("exc_read_rcblock does not support coupling.")
610      end if
611      ! It should be checked.
612      spin_dim=neh1
613      do block=1,2
614        ABI_MALLOC(buffer_dpc,(neh1))
615        if (block==1) spad=0
616        if (block==2) spad=neh1
617        do itp=1,spin_dim
618          !
619          read(funit, err=10, iomsg=errmsg) buffer_dpc(1:itp)
620          !
621          ! Fill the upper triangle if this node treats this column
622          col = itp+spad
623          if (col>=my_t1 .and. col<=my_t2) then
624            do it=1,itp
625              row = it + spad
626              hmat(row,col) = buffer_dpc(it)
627            end do
628            ! Force the diagonal to be real.
629            if (is_resonant .and.diago_is_real) hmat(col,col) = DBLE(hmat(col,col))
630          end if
631          !
632          ! Reconstruct the rows below the diagonal (diagonal part is not touched here).
633          row = itp + spad
634          if (is_resonant) then ! Hermitian
635            do it=1,itp-1
636              col = it + spad
637              if (col>=my_t1 .and. col<=my_t2) hmat(row,col) = CONJG(buffer_dpc(it))
638            end do
639          else  ! Coupling is symmetric.
640            do it=1,itp-1
641              col = it + spad
642              if (col>=my_t1 .and. col<=my_t2) hmat(row,col) = buffer_dpc(it)
643            end do
644          end if
645          !
646        end do ! itp
647        ABI_FREE(buffer_dpc)
648      end do ! block
649      !
650      ! B) Kx_{down up} = Kx_{up down}^H.
651      ! FIXME this part wont work if we have a different number of e-h pairs
652      spad=neh1
653      spin_dim=neh1
654      ABI_MALLOC(buffer_dpc,(neh1))
655      do itp=1,spin_dim
656        read(funit, err=10, iomsg=errmsg) buffer_dpc(1:spin_dim)
657        have_col = (spad+itp>=my_t1 .and. spad+itp<=my_t2)
658        if (have_col) hmat(1:spin_dim,spad+itp) = buffer_dpc(1:spin_dim)
659        ! Construct and store the lower block
660        if (is_resonant) then ! Hermitian
661          do it=1,spin_dim
662            have_row = (it>=my_t1 .and. it<=my_t2)
663            if (have_row) hmat(spad+itp,it) = CONJG(buffer_dpc(it))
664          end do
665        else ! Symmetric
666          do it=1,spin_dim
667            have_row = (it>=my_t1 .and. it<=my_t2)
668            if (have_row) hmat(spad+itp,it) = buffer_dpc(it)
669          end do
670        end if
671      end do
672      ABI_FREE(buffer_dpc)
673    end if
674 
675    close(funit)
676 
677  else
678 #ifdef HAVE_MPI_IO
679    if (is_resonant) then
680      call wrtout(std_out,". Reading resonant block from file "//TRIM(fname)//" using MPI-IO","COLL")
681    else
682      call wrtout(std_out,". Reading coupling block from file "//TRIM(fname)//" using MPI-IO","COLL")
683    end if
684 
685    amode=MPI_MODE_RDONLY
686    call MPI_FILE_OPEN(comm,fname,amode,MPI_INFO_NULL,mpifh,mpierr)
687    msg = " FILE_OPEN "//TRIM(fname)
688    ABI_CHECK_MPI(mpierr,msg)
689 
690    call MPI_FILE_GET_SIZE(mpifh,fsize,mpierr)
691    write(std_out,*)" file size is ",fsize
692    !
693    ! Skip the header and find the offset for reading the matrix.
694    call exc_skip_bshdr_mpio(mpifh,xmpio_collective,ehdr_offset)
695    !
696    ! Read my columns from file.
697    old_type=MPI_DOUBLE_COMPLEX
698    glob_sizes = (/hsize,hsize/); my_cols=(/my_t1,my_t2/)
699 
700    if (nsppol==1) then
701      call xmpio_create_coldistr_from_fpacked(glob_sizes,my_cols,old_type,ham_type,my_offpad,offset_err)
702    else
703      MSG_WARNING("nsppol==2 => calling fp3blocks")
704      write(std_out,*)"neh, hsize",neh1,neh2,hsize
705 
706 #if 1
707      nrec=neh1+2*neh2
708      ABI_MALLOC(bsize_frecord,(nrec))
709      bsize_frecord(1:neh1)           = (/(irec*xmpi_bsize_dpc, irec=1,neh1)/)
710      bsize_frecord(neh1+1:neh1+neh2) = (/(irec*xmpi_bsize_dpc, irec=1,neh2)/)
711      bsize_frecord(neh1+neh2+1:)     = neh1*xmpi_bsize_dpc
712      call xmpio_check_frmarkers(mpifh,ehdr_offset,xmpio_collective,nrec,bsize_frecord,ierr)
713      ABI_CHECK(ierr==0,"Error in Fortran markers")
714      ABI_FREE(bsize_frecord)
715      MSG_COMMENT("Marker check ok")
716      call xmpi_barrier(comm)
717 #endif
718 
719      block_sizes(:,1) = (/neh1,neh1/)
720      block_sizes(:,2) = (/neh2,neh2/)
721      block_sizes(:,3) = (/neh1,neh2/)
722      MSG_ERROR("fp3blocks is buggy")
723      call xmpio_create_coldistr_from_fp3blocks(glob_sizes,block_sizes,my_cols,old_type,ham_type,my_offpad,offset_err)
724    end if
725 
726    if (offset_err/=0) then
727      write(msg,"(3a)")&
728 &      "Global position index cannot be stored in a standard Fortran integer ",ch10,&
729 &      "Excitonic matrix cannot be read with a single MPI-IO call."
730      MSG_ERROR(msg)
731    end if
732    !
733    ! The offset used for reading.
734    my_offset = ehdr_offset + my_offpad
735    write(std_out,*)"my_offset= ",my_offset
736 
737    etype=MPI_BYTE
738    call MPI_FILE_SET_VIEW(mpifh, my_offset, etype, ham_type, 'native', MPI_INFO_NULL, mpierr)
739    ABI_CHECK_MPI(mpierr,"SET_VIEW")
740    !
741    ! Release the MPI filetype.
742    call MPI_TYPE_FREE(ham_type,mpierr)
743    ABI_CHECK_MPI(mpierr,"TYPE_FREE")
744 
745    my_nel = my_nt*hsize
746    call MPI_FILE_READ_ALL(mpifh, hmat, my_nel, MPI_DOUBLE_COMPLEX, status, mpierr)
747    ABI_CHECK_MPI(mpierr,"READ_ALL")
748 
749    !call MPI_GET_COUNT(status, MPI_DOUBLE_COMPLEX, ncount, mpierr)
750    !write(std_out,*)"count, my_nel ",ncount,my_nel
751    !ABI_CHECK_MPI(mpierr,"READ_ALL")
752    !
753    ! Close the file.
754    call MPI_FILE_CLOSE(mpifh, mpierr)
755    ABI_CHECK_MPI(mpierr,"FILE_CLOSE")
756    !
757    ! Use the symmetries of the block to reconstruct the local buffer.
758    ! Coupling does not require in-place symmetrization since it is symmetric.
759    if (is_resonant) then
760      do itp=my_t1,my_t2
761        if (itp+1<=hsize) hmat(itp+1:,itp) = DCONJG(hmat(itp+1:,itp)) ! Lower triangle using Hermiticity.
762        if (diago_is_real) hmat(itp,itp) = DBLE(hmat(itp,itp))        ! The diagonal is forced to be real when energies are real.
763      end do
764    end if
765 
766    call xmpi_barrier(comm)
767 #else
768    MSG_ERROR("MPI-IO support not enabled")
769 #endif
770  end if
771 
772 !BEGINDEBUG
773 !  if ( ANY(hmat==HUGE(zero)) ) then
774 !    write(std_out,*)"COUNT",COUNT(hmat==HUGE(zero))," hsize= ",hsize
775 !    MSG_ERROR("Something wrong in the reading")
776 !  end if
777 !ENDDEBUG
778 
779  call xmpi_barrier(comm)
780 
781  return
782 
783 ! Handle IO Error
784 10 continue
785  MSG_ERROR(errmsg)
786 
787 end subroutine exc_read_rcblock

m_bse_io/exc_skip_bshdr [ Functions ]

[ Top ] [ m_bse_io ] [ Functions ]

NAME

  exc_skip_bshdr

FUNCTION

   Skip the header of the (BSR|BSC) files storing the excitonic Hamiltonian. Fortran version.

INPUTS

  funt=Unit number.

OUTPUT

  ierr=Status error.

SIDE EFFECTS

  Skip the header.

PARENTS

      exc_build_block

CHILDREN

      c_f_pointer

SOURCE

242 subroutine exc_skip_bshdr(funt,ierr)
243 
244 
245 !This section has been created automatically by the script Abilint (TD).
246 !Do not modify the following lines by hand.
247 #undef ABI_FUNC
248 #define ABI_FUNC 'exc_skip_bshdr'
249 !End of the abilint section
250 
251  implicit none
252 
253 !Arguments ------------------------------------
254  integer,intent(in) :: funt
255  integer,intent(out) :: ierr
256 
257 !Local variables-------------------------------
258  character(len=500) :: errmsg
259 ! *************************************************************************
260 
261  call hdr_skip(funt,ierr)
262  if (ierr/=0) RETURN
263  read(funt, err=10, iomsg=errmsg)
264 
265  return
266 
267 ! Handle IO Error
268 10 continue
269  ierr = 0
270  MSG_WARNING(errmsg)
271 
272 end subroutine exc_skip_bshdr

m_bse_io/exc_skip_bshdr_mpio [ Functions ]

[ Top ] [ m_bse_io ] [ Functions ]

NAME

  exc_skip_bshdr_mpio

FUNCTION

   Skip the header of the (BSR|BSC) files storing the excitonic Hamiltonian. MPI-IO version.

INPUTS

  mpifh=MPI-IO file handler.
  at_option

SIDE EFFECTS

  ehdr_offset

PARENTS

      exc_build_block,m_bse_io,m_exc_diago

CHILDREN

      c_f_pointer

SOURCE

299 subroutine exc_skip_bshdr_mpio(mpifh,at_option,ehdr_offset)
300 
301 
302 !This section has been created automatically by the script Abilint (TD).
303 !Do not modify the following lines by hand.
304 #undef ABI_FUNC
305 #define ABI_FUNC 'exc_skip_bshdr_mpio'
306 !End of the abilint section
307 
308  implicit none
309 
310  !Arguments ------------------------------------
311  integer,intent(in) :: mpifh,at_option
312  integer(XMPI_OFFSET_KIND),intent(inout) :: ehdr_offset
313 
314 !Local variables ------------------------------
315 !scalars
316  integer :: fform,ierr
317 #ifdef HAVE_MPI_IO
318  integer(XMPI_OFFSET_KIND) :: fmarker
319 #endif
320 
321  ! *************************************************************************
322 
323  call hdr_mpio_skip(mpifh,fform,ehdr_offset)
324 
325 #ifdef HAVE_MPI_IO
326  call xmpio_read_frm(mpifh,ehdr_offset,at_option,fmarker,ierr)
327  !write(std_out,*)"fmarker last record ",fmarker
328 #else
329  MSG_ERROR("You should not be here")
330 #endif
331 
332 end subroutine exc_skip_bshdr_mpio

m_bse_io/exc_write_bshdr [ Functions ]

[ Top ] [ m_bse_io ] [ Functions ]

NAME

  exc_write_bshdr

FUNCTION

   Writes the header of the (BSR|BSC) files storing the excitonic Hamiltonian.

INPUTS

  funt=Fortran unit number.
  Bsp<excparam>=Structure storing the parameters of the run.
  Hdr<hdr_type>=The abinit header.

OUTPUT

  Only writing

PARENTS

      exc_build_block

CHILDREN

      c_f_pointer

SOURCE

104 subroutine exc_write_bshdr(funt,Bsp,Hdr)
105 
106 
107 !This section has been created automatically by the script Abilint (TD).
108 !Do not modify the following lines by hand.
109 #undef ABI_FUNC
110 #define ABI_FUNC 'exc_write_bshdr'
111 !End of the abilint section
112 
113  implicit none
114 
115  !Arguments ------------------------------------
116  integer,intent(in) :: funt
117  type(excparam),intent(in) :: BSp
118  type(hdr_type),intent(inout) :: Hdr
119 
120 !Local variables ------------------------------
121 !scalars
122  integer :: fform_1002 = 1002 ! TODO: change setup_bse so that Hdr_bse reflects the parameters of the run.
123  integer :: ierr
124  character(len=500) :: errmsg
125  ! *************************************************************************
126 
127  call hdr_fort_write(Hdr, funt, fform_1002, ierr)
128  ABI_CHECK(ierr == 0, "hdr_fort_write returned ierr != 0")
129  write(funt, err=10, iomsg=errmsg) BSp%nreh,BSp%nkbz
130 
131  return
132 
133 ! Handle IO Error
134 10 continue
135  MSG_ERROR(errmsg)
136 
137 end subroutine exc_write_bshdr

m_bse_io/exc_write_optme [ Functions ]

[ Top ] [ m_bse_io ] [ Functions ]

NAME

  exc_write_optme

FUNCTION

   Writes the optical matrix elements in the OME.nc file.
   Note that this is only available when NetCDF is available

INPUTS

  filname=filename used to write the optical matrix elements
  minb,maxb=minimum and max band index that have been calculated.
  nkbz=Number of points in the full Brillouin zone.
  nsppol=Number of independent spin polarizations.
  nq=Number of "small" q for optical limit.
  opt_cvk=Optical matrix elements to be written

OUTPUT

  ierr=return status of the writing process.
      => 0 if everything was ok
      => -1 if NetCDF is not available
      => 1 if NetCDF is available

PARENTS

CHILDREN

      c_f_pointer

SOURCE

1552 subroutine exc_write_optme(filname,minb,maxb,nkbz,nsppol,nq,opt_cvk,ierr)
1553 
1554 
1555 !This section has been created automatically by the script Abilint (TD).
1556 !Do not modify the following lines by hand.
1557 #undef ABI_FUNC
1558 #define ABI_FUNC 'exc_write_optme'
1559 !End of the abilint section
1560 
1561  implicit none
1562 
1563  !Arguments ------------------------------------
1564  integer,intent(in) :: minb,maxb,nkbz,nsppol,nq
1565  character(len=fnlen),intent(in) :: filname
1566  complex(dpc),intent(in) :: opt_cvk(minb:maxb,minb:maxb,nkbz,nsppol,nq)
1567  integer,intent(out) :: ierr
1568 
1569 !Local variables ------------------------------
1570 !scalars
1571  integer :: ncid,ncerr,cmplx_id,nband_id,nkbz_id,nsppol_id,nq_id,xyz_id
1572  integer :: minb_id,maxb_id,ome_id,iq,is,ik,ib,jb
1573  integer :: dimOME(6),dimKPT(2),dimQPT(2),dimSCA(0),start6(6),count6(6)
1574  real(dp) :: complex2(2)
1575  ! *************************************************************************
1576 
1577 #ifdef HAVE_NETCDF
1578  ierr = 1
1579 
1580 !1. Create netCDF file
1581  ncerr = nctk_open_create(ncid, filname, xmpi_comm_self)
1582  !MT aug 2013: keep the following on THREE lines in order to help abilint script
1583  NCF_CHECK_MSG(ncerr," create netcdf OME file")
1584 
1585 !2. Define dimensions
1586  NCF_CHECK(nf90_def_dim(ncid,"xyz",3,xyz_id))
1587  NCF_CHECK(nf90_def_dim(ncid,"cmplx",2,cmplx_id))
1588  NCF_CHECK(nf90_def_dim(ncid,"nband",maxb-minb+1,nband_id))
1589  NCF_CHECK(nf90_def_dim(ncid,"nkbz",nkbz,nkbz_id))
1590  NCF_CHECK(nf90_def_dim(ncid,"nq",nq,nq_id))
1591  NCF_CHECK(nf90_def_dim(ncid,"nsppol",nsppol,nsppol_id))
1592 
1593 !Dimensions for optical matrix elements
1594  dimOME = (/ cmplx_id, nband_id, nband_id, nkbz_id, nsppol_id, nq_id /)
1595 !Dimensions for kpoint positions
1596  dimKPT = (/ xyz_id, nkbz_id /)
1597 !Dimensions for qpoints for the optical limit
1598  dimQPT = (/ xyz_id, nq_id /)
1599 
1600 !3. Define variables
1601 
1602  call ab_define_var(ncid, dimOME, ome_id, NF90_DOUBLE,&
1603 & "OME", "Values of optical matrix elements","Tobedone")
1604 
1605  call ab_define_var(ncid, dimSCA, minb_id, NF90_INT,"minb",&
1606 & "Minimum band index for the optical matrix elements", "Dimensionless")
1607 
1608  call ab_define_var(ncid, dimSCA, maxb_id, NF90_INT,"maxb",&
1609 & "Maximum band index for the optical matrix elements", "Dimensionless")
1610 
1611 !4. End define mode
1612  NCF_CHECK(nf90_enddef(ncid))
1613 
1614 !5 Write scalars (minb and maxb)
1615 
1616  NCF_CHECK(nf90_put_var(ncid, minb_id, minb))
1617  NCF_CHECK(nf90_put_var(ncid, maxb_id, maxb))
1618 
1619 !6 Write optical matrix elements
1620 
1621  do iq=1,nq
1622    do is=1,nsppol
1623      do ik=1,nkbz
1624        do ib=minb,maxb
1625          do jb=minb,maxb
1626            start6 = (/ 1, ib-minb+1, jb-minb+1, ik, is, iq /)
1627            count6 = (/ 2, 1, 1, 1, 1, 1 /)
1628            complex2 = (/ REAL(opt_cvk(ib,jb,ik,is,iq)),AIMAG(opt_cvk(ib,jb,ik,is,iq)) /)
1629            ncerr = nf90_put_var(ncid, ome_id, complex2, start = start6, count = count6)
1630            NCF_CHECK_MSG(ncerr," write variable OME")
1631        end do
1632        end do
1633      end do
1634    end do
1635  end do
1636 
1637 !7 Close file
1638 
1639  NCF_CHECK(nf90_close(ncid))
1640  ierr = 0
1641 
1642 #else
1643  ierr = -1
1644 #endif
1645 
1646 end subroutine exc_write_optme

m_bse_io/offset_in_file [ Functions ]

[ Top ] [ m_bse_io ] [ Functions ]

NAME

  offset_in_file

FUNCTION

  Return the offset of the matrix element (row_glob,col_glob)
  size_glob(2) gives the number of row and column of the global matrix
  nsblocks is the number of sublocks, used for nsppol==2 (not used if 1)
  sub_block(2,2,nsblocks)= For each subblock the coordinates of the first and last element.

INPUTS

OUTPUT

PARENTS

SOURCE

1198 function offset_in_file(row_glob,col_glob,size_glob,nsblocks,sub_block,bsize_elm,bsize_frm)
1199 
1200 
1201 !This section has been created automatically by the script Abilint (TD).
1202 !Do not modify the following lines by hand.
1203 #undef ABI_FUNC
1204 #define ABI_FUNC 'offset_in_file'
1205 !End of the abilint section
1206 
1207  implicit none
1208 
1209  !Arguments ------------------------------------
1210  integer(XMPI_OFFSET_KIND) :: offset_in_file
1211  integer,intent(in) :: row_glob,col_glob,nsblocks,bsize_elm,bsize_frm
1212  integer,intent(in) :: size_glob(2),sub_block(2,2,nsblocks)
1213 
1214 !Local variables ------------------------------
1215 !scalars
1216  integer :: ii,jj,ijp_glob,swap
1217  integer(XMPI_OFFSET_KIND) :: my_offset
1218 
1219  ! *************************************************************************
1220 
1221  if (nsblocks==1) then
1222    ii = row_glob
1223    jj = col_glob
1224    if (ii>size_glob(1)/2) ii = ii - size_glob(1)/2 ! Wrap the index.
1225    if (jj>size_glob(2)/2) jj = jj - size_glob(2)/2
1226    if (jj<ii) then ! Exchange the indeces since the symmetric element is read.
1227      swap = jj
1228      jj   = ii
1229      ii   = swap
1230    end if
1231    ijp_glob = ii + jj*(jj-1)/2  ! Index for packed storage mode.
1232    my_offset = (ijp_glob-1)*bsize_elm + (jj-1)*2*bsize_frm
1233  else
1234    ABI_UNUSED(sub_block(1,1,1))
1235    MSG_ERROR("nsppol==2 not coded")
1236  end if
1237 
1238  offset_in_file = my_offset
1239 
1240 end function offset_in_file

m_bse_io/rrs_of_glob [ Functions ]

[ Top ] [ m_bse_io ] [ Functions ]

NAME

  rrs_of_glob

FUNCTION

   [+1,-1,0] if (row_glob,col_glob) belongs to the [ resonant, anti-resonant, (anti)coupling block ]

INPUTS

OUTPUT

PARENTS

SOURCE

1089 pure function rrs_of_glob(row_glob,col_glob,size_glob)
1090 
1091 
1092 !This section has been created automatically by the script Abilint (TD).
1093 !Do not modify the following lines by hand.
1094 #undef ABI_FUNC
1095 #define ABI_FUNC 'rrs_of_glob'
1096 !End of the abilint section
1097 
1098  implicit none
1099 
1100 !Arguments ------------------------------------
1101  integer :: rrs_of_glob
1102  integer,intent(in) :: row_glob,col_glob
1103  integer,intent(in) :: size_glob(2)
1104 
1105 !Local variables ------------------------------
1106 !scalars
1107  integer :: nreh1,nreh2
1108 
1109 ! *************************************************************************
1110 
1111  nreh1=size_glob(1)/2 ! Matrix is square and nreh1==nreh2 but oh well.
1112  nreh2=size_glob(2)/2
1113 
1114  if (row_glob<=nreh1 .and. col_glob<=nreh2) then
1115    rrs_of_glob=+1  ! Resonant.
1116  else if (row_glob >nreh1 .and. col_glob >nreh2) then
1117    rrs_of_glob=-1  ! anti-Resonant.
1118  else
1119    rrs_of_glob=0
1120  end if
1121 
1122 end function rrs_of_glob

m_hexc/exc_ham_ncwrite [ Functions ]

[ Top ] [ m_hexc ] [ Functions ]

NAME

 exc_ham_ncwrite

FUNCTION

  Writes the content of a hexc object to a NETCDF file
  according to the ETSF-IO specifications.

INPUTS

  ncid =NC file handle

OUTPUT

PARENTS

      m_hexc

CHILDREN

      c_f_pointer

SOURCE

1672 subroutine exc_ham_ncwrite(ncid,Kmesh,BSp,hsize,nreh,vcks2t,hreso,diag)
1673 
1674 
1675 !This section has been created automatically by the script Abilint (TD).
1676 !Do not modify the following lines by hand.
1677 #undef ABI_FUNC
1678 #define ABI_FUNC 'exc_ham_ncwrite'
1679 !End of the abilint section
1680 
1681  implicit none
1682 !Arguments ------------------------------------
1683 !scalars
1684  integer,intent(in) :: ncid
1685  integer,intent(in) :: hsize
1686  type(kmesh_t),intent(in) :: Kmesh
1687  type(excparam),intent(in) :: BSp
1688  integer,intent(in) :: nreh(BSp%nsppol)
1689  integer,target,intent(in) :: vcks2t(BSp%maxnbndv,BSp%maxnbndc,Kmesh%nbz,BSp%nsppol)
1690  complex(dpc),target,intent(in) :: hreso(hsize,hsize)
1691  complex(dpc),target,intent(in) :: diag(hsize)
1692 
1693 !Local variables-------------------------------
1694 #ifdef HAVE_NETCDF
1695  integer :: ncerr
1696  integer :: max_nreh, sum_nreh
1697  real(dp), ABI_CONTIGUOUS pointer :: r2vals(:,:),r3vals(:,:,:)
1698 
1699 ! *************************************************************************
1700 
1701  ! ==============================================
1702  ! === Write the dimensions specified by ETSF ===
1703  ! ==============================================
1704  max_nreh = MAXVAL(nreh)
1705  sum_nreh = SUM(nreh)
1706 
1707  ncerr = nctk_def_dims(ncid, [nctkdim_t("number_of_reduced_dimensions", 3), nctkdim_t("number_of_spins", bsp%nsppol),&
1708     nctkdim_t("number_of_kpoints", kmesh%nbz), nctkdim_t("max_number_of_valence_bands", bsp%maxnbndv),&
1709     nctkdim_t("max_number_of_conduction_bands", bsp%maxnbndc), nctkdim_t("max_number_of_transitions", max_nreh),&
1710     nctkdim_t("total_number_of_transitions", sum_nreh), nctkdim_t("cplex", 2)], defmode=.True.)
1711  NCF_CHECK(ncerr)
1712 
1713  ncerr = nctk_def_arrays(ncid, [&
1714    nctkarr_t('vcks2t', "i", 'max_number_of_valence_bands, max_number_of_conduction_bands, number_of_kpoints, number_of_spins'),&
1715    nctkarr_t('hamiltonian', "dp", 'cplex total_number_of_transitions total_number_of_transitions'),&
1716    nctkarr_t('diagonal', "dp", 'cplex total_number_of_transitions'),&
1717    nctkarr_t('lomo', "dp", 'number_of_spins'), &
1718    nctkarr_t('homo', "dp", 'number_of_spins'), &
1719    nctkarr_t('lumo', "dp", 'number_of_spins'), &
1720    nctkarr_t('humo', "dp", 'number_of_spins'), &
1721    nctkarr_t("reduced_coordinates_of_kpoints", "dp", "number_of_reduced_dimensions, number_of_kpoints"), &
1722    nctkarr_t("kpoint_weights", "dp", "number_of_kpoints") &
1723    ])
1724  NCF_CHECK(ncerr)
1725 
1726 ! Write data
1727  NCF_CHECK(nctk_set_datamode(ncid))
1728  NCF_CHECK(nf90_put_var(ncid, vid('vcks2t'), vcks2t))
1729  NCF_CHECK(nf90_put_var(ncid, vid('lomo'), BSp%lomo_spin))
1730  NCF_CHECK(nf90_put_var(ncid, vid('homo'), BSp%homo_spin))
1731  NCF_CHECK(nf90_put_var(ncid, vid('lumo'), BSp%lumo_spin))
1732  NCF_CHECK(nf90_put_var(ncid, vid('humo'), BSp%humo_spin))
1733 
1734  call c_f_pointer(c_loc(hreso(1,1)), r3vals, shape=[2, hsize, hsize])
1735  NCF_CHECK(nf90_put_var(ncid, vid("hamiltonian"), r3vals))
1736 
1737  call c_f_pointer(c_loc(diag(1)), r2vals, shape=[2, hsize])
1738  NCF_CHECK(nf90_put_var(ncid, vid("diag"), r2vals))
1739 
1740  ! Write K-points
1741  NCF_CHECK(nf90_put_var(ncid, vid("reduced_coordinates_of_kpoints"), kmesh%bz))
1742  NCF_CHECK(nf90_put_var(ncid, vid("kpoint_weights"), kmesh%wt))
1743  NCF_CHECK(ncerr)
1744 
1745 #else
1746  MSG_ERROR("netcdf support is not activated. ")
1747 #endif
1748 
1749 contains
1750  integer function vid(vname)
1751 
1752 !This section has been created automatically by the script Abilint (TD).
1753 !Do not modify the following lines by hand.
1754 #undef ABI_FUNC
1755 #define ABI_FUNC 'vid'
1756 !End of the abilint section
1757 
1758    character(len=*),intent(in) :: vname
1759    vid = nctk_idname(ncid, vname)
1760  end function vid
1761 
1762 end subroutine exc_ham_ncwrite