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