Line data Source code
1 : !!****m* ABINIT/m_bse_io
2 : !! NAME
3 : !! m_bse_io
4 : !!
5 : !! FUNCTION
6 : !! This module provides routines to read the Bethe-Salpeter Hamiltonian from file
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2008-2026 ABINIT group (MG)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
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, intrinsic :: iso_c_binding
25 : use defs_basis
26 : USE_MPI
27 : use m_xmpi
28 : use m_errors
29 : use m_abicore
30 : use netcdf
31 : use m_nctk
32 : use m_hdr
33 :
34 : use m_time, only : cwtime
35 : use m_fstrings, only : toupper
36 : use m_io_tools, only : open_file
37 : use m_numeric_tools, only : arth
38 : use m_special_funcs, only : gaussian
39 : use m_bs_defs, only : excparam
40 : use m_bz_mesh, only : kmesh_t
41 :
42 : implicit none
43 :
44 : #if defined HAVE_MPI1
45 : include 'mpif.h'
46 : #endif
47 :
48 : private
49 :
50 : public :: exc_read_rblock_fio ! Reads the entire resonant sub-block from file using Fortran IO.
51 : public :: exc_read_rcblock ! Reads a distributed sub-block of the excitonic Hamiltonian from file.
52 : public :: exc_fullh_from_blocks ! Initialize the specified sub-blocks of the *full* matrix (reso+anti-reso) from file.
53 : public :: rrs_of_glob ! [+1,-1,0] if (row_glob,col_glob) belongs to the [ resonant, anti-resonant, (anti)coupling block ]
54 : public :: ccs_of_glob ! [+1,-1,0] if (row_glob,col_glob) belongs to the [ coupling, anti-coupling, (anti)resonant block ]
55 : public :: offset_in_file ! Function used to describe the way the Hamiltonian is stored on disk.
56 : public :: exc_write_bshdr ! Writes the Header of the (BSR|BSC) files storing the excitonic Hamiltonian.
57 : public :: exc_read_bshdr ! Reads the Header of the (BSR|BSC) files.
58 : public :: exc_skip_bshdr ! Skip the Header of the (BSR|BSC) files. Fortran version.
59 : public :: exc_skip_bshdr_mpio ! Skip the Header of the (BSR|BSC) files. MPI-IO version.
60 : public :: exc_read_eigen ! Read selected energies and eigenvectors from the BSEIG file.
61 : 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.
62 : public :: exc_write_optme ! Writes the OME file storing the optical matrix elements
63 : public :: exc_ham_ncwrite ! Writes the hamiltonian in NETCDF format
64 : !!***
65 :
66 : CONTAINS !====================================================================
67 : !!***
68 :
69 : !----------------------------------------------------------------------
70 :
71 : !!****f* m_bse_io/exc_write_bshdr
72 : !! NAME
73 : !! exc_write_bshdr
74 : !!
75 : !! FUNCTION
76 : !! Writes the header of the (BSR|BSC) files storing the excitonic Hamiltonian.
77 : !!
78 : !! INPUTS
79 : !! funt=Fortran unit number.
80 : !! Bsp<excparam>=Structure storing the parameters of the run.
81 : !! Hdr<hdr_type>=The abinit header.
82 : !!
83 : !! OUTPUT
84 : !! Only writing
85 : !!
86 : !! SOURCE
87 :
88 22 : subroutine exc_write_bshdr(funt,Bsp,Hdr)
89 :
90 : !Arguments ------------------------------------
91 : integer,intent(in) :: funt
92 : type(excparam),intent(in) :: BSp
93 : type(hdr_type),intent(inout) :: Hdr
94 :
95 : !Local variables ------------------------------
96 : !scalars
97 : integer :: fform_1002 = 1002 ! TODO: change setup_bse so that Hdr_bse reflects the parameters of the run.
98 : integer :: ierr
99 : character(len=500) :: errmsg
100 : ! *************************************************************************
101 :
102 22 : call hdr%fort_write(funt, fform_1002, ierr)
103 22 : ABI_CHECK(ierr == 0, "hdr_fort_write returned ierr != 0")
104 22 : write(funt, err=10, iomsg=errmsg) BSp%nreh,BSp%nkbz
105 :
106 22 : return
107 :
108 : ! Handle IO Error
109 : 10 continue
110 0 : ABI_ERROR(errmsg)
111 :
112 : end subroutine exc_write_bshdr
113 : !!***
114 :
115 : !----------------------------------------------------------------------
116 :
117 : !!****f* m_bse_io/exc_read_bshdr
118 : !! NAME
119 : !! exc_read_bshdr
120 : !!
121 : !! FUNCTION
122 : !! Reads the header of the (BSR|BSC) files storing the excitonic Hamiltonian.
123 : !! and performs basilar consistency checks.
124 : !!
125 : !! INPUTS
126 : !! funt=Unit number.
127 : !! Bsp<excparam>=Structure storing the parameters of the run.
128 : !! Hdr<hdr_type>=The abinit header.
129 : !!
130 : !! OUTPUT
131 : !! fform=Integer defining the file format.
132 : !! ierr=Status error.
133 : !!
134 : !! SOURCE
135 :
136 51 : subroutine exc_read_bshdr(funt,Bsp,fform,ierr)
137 :
138 : !Arguments ------------------------------------
139 : integer,intent(in) :: funt
140 : integer,intent(out) :: fform,ierr
141 : type(excparam),intent(in) :: BSp
142 :
143 : !Local variables ------------------------------
144 : !scalars
145 : integer :: nkbz_read
146 : character(len=500) :: errmsg
147 51 : type(hdr_type) :: Hdr
148 : !arrays
149 102 : integer :: nreh_read(SIZE(BSp%nreh))
150 : ! *************************************************************************
151 :
152 51 : ierr=0
153 :
154 : ! Read the header and perform consistency checks.
155 51 : call hdr%fort_read(funt, fform, rewind=.True.)
156 51 : ABI_CHECK(fform /= 0, "hdr_fort_read returned fform == 0")
157 :
158 51 : read(funt, err=10, iomsg=errmsg) nreh_read, nkbz_read
159 51 : call Hdr%free()
160 :
161 104 : if (ANY(nreh_read /= BSp%nreh)) then
162 0 : call wrtout(std_out,"Wrong number of e-h transitions")
163 0 : ierr = ierr + 1
164 : end if
165 :
166 : return
167 :
168 0 : 10 ierr = 1
169 0 : ABI_WARNING(errmsg)
170 :
171 51 : end subroutine exc_read_bshdr
172 : !!***
173 :
174 : !----------------------------------------------------------------------
175 :
176 : !!****f* m_bse_io/exc_skip_bshdr
177 : !! NAME
178 : !! exc_skip_bshdr
179 : !!
180 : !! FUNCTION
181 : !! Skip the header of the (BSR|BSC) files storing the excitonic Hamiltonian. Fortran version.
182 : !!
183 : !! INPUTS
184 : !! funt=Unit number.
185 : !!
186 : !! OUTPUT
187 : !! ierr=Status error.
188 : !!
189 : !! SIDE EFFECTS
190 : !! Skip the header.
191 : !!
192 : !! SOURCE
193 :
194 22 : subroutine exc_skip_bshdr(funt,ierr)
195 :
196 : !Arguments ------------------------------------
197 : integer,intent(in) :: funt
198 : integer,intent(out) :: ierr
199 :
200 : !Local variables-------------------------------
201 : character(len=500) :: errmsg
202 : ! *************************************************************************
203 :
204 22 : call hdr_skip(funt,ierr)
205 44 : if (ierr/=0) RETURN
206 22 : read(funt, err=10, iomsg=errmsg)
207 :
208 22 : return
209 :
210 : ! Handle IO Error
211 : 10 continue
212 0 : ierr = 0
213 0 : ABI_WARNING(errmsg)
214 :
215 : end subroutine exc_skip_bshdr
216 : !!***
217 :
218 : !----------------------------------------------------------------------
219 :
220 : !!****f* m_bse_io/exc_skip_bshdr_mpio
221 : !! NAME
222 : !! exc_skip_bshdr_mpio
223 : !!
224 : !! FUNCTION
225 : !! Skip the header of the (BSR|BSC) files storing the excitonic Hamiltonian. MPI-IO version.
226 : !!
227 : !! INPUTS
228 : !! mpifh=MPI-IO file handler.
229 : !! at_option
230 : !!
231 : !! SIDE EFFECTS
232 : !! ehdr_offset
233 : !!
234 : !! SOURCE
235 :
236 0 : subroutine exc_skip_bshdr_mpio(mpifh,at_option,ehdr_offset)
237 :
238 : !Arguments ------------------------------------
239 : integer,intent(in) :: mpifh,at_option
240 : integer(XMPI_OFFSET_KIND),intent(inout) :: ehdr_offset
241 :
242 : !Local variables ------------------------------
243 : integer :: fform,ierr
244 : #ifdef HAVE_MPI_IO
245 : integer(XMPI_OFFSET_KIND) :: fmarker
246 : #endif
247 : ! *************************************************************************
248 :
249 0 : call hdr_mpio_skip(mpifh,fform,ehdr_offset)
250 :
251 : #ifdef HAVE_MPI_IO
252 0 : call xmpio_read_frm(mpifh,ehdr_offset,at_option,fmarker,ierr)
253 : !write(std_out,*)"fmarker last record ",fmarker
254 : #else
255 : ABI_ERROR("You should not be here")
256 : #endif
257 :
258 0 : end subroutine exc_skip_bshdr_mpio
259 : !!***
260 :
261 : !----------------------------------------------------------------------
262 :
263 : !!****f* m_bse_io/exc_read_eigen
264 : !! NAME
265 : !! exc_read_eigen
266 : !!
267 : !! FUNCTION
268 : !! Read selected energies and eigenvectors from the BSEIG file.
269 : !!
270 : !! INPUTS
271 : !! eig_fname=The name of the file storing the excitonic eigenvectors.
272 : !! hsize=Size of the Hamiltonian.
273 : !! nvec=Number of excitonic states to analyze.
274 : !! vec_idx(nvec)=List with the indices of the excitonic states sorted in ascending order.
275 : !! [Bsp]<excparam>=Structure storing the parameters of the run. If present the
276 : !! routine will perform additional consistency checks to make sure that
277 : !! the content of the file is consistent with the present run.
278 : !!
279 : !! OUTPUT
280 : !! [ene_list(nvec)]=Excitonic energies
281 : !! vec_list(hsize,nvec)=Excitonic eigenvectors.
282 : !!
283 : !! SOURCE
284 :
285 0 : subroutine exc_read_eigen(eig_fname,hsize,nvec,vec_idx,vec_list,ene_list,Bsp)
286 :
287 : !Arguments ------------------------------------
288 : !scalars
289 : integer,intent(in) :: nvec,hsize
290 : character(len=*),intent(in) :: eig_fname
291 : type(excparam),optional,intent(in) :: BSp
292 : ! arrays
293 : integer,intent(in) :: vec_idx(nvec)
294 : real(dp),optional,intent(out) :: ene_list(nvec)
295 : complex(dp),intent(out) :: vec_list(hsize,nvec)
296 :
297 : !Local variables ------------------------------
298 : !scalars
299 : integer :: eig_unt,hsize_read,neig_read,ll,vec
300 : character(len=500) :: msg,errmsg
301 : !arrays
302 : !real(dp),allocatable :: exc_ene(:)
303 0 : complex(dp),allocatable :: exc_ene_cplx(:)
304 : ! *************************************************************************
305 :
306 : ABI_UNUSED(BSp%nline)
307 :
308 0 : if (open_file(eig_fname,msg,newunit=eig_unt,form="unformatted",status="old",action="read") /= 0) then
309 0 : ABI_ERROR(msg)
310 : end if
311 :
312 0 : read(eig_unt, err=10, iomsg=errmsg)hsize_read,neig_read
313 : !write(std_out,*)hsize_read, neig_read
314 :
315 0 : if (hsize_read/=hsize) then
316 0 : write(msg,'(a,2(1x,i0))')" hsize_read/=hsize: ",hsize_read,hsize
317 0 : ABI_ERROR(msg)
318 : end if
319 :
320 : ! Read eigenvalues, ignore possibly small imaginary part.
321 0 : ABI_MALLOC(exc_ene_cplx,(neig_read))
322 0 : read(eig_unt, err=10, iomsg=errmsg) exc_ene_cplx
323 :
324 0 : if (PRESENT(ene_list)) then
325 0 : do vec=1,nvec
326 0 : ll = vec_idx(vec)
327 0 : ene_list(vec) = DBLE(exc_ene_cplx(ll))
328 : end do
329 : end if
330 0 : ABI_FREE(exc_ene_cplx)
331 :
332 0 : vec=1
333 0 : do ll=1,neig_read ! Read the selected excitons.
334 0 : if (ll==vec_idx(vec)) then
335 0 : read(eig_unt, err=10, iomsg=errmsg) vec_list(:,vec)
336 0 : if (vec==nvec) EXIT
337 0 : vec=vec+1
338 : else
339 0 : read(eig_unt, err=10, iomsg=errmsg)
340 : end if
341 : end do
342 :
343 0 : close(eig_unt, err=10, iomsg=errmsg)
344 :
345 0 : if (vec/=nvec) then
346 0 : write(msg,'(a,2(1x,i0))')" vec_idx is wrong, vec/=nvec ",vec,nvec+1
347 0 : ABI_ERROR(msg)
348 : end if
349 :
350 : return
351 :
352 : ! Handle IO-error
353 : 10 continue
354 0 : ABI_ERROR(errmsg)
355 :
356 0 : end subroutine exc_read_eigen
357 : !!***
358 :
359 : !----------------------------------------------------------------------
360 :
361 : !!****f* m_bse_io/exc_read_rcblock
362 : !! NAME
363 : !! exc_read_rcblock
364 : !!
365 : !! FUNCTION
366 : !! Reads the excitonic Hamiltonian from file
367 : !!
368 : !! INPUTS
369 : !! fname=File name.
370 : !! diago_is_real=.TRUE. if diagonal elements are real (used only if is_resonant==.TRUE.)
371 : !! nreh(nsppol)=Number of resonant transition for the two spins.
372 : !! is_resonant=Set to .TRUE. if the block is resonant.
373 : !! hsize=Dimension of the block.
374 : !! nsppol=2 for spin polarized systems. 1 otherwise.
375 : !! my_t1,my_t2=The first and the last columns of the matrix treated by this node.
376 : !! use_mpio=.TRUE. is MPI-IO routines are used.
377 : !! comm=MPI communicator.
378 : !!
379 : !! OUTPUT
380 : !! hmat(hsize,my_t1:my_t2)=The block read from file fname.
381 : !!
382 : !! TODO
383 : !! Remove Bsp
384 : !!
385 : !! SOURCE
386 :
387 30 : subroutine exc_read_rcblock(fname,Bsp,is_resonant,diago_is_real,nsppol,nreh,hsize,my_t1,my_t2,hmat,use_mpio,comm)
388 :
389 : !Arguments ------------------------------------
390 : !scalars
391 : integer,intent(in) :: comm,hsize,my_t1,my_t2,nsppol
392 : logical,intent(in) :: is_resonant,use_mpio,diago_is_real
393 : character(len=*),intent(in) :: fname
394 : type(excparam),intent(in) :: Bsp
395 : !arrays
396 : integer,intent(in) :: nreh(nsppol)
397 : complex(dp),intent(out) :: hmat(hsize,my_t1:my_t2)
398 :
399 : !Local variables ------------------------------
400 : !scalars
401 : integer,parameter :: master=0
402 : integer :: it,itp,funit,nproc,my_rank,neh1,neh2
403 : integer :: fform,my_nt
404 : integer :: row,col,block,spad,spin_dim,ierr,size_exp
405 : real(dp) :: cputime,walltime,gflops
406 : character(len=500) :: msg,errmsg
407 : !type(Hdr_type) :: bse_Hdr
408 : !arrays
409 30 : complex(dp),allocatable :: buffer_dpc(:)
410 : logical :: have_row,have_col
411 : #ifdef HAVE_MPI_IO
412 : integer :: mpierr,mpifh,ham_type,my_nel,old_type,etype,offset_err,amode
413 : integer :: irec,nrec !,ncount
414 : integer(XMPI_OFFSET_KIND) :: ehdr_offset,my_offset,my_offpad,fsize
415 30 : integer(XMPI_OFFSET_KIND),allocatable :: bsize_frecord(:)
416 : integer :: glob_sizes(2),my_cols(2), block_sizes(2,3), status(MPI_STATUS_SIZE)
417 : #endif
418 : !************************************************************************
419 :
420 30 : nproc = xmpi_comm_size(comm); my_rank= xmpi_comm_rank(comm)
421 :
422 : ! Compute the (Expected) size of the hamiltonian.
423 30 : neh1 = nreh(1)
424 30 : neh2 = neh1; if (nsppol==2) neh2=nreh(2)
425 :
426 30 : size_exp=neh1; if (nsppol==2) size_exp=SUM(nreh)
427 :
428 30 : ABI_CHECK(hsize==size_exp,"Wrong hsize")
429 30 : if (neh1/=neh2) then
430 0 : ABI_ERROR("BSE code does not support different number of transitions for the two spin channels")
431 : end if
432 :
433 30 : my_nt = my_t2-my_t1+1
434 : !hmat = HUGE(zero)
435 :
436 30 : if (.not.use_mpio) then
437 :
438 30 : if (is_resonant) then
439 30 : call wrtout(std_out,". Reading resonant block from file: "//TRIM(fname)//" using Fortran-IO")
440 : else
441 0 : call wrtout(std_out,". Reading coupling block from file: "//TRIM(fname)//" using Fortran-IO")
442 : end if
443 :
444 30 : if (open_file(fname,msg,newunit=funit,form="unformatted",status="old",action="read") /= 0) then
445 0 : ABI_ERROR(msg)
446 : end if
447 : !
448 : ! Read the header and perform consistency checks.
449 30 : call exc_read_bshdr(funit,Bsp,fform,ierr)
450 30 : ABI_CHECK(ierr==0,"Wrong BSE header")
451 : !
452 : ! Construct full excitonic Hamiltonian using symmetries.
453 30 : if (nsppol==1) then
454 30 : call cwtime(cputime,walltime,gflops,"start")
455 : !
456 90 : ABI_MALLOC(buffer_dpc,(neh1))
457 8990 : do itp=1,hsize
458 8960 : read(funit, err=10, iomsg=errmsg) buffer_dpc(1:itp)
459 : !
460 : ! Fill the upper triangle if I have this column.
461 8960 : if (itp>=my_t1 .and. itp<=my_t2) then
462 1704832 : do it=1,itp
463 1704832 : hmat(it,itp) = buffer_dpc(it)
464 : end do
465 : ! Force the diagonal to be real.
466 7424 : if (is_resonant .and.diago_is_real) hmat(itp,itp) = DBLE(hmat(itp,itp))
467 : end if
468 : !
469 : ! Reconstruct the rows below the diagonal (diagonal part is not touched here).
470 8990 : if (is_resonant) then ! Use Hermiticity.
471 1993088 : do it=1,itp-1
472 1993088 : if (it>=my_t1 .and. it<=my_t2) hmat(itp,it) = CONJG(buffer_dpc(it))
473 : end do
474 : else ! Coupling is symmetric.
475 0 : do it=1,itp-1
476 0 : if (it>=my_t1 .and. it<=my_t2) hmat(itp,it) = buffer_dpc(it)
477 : end do
478 : end if
479 : !
480 : end do ! itp
481 30 : ABI_FREE(buffer_dpc)
482 :
483 30 : call cwtime(cputime,walltime,gflops,"stop")
484 30 : write(msg,'(2(a,f9.1),a)')" Fortran-IO completed. cpu_time ",cputime,"[s], walltime ",walltime," [s]"
485 90 : call wrtout(std_out, msg, do_flush=.True.)
486 : else
487 : ! Spin polarized case.
488 : !
489 : ! The file contains
490 : ! A) The up-up and the down-down block in packed form
491 : ! (only the upper triangle is stored since the blocks are Hermitian)
492 : ! B) The entire up-down exchange block (no symmetry here)
493 : !
494 : ! A) Construct resonant blocks from the upper triangles stored on file.
495 : ! FIXME this part won't work if we have a different number of e-h pairs
496 0 : if (.not.is_resonant) then
497 0 : ABI_ERROR("exc_read_rcblock does not support coupling.")
498 : end if
499 : ! It should be checked.
500 0 : spin_dim=neh1
501 0 : do block=1,2
502 0 : ABI_MALLOC(buffer_dpc,(neh1))
503 0 : if (block==1) spad=0
504 0 : if (block==2) spad=neh1
505 0 : do itp=1,spin_dim
506 : !
507 0 : read(funit, err=10, iomsg=errmsg) buffer_dpc(1:itp)
508 : !
509 : ! Fill the upper triangle if this node treats this column
510 0 : col = itp+spad
511 0 : if (col>=my_t1 .and. col<=my_t2) then
512 0 : do it=1,itp
513 0 : row = it + spad
514 0 : hmat(row,col) = buffer_dpc(it)
515 : end do
516 : ! Force the diagonal to be real.
517 0 : if (is_resonant .and.diago_is_real) hmat(col,col) = DBLE(hmat(col,col))
518 : end if
519 : !
520 : ! Reconstruct the rows below the diagonal (diagonal part is not touched here).
521 0 : row = itp + spad
522 0 : if (is_resonant) then ! Hermitian
523 0 : do it=1,itp-1
524 0 : col = it + spad
525 0 : if (col>=my_t1 .and. col<=my_t2) hmat(row,col) = CONJG(buffer_dpc(it))
526 : end do
527 : else ! Coupling is symmetric.
528 0 : do it=1,itp-1
529 0 : col = it + spad
530 0 : if (col>=my_t1 .and. col<=my_t2) hmat(row,col) = buffer_dpc(it)
531 : end do
532 : end if
533 : !
534 : end do ! itp
535 0 : ABI_FREE(buffer_dpc)
536 : end do ! block
537 : !
538 : ! B) Kx_{down up} = Kx_{up down}^H.
539 : ! FIXME this part won't work if we have a different number of e-h pairs
540 0 : spad=neh1
541 0 : spin_dim=neh1
542 0 : ABI_MALLOC(buffer_dpc,(neh1))
543 0 : do itp=1,spin_dim
544 0 : read(funit, err=10, iomsg=errmsg) buffer_dpc(1:spin_dim)
545 0 : have_col = (spad+itp>=my_t1 .and. spad+itp<=my_t2)
546 0 : if (have_col) hmat(1:spin_dim,spad+itp) = buffer_dpc(1:spin_dim)
547 : ! Construct and store the lower block
548 0 : if (is_resonant) then ! Hermitian
549 0 : do it=1,spin_dim
550 0 : have_row = (it>=my_t1 .and. it<=my_t2)
551 0 : if (have_row) hmat(spad+itp,it) = CONJG(buffer_dpc(it))
552 : end do
553 : else ! Symmetric
554 0 : do it=1,spin_dim
555 0 : have_row = (it>=my_t1 .and. it<=my_t2)
556 0 : if (have_row) hmat(spad+itp,it) = buffer_dpc(it)
557 : end do
558 : end if
559 : end do
560 0 : ABI_FREE(buffer_dpc)
561 : end if
562 :
563 30 : close(funit)
564 :
565 : else
566 : #ifdef HAVE_MPI_IO
567 0 : if (is_resonant) then
568 0 : call wrtout(std_out,". Reading resonant block from file "//TRIM(fname)//" using MPI-IO")
569 : else
570 0 : call wrtout(std_out,". Reading coupling block from file "//TRIM(fname)//" using MPI-IO")
571 : end if
572 :
573 0 : amode=MPI_MODE_RDONLY
574 0 : call MPI_FILE_OPEN(comm,fname,amode,MPI_INFO_NULL,mpifh,mpierr)
575 0 : msg = " FILE_OPEN "//TRIM(fname)
576 0 : ABI_CHECK_MPI(mpierr,msg)
577 :
578 0 : call MPI_FILE_GET_SIZE(mpifh,fsize,mpierr)
579 : !write(std_out,*)" file size is ",fsize
580 : !
581 : ! Skip the header and find the offset for reading the matrix.
582 0 : call exc_skip_bshdr_mpio(mpifh,xmpio_collective,ehdr_offset)
583 : !
584 : ! Read my columns from file.
585 0 : old_type=MPI_DOUBLE_COMPLEX
586 0 : glob_sizes = (/hsize,hsize/); my_cols=(/my_t1,my_t2/)
587 :
588 0 : if (nsppol==1) then
589 0 : call xmpio_create_coldistr_from_fpacked(glob_sizes,my_cols,old_type,ham_type,my_offpad,offset_err)
590 : else
591 0 : ABI_WARNING("nsppol==2 => calling fp3blocks")
592 0 : write(std_out,*)"neh, hsize",neh1,neh2,hsize
593 :
594 0 : nrec=neh1+2*neh2
595 0 : ABI_MALLOC(bsize_frecord,(nrec))
596 0 : bsize_frecord(1:neh1) = (/(irec*xmpi_bsize_dpc, irec=1,neh1)/)
597 0 : bsize_frecord(neh1+1:neh1+neh2) = (/(irec*xmpi_bsize_dpc, irec=1,neh2)/)
598 0 : bsize_frecord(neh1+neh2+1:) = neh1*xmpi_bsize_dpc
599 0 : call xmpio_check_frmarkers(mpifh,ehdr_offset,xmpio_collective,nrec,bsize_frecord,ierr)
600 0 : ABI_CHECK(ierr==0,"Error in Fortran markers")
601 0 : ABI_FREE(bsize_frecord)
602 0 : ABI_COMMENT("Marker check ok")
603 0 : call xmpi_barrier(comm)
604 :
605 0 : block_sizes(:,1) = (/neh1,neh1/)
606 0 : block_sizes(:,2) = (/neh2,neh2/)
607 0 : block_sizes(:,3) = (/neh1,neh2/)
608 0 : ABI_ERROR("fp3blocks is buggy")
609 0 : call xmpio_create_coldistr_from_fp3blocks(glob_sizes,block_sizes,my_cols,old_type,ham_type,my_offpad,offset_err)
610 : end if
611 :
612 0 : if (offset_err/=0) then
613 : write(msg,"(3a)")&
614 0 : "Global position index cannot be stored in a standard Fortran integer ",ch10,&
615 0 : "Excitonic matrix cannot be read with a single MPI-IO call."
616 0 : ABI_ERROR(msg)
617 : end if
618 : !
619 : ! The offset used for reading.
620 0 : my_offset = ehdr_offset + my_offpad
621 0 : write(std_out,*)"my_offset= ",my_offset
622 :
623 0 : etype=MPI_BYTE
624 0 : call MPI_FILE_SET_VIEW(mpifh, my_offset, etype, ham_type, 'native', MPI_INFO_NULL, mpierr)
625 0 : ABI_CHECK_MPI(mpierr,"SET_VIEW")
626 :
627 : ! Release the MPI filetype.
628 0 : call MPI_TYPE_FREE(ham_type,mpierr)
629 0 : ABI_CHECK_MPI(mpierr,"TYPE_FREE")
630 :
631 0 : my_nel = my_nt*hsize
632 0 : call MPI_FILE_READ_ALL(mpifh, hmat, my_nel, MPI_DOUBLE_COMPLEX, status, mpierr)
633 0 : ABI_CHECK_MPI(mpierr,"READ_ALL")
634 :
635 : !call MPI_GET_COUNT(status, MPI_DOUBLE_COMPLEX, ncount, mpierr)
636 : !write(std_out,*)"count, my_nel ",ncount,my_nel
637 : !ABI_CHECK_MPI(mpierr,"READ_ALL")
638 : !
639 : ! Close the file.
640 0 : call MPI_FILE_CLOSE(mpifh, mpierr)
641 0 : ABI_CHECK_MPI(mpierr,"FILE_CLOSE")
642 : !
643 : ! Use the symmetries of the block to reconstruct the local buffer.
644 : ! Coupling does not require in-place symmetrization since it is symmetric.
645 0 : if (is_resonant) then
646 0 : do itp=my_t1,my_t2
647 0 : if (itp+1<=hsize) hmat(itp+1:,itp) = DCONJG(hmat(itp+1:,itp)) ! Lower triangle using Hermiticity.
648 0 : if (diago_is_real) hmat(itp,itp) = DBLE(hmat(itp,itp)) ! The diagonal is forced to be real when energies are real.
649 : end do
650 : end if
651 :
652 0 : call xmpi_barrier(comm)
653 : #else
654 : ABI_ERROR("MPI-IO support not enabled")
655 : #endif
656 : end if
657 :
658 : !BEGINDEBUG
659 : ! if ( ANY(hmat==HUGE(zero)) ) then
660 : ! write(std_out,*)"COUNT",COUNT(hmat==HUGE(zero))," hsize= ",hsize
661 : ! ABI_ERROR("Something wrong in the reading")
662 : ! end if
663 : !ENDDEBUG
664 :
665 30 : call xmpi_barrier(comm)
666 :
667 : return
668 :
669 : ! Handle IO Error
670 : 10 continue
671 0 : ABI_ERROR(errmsg)
672 :
673 0 : end subroutine exc_read_rcblock
674 : !!***
675 :
676 : !----------------------------------------------------------------------
677 :
678 : !!****f* m_bse_io/exc_fullh_from_blocks
679 : !! NAME
680 : !! exc_fullh_from_blocks
681 : !!
682 : !! FUNCTION
683 : !! Construct the matrix F H
684 : !!
685 : !! INPUTS
686 : !! funt
687 : !! nsppol
688 : !! block_type
689 : !! "Resonant"
690 : !! "Coupling"
691 : !! row_sign
692 : !! -1 to read ( R C )
693 : !! (-C* -R*)
694 : !!
695 : !! +1 to read ( R C )
696 : !! ( C* R*)
697 : !! diago_is_real=Used when block_type=resonat to specify whether the diagonal matrix elements are
698 : !! real or complex (when QP linewidth are included)
699 : !! nreh(nsppol)
700 : !! exc_size=Size of the full excitonic Hamiltonian.
701 : !!
702 : !! SIDE EFFECTS
703 : !! exc_ham(exc_size,exc_size)
704 : !!
705 : !! NOTES
706 : !!
707 : !! SOURCE
708 :
709 2 : subroutine exc_fullh_from_blocks(funt,block_type,nsppol,row_sign,diago_is_real,nreh,exc_size,exc_ham)
710 :
711 : !Arguments ------------------------------------
712 : !scalars
713 : integer,intent(in) :: funt,exc_size,nsppol,row_sign
714 : logical,intent(in) :: diago_is_real
715 : character(len=*),intent(in) :: block_type
716 : !arrays
717 : integer,intent(in) :: nreh(nsppol)
718 : complex(dp),intent(inout) :: exc_ham(exc_size,exc_size)
719 :
720 : !Local variables-------------------------------
721 : !scalars
722 : integer :: it,itp,szbuf,neh,pad_c1,pad_r1,spin_dim,spad_r,spad_c
723 : integer :: block,spad,row1,col1,row2,col2,spin_stride,ierr
724 : complex(dp) :: cttp
725 : character(len=500) :: errmsg
726 : !arrays
727 2 : complex(dp),allocatable :: cbuff_dpc(:)
728 : ! *********************************************************************
729 :
730 2 : szbuf=exc_size ! FIXME oversized!
731 2 : neh = nreh(1)
732 :
733 2 : if (nsppol==2) then
734 0 : ABI_WARNING("nsppol==2 is very experimental")
735 : end if
736 4 : if (ANY(nreh(1)/=nreh)) then
737 0 : ABI_ERROR(" different nreh are not supported")
738 : end if
739 :
740 6 : ABI_MALLOC_OR_DIE(cbuff_dpc,(exc_size), ierr)
741 : !
742 : ! The two cases nsppol==1,2 can be merged but we keep them
743 : ! separated to keep to code readable.
744 :
745 5 : SELECT CASE (toupper(block_type))
746 : CASE ("RESONANT")
747 :
748 1 : if (nsppol==1) then
749 : !
750 : ! Construct resonant block from the upper triangle stored on file.
751 385 : neh = nreh(1)
752 385 : do itp=1,neh
753 384 : read(funt, err=10, iomsg=errmsg) cbuff_dpc(1:itp)
754 73920 : do it=1,itp-1 ! The diagonal is treated below.
755 73536 : cttp = cbuff_dpc(it)
756 73536 : exc_ham(it ,itp) = cttp ! R
757 73536 : exc_ham(itp ,it ) = CONJG(cttp) ! row_sign R*
758 73536 : exc_ham(neh+it ,neh+itp) = row_sign*CONJG(cttp)
759 73920 : exc_ham(neh+itp,neh+it ) = row_sign*cttp
760 : end do
761 385 : if (diago_is_real) then
762 384 : exc_ham(itp ,itp) = DBLE(cbuff_dpc(itp))
763 384 : exc_ham(neh+itp,neh+itp) = row_sign*DBLE(cbuff_dpc(itp))
764 : else
765 0 : exc_ham(itp,itp) = cbuff_dpc(itp)
766 0 : exc_ham(neh+itp,neh+itp) = row_sign*(CONJG(cbuff_dpc(itp)))
767 : end if
768 : end do
769 : !
770 : !
771 : else
772 : ! FIXME this part won't work if we have a different number of e-h pairs
773 0 : ABI_CHECK(ALL(nreh==nreh(1)),"Different number of transitions")
774 : ! The file contains
775 : ! A) The up-up and the down-down block in packed form
776 : ! (only the upper triangle is stored since these blocks are Hermitian)
777 : ! B) The entire up-down exchange block (no symmetry here)
778 : !
779 : ! The resonant block is given by:
780 : ! | (v'c' up) | (v'c' dwn) |
781 : ! ------------------------------ where v_{-+} = v_{+-}^H when the momentum of the photon is neglected.
782 : ! | [T-W+v]++ | v+- | (vc up) Note that v_{+-} is not Hermitian due to the presence of different spins.
783 : ! R = ------------------------------ Actually it reduces to a Hermitian matrix when the system is not spin polarized.
784 : ! | v-+ | [T-W+v]-- | (vc dwn) [T-W+v] is Hermitian provided the the QP energies are purely real.
785 : ! ------------------------------
786 : !
787 : ! *) Fill the diagonal blocks.
788 : ! only the upper triangle is stored on file.
789 : ! row1,col1 refer to the resonant block.
790 : ! row2,col2 refer to the anti-resonant block.
791 0 : do block=1,2
792 0 : if (block==1) then
793 0 : spad=0
794 0 : spin_stride=SUM(nreh)
795 : end if
796 0 : if (block==2) then
797 0 : spad=nreh(1)
798 0 : spin_stride=2*nreh(1) + nreh(2)
799 : end if
800 0 : do itp=1,nreh(block)
801 0 : read(funt, err=10, iomsg=errmsg) cbuff_dpc(1:itp)
802 0 : col1 = itp+spad
803 0 : col2 = itp+spin_stride
804 0 : do it=1,itp-1
805 0 : cttp = cbuff_dpc(it)
806 0 : row1 = it + spad
807 0 : row2 = it + spin_stride
808 0 : exc_ham(row1,col1) = cttp ! [T-W+v]
809 0 : exc_ham(col1,row1) = CONJG(cttp)
810 0 : exc_ham(row2,col2) = row_sign*CONJG(cttp) ! row_sign [T-W+v]*
811 0 : exc_ham(col2,row2) = row_sign*cttp
812 : end do
813 0 : if (diago_is_real) then
814 0 : exc_ham(col1,col1) = DBLE(cbuff_dpc(itp))
815 0 : exc_ham(col2,col2) = row_sign*DBLE(cbuff_dpc(itp))
816 : else
817 0 : exc_ham(col1,col1) = cbuff_dpc(itp)
818 0 : exc_ham(col2,col2) = row_sign*CONJG(cbuff_dpc(itp))
819 : end if
820 : end do
821 : end do
822 : !
823 : ! Read v+- and reconstruct resonant and anti-resonat blocks.
824 : ! TODO recheck this
825 0 : pad_r1=SUM(nreh)
826 0 : pad_c1=2*nreh(1) + nreh(2)
827 :
828 0 : spin_dim=nreh(1)
829 0 : do itp=1,spin_dim
830 0 : read(funt, err=10, iomsg=errmsg) cbuff_dpc(1:spin_dim)
831 0 : exc_ham(1:spin_dim,nreh(1)+itp) = cbuff_dpc(1:spin_dim) ! upper reso
832 0 : exc_ham(1+pad_r1:pad_r1+spin_dim,pad_c1+itp) = row_sign*CONJG(cbuff_dpc(1:spin_dim)) ! upper anti-reso
833 0 : col1 = itp+nreh(1)
834 0 : col2 = itp+(2*nreh(1) + nreh(2))
835 0 : do it=1,spin_dim
836 0 : cttp = cbuff_dpc(it)
837 0 : row1 = it
838 0 : exc_ham(col1,row1) = CONJG(cttp) ! lower reso.
839 0 : row2 = it + SUM(nreh)
840 0 : exc_ham(col2,row2) = row_sign*cttp ! lower anti-reso.
841 : end do
842 : end do
843 : end if
844 :
845 : CASE ("COUPLING")
846 : !
847 1 : if (nsppol==1) then
848 : !
849 385 : do itp=1,neh
850 384 : read(funt, err=10, iomsg=errmsg) cbuff_dpc(1:itp)
851 74305 : do it=1,itp
852 73920 : cttp = cbuff_dpc(it)
853 73920 : exc_ham(it ,neh+itp) = cttp
854 73920 : exc_ham(itp ,neh+it ) = cttp
855 73920 : exc_ham(neh+it ,itp ) = row_sign*CONJG(cttp)
856 74304 : exc_ham(neh+itp,it ) = row_sign*CONJG(cttp)
857 : end do
858 : end do
859 : !
860 : else
861 : ! The coupling block is given by:
862 : ! | (c'v' up) | (c'v dwn) |
863 : ! ----------------------------------- where v_{-+} = v_{+-}^t when the momentum of the photon is neglected.
864 : ! | [-W+v]++ | v+- | (vc up) The entire matrix v_{+-} is stored on file.
865 : ! C = -----------------------------------
866 : ! | v-+ | [-W+v]-- | (vc dwn)
867 : ! -----------------------------------
868 : !
869 : ! *) Fill blocks that are diagonal in spin coordinates.
870 : ! row1,col1 refer to the resonat block.
871 : ! row2,col2 refer to the anti-resonant block.
872 0 : do block=1,2
873 0 : if (block==1) then
874 0 : spad_r=0
875 0 : spad_c=SUM(nreh)
876 : end if
877 0 : if (block==2) then
878 0 : spad_r=nreh(1)
879 0 : spad_c=2*nreh(1)+nreh(2)
880 : end if
881 0 : do itp=1,nreh(block)
882 0 : read(funt, err=10, iomsg=errmsg) cbuff_dpc(1:itp)
883 0 : col1 = itp+spad_c
884 0 : row2 = itp+spad_r
885 0 : do it=1,itp-1
886 0 : cttp = cbuff_dpc(it)
887 0 : row1 = it + spad_r
888 0 : col2 = it + spad_c
889 0 : exc_ham(row1,col1) = cttp ! upper coupling
890 0 : exc_ham(row2,col2) = cttp ! lower coupling (symmetric)
891 0 : exc_ham(col1,row1) = row_sign*CONJG(cttp) ! lower anti-coupling
892 0 : exc_ham(col2,row2) = row_sign*CONJG(cttp) ! upper anti-couling
893 : end do
894 : ! TODO recheck this
895 0 : row1 = itp+spad_r
896 0 : exc_ham(row1,col1) = cbuff_dpc(itp) ! Diagonals of the block.
897 0 : col2 = itp+spad_c
898 0 : exc_ham(col2,row2) = row_sign*CONJG(cbuff_dpc(itp))
899 : end do
900 : end do
901 : !
902 : ! Read Full v+- and reconstruct resonant and anti-resonat blocks.
903 : ! TODO recheck this
904 0 : spad=2*nreh(1) + nreh(2)
905 0 : pad_r1=SUM(nreh)
906 0 : pad_c1=2*nreh(1) + nreh(2)
907 :
908 0 : spin_dim=nreh(1)
909 0 : do itp=1,spin_dim
910 0 : read(funt, err=10, iomsg=errmsg) cbuff_dpc(1:spin_dim)
911 0 : exc_ham(1:spin_dim,spad+itp) = cbuff_dpc(1:spin_dim) ! upper block reso
912 0 : exc_ham(1+pad_r1:pad_r1+spin_dim,nreh(1)+itp) = row_sign*CONJG(cbuff_dpc(1:spin_dim)) ! upper block anti-reso
913 0 : col1 = itp+spad
914 0 : row2 = itp+nreh(1)
915 0 : do it=1,spin_dim
916 0 : cttp = cbuff_dpc(it)
917 0 : row1 = it
918 0 : exc_ham(col1,row1) = row_sign*CONJG(cttp) ! lower anti-reso.
919 0 : col2 = it + SUM(nreh)
920 0 : exc_ham(row2,col2) = cttp ! lower reso.
921 : end do
922 : end do
923 : !
924 : end if
925 :
926 : CASE DEFAULT
927 2 : ABI_ERROR("Unknown block_type: "//TRIM(block_type))
928 : END SELECT
929 :
930 2 : ABI_FREE(cbuff_dpc)
931 :
932 : return
933 :
934 : ! Handle IO Error
935 : 10 continue
936 0 : ABI_ERROR(errmsg)
937 :
938 0 : end subroutine exc_fullh_from_blocks
939 : !!***
940 :
941 : !----------------------------------------------------------------------
942 :
943 : !!****f* m_bse_io/rrs_of_glob
944 : !! NAME
945 : !! rrs_of_glob
946 : !!
947 : !! FUNCTION
948 : !! [+1,-1,0] if (row_glob,col_glob) belongs to the [ resonant, anti-resonant, (anti)coupling block ]
949 : !!
950 : !! INPUTS
951 : !!
952 : !! OUTPUT
953 : !!
954 : !! SOURCE
955 :
956 0 : pure function rrs_of_glob(row_glob,col_glob,size_glob)
957 :
958 : !Arguments ------------------------------------
959 : integer :: rrs_of_glob
960 : integer,intent(in) :: row_glob,col_glob
961 : integer,intent(in) :: size_glob(2)
962 :
963 : !Local variables ------------------------------
964 : integer :: nreh1,nreh2
965 : ! *************************************************************************
966 :
967 0 : nreh1=size_glob(1)/2 ! Matrix is square and nreh1==nreh2 but oh well.
968 0 : nreh2=size_glob(2)/2
969 :
970 0 : if (row_glob<=nreh1 .and. col_glob<=nreh2) then
971 : rrs_of_glob=+1 ! Resonant.
972 0 : else if (row_glob >nreh1 .and. col_glob >nreh2) then
973 : rrs_of_glob=-1 ! anti-Resonant.
974 : else
975 0 : rrs_of_glob=0
976 : end if
977 :
978 0 : end function rrs_of_glob
979 : !!***
980 :
981 : !----------------------------------------------------------------------
982 :
983 : !!****f* m_bse_io/ccs_of_glob
984 : !! NAME
985 : !!! ccs_of_glob
986 : !!
987 : !! FUNCTION
988 : !! [+1,-1,0] if (row_glob,col_glob) belongs to the [ coupling, anti-coupling, (anti)resonant block ]
989 : !!
990 : !! INPUTS
991 : !!
992 : !! OUTPUT
993 : !!
994 : !! SOURCE
995 :
996 0 : pure function ccs_of_glob(row_glob,col_glob,size_glob)
997 :
998 : !Arguments ------------------------------------
999 : integer :: ccs_of_glob
1000 : integer,intent(in) :: row_glob,col_glob
1001 : integer,intent(in) :: size_glob(2)
1002 :
1003 : !Local variables ------------------------------
1004 : integer :: nreh1,nreh2
1005 : ! *************************************************************************
1006 :
1007 0 : nreh1=size_glob(1)/2 ! Matrix is square and nreh1==nreh2 but oh well.
1008 0 : nreh2=size_glob(2)/2
1009 :
1010 0 : if (row_glob<=nreh1 .and. col_glob >nreh2) then ! Coupling.
1011 : ccs_of_glob = +1
1012 0 : else if (row_glob >nreh1 .and. col_glob<=nreh2) then ! anti-Coupling
1013 : ccs_of_glob = -1
1014 : else
1015 0 : ccs_of_glob = 0
1016 : end if
1017 :
1018 0 : end function ccs_of_glob
1019 : !!***
1020 :
1021 : !----------------------------------------------------------------------
1022 :
1023 : !!****f* m_bse_io/offset_in_file
1024 : !! NAME
1025 : !! offset_in_file
1026 : !!
1027 : !! FUNCTION
1028 : !! Return the offset of the matrix element (row_glob,col_glob)
1029 : !! size_glob(2) gives the number of row and column of the global matrix
1030 : !! nsblocks is the number of sublocks, used for nsppol==2 (not used if 1)
1031 : !! sub_block(2,2,nsblocks)= For each subblock the coordinates of the first and last element.
1032 : !!
1033 : !! INPUTS
1034 : !!
1035 : !! OUTPUT
1036 : !!
1037 : !! SOURCE
1038 :
1039 0 : function offset_in_file(row_glob,col_glob,size_glob,nsblocks,sub_block,bsize_elm,bsize_frm)
1040 :
1041 : !Arguments ------------------------------------
1042 : integer(XMPI_OFFSET_KIND) :: offset_in_file
1043 : integer,intent(in) :: row_glob,col_glob,nsblocks,bsize_elm,bsize_frm
1044 : integer,intent(in) :: size_glob(2),sub_block(2,2,nsblocks)
1045 :
1046 : !Local variables ------------------------------
1047 : !scalars
1048 : integer :: ii,jj,ijp_glob,swap
1049 : integer(XMPI_OFFSET_KIND) :: my_offset
1050 : ! *************************************************************************
1051 :
1052 0 : if (nsblocks==1) then
1053 0 : ii = row_glob
1054 0 : jj = col_glob
1055 0 : if (ii>size_glob(1)/2) ii = ii - size_glob(1)/2 ! Wrap the index.
1056 0 : if (jj>size_glob(2)/2) jj = jj - size_glob(2)/2
1057 0 : if (jj<ii) then ! Exchange the indices since the symmetric element is read.
1058 0 : swap = jj
1059 0 : jj = ii
1060 0 : ii = swap
1061 : end if
1062 0 : ijp_glob = ii + jj*(jj-1)/2 ! Index for packed storage mode.
1063 0 : my_offset = (ijp_glob-1)*bsize_elm + (jj-1)*2*bsize_frm
1064 : else
1065 : ABI_UNUSED(sub_block(1,1,1))
1066 0 : ABI_ERROR("nsppol==2 not coded")
1067 : end if
1068 :
1069 0 : offset_in_file = my_offset
1070 :
1071 0 : end function offset_in_file
1072 : !!***
1073 :
1074 : !----------------------------------------------------------------------
1075 :
1076 : !!****f* m_bse_io/exc_read_rblock_fio
1077 : !! NAME
1078 : !! exc_read_rblock_fio
1079 : !!
1080 : !! FUNCTION
1081 : !! Reads the resonant block from file using Fortran IO.
1082 : !!
1083 : !! INPUTS
1084 : !! funt=Fortran unit number.
1085 : !! nsppol=Number of spins
1086 : !! exc_size=Size of the resonant bock.
1087 : !! diago_is_real=.TRUE. if diagonal elements are real.
1088 : !! nreh(nsppol)=Number of resonant transitions for each spin.
1089 : !!
1090 : !! OUTPUT
1091 : !! ierr=Status error
1092 : !! exc_mat(exc_size,exc_size)=The resonant block.
1093 : !!
1094 : !! SOURCE
1095 :
1096 19 : subroutine exc_read_rblock_fio(funt,diago_is_real,nsppol,nreh,exc_size,exc_mat,ierr)
1097 :
1098 : !Arguments ------------------------------------
1099 : !scalars
1100 : integer,intent(in) :: funt,nsppol,exc_size
1101 : logical,intent(in) :: diago_is_real
1102 : integer,intent(out) :: ierr
1103 : !arrays
1104 : integer,intent(in) :: nreh(nsppol)
1105 : complex(dp),intent(out) :: exc_mat(exc_size,exc_size)
1106 :
1107 : !Local variables ------------------------------
1108 : !scalars
1109 : integer :: itp,it,block,col,row,spad
1110 : complex(dp) :: ctemp
1111 : character(len=500) :: errmsg
1112 : !arrays
1113 19 : complex(dp),allocatable :: cbuff_dpc(:)
1114 : ! *************************************************************************
1115 :
1116 19 : ierr=0
1117 :
1118 : ! Construct full resonant block using Hermiticity. File is always in double precision.
1119 57 : ABI_MALLOC(cbuff_dpc,(exc_size))
1120 :
1121 19 : if (nsppol==1) then ! Construct resonant block from the upper triangle stored on file.
1122 : !
1123 2801 : do itp=1,nreh(1)
1124 2784 : read(funt, err=10, iomsg=errmsg) cbuff_dpc(1:itp)
1125 356208 : do it=1,itp-1 ! Diagonal is treated below.
1126 353424 : ctemp = cbuff_dpc(it)
1127 353424 : exc_mat(it,itp) = ctemp
1128 356208 : exc_mat(itp,it) = CONJG(ctemp)
1129 : end do
1130 2801 : if (diago_is_real) then
1131 2784 : exc_mat(itp,itp) = DBLE(cbuff_dpc(itp))
1132 : else
1133 0 : exc_mat(itp,itp) = cbuff_dpc(itp)
1134 : end if
1135 : end do
1136 : !
1137 : else
1138 : ! The file contains
1139 : ! A) The up-up and the down-down block in packed form
1140 : ! (only the upper triangle is stored since these blocks are Hermitian)
1141 : ! B) The entire up-down exchange block (no symmetry here)
1142 : !
1143 : ! A) Construct resonant blocks from the upper triangles stored on file.
1144 : ! FIXME this part won't work if we have a different number of e-h pairs
1145 : !ABI_CHECK(ALL(nreh==nreh(1)),"Different number of transitions")
1146 :
1147 6 : do block=1,2
1148 4 : if (block==1) spad=0
1149 2 : if (block==2) spad=nreh(1)
1150 1542 : do itp=1,nreh(block)
1151 1536 : read(funt, err=10, iomsg=errmsg) cbuff_dpc(1:itp)
1152 1536 : col = itp+spad
1153 295680 : do it=1,itp-1 ! diagonal is treated below.
1154 294144 : ctemp = cbuff_dpc(it)
1155 294144 : row = it + spad
1156 294144 : exc_mat(row,col) = ctemp
1157 295680 : exc_mat(col,row) = CONJG(ctemp)
1158 : end do
1159 1540 : if (diago_is_real) then
1160 1536 : exc_mat(col,col) = DBLE(cbuff_dpc(itp))
1161 : else
1162 0 : exc_mat(col,col) = cbuff_dpc(itp)
1163 : end if
1164 : end do
1165 : end do
1166 : !
1167 : ! read v+- that is a matrix with shape nreh(1) X nreh(2)
1168 2 : spad=nreh(1)
1169 770 : do itp=1,nreh(2)
1170 768 : read(funt, err=10, iomsg=errmsg) cbuff_dpc(1:nreh(1))
1171 295682 : exc_mat(1:nreh(1),spad+itp) = cbuff_dpc(1:nreh(1))
1172 : end do
1173 :
1174 : end if
1175 :
1176 19 : ABI_FREE(cbuff_dpc)
1177 :
1178 : return
1179 :
1180 : ! Raise the error.
1181 : 10 continue
1182 0 : ierr = 1
1183 0 : ABI_WARNING(errmsg)
1184 :
1185 0 : end subroutine exc_read_rblock_fio
1186 : !!***
1187 :
1188 : !----------------------------------------------------------------------
1189 :
1190 : !!****f* m_bse_io/exc_amplitude
1191 : !! NAME
1192 : !! exc_amplitude
1193 : !!
1194 : !! FUNCTION
1195 : !! Calculate the amplitude function of the excitonic eigenstate |exc_vec\>
1196 : !! F(w) = \sum_t |<t|exc_vec>|^2 \delta(ww- ene_t) where the sum over t is done
1197 : !! of the full set of transitions used to construct the BS Hamiltoniam.
1198 : !!
1199 : !! INPUTS
1200 : !! Bsp<excparam>=Structure storing the parameters of the run.
1201 : !! eig_fname=The name of the file storing the excitonic eigenvectors.
1202 : !! nvec=Number of excitonic states to analyze.
1203 : !! vec_idx(nvec)=List with the indices of the excitonic states sorted in ascending order.
1204 : !! out_fname=The name of the file where the results are written.
1205 : !!
1206 : !! OUTPUT
1207 : !! Only writing.
1208 : !!
1209 : !! SOURCE
1210 :
1211 0 : subroutine exc_amplitude(Bsp,eig_fname,nvec,vec_idx,out_fname)
1212 :
1213 : !Arguments ------------------------------------
1214 : !scalars
1215 : integer,intent(in) :: nvec
1216 : character(len=*),intent(in) :: eig_fname,out_fname
1217 : type(excparam),intent(in) :: BSp
1218 : ! arrays
1219 : integer,intent(in) :: vec_idx(nvec)
1220 :
1221 : !Local variables ------------------------------
1222 : !scalars
1223 : integer :: vec,art_idx,ierr
1224 : integer :: spin,iw,it,nw,pos_w,neg_w,out_unt,rt_idx,hsize
1225 : real(dp) :: ene_rt,ampl_eh,ampl_he,xx,stdev,w_max,w_min,step
1226 : character(len=500) :: msg
1227 : !arrays
1228 0 : real(dp),allocatable :: wmesh(:),amplitude(:),ene_list(:)
1229 0 : complex(dp),allocatable :: vec_list(:,:)
1230 : ! *************************************************************************
1231 :
1232 : ! Setup of the frequency mesh for F(w).
1233 0 : w_min=greatest_real; w_max=smallest_real
1234 0 : do spin=1,BSp%nsppol
1235 0 : do it=1,BSp%nreh(spin)
1236 0 : ene_rt = Bsp%Trans(it,spin)%en
1237 0 : w_min = MIN(w_min,ene_rt)
1238 0 : w_max = MAX(w_max,ene_rt)
1239 : end do
1240 : end do
1241 :
1242 0 : step = Bsp%domega
1243 0 : if (Bsp%use_coupling==0) then
1244 0 : nw = (w_max - w_min)/step + 1
1245 0 : ABI_MALLOC(wmesh,(nw))
1246 0 : wmesh = arth(w_min,step,nw)
1247 : else
1248 : ! Both positive and negative frequencies are needed.
1249 0 : pos_w = (w_max - w_min)/step + 1
1250 0 : neg_w = pos_w; if (ABS(w_min) < tol6) neg_w=neg_w-1 ! zero should not included twice.
1251 0 : nw = neg_w + pos_w
1252 0 : ABI_MALLOC(wmesh,(nw))
1253 0 : wmesh(1:neg_w) = arth(-w_max,step,neg_w)
1254 0 : wmesh(neg_w+1:) = arth( w_min,step,pos_w)
1255 : end if
1256 : !
1257 : ! Read selected eigenvectors.
1258 0 : hsize = SUM(Bsp%nreh); if (Bsp%use_coupling>0) hsize=2*hsize
1259 :
1260 0 : ABI_MALLOC(ene_list,(nvec))
1261 0 : ABI_MALLOC_OR_DIE(vec_list,(hsize,nvec), ierr)
1262 :
1263 0 : call exc_read_eigen(eig_fname,hsize,nvec,vec_idx,vec_list,ene_list,Bsp=Bsp)
1264 :
1265 0 : ABI_FREE(ene_list)
1266 :
1267 0 : if (open_file(out_fname,msg,newunit=out_unt,form="formatted",action="write") /= 0) then
1268 0 : ABI_ERROR(msg)
1269 : end if
1270 :
1271 0 : write(out_unt,*)"# Amplitude functions F(w) = \sum_t |<t|exc_vec>|^2 \delta(ww- ene_t), w is given in eV. "
1272 0 : write(out_unt,*)"# Number of excitonic eigenvectors analyzed: ",nvec
1273 :
1274 0 : ABI_MALLOC(amplitude,(nw))
1275 0 : stdev = BSp%broad ! Broadening for the gaussian.
1276 :
1277 0 : do vec=1,nvec
1278 : !
1279 : ! amplitude(ww) = \sum_t |<t|exc_vec>|^2 \delta(ww- ene_t)
1280 0 : amplitude = zero
1281 0 : do spin=1,BSp%nsppol
1282 0 : do it=1,BSp%nreh(spin)
1283 0 : ene_rt = Bsp%Trans(it,spin)%en
1284 0 : rt_idx = it + (spin-1)*Bsp%nreh(1)
1285 0 : ampl_eh = (ABS(vec_list(rt_idx,vec)))**2
1286 0 : if (Bsp%use_coupling>0) then ! Need the index and the amplitude of the anti-resonant component.
1287 0 : art_idx = it + (spin-1)*Bsp%nreh(1) + SUM(Bsp%nreh)
1288 0 : ampl_he = (ABS(vec_list(art_idx,vec)))**2
1289 : end if
1290 : !
1291 0 : do iw=1,nw ! Accumulate
1292 0 : xx = wmesh(iw) - ene_rt
1293 0 : amplitude(iw) = amplitude(iw) + ampl_eh * gaussian(xx, stdev)
1294 0 : if (Bsp%use_coupling>0) then
1295 0 : xx = wmesh(iw) + ene_rt
1296 0 : amplitude(iw) = amplitude(iw) + ampl_he * gaussian(xx, stdev)
1297 : end if
1298 : end do
1299 : !
1300 : end do
1301 : end do
1302 : !
1303 : ! Write results
1304 0 : write(out_unt,*)"# Amplitude function F(w) for exc_vector index ",vec_idx(vec)
1305 0 : do iw=1,nw
1306 0 : write(out_unt,*)wmesh(iw)*Ha_eV,amplitude(iw)
1307 : end do
1308 0 : write(out_unt,*)"#"
1309 : end do
1310 :
1311 0 : close(out_unt)
1312 :
1313 0 : ABI_FREE(amplitude)
1314 0 : ABI_FREE(wmesh)
1315 0 : ABI_FREE(vec_list)
1316 :
1317 0 : end subroutine exc_amplitude
1318 : !!***
1319 :
1320 : !----------------------------------------------------------------------
1321 :
1322 : !!****f* m_bse_io/exc_write_optme
1323 : !! NAME
1324 : !! exc_write_optme
1325 : !!
1326 : !! FUNCTION
1327 : !! Writes the optical matrix elements in the OME.nc file.
1328 : !! Note that this is only available when NetCDF is available
1329 : !!
1330 : !! INPUTS
1331 : !! filename=filename used to write the optical matrix elements
1332 : !! minb,maxb=minimum and max band index that have been calculated.
1333 : !! nkbz=Number of points in the full Brillouin zone.
1334 : !! nsppol=Number of independent spin polarizations.
1335 : !! nq=Number of "small" q for optical limit.
1336 : !! opt_cvk=Optical matrix elements to be written
1337 : !!
1338 : !! OUTPUT
1339 : !! ierr=return status of the writing process.
1340 : !! => 0 if everything was ok
1341 : !! => -1 if NetCDF is not available
1342 : !! => 1 if NetCDF is available
1343 : !!
1344 : !! SOURCE
1345 :
1346 0 : subroutine exc_write_optme(filename,minb,maxb,nkbz,nsppol,nq,opt_cvk,ierr)
1347 :
1348 : !Arguments ------------------------------------
1349 : integer,intent(in) :: minb,maxb,nkbz,nsppol,nq
1350 : character(len=fnlen),intent(in) :: filename
1351 : complex(dp),intent(in) :: opt_cvk(minb:maxb,minb:maxb,nkbz,nsppol,nq)
1352 : integer,intent(out) :: ierr
1353 :
1354 : !Local variables ------------------------------
1355 : !scalars
1356 : integer :: ncid,cmplx_id,nband_id,nkbz_id,nsppol_id,nq_id,xyz_id
1357 : integer :: minb_id,maxb_id,ome_id,iq,is,ik,ib,jb
1358 : integer :: dimOME(6),dimKPT(2),dimQPT(2),dimSCA(0),start6(6),count6(6)
1359 : real(dp) :: complex2(2)
1360 : ! *************************************************************************
1361 :
1362 0 : ierr = 1
1363 :
1364 : !1. Create netCDF file
1365 0 : NCF_CHECK_MSG(nctk_open_create(ncid, filename, xmpi_comm_self), " create netcdf OME file")
1366 :
1367 : !2. Define dimensions
1368 0 : NCF_CHECK(nf90_def_dim(ncid,"xyz",3,xyz_id))
1369 0 : NCF_CHECK(nf90_def_dim(ncid,"cmplx",2,cmplx_id))
1370 0 : NCF_CHECK(nf90_def_dim(ncid,"nband",maxb-minb+1,nband_id))
1371 0 : NCF_CHECK(nf90_def_dim(ncid,"nkbz",nkbz,nkbz_id))
1372 0 : NCF_CHECK(nf90_def_dim(ncid,"nq",nq,nq_id))
1373 0 : NCF_CHECK(nf90_def_dim(ncid,"nsppol",nsppol,nsppol_id))
1374 :
1375 : !Dimensions for optical matrix elements
1376 0 : dimOME = (/ cmplx_id, nband_id, nband_id, nkbz_id, nsppol_id, nq_id /)
1377 : !Dimensions for kpoint positions
1378 0 : dimKPT = (/ xyz_id, nkbz_id /)
1379 : !Dimensions for qpoints for the optical limit
1380 0 : dimQPT = (/ xyz_id, nq_id /)
1381 :
1382 : !3. Define variables
1383 :
1384 : call ab_define_var(ncid, dimOME, ome_id, NF90_DOUBLE,&
1385 0 : "OME", "Values of optical matrix elements","Tobedone")
1386 :
1387 : call ab_define_var(ncid, dimSCA, minb_id, NF90_INT,"minb",&
1388 0 : "Minimum band index for the optical matrix elements", "Dimensionless")
1389 :
1390 : call ab_define_var(ncid, dimSCA, maxb_id, NF90_INT,"maxb",&
1391 0 : "Maximum band index for the optical matrix elements", "Dimensionless")
1392 :
1393 : !4. End define mode
1394 0 : NCF_CHECK(nf90_enddef(ncid))
1395 :
1396 : !5 Write scalars (minb and maxb)
1397 :
1398 0 : NCF_CHECK(nf90_put_var(ncid, minb_id, minb))
1399 0 : NCF_CHECK(nf90_put_var(ncid, maxb_id, maxb))
1400 :
1401 : !6 Write optical matrix elements
1402 :
1403 0 : do iq=1,nq
1404 0 : do is=1,nsppol
1405 0 : do ik=1,nkbz
1406 0 : do ib=minb,maxb
1407 0 : do jb=minb,maxb
1408 0 : start6 = (/ 1, ib-minb+1, jb-minb+1, ik, is, iq /)
1409 0 : count6 = (/ 2, 1, 1, 1, 1, 1 /)
1410 0 : complex2 = (/ REAL(opt_cvk(ib,jb,ik,is,iq)),AIMAG(opt_cvk(ib,jb,ik,is,iq)) /)
1411 0 : NCF_CHECK(nf90_put_var(ncid, ome_id, complex2, start = start6, count = count6))
1412 : end do
1413 : end do
1414 : end do
1415 : end do
1416 : end do
1417 :
1418 : !7 Close file
1419 0 : NCF_CHECK(nf90_close(ncid))
1420 0 : ierr = 0
1421 :
1422 0 : end subroutine exc_write_optme
1423 : !!***
1424 :
1425 : !----------------------------------------------------------------------
1426 :
1427 : !!****f* m_hexc/exc_ham_ncwrite
1428 : !! NAME
1429 : !! exc_ham_ncwrite
1430 : !!
1431 : !! FUNCTION
1432 : !! Writes the content of a hexc object to a NETCDF file
1433 : !! according to the ETSF-IO specifications.
1434 : !!
1435 : !! INPUTS
1436 : !! ncid =NC file handle
1437 : !!
1438 : !! OUTPUT
1439 : !!
1440 : !! SOURCE
1441 :
1442 0 : subroutine exc_ham_ncwrite(ncid,Kmesh,BSp,hsize,nreh,vcks2t,hreso,diag)
1443 :
1444 : !Arguments ------------------------------------
1445 : !scalars
1446 : integer,intent(in) :: ncid
1447 : integer,intent(in) :: hsize
1448 : type(kmesh_t),intent(in) :: Kmesh
1449 : type(excparam),intent(in) :: BSp
1450 : integer,intent(in) :: nreh(BSp%nsppol)
1451 : integer,target,intent(in) :: vcks2t(BSp%maxnbndv,BSp%maxnbndc,Kmesh%nbz,BSp%nsppol)
1452 : complex(dp),target,intent(in) :: hreso(hsize,hsize)
1453 : complex(dp),target,intent(in) :: diag(hsize)
1454 :
1455 : !Local variables-------------------------------
1456 : integer :: ncerr
1457 : integer :: max_nreh, sum_nreh
1458 0 : real(dp), contiguous, pointer :: r2vals(:,:),r3vals(:,:,:)
1459 : ! *************************************************************************
1460 :
1461 : ! ==============================================
1462 : ! === Write the dimensions specified by ETSF ===
1463 : ! ==============================================
1464 0 : max_nreh = MAXVAL(nreh)
1465 0 : sum_nreh = SUM(nreh)
1466 :
1467 : ncerr = nctk_def_dims(ncid, [nctkdim_t("number_of_reduced_dimensions", 3), nctkdim_t("number_of_spins", bsp%nsppol),&
1468 : nctkdim_t("number_of_kpoints", kmesh%nbz), nctkdim_t("max_number_of_valence_bands", bsp%maxnbndv),&
1469 : nctkdim_t("max_number_of_conduction_bands", bsp%maxnbndc), nctkdim_t("max_number_of_transitions", max_nreh),&
1470 0 : nctkdim_t("total_number_of_transitions", sum_nreh), nctkdim_t("cplex", 2)], defmode=.True.)
1471 0 : NCF_CHECK(ncerr)
1472 :
1473 : ncerr = nctk_def_arrays(ncid, [&
1474 : nctkarr_t('vcks2t', "i", 'max_number_of_valence_bands, max_number_of_conduction_bands, number_of_kpoints, number_of_spins'),&
1475 : nctkarr_t('hamiltonian', "dp", 'cplex total_number_of_transitions total_number_of_transitions'),&
1476 : nctkarr_t('diagonal', "dp", 'cplex total_number_of_transitions'),&
1477 : nctkarr_t('lomo', "dp", 'number_of_spins'), &
1478 : nctkarr_t('homo', "dp", 'number_of_spins'), &
1479 : nctkarr_t('lumo', "dp", 'number_of_spins'), &
1480 : nctkarr_t('humo', "dp", 'number_of_spins'), &
1481 : nctkarr_t("reduced_coordinates_of_kpoints", "dp", "number_of_reduced_dimensions, number_of_kpoints"), &
1482 : nctkarr_t("kpoint_weights", "dp", "number_of_kpoints") &
1483 0 : ])
1484 0 : NCF_CHECK(ncerr)
1485 :
1486 : ! Write data
1487 0 : NCF_CHECK(nctk_set_datamode(ncid))
1488 0 : NCF_CHECK(nf90_put_var(ncid, vid('vcks2t'), vcks2t))
1489 0 : NCF_CHECK(nf90_put_var(ncid, vid('lomo'), BSp%lomo_spin))
1490 0 : NCF_CHECK(nf90_put_var(ncid, vid('homo'), BSp%homo_spin))
1491 0 : NCF_CHECK(nf90_put_var(ncid, vid('lumo'), BSp%lumo_spin))
1492 0 : NCF_CHECK(nf90_put_var(ncid, vid('humo'), BSp%humo_spin))
1493 :
1494 0 : call c_f_pointer(c_loc(hreso(1,1)), r3vals, shape=[2, hsize, hsize])
1495 0 : NCF_CHECK(nf90_put_var(ncid, vid("hamiltonian"), r3vals))
1496 :
1497 0 : call c_f_pointer(c_loc(diag(1)), r2vals, shape=[2, hsize])
1498 0 : NCF_CHECK(nf90_put_var(ncid, vid("diag"), r2vals))
1499 :
1500 : ! Write K-points
1501 0 : NCF_CHECK(nf90_put_var(ncid, vid("reduced_coordinates_of_kpoints"), kmesh%bz))
1502 0 : NCF_CHECK(nf90_put_var(ncid, vid("kpoint_weights"), kmesh%wt))
1503 0 : NCF_CHECK(ncerr)
1504 :
1505 : contains
1506 0 : integer function vid(vname)
1507 : character(len=*),intent(in) :: vname
1508 0 : vid = nctk_idname(ncid, vname)
1509 : end function vid
1510 :
1511 : end subroutine exc_ham_ncwrite
1512 : !!***
1513 :
1514 : !----------------------------------------------------------------------
1515 :
1516 : END MODULE m_bse_io
1517 : !!***
|
