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