TABLE OF CONTENTS
- ABINIT/m_bse_io
- m_bse_io/ccs_of_glob
- m_bse_io/exc_amplitude
- m_bse_io/exc_fullh_from_blocks
- m_bse_io/exc_read_bshdr
- m_bse_io/exc_read_eigen
- m_bse_io/exc_read_rblock_fio
- m_bse_io/exc_read_rcblock
- m_bse_io/exc_skip_bshdr
- m_bse_io/exc_skip_bshdr_mpio
- m_bse_io/exc_write_bshdr
- m_bse_io/exc_write_optme
- m_bse_io/offset_in_file
- m_bse_io/rrs_of_glob
- m_hexc/exc_ham_ncwrite
ABINIT/m_bse_io [ Modules ]
NAME
m_bse_io
FUNCTION
This module provides routines to read the Bethe-Salpeter Hamiltonian from file
COPYRIGHT
Copyright (C) 2008-2024 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 .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 MODULE m_bse_io 23 24 use defs_basis 25 use m_xmpi 26 use m_errors 27 use m_abicore 28 #if defined HAVE_MPI2 29 use mpi 30 #endif 31 use netcdf 32 use m_nctk 33 use, intrinsic :: iso_c_binding 34 use m_hdr 35 36 use m_time, only : cwtime 37 use m_fstrings, only : toupper 38 use m_io_tools, only : open_file 39 use m_numeric_tools, only : arth 40 use m_special_funcs, only : gaussian 41 use m_bs_defs, only : excparam 42 use m_bz_mesh, only : kmesh_t 43 44 implicit none 45 46 #if defined HAVE_MPI1 47 include 'mpif.h' 48 #endif 49 50 private 51 52 public :: exc_read_rblock_fio ! Reads the entire resonant sub-block from file using Fortran IO. 53 public :: exc_read_rcblock ! Reads a distributed sub-block of the excitonic Hamiltonian from file. 54 public :: exc_fullh_from_blocks ! Initialize the specified sub-blocks of the *full* matrix (reso+anti-reso) from file. 55 public :: rrs_of_glob ! [+1,-1,0] if (row_glob,col_glob) belongs to the [ resonant, anti-resonant, (anti)coupling block ] 56 public :: ccs_of_glob ! [+1,-1,0] if (row_glob,col_glob) belongs to the [ coupling, anti-coupling, (anti)resonant block ] 57 public :: offset_in_file ! Function used to describe the way the Hamiltonian is stored on disk. 58 public :: exc_write_bshdr ! Writes the Header of the (BSR|BSC) files storing the excitonic Hamiltonian. 59 public :: exc_read_bshdr ! Reads the Header of the (BSR|BSC) files. 60 public :: exc_skip_bshdr ! Skip the Header of the (BSR|BSC) files. Fortran version. 61 public :: exc_skip_bshdr_mpio ! Skip the Header of the (BSR|BSC) files. MPI-IO version. 62 public :: exc_read_eigen ! Read selected energies and eigenvectors from the BSEIG file. 63 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. 64 public :: exc_write_optme ! Writes the OME file storing the optical matrix elements 65 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
SOURCE
1015 pure function ccs_of_glob(row_glob,col_glob,size_glob) 1016 1017 !Arguments ------------------------------------ 1018 integer :: ccs_of_glob 1019 integer,intent(in) :: row_glob,col_glob 1020 integer,intent(in) :: size_glob(2) 1021 1022 !Local variables ------------------------------ 1023 !scalars 1024 integer :: nreh1,nreh2 1025 1026 ! ************************************************************************* 1027 1028 nreh1=size_glob(1)/2 ! Matrix is square and nreh1==nreh2 but oh well. 1029 nreh2=size_glob(2)/2 1030 1031 if (row_glob<=nreh1 .and. col_glob >nreh2) then ! Coupling. 1032 ccs_of_glob = +1 1033 else if (row_glob >nreh1 .and. col_glob<=nreh2) then ! anti-Coupling 1034 ccs_of_glob = -1 1035 else 1036 ccs_of_glob = 0 1037 end if 1038 1039 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 indices of the excitonic states sorted in ascending order. out_fname=The name of the file where the results are written.
OUTPUT
Only writing.
SOURCE
1234 subroutine exc_amplitude(Bsp,eig_fname,nvec,vec_idx,out_fname) 1235 1236 !Arguments ------------------------------------ 1237 !scalars 1238 integer,intent(in) :: nvec 1239 character(len=*),intent(in) :: eig_fname,out_fname 1240 type(excparam),intent(in) :: BSp 1241 ! arrays 1242 integer,intent(in) :: vec_idx(nvec) 1243 1244 !Local variables ------------------------------ 1245 !scalars 1246 integer :: vec,art_idx,ierr 1247 integer :: spin,iw,it,nw,pos_w,neg_w,out_unt,rt_idx,hsize 1248 real(dp) :: ene_rt,ampl_eh,ampl_he,xx,stdev,w_max,w_min,step 1249 character(len=500) :: msg 1250 !arrays 1251 real(dp),allocatable :: wmesh(:),amplitude(:),ene_list(:) 1252 complex(dpc),allocatable :: vec_list(:,:) 1253 1254 ! ************************************************************************* 1255 1256 ! Setup of the frequency mesh for F(w). 1257 w_min=greatest_real; w_max=smallest_real 1258 do spin=1,BSp%nsppol 1259 do it=1,BSp%nreh(spin) 1260 ene_rt = Bsp%Trans(it,spin)%en 1261 w_min = MIN(w_min,ene_rt) 1262 w_max = MAX(w_max,ene_rt) 1263 end do 1264 end do 1265 1266 step = Bsp%domega 1267 if (Bsp%use_coupling==0) then 1268 nw = (w_max - w_min)/step + 1 1269 ABI_MALLOC(wmesh,(nw)) 1270 wmesh = arth(w_min,step,nw) 1271 else 1272 ! Both positive and negative frequencies are needed. 1273 pos_w = (w_max - w_min)/step + 1 1274 neg_w = pos_w; if (ABS(w_min) < tol6) neg_w=neg_w-1 ! zero should not included twice. 1275 nw = neg_w + pos_w 1276 ABI_MALLOC(wmesh,(nw)) 1277 wmesh(1:neg_w) = arth(-w_max,step,neg_w) 1278 wmesh(neg_w+1:) = arth( w_min,step,pos_w) 1279 end if 1280 ! 1281 ! Read selected eigenvectors. 1282 hsize = SUM(Bsp%nreh); if (Bsp%use_coupling>0) hsize=2*hsize 1283 1284 ABI_MALLOC(ene_list,(nvec)) 1285 ABI_MALLOC_OR_DIE(vec_list,(hsize,nvec), ierr) 1286 1287 call exc_read_eigen(eig_fname,hsize,nvec,vec_idx,vec_list,ene_list,Bsp=Bsp) 1288 1289 ABI_FREE(ene_list) 1290 1291 if (open_file(out_fname,msg,newunit=out_unt,form="formatted",action="write") /= 0) then 1292 ABI_ERROR(msg) 1293 end if 1294 1295 write(out_unt,*)"# Amplitude functions F(w) = \sum_t |<t|exc_vec>|^2 \delta(ww- ene_t), w is given in eV. " 1296 write(out_unt,*)"# Number of excitonic eigenvectors analyzed: ",nvec 1297 1298 ABI_MALLOC(amplitude,(nw)) 1299 stdev = BSp%broad ! Broadening for the gaussian. 1300 1301 do vec=1,nvec 1302 ! 1303 ! amplitude(ww) = \sum_t |<t|exc_vec>|^2 \delta(ww- ene_t) 1304 amplitude = zero 1305 do spin=1,BSp%nsppol 1306 do it=1,BSp%nreh(spin) 1307 ene_rt = Bsp%Trans(it,spin)%en 1308 rt_idx = it + (spin-1)*Bsp%nreh(1) 1309 ampl_eh = (ABS(vec_list(rt_idx,vec)))**2 1310 if (Bsp%use_coupling>0) then ! Need the index and the amplitude of the anti-resonant component. 1311 art_idx = it + (spin-1)*Bsp%nreh(1) + SUM(Bsp%nreh) 1312 ampl_he = (ABS(vec_list(art_idx,vec)))**2 1313 end if 1314 ! 1315 do iw=1,nw ! Accumulate 1316 xx = wmesh(iw) - ene_rt 1317 amplitude(iw) = amplitude(iw) + ampl_eh * gaussian(xx, stdev) 1318 if (Bsp%use_coupling>0) then 1319 xx = wmesh(iw) + ene_rt 1320 amplitude(iw) = amplitude(iw) + ampl_he * gaussian(xx, stdev) 1321 end if 1322 end do 1323 ! 1324 end do 1325 end do 1326 ! 1327 ! Write results 1328 write(out_unt,*)"# Amplitude function F(w) for exc_vector index ",vec_idx(vec) 1329 do iw=1,nw 1330 write(out_unt,*)wmesh(iw)*Ha_eV,amplitude(iw) 1331 end do 1332 write(out_unt,*)"#" 1333 end do 1334 1335 close(out_unt) 1336 1337 ABI_FREE(amplitude) 1338 ABI_FREE(wmesh) 1339 ABI_FREE(vec_list) 1340 1341 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
SOURCE
725 subroutine exc_fullh_from_blocks(funt,block_type,nsppol,row_sign,diago_is_real,nreh,exc_size,exc_ham) 726 727 !Arguments ------------------------------------ 728 !scalars 729 integer,intent(in) :: funt,exc_size,nsppol,row_sign 730 logical,intent(in) :: diago_is_real 731 character(len=*),intent(in) :: block_type 732 !arrays 733 integer,intent(in) :: nreh(nsppol) 734 complex(dpc),intent(inout) :: exc_ham(exc_size,exc_size) 735 736 !Local variables------------------------------- 737 !scalars 738 integer :: it,itp,szbuf,neh,pad_c1,pad_r1,spin_dim,spad_r,spad_c 739 integer :: block,spad,row1,col1,row2,col2,spin_stride,ierr 740 complex(dpc) :: cttp 741 character(len=500) :: errmsg 742 !arrays 743 complex(dpc),allocatable :: cbuff_dpc(:) 744 745 ! ********************************************************************* 746 747 szbuf=exc_size ! FIXME oversized! 748 neh = nreh(1) 749 750 if (nsppol==2) then 751 ABI_WARNING("nsppol==2 is very experimental") 752 end if 753 if (ANY(nreh(1)/=nreh)) then 754 ABI_ERROR(" different nreh are not supported") 755 end if 756 757 ABI_MALLOC_OR_DIE(cbuff_dpc,(exc_size), ierr) 758 ! 759 ! The two cases nsppol==1,2 can be merged but we keep them 760 ! separated to keep to code readable. 761 762 SELECT CASE (toupper(block_type)) 763 CASE ("RESONANT") 764 765 if (nsppol==1) then 766 ! 767 ! Construct resonant block from the upper triangle stored on file. 768 neh = nreh(1) 769 do itp=1,neh 770 read(funt, err=10, iomsg=errmsg) cbuff_dpc(1:itp) 771 do it=1,itp-1 ! The diagonal is treated below. 772 cttp = cbuff_dpc(it) 773 exc_ham(it ,itp) = cttp ! R 774 exc_ham(itp ,it ) = CONJG(cttp) ! row_sign R* 775 exc_ham(neh+it ,neh+itp) = row_sign*CONJG(cttp) 776 exc_ham(neh+itp,neh+it ) = row_sign*cttp 777 end do 778 if (diago_is_real) then 779 exc_ham(itp ,itp) = DBLE(cbuff_dpc(itp)) 780 exc_ham(neh+itp,neh+itp) = row_sign*DBLE(cbuff_dpc(itp)) 781 else 782 exc_ham(itp,itp) = cbuff_dpc(itp) 783 exc_ham(neh+itp,neh+itp) = row_sign*(CONJG(cbuff_dpc(itp))) 784 end if 785 end do 786 ! 787 ! 788 else 789 ! FIXME this part wont work if we have a different number of e-h pairs 790 ABI_CHECK(ALL(nreh==nreh(1)),"Different number of transitions") 791 ! The file contains 792 ! A) The up-up and the down-down block in packed form 793 ! (only the upper triangle is stored since these blocks are Hermitian) 794 ! B) The entire up-down exchange block (no symmetry here) 795 ! 796 ! The resonant block is given by: 797 ! | (v'c' up) | (v'c' dwn) | 798 ! ------------------------------ where v_{-+} = v_{+-}^H when the momentum of the photon is neglected. 799 ! | [T-W+v]++ | v+- | (vc up) Note that v_{+-} is not Hermitian due to the presence of different spins. 800 ! R = ------------------------------ Actually it reduces to a Hermitian matrix when the system is not spin polarized. 801 ! | v-+ | [T-W+v]-- | (vc dwn) [T-W+v] is Hermitian provided the the QP energies are purely real. 802 ! ------------------------------ 803 ! 804 ! *) Fill the diagonal blocks. 805 ! only the upper triangle is stored on file. 806 ! row1,col1 refer to the resonant block. 807 ! row2,col2 refer to the anti-resonant block. 808 do block=1,2 809 if (block==1) then 810 spad=0 811 spin_stride=SUM(nreh) 812 end if 813 if (block==2) then 814 spad=nreh(1) 815 spin_stride=2*nreh(1) + nreh(2) 816 end if 817 do itp=1,nreh(block) 818 read(funt, err=10, iomsg=errmsg) cbuff_dpc(1:itp) 819 col1 = itp+spad 820 col2 = itp+spin_stride 821 do it=1,itp-1 822 cttp = cbuff_dpc(it) 823 row1 = it + spad 824 row2 = it + spin_stride 825 exc_ham(row1,col1) = cttp ! [T-W+v] 826 exc_ham(col1,row1) = CONJG(cttp) 827 exc_ham(row2,col2) = row_sign*CONJG(cttp) ! row_sign [T-W+v]* 828 exc_ham(col2,row2) = row_sign*cttp 829 end do 830 if (diago_is_real) then 831 exc_ham(col1,col1) = DBLE(cbuff_dpc(itp)) 832 exc_ham(col2,col2) = row_sign*DBLE(cbuff_dpc(itp)) 833 else 834 exc_ham(col1,col1) = cbuff_dpc(itp) 835 exc_ham(col2,col2) = row_sign*CONJG(cbuff_dpc(itp)) 836 end if 837 end do 838 end do 839 ! 840 ! Read v+- and reconstruct resonant and anti-resonat blocks. 841 ! TODO recheck this 842 pad_r1=SUM(nreh) 843 pad_c1=2*nreh(1) + nreh(2) 844 845 spin_dim=nreh(1) 846 do itp=1,spin_dim 847 read(funt, err=10, iomsg=errmsg) cbuff_dpc(1:spin_dim) 848 exc_ham(1:spin_dim,nreh(1)+itp) = cbuff_dpc(1:spin_dim) ! upper reso 849 exc_ham(1+pad_r1:pad_r1+spin_dim,pad_c1+itp) = row_sign*CONJG(cbuff_dpc(1:spin_dim)) ! upper anti-reso 850 col1 = itp+nreh(1) 851 col2 = itp+(2*nreh(1) + nreh(2)) 852 do it=1,spin_dim 853 cttp = cbuff_dpc(it) 854 row1 = it 855 exc_ham(col1,row1) = CONJG(cttp) ! lower reso. 856 row2 = it + SUM(nreh) 857 exc_ham(col2,row2) = row_sign*cttp ! lower anti-reso. 858 end do 859 end do 860 end if 861 862 CASE ("COUPLING") 863 ! 864 if (nsppol==1) then 865 ! 866 do itp=1,neh 867 read(funt, err=10, iomsg=errmsg) cbuff_dpc(1:itp) 868 do it=1,itp 869 cttp = cbuff_dpc(it) 870 exc_ham(it ,neh+itp) = cttp 871 exc_ham(itp ,neh+it ) = cttp 872 exc_ham(neh+it ,itp ) = row_sign*CONJG(cttp) 873 exc_ham(neh+itp,it ) = row_sign*CONJG(cttp) 874 end do 875 end do 876 ! 877 else 878 ! The coupling block is given by: 879 ! | (c'v' up) | (c'v dwn) | 880 ! ----------------------------------- where v_{-+} = v_{+-}^t when the momentum of the photon is neglected. 881 ! | [-W+v]++ | v+- | (vc up) The entire matrix v_{+-} is stored on file. 882 ! C = ----------------------------------- 883 ! | v-+ | [-W+v]-- | (vc dwn) 884 ! ----------------------------------- 885 ! 886 ! *) Fill blocks that are diagonal in spin coordinates. 887 ! row1,col1 refer to the resonat block. 888 ! row2,col2 refer to the anti-resonant block. 889 do block=1,2 890 if (block==1) then 891 spad_r=0 892 spad_c=SUM(nreh) 893 end if 894 if (block==2) then 895 spad_r=nreh(1) 896 spad_c=2*nreh(1)+nreh(2) 897 end if 898 do itp=1,nreh(block) 899 read(funt, err=10, iomsg=errmsg) cbuff_dpc(1:itp) 900 col1 = itp+spad_c 901 row2 = itp+spad_r 902 do it=1,itp-1 903 cttp = cbuff_dpc(it) 904 row1 = it + spad_r 905 col2 = it + spad_c 906 exc_ham(row1,col1) = cttp ! upper coupling 907 exc_ham(row2,col2) = cttp ! lower coupling (symmetric) 908 exc_ham(col1,row1) = row_sign*CONJG(cttp) ! lower anti-coupling 909 exc_ham(col2,row2) = row_sign*CONJG(cttp) ! upper anti-couling 910 end do 911 ! TODO recheck this 912 row1 = itp+spad_r 913 exc_ham(row1,col1) = cbuff_dpc(itp) ! Diagonals of the block. 914 col2 = itp+spad_c 915 exc_ham(col2,row2) = row_sign*CONJG(cbuff_dpc(itp)) 916 end do 917 end do 918 ! 919 ! Read Full v+- and reconstruct resonant and anti-resonat blocks. 920 ! TODO recheck this 921 spad=2*nreh(1) + nreh(2) 922 pad_r1=SUM(nreh) 923 pad_c1=2*nreh(1) + nreh(2) 924 925 spin_dim=nreh(1) 926 do itp=1,spin_dim 927 read(funt, err=10, iomsg=errmsg) cbuff_dpc(1:spin_dim) 928 exc_ham(1:spin_dim,spad+itp) = cbuff_dpc(1:spin_dim) ! upper block reso 929 exc_ham(1+pad_r1:pad_r1+spin_dim,nreh(1)+itp) = row_sign*CONJG(cbuff_dpc(1:spin_dim)) ! upper block anti-reso 930 col1 = itp+spad 931 row2 = itp+nreh(1) 932 do it=1,spin_dim 933 cttp = cbuff_dpc(it) 934 row1 = it 935 exc_ham(col1,row1) = row_sign*CONJG(cttp) ! lower anti-reso. 936 col2 = it + SUM(nreh) 937 exc_ham(row2,col2) = cttp ! lower reso. 938 end do 939 end do 940 ! 941 end if 942 943 CASE DEFAULT 944 ABI_ERROR("Unknown block_type: "//TRIM(block_type)) 945 END SELECT 946 947 ABI_FREE(cbuff_dpc) 948 949 return 950 951 ! Handle IO Error 952 10 continue 953 ABI_ERROR(errmsg) 954 955 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.
SOURCE
138 subroutine exc_read_bshdr(funt,Bsp,fform,ierr) 139 140 !Arguments ------------------------------------ 141 integer,intent(in) :: funt 142 integer,intent(out) :: fform,ierr 143 type(excparam),intent(in) :: BSp 144 145 !Local variables ------------------------------ 146 !scalars 147 integer :: nkbz_read 148 character(len=500) :: errmsg 149 type(hdr_type) :: Hdr 150 !arrays 151 integer :: nreh_read(SIZE(BSp%nreh)) 152 153 ! ************************************************************************* 154 155 ierr=0 156 157 ! Read the header and perform consistency checks. 158 call hdr_fort_read(hdr, funt, fform, rewind=.True.) 159 ABI_CHECK(fform /= 0, "hdr_fort_read returned fform == 0") 160 161 read(funt, err=10, iomsg=errmsg) nreh_read, nkbz_read 162 163 call Hdr%free() 164 165 if (ANY(nreh_read/=BSp%nreh)) then 166 call wrtout(std_out,"Wrong number of e-h transitions","COLL") 167 ierr = ierr + 1 168 end if 169 170 return 171 172 10 ierr = 1 173 ABI_WARNING(errmsg) 174 175 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 indices 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.
SOURCE
291 subroutine exc_read_eigen(eig_fname,hsize,nvec,vec_idx,vec_list,ene_list,Bsp) 292 293 !Arguments ------------------------------------ 294 !scalars 295 integer,intent(in) :: nvec,hsize 296 character(len=*),intent(in) :: eig_fname 297 type(excparam),optional,intent(in) :: BSp 298 ! arrays 299 integer,intent(in) :: vec_idx(nvec) 300 real(dp),optional,intent(out) :: ene_list(nvec) 301 complex(dpc),intent(out) :: vec_list(hsize,nvec) 302 303 !Local variables ------------------------------ 304 !scalars 305 integer :: eig_unt,hsize_read,neig_read,ll,vec 306 character(len=500) :: msg,errmsg 307 !arrays 308 !real(dp),allocatable :: exc_ene(:) 309 complex(dpc),allocatable :: exc_ene_cplx(:) 310 311 ! ************************************************************************* 312 313 ABI_UNUSED(BSp%nline) 314 315 if (open_file(eig_fname,msg,newunit=eig_unt,form="unformatted",status="old",action="read") /= 0) then 316 ABI_ERROR(msg) 317 end if 318 319 read(eig_unt, err=10, iomsg=errmsg)hsize_read,neig_read 320 !write(std_out,*)hsize_read, neig_read 321 322 if (hsize_read/=hsize) then 323 write(msg,'(a,2(1x,i0))')" hsize_read/=hsize: ",hsize_read,hsize 324 ABI_ERROR(msg) 325 end if 326 327 ! Read eigenvalues, ignore possibly small imaginary part. 328 ABI_MALLOC(exc_ene_cplx,(neig_read)) 329 read(eig_unt, err=10, iomsg=errmsg) exc_ene_cplx 330 331 if (PRESENT(ene_list)) then 332 do vec=1,nvec 333 ll = vec_idx(vec) 334 ene_list(vec) = DBLE(exc_ene_cplx(ll)) 335 end do 336 end if 337 ABI_FREE(exc_ene_cplx) 338 339 vec=1 340 do ll=1,neig_read ! Read the selected excitons. 341 if (ll==vec_idx(vec)) then 342 read(eig_unt, err=10, iomsg=errmsg) vec_list(:,vec) 343 if (vec==nvec) EXIT 344 vec=vec+1 345 else 346 read(eig_unt, err=10, iomsg=errmsg) 347 end if 348 end do 349 350 close(eig_unt, err=10, iomsg=errmsg) 351 352 if (vec/=nvec) then 353 write(msg,'(a,2(1x,i0))')" vec_idx is wrong, vec/=nvec ",vec,nvec+1 354 ABI_ERROR(msg) 355 end if 356 357 return 358 359 ! Handle IO-error 360 10 continue 361 ABI_ERROR(errmsg) 362 363 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.
SOURCE
1118 subroutine exc_read_rblock_fio(funt,diago_is_real,nsppol,nreh,exc_size,exc_mat,ierr) 1119 1120 !Arguments ------------------------------------ 1121 !scalars 1122 integer,intent(in) :: funt,nsppol,exc_size 1123 logical,intent(in) :: diago_is_real 1124 integer,intent(out) :: ierr 1125 !arrays 1126 integer,intent(in) :: nreh(nsppol) 1127 complex(dpc),intent(out) :: exc_mat(exc_size,exc_size) 1128 1129 !Local variables ------------------------------ 1130 !scalars 1131 integer :: itp,it,block,col,row,spad 1132 complex(dpc) :: ctemp 1133 character(len=500) :: errmsg 1134 !arrays 1135 complex(dpc),allocatable :: cbuff_dpc(:) 1136 1137 ! ************************************************************************* 1138 1139 ierr=0 1140 1141 ! Construct full resonant block using Hermiticity. File is always in double precision. 1142 ABI_MALLOC(cbuff_dpc,(exc_size)) 1143 1144 if (nsppol==1) then ! Construct resonant block from the upper triangle stored on file. 1145 ! 1146 do itp=1,nreh(1) 1147 read(funt, err=10, iomsg=errmsg) cbuff_dpc(1:itp) 1148 do it=1,itp-1 ! Diagonal is treated below. 1149 ctemp = cbuff_dpc(it) 1150 exc_mat(it,itp) = ctemp 1151 exc_mat(itp,it) = CONJG(ctemp) 1152 end do 1153 if (diago_is_real) then 1154 exc_mat(itp,itp) = DBLE(cbuff_dpc(itp)) 1155 else 1156 exc_mat(itp,itp) = cbuff_dpc(itp) 1157 end if 1158 end do 1159 ! 1160 else 1161 ! The file contains 1162 ! A) The up-up and the down-down block in packed form 1163 ! (only the upper triangle is stored since these blocks are Hermitian) 1164 ! B) The entire up-down exchange block (no symmetry here) 1165 ! 1166 ! A) Construct resonant blocks from the upper triangles stored on file. 1167 ! FIXME this part wont work if we have a different number of e-h pairs 1168 !ABI_CHECK(ALL(nreh==nreh(1)),"Different number of transitions") 1169 1170 do block=1,2 1171 if (block==1) spad=0 1172 if (block==2) spad=nreh(1) 1173 do itp=1,nreh(block) 1174 read(funt, err=10, iomsg=errmsg) cbuff_dpc(1:itp) 1175 col = itp+spad 1176 do it=1,itp-1 ! diagonal is treated below. 1177 ctemp = cbuff_dpc(it) 1178 row = it + spad 1179 exc_mat(row,col) = ctemp 1180 exc_mat(col,row) = CONJG(ctemp) 1181 end do 1182 if (diago_is_real) then 1183 exc_mat(col,col) = DBLE(cbuff_dpc(itp)) 1184 else 1185 exc_mat(col,col) = cbuff_dpc(itp) 1186 end if 1187 end do 1188 end do 1189 ! 1190 ! read v+- that is a matrix with shape nreh(1) X nreh(2) 1191 spad=nreh(1) 1192 do itp=1,nreh(2) 1193 read(funt, err=10, iomsg=errmsg) cbuff_dpc(1:nreh(1)) 1194 exc_mat(1:nreh(1),spad+itp) = cbuff_dpc(1:nreh(1)) 1195 end do 1196 1197 end if 1198 1199 ABI_FREE(cbuff_dpc) 1200 1201 return 1202 1203 ! Raise the error. 1204 10 continue 1205 ierr = 1 1206 ABI_WARNING(errmsg) 1207 1208 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
SOURCE
394 subroutine exc_read_rcblock(fname,Bsp,is_resonant,diago_is_real,nsppol,nreh,hsize,my_t1,my_t2,hmat,use_mpio,comm) 395 396 !Arguments ------------------------------------ 397 !scalars 398 integer,intent(in) :: comm,hsize,my_t1,my_t2,nsppol 399 logical,intent(in) :: is_resonant,use_mpio,diago_is_real 400 character(len=*),intent(in) :: fname 401 type(excparam),intent(in) :: Bsp 402 !arrays 403 integer,intent(in) :: nreh(nsppol) 404 complex(dpc),intent(out) :: hmat(hsize,my_t1:my_t2) 405 406 !Local variables ------------------------------ 407 !scalars 408 integer,parameter :: master=0 409 integer :: it,itp,funit,nproc,my_rank,neh1,neh2 410 integer :: fform,my_nt 411 integer :: row,col,block,spad,spin_dim,ierr,size_exp 412 real(dp) :: cputime,walltime,gflops 413 character(len=500) :: msg,errmsg 414 !type(Hdr_type) :: bse_Hdr 415 !arrays 416 complex(dpc),allocatable :: buffer_dpc(:) 417 logical :: have_row,have_col 418 #ifdef HAVE_MPI_IO 419 integer :: mpierr,mpifh,ham_type,my_nel,old_type,etype,offset_err,amode 420 integer :: irec,nrec !,ncount 421 integer(XMPI_OFFSET_KIND) :: ehdr_offset,my_offset,my_offpad,fsize 422 integer(XMPI_OFFSET_KIND),allocatable :: bsize_frecord(:) 423 integer :: glob_sizes(2),my_cols(2) 424 integer :: block_sizes(2,3) 425 integer :: status(MPI_STATUS_SIZE) 426 #endif 427 428 !************************************************************************ 429 430 nproc = xmpi_comm_size(comm); my_rank= xmpi_comm_rank(comm) 431 432 ! Compute the (Expected) size of the hamiltonian. 433 neh1 = nreh(1) 434 neh2 = neh1; if (nsppol==2) neh2=nreh(2) 435 436 size_exp=neh1; if (nsppol==2) size_exp=SUM(nreh) 437 438 ABI_CHECK(hsize==size_exp,"Wrong hsize") 439 if (neh1/=neh2) then 440 msg = "BSE code does not support different number of transitions for the two spin channels" 441 ABI_ERROR(msg) 442 end if 443 444 my_nt = my_t2-my_t1+1 445 446 !BEGINDEBUG 447 ! hmat = HUGE(zero) 448 !ENDDEBUG 449 450 if (.not.use_mpio) then 451 452 if (is_resonant) then 453 call wrtout(std_out,". Reading resonant block from file: "//TRIM(fname)//" using Fortran-IO","COLL") 454 else 455 call wrtout(std_out,". Reading coupling block from file: "//TRIM(fname)//" using Fortran-IO","COLL") 456 end if 457 458 if (open_file(fname,msg,newunit=funit,form="unformatted",status="old",action="read") /= 0) then 459 ABI_ERROR(msg) 460 end if 461 ! 462 ! Read the header and perform consistency checks. 463 call exc_read_bshdr(funit,Bsp,fform,ierr) 464 ABI_CHECK(ierr==0,"Wrong BSE header") 465 ! 466 ! Construct full excitonic Hamiltonian using symmetries. 467 if (nsppol==1) then 468 call cwtime(cputime,walltime,gflops,"start") 469 ! 470 ABI_MALLOC(buffer_dpc,(neh1)) 471 do itp=1,hsize 472 read(funit, err=10, iomsg=errmsg) buffer_dpc(1:itp) 473 ! 474 ! Fill the upper triangle if I have this column. 475 if (itp>=my_t1 .and. itp<=my_t2) then 476 do it=1,itp 477 hmat(it,itp) = buffer_dpc(it) 478 end do 479 ! Force the diagonal to be real. 480 if (is_resonant .and.diago_is_real) hmat(itp,itp) = DBLE(hmat(itp,itp)) 481 end if 482 ! 483 ! Reconstruct the rows below the diagonal (diagonal part is not touched here). 484 if (is_resonant) then ! Use Hermiticity. 485 do it=1,itp-1 486 if (it>=my_t1 .and. it<=my_t2) hmat(itp,it) = CONJG(buffer_dpc(it)) 487 end do 488 else ! Coupling is symmetric. 489 do it=1,itp-1 490 if (it>=my_t1 .and. it<=my_t2) hmat(itp,it) = buffer_dpc(it) 491 end do 492 end if 493 ! 494 end do ! itp 495 ABI_FREE(buffer_dpc) 496 497 call cwtime(cputime,walltime,gflops,"stop") 498 write(msg,'(2(a,f9.1),a)')" Fortran-IO completed. cpu_time ",cputime,"[s], walltime ",walltime," [s]" 499 call wrtout(std_out,msg,"COLL", do_flush=.True.) 500 else 501 ! Spin polarized case. 502 ! 503 ! The file contains 504 ! A) The up-up and the down-down block in packed form 505 ! (only the upper triangle is stored since the blocks are Hermitian) 506 ! B) The entire up-down exchange block (no symmetry here) 507 ! 508 ! A) Construct resonant blocks from the upper triangles stored on file. 509 ! FIXME this part wont work if we have a different number of e-h pairs 510 if (.not.is_resonant) then 511 ABI_ERROR("exc_read_rcblock does not support coupling.") 512 end if 513 ! It should be checked. 514 spin_dim=neh1 515 do block=1,2 516 ABI_MALLOC(buffer_dpc,(neh1)) 517 if (block==1) spad=0 518 if (block==2) spad=neh1 519 do itp=1,spin_dim 520 ! 521 read(funit, err=10, iomsg=errmsg) buffer_dpc(1:itp) 522 ! 523 ! Fill the upper triangle if this node treats this column 524 col = itp+spad 525 if (col>=my_t1 .and. col<=my_t2) then 526 do it=1,itp 527 row = it + spad 528 hmat(row,col) = buffer_dpc(it) 529 end do 530 ! Force the diagonal to be real. 531 if (is_resonant .and.diago_is_real) hmat(col,col) = DBLE(hmat(col,col)) 532 end if 533 ! 534 ! Reconstruct the rows below the diagonal (diagonal part is not touched here). 535 row = itp + spad 536 if (is_resonant) then ! Hermitian 537 do it=1,itp-1 538 col = it + spad 539 if (col>=my_t1 .and. col<=my_t2) hmat(row,col) = CONJG(buffer_dpc(it)) 540 end do 541 else ! Coupling is symmetric. 542 do it=1,itp-1 543 col = it + spad 544 if (col>=my_t1 .and. col<=my_t2) hmat(row,col) = buffer_dpc(it) 545 end do 546 end if 547 ! 548 end do ! itp 549 ABI_FREE(buffer_dpc) 550 end do ! block 551 ! 552 ! B) Kx_{down up} = Kx_{up down}^H. 553 ! FIXME this part wont work if we have a different number of e-h pairs 554 spad=neh1 555 spin_dim=neh1 556 ABI_MALLOC(buffer_dpc,(neh1)) 557 do itp=1,spin_dim 558 read(funit, err=10, iomsg=errmsg) buffer_dpc(1:spin_dim) 559 have_col = (spad+itp>=my_t1 .and. spad+itp<=my_t2) 560 if (have_col) hmat(1:spin_dim,spad+itp) = buffer_dpc(1:spin_dim) 561 ! Construct and store the lower block 562 if (is_resonant) then ! Hermitian 563 do it=1,spin_dim 564 have_row = (it>=my_t1 .and. it<=my_t2) 565 if (have_row) hmat(spad+itp,it) = CONJG(buffer_dpc(it)) 566 end do 567 else ! Symmetric 568 do it=1,spin_dim 569 have_row = (it>=my_t1 .and. it<=my_t2) 570 if (have_row) hmat(spad+itp,it) = buffer_dpc(it) 571 end do 572 end if 573 end do 574 ABI_FREE(buffer_dpc) 575 end if 576 577 close(funit) 578 579 else 580 #ifdef HAVE_MPI_IO 581 if (is_resonant) then 582 call wrtout(std_out,". Reading resonant block from file "//TRIM(fname)//" using MPI-IO","COLL") 583 else 584 call wrtout(std_out,". Reading coupling block from file "//TRIM(fname)//" using MPI-IO","COLL") 585 end if 586 587 amode=MPI_MODE_RDONLY 588 call MPI_FILE_OPEN(comm,fname,amode,MPI_INFO_NULL,mpifh,mpierr) 589 msg = " FILE_OPEN "//TRIM(fname) 590 ABI_CHECK_MPI(mpierr,msg) 591 592 call MPI_FILE_GET_SIZE(mpifh,fsize,mpierr) 593 write(std_out,*)" file size is ",fsize 594 ! 595 ! Skip the header and find the offset for reading the matrix. 596 call exc_skip_bshdr_mpio(mpifh,xmpio_collective,ehdr_offset) 597 ! 598 ! Read my columns from file. 599 old_type=MPI_DOUBLE_COMPLEX 600 glob_sizes = (/hsize,hsize/); my_cols=(/my_t1,my_t2/) 601 602 if (nsppol==1) then 603 call xmpio_create_coldistr_from_fpacked(glob_sizes,my_cols,old_type,ham_type,my_offpad,offset_err) 604 else 605 ABI_WARNING("nsppol==2 => calling fp3blocks") 606 write(std_out,*)"neh, hsize",neh1,neh2,hsize 607 608 #if 1 609 nrec=neh1+2*neh2 610 ABI_MALLOC(bsize_frecord,(nrec)) 611 bsize_frecord(1:neh1) = (/(irec*xmpi_bsize_dpc, irec=1,neh1)/) 612 bsize_frecord(neh1+1:neh1+neh2) = (/(irec*xmpi_bsize_dpc, irec=1,neh2)/) 613 bsize_frecord(neh1+neh2+1:) = neh1*xmpi_bsize_dpc 614 call xmpio_check_frmarkers(mpifh,ehdr_offset,xmpio_collective,nrec,bsize_frecord,ierr) 615 ABI_CHECK(ierr==0,"Error in Fortran markers") 616 ABI_FREE(bsize_frecord) 617 ABI_COMMENT("Marker check ok") 618 call xmpi_barrier(comm) 619 #endif 620 621 block_sizes(:,1) = (/neh1,neh1/) 622 block_sizes(:,2) = (/neh2,neh2/) 623 block_sizes(:,3) = (/neh1,neh2/) 624 ABI_ERROR("fp3blocks is buggy") 625 call xmpio_create_coldistr_from_fp3blocks(glob_sizes,block_sizes,my_cols,old_type,ham_type,my_offpad,offset_err) 626 end if 627 628 if (offset_err/=0) then 629 write(msg,"(3a)")& 630 & "Global position index cannot be stored in a standard Fortran integer ",ch10,& 631 & "Excitonic matrix cannot be read with a single MPI-IO call." 632 ABI_ERROR(msg) 633 end if 634 ! 635 ! The offset used for reading. 636 my_offset = ehdr_offset + my_offpad 637 write(std_out,*)"my_offset= ",my_offset 638 639 etype=MPI_BYTE 640 call MPI_FILE_SET_VIEW(mpifh, my_offset, etype, ham_type, 'native', MPI_INFO_NULL, mpierr) 641 ABI_CHECK_MPI(mpierr,"SET_VIEW") 642 ! 643 ! Release the MPI filetype. 644 call MPI_TYPE_FREE(ham_type,mpierr) 645 ABI_CHECK_MPI(mpierr,"TYPE_FREE") 646 647 my_nel = my_nt*hsize 648 call MPI_FILE_READ_ALL(mpifh, hmat, my_nel, MPI_DOUBLE_COMPLEX, status, mpierr) 649 ABI_CHECK_MPI(mpierr,"READ_ALL") 650 651 !call MPI_GET_COUNT(status, MPI_DOUBLE_COMPLEX, ncount, mpierr) 652 !write(std_out,*)"count, my_nel ",ncount,my_nel 653 !ABI_CHECK_MPI(mpierr,"READ_ALL") 654 ! 655 ! Close the file. 656 call MPI_FILE_CLOSE(mpifh, mpierr) 657 ABI_CHECK_MPI(mpierr,"FILE_CLOSE") 658 ! 659 ! Use the symmetries of the block to reconstruct the local buffer. 660 ! Coupling does not require in-place symmetrization since it is symmetric. 661 if (is_resonant) then 662 do itp=my_t1,my_t2 663 if (itp+1<=hsize) hmat(itp+1:,itp) = DCONJG(hmat(itp+1:,itp)) ! Lower triangle using Hermiticity. 664 if (diago_is_real) hmat(itp,itp) = DBLE(hmat(itp,itp)) ! The diagonal is forced to be real when energies are real. 665 end do 666 end if 667 668 call xmpi_barrier(comm) 669 #else 670 ABI_ERROR("MPI-IO support not enabled") 671 #endif 672 end if 673 674 !BEGINDEBUG 675 ! if ( ANY(hmat==HUGE(zero)) ) then 676 ! write(std_out,*)"COUNT",COUNT(hmat==HUGE(zero))," hsize= ",hsize 677 ! ABI_ERROR("Something wrong in the reading") 678 ! end if 679 !ENDDEBUG 680 681 call xmpi_barrier(comm) 682 683 return 684 685 ! Handle IO Error 686 10 continue 687 ABI_ERROR(errmsg) 688 689 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.
SOURCE
198 subroutine exc_skip_bshdr(funt,ierr) 199 200 !Arguments ------------------------------------ 201 integer,intent(in) :: funt 202 integer,intent(out) :: ierr 203 204 !Local variables------------------------------- 205 character(len=500) :: errmsg 206 ! ************************************************************************* 207 208 call hdr_skip(funt,ierr) 209 if (ierr/=0) RETURN 210 read(funt, err=10, iomsg=errmsg) 211 212 return 213 214 ! Handle IO Error 215 10 continue 216 ierr = 0 217 ABI_WARNING(errmsg) 218 219 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
SOURCE
240 subroutine exc_skip_bshdr_mpio(mpifh,at_option,ehdr_offset) 241 242 !Arguments ------------------------------------ 243 integer,intent(in) :: mpifh,at_option 244 integer(XMPI_OFFSET_KIND),intent(inout) :: ehdr_offset 245 246 !Local variables ------------------------------ 247 !scalars 248 integer :: fform,ierr 249 #ifdef HAVE_MPI_IO 250 integer(XMPI_OFFSET_KIND) :: fmarker 251 #endif 252 253 ! ************************************************************************* 254 255 call hdr_mpio_skip(mpifh,fform,ehdr_offset) 256 257 #ifdef HAVE_MPI_IO 258 call xmpio_read_frm(mpifh,ehdr_offset,at_option,fmarker,ierr) 259 !write(std_out,*)"fmarker last record ",fmarker 260 #else 261 ABI_ERROR("You should not be here") 262 #endif 263 264 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
SOURCE
90 subroutine exc_write_bshdr(funt,Bsp,Hdr) 91 92 !Arguments ------------------------------------ 93 integer,intent(in) :: funt 94 type(excparam),intent(in) :: BSp 95 type(hdr_type),intent(inout) :: Hdr 96 97 !Local variables ------------------------------ 98 !scalars 99 integer :: fform_1002 = 1002 ! TODO: change setup_bse so that Hdr_bse reflects the parameters of the run. 100 integer :: ierr 101 character(len=500) :: errmsg 102 ! ************************************************************************* 103 104 call hdr%fort_write(funt, fform_1002, ierr) 105 ABI_CHECK(ierr == 0, "hdr_fort_write returned ierr != 0") 106 write(funt, err=10, iomsg=errmsg) BSp%nreh,BSp%nkbz 107 108 return 109 110 ! Handle IO Error 111 10 continue 112 ABI_ERROR(errmsg) 113 114 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
SOURCE
1370 subroutine exc_write_optme(filname,minb,maxb,nkbz,nsppol,nq,opt_cvk,ierr) 1371 1372 !Arguments ------------------------------------ 1373 integer,intent(in) :: minb,maxb,nkbz,nsppol,nq 1374 character(len=fnlen),intent(in) :: filname 1375 complex(dpc),intent(in) :: opt_cvk(minb:maxb,minb:maxb,nkbz,nsppol,nq) 1376 integer,intent(out) :: ierr 1377 1378 !Local variables ------------------------------ 1379 !scalars 1380 integer :: ncid,cmplx_id,nband_id,nkbz_id,nsppol_id,nq_id,xyz_id 1381 integer :: minb_id,maxb_id,ome_id,iq,is,ik,ib,jb 1382 integer :: dimOME(6),dimKPT(2),dimQPT(2),dimSCA(0),start6(6),count6(6) 1383 real(dp) :: complex2(2) 1384 ! ************************************************************************* 1385 1386 ierr = 1 1387 1388 !1. Create netCDF file 1389 NCF_CHECK_MSG(nctk_open_create(ncid, filname, xmpi_comm_self), " create netcdf OME file") 1390 1391 !2. Define dimensions 1392 NCF_CHECK(nf90_def_dim(ncid,"xyz",3,xyz_id)) 1393 NCF_CHECK(nf90_def_dim(ncid,"cmplx",2,cmplx_id)) 1394 NCF_CHECK(nf90_def_dim(ncid,"nband",maxb-minb+1,nband_id)) 1395 NCF_CHECK(nf90_def_dim(ncid,"nkbz",nkbz,nkbz_id)) 1396 NCF_CHECK(nf90_def_dim(ncid,"nq",nq,nq_id)) 1397 NCF_CHECK(nf90_def_dim(ncid,"nsppol",nsppol,nsppol_id)) 1398 1399 !Dimensions for optical matrix elements 1400 dimOME = (/ cmplx_id, nband_id, nband_id, nkbz_id, nsppol_id, nq_id /) 1401 !Dimensions for kpoint positions 1402 dimKPT = (/ xyz_id, nkbz_id /) 1403 !Dimensions for qpoints for the optical limit 1404 dimQPT = (/ xyz_id, nq_id /) 1405 1406 !3. Define variables 1407 1408 call ab_define_var(ncid, dimOME, ome_id, NF90_DOUBLE,& 1409 "OME", "Values of optical matrix elements","Tobedone") 1410 1411 call ab_define_var(ncid, dimSCA, minb_id, NF90_INT,"minb",& 1412 "Minimum band index for the optical matrix elements", "Dimensionless") 1413 1414 call ab_define_var(ncid, dimSCA, maxb_id, NF90_INT,"maxb",& 1415 "Maximum band index for the optical matrix elements", "Dimensionless") 1416 1417 !4. End define mode 1418 NCF_CHECK(nf90_enddef(ncid)) 1419 1420 !5 Write scalars (minb and maxb) 1421 1422 NCF_CHECK(nf90_put_var(ncid, minb_id, minb)) 1423 NCF_CHECK(nf90_put_var(ncid, maxb_id, maxb)) 1424 1425 !6 Write optical matrix elements 1426 1427 do iq=1,nq 1428 do is=1,nsppol 1429 do ik=1,nkbz 1430 do ib=minb,maxb 1431 do jb=minb,maxb 1432 start6 = (/ 1, ib-minb+1, jb-minb+1, ik, is, iq /) 1433 count6 = (/ 2, 1, 1, 1, 1, 1 /) 1434 complex2 = (/ REAL(opt_cvk(ib,jb,ik,is,iq)),AIMAG(opt_cvk(ib,jb,ik,is,iq)) /) 1435 NCF_CHECK(nf90_put_var(ncid, ome_id, complex2, start = start6, count = count6)) 1436 end do 1437 end do 1438 end do 1439 end do 1440 end do 1441 1442 !7 Close file 1443 NCF_CHECK(nf90_close(ncid)) 1444 ierr = 0 1445 1446 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
SOURCE
1060 function offset_in_file(row_glob,col_glob,size_glob,nsblocks,sub_block,bsize_elm,bsize_frm) 1061 1062 !Arguments ------------------------------------ 1063 integer(XMPI_OFFSET_KIND) :: offset_in_file 1064 integer,intent(in) :: row_glob,col_glob,nsblocks,bsize_elm,bsize_frm 1065 integer,intent(in) :: size_glob(2),sub_block(2,2,nsblocks) 1066 1067 !Local variables ------------------------------ 1068 !scalars 1069 integer :: ii,jj,ijp_glob,swap 1070 integer(XMPI_OFFSET_KIND) :: my_offset 1071 1072 ! ************************************************************************* 1073 1074 if (nsblocks==1) then 1075 ii = row_glob 1076 jj = col_glob 1077 if (ii>size_glob(1)/2) ii = ii - size_glob(1)/2 ! Wrap the index. 1078 if (jj>size_glob(2)/2) jj = jj - size_glob(2)/2 1079 if (jj<ii) then ! Exchange the indices since the symmetric element is read. 1080 swap = jj 1081 jj = ii 1082 ii = swap 1083 end if 1084 ijp_glob = ii + jj*(jj-1)/2 ! Index for packed storage mode. 1085 my_offset = (ijp_glob-1)*bsize_elm + (jj-1)*2*bsize_frm 1086 else 1087 ABI_UNUSED(sub_block(1,1,1)) 1088 ABI_ERROR("nsppol==2 not coded") 1089 end if 1090 1091 offset_in_file = my_offset 1092 1093 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
SOURCE
973 pure function rrs_of_glob(row_glob,col_glob,size_glob) 974 975 !Arguments ------------------------------------ 976 integer :: rrs_of_glob 977 integer,intent(in) :: row_glob,col_glob 978 integer,intent(in) :: size_glob(2) 979 980 !Local variables ------------------------------ 981 !scalars 982 integer :: nreh1,nreh2 983 984 ! ************************************************************************* 985 986 nreh1=size_glob(1)/2 ! Matrix is square and nreh1==nreh2 but oh well. 987 nreh2=size_glob(2)/2 988 989 if (row_glob<=nreh1 .and. col_glob<=nreh2) then 990 rrs_of_glob=+1 ! Resonant. 991 else if (row_glob >nreh1 .and. col_glob >nreh2) then 992 rrs_of_glob=-1 ! anti-Resonant. 993 else 994 rrs_of_glob=0 995 end if 996 997 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
SOURCE
1466 subroutine exc_ham_ncwrite(ncid,Kmesh,BSp,hsize,nreh,vcks2t,hreso,diag) 1467 1468 !Arguments ------------------------------------ 1469 !scalars 1470 integer,intent(in) :: ncid 1471 integer,intent(in) :: hsize 1472 type(kmesh_t),intent(in) :: Kmesh 1473 type(excparam),intent(in) :: BSp 1474 integer,intent(in) :: nreh(BSp%nsppol) 1475 integer,target,intent(in) :: vcks2t(BSp%maxnbndv,BSp%maxnbndc,Kmesh%nbz,BSp%nsppol) 1476 complex(dpc),target,intent(in) :: hreso(hsize,hsize) 1477 complex(dpc),target,intent(in) :: diag(hsize) 1478 1479 !Local variables------------------------------- 1480 integer :: ncerr 1481 integer :: max_nreh, sum_nreh 1482 real(dp), ABI_CONTIGUOUS pointer :: r2vals(:,:),r3vals(:,:,:) 1483 1484 ! ************************************************************************* 1485 1486 ! ============================================== 1487 ! === Write the dimensions specified by ETSF === 1488 ! ============================================== 1489 max_nreh = MAXVAL(nreh) 1490 sum_nreh = SUM(nreh) 1491 1492 ncerr = nctk_def_dims(ncid, [nctkdim_t("number_of_reduced_dimensions", 3), nctkdim_t("number_of_spins", bsp%nsppol),& 1493 nctkdim_t("number_of_kpoints", kmesh%nbz), nctkdim_t("max_number_of_valence_bands", bsp%maxnbndv),& 1494 nctkdim_t("max_number_of_conduction_bands", bsp%maxnbndc), nctkdim_t("max_number_of_transitions", max_nreh),& 1495 nctkdim_t("total_number_of_transitions", sum_nreh), nctkdim_t("cplex", 2)], defmode=.True.) 1496 NCF_CHECK(ncerr) 1497 1498 ncerr = nctk_def_arrays(ncid, [& 1499 nctkarr_t('vcks2t', "i", 'max_number_of_valence_bands, max_number_of_conduction_bands, number_of_kpoints, number_of_spins'),& 1500 nctkarr_t('hamiltonian', "dp", 'cplex total_number_of_transitions total_number_of_transitions'),& 1501 nctkarr_t('diagonal', "dp", 'cplex total_number_of_transitions'),& 1502 nctkarr_t('lomo', "dp", 'number_of_spins'), & 1503 nctkarr_t('homo', "dp", 'number_of_spins'), & 1504 nctkarr_t('lumo', "dp", 'number_of_spins'), & 1505 nctkarr_t('humo', "dp", 'number_of_spins'), & 1506 nctkarr_t("reduced_coordinates_of_kpoints", "dp", "number_of_reduced_dimensions, number_of_kpoints"), & 1507 nctkarr_t("kpoint_weights", "dp", "number_of_kpoints") & 1508 ]) 1509 NCF_CHECK(ncerr) 1510 1511 ! Write data 1512 NCF_CHECK(nctk_set_datamode(ncid)) 1513 NCF_CHECK(nf90_put_var(ncid, vid('vcks2t'), vcks2t)) 1514 NCF_CHECK(nf90_put_var(ncid, vid('lomo'), BSp%lomo_spin)) 1515 NCF_CHECK(nf90_put_var(ncid, vid('homo'), BSp%homo_spin)) 1516 NCF_CHECK(nf90_put_var(ncid, vid('lumo'), BSp%lumo_spin)) 1517 NCF_CHECK(nf90_put_var(ncid, vid('humo'), BSp%humo_spin)) 1518 1519 call c_f_pointer(c_loc(hreso(1,1)), r3vals, shape=[2, hsize, hsize]) 1520 NCF_CHECK(nf90_put_var(ncid, vid("hamiltonian"), r3vals)) 1521 1522 call c_f_pointer(c_loc(diag(1)), r2vals, shape=[2, hsize]) 1523 NCF_CHECK(nf90_put_var(ncid, vid("diag"), r2vals)) 1524 1525 ! Write K-points 1526 NCF_CHECK(nf90_put_var(ncid, vid("reduced_coordinates_of_kpoints"), kmesh%bz)) 1527 NCF_CHECK(nf90_put_var(ncid, vid("kpoint_weights"), kmesh%wt)) 1528 NCF_CHECK(ncerr) 1529 1530 contains 1531 integer function vid(vname) 1532 character(len=*),intent(in) :: vname 1533 vid = nctk_idname(ncid, vname) 1534 end function vid 1535 1536 end subroutine exc_ham_ncwrite