Line data Source code
1 : !!****m* ABINIT/m_hexc
2 : !! NAME
3 : !! m_hexc
4 : !!
5 : !! FUNCTION
6 : !! Module for excitonic hamiltonian for Haydock
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2014-2026 ABINIT group (M.Giantomassi, Y. Gillet)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
16 : #if defined HAVE_CONFIG_H
17 : #include "config.h"
18 : #endif
19 :
20 : #include "abi_common.h"
21 :
22 : MODULE m_hexc
23 :
24 : use defs_basis
25 : use m_abicore
26 : use m_bs_defs
27 : use m_xmpi
28 : use m_errors
29 : use m_haydock_io
30 : use m_linalg_interfaces
31 : use netcdf
32 : use m_nctk
33 :
34 : use m_time, only : timab
35 : use m_fstrings, only : indent, strcat, sjoin, itoa
36 : use defs_datatypes, only : pseudopotential_type
37 : use m_hide_blas, only : xdotc, xgemv
38 : use m_numeric_tools, only : print_arr, symmetrize, hermitianize, wrap2_pmhalf
39 : use m_crystal, only : crystal_t
40 : use m_bz_mesh, only : kmesh_t, findqg0
41 : use m_double_grid, only : double_grid_t, get_kpt_from_indices_coarse, compute_corresp
42 : use m_wfd, only : wfdgw_t
43 : use m_bse_io, only : exc_read_rcblock, exc_write_optme, exc_ham_ncwrite
44 : use m_pawtab, only : pawtab_type
45 : use m_vcoul, only : vcoul_t
46 : use m_bseinterp, only : interpolator_t
47 : use m_ebands, only : ebands_t
48 :
49 : implicit none
50 :
51 : private
52 : !!***
53 :
54 : !!****t* m_haydock/hexc_t
55 : !! NAME
56 : !! hexc_t
57 : !!
58 : !! FUNCTION
59 : !! Store the excitonic hamiltonian and other related information
60 : !!
61 : !! SOURCE
62 :
63 : type, public :: hexc_t
64 :
65 : !scalars
66 : integer :: comm
67 : ! MPI communicator
68 :
69 : integer :: hsize_coarse
70 : ! Size of the coarse hamiltonian
71 :
72 : integer :: hsize
73 : ! Size of the hamiltonian and the kets
74 : ! (= hsize_coarse without interpolation, = hsize_dense with interpolation)
75 :
76 : integer :: nbz
77 : ! Number of kpoints for the full problem
78 : ! (= nbz_coarse without interpolation, = nbz_dense with interpolation)
79 :
80 : integer :: my_t1
81 : ! Lower limit of MPI paral
82 :
83 : integer :: my_t2
84 : ! Upper limit of MPI paral
85 :
86 : integer :: my_nt
87 : ! Number of transitions treat by node
88 : ! = my_t2 - my_t1 + 1
89 :
90 : integer :: nbnd_coarse
91 : ! Product of number of bands conduction X valence
92 :
93 : ! Pointers to data that are already in memory
94 : type(excparam),pointer :: bsp => null()
95 : ! parameters for BS
96 :
97 : type(excfiles),pointer :: bs_files => null()
98 : ! files for BSE
99 :
100 : type(crystal_t),pointer :: crystal => null()
101 : ! crystal info
102 :
103 : type(kmesh_t),pointer :: kmesh_coarse => null()
104 : ! kmesh of the coarse mesh
105 :
106 : type(kmesh_t),pointer :: kmesh => null()
107 : ! kmesh of the full problem
108 :
109 : type(ebands_t),pointer :: ks_bst => null()
110 : type(ebands_t),pointer :: qp_bst => null()
111 : ! band structures of the full problem
112 :
113 : type(wfdgw_t),pointer :: wfd_coarse => null()
114 : ! Wfd of the coarse problem
115 :
116 : type(wfdgw_t),pointer :: wfd => null()
117 : ! wfd of the full problem
118 :
119 : !arrays
120 : complex(dp),allocatable :: hreso(:,:)
121 : ! Resonant part of the hamiltonian
122 :
123 : complex(dp),allocatable :: hcoup(:,:)
124 : ! Coupling part of the hamiltonian
125 :
126 : complex(dp),allocatable :: diag_coarse(:)
127 : ! Diagonal part of the hamiltonian with transition energies
128 :
129 : contains
130 : procedure :: init => hexc_init ! Construct the object
131 : procedure :: free => hexc_free ! Free memory
132 : procedure :: build_hinterp => hexc_build_hinterp ! Interpolate the Hamiltonian and store it in memory
133 : procedure :: matmul_tda => hexc_matmul_tda ! Matrix-vector multiplication (TDA)
134 : procedure :: matmul_full => hexc_matmul_full ! Matrix-vector multiplication (TDA + Coupling)
135 : procedure :: matmul_elphon => hexc_matmul_elphon ! Matrix-vector multiplication (TDA + elphon)
136 : end type hexc_t
137 : !!***
138 :
139 : !----------------------------------------------------------------------
140 :
141 : !!****t* m_hexc/hexc_interp_t
142 : !! NAME
143 : !! hexc_interp_t
144 : !!
145 : !! FUNCTION
146 : !! Store information about interpolation of excitonic hamiltonian
147 : !!
148 : !! SOURCE
149 :
150 : type,public :: hexc_interp_t
151 :
152 : !scalars
153 : integer :: hsize_dense
154 : ! Size of the dense hamiltonian
155 :
156 : real(dp) :: m3_width
157 : ! Width of the region where M3 is applied instead of M1
158 :
159 : type(interpolator_t) :: interpolator
160 : ! Interpolator containing overlaps and interpolation info
161 :
162 : ! Pointers to datatypes that are already in memory
163 : type(kmesh_t),pointer :: kmesh_dense => null()
164 : ! kmesh of the dense mesh
165 :
166 : type(vcoul_t),pointer :: vcp_dense => null()
167 : ! coulomb interaction on the dense mesh
168 :
169 : !arrays
170 : integer,allocatable :: kdense2div(:)
171 : ! kdense2div(nbz_dense)
172 : ! Index of kpoint -> Index of division
173 :
174 : integer,allocatable :: div2kdense(:,:)
175 : ! div2kdense(nbz_coarse,ndiv)
176 : ! Index of kpoint coarse + Index of division -> Index of kdense
177 :
178 : complex(dp),allocatable :: diag_dense(:)
179 : ! diag_dense(hsize_dense)
180 : ! Diagonal part of the dense hamiltonian
181 :
182 : complex(dp),allocatable :: hinterp(:,:)
183 : ! hinterp(hsize_dense,hsize_dense)
184 : ! Interpolated hamiltonian
185 :
186 : complex(dp),allocatable :: all_hmat(:,:)
187 : ! all_hmat,(hsize,hsize))
188 : ! Coarse excitonic matrix in a format suitable for interpolation in k-space
189 :
190 : complex(dp),allocatable :: all_acoeffs(:,:)
191 : ! all_acoeffs(hsize,hsize))
192 : ! a coefficients in a format suitable for interpolation in k-space
193 :
194 : complex(dp),allocatable :: all_bcoeffs(:,:)
195 : ! all_bcoeffs(hsize,hsize))
196 : ! b coefficients in a format suitable for interpolation in k-space
197 :
198 : complex(dp),allocatable :: all_ccoeffs(:,:)
199 : ! all_ccoeffs(hsize,hsize))
200 : ! c coefficients in a format suitable for interpolation in k-space
201 :
202 : contains
203 : procedure :: init => hexc_interp_init ! Construct the object for interpolated ham
204 : procedure :: free => hexc_interp_free ! Free memory for interpolated ham
205 : end type hexc_interp_t
206 : !!***
207 :
208 : !----------------------------------------------------------------------
209 :
210 : CONTAINS !=======================================================================
211 : !!***
212 :
213 : !!****f* m_hexc/hexc_init
214 : !! NAME
215 : !! hexc_init
216 : !!
217 : !! FUNCTION
218 : !! Construct the hexc object
219 : !!
220 : !! INPUTS
221 : !! BSp<excparam>=Parameters of BS
222 : !! BS_files<excparam>=Files for BS
223 : !! Cryst<crystal_t>=Info on the crystalline structure
224 : !! Kmesh_coarse<kmesh_t>=Kmesh info
225 : !! Wfd_coarse<wfdgw_t>=Wavefunction descriptor
226 : !! KS_BSt<ebands_t>=Kohn-Sham band structure
227 : !! QP_BSt<ebands_t>=Quasi-Particle band structure
228 : !! comm=communicator
229 : !!
230 : !! OUTPUT
231 : !! hexc<hexc_t>=Excitonic Hamiltonian
232 : !!
233 : !! SOURCE
234 :
235 23 : subroutine hexc_init(hexc, BSp, BS_files, Cryst, Kmesh_coarse, Wfd_coarse, KS_BSt, QP_BSt, comm)
236 :
237 : !Arguments ---------------------------
238 : !scalars
239 : class(hexc_t),intent(inout) :: hexc
240 : integer,intent(in) :: comm
241 : type(excparam),intent(in),target :: BSp
242 : type(excfiles),intent(in),target :: BS_files
243 : type(crystal_t),intent(in),target :: Cryst
244 : type(kmesh_t),intent(in),target :: Kmesh_coarse
245 : type(wfdgw_t),intent(in),target :: Wfd_coarse
246 : type(ebands_t),intent(in),target :: KS_BSt, QP_BSt
247 :
248 : !Local variables ---------------------
249 : !scalars
250 : integer :: ierr,ncid,ncerr, max_r, max_c, hsize, spin, spad, itt ! For diagonal !
251 : logical :: is_resonant, diago_is_real, use_mpio=.FALSE.
252 : character(len=fnlen) :: hreso_fname, hcoup_fname
253 : !character(len=500) :: msg
254 : !arrays
255 23 : complex(dp),allocatable :: test(:,:)
256 : !*****************************************************************************
257 :
258 23 : hexc%bsp => BSp
259 23 : hexc%bs_files => BS_files
260 23 : hexc%crystal => Cryst
261 23 : hexc%kmesh_coarse => Kmesh_coarse
262 :
263 23 : hexc%comm = comm
264 46 : hsize = SUM(BSp%nreh)
265 :
266 23 : hexc%ks_bst => KS_BSt
267 23 : hexc%qp_bst => QP_BSt
268 23 : hexc%wfd => Wfd_coarse
269 23 : hexc%wfd_coarse => Wfd_coarse
270 23 : hexc%kmesh => Kmesh_coarse
271 23 : hexc%hsize_coarse = hsize
272 23 : hexc%hsize = hsize
273 23 : hexc%nbz = Kmesh_coarse%nbz
274 :
275 23 : hexc%nbnd_coarse = BSp%maxnbndv*BSp%maxnbndc
276 :
277 : ! Divide the columns of the Hamiltonian among the nodes.
278 23 : call xmpi_split_work(hsize,comm,hexc%my_t1,hexc%my_t2)
279 :
280 23 : hexc%my_nt = hexc%my_t2 - hexc%my_t1 + 1
281 23 : ABI_CHECK(hexc%my_nt>0,"found processor with 0 rows")
282 :
283 92 : ABI_MALLOC_OR_DIE(hexc%hreso,(hsize,hexc%my_t1:hexc%my_t2), ierr)
284 :
285 : ! Read the resonant block from file.
286 23 : if (BS_files%in_hreso /= BSE_NOFILE) then
287 5 : hreso_fname = BS_files%in_hreso
288 : else
289 18 : hreso_fname = BS_files%out_hreso
290 : end if
291 :
292 23 : is_resonant=.TRUE.; diago_is_real=(.not.BSp%have_complex_ene)
293 : call exc_read_rcblock(hreso_fname,Bsp,is_resonant,diago_is_real,BSp%nsppol,BSp%nreh,hsize,&
294 23 : & hexc%my_t1,hexc%my_t2,hexc%hreso,use_mpio,comm)
295 :
296 : !BEGIN DEBUG
297 23 : if (use_mpio) then
298 0 : ABI_WARNING("Testing MPI-IO routines")
299 0 : ABI_MALLOC_OR_DIE(test,(hsize,hexc%my_t1:hexc%my_t2), ierr)
300 0 : diago_is_real=(.not.BSp%have_complex_ene)
301 : call exc_read_rcblock(hreso_fname,Bsp,is_resonant,diago_is_real,Bsp%nsppol,Bsp%nreh,hsize,&
302 0 : hexc%my_t1,hexc%my_t2,test,.FALSE.,comm)
303 0 : test = test-hexc%hreso
304 0 : write(std_out,*)"DEBUG: Diff MPI-IO - Fortran ",MAXVAL(ABS(test))
305 0 : max_r=20; max_c=10
306 0 : write(std_out,*)" **** Testing resonant block **** "
307 0 : call print_arr([std_out], test,max_r=max_r,max_c=max_c)
308 0 : if (BSp%nsppol==2) then
309 0 : write(std_out,*)" **** D down down ****"
310 0 : call print_arr([std_out], test(hsize/2+1:,hsize/2+1:),max_r=max_r,max_c=max_c)
311 0 : write(std_out,*)" **** V up down ****"
312 0 : call print_arr([std_out], test(1:hsize/2,hsize/2+1:),max_r=max_r,max_c=max_c)
313 0 : write(std_out,*)" **** V down up ****"
314 0 : call print_arr([std_out], test(hsize/2+1:,1:hsize/2),max_r=max_r,max_c=max_c)
315 : end if
316 0 : ABI_FREE(test)
317 : end if
318 : !END DEBUG
319 :
320 : !
321 : ! Read coupling block.
322 23 : if (BSp%use_coupling>0) then
323 0 : ABI_CHECK(.not. Bsp%use_interp,"interpolation with coupling not coded!")
324 0 : if (BS_files%in_hcoup /= BSE_NOFILE) then
325 0 : hcoup_fname = BS_files%in_hcoup
326 : else
327 0 : hcoup_fname = BS_files%out_hcoup
328 : end if
329 :
330 0 : ABI_MALLOC_OR_DIE(hexc%hcoup,(hsize,hexc%my_t1:hexc%my_t2), ierr)
331 0 : is_resonant=.FALSE.; diago_is_real=.FALSE.
332 : call exc_read_rcblock(hcoup_fname,Bsp,is_resonant,diago_is_real,BSp%nsppol,BSp%nreh,hsize,&
333 0 : hexc%my_t1,hexc%my_t2,hexc%hcoup,use_mpio,comm)
334 : !call symmetrize(hcoup,"ALL")
335 :
336 0 : if (use_mpio) then
337 0 : ABI_WARNING("Testing MPI-IO routines")
338 0 : ABI_MALLOC_OR_DIE(test,(hsize,hexc%my_t1:hexc%my_t2), ierr)
339 0 : diago_is_real=.FALSE.
340 0 : call exc_read_rcblock(hcoup_fname,Bsp,is_resonant,diago_is_real,BSp%nsppol,Bsp%nreh,hsize,hexc%my_t1,hexc%my_t2,test,.FALSE.,comm)
341 0 : test = test-hexc%hcoup
342 0 : write(std_out,*)"DEBUG: Diff MPI-IO - Fortran ",MAXVAL(ABS(test))
343 0 : max_r=20; max_c=10
344 0 : write(std_out,*)" **** Testing coupling block **** "
345 0 : call print_arr([std_out], test,max_r=max_r,max_c=max_c)
346 0 : if (BSp%nsppol==2) then
347 0 : write(std_out,*)" **** D down down ****"
348 0 : call print_arr([std_out], test(hsize/2+1:,hsize/2+1:),max_r=max_r,max_c=max_c)
349 0 : write(std_out,*)" **** V up down ****"
350 0 : call print_arr([std_out], test(1:hsize/2,hsize/2+1:),max_r=max_r,max_c=max_c)
351 0 : write(std_out,*)" **** V down up ****"
352 0 : call print_arr([std_out], test(hsize/2+1:,1:hsize/2),max_r=max_r,max_c=max_c)
353 : end if
354 0 : ABI_FREE(test)
355 : end if
356 : end if
357 :
358 23 : if(BSp%prt_ncham .or. BSp%use_interp) then
359 : ! I want to store the diagonal part for future use (printing or interpolation) !
360 12 : ABI_MALLOC(hexc%diag_coarse,(hexc%hsize_coarse))
361 4 : spad=0
362 8 : do spin=1,BSp%nsppol
363 4 : if(spin==2) spad=BSp%nreh(1)
364 392 : do itt=1,BSp%nreh(spin) ! 1 is for spin 1
365 388 : hexc%diag_coarse(spad+itt) = Bsp%Trans(itt,spin)%en
366 : end do
367 : end do
368 :
369 4 : if (BSp%prt_ncham) then
370 0 : ncerr = nctk_open_create(ncid, trim(hexc%BS_files%out_basename)//"_HEXC.nc", xmpi_comm_self)
371 0 : NCF_CHECK_MSG(ncerr, "Creating HEXC file")
372 : call exc_ham_ncwrite(ncid, hexc%Kmesh_coarse, hexc%BSp, hexc%hsize_coarse, hexc%BSp%nreh, &
373 0 : & hexc%BSp%vcks2t,hexc%hreso,hexc%diag_coarse)
374 0 : NCF_CHECK(nf90_close(ncid))
375 : end if
376 : end if
377 :
378 23 : end subroutine hexc_init
379 : !!***
380 :
381 : !-------------------------------------------------------------------
382 :
383 : !!****f* m_hexc/hexc_interp_init
384 : !! NAME
385 : !! hexc_interp_init
386 : !!
387 : !! FUNCTION
388 : !! Construct the hexc_interp object
389 : !!
390 : !! INPUTS
391 : !! hexc<hexc_t>=Excitonic hamiltonian
392 : !! Kmesh_dense<kmesh_t>=Kmesh info
393 : !! Vcp_dense<vcoul_t>=Dense mesh info about coulomb
394 : !! double_grid<double_grid_t>=Link between dense and coarse mesh
395 : !! Wfd_dense<wfdgw_t>=Wavefunction descriptor
396 : !! KS_BSt_dense<ebands_t>=Kohn-Sham band structure
397 : !! QP_BSt_dense<ebands_t>=Quasi-Particle band structure
398 : !! Psps <type(pseudopotential_type)>=variables related to pseudopotentials.
399 : !! Pawtab(Cryst%ntypat*usepaw)<pawtab_type>=PAW tabulated starting data.
400 : !! comm=communicator
401 : !!
402 : !! OUTPUT
403 : !! hexc_i<hexc_interp_t>=Interpolated excitonic hamiltonian
404 : !!
405 : !! SIDE EFFECTS
406 : !! hexc
407 : !! Will be modified so that the size of the problem is full interpolated hamiltonian
408 : !! Wfd, Wfd_dense
409 : !! The memory might be modified by computing wavefunctions
410 : !!
411 : !! SOURCE
412 :
413 4 : subroutine hexc_interp_init(hexc_i, hexc, m3_width, method, Kmesh_dense, Vcp_dense, &
414 4 : double_grid, Wfd_dense, KS_BSt_dense, QP_BSt_dense, Psps, Pawtab)
415 :
416 : !Arguments ---------------------------
417 : !scalars
418 : class(hexc_interp_t),intent(inout) :: hexc_i
419 : integer,intent(in) :: method
420 : real(dp),intent(in) :: m3_width
421 : type(hexc_t),intent(inout) :: hexc
422 : type(double_grid_t),intent(in),target :: double_grid
423 : type(wfdgw_t),intent(inout),target :: Wfd_dense !, Wfd
424 : type(kmesh_t),intent(in),target :: Kmesh_dense
425 : type(pseudopotential_type),intent(in) :: Psps
426 : type(vcoul_t),intent(in),target :: Vcp_dense
427 : type(ebands_t),intent(in),target :: KS_BSt_dense, QP_BSt_dense
428 : !arrays
429 : type(pawtab_type),intent(in) :: Pawtab(hexc%crystal%ntypat*hexc%Wfd_coarse%usepaw)
430 :
431 : !Local variables ---------------------
432 : !scalars
433 : integer,parameter :: spin1 = 1
434 : integer :: nsppol,ierr,ii,itt,nproc, my_rank,hsize
435 : real(dp),parameter :: threshold = 0.1_dp
436 4 : type(excparam) :: BSp
437 : logical :: is_resonant, diago_is_real, use_mpio
438 : !arrays
439 : character(len=fnlen) :: tmpfname, hreso_fname
440 : !*****************************************************************************
441 :
442 4 : BSp = hexc%bsp
443 4 : ABI_CHECK(BSp%nsppol == 1,"nsppol > 1 not implemented yet")
444 :
445 4 : hsize = hexc%hsize_coarse
446 :
447 4 : nproc = xmpi_comm_size(hexc%comm); my_rank= xmpi_comm_rank(hexc%comm)
448 4 : nsppol = hexc%Bsp%nsppol
449 :
450 4 : ABI_CHECK(nproc == 1,"Parallelization not available in interpolation")
451 :
452 4 : hexc_i%m3_width = m3_width
453 :
454 4 : hexc_i%kmesh_dense => Kmesh_dense
455 4 : hexc_i%vcp_dense => Vcp_dense
456 8 : hexc_i%hsize_dense = SUM(BSp%nreh_interp)
457 4 : hexc%hsize = hexc_i%hsize_dense
458 4 : hexc%nbz = Kmesh_dense%nbz
459 :
460 4 : hexc%ks_bst => KS_BSt_dense
461 4 : hexc%qp_bst => QP_BSt_dense
462 4 : hexc%wfd => Wfd_dense
463 4 : hexc%kmesh => Kmesh_dense
464 :
465 : ! No parallelization !
466 4 : hexc%my_t1 = 1
467 4 : hexc%my_t2 = hexc_i%hsize_dense
468 :
469 : ! Initialize the interpolator
470 : call hexc_i%interpolator%init(double_grid, Wfd_dense, hexc%Wfd_coarse, Kmesh_dense, &
471 4 : hexc%Kmesh_coarse, hexc%BSp, hexc%crystal, Psps, Pawtab, method)
472 :
473 4 : if (BSp%sum_overlaps) call hexc_i%interpolator%normalize()
474 :
475 12 : ABI_MALLOC(hexc_i%kdense2div,(double_grid%nbz_dense))
476 16 : ABI_MALLOC(hexc_i%div2kdense,(double_grid%nbz_coarse,double_grid%ndiv))
477 :
478 4 : call compute_corresp(double_grid,hexc_i%div2kdense,hexc_i%kdense2div)
479 :
480 4 : if (any(BSp%interp_mode == [2,3])) then
481 : ! Read a, b, c coefficient matrices from file.
482 : ! For the time being, we read the full matrix in a temporary array, and
483 : ! then we store the data in a form suitable for the interpolation.
484 2 : is_resonant=.TRUE.; diago_is_real=(.not.BSp%have_complex_ene); use_mpio=.FALSE.
485 :
486 2 : if (hexc%bs_files%in_hreso /= BSE_NOFILE) then
487 2 : hreso_fname = hexc%bs_files%in_hreso
488 : else
489 0 : hreso_fname = hexc%bs_files%out_hreso
490 : end if
491 :
492 2 : tmpfname = hreso_fname; ii = LEN_TRIM(hreso_fname)
493 :
494 : ! TODO: Write new IO routines to read MPI-distributed data in a format suitable for the interpolation
495 2 : tmpfname(ii-2:ii+1) = 'ABSR'
496 8 : ABI_MALLOC_OR_DIE(hexc_i%all_acoeffs,(hsize,hsize), ierr)
497 : call exc_read_rcblock(tmpfname,Bsp,is_resonant,diago_is_real,nsppol,BSp%nreh,hsize,1,hsize,&
498 2 : hexc_i%all_acoeffs,use_mpio,hexc%comm)
499 :
500 2 : tmpfname(ii-2:ii+1) = 'BBSR'
501 6 : ABI_MALLOC_OR_DIE(hexc_i%all_bcoeffs,(hsize,hsize), ierr)
502 : call exc_read_rcblock(tmpfname,Bsp,is_resonant,diago_is_real,nsppol,BSp%nreh,hsize,1,hsize,&
503 2 : hexc_i%all_bcoeffs,use_mpio,hexc%comm)
504 :
505 2 : tmpfname(ii-2:ii+1) = 'CBSR'
506 6 : ABI_MALLOC_OR_DIE(hexc_i%all_ccoeffs,(hsize,hsize), ierr)
507 : call exc_read_rcblock(tmpfname,Bsp,is_resonant,diago_is_real,nsppol,BSp%nreh,hsize,1,hsize,&
508 2 : hexc_i%all_ccoeffs,use_mpio,hexc%comm)
509 : end if
510 :
511 : ! Compute overlaps & compute all hmat
512 16 : ABI_MALLOC_OR_DIE(hexc_i%all_hmat,(hsize,hsize), ierr)
513 :
514 37252 : hexc_i%all_hmat(:,:) = hexc%hreso(:,:)
515 :
516 388 : do itt=1,hsize
517 388 : hexc_i%all_hmat(itt,itt) = hexc_i%all_hmat(itt,itt) - hexc%diag_coarse(itt)
518 : end do
519 :
520 : ! I don't need the diag_coarse any more
521 4 : ABI_FREE(hexc%diag_coarse)
522 :
523 : ! Compute diagonal part of the dense Ham
524 12 : ABI_MALLOC(hexc_i%diag_dense,(hexc_i%hsize_dense))
525 3076 : do itt=1,BSp%nreh_interp(spin1) ! 1 is for spin 1
526 3076 : hexc_i%diag_dense(itt) = Bsp%Trans_interp(itt,spin1)%en
527 : end do
528 :
529 4 : end subroutine hexc_interp_init
530 : !!***
531 :
532 : !-------------------------------------------------------------------
533 :
534 : !!****f* m_hexc/hexc_build_hinterp
535 : !! NAME
536 : !! hexc_build_hinterp
537 : !!
538 : !! FUNCTION
539 : !! Pre-compute interpolated hamiltonian and store it in memory
540 : !!
541 : !! SIDE EFFECTS
542 : !! hexc, hexc_i
543 : !! Pre-compute info to save CPU time when computing matmul
544 : !!
545 : !! SOURCE
546 :
547 2 : subroutine hexc_build_hinterp(hexc, hexc_i)
548 :
549 : !Arguments ---------------------------
550 : class(hexc_t),intent(inout) :: hexc
551 : class(hexc_interp_t),intent(inout) :: hexc_i
552 :
553 : !Local variables ---------------------
554 : integer :: ierr,ncerr,ncid
555 : character(len=500) :: msg
556 : !*****************************************************************************
557 :
558 2 : write(msg,"(a,f8.1,a)")"Memory needed for hinterp = ",one*(hexc_i%hsize_dense**2)*2*dp*b2Mb," Mb"
559 2 : call wrtout(std_out,msg,"COLL")
560 :
561 8 : ABI_MALLOC_OR_DIE(hexc_i%hinterp,(hexc_i%hsize_dense,hexc_i%hsize_dense), ierr)
562 :
563 : call hexc_compute_hinterp(hexc%BSp, hexc%hsize_coarse, hexc_i%hsize_dense, hexc_i%all_hmat, &
564 : hexc_i%interpolator%double_grid,hexc%nbnd_coarse, hexc_i%interpolator, &
565 : hexc_i%kdense2div, hexc_i%all_acoeffs,hexc_i%all_bcoeffs, hexc_i%all_ccoeffs, &
566 2 : hexc_i%Kmesh_dense, hexc_i%Vcp_dense, hexc%crystal%gmet, hexc_i%hinterp, hexc_i%m3_width)
567 :
568 2 : ABI_SFREE(hexc_i%all_acoeffs)
569 2 : ABI_SFREE(hexc_i%all_bcoeffs)
570 2 : ABI_SFREE(hexc_i%all_ccoeffs)
571 :
572 :
573 2 : if (hexc%BSp%prt_ncham) then
574 0 : ABI_COMMENT("Printing HEXC_I.nc file")
575 0 : ncerr = nctk_open_create(ncid, trim(hexc%BS_files%out_basename)//"_HEXC_I.nc", xmpi_comm_self)
576 0 : NCF_CHECK_MSG(ncerr, "Creating HEXC_I file")
577 : call exc_ham_ncwrite(ncid, hexc_i%Kmesh_dense, hexc%BSp, hexc_i%hsize_dense, hexc%BSp%nreh_interp, &
578 0 : hexc%BSp%vcks2t_interp, hexc_i%hinterp, hexc_i%diag_dense)
579 0 : NCF_CHECK(nf90_close(ncid))
580 : end if
581 :
582 2 : end subroutine hexc_build_hinterp
583 : !!***
584 :
585 : !-------------------------------------------------------------------
586 :
587 : !!****f* m_hexc/hexc_compute_subhinterp
588 : !! NAME
589 : !! hexc_compute_subhinterp
590 : !!
591 : !! FUNCTION
592 : !! Compute the interpolation for work_coeffs in dense mesh
593 : !!
594 : !! INPUTS
595 : !! BSp<excparam>=Parameters for BS run
596 : !! grid<double_grid_t>=Double grid info
597 : !! nbnd_coarse=Number of bands (= nbndv * nbndc)
598 : !! interpolator<interpolator_t>=Interpolation info
599 : !! kdense2div=Mapping between dense point and coarse point
600 : !! work_coeffs=Coefficients to be interpolated
601 : !! ikp_dense=Current kpoint
602 : !! overlaps=Wavefunction overlaps
603 : !!
604 : !! OUTPUT
605 : !! Cmat(nbnd_coarse) = Interpolated coefficients
606 : !!
607 : !! SOURCE
608 :
609 1397760 : subroutine hexc_compute_subhinterp(BSp,grid,nbnd_coarse,&
610 1397760 : interpolator,kdense2div,work_coeffs,Cmat,ikp_dense,overlaps)
611 :
612 : !Arguments ------------------------------------
613 : !scalars
614 : type(excparam),intent(in) :: BSp
615 : integer,intent(in) :: nbnd_coarse
616 : integer,intent(in) :: ikp_dense
617 : type(double_grid_t),intent(in) :: grid
618 : type(interpolator_t),target,intent(inout) :: interpolator
619 : !arrays
620 : integer,intent(in) :: kdense2div(grid%nbz_dense)
621 : complex(gwp),intent(in) :: overlaps(interpolator%mband_coarse,interpolator%mband_dense,interpolator%nvert)
622 : complex(dp),intent(in) :: work_coeffs(nbnd_coarse,interpolator%nvert)
623 : complex(dp),intent(out) :: Cmat(nbnd_coarse)
624 :
625 : !Local variables ------------------------------
626 : !scalars
627 : integer,parameter :: spin1=1,spin2=1
628 : integer :: iv1,ic1
629 : integer :: icp,ivp,idivp,ibndp_coarse,ibndp_coarse1,ineighbourp
630 : integer :: indwithnb
631 : integer :: lumo2,lomo2,humo2,homo2
632 : complex(dp) :: tmp_val, tmp2, tmp4
633 : !arrays
634 1397760 : complex(dp),contiguous, pointer :: btemp(:),ctemp(:)
635 : !*********************************************************************
636 :
637 1397760 : btemp => interpolator%btemp
638 1397760 : ctemp => interpolator%ctemp
639 :
640 135582720 : btemp = czero
641 135582720 : ctemp = czero
642 :
643 1397760 : lumo2 = BSp%lumo_spin(spin2)
644 1397760 : lomo2 = BSp%lomo_spin(spin2)
645 1397760 : humo2 = BSp%humo_spin(spin2)
646 1397760 : homo2 = BSp%homo_spin(spin2)
647 :
648 18170880 : Cmat = czero
649 :
650 1397760 : idivp = kdense2div(ikp_dense)
651 :
652 12579840 : do ineighbourp = 1,interpolator%nvert
653 :
654 : btemp(((ineighbourp-1)*nbnd_coarse+1):(ineighbourp*nbnd_coarse)) = &
655 146764800 : interpolator%interp_factors(ineighbourp,idivp)*work_coeffs(:,ineighbourp)
656 :
657 : end do !ineighbourp
658 :
659 : ! Loop over the (c', v') part of the right transition
660 5591040 : do ivp = lomo2,homo2
661 22364160 : do icp = lumo2,humo2
662 :
663 16773120 : ibndp_coarse = (ivp-lomo2)*BSp%maxnbndc+(icp-lumo2+1)
664 : ! Now we now it_dense, and itp_dense
665 :
666 150958080 : do ineighbourp = 1,interpolator%nvert
667 553512960 : do iv1 = lomo2, homo2
668 : ! BSp%lumo_spin(spin2),BSp%humo_spin(spin2)
669 402554880 : tmp4 = overlaps(iv1,ivp,ineighbourp)
670 :
671 2146959360 : do ic1 = lumo2,humo2
672 1610219520 : tmp2 = GWPC_CONJG(overlaps(ic1,icp,ineighbourp))
673 : ! BSp%lomo_spin(spin2),BSp%homo_spin(spin2)
674 :
675 1610219520 : ibndp_coarse1 = (iv1-lomo2)*BSp%maxnbndc+(ic1-lumo2+1)
676 1610219520 : indwithnb = (ineighbourp-1)*nbnd_coarse+ibndp_coarse1
677 :
678 2012774400 : ctemp(indwithnb) = tmp4 * tmp2
679 : end do ! iv1
680 : end do ! ic1
681 : end do ! ineighbourp
682 :
683 16773120 : tmp_val = xdotc(interpolator%nvert*nbnd_coarse,ctemp,1,btemp,1)
684 : !tmp_val = DOT_PRODUCT(ctemp,btemp)
685 :
686 20966400 : Cmat(ibndp_coarse) = tmp_val
687 : end do !ivp
688 : end do !icp
689 :
690 1397760 : nullify(btemp)
691 1397760 : nullify(ctemp)
692 :
693 1397760 : end subroutine hexc_compute_subhinterp
694 : !!***
695 :
696 : !----------------------------------------------------------------------
697 :
698 : !!****f* m_hexc/hexc_compute_hinterp
699 : !! NAME
700 : !! hexc_compute_hinterp
701 : !!
702 : !! FUNCTION
703 : !! Compute interpolated matrix elements for methods 2 and 3
704 : !!
705 : !! INPUTS
706 : !! BSp<type(excparam)=The parameter for the Bethe-Salpeter run.
707 : !! hsize_coarse=Size of the coarse Hamiltonian
708 : !! hsize_dense=Size of the dense Hamiltonian
709 : !! hmat(hsize_coarse,hsize_coarse,8)=Excitonic matrix
710 : !! grid<double_grid_t> = Correspondence between coarse and dense k-mesh.
711 : !! nbnd_coarse = Total number of bands
712 : !! interpolator<interpolator_t> = Interpolator
713 : !! kdense2div = Mapping kdense2div
714 : !! acoeffs, bcoeffs, ccoeffs = decomposition "W = a/q^2 + b/q + c"
715 : !! Kmesh_dense<type(kmesh_t)>=The list of k-points in the BZ, IBZ and symmetry tables.
716 : !! Vcp_dense<vcoul_t>=Coulomb interation in G-space on the dense Q-mesh
717 : !! gmet(3,3)=Metric tensor in G-space
718 : !!
719 : !! OUTPUT
720 : !! hinterp = Interpolated hamiltonian
721 : !!
722 : !! SOURCE
723 :
724 2 : subroutine hexc_compute_hinterp(BSp,hsize_coarse,hsize_dense,hmat,grid,nbnd_coarse,&
725 2 : & interpolator,kdense2div,acoeffs,bcoeffs,ccoeffs,Kmesh_dense,Vcp_dense,gmet,hinterp,&
726 : & m3_width)
727 :
728 : !Arguments ------------------------------------
729 : !scalars
730 : integer,intent(in) :: hsize_coarse,hsize_dense,nbnd_coarse !,ntrans
731 : real(dp),intent(in) :: m3_width
732 : type(excparam),intent(in) :: BSp
733 : type(double_grid_t),intent(in) :: grid
734 : type(vcoul_t),intent(in) :: Vcp_dense
735 : type(kmesh_t),intent(in) :: Kmesh_dense
736 : type(interpolator_t),target,intent(inout) :: interpolator
737 : !arrays
738 : integer,intent(in) :: kdense2div(grid%nbz_dense)
739 : real(dp),intent(in) :: gmet(3,3)
740 : complex(dp),intent(in) :: hmat(hsize_coarse,hsize_coarse)
741 : complex(dp),intent(in) :: acoeffs(hsize_coarse,hsize_coarse)
742 : complex(dp),intent(in) :: bcoeffs(hsize_coarse,hsize_coarse)
743 : complex(dp),intent(in) :: ccoeffs(hsize_coarse,hsize_coarse)
744 : complex(dp),intent(out) :: hinterp(hsize_dense,hsize_dense)
745 :
746 : !Local variables ------------------------------
747 : !scalars
748 : integer,parameter :: spin1=1, spin2=1
749 : integer :: ic,iv,iv1,ic1,ik_dense,ik_coarse,it_coarse,it_dense,idiv,ibnd_coarse,ibnd_coarse1,ineighbour
750 : integer :: icp,ivp,ikp_dense,ikp_coarse,itp_coarse,itp_dense,idivp,ibndp_coarse,ibndp_coarse1,ineighbourp,itp_coarse1
751 : integer :: itc,it_dense1,indwithnb, corresp_ind, limitnbz, inb,ierr
752 : real(dp) :: factor,vc_sqrt_qbz,qnorm
753 : complex(dp) :: term
754 : logical :: newway, use_herm
755 : !arrays
756 : real(dp) :: kmkp(3),q2(3),shift(3),qinred(3),tsec(2)
757 2 : complex(dp),allocatable :: Cmat(:,:,:) !Temp matrices for optimized version
758 2 : complex(dp),allocatable :: tmp_Cmat(:)
759 2 : complex(dp),allocatable :: work_coeffs(:,:)
760 2 : integer,allocatable :: band2it(:)
761 2 : complex(dp),contiguous, pointer :: btemp(:),ctemp(:)
762 : !************************************************************************
763 :
764 2 : call timab(696,1,tsec)
765 :
766 2 : newway = .True.
767 2 : use_herm = .True.
768 :
769 2 : if (any(BSp%interp_mode == [2,3])) then
770 2 : if (Vcp_dense%mode /= 'CRYSTAL' .and. Vcp_dense%mode /= 'AUXILIARY_FUNCTION') then
771 0 : ABI_ERROR('Vcp_dense%mode not implemented yet !')
772 : end if
773 : end if
774 :
775 2 : if (BSp%nsppol > 1) then
776 0 : ABI_ERROR("nsppol > 1 not yet implemented")
777 : end if
778 :
779 2 : factor = one/grid%ndiv
780 :
781 1181186 : hinterp = czero; term = czero
782 :
783 10 : ABI_MALLOC_OR_DIE(Cmat,(nbnd_coarse,nbnd_coarse,interpolator%nvert), ierr)
784 2514 : Cmat = czero
785 :
786 6 : ABI_MALLOC(band2it,(nbnd_coarse))
787 2 : call interpolator%alloc_work(nbnd_coarse*interpolator%nvert)
788 :
789 2 : btemp => interpolator%btemp
790 2 : ctemp => interpolator%ctemp
791 :
792 : if(newway) then
793 6 : ABI_MALLOC(tmp_Cmat,(nbnd_coarse))
794 8 : ABI_MALLOC(work_coeffs,(nbnd_coarse,interpolator%nvert))
795 : end if
796 :
797 130 : do ik_dense = 1,grid%nbz_dense
798 128 : write(std_out,*) "Kdense = ",ik_dense,"/",grid%nbz_dense
799 128 : ik_coarse = grid%dense_to_coarse(ik_dense)
800 : if(use_herm) then
801 : limitnbz = ik_dense
802 : else
803 : limitnbz = grid%nbz_dense
804 : end if
805 :
806 4290 : do ikp_dense = 1,limitnbz
807 4160 : ikp_coarse = grid%dense_to_coarse(ikp_dense)
808 :
809 16640 : do iv1 = BSp%lomo_spin(spin2),BSp%homo_spin(spin2)
810 66560 : do ic1 = BSp%lumo_spin(spin2),BSp%humo_spin(spin2)
811 49920 : itp_coarse1 = BSp%vcks2t(iv1,ic1,ikp_coarse,spin2)
812 49920 : ibndp_coarse1 = (iv1-BSp%lomo_spin(spin2))*BSp%maxnbndc+(ic1-BSp%lumo_spin(spin2)+1)
813 :
814 62400 : band2it(ibndp_coarse1) = itp_coarse1
815 : end do
816 : end do
817 :
818 :
819 4160 : if (any(BSp%interp_mode == [2,3])) then
820 : ! Check if we are along the diagonal
821 16640 : kmkp = Kmesh_dense%bz(:,ik_dense) - Kmesh_dense%bz(:,ikp_dense)
822 :
823 16640 : call wrap2_pmhalf(kmkp(:),q2(:),shift(:))
824 104000 : qinred = MATMUL(grid%kptrlatt_coarse,q2)
825 :
826 : ! We are outside the diagonal
827 16640 : if (BSp%interp_mode==3 .and. ANY((ABS(qinred)-tol7) > m3_width)) cycle
828 :
829 66560 : qnorm = two_pi*SQRT(DOT_PRODUCT(q2,MATMUL(gmet,q2)))
830 :
831 5696 : if(ALL(ABS(q2(:)) < 1.e-3)) then
832 128 : vc_sqrt_qbz = SQRT(Vcp_dense%i_sz)
833 : else
834 4032 : vc_sqrt_qbz = SQRT(four_pi/qnorm**2)
835 : end if
836 :
837 : !!DEBUG CHK !
838 : !!COMPUTE Qpoint
839 : !call findqg0(iq_bz,g0,kmkp,Qmesh_dense%nbz,Qmesh_dense%bz,BSp%mG0)
840 :
841 : !! * Get iq_ibz, and symmetries from iq_bz
842 : !call qmesh_dense%get_BZ_item(Qmesh_dense,iq_bz,qbz,iq_ibz,isym_q,itim_q)
843 :
844 : !if(iq_ibz > 1 .and. ABS(vc_sqrt_qbz - Vcp_dense%vc_sqrt(1,iq_ibz)) > 1.e-3) then
845 : ! write(*,*) "vc_sqrt_qbz = ",vc_sqrt_qbz
846 : ! write(*,*) "Vcp_dense%vc_sqrt(1,iq_ibz) = ",Vcp_dense%vc_sqrt(1,iq_ibz)
847 : ! ABI_ERROR("vcp are not the same !")
848 : !else if(iq_ibz == 1 .and. ABS(vc_sqrt_qbz - SQRT(Vcp_dense%i_sz)) > 1.e-3) then
849 : ! write(*,*) "vc_sqrt_qbz = ",vc_sqrt_qbz
850 : ! write(*,*) "SQRT(Vcp_dense%i_sz) = ",SQRT(Vcp_dense%i_sz)
851 : ! ABI_ERROR("vcp are not the same !")
852 : !end if
853 : !!END DEBUG CHK !
854 : end if
855 :
856 : if(newway) then
857 :
858 5229120 : Cmat = czero
859 :
860 436800 : work_coeffs = czero
861 :
862 37440 : do ineighbour = 1,interpolator%nvert
863 :
864 : ! Loop over the (c, v) part of the left transition
865 137280 : do iv = BSp%lomo_spin(spin1),BSp%homo_spin(spin1)
866 532480 : do ic = BSp%lumo_spin(spin1),BSp%humo_spin(spin1)
867 :
868 399360 : it_dense = BSp%vcks2t_interp(iv,ic,ik_dense,spin1)
869 399360 : it_coarse = BSp%vcks2t(iv,ic,ik_coarse,spin1)
870 399360 : ibnd_coarse = (iv-BSp%lomo_spin(spin1))*BSp%maxnbndc+(ic-BSp%lumo_spin(spin1)+1)
871 :
872 399360 : itc = interpolator%corresp(it_coarse,ineighbour,spin1)
873 :
874 399360 : if (any(BSp%interp_mode == [1,3,4])) then
875 :
876 : !work_coeffs(:,:) = hmat(itc,band2it(:),:)
877 1797120 : do inb = 1,interpolator%nvert
878 40135680 : work_coeffs(:,inb) = hmat(itc,interpolator%corresp(band2it(:),inb,spin2))
879 : end do
880 :
881 : call hexc_compute_subhinterp(BSp,grid,nbnd_coarse,&
882 : & interpolator,kdense2div,&
883 : & work_coeffs,tmp_Cmat,ikp_dense,&
884 199680 : & interpolator%overlaps(:,:,:,ikp_dense,spin2))
885 :
886 199680 : if(any(BSp%interp_mode == [1,4])) then
887 0 : Cmat(ibnd_coarse,:,ineighbour) = tmp_Cmat
888 199680 : else if (BSp%interp_mode == 3) then
889 2595840 : Cmat(ibnd_coarse,:,ineighbour) = -tmp_Cmat
890 : end if
891 : end if
892 :
893 :
894 399360 : if (any(BSp%interp_mode == [2,3])) then
895 : !work_coeffs(:,:) = acoeffs(itc,band2it(:),:)
896 3594240 : do inb = 1,interpolator%nvert
897 80271360 : work_coeffs(:,inb) = acoeffs(itc,interpolator%corresp(band2it(:),inb,spin2))
898 : end do
899 :
900 : call hexc_compute_subhinterp(BSp,grid,nbnd_coarse,&
901 : & interpolator,kdense2div,&
902 : & work_coeffs,tmp_Cmat,ikp_dense,&
903 399360 : & interpolator%overlaps(:,:,:,ikp_dense,spin2))
904 :
905 5191680 : tmp_Cmat = tmp_Cmat * (vc_sqrt_qbz**2)
906 5191680 : Cmat(ibnd_coarse,:,ineighbour) = Cmat(ibnd_coarse,:,ineighbour) + tmp_Cmat
907 : end if
908 :
909 :
910 399360 : if (any(BSp%interp_mode == [2,3])) then
911 : !work_coeffs(:,:) = bcoeffs(itc,band2it(:),:)
912 3594240 : do inb = 1,interpolator%nvert
913 80271360 : work_coeffs(:,inb) = bcoeffs(itc,interpolator%corresp(band2it(:),inb,spin2))
914 : end do
915 :
916 : call hexc_compute_subhinterp(BSp,grid,nbnd_coarse,&
917 : & interpolator,kdense2div,&
918 : & work_coeffs,tmp_Cmat,ikp_dense,&
919 399360 : & interpolator%overlaps(:,:,:,ikp_dense,spin2))
920 :
921 5191680 : tmp_Cmat = tmp_Cmat * (vc_sqrt_qbz)
922 5191680 : Cmat(ibnd_coarse,:,ineighbour) = Cmat(ibnd_coarse,:,ineighbour) + tmp_Cmat
923 : end if
924 :
925 :
926 499200 : if (any(BSp%interp_mode == [2,3])) then
927 : !work_coeffs(:,:) = ccoeffs(itc,band2it(:),:)
928 3594240 : do inb = 1,interpolator%nvert
929 80271360 : work_coeffs(:,inb) = ccoeffs(itc,interpolator%corresp(band2it(:),inb,spin2))
930 : end do
931 :
932 : call hexc_compute_subhinterp(BSp,grid,nbnd_coarse,&
933 : & interpolator,kdense2div,&
934 : & work_coeffs,tmp_Cmat,ikp_dense,&
935 399360 : & interpolator%overlaps(:,:,:,ikp_dense,spin2))
936 :
937 5191680 : Cmat(ibnd_coarse,:,ineighbour) = Cmat(ibnd_coarse,:,ineighbour) + tmp_Cmat
938 : end if
939 : end do ! ic
940 : end do ! iv
941 : end do ! ineighbour
942 :
943 : else
944 : ! Loop over the (c, v) part of the left transition
945 : do iv = BSp%lomo_spin(spin1),BSp%homo_spin(spin1)
946 : do ic = BSp%lumo_spin(spin1),BSp%humo_spin(spin1)
947 :
948 : it_dense = BSp%vcks2t_interp(iv,ic,ik_dense,spin1)
949 : it_coarse = BSp%vcks2t(iv,ic,ik_coarse,spin1)
950 : ibnd_coarse = (iv-BSp%lomo_spin(spin1))*BSp%maxnbndc+(ic-BSp%lumo_spin(spin1)+1)
951 :
952 : ! Loop over the (c', v') part of the right transition
953 : do ivp = BSp%lomo_spin(spin2),BSp%homo_spin(spin2)
954 : do icp = BSp%lumo_spin(spin2),BSp%humo_spin(spin2)
955 :
956 : itp_dense = BSp%vcks2t_interp(ivp,icp,ikp_dense,spin2)
957 : itp_coarse = BSp%vcks2t(ivp,icp,ikp_coarse,spin2)
958 : ibndp_coarse = (ivp-Bsp%lomo_spin(spin2))*BSp%maxnbndc+(icp-BSp%lumo_spin(spin2)+1)
959 : ! Now we now it_dense, and itp_dense
960 :
961 : idivp = kdense2div(ikp_dense)
962 :
963 : btemp = czero; ctemp = czero
964 :
965 : ! MG TODO: This way of looping is not optimal
966 : do ineighbour = 1,interpolator%nvert
967 : itc = interpolator%corresp(it_coarse,ineighbour,spin1)
968 :
969 : do ineighbourp = 1,interpolator%nvert
970 :
971 : do iv1 = BSp%lomo_spin(spin2),BSp%homo_spin(spin2)
972 : do ic1 = BSp%lumo_spin(spin2),BSp%humo_spin(spin2)
973 :
974 : ibndp_coarse1 = (iv1-BSp%lomo_spin(spin2))*BSp%maxnbndc+(ic1-BSp%lumo_spin(spin2)+1)
975 : indwithnb = (ineighbourp-1)*nbnd_coarse+ibndp_coarse1
976 : itp_coarse1 = BSp%vcks2t(iv1,ic1,ikp_coarse,spin2)
977 : corresp_ind = interpolator%corresp(itp_coarse1,ineighbourp,spin2)
978 :
979 : select case (BSp%interp_mode)
980 : case (1,4)
981 : interpolator%btemp(indwithnb) = hmat(itc,corresp_ind)
982 : case (2)
983 : interpolator%btemp(indwithnb) = acoeffs(itc,corresp_ind)*(vc_sqrt_qbz**2) &
984 : & + bcoeffs(itc,corresp_ind)*(vc_sqrt_qbz) &
985 : & + ccoeffs(itc,corresp_ind)
986 : case (3)
987 : ! Diff between divergence and hmat
988 : interpolator%btemp(indwithnb) = acoeffs(itc,corresp_ind)*(vc_sqrt_qbz**2) &
989 : & + bcoeffs(itc,corresp_ind)*(vc_sqrt_qbz) &
990 : & + ccoeffs(itc,corresp_ind) &
991 : & - hmat(itc,corresp_ind)
992 : case default
993 : ABI_ERROR("Wrong Bsp%interp_mode")
994 : end select
995 :
996 : ctemp(indwithnb) = &
997 : & interpolator%overlaps(iv1,ivp,ineighbourp,ikp_dense,spin2) &
998 : & * GWPC_CONJG(interpolator%overlaps(ic1,icp,ineighbourp,ikp_dense,spin2)) &
999 : & *interpolator%interp_factors(ineighbourp,idivp)
1000 : end do ! ic1
1001 : end do !iv1
1002 :
1003 : end do !ineighbourp
1004 : Cmat(ibnd_coarse,ibndp_coarse,ineighbour) = xdotc(interpolator%nvert*nbnd_coarse,&
1005 : & ctemp,1,btemp,1)
1006 : end do !ineighbour
1007 :
1008 : end do !icp
1009 : end do !ivp
1010 :
1011 : end do !ic
1012 : end do !iv
1013 :
1014 : end if
1015 :
1016 16768 : do iv = BSp%lomo_spin(spin1),BSp%homo_spin(spin1)
1017 66560 : do ic = BSp%lumo_spin(spin1),BSp%humo_spin(spin1)
1018 49920 : it_dense = BSp%vcks2t_interp(iv,ic,ik_dense,spin1)
1019 49920 : it_coarse = BSp%vcks2t(iv,ic,ik_coarse,spin1)
1020 49920 : ibnd_coarse = (iv-BSp%lomo_spin(spin1))*BSp%maxnbndc+(ic-BSp%lumo_spin(spin1)+1)
1021 :
1022 49920 : idiv = kdense2div(ik_dense)
1023 :
1024 212160 : do ivp = BSp%lomo_spin(spin2),BSp%homo_spin(spin2)
1025 798720 : do icp = BSp%lumo_spin(spin2),BSp%humo_spin(spin2)
1026 599040 : itp_dense = BSp%vcks2t_interp(ivp,icp,ikp_dense,spin2)
1027 599040 : itp_coarse = BSp%vcks2t(ivp,icp,ikp_coarse,spin2)
1028 599040 : ibndp_coarse = (ivp-Bsp%lomo_spin(spin2))*BSp%maxnbndc+(icp-BSp%lumo_spin(spin2)+1)
1029 :
1030 : ! Hermicity
1031 : if(use_herm .and. it_dense < itp_dense) then
1032 : continue
1033 : end if
1034 :
1035 : !btemp = czero; ctemp = czero
1036 :
1037 5391360 : do ineighbour = 1,interpolator%nvert
1038 19768320 : do iv1 = BSp%lomo_spin(spin1),BSp%homo_spin(spin1)
1039 76677120 : do ic1 = BSp%lumo_spin(spin1),BSp%humo_spin(spin1)
1040 57507840 : ibnd_coarse1 = (iv1-BSp%lomo_spin(spin1))*BSp%maxnbndc+(ic1-BSp%lumo_spin(spin1)+1)
1041 57507840 : it_dense1 = BSp%vcks2t_interp(iv1,ic1,ik_dense,spin1)
1042 57507840 : indwithnb = (ineighbour-1)*nbnd_coarse+ibnd_coarse1
1043 :
1044 57507840 : btemp(indwithnb) = Cmat(ibnd_coarse1,ibndp_coarse,ineighbour)
1045 :
1046 : ctemp(indwithnb) = GWPC_CONJG(interpolator%overlaps(iv1,iv,ineighbour,ik_dense,spin1)) &
1047 : & *interpolator%overlaps(ic1,ic,ineighbour,ik_dense,spin1) &
1048 71884800 : & *interpolator%interp_factors(ineighbour,idiv)
1049 : end do !ic1
1050 : end do !iv1
1051 : end do !ineighbour
1052 :
1053 : ! Save interpolated value.
1054 748800 : hinterp(it_dense,itp_dense) = xdotc(interpolator%nvert*nbnd_coarse,ctemp,1,btemp,1)
1055 : !DOT_PRODUCT(ctemp,btemp)
1056 :
1057 : end do !icp
1058 : end do !ivp
1059 :
1060 : end do !ic
1061 : end do !iv
1062 :
1063 : end do !ikp
1064 : end do !ik
1065 :
1066 : ! Enforce hermiticity
1067 : if(use_herm) then
1068 1538 : do itp_dense = 1,BSp%nreh_interp(spin2)
1069 592130 : do it_dense = itp_dense,BSp%nreh_interp(spin1)
1070 592128 : if(it_dense == itp_dense) then
1071 1536 : hinterp(itp_dense,it_dense) = DBLE(hinterp(it_dense,itp_dense))
1072 : else
1073 589056 : hinterp(itp_dense,it_dense) = CONJG(hinterp(it_dense,itp_dense))
1074 : end if
1075 : end do
1076 : end do
1077 : end if
1078 :
1079 2 : ABI_FREE(Cmat)
1080 2 : ABI_FREE(band2it)
1081 :
1082 2 : ABI_SFREE(tmp_Cmat)
1083 2 : ABI_SFREE(work_coeffs)
1084 :
1085 2 : nullify(ctemp)
1086 2 : nullify(btemp)
1087 2 : call interpolator%int_free()
1088 :
1089 1181186 : hinterp = hinterp*factor
1090 :
1091 2 : call timab(696,2,tsec)
1092 :
1093 2 : end subroutine hexc_compute_hinterp
1094 : !!***
1095 :
1096 : !----------------------------------------------------------------------
1097 :
1098 : !!****f* m_hexc/hexc_free
1099 : !! NAME
1100 : !! hexc_free
1101 : !!
1102 : !! FUNCTION
1103 : !! Destroy the interpolator object in memory
1104 : !!
1105 : !! SOURCE
1106 :
1107 23 : subroutine hexc_free(hexc)
1108 :
1109 : !Arguments ---------------------------
1110 : class(hexc_t),intent(inout) :: hexc
1111 : !*****************************************************************************
1112 :
1113 23 : if (associated(hexc%bsp)) then
1114 23 : nullify(hexc%bsp)
1115 : end if
1116 :
1117 23 : if (associated(hexc%crystal)) then
1118 23 : nullify(hexc%crystal)
1119 : end if
1120 :
1121 23 : if (associated(hexc%kmesh_coarse)) then
1122 23 : nullify(hexc%kmesh_coarse)
1123 : end if
1124 :
1125 23 : ABI_SFREE(hexc%hreso)
1126 23 : ABI_SFREE(hexc%hcoup)
1127 23 : ABI_SFREE(hexc%diag_coarse)
1128 :
1129 23 : end subroutine hexc_free
1130 : !!***
1131 :
1132 : !-------------------------------------------------------------------
1133 :
1134 : !!****f* m_hexc/hexc_interp_free
1135 : !! NAME
1136 : !! hexc_interp_free
1137 : !!
1138 : !! FUNCTION
1139 : !! Free dynamic memory.
1140 : !!
1141 : !! SOURCE
1142 :
1143 23 : subroutine hexc_interp_free(hexc_i)
1144 :
1145 : !Arguments ---------------------------
1146 : class(hexc_interp_t),intent(inout) :: hexc_i
1147 : !*****************************************************************************
1148 :
1149 23 : ABI_SFREE(hexc_i%kdense2div)
1150 23 : ABI_SFREE(hexc_i%div2kdense)
1151 23 : ABI_SFREE(hexc_i%diag_dense)
1152 23 : ABI_SFREE(hexc_i%hinterp)
1153 23 : ABI_SFREE(hexc_i%all_hmat)
1154 23 : ABI_SFREE(hexc_i%all_acoeffs)
1155 23 : ABI_SFREE(hexc_i%all_bcoeffs)
1156 23 : ABI_SFREE(hexc_i%all_ccoeffs)
1157 :
1158 23 : if (associated(hexc_i%kmesh_dense)) then
1159 4 : nullify(hexc_i%kmesh_dense)
1160 : end if
1161 :
1162 23 : if (associated(hexc_i%vcp_dense)) then
1163 4 : nullify(hexc_i%vcp_dense)
1164 : end if
1165 :
1166 23 : call hexc_i%interpolator%free()
1167 :
1168 23 : end subroutine hexc_interp_free
1169 : !!***
1170 :
1171 : !----------------------------------------------------------------------
1172 :
1173 : !!****f* m_hexc/hexc_interp_matmul
1174 : !! NAME
1175 : !! hexc_interp_matmul
1176 : !!
1177 : !! FUNCTION
1178 : !! Compute matrix-vector product Hmat * phi by interpolating coarse Hmat
1179 : !!
1180 : !! INPUTS
1181 : !! BSp<type(excparam)>=Parameters defining the BS calculation
1182 : !! hsize_coarse = Size of the coarse hamiltonian
1183 : !! hsize_dense = Size of the dense hamiltonian
1184 : !! hmat(hsize_coarse,hsize_coarse,8) = coarse hamiltonian
1185 : !! phi(hsize_dense) = ket on which apply the matrix
1186 : !! grid <double_grid_t> = Correspondence between coarse and dense k-mesh.
1187 : !! nbnd_coarse = Total number of bands
1188 : !! interpolator<interpolator_t> = Interpolator
1189 : !! div2kdense = Mapping from coarse and division -> dense mesh
1190 : !! kdense2div = Mapping from dense mesh -> division
1191 : !!
1192 : !! OUTPUT
1193 : !! hphi(hsize_dense) = Interp(hmat)*phi
1194 : !!
1195 : !! SOURCE
1196 :
1197 3259 : subroutine hexc_interp_matmul(BSp,hsize_coarse,hsize_dense,hmat,phi,hphi,grid,&
1198 3259 : nbnd_coarse,interpolator,div2kdense,kdense2div)
1199 :
1200 : !Arguments ------------------------------------
1201 : !scalars
1202 : integer,intent(in) :: hsize_coarse,hsize_dense,nbnd_coarse !,ntrans
1203 : type(excparam),intent(in) :: BSp
1204 : type(double_grid_t),intent(in) :: grid
1205 : type(interpolator_t),intent(in) :: interpolator
1206 : !arrays
1207 : integer,intent(in) :: div2kdense(grid%nbz_coarse,grid%ndiv), kdense2div(grid%nbz_dense)
1208 : complex(dp),intent(in) :: phi(hsize_dense)
1209 : complex(dp),intent(in) :: hmat(hsize_coarse,hsize_coarse)
1210 : complex(dp),intent(inout) :: hphi(hsize_dense)
1211 :
1212 : !Local variables ------------------------------
1213 : !scalars
1214 : integer :: itt,ik_dense,ik_coarse,it_coarse, ic,iv,iv1,ic1, ibnd_coarse
1215 : integer :: ibnd_coarse1, ineighbour,idense,ikpt
1216 : integer :: my_k1,my_k2,ind_with_nb,is, is1
1217 : real(dp) :: factor
1218 : complex(dp) :: tmp
1219 : logical,parameter :: use_blas=.True.
1220 : !arrays
1221 6518 : integer :: allindices(nbnd_coarse)
1222 6518 : complex(dp) :: allp(hsize_coarse,interpolator%nvert), test(hsize_coarse)
1223 6518 : complex(dp) :: ophi(grid%nbz_dense,interpolator%nvert,nbnd_coarse)
1224 3259 : complex(dp),allocatable :: b(:), c(:),A(:,:), tmp_array(:), tmp_array2(:,:)
1225 : !************************************************************************
1226 :
1227 3259 : factor = one/grid%ndiv
1228 : !hphi = czero
1229 :
1230 : ! Outer index : k point in the dense zone
1231 : ! Sum over vc
1232 : ! Index of result : k point in the dense zone, v2,c2,neighbour
1233 :
1234 : ! Parallelization on nbz in the coarse mesh !
1235 3259 : my_k1 = 1
1236 3259 : my_k2 = grid%nbz_coarse
1237 :
1238 13036 : ABI_MALLOC(A,(interpolator%nvert*nbnd_coarse,nbnd_coarse))
1239 9777 : ABI_MALLOC(b,(nbnd_coarse))
1240 9777 : ABI_MALLOC(c,(interpolator%nvert*nbnd_coarse))
1241 :
1242 14673866 : c = czero; ophi = czero
1243 :
1244 211835 : do ik_dense = 1,grid%nbz_dense
1245 : ! if( ik_dense is not in my set of k-points)
1246 : ! ! continue
1247 : !
1248 417152 : do is1 = 1, BSp%nsppol
1249 1042880 : do iv1 = BSp%lomo_spin(is1),Bsp%homo_spin(is1)
1250 3337216 : do ic1 = BSp%lumo_spin(is1),Bsp%humo_spin(is1)
1251 2502912 : ibnd_coarse = (iv1-BSp%lomo_spin(is1))*BSp%maxnbndc+(ic1-BSp%lumo_spin(is1)+1)
1252 2502912 : itt = BSp%vcks2t_interp(iv1,ic1,ik_dense,is1)
1253 3128640 : allindices(ibnd_coarse) = itt
1254 : end do !ic1
1255 : end do !iv1
1256 : end do !is1
1257 :
1258 2711488 : b(:) = phi(allindices(:))
1259 :
1260 417152 : do is = 1, BSp%nsppol
1261 1042880 : do iv = BSp%lomo_spin(is),Bsp%homo_spin(is)
1262 3337216 : do ic = BSp%lumo_spin(is),Bsp%humo_spin(is)
1263 2502912 : ibnd_coarse = (iv-BSp%lomo_spin(is))*BSp%maxnbndc+(ic-BSp%lumo_spin(is)+1)
1264 2502912 : idense = Bsp%vcks2t_interp(iv,ic,ik_dense,is)
1265 :
1266 17313600 : do ineighbour = 1,interpolator%nvert
1267 14184960 : ind_with_nb = (ineighbour-1)*(nbnd_coarse)+ibnd_coarse
1268 :
1269 : !A(ind_with_nb,:) = overlaps(allindices(:),ibnd_coarse,ineighbour)
1270 :
1271 : ! Should be optimized !!!
1272 59242752 : do iv1 = BSp%lomo_spin(is),Bsp%homo_spin(is)
1273 226959360 : do ic1 = BSp%lumo_spin(is),Bsp%humo_spin(is)
1274 170219520 : ibnd_coarse1 = (iv1-BSp%lomo_spin(is))*BSp%maxnbndc+(ic1-BSp%lumo_spin(is)+1)
1275 : A(ind_with_nb,ibnd_coarse1) = GWPC_CONJG(interpolator%overlaps(iv,iv1,ineighbour,ik_dense,is)) &
1276 212774400 : & *interpolator%overlaps(ic,ic1,ineighbour,ik_dense,is)
1277 : end do !ic1
1278 : end do !iv1
1279 : end do !ineighbour
1280 : end do !ic
1281 : end do !iv
1282 : end do !is
1283 :
1284 : if(use_blas) then
1285 208576 : call xgemv('N',interpolator%nvert*nbnd_coarse,nbnd_coarse,cone,A,interpolator%nvert*nbnd_coarse,b,1,czero,c,1)
1286 : else
1287 : c = MATMUL(A,b)
1288 : end if
1289 :
1290 420411 : do is = 1, BSp%nsppol
1291 1042880 : do iv = BSp%lomo_spin(is),BSp%homo_spin(is)
1292 3337216 : do ic = BSp%lumo_spin(is),BSp%humo_spin(is)
1293 2502912 : ibnd_coarse = (iv-BSp%lomo_spin(is))*BSp%maxnbndc+(ic-BSp%lumo_spin(is)+1)
1294 17313600 : do ineighbour = 1,interpolator%nvert
1295 14184960 : ind_with_nb = (ineighbour-1)*(nbnd_coarse)+ibnd_coarse
1296 16687872 : ophi(ik_dense,ineighbour,ibnd_coarse) = c(ind_with_nb)
1297 : end do !ineighbour
1298 : end do !ic
1299 : end do !iv
1300 : end do !is
1301 :
1302 : end do !ik_dense
1303 :
1304 3259 : ABI_FREE(A)
1305 3259 : ABI_FREE(b)
1306 3259 : ABI_FREE(c)
1307 :
1308 : !call xmpi_sum_(ophi,comm,ierr)
1309 :
1310 : ! Outer index : k,v,c in the coarse zone, ineighbour
1311 : ! Sum over all k-dense relative to one coarse point
1312 : ! Index of result : k,v,c in the coarse zone, ineighbour
1313 :
1314 9777 : ABI_MALLOC(b,(grid%ndiv))
1315 6518 : ABI_MALLOC(c,(grid%ndiv))
1316 :
1317 1794849 : allp = czero
1318 :
1319 6518 : do is = 1, BSp%nsppol
1320 24988 : do ineighbour = 1,interpolator%nvert
1321 :
1322 1794849 : do it_coarse = 1, BSp%nreh(is)
1323 : ibnd_coarse = (Bsp%trans(it_coarse,is)%v-BSp%lomo_spin(is))*BSp%maxnbndc+&
1324 1773120 : (BSp%Trans(it_coarse,is)%c-BSp%lumo_spin(is)+1)
1325 1773120 : ik_coarse = BSp%trans(it_coarse,is)%k
1326 : !b(:) = interp_factors(it_coarse,ineighbour,:)
1327 15958080 : b(:) = interpolator%interp_factors(ineighbour,:)
1328 : !c(:) = ophi(indices(it_coarse,:),ineighbour,ibnd_coarse)
1329 15958080 : c(:) = ophi(div2kdense(ik_coarse,:),ineighbour,ibnd_coarse)
1330 15958080 : tmp = DOT_PRODUCT(b,c)
1331 1791590 : allp(it_coarse,ineighbour) = tmp
1332 : end do
1333 :
1334 : end do
1335 : end do
1336 :
1337 : !call xmpi_sum_(allp,comm,ierr)
1338 :
1339 3259 : ABI_FREE(b)
1340 3259 : ABI_FREE(c)
1341 :
1342 9777 : ABI_MALLOC(tmp_array,(hsize_coarse))
1343 13036 : ABI_MALLOC(tmp_array2,(hsize_coarse,hsize_coarse))
1344 316123 : tmp_array(:) = czero
1345 30351067 : tmp_array2(:,:) = czero
1346 :
1347 316123 : test = czero
1348 :
1349 : ! Second step : Multiplication by hmat
1350 21729 : do ineighbour = 1,interpolator%nvert
1351 173802700 : tmp_array2 = factor*hmat(:,interpolator%corresp(:,ineighbour,1)) ! 1 is for spin
1352 3259 : if(use_blas) then
1353 : !call xgemv('N',hsize_coarse,hsize_coarse,cone,factor*(hmat(:,:,ineighbour)),hsize_coarse,allp(:,ineighbour),1,czero,tmp_array,1)
1354 : !tmp_array2 = hmat(:,:,ineighbour)
1355 : !tmp_array2 = factor*tmp_array2
1356 18470 : call xgemv('N',hsize_coarse,hsize_coarse,cone,tmp_array2,hsize_coarse,allp(:,ineighbour),1,czero,tmp_array,1)
1357 1791590 : test = test + tmp_array
1358 : else
1359 : !test = test+MATMUL(factor*(hmat(:,interpolator%corresp(:,ineighbour,1))),allp(:,ineighbour))
1360 : test = test+MATMUL(tmp_array2,allp(:,ineighbour))
1361 : end if
1362 : end do
1363 :
1364 3259 : ABI_FREE(tmp_array)
1365 3259 : ABI_FREE(tmp_array2)
1366 :
1367 : ! Outer index : ineighbour
1368 : ! Sum over all v c
1369 : ! Index of result : ineighbour, k_dense, v,c
1370 13036 : ABI_MALLOC(A,(nbnd_coarse,nbnd_coarse))
1371 6518 : ABI_MALLOC(b,(nbnd_coarse))
1372 6518 : ABI_MALLOC(c,(nbnd_coarse))
1373 42367 : c = czero
1374 :
1375 21729 : do ineighbour = 1,interpolator%nvert
1376 1203809 : do ik_dense = 1,grid%nbz_dense
1377 :
1378 2364160 : do is1 = 1, Bsp%nsppol
1379 5910400 : do iv1 = Bsp%lomo_spin(is1),Bsp%homo_spin(is1)
1380 18913280 : do ic1 = BSp%lumo_spin(is1), Bsp%humo_spin(is1)
1381 14184960 : ibnd_coarse = (iv1-BSp%lomo_spin(is1))*BSp%maxnbndc+(ic1-BSp%lumo_spin(is1)+1)
1382 :
1383 14184960 : ik_coarse = grid%dense_to_coarse(ik_dense)
1384 14184960 : itt = BSp%vcks2t(iv1,ic1,ik_coarse,is1)
1385 17731200 : b(ibnd_coarse) = test(interpolator%corresp(itt,ineighbour,is1))
1386 : end do ! ic1
1387 : end do ! iv1
1388 : end do ! is1
1389 :
1390 2364160 : do is = 1, BSp%nsppol
1391 5910400 : do iv = BSp%lomo_spin(is),Bsp%homo_spin(is)
1392 18913280 : do ic = BSp%lumo_spin(is),BSp%humo_spin(is)
1393 14184960 : ibnd_coarse = (iv-BSp%lomo_spin(is))*Bsp%maxnbndc+(ic-BSp%lumo_spin(is)+1)
1394 14184960 : idense = BSp%vcks2t_interp(iv,ic,ik_dense,is)
1395 :
1396 : !A(ibnd_coarse,:) = CONJG(overlaps(idense,:,ineighbour))
1397 :
1398 : ! Should be optimized !!!
1399 60286080 : do iv1 = BSp%lomo_spin(is),Bsp%homo_spin(is)
1400 226959360 : do ic1 = BSp%lumo_spin(is),Bsp%humo_spin(is)
1401 170219520 : ibnd_coarse1 = (iv1-BSp%lomo_spin(is))*BSp%maxnbndc+(ic1-BSp%lumo_spin(is)+1)
1402 : A(ibnd_coarse,ibnd_coarse1) = (interpolator%overlaps(iv1,iv,ineighbour,ik_dense,is)) &
1403 212774400 : *GWPC_CONJG(interpolator%overlaps(ic1,ic,ineighbour,ik_dense,is))
1404 : end do !ic1
1405 : end do !iv1
1406 : end do ! ic
1407 : end do !iv
1408 : end do !is
1409 :
1410 : if(use_blas) then
1411 1182080 : call xgemv('N',nbnd_coarse,nbnd_coarse,cone,A,nbnd_coarse,b,1,czero,c,1)
1412 : else
1413 : c = MATMUL(A,b)
1414 : end if
1415 :
1416 2382630 : do is = 1, BSp%nsppol
1417 5910400 : do iv = BSp%lomo_spin(is),Bsp%homo_spin(is)
1418 18913280 : do ic = BSp%lumo_spin(is),BSp%humo_spin(is)
1419 14184960 : ibnd_coarse = (iv-BSp%lomo_spin(is))*BSp%maxnbndc+(ic-BSp%lumo_spin(is)+1)
1420 14184960 : idense = Bsp%vcks2t_interp(iv,ic,ik_dense,is)
1421 : !ophi(ik_dense,ineighbour,ibnd_coarse) = c(idense)
1422 17731200 : ophi(ik_dense,ineighbour,ibnd_coarse) = c(ibnd_coarse)
1423 : end do
1424 : end do
1425 : end do
1426 :
1427 : end do ! ik_dense
1428 : end do ! ineighbour
1429 :
1430 : !call xmpi_sum_(ophi,comm,ierr)
1431 :
1432 3259 : ABI_FREE(A)
1433 3259 : ABI_FREE(b)
1434 3259 : ABI_FREE(c)
1435 :
1436 : ! Outer indices : it_dense
1437 : ! Sum over neighbours
1438 : ! Index of result : it_dense (ik,ic,iv dense)
1439 :
1440 9777 : ABI_MALLOC(b,(interpolator%nvert))
1441 6518 : ABI_MALLOC(c,(interpolator%nvert))
1442 :
1443 6518 : do is = 1, BSp%nsppol
1444 2509430 : do itt = 1,BSp%nreh_interp(is)
1445 : ! From itt -> ik_ibz,ic,iv
1446 2502912 : ik_dense = BSp%Trans_interp(itt,is)%k
1447 2502912 : ic = BSp%Trans_interp(itt,is)%c
1448 2502912 : iv = BSp%Trans_interp(itt,is)%v
1449 :
1450 : ! From ik_ibz in the dense mesh -> indices_dense
1451 2502912 : ik_coarse = grid%dense_to_coarse(ik_dense)
1452 2502912 : it_coarse = BSp%vcks2t(iv,ic,ik_coarse,is)
1453 :
1454 2502912 : ibnd_coarse = (iv-BSp%lomo_spin(is))*BSp%maxnbndc+(ic-BSp%lumo_spin(is)+1)
1455 :
1456 2502912 : ikpt = kdense2div(ik_dense)
1457 : !ikpt = -1
1458 : !do ix = 1,grid%ndiv
1459 : ! if (indices(it_coarse,ix) == ik_dense) then
1460 : ! ikpt = ix
1461 : ! exit
1462 : ! end if
1463 : !end do
1464 : !ABI_CHECK(ikpt/=-1,"Cannot find ik_dense")
1465 :
1466 : !b = interp_factors(it_coarse,:,ikpt)
1467 19190784 : b = interpolator%interp_factors(:,ikpt)
1468 19190784 : c = ophi(ik_dense,:,ibnd_coarse)
1469 :
1470 2506171 : hphi(itt) = hphi(itt) + xdotc(interpolator%nvert, b, 1, c, 1)
1471 : end do
1472 : end do
1473 :
1474 3259 : ABI_FREE(b)
1475 3259 : ABI_FREE(c)
1476 :
1477 3259 : end subroutine hexc_interp_matmul
1478 : !!***
1479 :
1480 : !-------------------------------------------------------------------
1481 :
1482 : !!****f* m_hexc/hexc_matmul_tda
1483 : !! NAME
1484 : !! hexc_matmul_tda
1485 : !!
1486 : !! FUNCTION
1487 : !! Compute H_exc |\psi>
1488 : !!
1489 : !! INPUTS
1490 : !! hexc_i<hexc_interp_t> = Interpolated excitonic hamiltonian
1491 : !! phi = Input ket
1492 : !!
1493 : !! OUTPUT
1494 : !! hphi = hreso * phi
1495 : !!
1496 : !! SOURCE
1497 :
1498 9793 : subroutine hexc_matmul_tda(hexc, hexc_i, phi, hphi)
1499 :
1500 : !Arguments ---------------------------
1501 : class(hexc_t),intent(in) :: hexc
1502 : type(hexc_interp_t),intent(in) :: hexc_i
1503 : complex(dp),intent(in) :: phi(hexc%hsize)
1504 : complex(dp),intent(out) :: hphi(hexc%hsize)
1505 :
1506 : !Local variables ---------------------
1507 : integer :: ierr
1508 : real(dp) :: tsec(2)
1509 : !*****************************************************************************
1510 :
1511 9793 : call timab(697,1,tsec)
1512 :
1513 9793 : if (hexc%BSp%use_interp) then
1514 3352840 : hphi = hexc_i%diag_dense * phi
1515 :
1516 4360 : if (any(hexc%BSp%interp_mode == [2,3,4])) then
1517 : ! hphi = hphi + MATMUL(hinterp,phi)
1518 2211 : call xgemv('N',hexc_i%hsize_dense,hexc_i%hsize_dense,cone,hexc_i%hinterp,hexc_i%hsize_dense,phi,1,cone,hphi,1)
1519 2211 : if (any(hexc%BSp%interp_mode == [2,4])) then
1520 1101 : call timab(697,2,tsec)
1521 1101 : return ! We are done
1522 : end if
1523 : end if
1524 :
1525 : call hexc_interp_matmul(hexc%bsp, hexc%hsize_coarse, hexc_i%hsize_dense, hexc_i%all_hmat, phi, hphi, &
1526 3259 : hexc_i%interpolator%double_grid,hexc%nbnd_coarse, hexc_i%interpolator, hexc_i%div2kdense, hexc_i%kdense2div)
1527 :
1528 : else ! No interpolation
1529 5433 : call xgemv('N',hexc%hsize,hexc%my_nt,cone,hexc%hreso,hexc%hsize,phi,1,czero,hphi,1)
1530 5433 : call xmpi_sum(hphi,hexc%comm,ierr)
1531 : end if
1532 :
1533 8692 : call timab(697,2,tsec)
1534 :
1535 : end subroutine hexc_matmul_tda
1536 : !!***
1537 :
1538 : !-------------------------------------------------------------------
1539 :
1540 : !!****f* m_hexc/hexc_matmul_elphon
1541 : !! NAME
1542 : !! hexc_matmul_elphon
1543 : !!
1544 : !! FUNCTION
1545 : !! Compute H | \psi > + E_{elphon} | \psi >
1546 : !!
1547 : !! INPUTS
1548 : !! hexc<hexc_t> = Excitonic hamiltonian
1549 : !! phi = Input ket
1550 : !! ep_renorm = vector with electron-phonon renorms
1551 : !! op = 'N' for H | psi >, 'C' for H^\dagger | psi >
1552 : !!
1553 : !! OUTPUT
1554 : !! hphi = hreso * phi + ep_renorm * phi
1555 : !!
1556 : !! SOURCE
1557 :
1558 6432 : subroutine hexc_matmul_elphon(hexc, phi, hphi, op, ep_renorm)
1559 :
1560 : !Arguments ---------------------------
1561 : class(hexc_t),intent(in) :: hexc
1562 : character,intent(in) :: op
1563 : complex(dp),intent(in) :: phi(hexc%my_nt)
1564 : complex(dp),intent(out) :: hphi(hexc%hsize)
1565 : complex(dp),intent(in) :: ep_renorm(hexc%hsize)
1566 :
1567 : !Local variables ---------------------
1568 : integer :: ierr
1569 : real(dp) :: tsec(2)
1570 : !*****************************************************************************
1571 :
1572 6432 : call timab(697,1,tsec)
1573 :
1574 6432 : if(hexc%BSp%use_interp) then
1575 0 : ABI_ERROR('Not yet implemented with interpolation !')
1576 : else ! No interpolation
1577 : ! As our matrix is hermitian (hreso), we should always use 'N' here (it is stored column-wise !)
1578 6432 : call xgemv('N',hexc%hsize,hexc%my_nt,cone,hexc%hreso,hexc%hsize,phi,1,czero,hphi,1)
1579 :
1580 : !!! ep_renorm is stored on each cpu
1581 6432 : if (op == 'N') then
1582 311952 : hphi(hexc%my_t1:hexc%my_t2) = hphi(hexc%my_t1:hexc%my_t2) + ep_renorm(hexc%my_t1:hexc%my_t2) * phi
1583 3216 : else if(op == 'C') then
1584 311952 : hphi(hexc%my_t1:hexc%my_t2) = hphi(hexc%my_t1:hexc%my_t2) + CONJG(ep_renorm(hexc%my_t1:hexc%my_t2)) * phi
1585 : end if
1586 6432 : call xmpi_sum(hphi,hexc%comm,ierr)
1587 : end if
1588 :
1589 6432 : call timab(697,2,tsec)
1590 :
1591 6432 : end subroutine hexc_matmul_elphon
1592 : !!***
1593 :
1594 : !-------------------------------------------------------------------
1595 :
1596 : !!****f* m_hexc/hexc_matmul_full
1597 : !! NAME
1598 : !! hexc_matmul_full
1599 : !!
1600 : !! FUNCTION
1601 : !! Compute H | \psi >
1602 : !!
1603 : !! INPUTS
1604 : !! hexc<hexc_t> = Excitonic hamiltonian
1605 : !! hexc_i<hexc_interp_t> = Interpolated hamiltonian
1606 : !! phi = Input ket
1607 : !! parity = -1 or +1 parameter
1608 : !!
1609 : !! OUTPUT
1610 : !! hphi = hreso * phi + parity * hcoup * CONJ(phi)
1611 : !!
1612 : !! SOURCE
1613 :
1614 0 : subroutine hexc_matmul_full(hexc, hexc_i, phi, hphi, parity)
1615 :
1616 : !Arguments ---------------------------
1617 : class(hexc_t),intent(in) :: hexc
1618 : type(hexc_interp_t),intent(in) :: hexc_i
1619 : complex(dp),intent(in) :: phi(hexc%hsize)
1620 : complex(dp),intent(out) :: hphi(hexc%hsize)
1621 : integer,intent(in) :: parity
1622 :
1623 : !Local variables ---------------------
1624 : real(dp) :: tsec(2)
1625 : !*****************************************************************************
1626 :
1627 0 : call timab(697,1,tsec)
1628 :
1629 : ABI_UNUSED(hexc_i%hsize_dense)
1630 :
1631 0 : if(hexc%BSp%use_interp) then
1632 0 : ABI_ERROR("Coupling is not yet implemented with interpolation")
1633 : else
1634 : ! No interpolation
1635 0 : hphi = MATMUL(hexc%hreso,phi) + parity * MATMUL(hexc%hcoup,CONJG(phi))
1636 : end if
1637 :
1638 0 : call timab(697,2,tsec)
1639 :
1640 0 : end subroutine hexc_matmul_full
1641 : !!***
1642 :
1643 : !-------------------------------------------------------------------
1644 :
1645 0 : end module m_hexc
1646 : !!***
|
