Branch data 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 : : !!***
|
