Branch data Line data Source code
1 : : !!****m* ABINIT/m_exc_diago
2 : : !! NAME
3 : : !! m_exc_diago
4 : : !!
5 : : !! FUNCTION
6 : : !!
7 : : !! COPYRIGHT
8 : : !! Copyright (C) 2009-2026 ABINIT and EXC groups (L.Reining, V.Olevano, F.Sottile, S.Albrecht, G.Onida, M.Giantomassi)
9 : : !! This file is distributed under the terms of the
10 : : !! GNU General Public License, see ~abinit/COPYING
11 : : !! or http://www.gnu.org/copyleft/gpl.txt .
12 : : !!
13 : : !! SOURCE
14 : :
15 : : #if defined HAVE_CONFIG_H
16 : : #include "config.h"
17 : : #endif
18 : :
19 : : #include "abi_common.h"
20 : :
21 : : MODULE m_exc_diago
22 : :
23 : : use defs_basis
24 : : use m_bs_defs
25 : : use m_abicore
26 : : use m_errors
27 : : USE_MPI
28 : : use m_xmpi
29 : : use m_hdr
30 : : use m_sort
31 : : use m_slk
32 : :
33 : : use defs_datatypes, only : pseudopotential_type
34 : : use m_io_tools, only : open_file
35 : : use m_fstrings, only : int2char4
36 : : use m_numeric_tools, only : print_arr, hermitianize
37 : : use m_crystal, only : crystal_t
38 : : use m_kpts, only : listkk
39 : : use m_bz_mesh, only : kmesh_t
40 : : use m_ebands, only : ebands_t
41 : : use m_eprenorms, only : eprenorms_t
42 : : use m_wfd, only : wfdgw_t
43 : : use m_paw_hr, only : pawhur_t
44 : : use m_pawtab, only : pawtab_type
45 : : use m_exc_itdiago, only : exc_iterative_diago
46 : : use m_hide_lapack, only : xheev, xheevx, xgeev, xhegvx, xginv, xhdp_invert, xhegv
47 : : use m_hide_blas, only : xdotc, xgemm
48 : : use m_bse_io, only : exc_fullh_from_blocks, offset_in_file, rrs_of_glob, ccs_of_glob, &
49 : : & exc_read_bshdr, exc_skip_bshdr_mpio, exc_read_rblock_fio
50 : : use m_exc_spectra, only : build_spectra
51 : :
52 : : implicit none
53 : :
54 : : private
55 : :
56 : : #if defined HAVE_MPI1
57 : : include 'mpif.h'
58 : : #endif
59 : :
60 : : !#define DEV_MG_DEBUG_THIS
61 : :
62 : : public :: exc_diago_driver ! Driver routine for the direct diagonalization of the BSE Hamiltonian (main entry point)
63 : : !!***
64 : :
65 : : contains
66 : :
67 : : !!****f* m_exc_diago/exc_diago_driver
68 : : !! NAME
69 : : !! exc_diago_driver
70 : : !!
71 : : !! FUNCTION
72 : : !! Driver routine for the direct diagonalization of the Hermitian excitonic Hamiltonian.
73 : : !!
74 : : !! INPUTS
75 : : !! neh=Rank of the resonant block of the Hamiltonian.
76 : : !! BS_files<excfiles>=Datatype storing names and files used in the Bethe-Salpeter code.
77 : : !! %exh=Name of the file storing the excitonic resonant part.
78 : : !!
79 : : !! OUTPUT
80 : : !! Eigenvalues and eigenvectors are written on file.
81 : : !!
82 : : !! SOURCE
83 : :
84 : 6 : subroutine exc_diago_driver(Wfd,Bsp,BS_files,KS_BSt,QP_BSt,Cryst,Kmesh,Psps,&
85 : 6 : & Pawtab,Hur,Hdr_bse,drude_plsmf,Epren)
86 : :
87 : : !Arguments ------------------------------------
88 : : !scalars
89 : : real(dp),intent(in) :: drude_plsmf
90 : : type(excparam),intent(in) :: BSp
91 : : type(excfiles),intent(in) :: BS_files
92 : : type(Hdr_type),intent(in) :: Hdr_bse
93 : : type(crystal_t),intent(in) :: Cryst
94 : : type(pseudopotential_type),intent(in) :: Psps
95 : : type(kmesh_t),intent(in) :: Kmesh
96 : : type(ebands_t),intent(in) :: KS_BSt,QP_BSt
97 : : type(wfdgw_t),intent(inout) :: Wfd
98 : : type(eprenorms_t),intent(in) :: Epren
99 : : !arrays
100 : : type(pawtab_type),intent(in) :: Pawtab(Cryst%ntypat*Wfd%usepaw)
101 : : type(pawhur_t),intent(in) :: Hur(Cryst%natom*Wfd%usepaw)
102 : :
103 : : !Local variables ------------------------------
104 : : !scalars
105 : : integer :: my_rank,master,comm,prtvol
106 : : complex(dp) :: exc_gap,gw_gap
107 : : logical :: eval_eigenstates
108 : : character(len=500) :: msg
109 : : !arrays
110 : 12 : real(dp) :: gaps(3,QP_BSt%nsppol)
111 : : !************************************************************************
112 : :
113 : : DBG_ENTER("COLL")
114 : :
115 : 6 : comm = Wfd%comm
116 : 6 : my_rank = Wfd%my_rank
117 : 6 : master = Wfd%master
118 : 6 : prtvol = Wfd%prtvol
119 : :
120 [ - + ]: 6 : if (BSp%have_complex_ene) then
121 : 0 : ABI_ERROR("Complex energies are not supported yet")
122 : : end if
123 : : !
124 : : ! This trick is needed to restart a CG run, use DDIAGO to calculate the spectra reusing an old BSEIG file.
125 [ - + - - ]: 6 : eval_eigenstates = (BS_files%in_eig == BSE_NOFILE) .or. (Bsp%algorithm == BSE_ALGO_CG)
126 : :
127 : : if (eval_eigenstates) then
128 : : !
129 : 11 : select case (BSp%algorithm)
130 : : case (BSE_ALGO_DDIAGO)
131 [ + + ]: 5 : if (BSp%use_coupling==0) then
132 : 4 : call exc_diago_resonant(BSp,BS_files,Hdr_bse,prtvol,comm)
133 [ + + ]: 4 : if(Bsp%do_ep_renorm) then
134 : 1 : call exc_diago_resonant(BSp,BS_files,Hdr_bse,prtvol,comm,Epren=Epren,Kmesh=Kmesh,Cryst=Cryst,elph_lifetime=.TRUE.)
135 : : end if
136 : : else
137 [ - + ]: 1 : if (Bsp%have_complex_ene) then
138 : : ! Solve Hv = ev with generic complex matrix.
139 : 0 : call exc_diago_coupling(BSp,BS_files,Hdr_bse,prtvol,comm)
140 : : else
141 : : ! Solve generalized eigenvalue problem F Hbar with Hbar Hermitian definitive positive matrix.
142 : 1 : call exc_diago_coupling_hegv(BSp,BS_files,Hdr_bse,prtvol,comm)
143 : : end if
144 : : end if
145 : :
146 : : case (BSE_ALGO_CG)
147 [ + - ]: 1 : if (BSp%use_coupling==0) then
148 : 1 : call exc_iterative_diago(Bsp,BS_files,Hdr_bse,prtvol,comm)
149 : : else
150 : 0 : ABI_ERROR("CG + coupling not coded")
151 : : end if
152 : :
153 : : case default
154 : 0 : write(msg,'(a,i0)')" Wrong value for Bsp%algorithm: ",Bsp%algorithm
155 [ + + - ]: 6 : ABI_ERROR(msg)
156 : : end select
157 : : !
158 [ + - ]: 6 : if (my_rank==master) then
159 : 6 : call QP_BSt%report_gap(header="QP bands",unit=std_out,gaps=gaps)
160 [ + - - + : 19 : gw_gap = MINVAL(gaps(2,:))
- - + + ]
161 : 6 : call exc_print_eig(BSp,BS_files%out_eig,gw_gap,exc_gap)
162 : : end if
163 : 6 : call xmpi_barrier(comm)
164 : : !
165 : : end if
166 : :
167 : 6 : call build_spectra(BSp,BS_files,Cryst,Kmesh,KS_BSt,QP_BSt,Psps,Pawtab,Wfd,Hur,drude_plsmf,comm)
168 : :
169 : : ! Electron-phonon renormalization !
170 [ + + + + : 6 : if (BSp%algorithm == BSE_ALGO_DDIAGO .and. BSp%use_coupling == 0 .and. BSp%do_ep_renorm) then
+ + ]
171 : 1 : call build_spectra(BSp,BS_files,Cryst,Kmesh,KS_BSt,QP_BSt,Psps,Pawtab,Wfd,Hur,drude_plsmf,comm,Epren=Epren)
172 : : end if
173 : :
174 : : DBG_EXIT("COLL")
175 : :
176 : 6 : end subroutine exc_diago_driver
177 : : !!***
178 : :
179 : : !----------------------------------------------------------------------
180 : :
181 : : !!****f* m_exc_diago/exc_diago_resonant
182 : : !! NAME
183 : : !! exc_diago_resonant
184 : : !!
185 : : !! FUNCTION
186 : : !! Calculates eigenvalues and eigenvectors of the Hermitian excitonic Hamiltonian (coupling is neglected).
187 : : !!
188 : : !! INPUTS
189 : : !! Bsp
190 : : !! BS_files<excfiles>=Datatype storing names and files used in the Bethe-Salpeter code.
191 : : !! comm=MPI communicator.
192 : : !! bseig_fname=The name of the output file
193 : : !! prtvol=Verbosity level.
194 : : !!
195 : : !! OUTPUT
196 : : !! Eigenvalues and eigenvectors are written on file bseig_fname
197 : : !!
198 : : !! SOURCE
199 : :
200 : 5 : subroutine exc_diago_resonant(Bsp,BS_files,Hdr_bse,prtvol,comm,Epren,Kmesh,Cryst,elph_lifetime)
201 : :
202 : : !Arguments ------------------------------------
203 : : !scalars
204 : : integer,intent(in) :: comm,prtvol
205 : : logical,optional,intent(in) :: elph_lifetime
206 : : type(excparam),intent(in) :: BSp
207 : : type(excfiles),intent(in) :: BS_files
208 : : type(Hdr_type),intent(in) :: Hdr_bse
209 : : type(eprenorms_t),optional,intent(in) :: Epren
210 : : type(kmesh_t),optional,intent(in) :: Kmesh
211 : : type(crystal_t),optional,intent(in) :: Cryst
212 : :
213 : : !Local variables ------------------------------
214 : : !scalars
215 : : integer,parameter :: master=0
216 : : integer :: ii,it,mi,hreso_unt,eig_unt,exc_size,neh1,neh2,j
217 : : integer :: nsppol,il,iu,mene_found,nstates
218 : : integer :: nprocs,my_rank,fform,nene_printed,ierr
219 : : real(dp) :: exc_gap,exc_maxene,abstol
220 : : real(dp) :: vl,vu
221 : : logical :: use_scalapack,do_full_diago,diagonal_is_real
222 : : character(len=500) :: msg
223 : : character(len=fnlen) :: hreso_fname,bseig_fname
224 : : !arrays
225 : 5 : real(dp),allocatable :: exc_ene(:)
226 : 5 : complex(dp),allocatable :: exc_mat(:,:),exc_vec(:,:)
227 : : #if defined HAVE_LINALG_SCALAPACK && defined HAVE_MPI_IO
228 : : integer :: amode,mpi_fh,istwf_k,tbloc
229 : : integer(XMPI_OFFSET_KIND) :: ehdr_offset,fmarker
230 : : integer :: block_sizes(2,3),array_of_sizes(2),gsub(2,2)
231 : : logical,parameter :: is_fortran_file=.TRUE.
232 : : real(dp),external :: PDLAMCH
233 : : type(slkmat_dp_t) :: Slk_mat,Slk_vec
234 : : type(slk_processor_t) :: Slk_processor
235 : : #endif
236 : :
237 : : integer :: ik, ic, iv, isppol, ireh, ep_ik, itemp
238 : : complex(dp) :: en
239 : :
240 : : real(dp) :: dksqmax
241 : 5 : integer,allocatable :: bs2eph(:,:)
242 : : integer :: sppoldbl, timrev
243 : : logical :: do_ep_renorm, do_ep_lifetime
244 : : integer :: ntemp
245 : : character(len=4) :: ts
246 : 5 : complex(dp),allocatable :: exc_vl(:,:),exc_ene_c(:)
247 : : complex(dp) :: ctemp
248 : : !! complex(dp),allocatable :: ovlp(:,:)
249 : : !************************************************************************
250 : :
251 : : DBG_ENTER("PERS")
252 : :
253 : 5 : nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
254 : :
255 [ - + ]: 5 : if (BSp%have_complex_ene) then ! QP lifetimes are not included
256 : 0 : ABI_ERROR("complex energies not coded yet")
257 : : end if
258 : :
259 [ + + + - : 11 : if (ANY(Bsp%nreh/=Bsp%nreh(1))) then
- + ]
260 : 0 : write(std_out,*)" Bsp%nreh: ",Bsp%nreh
261 : 0 : write(msg,'(a)')" BSE code does not support different number of transitions for the two spin channels"
262 : 0 : ABI_WARNING(msg)
263 : : end if
264 : :
265 : 5 : nsppol = Hdr_bse%nsppol
266 [ + + ]: 11 : exc_size = SUM(BSp%nreh)
267 : 5 : nstates = BSp%nstates; do_full_diago=(Bsp%nstates==exc_size)
268 : :
269 : 5 : neh1 = Bsp%nreh(1); neh2 = neh1
270 : : if (Hdr_bse%nsppol==2) neh2 = Bsp%nreh(2)
271 : :
272 : : ! Scalapack is disabled due to portability issues in slk_read
273 : : ! This part should be rewritten with hdf5 + mpi-io
274 : :
275 : : use_scalapack = .FALSE.
276 : : !#if defined HAVE_LINALG_SCALAPACK
277 : : ! use_scalapack = (nprocs > 1)
278 : : !#endif
279 : : if (use_scalapack .and. nsppol == 2) then
280 : : use_scalapack = .False.
281 : : ABI_WARNING("Scalapack with nsppol==2 not yet available. Using sequential version")
282 : : end if
283 : :
284 [ - + ]: 5 : if (.not.use_scalapack .and. my_rank/=master) GOTO 10 ! Inversion is done by master only.
285 : :
286 [ - + ]: 5 : nene_printed = MIN(32*nsppol,nstates); if (prtvol>10) nene_printed = nstates
287 : :
288 [ + + ]: 5 : if (BS_files%in_hreso /= BSE_NOFILE) then
289 : 1 : hreso_fname = BS_files%in_hreso
290 : : else
291 : 4 : hreso_fname = BS_files%out_hreso
292 : : end if
293 : :
294 : 5 : bseig_fname = BS_files%out_eig
295 [ - + ]: 5 : if (BS_files%in_eig /= BSE_NOFILE) then
296 : 0 : ABI_ERROR("BS_files%in_eig is defined!")
297 : : end if
298 : :
299 : 5 : write(msg,'(a,i0)')' Direct diagonalization of the resonant excitonic Hamiltonian, Matrix size= ',exc_size
300 [ + + ]: 15 : call wrtout([std_out, ab_out], msg)
301 : :
302 [ - + + - : 15 : ABI_MALLOC_OR_DIE(exc_ene,(exc_size), ierr)
+ - + - +
- - + -
+ ]
303 [ - + + - : 15 : ABI_MALLOC_OR_DIE(exc_ene_c,(exc_size), ierr)
+ - + - +
- - + -
+ ]
304 : :
305 : 5 : do_ep_renorm = .FALSE.
306 : 5 : ntemp = 1
307 : 5 : do_ep_lifetime = .FALSE.
308 [ + + + + ]: 5 : if(BSp%do_ep_renorm .and. present(Epren)) then
309 : 1 : do_ep_renorm = .TRUE.
310 : 1 : ntemp = Epren%ntemp
311 [ + - ]: 1 : if(present(elph_lifetime)) then
312 : 1 : do_ep_lifetime = elph_lifetime
313 : : end if
314 : : end if
315 : :
316 : : if (do_ep_renorm) then
317 [ - + ]: 1 : ABI_CHECK(nsppol == 1, "Nsppol == 2 not supported with elphon renormalizations")
318 : : end if
319 : :
320 : : SELECT CASE (use_scalapack)
321 : : CASE (.FALSE.)
322 : :
323 : 5 : write(msg,'(a)')". Using LAPACK sequential version. "
324 [ + + ]: 15 : call wrtout([std_out, ab_out], msg)
325 : 5 : write(msg,'(a,f8.1,a)')' Allocating excitonic eigenvalues. Memory required: ', exc_size*dp*b2Mb,' Mb. '
326 : 5 : call wrtout(std_out, msg)
327 : :
328 : 5 : write(msg,'(a,f8.1,a)')' Allocating excitonic hamiltonian. Memory required: ',exc_size**2*dp*b2Mb,' Mb.'
329 : 5 : call wrtout(std_out, msg, do_flush=.True.)
330 : :
331 [ - + + - : 20 : ABI_MALLOC_OR_DIE(exc_mat,(exc_size,exc_size), ierr)
- + + - +
- + - + -
- + - + ]
332 : :
333 [ + + ]: 5 : if (do_ep_renorm) then
334 [ - + + - : 3 : ABI_MALLOC_OR_DIE(exc_vl,(exc_size,exc_size),ierr)
- + + - +
- - + -
+ ]
335 : : end if
336 : : !exc_mat = HUGE(zero)
337 : : !
338 : : ! Read data from file.
339 [ - + ]: 5 : if (open_file(hreso_fname,msg,newunit=hreso_unt,form="unformatted",status="old",action="read") /= 0) then
340 : 0 : ABI_ERROR(msg)
341 : : end if
342 : : !
343 : : ! Read the header and perform consistency checks.
344 : 5 : call exc_read_bshdr(hreso_unt,Bsp,fform,ierr)
345 [ - + ]: 5 : ABI_CHECK(ierr==0,"Fatal error, cannot continue")
346 : : !
347 : : ! Construct full resonant block using Hermiticity.
348 : 5 : diagonal_is_real = .not.Bsp%have_complex_ene
349 : 5 : call exc_read_rblock_fio(hreso_unt,diagonal_is_real,nsppol,Bsp%nreh,exc_size,exc_mat,ierr)
350 [ - + ]: 5 : ABI_CHECK(ierr==0,"Fatal error, cannot continue")
351 : :
352 : 5 : close(hreso_unt)
353 : :
354 [ + + ]: 5 : if (do_ep_renorm) then
355 : 1 : write(std_out,'(a)') "Mapping kpts from bse to eph"
356 : 1 : sppoldbl = 1 !; if (any(Cryst%symafm == -1) .and. Epren%nsppol == 1) nsppoldbl=2
357 [ - + + - : 3 : ABI_MALLOC(bs2eph, (Kmesh%nbz*sppoldbl, 6))
+ - + - -
+ - + ]
358 : 1 : timrev = 1
359 : : call listkk(dksqmax, Cryst%gmet, bs2eph, Epren%kpts, Kmesh%bz, Epren%nkpt, Kmesh%nbz, Cryst%nsym, &
360 [ + - ]: 1 : sppoldbl, Cryst%symafm, Cryst%symrel, timrev, xmpi_comm_self, use_symrec=.False.)
361 : : end if
362 : :
363 [ + + ]: 19 : do itemp = 1, ntemp
364 : :
365 : : !TODO should find a way not to read again and again !
366 : : ! but without storing it twice !!!
367 : : !exc_mat = HUGE(zero)
368 : : !
369 : : ! Read data from file.
370 [ - + ]: 14 : if (open_file(hreso_fname,msg,newunit=hreso_unt,form="unformatted",status="old",action="read") /= 0) then
371 : 0 : ABI_ERROR(msg)
372 : : end if
373 : : !
374 : : ! Read the header and perform consistency checks.
375 : 14 : call exc_read_bshdr(hreso_unt,Bsp,fform,ierr)
376 [ - + ]: 14 : ABI_CHECK(ierr==0,"Fatal error, cannot continue")
377 : : !
378 : : ! Construct full resonant block using Hermiticity.
379 : 14 : diagonal_is_real = .not.Bsp%have_complex_ene
380 : 14 : call exc_read_rblock_fio(hreso_unt,diagonal_is_real,nsppol,Bsp%nreh,exc_size,exc_mat,ierr)
381 [ - + ]: 14 : ABI_CHECK(ierr==0,"Fatal error, cannot continue")
382 : :
383 : 14 : close(hreso_unt)
384 : :
385 : 14 : bseig_fname = BS_files%out_eig
386 : :
387 [ + + ]: 14 : if (do_ep_renorm) then
388 : 10 : write(std_out,'(a,i4)') "Will perform elphon renormalization for itemp = ",itemp
389 : :
390 : 10 : call int2char4(itemp,ts)
391 : :
392 [ + - ]: 10 : bseig_fname = TRIM(BS_files%out_eig) // TRIM("_T") // ts
393 : :
394 : : ! Should patch the diagonal of exc_mat
395 : :
396 [ + + ]: 20 : do isppol = 1, BSp%nsppol
397 [ + + ]: 980 : do ireh = 1, BSp%nreh(isppol)
398 : 960 : ic = BSp%Trans(ireh,isppol)%c
399 : 960 : iv = BSp%Trans(ireh,isppol)%v
400 : 960 : ik = BSp%Trans(ireh,isppol)%k ! In the full bz
401 : 960 : en = BSp%Trans(ireh,isppol)%en
402 : :
403 : 960 : ep_ik = bs2eph(ik,1)
404 : :
405 : : !TODO support multiple spins !
406 [ - + ]: 960 : if(ABS(en - (Epren%eigens(ic,ep_ik,isppol)-Epren%eigens(iv,ep_ik,isppol)+BSp%mbpt_sciss)) > tol3) then
407 : 0 : ABI_ERROR("Eigen from the transition does not correspond to the EP file !")
408 : : end if
409 : 960 : exc_mat(ireh,ireh) = exc_mat(ireh,ireh) + (Epren%renorms(1,ic,ik,isppol,itemp) - Epren%renorms(1,iv,ik,isppol,itemp))
410 : :
411 : : ! Add lifetime
412 [ + - ]: 970 : if(do_ep_lifetime) then
413 : : exc_mat(ireh,ireh) = exc_mat(ireh,ireh) - j_dpc*(Epren%linewidth(1,ic,ik,isppol,itemp) &
414 : 960 : & + Epren%linewidth(1,iv,ik,isppol,itemp))
415 : : end if
416 : :
417 : : end do
418 : : end do
419 : :
420 : : end if
421 : :
422 [ + + ]: 14 : if (do_full_diago) then
423 [ + + ]: 13 : if(do_ep_renorm) then
424 : 10 : call wrtout(std_out," Full diagonalization with XGEEV... ")
425 [ - + + - : 30 : ABI_MALLOC(exc_vec,(exc_size,exc_size))
- + + - -
+ - + ]
426 : 10 : call xgeev('V','V',exc_size,exc_mat,exc_size,exc_ene_c,exc_vl,exc_size,exc_vec,exc_size)
427 [ + + + + ]: 93130 : exc_mat(:,1:nstates) = exc_vec
428 [ - + ]: 10 : ABI_FREE(exc_vec)
429 : : else
430 : 3 : call wrtout(std_out," Full diagonalization with XHEEV... ")
431 : 3 : call xheev("Vectors","Upper",exc_size,exc_mat,exc_ene)
432 [ + + ]: 1251 : exc_ene_c(:) = exc_ene(:)
433 : : end if
434 : : else
435 : 1 : call wrtout(std_out," Partial diagonalization with XHEEVX... ")
436 : 1 : abstol=zero; il=1; iu=nstates
437 [ - + + - : 4 : ABI_MALLOC_OR_DIE(exc_vec,(exc_size,nstates),ierr)
- + + - +
- + - + -
- + - + ]
438 : 1 : call xheevx("Vectors","Index","Upper",exc_size,exc_mat,vl,vu,il,iu,abstol,mene_found,exc_ene,exc_vec,exc_size)
439 [ + + + + ]: 3851 : exc_mat(:,1:nstates) = exc_vec
440 [ + + ]: 385 : exc_ene_c(:) = exc_ene(:)
441 [ - + ]: 1 : ABI_FREE(exc_vec)
442 : : end if
443 : : !
444 : : ! ==============================================
445 : : ! === Now exc_mat contains the eigenvectors ====
446 : : ! ==============================================
447 : :
448 : : ! * Write the final results.
449 [ + - ]: 14 : call wrtout(std_out,' Writing eigenvalues and eigenvectors to file: '//TRIM(bseig_fname))
450 : :
451 [ - + ]: 14 : if (open_file(bseig_fname,msg,newunit=eig_unt,form="unformatted",action="write") /= 0) then
452 : 0 : ABI_ERROR(msg)
453 : : end if
454 : :
455 : : !!! !DBYG
456 : : !!! !Compute overlap matrix
457 : : !!! ABI_MALLOC(ovlp,(exc_size,exc_size))
458 : : !!! do mi=1,nstates
459 : : !!! do ireh=1,nstates
460 : : !!! ovlp(mi,ireh) = xdotc(exc_size,exc_vl(:,mi),1,exc_mat(:,ireh),1)
461 : : !!! if(mi==ireh) then
462 : : !!! !if(ABS(ovlp(mi,ireh)) < 0.999) then
463 : : !!! ! write(*,*) "it,itp = ",mi,ireh,"ovlp = ",ovlp(mi,ireh)
464 : : !!! !end if
465 : : !!! else
466 : : !!! if(ABS(ovlp(mi,ireh)) > 0.001) then
467 : : !!! write(*,*) "it,itp = ",mi,ireh,"ovlp = ",ovlp(mi,ireh)
468 : : !!! end if
469 : : !!! end if
470 : : !!! end do
471 : : !!! end do
472 : : !!! !call xgemm("C","N",exc_size,nstates,nstates,cone,exc_vl,exc_size,exc_mat,exc_size,czero,ovlp,nstates)
473 : :
474 : : !!! write(777,*) ovlp
475 : : !!! ABI_FREE(ovlp)
476 : : !!! !ENDDBYG
477 : :
478 : : !% fform = 1002 ! FIXME
479 : : !% call hdr_io_int(fform,Hdr_bse,2,eig_unt)
480 : :
481 : 14 : write(eig_unt) do_ep_lifetime
482 : 14 : write(eig_unt) exc_size, nstates
483 : 14 : write(eig_unt) exc_ene_c(1:nstates)
484 [ + + ]: 2232 : do mi=1,nstates
485 : 2218 : write(eig_unt) exc_mat(1:exc_size,mi)
486 [ + + ]: 2232 : if(do_ep_lifetime) then
487 : 960 : write(eig_unt) exc_vl(1:exc_size,mi)
488 : : end if
489 : : end do
490 : :
491 : 47 : close(eig_unt)
492 : :
493 : : end do ! itemp
494 : :
495 [ - + ]: 5 : ABI_FREE(exc_mat)
496 [ + + ]: 15 : if (do_ep_renorm) then
497 [ - + ]: 1 : ABI_FREE(exc_vl)
498 [ - + ]: 1 : ABI_FREE(bs2eph)
499 : : end if
500 : :
501 : : CASE (.TRUE.)
502 : :
503 : : #if defined HAVE_LINALG_SCALAPACK && defined HAVE_MPI_IO
504 : : if (nsppol==2) then
505 : : ABI_WARNING("nsppol==2 + scalapack not coded yet")
506 : : end if
507 : :
508 : : istwf_k=1; tbloc=50
509 : : write(msg,'(2(a,i0))')". Using scaLAPACK version with nprocs= ",nprocs,"; block size= ",tbloc
510 : : call wrtout([std_out, ab_out], msg)
511 : :
512 : : write(msg,'(a,f8.1,a)')' Allocating excitonic eigenvalues. Memory required: ',exc_size*dp*b2Mb,' Mb. '
513 : : call wrtout(std_out, msg)
514 : : !
515 : : ! Init scaLAPACK environment.
516 : : call Slk_processor%init(comm)
517 : : !
518 : : ! Init scaLAPACK matrices
519 : : call Slk_mat%init(exc_size,exc_size,Slk_processor,istwf_k)
520 : : call Slk_vec%init(exc_size,exc_size,Slk_processor,istwf_k)
521 : : !
522 : : ! Open the file with MPI-IO and skip the record.
523 : : amode=MPI_MODE_RDONLY
524 : :
525 : : call MPI_FILE_OPEN(comm, hreso_fname, amode, MPI_INFO_NULL, mpi_fh, ierr)
526 : : ABI_CHECK_MPI(ierr,"MPI_IO error opening file: "//TRIM(hreso_fname))
527 : :
528 : : ! Skip the header and find the offset for reading the matrix.
529 : : call exc_skip_bshdr_mpio(mpi_fh,xmpio_collective,ehdr_offset)
530 : : !
531 : : ! Read scaLAPACK matrix from the file.
532 : : if (nsppol==1) then
533 : : call slk_read(Slk_mat,"Upper","Hermitian",is_fortran_file,mpi_fh=mpi_fh,offset=ehdr_offset)
534 : : else
535 : : array_of_sizes = (/exc_size,exc_size/)
536 : : block_sizes(:,1) = (/neh1,neh1/)
537 : : block_sizes(:,2) = (/neh2,neh2/)
538 : : block_sizes(:,3) = (/neh1,neh2/)
539 : : ABI_ERROR("Not tested")
540 : : !call slk_read_from_blocks(Slk_mat,array_of_sizes,block_sizes,is_fortran_file,mpi_fh=mpi_fh,offset=ehdr_offset)
541 : : end if
542 : :
543 : : call MPI_FILE_CLOSE(mpi_fh, ierr)
544 : : ABI_CHECK_MPI(ierr,"FILE_CLOSE")
545 : :
546 : : if (do_full_diago) then
547 : : call wrtout(std_out," Performing full diagonalization with scaLAPACK...")
548 : : call slk_mat%heev("Vectors","Upper",Slk_vec,exc_ene)
549 : : else
550 : : call wrtout(std_out," Performing partial diagonalization with scaLAPACK...")
551 : : il=1; iu=nstates; abstol=zero !ABSTOL = PDLAMCH(comm,'U')
552 : : call slk_mat%pzheevx("Vectors","Index","Upper",vl,vu,il,iu,abstol,Slk_vec,mene_found,exc_ene)
553 : : end if
554 : :
555 : : exc_ene_c(:) = exc_ene(:)
556 : :
557 : : call Slk_mat%free()
558 : :
559 : : call wrtout(std_out,' Writing eigenvalues/vectors to file: '//TRIM(bseig_fname), do_flush=.True.)
560 : :
561 : : ! Write distributed matrix on file bseig_fname with Fortran records.
562 : : if (my_rank==master) then ! Write exc eigenvalues. Vectors will be appended in slk_write.
563 : : if (open_file(bseig_fname,msg,newunit=eig_unt,form="unformatted",action="write") /= 0) then
564 : : ABI_ERROR(msg)
565 : : end if
566 : : write(eig_unt) exc_size, nstates
567 : : write(eig_unt) exc_ene_c(1:nstates)
568 : : close(eig_unt)
569 : : end if
570 : :
571 : : call xmpi_barrier(comm)
572 : : !
573 : : ! Open the file with MPI-IO and skip the record.
574 : : amode=MPI_MODE_RDWR
575 : :
576 : : call MPI_FILE_OPEN(comm, bseig_fname, amode, MPI_INFO_NULL, mpi_fh, ierr)
577 : : ABI_CHECK_MPI(ierr,"MPI_IO error opening file: "//TRIM(hreso_fname))
578 : :
579 : : !call MPI_FILE_SYNC(mpi_fh,ierr)
580 : :
581 : : ehdr_offset = 0
582 : : call xmpio_read_frm(mpi_fh,ehdr_offset,xmpio_collective,fmarker,ierr)
583 : : call xmpio_read_frm(mpi_fh,ehdr_offset,xmpio_collective,fmarker,ierr)
584 : :
585 : : write(std_out,*)"Writing nstates ",nstates
586 : : gsub(:,1) = (/1,1/)
587 : : gsub(:,2) = (/exc_size,nstates/)
588 : : call slk_write(Slk_vec,"All",is_fortran_file,mpi_fh=mpi_fh,offset=ehdr_offset,glob_subarray=gsub)
589 : :
590 : : call MPI_FILE_CLOSE(mpi_fh, ierr)
591 : : ABI_CHECK_MPI(ierr,"FILE_CLOSE")
592 : :
593 : : call Slk_vec%free()
594 : : call Slk_processor%free()
595 : : call xmpi_barrier(comm)
596 : : #else
597 : : ABI_BUG("You should not be here!")
598 : : #endif
599 : :
600 : : END SELECT
601 : :
602 : : ! Order the eigenvalues
603 [ + + ]: 1354 : do ii=nstates,2,-1
604 [ + + ]: 378583 : do j=1,ii-1
605 [ + + ]: 378578 : if (DBLE(exc_ene_c(j)) > DBLE(exc_ene_c(j+1))) then
606 : 3583 : ctemp = exc_ene_c(j)
607 : 3583 : exc_ene_c(j) = exc_ene_c(j+1)
608 : 3583 : exc_ene_c(j+1) = ctemp
609 : : end if
610 : : end do
611 : : end do
612 : :
613 : :
614 : 5 : write(msg,'(a,i4)')' Excitonic eigenvalues in eV up to n= ',nene_printed
615 [ + + ]: 15 : call wrtout([std_out, ab_out], msg)
616 : :
617 [ + + ]: 27 : do it=0,(nene_printed-1)/8
618 [ + - + + ]: 192 : write(msg,'(8f10.5)') ( DBLE(exc_ene_c(ii))*Ha_eV, ii=1+it*8,MIN(it*8+8,nene_printed) )
619 [ + + ]: 71 : call wrtout([std_out, ab_out], msg)
620 : : end do
621 : :
622 [ + - - + : 1364 : exc_gap = MINVAL(DBLE(exc_ene_c(1:nstates)))
- - + + ]
623 [ + - - + : 1364 : exc_maxene = MAXVAL(DBLE(exc_ene_c(1:nstates)))
- - + + ]
624 : :
625 : 5 : write(msg,'(a,2(a,f7.2,2a),a)')ch10,&
626 : 5 : " First excitonic eigenvalue= ",exc_gap*Ha_eV, " [eV]",ch10,&
627 : 10 : " Last excitonic eigenvalue= ",exc_maxene*Ha_eV," [eV]",ch10,ch10
628 [ + + ]: 15 : call wrtout([std_out, ab_out], msg, do_flush=.True.)
629 : :
630 [ - + ]: 5 : ABI_FREE(exc_ene_c)
631 [ - + ]: 5 : ABI_FREE(exc_ene)
632 : :
633 : 5 : 10 call xmpi_barrier(comm)
634 : :
635 : : DBG_EXIT("PERS")
636 : :
637 [ - + - + ]: 5 : end subroutine exc_diago_resonant
638 : : !!***
639 : :
640 : : !----------------------------------------------------------------------
641 : :
642 : : !!****f* m_exc_diago/exc_print_eig
643 : : !! NAME
644 : : !! exc_print_eig
645 : : !!
646 : : !! FUNCTION
647 : : !! Print excitonic eigenvalues on std_out and ab_out.
648 : : !!
649 : : !! INPUTS
650 : : !! gw_gap=GW direct gap.
651 : : !! bseig_fname=The name of file containing eigenvalues and eigenvectors
652 : : !!
653 : : !! OUTPUT
654 : : !! exc_gap=Excitonic direct gap.
655 : : !! Additional info on the Excitonic spectrum are reported on standard output.
656 : : !!
657 : : !! SOURCE
658 : :
659 : 6 : subroutine exc_print_eig(BSp,bseig_fname,gw_gap,exc_gap)
660 : :
661 : : !Arguments ------------------------------------
662 : : !scalars
663 : : complex(dp),intent(in) :: gw_gap
664 : : complex(dp),intent(out) :: exc_gap
665 : : character(len=*),intent(in) :: bseig_fname
666 : : type(excparam),intent(in) :: BSp
667 : :
668 : : !Local variables ------------------------------
669 : : !scalars
670 : : integer :: nstates_read,ii,j,k,eig_unt,ieig,hsize_exp, hsize_read !,nstates
671 : : complex(dp) :: bind_energy,ctemp
672 : : character(len=500) :: msg
673 : : !type(Hdr_type) :: tmp_Hdr
674 : : !arrays
675 : 6 : integer,allocatable :: iperm(:)
676 : 6 : real(dp),allocatable :: exc_rene(:)
677 : 6 : complex(dp),allocatable :: exc_cene(:)
678 : : !************************************************************************
679 : :
680 : 6 : exc_gap = czero
681 : :
682 [ - + ]: 6 : if (open_file(bseig_fname,msg,newunit=eig_unt,form="unformatted",status="old",action="read") /= 0) then
683 : 0 : ABI_ERROR(msg)
684 : : end if
685 : :
686 : 6 : read(eig_unt) ! do_ep_lifetime
687 : 6 : read(eig_unt) hsize_read, nstates_read
688 : :
689 [ + + + + ]: 12 : if (BSp%use_coupling==0) hsize_exp = SUM(Bsp%nreh)
690 [ + + + + ]: 7 : if (BSp%use_coupling>0) hsize_exp = 2*SUM(Bsp%nreh)
691 : :
692 [ - + ]: 6 : if (hsize_exp /= hsize_read) then
693 : 0 : write(msg,'(2(a,i0))')" Wrong dimension: read: ",hsize_read," expected= ",hsize_exp
694 : 0 : ABI_ERROR(msg)
695 : : end if
696 : :
697 [ - + + - : 18 : ABI_MALLOC(exc_cene,(nstates_read))
+ - + - -
+ - + ]
698 : 6 : read(eig_unt) exc_cene(:)
699 : :
700 [ - + + - : 18 : ABI_MALLOC(exc_rene,(nstates_read))
+ - + - -
+ - + ]
701 [ + - - + : 2048 : exc_rene = DBLE(exc_cene)
- - + + ]
702 : :
703 [ - + + - : 18 : ABI_MALLOC(iperm,(nstates_read))
+ - + - -
+ - + ]
704 [ + - + + : 6126 : iperm = (/(ii, ii=1,nstates_read)/)
+ - + + +
- - + - -
+ + ]
705 : :
706 : 6 : call sort_dp(nstates_read,exc_rene,iperm,tol6)
707 : :
708 [ - + ]: 6 : ABI_FREE(exc_rene)
709 [ - + ]: 6 : ABI_FREE(iperm)
710 : :
711 : : ! put in ascending order
712 [ + + ]: 2036 : do ii=nstates_read,2,-1
713 [ + + ]: 669278 : do j=1,ii-1
714 [ + - ]: 669272 : if (DBLE(exc_cene(j)) > DBLE(exc_cene(j+1))) then
715 : 0 : ctemp = exc_cene(j)
716 : 0 : exc_cene(j) = exc_cene(j+1)
717 : 0 : exc_cene(j+1) = ctemp
718 : : end if
719 : : end do
720 : : end do
721 : :
722 : 6 : exc_gap = DCMPLX(ABS(DBLE(exc_cene(1))),AIMAG(exc_cene(1)))
723 : :
724 [ + + ]: 2042 : do ii=1,nstates_read
725 [ + + ]: 2042 : if (ABS(DBLE(exc_cene(ii))) < DBLE(exc_gap)) then
726 : 383 : exc_gap = DCMPLX(ABS(DBLE(exc_cene(ii))),AIMAG(exc_cene(ii)))
727 : : end if
728 : : end do
729 : :
730 : 6 : bind_energy = gw_gap - exc_gap
731 : :
732 : : write(msg,"(3(a,2f6.2,2a))")&
733 : 6 : " GW direct gap ",gw_gap*Ha_eV, " [eV] ",ch10,&
734 : 6 : " EXC direct gap ",exc_gap*Ha_eV, " [eV] ",ch10,&
735 : 12 : " EXC binding energy ",bind_energy*Ha_eV," [eV] ",ch10
736 [ + + ]: 18 : call wrtout([std_out, ab_out], msg)
737 : :
738 : 6 : msg=' Excitonic eigenvalues up to the GW energy gap [eV]'
739 [ + + ]: 18 : call wrtout([std_out, ab_out], msg)
740 : :
741 [ - + ]: 390 : do ii=1,nstates_read
742 [ + + ]: 390 : if (DBLE(exc_cene(ii)) > zero) EXIT
743 : : end do
744 : :
745 [ + + ]: 96 : do j=ii,nstates_read
746 [ + + ]: 96 : if (DBLE(exc_cene(j)) > DBLE(gw_gap)) EXIT
747 : : end do
748 : 6 : j=j-1
749 : :
750 [ + + ]: 96 : do ieig=ii,j
751 : 90 : write(msg,'(i3,a,2f6.2,a)')ieig," (",exc_cene(ieig)*Ha_eV,")"
752 [ + + ]: 276 : call wrtout([std_out, ab_out], msg)
753 : : end do
754 : :
755 : 6 : ii=ii-1
756 [ + + ]: 22 : do j=ii,1,-1
757 [ + + ]: 22 : if (ABS(DBLE(exc_cene(j))) > DBLE(gw_gap)) EXIT
758 : : end do
759 : 6 : j=j+1
760 : :
761 : : ! This coding is not portable, write to ab_out has been disabled.
762 [ + + ]: 6 : if (ii>0) then
763 [ + + ]: 17 : do k=ii,j,-1
764 : 16 : write(msg,'(i3,a,2f6.2,a)')k," (",exc_cene(k)*Ha_eV,")"
765 : 17 : call wrtout(std_out, msg)
766 : : end do
767 : : end if
768 : :
769 : 6 : ABI_FREE(exc_cene)
770 : :
771 : 6 : close(eig_unt)
772 : :
773 : 6 : end subroutine exc_print_eig
774 : : !!***
775 : :
776 : : !----------------------------------------------------------------------
777 : :
778 : : !!****f* m_exc_diago/exc_diago_coupling
779 : : !! NAME
780 : : !! exc_diago_coupling
781 : : !!
782 : : !! FUNCTION
783 : : !! Calculate excitonic eigenvalues and eigenvectors by performing a direct diagonalization.
784 : : !! of the non Hermitian excitonic Hamiltoninan (resonant + coupling).
785 : : !!
786 : : !! INPUTS
787 : : !! bseig_fname=The name of the output file.
788 : : !! Bsp
789 : : !! neh=Rank of the resonant block of the Hamiltoninan (equal to the rank of the coupling part)
790 : : !! comm=MPI communicator.
791 : : !! BS_files<excfiles>=Datatype storing names and files used in the Bethe-Salpeter code.
792 : : !!
793 : : !! OUTPUT
794 : : !! Excitonic eigenvectors and eigenvalues are written on file BS_files%out_eig.
795 : : !!
796 : : !! SOURCE
797 : :
798 : 0 : subroutine exc_diago_coupling(Bsp,BS_files,Hdr_bse,prtvol,comm)
799 : :
800 : : !Arguments ------------------------------------
801 : : !scalars
802 : : integer,intent(in) :: comm,prtvol
803 : : type(excfiles),intent(in) :: BS_files
804 : : type(excparam),intent(in) :: BSp
805 : : type(Hdr_type),intent(in) :: Hdr_bse
806 : :
807 : : !Local variables ------------------------------
808 : : !scalars
809 : : integer,parameter :: master=0,ldvl=1
810 : : integer(i8b) :: bsize_ham
811 : : integer :: ii,exc_size,hreso_unt,hcoup_unt,eig_unt,nsppol,nstates
812 : : integer :: bsz,block,bs1,bs2,jj, fform,row_sign, mi,it,nprocs,my_rank !itp
813 : : integer :: nene_printed,ierr
814 : : real(dp) :: exc_gap,exc_maxene,temp
815 : : logical :: diago_is_real,do_full_diago, do_ep_lifetime
816 : : character(len=500) :: msg
817 : : character(len=fnlen) :: hreso_fname,hcoup_fname,bseig_fname
818 : : !arrays
819 : 0 : complex(dp),allocatable :: exc_ham(:,:),exc_rvect(:,:),exc_ene(:),ovlp(:,:), cbuff(:,:)
820 : : complex(dp) :: vl_dpc(ldvl,1)
821 : : !************************************************************************
822 : :
823 : 0 : my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm)
824 : :
825 : 0 : nsppol = Hdr_bse%nsppol
826 [ # # ]: 0 : if (nsppol==2) then
827 : 0 : ABI_WARNING("nsppol==2 with coupling is still under development")
828 : : end if
829 : :
830 [ # # ]: 0 : if (nprocs > 1) then
831 : 0 : ABI_WARNING("Scalapack does not provide ZGEEV, diagonalization is done in sequential!")
832 : : end if
833 : :
834 [ # # ]: 0 : exc_size = 2*SUM(BSp%nreh)
835 : 0 : nstates = BSp%nstates
836 : 0 : do_full_diago = (exc_size==nstates)
837 [ # # ]: 0 : ABI_CHECK(do_full_diago,"Partial diago not coded yet")
838 : :
839 : 0 : bseig_fname = BS_files%out_eig
840 [ # # ]: 0 : if (BS_files%in_eig /= BSE_NOFILE) then
841 : 0 : ABI_ERROR("BS_files%in_eig is defined!")
842 : : end if
843 : : !
844 : : ! Only master performs the diagonalization since ScaLAPACK does not provide the parallel version of ZGEEV.
845 [ # # ]: 0 : if (my_rank/=master) GOTO 10
846 : :
847 : 0 : write(msg,'(a,i0)')' Direct diagonalization of the full excitonic Hamiltonian, Matrix size= ',exc_size
848 [ # # ]: 0 : call wrtout([std_out, ab_out], msg)
849 : :
850 : 0 : bsize_ham = 2*dp*exc_size**2
851 : 0 : write(msg,'(a,f9.2,a)')' Allocating full excitonic Hamiltonian. Memory requested: ',bsize_ham*b2Gb,' Gb. '
852 : 0 : call wrtout(std_out, msg)
853 : :
854 [ # # # # : 0 : ABI_MALLOC_OR_DIE(exc_ham,(exc_size,exc_size), ierr)
# # # # #
# # # # #
# # # # ]
855 : :
856 : : write(msg,'(3a,f8.1,3a,f8.1,a)')&
857 : 0 : ' Allocating excitonic eigenvalues and eigenvectors. ',ch10,&
858 : 0 : ' Memory-space requested: ',2*dp*exc_size*b2Gb,' Gb. ',ch10,&
859 : 0 : ' Memory-space requested: ',bsize_ham*b2Gb,' Gb. '
860 : 0 : call wrtout(std_out, msg)
861 : :
862 [ # # # # : 0 : ABI_MALLOC_OR_DIE(exc_ene,(exc_size), ierr)
# # # # #
# # # #
# ]
863 : :
864 [ # # ]: 0 : if (BS_files%in_hreso /= BSE_NOFILE) then
865 : 0 : hreso_fname = BS_files%in_hreso
866 : : else
867 : 0 : hreso_fname = BS_files%out_hreso
868 : : end if
869 : :
870 [ # # ]: 0 : call wrtout(std_out,' Reading resonant excitonic Hamiltonian from '//TRIM(hreso_fname))
871 : :
872 [ # # ]: 0 : if (open_file(hreso_fname,msg,newunit=hreso_unt,form="unformatted",status="old",action="read") /= 0) then
873 : 0 : ABI_ERROR(msg)
874 : : end if
875 : : !
876 : : ! Read the header and perform consistency checks.
877 : 0 : call exc_read_bshdr(hreso_unt,Bsp,fform,ierr)
878 [ # # ]: 0 : ABI_CHECK(ierr==0,"Wrong header")
879 : : !
880 : : ! Construct resonant and anti-resonant part of the excitonic Hamiltonian using Hermiticity. File is always in double precision.
881 : : ! Fill exc_ham with ( R 0 )
882 : : ! ( 0 -R*)
883 : : !BEGINDEBUG
884 [ # # # # ]: 0 : exc_ham = HUGE(one)
885 : : !ENDDEBUG
886 : :
887 : 0 : row_sign=-1; diago_is_real=(.not.BSp%have_complex_ene)
888 : 0 : call exc_fullh_from_blocks(hreso_unt,"Resonant",nsppol,row_sign,diago_is_real,BSp%nreh,exc_size,exc_ham)
889 : 0 : close(hreso_unt)
890 : :
891 [ # # ]: 0 : if (BS_files%in_hcoup /= BSE_NOFILE) then
892 : 0 : hcoup_fname = BS_files%in_hcoup
893 : : else
894 : 0 : hcoup_fname = BS_files%out_hcoup
895 : : end if
896 : :
897 [ # # ]: 0 : call wrtout(std_out,' Reading coupling excitonic Hamiltonian from '//TRIM(hcoup_fname))
898 [ # # ]: 0 : if (open_file(hcoup_fname,msg,newunit=hcoup_unt,form="unformatted",status="old",action="read") /= 0) then
899 : 0 : ABI_ERROR(msg)
900 : : end if
901 : : !
902 : : ! Read the header and perform consistency checks.
903 : 0 : call exc_read_bshdr(hcoup_unt,Bsp,fform,ierr)
904 [ # # ]: 0 : ABI_CHECK(ierr==0,"Wrong header")
905 : : !
906 : : ! Fill exc_ham with ( 0 C) to have ( R C )
907 : : ! (-C* 0) (-C* -R*)
908 : 0 : row_sign=-1; diago_is_real=(.not.BSp%have_complex_ene) ! not used here
909 : 0 : call exc_fullh_from_blocks(hcoup_unt,"Coupling",nsppol,row_sign,diago_is_real,BSp%nreh,exc_size,exc_ham)
910 : :
911 : : !BEGINDEBUG
912 [ # # # # : 0 : if (ANY(exc_ham==HUGE(one))) then
# # # # ]
913 [ # # # # : 0 : write(msg,'(a,2(1x,i0))')"There is a bug in exc_fullh_from_blocks",COUNT(exc_ham==HUGE(one)),exc_size**2
# # ]
914 : 0 : ABI_WARNING(msg)
915 : 0 : bsz = Bsp%nreh(1)
916 [ # # # # : 0 : ABI_MALLOC(cbuff,(bsz,bsz))
# # # # #
# # # # #
# # ]
917 : 0 : block=0
918 [ # # ]: 0 : do jj=1,2*nsppol
919 [ # # ]: 0 : do ii=1,2*nsppol
920 : 0 : block=block+1
921 : 0 : bs1 = (ii-1)*bsz+1
922 : 0 : bs2 = (jj-1)*bsz+1
923 [ # # # # : 0 : cbuff = exc_ham(bs1:bs1+bsz-1,bs2:bs2+bsz-1)
# # # # #
# # # ]
924 [ # # # # : 0 : if (ANY(cbuff==HUGE(one))) then
# # # # ]
925 [ # # # # : 0 : write(std_out,*)" for block ",ii,jj," found ",COUNT(cbuff==HUGE(one))," wrong entries"
# # ]
926 : : end if
927 : : end do
928 : : end do
929 : :
930 [ # # ]: 0 : ABI_FREE(cbuff)
931 : 0 : ABI_ERROR("Cannot continue")
932 : : end if
933 : : !ENDDEBUG
934 : :
935 : 0 : close(hcoup_unt)
936 : : !
937 : : ! ======================================================
938 : : ! ==== Calculate right eigenvectors and eigenvalues ====
939 : : ! ======================================================
940 [ # # # # : 0 : ABI_MALLOC_OR_DIE(exc_rvect,(exc_size,exc_size), ierr)
# # # # #
# # # #
# ]
941 : :
942 [ # # ]: 0 : if (do_full_diago) then
943 : 0 : call wrtout(std_out,"Complete direct diagonalization with xgeev...")
944 : 0 : call xgeev("No_left_eigen","Vectors",exc_size,exc_ham,exc_size,exc_ene,vl_dpc,ldvl,exc_rvect,exc_size)
945 : : else
946 : 0 : ABI_ERROR("Not implemented error")
947 : : end if
948 : :
949 [ # # ]: 0 : ABI_FREE(exc_ham)
950 : :
951 [ # # # # : 0 : exc_gap = MINVAL(ABS(DBLE (exc_ene(1:nstates))))
# # # # ]
952 [ # # # # : 0 : exc_maxene = MAXVAL(ABS(DBLE (exc_ene(1:nstates))))
# # # # ]
953 [ # # # # : 0 : temp = MAXVAL(ABS(AIMAG(exc_ene(1:nstates))))
# # # # ]
954 : :
955 : : write(msg,'(2(a,f7.2,2a),a,es9.2,2a)')&
956 : 0 : " First excitonic eigenvalue: ",exc_gap*Ha_eV, " [eV].",ch10,&
957 : 0 : " Last excitonic eigenvalue: ",exc_maxene*Ha_eV," [eV].",ch10,&
958 : 0 : " Largest imaginary part: ",temp*Ha_eV, " [eV] ",ch10
959 [ # # ]: 0 : call wrtout([std_out, ab_out], msg)
960 : :
961 [ # # ]: 0 : nene_printed = MIN(32*nsppol,nstates); if (prtvol>10) nene_printed = nstates
962 : :
963 : : ! This is not portable as the the eigenvalues calculated by ZGEEV are not sorted.
964 : : ! Even two subsequent calculations with the same input on the same machine
965 : : ! might produce different orderings. Might sort the eigenvalues though, just for printing.
966 : :
967 : 0 : write(msg,'(a,i0)')' Complex excitonic eigenvalues in eV up to n= ',nene_printed
968 : 0 : call wrtout(std_out,msg)
969 : :
970 [ # # ]: 0 : do it=0,(nene_printed-1)/4
971 [ # # # # ]: 0 : write(msg,'(8f10.5)') ( exc_ene(ii)*Ha_eV, ii=1+it*4,MIN(it*4+4,nene_printed) )
972 : 0 : call wrtout(std_out,msg)
973 : : end do
974 : :
975 [ # # ]: 0 : call wrtout(std_out,ch10//" Writing eigenvalues and eigenvectors on file "//TRIM(bseig_fname))
976 : :
977 [ # # ]: 0 : if (open_file(bseig_fname,msg,newunit=eig_unt,form="unformatted",action="write") /= 0) then
978 : 0 : ABI_ERROR(msg)
979 : : end if
980 : :
981 : : !YG : new version with lifetime
982 : 0 : do_ep_lifetime = .FALSE.
983 : 0 : write(eig_unt) do_ep_lifetime
984 : :
985 : 0 : write(eig_unt)exc_size,nstates
986 [ # # ]: 0 : write(eig_unt)CMPLX(exc_ene(1:nstates),kind=dp)
987 [ # # ]: 0 : do mi=1,nstates
988 : 0 : write(eig_unt) exc_rvect(:,mi)
989 : : end do
990 : :
991 [ # # ]: 0 : ABI_FREE(exc_ene)
992 : :
993 [ # # # # : 0 : ABI_MALLOC_OR_DIE(ovlp,(nstates,nstates), ierr)
# # # # #
# # # # #
# # # # ]
994 : :
995 : 0 : call wrtout(std_out,' Calculating overlap matrix... ')
996 : :
997 : : !do itp=1,nstates
998 : : ! do it=1,nstates
999 : : ! ovlp(it,itp) = xdotc(exc_size,exc_rvect(:,it),1,exc_rvect(:,itp),1)
1000 : : ! end do
1001 : : !end do
1002 : 0 : call xgemm("C","N",exc_size,nstates,nstates,cone,exc_rvect,exc_size,exc_rvect,exc_size,czero,ovlp,nstates)
1003 [ # # ]: 0 : ABI_FREE(exc_rvect)
1004 : :
1005 : 0 : call wrtout(std_out," Inverting overlap matrix... ")
1006 : :
1007 : : ! Version for generic complex matrix.
1008 : : !call xginv(ovlp,exc_size)
1009 : :
1010 : : ! The overlap is Hermitian definite positive.
1011 : 0 : call xhdp_invert("Upper",ovlp,nstates)
1012 : 0 : call hermitianize(ovlp,"Upper")
1013 : :
1014 [ # # ]: 0 : call wrtout(std_out,' Writing overlap matrix S^-1 on file: '//TRIM(bseig_fname))
1015 : :
1016 [ # # ]: 0 : do it=1,nstates
1017 [ # # ]: 0 : write(eig_unt) CMPLX(ovlp(:,it),kind=dp)
1018 : : end do
1019 : :
1020 [ # # ]: 0 : ABI_FREE(ovlp)
1021 : :
1022 : 0 : close(eig_unt)
1023 : :
1024 : 0 : 10 call xmpi_barrier(comm)
1025 : :
1026 [ # # ]: 0 : end subroutine exc_diago_coupling
1027 : : !!***
1028 : :
1029 : : !----------------------------------------------------------------------
1030 : :
1031 : : !!****f* m_exc_diago/exc_diago_coupling_hegv
1032 : : !! NAME
1033 : : !! exc_diago_coupling_hegv
1034 : : !!
1035 : : !! FUNCTION
1036 : : !! Calculate excitonic eigenvalues and eigenvectors by performing a direct diagonalization.
1037 : : !! of the non Hermitian excitonic Hamiltonian (resonant + coupling).
1038 : : !!
1039 : : !! INPUTS
1040 : : !! bseig_fname=The name of the output file.
1041 : : !! Bsp
1042 : : !! neh=Rank of the resonant block of the Hamiltoninan (equal to the rank of the coupling part)
1043 : : !! comm=MPI communicator.
1044 : : !! BS_files<excfiles>=Datatype storing names and files used in the Bethe-Salpeter code.
1045 : : !!
1046 : : !! OUTPUT
1047 : : !! Excitonic eigenvectors and eigenvalues are written to file BS_files%out_eig.
1048 : : !!
1049 : : !! SOURCE
1050 : :
1051 : 1 : subroutine exc_diago_coupling_hegv(Bsp,BS_files,Hdr_bse,prtvol,comm)
1052 : :
1053 : : !Arguments ------------------------------------
1054 : : !scalars
1055 : : integer,intent(in) :: comm,prtvol
1056 : : type(excparam),intent(in) :: BSp
1057 : : type(excfiles),intent(in) :: BS_files
1058 : : type(Hdr_type),intent(in) :: Hdr_bse
1059 : :
1060 : : !Local variables ------------------------------
1061 : : !scalars
1062 : : integer,parameter :: master=0
1063 : : integer(i8b) :: bsize_ham
1064 : : integer :: itype,il,iu,spin,row1,row2,pad_r1,pad_r2,neh1,neh2
1065 : : integer :: ii,exc_size,hreso_unt,hcoup_unt,eig_unt
1066 : : integer :: fform,neig_found,nstates, mi,it,nprocs,my_rank
1067 : : integer :: nene_printed,nsppol,row_sign,ierr
1068 : : real(dp) :: exc_gap,exc_maxene,abstol,vl,vu
1069 : : character(len=500) :: msg
1070 : : character(len=fnlen) :: reso_fname,coup_fname,bseig_fname
1071 : : logical :: use_scalapack,do_full_diago,diago_is_real, do_ep_lifetime
1072 : : !arrays
1073 : 1 : real(dp),allocatable :: exc_ene(:) !,test_ene(:)
1074 : 1 : complex(dp),allocatable :: exc_ham(:,:),exc_rvect(:,:),fmat(:,:),ovlp(:,:)
1075 : : #if defined HAVE_LINALG_SCALAPACK && defined HAVE_MPI_IO
1076 : : integer,parameter :: istwfk1=1
1077 : : integer :: amode,mpi_fh,tbloc,mene_found,mpi_err,my_nel,nsblocks
1078 : : integer :: iloc,jloc,iglob,jglob,etype,slk_mask_type,offset_err,el,rrs_kind,ccs_kind
1079 : : !integer :: max_r,max_c
1080 : : integer(XMPI_OFFSET_KIND) :: ehdr_offset,fmarker,my_offset
1081 : : integer :: gsub(2,2)
1082 : : logical,parameter :: is_fortran_file=.TRUE.
1083 : : complex(dp) :: ctmp
1084 : : integer,allocatable :: sub_block(:,:,:)
1085 : : integer,pointer :: myel2loc(:,:)
1086 : : complex(dp),allocatable :: tmp_cbuffer(:)
1087 : : character(50) :: uplo
1088 : : real(dp),external :: PDLAMCH
1089 : : type(slkmat_dp_t) :: Slk_F,Slk_Hbar,Slk_vec,Slk_ovlp !,Slk_tmp
1090 : : type(slk_processor_t) :: Slk_processor
1091 : : #endif
1092 : : !************************************************************************
1093 : :
1094 : 1 : my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm)
1095 : :
1096 : 1 : nsppol = Hdr_bse%nsppol
1097 [ - + ]: 1 : if (nsppol==2) then
1098 : 0 : ABI_WARNING("nsppol==2 is still under development!")
1099 : : end if
1100 : :
1101 : 1 : neh1 = BSp%nreh(1); neh2=neh1
1102 : : if (nsppol==2) neh2 = BSp%nreh(2)
1103 : :
1104 [ + + ]: 2 : exc_size = 2*SUM(Bsp%nreh)
1105 : 1 : nstates = Bsp%nstates
1106 : 1 : do_full_diago=(nstates==exc_size)
1107 : :
1108 : 1 : write(msg,'(a,i0)')'. Direct diagonalization of the full excitonic Hamiltonian, Matrix size= ',exc_size
1109 [ + + ]: 3 : call wrtout([std_out, ab_out], msg)
1110 : :
1111 : 1 : bseig_fname = BS_files%out_eig
1112 [ - + ]: 1 : if (BS_files%in_eig /= BSE_NOFILE) then
1113 : 0 : ABI_ERROR("BS_files%in_eig is defined!")
1114 : : end if
1115 : :
1116 [ + - ]: 1 : if (BS_files%in_hreso /= BSE_NOFILE) then
1117 : 1 : reso_fname = BS_files%in_hreso
1118 : : else
1119 : 0 : reso_fname = BS_files%out_hreso
1120 : : end if
1121 [ + - ]: 1 : call wrtout(std_out,' Reading resonant excitonic Hamiltonian from '//TRIM(reso_fname))
1122 : :
1123 [ - + ]: 1 : if (BS_files%in_hcoup /= BSE_NOFILE) then
1124 : 0 : coup_fname = BS_files%in_hcoup
1125 : : else
1126 : 1 : coup_fname = BS_files%out_hcoup
1127 : : end if
1128 [ + - ]: 1 : call wrtout(std_out,' Reading coupling excitonic Hamiltonian from '//TRIM(coup_fname))
1129 : :
1130 : : ! TODO: Reintegrate SCALAPACK: use new format
1131 : : ! --- !ERROR
1132 : : ! src_file: m_exc_diago.F90
1133 : : ! src_line: 1398
1134 : : ! mpi_rank: 0
1135 : : ! message: |
1136 : : ! SET_VIEW
1137 : : ! Other I/O error , error stack:
1138 : : ! ADIO_Set_view(48): **iobadoverlap displacements of filetype must be in a monotonically nondecreasing order
1139 : : ! ...
1140 : :
1141 : 1 : use_scalapack = .FALSE.
1142 : : #ifdef HAVE_LINALG_SCALAPACK
1143 : : ! This is alway false. I use this trick so that the second case below is always compiled
1144 : : ! to avoid regressions.
1145 : : use_scalapack = nprocs > 1 .and. nsppol > 5
1146 : : #endif
1147 : : !use_scalapack = .FALSE.
1148 : : !use_scalapack = .TRUE.
1149 : :
1150 [ - + ]: 1 : if (.not.use_scalapack .and. my_rank/=master) GOTO 10
1151 : :
1152 [ - + + - : 3 : ABI_MALLOC_OR_DIE(exc_ene,(exc_size), ierr)
+ - + - +
- - + -
+ ]
1153 : :
1154 : : SELECT CASE (use_scalapack)
1155 : :
1156 : : CASE (.FALSE.)
1157 : 1 : write(msg,'(a)')". Using LAPACK sequential version to solve FHv = ev with H positive definite. "
1158 [ + + ]: 3 : call wrtout([std_out, ab_out], msg)
1159 : :
1160 : 1 : bsize_ham = 2*dp*exc_size**2
1161 : 1 : write(msg,'(a,f9.2,a)')' Allocating full excitonic Hamiltonian. Memory requested: ',2*bsize_ham*b2Gb,' Gb. '
1162 : 1 : call wrtout(std_out, msg)
1163 : :
1164 [ - + + - : 4 : ABI_MALLOC_OR_DIE(exc_ham,(exc_size,exc_size), ierr)
- + + - +
- + - + -
- + - + ]
1165 [ - + + - : 3 : ABI_MALLOC_OR_DIE(fmat,(exc_size,exc_size), ierr)
- + + - +
- - + -
+ ]
1166 : :
1167 : : write(msg,'(3a,f8.1,3a,f8.1,a)')&
1168 : 1 : ' Allocating excitonic eigenvalues and eigenvectors. ',ch10,&
1169 : 1 : ' Memory-space requested: ',2*dp*exc_size*b2Gb,' Gb. ',ch10,&
1170 : 2 : ' Memory-space requested: ',bsize_ham*b2Gb,' Gb. '
1171 : 1 : call wrtout(std_out, msg)
1172 : :
1173 [ - + ]: 1 : if (open_file(reso_fname,msg,newunit=hreso_unt,form="unformatted",status="old",action="read") /= 0) then
1174 : 0 : ABI_ERROR(msg)
1175 : : end if
1176 : : !
1177 : : ! Read the header and perform consistency checks.
1178 : 1 : call exc_read_bshdr(hreso_unt,Bsp,fform,ierr)
1179 [ - + ]: 1 : ABI_CHECK(ierr==0,"Wrong header")
1180 : : !
1181 : : ! Construct Hbar = ( R C )
1182 : : ! ( C* R*)
1183 : : !
1184 : 1 : row_sign=+1; diago_is_real=(.not.BSp%have_complex_ene)
1185 : 1 : call exc_fullh_from_blocks(hreso_unt,"Resonant",nsppol,row_sign,diago_is_real,Bsp%nreh,exc_size,exc_ham)
1186 : 1 : close(hreso_unt)
1187 : :
1188 [ - + ]: 1 : if (open_file(coup_fname,msg,newunit=hcoup_unt,form="unformatted",status="old",action="read") /= 0) then
1189 : 0 : ABI_ERROR(msg)
1190 : : end if
1191 : : !
1192 : : ! Read the header and perform consistency checks.
1193 : 1 : call exc_read_bshdr(hcoup_unt,Bsp,fform,ierr)
1194 [ - + ]: 1 : ABI_CHECK(ierr==0,"Wrong header")
1195 : :
1196 : 1 : row_sign=+1; diago_is_real=(.not.BSp%have_complex_ene) ! not used here.
1197 : 1 : call exc_fullh_from_blocks(hcoup_unt,"Coupling",nsppol,row_sign,diago_is_real,Bsp%nreh,exc_size,exc_ham)
1198 : 1 : close(hcoup_unt)
1199 : :
1200 : : !#ifdef DEV_MG_DEBUG_THIS
1201 : : !write(666)exc_ham
1202 : : !#endif
1203 : : !
1204 : : ! Fill fmat = (1 0)
1205 : : ! (0 -1)
1206 [ + + + + ]: 590593 : fmat = czero
1207 [ + + ]: 2 : do spin=1,nsppol
1208 : 1 : pad_r1 = (spin-1)*Bsp%nreh(1)
1209 [ + + ]: 2 : pad_r2 = SUM(Bsp%nreh)
1210 [ - + ]: 1 : if (spin==2) pad_r2 = pad_r2 + Bsp%nreh(1)
1211 [ + + ]: 386 : do it=1,Bsp%nreh(spin)
1212 : 384 : row1 = it + pad_r1
1213 : 384 : row2 = it + pad_r2
1214 : 384 : fmat(row1,row1) = cone
1215 : 385 : fmat(row2,row2) = -cone
1216 : : end do
1217 : : end do
1218 : : !
1219 : : ! ==================================================
1220 : : ! ==== Solve generalized EV problem F H u = e u ====
1221 : : ! ==================================================
1222 : : ! The eigenvectors Z are normalized as follows: if ITYPE = 1 or 2, Z**T*B*Z = I; if ITYPE = 3, Z**T*inv(B)*Z = I.
1223 : : !
1224 : 1 : itype=2
1225 [ + - ]: 1 : if (do_full_diago) then
1226 : 1 : call wrtout(std_out," Full diagonalization with XHEGV... ")
1227 : 1 : call xhegv(itype,"Vectors","Upper",exc_size,fmat,exc_ham,exc_ene)
1228 : : else
1229 : 0 : call wrtout(std_out," Partial diagonalization with XHEGVX... ")
1230 [ # # # # : 0 : ABI_MALLOC_OR_DIE(exc_rvect,(exc_size,nstates), ierr)
# # # # #
# # # # #
# # # # ]
1231 : 0 : il=1; iu=1; abstol=zero
1232 : 0 : call xhegvx(itype,"Vectors","All","Upper",exc_size,fmat,exc_ham,vl,vu,il,iu,abstol,neig_found,exc_ene,exc_rvect,exc_size)
1233 : : end if
1234 : :
1235 [ - + ]: 1 : ABI_FREE(exc_ham)
1236 : :
1237 [ + - ]: 1 : if (do_full_diago) then
1238 [ - + + - : 4 : ABI_MALLOC(exc_rvect,(exc_size,nstates))
- + + - +
- + - - +
- + - + ]
1239 [ - + + + : 590594 : exc_rvect = fmat(:,1:nstates)
+ + ]
1240 : : end if
1241 : :
1242 [ - + ]: 1 : ABI_FREE(fmat)
1243 : :
1244 [ + - ]: 1 : call wrtout(std_out," Writing eigenvalues and eigenvectors on file: "//TRIM(bseig_fname))
1245 : :
1246 [ - + ]: 1 : if (open_file(bseig_fname,msg,newunit=eig_unt,form="unformatted",action="write") /= 0) then
1247 : 0 : ABI_ERROR(msg)
1248 : : end if
1249 : :
1250 : 1 : do_ep_lifetime = .FALSE.
1251 : 1 : write(eig_unt) do_ep_lifetime
1252 : 1 : write(eig_unt) exc_size, nstates
1253 [ + + ]: 769 : write(eig_unt) CMPLX(exc_ene(1:nstates),kind=dp)
1254 [ + + ]: 769 : do mi=1,nstates
1255 [ + + ]: 590593 : write(eig_unt) CMPLX(exc_rvect(:,mi),kind=dp)
1256 : : end do
1257 : :
1258 : : !#ifdef DEV_MG_DEBUG_THIS
1259 : : ! write(888)exc_rvect
1260 : : ! write(888)exc_ene
1261 : : !#endif
1262 : :
1263 [ - + + - : 4 : ABI_MALLOC_OR_DIE(ovlp,(nstates,nstates), ierr)
- + + - +
- + - + -
- + - + ]
1264 : :
1265 : 1 : call wrtout(std_out,' Calculating overlap matrix...')
1266 : :
1267 : 1 : call xgemm("C","N",exc_size,nstates,nstates,cone,exc_rvect,exc_size,exc_rvect,exc_size,czero,ovlp,nstates)
1268 [ - + ]: 1 : ABI_FREE(exc_rvect)
1269 : :
1270 : : !#ifdef DEV_MG_DEBUG_THIS
1271 : : !write(667)ovlp
1272 : : !#endif
1273 : :
1274 : 1 : call wrtout(std_out," Inverting overlap matrix... ")
1275 : : !
1276 : : ! The overlap is Hermitian definite positive.
1277 : 1 : call xhdp_invert("Upper",ovlp,nstates)
1278 : 1 : call hermitianize(ovlp,"Upper")
1279 : :
1280 : : ! Version for generic complex matrix.
1281 : : !call xginv(ovlp,nstates)
1282 : :
1283 : : !#ifdef DEV_MG_DEBUG_THIS
1284 : : !write(668,*)ovlp
1285 : : !#endif
1286 : :
1287 [ + - ]: 1 : call wrtout(std_out,' Writing overlap matrix O^-1 on file: '//TRIM(bseig_fname))
1288 [ + + ]: 769 : do it=1,nstates
1289 : 769 : write(eig_unt) ovlp(:,it)
1290 : : end do
1291 : :
1292 [ - + ]: 1 : ABI_FREE(ovlp)
1293 : 3 : close(eig_unt)
1294 : :
1295 : : CASE (.TRUE.)
1296 : :
1297 : : #if defined HAVE_LINALG_SCALAPACK && defined HAVE_MPI_IO
1298 : : !
1299 : : ! Init scaLAPACK matrix Hbar = ( R C )
1300 : : ! ( C* R*)
1301 : : ! Battle plan:
1302 : : ! Here the reading is complicated by the fact that R and C are stored on two different files
1303 : : ! and moreover the matrices are in packed storage mode.
1304 : : ! For initializing the local part of the resonant and anti-resonant block we have to allocate
1305 : : ! a temporary buffer. Then we read the buffer from file and the corresponding elements of the
1306 : : ! scaLAPACK matrix are initialized taking into account the symmetries of R (Hermitian)
1307 : : ! The same procedure is used to read the coupling and the anti-coupling part (Symmetric).
1308 : : !
1309 : : tbloc=50
1310 : : write(msg,'(2(a,i0))')". Using MPI-IO + scaLAPACK version with nprocs= ",nprocs,"; block size= ",tbloc
1311 : : call wrtout([std_out, ab_out], msg, do_flush=.True.)
1312 : : !
1313 : : ! Init scaLAPACK environment.
1314 : : call Slk_processor%init(comm)
1315 : : !
1316 : : ! Open the Resonant file with MPI-IO and skip the record.
1317 : : amode=MPI_MODE_RDONLY
1318 : :
1319 : : call MPI_FILE_OPEN(comm, reso_fname, amode, MPI_INFO_NULL, mpi_fh, mpi_err)
1320 : : msg = " MPI_IO error opening file: "//TRIM(reso_fname)
1321 : : ABI_CHECK_MPI(mpi_err,msg)
1322 : : !
1323 : : ! Skip the header and find the offset for reading the matrix.
1324 : : call exc_skip_bshdr_mpio(mpi_fh,xmpio_collective,ehdr_offset)
1325 : : !
1326 : : ! Read = ( R - )
1327 : : ! ( - R*)
1328 : : call Slk_Hbar%init(exc_size,exc_size,Slk_processor,istwfk1)
1329 : :
1330 : : nullify(myel2loc)
1331 : : nsblocks=nsppol
1332 : : ABI_MALLOC(sub_block,(2,2,nsblocks))
1333 : : ABI_CHECK(nsppol==1,"nsppol==2 not coded yet")
1334 : :
1335 : : call slk_single_fview_read_mask(Slk_Hbar,rrs_of_glob,offset_in_file,nsblocks,sub_block,my_nel,myel2loc,etype,slk_mask_type,&
1336 : : & offset_err,is_fortran_file)
1337 : :
1338 : : if (offset_err/=0) then
1339 : : write(msg,"(3a)")&
1340 : : & " Global position index cannot be stored in a standard Fortran integer ",ch10,&
1341 : : & " Excitonic matrix cannot be read with a single MPI-IO call."
1342 : : ABI_ERROR(msg)
1343 : : end if
1344 : :
1345 : : ! Shift the offset because the view starts at the fist matrix element!
1346 : : ! TODO should rationalize the treatment of the offset
1347 : : my_offset = ehdr_offset + xmpio_bsize_frm
1348 : : call MPI_FILE_SET_VIEW(mpi_fh, my_offset, etype, slk_mask_type, 'native', MPI_INFO_NULL, mpi_err)
1349 : : ABI_CHECK_MPI(mpi_err,"SET_VIEW")
1350 : :
1351 : : call MPI_TYPE_FREE(slk_mask_type,mpi_err)
1352 : : ABI_CHECK_MPI(mpi_err,"MPI_TYPE_FREE")
1353 : : !
1354 : : ! Read my portion of the R,-R* sublocks and store the values in a temporary buffer.
1355 : : ABI_MALLOC_OR_DIE(tmp_cbuffer,(my_nel), ierr)
1356 : :
1357 : : call xmpi_barrier(comm)
1358 : :
1359 : : call MPI_FILE_READ_ALL(mpi_fh, tmp_cbuffer, my_nel, MPI_DOUBLE_COMPLEX, MPI_STATUS_IGNORE, mpi_err)
1360 : : ABI_CHECK_MPI(mpi_err,"READ_ALL")
1361 : : !
1362 : : ! Symmetrize my Resonant part.
1363 : : do el=1,my_nel
1364 : : iloc = myel2loc(1,el)
1365 : : jloc = myel2loc(2,el)
1366 : : call Slk_Hbar%loc2glob(iloc,jloc,iglob,jglob)
1367 : : ctmp = tmp_cbuffer(el)
1368 : : if (iglob==jglob.and..not.Bsp%have_complex_ene) ctmp = DBLE(ctmp) ! Force the diagonal to be real.
1369 : : rrs_kind = rrs_of_glob(iglob,jglob,Slk_Hbar%size_global)
1370 : : if (rrs_kind==1.and.jglob<iglob) then ! Lower resonant
1371 : : ctmp = DCONJG(ctmp)
1372 : : else if (rrs_kind==-1.and.jglob>=iglob) then ! Lower Anti-resonant (Diagonal is included).
1373 : : ctmp = DCONJG(ctmp)
1374 : : end if
1375 : : Slk_Hbar%buffer_cplx(iloc,jloc) = ctmp
1376 : : end do
1377 : :
1378 : : ABI_FREE(tmp_cbuffer)
1379 : : ABI_FREE(myel2loc)
1380 : :
1381 : : call MPI_FILE_CLOSE(mpi_fh, mpi_err)
1382 : : ABI_CHECK_MPI(mpi_err,"FILE_CLOSE")
1383 : : !
1384 : : ! Read = ( - C)
1385 : : ! (-C* -)
1386 : : !
1387 : : call MPI_FILE_OPEN(comm, coup_fname, amode, MPI_INFO_NULL, mpi_fh, mpi_err)
1388 : : msg = " MPI_IO error opening file: "//TRIM(coup_fname)
1389 : : ABI_CHECK_MPI(mpi_err,msg)
1390 : : !
1391 : : ! Skip the header and find the offset for reading the matrix.
1392 : : call exc_skip_bshdr_mpio(mpi_fh,xmpio_collective,ehdr_offset)
1393 : :
1394 : : nullify(myel2loc)
1395 : : call slk_single_fview_read_mask(Slk_Hbar,ccs_of_glob,offset_in_file,nsblocks,sub_block,my_nel,myel2loc,etype,slk_mask_type,&
1396 : : & offset_err,is_fortran_file)
1397 : :
1398 : : ABI_FREE(sub_block)
1399 : :
1400 : : if (offset_err/=0) then
1401 : : write(msg,"(3a)")&
1402 : : & " Global position index cannot be stored in a standard Fortran integer ",ch10,&
1403 : : & " Excitonic matrix cannot be read with a single MPI-IO call."
1404 : : ABI_ERROR(msg)
1405 : : end if
1406 : : !
1407 : : ! Shift the offset because the view starts at the fist matrix element!
1408 : : ! TODO should rationalize the treatment of the offset so that the client code
1409 : : ! will automatically receive my_offset.
1410 : : my_offset = ehdr_offset + xmpio_bsize_frm
1411 : : call MPI_FILE_SET_VIEW(mpi_fh, my_offset, etype, slk_mask_type, 'native', MPI_INFO_NULL, mpi_err)
1412 : : ABI_CHECK_MPI(mpi_err,"SET_VIEW")
1413 : :
1414 : : call MPI_TYPE_FREE(slk_mask_type,mpi_err)
1415 : : ABI_CHECK_MPI(mpi_err,"MPI_TYPE_FREE")
1416 : : !
1417 : : ! Read my portion of the C-C* blocks and store the values in a temporary buffer.
1418 : : ABI_MALLOC_OR_DIE(tmp_cbuffer,(my_nel), ierr)
1419 : :
1420 : : call MPI_FILE_READ_ALL(mpi_fh, tmp_cbuffer, my_nel, MPI_DOUBLE_COMPLEX, MPI_STATUS_IGNORE, mpi_err)
1421 : : ABI_CHECK_MPI(mpi_err,"READ_ALL")
1422 : : !
1423 : : ! Symmetrize my coupling part.
1424 : : ! Coupling block is symmetric => No symmetrization of the lower triangle.
1425 : : do el=1,my_nel
1426 : : iloc = myel2loc(1,el)
1427 : : jloc = myel2loc(2,el)
1428 : : call Slk_Hbar%loc2glob(iloc, jloc, iglob, jglob)
1429 : : ccs_kind = ccs_of_glob(iglob,jglob,Slk_Hbar%size_global)
1430 : : ctmp = tmp_cbuffer(el)
1431 : : if (ccs_kind==-1) ctmp = DCONJG(ctmp) ! Anti-coupling (Diagonal is included).
1432 : : Slk_Hbar%buffer_cplx(iloc,jloc) = ctmp
1433 : : end do
1434 : :
1435 : : ABI_FREE(tmp_cbuffer)
1436 : : ABI_FREE(myel2loc)
1437 : :
1438 : : !max_r=20; max_c=10
1439 : : !call print_arr(Slk_Hbar%buffer_cplx,max_r=max_r,max_c=max_c,unit=std_out)
1440 : :
1441 : : !#ifdef DEV_MG_DEBUG_THIS
1442 : : ! ABI_MALLOC(exc_ham,(exc_size,exc_size))
1443 : : ! read(666)exc_ham
1444 : : !
1445 : : ! write(std_out,*)"Error Hbar: ",MAXVAL(ABS(exc_ham-Slk_Hbar%buffer_cplx))
1446 : : ! ABI_FREE(exc_ham)
1447 : : !#endif
1448 : :
1449 : : call MPI_FILE_CLOSE(mpi_fh, mpi_err)
1450 : : ABI_CHECK_MPI(mpi_err,"FILE_CLOSE")
1451 : : !
1452 : : ! Init scaLAPACK matrix F
1453 : : call Slk_F%init(exc_size,exc_size,Slk_processor,istwfk1)
1454 : : !
1455 : : ! Global F = (1 0)
1456 : : ! (0 -1)
1457 : : do jloc=1,Slk_F%size_local(2)
1458 : : do iloc=1,Slk_F%size_local(1)
1459 : : call Slk_F%loc2glob(iloc, jloc, iglob, jglob)
1460 : : if (iglob==jglob) then
1461 : : if (iglob<=SUM(Bsp%nreh)) then
1462 : : Slk_F%buffer_cplx(iloc,jloc) = cone
1463 : : else
1464 : : Slk_F%buffer_cplx(iloc,jloc) = -cone
1465 : : end if
1466 : : else
1467 : : Slk_F%buffer_cplx(iloc,jloc) = czero
1468 : : end if
1469 : : end do
1470 : : end do
1471 : : !
1472 : : ! ===========================================================
1473 : : ! ==== Solve generalized EV problem H u = F Hbar u = e u ====
1474 : : ! ===========================================================
1475 : : call Slk_vec%init(exc_size,exc_size,Slk_processor,istwfk1)
1476 : : !
1477 : : itype=2; vl=1; vu=1; il=1; iu=nstates
1478 : : abstol=zero !ABSTOL = PDLAMCH(comm,'U')
1479 : :
1480 : : !#if 1
1481 : : if (do_full_diago) then
1482 : : call slk_F%pzhegvx(itype,"Vectors","All","Upper",Slk_Hbar,vl,vu,il,iu,abstol,Slk_vec,mene_found,exc_ene)
1483 : : else
1484 : : ABI_WARNING("Partial diago is still under testing")
1485 : : call slk_F%pzhegvx(itype,"Vectors","Index","Upper",Slk_Hbar,vl,vu,il,iu,abstol,Slk_vec,mene_found,exc_ene)
1486 : : end if
1487 : : !#else
1488 : : ! call xhegv(itype,"Vectors","Upper",exc_size,Slk_F%buffer_cplx,Slk_Hbar%buffer_cplx,exc_ene)
1489 : : ! Slk_vec%buffer_cplx = Slk_F%buffer_cplx
1490 : : !#endif
1491 : :
1492 : : !#ifdef DEV_MG_DEBUG_THIS
1493 : : ! if (PRODUCT(Slk_Hbar%size_local) /= exc_size**2) then
1494 : : ! ABI_ERROR("Wrong size")
1495 : : ! end if
1496 : : !
1497 : : ! ABI_MALLOC(exc_ham,(exc_size,exc_size))
1498 : : ! read(888)exc_ham
1499 : : !
1500 : : ! write(std_out,*)"Error rvec: ",MAXVAL(ABS(exc_ham-Slk_vec%buffer_cplx))
1501 : : ! ABI_FREE(exc_ham)
1502 : : !
1503 : : ! ABI_MALLOC(test_ene,(exc_size))
1504 : : ! read(888)test_ene
1505 : : ! write(std_out,*)"Error ene: ",MAXVAL(ABS(exc_ene-test_ene))
1506 : : ! ABI_FREE(test_ene)
1507 : : !#endif
1508 : :
1509 : : call Slk_F%free()
1510 : : call Slk_Hbar%free()
1511 : :
1512 : : call wrtout(std_out,ch10//" Writing eigenvalues and eigenvectors on file: "//TRIM(bseig_fname))
1513 : : !
1514 : : ! Open the file with Fortran-IO to write the Header.
1515 : : if (my_rank==master) then
1516 : : if (open_file(bseig_fname,msg,newunit=eig_unt,form="unformatted",action="write") /= 0) then
1517 : : ABI_ERROR(msg)
1518 : : end if
1519 : :
1520 : : write(eig_unt) exc_size,nstates
1521 : : write(eig_unt) CMPLX(exc_ene(1:nstates),kind=dp)
1522 : : !do mi=1,exc_size
1523 : : ! write(eig_unt) CMPLX(exc_rvect(:,mi),kind=dp)
1524 : : !end do
1525 : : close(eig_unt)
1526 : : end if
1527 : : !
1528 : : ! Open the file with MPI-IO and write the distributed eigevectors.
1529 : : call xmpi_barrier(comm)
1530 : : amode=MPI_MODE_RDWR
1531 : : call MPI_FILE_OPEN(comm, bseig_fname, amode, MPI_INFO_NULL, mpi_fh, mpi_err)
1532 : : ABI_CHECK_MPI(mpi_err,"FILE_OPEN: "//TRIM(bseig_fname))
1533 : : !
1534 : : ! Skip the header and find the offset for writing the matrix.
1535 : : ehdr_offset=0
1536 : : !call hdr_mpio_skip(mpi_fh,fform,ehdr_offset)
1537 : :
1538 : : call xmpio_read_frm(mpi_fh,ehdr_offset,xmpio_collective,fmarker,mpi_err)
1539 : : write(std_out,*)" fmarker1 = ",fmarker
1540 : : call xmpio_read_frm(mpi_fh,ehdr_offset,xmpio_collective,fmarker,mpi_err)
1541 : : write(std_out,*)" fmarker2 = ",fmarker
1542 : :
1543 : : write(std_out,*)" Writing nstates ",nstates
1544 : : gsub(:,1) = (/1,1/)
1545 : : gsub(:,2) = (/exc_size,nstates/)
1546 : : call slk_write(Slk_vec,"All",is_fortran_file,mpi_fh=mpi_fh,offset=ehdr_offset,glob_subarray=gsub)
1547 : :
1548 : : call wrtout(std_out,' Calculating overlap matrix... ')
1549 : : if (.not.do_full_diago) then
1550 : : ABI_ERROR(" Init of Slk_ovlp is wrong")
1551 : : end if
1552 : :
1553 : : call Slk_ovlp%init(exc_size,exc_size,Slk_processor,istwfk1)
1554 : :
1555 : : ! Calculate the overlap matrix.
1556 : : ! FIXME
1557 : : ! The ESLL manual says that "matrices matrix1 and matrix2 must have no common elements;
1558 : : ! otherwise, results are unpredictable."
1559 : : ! However the official scaLAPACK documentation does not report this (severe) limitation.
1560 : :
1561 : : !call Slk_tmp%init(exc_size,exc_size,Slk_processor,istwfk1)
1562 : : !Slk_tmp%buffer_cplx = Slk_vec%buffer_cplx
1563 : : !call slk_pgemm("C","N",Slk_tmp,cone,Slk_vec,czero,Slk_ovlp)
1564 : : !call Slk_tmp%free()
1565 : :
1566 : : call slk_pgemm("C","N",Slk_vec,cone,Slk_vec,czero,Slk_ovlp)
1567 : :
1568 : : !#ifdef DEV_MG_DEBUG_THIS
1569 : : ! ABI_MALLOC(exc_ham,(exc_size,exc_size))
1570 : : ! read(667)exc_ham
1571 : : !
1572 : : ! write(std_out,*)"Error Ovlp: ",MAXVAL(ABS(exc_ham-Slk_ovlp%buffer_cplx))
1573 : : ! !Slk_ovlp%buffer_cplx = exc_ham
1574 : : !#endif
1575 : :
1576 : : !max_r=20; max_c=10
1577 : : !call print_arr(Slk_ovlp%buffer_cplx,max_r=max_r,max_c=max_c,unit=std_out)
1578 : :
1579 : : call Slk_vec%free()
1580 : :
1581 : : call wrtout(std_out," Inverting overlap matrix... ")
1582 : : uplo="Upper"
1583 : :
1584 : : !#if 0
1585 : : !!DEBUG
1586 : : ! call xhdp_invert(uplo,Slk_ovlp%buffer_cplx,exc_size)
1587 : : !
1588 : : ! !call slk_symmetrize(Slk_ovlp,uplo,"Hermitian")
1589 : : ! call hermitianize(Slk_ovlp%buffer_cplx,uplo)
1590 : : !
1591 : : ! exc_ham = MATMUL(exc_ham,Slk_ovlp%buffer_cplx)
1592 : : ! do it=1,exc_size
1593 : : ! exc_ham(it,it) = exc_ham(it,it) - cone
1594 : : ! end do
1595 : : !
1596 : : ! write(std_out,*)"Error Inversion: ",MAXVAL(ABS(exc_ham))
1597 : : ! ABI_FREE(exc_ham)
1598 : : !!END DEBUG
1599 : : !
1600 : : !#else
1601 : : ! call Slk_ovlp%hpd_invert(uplo)
1602 : : ! call hermitianize(Slk_ovlp%buffer_cplx,uplo)
1603 : : ! !call slk_symmetrize(Slk_ovlp,uplo,"Hermitian")
1604 : :
1605 : : call Slk_ovlp%invert() ! Version for generic complex matrix.
1606 : : !#endif
1607 : :
1608 : : !#ifdef DEV_MG_DEBUG_THIS
1609 : : ! ABI_MALLOC(exc_ham,(exc_size,exc_size))
1610 : : ! read(668)exc_ham
1611 : : ! write(std_out,*)"Error in Inv Ovlp: ",MAXVAL(ABS(exc_ham-Slk_ovlp%buffer_cplx))
1612 : : !
1613 : : ! !exc_ham = exc_ham-Slk_ovlp%buffer_cplx
1614 : : ! !do it=1,exc_size
1615 : : ! ! if ( MAXVAL(ABS(exc_ham(:,it))) > 0.1 ) write(std_out,*)"it: ",it,exc_ham(:,it)
1616 : : ! !end do
1617 : : !
1618 : : ! !Slk_ovlp%buffer_cplx = exc_ham
1619 : : ! ABI_FREE(exc_ham)
1620 : : !
1621 : : ! !write(std_out,*)"MAX ERR",MAXVAL(ABS(Slk_ovlp%buffer_cplx - TRANSPOSE(DCONJG(Slk_ovlp%buffer_cplx))))
1622 : : !#endif
1623 : :
1624 : : call wrtout(std_out,' Writing overlap matrix S^-1 on file: '//TRIM(bseig_fname))
1625 : :
1626 : : call slk_write(Slk_ovlp,"All",is_fortran_file,mpi_fh=mpi_fh,offset=ehdr_offset)
1627 : :
1628 : : call MPI_FILE_CLOSE(mpi_fh, mpi_err)
1629 : : ABI_CHECK_MPI(mpi_err,"FILE_CLOSE")
1630 : :
1631 : : call Slk_ovlp%free()
1632 : : call Slk_processor%free()
1633 : : #else
1634 : : ABI_BUG("You should not be here!")
1635 : : #endif
1636 : :
1637 : : END SELECT
1638 : :
1639 [ + - - + : 770 : exc_gap = MINVAL(ABS(exc_ene(1:nstates)))
- - + + ]
1640 [ + - - + : 770 : exc_maxene = MAXVAL(ABS(exc_ene(1:nstates)))
- - + + ]
1641 : :
1642 : : write(msg,'(2(a,f7.2,2a))')&
1643 : 1 : " First excitonic eigenvalue: ",exc_gap*Ha_eV, " [eV].",ch10,&
1644 : 2 : " Last excitonic eigenvalue: ",exc_maxene*Ha_eV," [eV].",ch10
1645 [ + + ]: 3 : call wrtout([std_out, ab_out], msg)
1646 : :
1647 [ - + ]: 1 : nene_printed = MIN(32*nsppol,nstates); if (prtvol>10) nene_printed = nstates
1648 : 1 : write(msg,'(a,i0)')' Complex excitonic eigenvalues in eV up to n= ',nene_printed
1649 : 1 : call wrtout(std_out, msg)
1650 : :
1651 [ + + ]: 9 : do it=0,(nene_printed-1)/4
1652 [ + - + + ]: 40 : write(msg,'(4f10.5)') ( exc_ene(ii)*Ha_eV, ii=1+it*4,MIN(it*4+4,nene_printed) )
1653 : 9 : call wrtout(std_out, msg)
1654 : : end do
1655 : :
1656 [ - + ]: 1 : ABI_FREE(exc_ene)
1657 : :
1658 : 1 : 10 call xmpi_barrier(comm)
1659 : :
1660 [ - + ]: 1 : end subroutine exc_diago_coupling_hegv
1661 : : !!***
1662 : :
1663 : : !----------------------------------------------------------------------
1664 : :
1665 : : end module m_exc_diago
1666 : : !!***
|
