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 : : ! v(q) = 4pi/q^2 for q /= 0 and i_sz for q = 0
771 [ + - + + : 6 : if (all(Vcp_dense%mode /= [character(len=len(Vcp_dense%mode)) :: 'CRYSTAL', 'AUXILIARY_FUNCTION', 'AUX_GB'])) then
- + ]
772 : 0 : ABI_ERROR('Vcp_dense%mode not implemented yet !')
773 : : end if
774 : : end if
775 : :
776 [ - + ]: 2 : if (BSp%nsppol > 1) then
777 : 0 : ABI_ERROR("nsppol > 1 not yet implemented")
778 : : end if
779 : :
780 : 2 : factor = one/grid%ndiv
781 : :
782 [ + + + + ]: 1181186 : hinterp = czero; term = czero
783 : :
784 [ - + + - : 10 : ABI_MALLOC_OR_DIE(Cmat,(nbnd_coarse,nbnd_coarse,interpolator%nvert), ierr)
- + + - -
+ + - + -
+ - + - -
+ - + ]
785 [ + + + + : 2514 : Cmat = czero
+ + ]
786 : :
787 [ - + + - : 6 : ABI_MALLOC(band2it,(nbnd_coarse))
+ - + - -
+ - + ]
788 : 2 : call interpolator%alloc_work(nbnd_coarse*interpolator%nvert)
789 : :
790 : 2 : btemp => interpolator%btemp
791 : 2 : ctemp => interpolator%ctemp
792 : :
793 : : if(newway) then
794 [ - + + - : 6 : ABI_MALLOC(tmp_Cmat,(nbnd_coarse))
+ - + - -
+ - + ]
795 [ - + + - : 8 : ABI_MALLOC(work_coeffs,(nbnd_coarse,interpolator%nvert))
- + + - +
- + - - +
- + ]
796 : : end if
797 : :
798 [ + + ]: 130 : do ik_dense = 1,grid%nbz_dense
799 : 128 : write(std_out,*) "Kdense = ",ik_dense,"/",grid%nbz_dense
800 : 128 : ik_coarse = grid%dense_to_coarse(ik_dense)
801 : : if(use_herm) then
802 : : limitnbz = ik_dense
803 : : else
804 : : limitnbz = grid%nbz_dense
805 : : end if
806 : :
807 [ + + ]: 4290 : do ikp_dense = 1,limitnbz
808 : 4160 : ikp_coarse = grid%dense_to_coarse(ikp_dense)
809 : :
810 [ + + ]: 16640 : do iv1 = BSp%lomo_spin(spin2),BSp%homo_spin(spin2)
811 [ + + ]: 66560 : do ic1 = BSp%lumo_spin(spin2),BSp%humo_spin(spin2)
812 : 49920 : itp_coarse1 = BSp%vcks2t(iv1,ic1,ikp_coarse,spin2)
813 : 49920 : ibndp_coarse1 = (iv1-BSp%lomo_spin(spin2))*BSp%maxnbndc+(ic1-BSp%lumo_spin(spin2)+1)
814 : :
815 : 62400 : band2it(ibndp_coarse1) = itp_coarse1
816 : : end do
817 : : end do
818 : :
819 : :
820 [ + - ]: 4160 : if (any(BSp%interp_mode == [2,3])) then
821 : : ! Check if we are along the diagonal
822 [ + + ]: 16640 : kmkp = Kmesh_dense%bz(:,ik_dense) - Kmesh_dense%bz(:,ikp_dense)
823 : :
824 [ + + ]: 16640 : call wrap2_pmhalf(kmkp(:),q2(:),shift(:))
825 [ + + + + : 104000 : qinred = MATMUL(grid%kptrlatt_coarse,q2)
+ + + + ]
826 : :
827 : : ! We are outside the diagonal
828 [ + + + - : 16640 : if (BSp%interp_mode==3 .and. ANY((ABS(qinred)-tol7) > m3_width)) cycle
+ + + - ]
829 : :
830 [ + + + + : 66560 : qnorm = two_pi*SQRT(DOT_PRODUCT(q2,MATMUL(gmet,q2)))
+ + ]
831 : :
832 [ + + + + : 5696 : if(ALL(ABS(q2(:)) < 1.e-3)) then
+ + ]
833 : 128 : vc_sqrt_qbz = SQRT(Vcp_dense%i_sz)
834 : : else
835 : 4032 : vc_sqrt_qbz = SQRT(four_pi/qnorm**2)
836 : : end if
837 : :
838 : : !!DEBUG CHK !
839 : : !!COMPUTE Qpoint
840 : : !call findqg0(iq_bz,g0,kmkp,Qmesh_dense%nbz,Qmesh_dense%bz,BSp%mG0)
841 : :
842 : : !! * Get iq_ibz, and symmetries from iq_bz
843 : : !call qmesh_dense%get_BZ_item(Qmesh_dense,iq_bz,qbz,iq_ibz,isym_q,itim_q)
844 : :
845 : : !if(iq_ibz > 1 .and. ABS(vc_sqrt_qbz - Vcp_dense%vc_sqrt(1,iq_ibz)) > 1.e-3) then
846 : : ! write(*,*) "vc_sqrt_qbz = ",vc_sqrt_qbz
847 : : ! write(*,*) "Vcp_dense%vc_sqrt(1,iq_ibz) = ",Vcp_dense%vc_sqrt(1,iq_ibz)
848 : : ! ABI_ERROR("vcp are not the same !")
849 : : !else if(iq_ibz == 1 .and. ABS(vc_sqrt_qbz - SQRT(Vcp_dense%i_sz)) > 1.e-3) then
850 : : ! write(*,*) "vc_sqrt_qbz = ",vc_sqrt_qbz
851 : : ! write(*,*) "SQRT(Vcp_dense%i_sz) = ",SQRT(Vcp_dense%i_sz)
852 : : ! ABI_ERROR("vcp are not the same !")
853 : : !end if
854 : : !!END DEBUG CHK !
855 : : end if
856 : :
857 : : if(newway) then
858 : :
859 [ + + + + : 5229120 : Cmat = czero
+ + ]
860 : :
861 [ + + + + ]: 436800 : work_coeffs = czero
862 : :
863 [ + + ]: 37440 : do ineighbour = 1,interpolator%nvert
864 : :
865 : : ! Loop over the (c, v) part of the left transition
866 [ + + ]: 137280 : do iv = BSp%lomo_spin(spin1),BSp%homo_spin(spin1)
867 [ + + ]: 532480 : do ic = BSp%lumo_spin(spin1),BSp%humo_spin(spin1)
868 : :
869 : 399360 : it_dense = BSp%vcks2t_interp(iv,ic,ik_dense,spin1)
870 : 399360 : it_coarse = BSp%vcks2t(iv,ic,ik_coarse,spin1)
871 : 399360 : ibnd_coarse = (iv-BSp%lomo_spin(spin1))*BSp%maxnbndc+(ic-BSp%lumo_spin(spin1)+1)
872 : :
873 : 399360 : itc = interpolator%corresp(it_coarse,ineighbour,spin1)
874 : :
875 [ + + - + ]: 399360 : if (any(BSp%interp_mode == [1,3,4])) then
876 : :
877 : : !work_coeffs(:,:) = hmat(itc,band2it(:),:)
878 [ + + ]: 1797120 : do inb = 1,interpolator%nvert
879 [ + - + + : 40135680 : work_coeffs(:,inb) = hmat(itc,interpolator%corresp(band2it(:),inb,spin2))
+ + ]
880 : : end do
881 : :
882 : : call hexc_compute_subhinterp(BSp,grid,nbnd_coarse,&
883 : : & interpolator,kdense2div,&
884 : : & work_coeffs,tmp_Cmat,ikp_dense,&
885 : 199680 : & interpolator%overlaps(:,:,:,ikp_dense,spin2))
886 : :
887 [ - + ]: 199680 : if(any(BSp%interp_mode == [1,4])) then
888 [ # # ]: 0 : Cmat(ibnd_coarse,:,ineighbour) = tmp_Cmat
889 [ + - ]: 199680 : else if (BSp%interp_mode == 3) then
890 [ + + ]: 2595840 : Cmat(ibnd_coarse,:,ineighbour) = -tmp_Cmat
891 : : end if
892 : : end if
893 : :
894 : :
895 [ + - ]: 399360 : if (any(BSp%interp_mode == [2,3])) then
896 : : !work_coeffs(:,:) = acoeffs(itc,band2it(:),:)
897 [ + + ]: 3594240 : do inb = 1,interpolator%nvert
898 [ + - + + : 80271360 : work_coeffs(:,inb) = acoeffs(itc,interpolator%corresp(band2it(:),inb,spin2))
+ + ]
899 : : end do
900 : :
901 : : call hexc_compute_subhinterp(BSp,grid,nbnd_coarse,&
902 : : & interpolator,kdense2div,&
903 : : & work_coeffs,tmp_Cmat,ikp_dense,&
904 : 399360 : & interpolator%overlaps(:,:,:,ikp_dense,spin2))
905 : :
906 [ + + ]: 5191680 : tmp_Cmat = tmp_Cmat * (vc_sqrt_qbz**2)
907 [ + + ]: 5191680 : Cmat(ibnd_coarse,:,ineighbour) = Cmat(ibnd_coarse,:,ineighbour) + tmp_Cmat
908 : : end if
909 : :
910 : :
911 [ + - ]: 399360 : if (any(BSp%interp_mode == [2,3])) then
912 : : !work_coeffs(:,:) = bcoeffs(itc,band2it(:),:)
913 [ + + ]: 3594240 : do inb = 1,interpolator%nvert
914 [ + - + + : 80271360 : work_coeffs(:,inb) = bcoeffs(itc,interpolator%corresp(band2it(:),inb,spin2))
+ + ]
915 : : end do
916 : :
917 : : call hexc_compute_subhinterp(BSp,grid,nbnd_coarse,&
918 : : & interpolator,kdense2div,&
919 : : & work_coeffs,tmp_Cmat,ikp_dense,&
920 : 399360 : & interpolator%overlaps(:,:,:,ikp_dense,spin2))
921 : :
922 [ + + ]: 5191680 : tmp_Cmat = tmp_Cmat * (vc_sqrt_qbz)
923 [ + + ]: 5191680 : Cmat(ibnd_coarse,:,ineighbour) = Cmat(ibnd_coarse,:,ineighbour) + tmp_Cmat
924 : : end if
925 : :
926 : :
927 [ + - ]: 499200 : if (any(BSp%interp_mode == [2,3])) then
928 : : !work_coeffs(:,:) = ccoeffs(itc,band2it(:),:)
929 [ + + ]: 3594240 : do inb = 1,interpolator%nvert
930 [ + - + + : 80271360 : work_coeffs(:,inb) = ccoeffs(itc,interpolator%corresp(band2it(:),inb,spin2))
+ + ]
931 : : end do
932 : :
933 : : call hexc_compute_subhinterp(BSp,grid,nbnd_coarse,&
934 : : & interpolator,kdense2div,&
935 : : & work_coeffs,tmp_Cmat,ikp_dense,&
936 : 399360 : & interpolator%overlaps(:,:,:,ikp_dense,spin2))
937 : :
938 [ + + ]: 5191680 : Cmat(ibnd_coarse,:,ineighbour) = Cmat(ibnd_coarse,:,ineighbour) + tmp_Cmat
939 : : end if
940 : : end do ! ic
941 : : end do ! iv
942 : : end do ! ineighbour
943 : :
944 : : else
945 : : ! Loop over the (c, v) part of the left transition
946 : : do iv = BSp%lomo_spin(spin1),BSp%homo_spin(spin1)
947 : : do ic = BSp%lumo_spin(spin1),BSp%humo_spin(spin1)
948 : :
949 : : it_dense = BSp%vcks2t_interp(iv,ic,ik_dense,spin1)
950 : : it_coarse = BSp%vcks2t(iv,ic,ik_coarse,spin1)
951 : : ibnd_coarse = (iv-BSp%lomo_spin(spin1))*BSp%maxnbndc+(ic-BSp%lumo_spin(spin1)+1)
952 : :
953 : : ! Loop over the (c', v') part of the right transition
954 : : do ivp = BSp%lomo_spin(spin2),BSp%homo_spin(spin2)
955 : : do icp = BSp%lumo_spin(spin2),BSp%humo_spin(spin2)
956 : :
957 : : itp_dense = BSp%vcks2t_interp(ivp,icp,ikp_dense,spin2)
958 : : itp_coarse = BSp%vcks2t(ivp,icp,ikp_coarse,spin2)
959 : : ibndp_coarse = (ivp-Bsp%lomo_spin(spin2))*BSp%maxnbndc+(icp-BSp%lumo_spin(spin2)+1)
960 : : ! Now we now it_dense, and itp_dense
961 : :
962 : : idivp = kdense2div(ikp_dense)
963 : :
964 : : btemp = czero; ctemp = czero
965 : :
966 : : ! MG TODO: This way of looping is not optimal
967 : : do ineighbour = 1,interpolator%nvert
968 : : itc = interpolator%corresp(it_coarse,ineighbour,spin1)
969 : :
970 : : do ineighbourp = 1,interpolator%nvert
971 : :
972 : : do iv1 = BSp%lomo_spin(spin2),BSp%homo_spin(spin2)
973 : : do ic1 = BSp%lumo_spin(spin2),BSp%humo_spin(spin2)
974 : :
975 : : ibndp_coarse1 = (iv1-BSp%lomo_spin(spin2))*BSp%maxnbndc+(ic1-BSp%lumo_spin(spin2)+1)
976 : : indwithnb = (ineighbourp-1)*nbnd_coarse+ibndp_coarse1
977 : : itp_coarse1 = BSp%vcks2t(iv1,ic1,ikp_coarse,spin2)
978 : : corresp_ind = interpolator%corresp(itp_coarse1,ineighbourp,spin2)
979 : :
980 : : select case (BSp%interp_mode)
981 : : case (1,4)
982 : : interpolator%btemp(indwithnb) = hmat(itc,corresp_ind)
983 : : case (2)
984 : : interpolator%btemp(indwithnb) = acoeffs(itc,corresp_ind)*(vc_sqrt_qbz**2) &
985 : : & + bcoeffs(itc,corresp_ind)*(vc_sqrt_qbz) &
986 : : & + ccoeffs(itc,corresp_ind)
987 : : case (3)
988 : : ! Diff between divergence and hmat
989 : : interpolator%btemp(indwithnb) = acoeffs(itc,corresp_ind)*(vc_sqrt_qbz**2) &
990 : : & + bcoeffs(itc,corresp_ind)*(vc_sqrt_qbz) &
991 : : & + ccoeffs(itc,corresp_ind) &
992 : : & - hmat(itc,corresp_ind)
993 : : case default
994 : : ABI_ERROR("Wrong Bsp%interp_mode")
995 : : end select
996 : :
997 : : ctemp(indwithnb) = &
998 : : & interpolator%overlaps(iv1,ivp,ineighbourp,ikp_dense,spin2) &
999 : : & * GWPC_CONJG(interpolator%overlaps(ic1,icp,ineighbourp,ikp_dense,spin2)) &
1000 : : & *interpolator%interp_factors(ineighbourp,idivp)
1001 : : end do ! ic1
1002 : : end do !iv1
1003 : :
1004 : : end do !ineighbourp
1005 : : Cmat(ibnd_coarse,ibndp_coarse,ineighbour) = xdotc(interpolator%nvert*nbnd_coarse,&
1006 : : & ctemp,1,btemp,1)
1007 : : end do !ineighbour
1008 : :
1009 : : end do !icp
1010 : : end do !ivp
1011 : :
1012 : : end do !ic
1013 : : end do !iv
1014 : :
1015 : : end if
1016 : :
1017 [ + + ]: 16768 : do iv = BSp%lomo_spin(spin1),BSp%homo_spin(spin1)
1018 [ + + ]: 66560 : do ic = BSp%lumo_spin(spin1),BSp%humo_spin(spin1)
1019 : 49920 : it_dense = BSp%vcks2t_interp(iv,ic,ik_dense,spin1)
1020 : 49920 : it_coarse = BSp%vcks2t(iv,ic,ik_coarse,spin1)
1021 : 49920 : ibnd_coarse = (iv-BSp%lomo_spin(spin1))*BSp%maxnbndc+(ic-BSp%lumo_spin(spin1)+1)
1022 : :
1023 : 49920 : idiv = kdense2div(ik_dense)
1024 : :
1025 [ + + ]: 212160 : do ivp = BSp%lomo_spin(spin2),BSp%homo_spin(spin2)
1026 [ + + ]: 798720 : do icp = BSp%lumo_spin(spin2),BSp%humo_spin(spin2)
1027 : 599040 : itp_dense = BSp%vcks2t_interp(ivp,icp,ikp_dense,spin2)
1028 : 599040 : itp_coarse = BSp%vcks2t(ivp,icp,ikp_coarse,spin2)
1029 : 599040 : ibndp_coarse = (ivp-Bsp%lomo_spin(spin2))*BSp%maxnbndc+(icp-BSp%lumo_spin(spin2)+1)
1030 : :
1031 : : ! Hermicity
1032 : : if(use_herm .and. it_dense < itp_dense) then
1033 : : continue
1034 : : end if
1035 : :
1036 : : !btemp = czero; ctemp = czero
1037 : :
1038 [ + + ]: 5391360 : do ineighbour = 1,interpolator%nvert
1039 [ + + ]: 19768320 : do iv1 = BSp%lomo_spin(spin1),BSp%homo_spin(spin1)
1040 [ + + ]: 76677120 : do ic1 = BSp%lumo_spin(spin1),BSp%humo_spin(spin1)
1041 : 57507840 : ibnd_coarse1 = (iv1-BSp%lomo_spin(spin1))*BSp%maxnbndc+(ic1-BSp%lumo_spin(spin1)+1)
1042 : 57507840 : it_dense1 = BSp%vcks2t_interp(iv1,ic1,ik_dense,spin1)
1043 : 57507840 : indwithnb = (ineighbour-1)*nbnd_coarse+ibnd_coarse1
1044 : :
1045 : 57507840 : btemp(indwithnb) = Cmat(ibnd_coarse1,ibndp_coarse,ineighbour)
1046 : :
1047 : : ctemp(indwithnb) = GWPC_CONJG(interpolator%overlaps(iv1,iv,ineighbour,ik_dense,spin1)) &
1048 : : & *interpolator%overlaps(ic1,ic,ineighbour,ik_dense,spin1) &
1049 : 71884800 : & *interpolator%interp_factors(ineighbour,idiv)
1050 : : end do !ic1
1051 : : end do !iv1
1052 : : end do !ineighbour
1053 : :
1054 : : ! Save interpolated value.
1055 : 748800 : hinterp(it_dense,itp_dense) = xdotc(interpolator%nvert*nbnd_coarse,ctemp,1,btemp,1)
1056 : : !DOT_PRODUCT(ctemp,btemp)
1057 : :
1058 : : end do !icp
1059 : : end do !ivp
1060 : :
1061 : : end do !ic
1062 : : end do !iv
1063 : :
1064 : : end do !ikp
1065 : : end do !ik
1066 : :
1067 : : ! Enforce hermiticity
1068 : : if(use_herm) then
1069 [ + + ]: 1538 : do itp_dense = 1,BSp%nreh_interp(spin2)
1070 [ + + ]: 592130 : do it_dense = itp_dense,BSp%nreh_interp(spin1)
1071 [ + + ]: 592128 : if(it_dense == itp_dense) then
1072 : 1536 : hinterp(itp_dense,it_dense) = DBLE(hinterp(it_dense,itp_dense))
1073 : : else
1074 : 589056 : hinterp(itp_dense,it_dense) = CONJG(hinterp(it_dense,itp_dense))
1075 : : end if
1076 : : end do
1077 : : end do
1078 : : end if
1079 : :
1080 [ - + ]: 2 : ABI_FREE(Cmat)
1081 [ - + ]: 2 : ABI_FREE(band2it)
1082 : :
1083 [ + - ]: 2 : ABI_SFREE(tmp_Cmat)
1084 [ + - ]: 2 : ABI_SFREE(work_coeffs)
1085 : :
1086 : 2 : nullify(ctemp)
1087 : 2 : nullify(btemp)
1088 : 2 : call interpolator%int_free()
1089 : :
1090 [ + + + + ]: 1181186 : hinterp = hinterp*factor
1091 : :
1092 : 2 : call timab(696,2,tsec)
1093 : :
1094 : 2 : end subroutine hexc_compute_hinterp
1095 : : !!***
1096 : :
1097 : : !----------------------------------------------------------------------
1098 : :
1099 : : !!****f* m_hexc/hexc_free
1100 : : !! NAME
1101 : : !! hexc_free
1102 : : !!
1103 : : !! FUNCTION
1104 : : !! Destroy the interpolator object in memory
1105 : : !!
1106 : : !! SOURCE
1107 : :
1108 : 23 : subroutine hexc_free(hexc)
1109 : :
1110 : : !Arguments ---------------------------
1111 : : class(hexc_t),intent(inout) :: hexc
1112 : : !*****************************************************************************
1113 : :
1114 [ + - ]: 23 : if (associated(hexc%bsp)) then
1115 : 23 : nullify(hexc%bsp)
1116 : : end if
1117 : :
1118 [ + - ]: 23 : if (associated(hexc%crystal)) then
1119 : 23 : nullify(hexc%crystal)
1120 : : end if
1121 : :
1122 [ + - ]: 23 : if (associated(hexc%kmesh_coarse)) then
1123 : 23 : nullify(hexc%kmesh_coarse)
1124 : : end if
1125 : :
1126 [ + - ]: 23 : ABI_SFREE(hexc%hreso)
1127 [ - + ]: 23 : ABI_SFREE(hexc%hcoup)
1128 [ - + ]: 23 : ABI_SFREE(hexc%diag_coarse)
1129 : :
1130 : 23 : end subroutine hexc_free
1131 : : !!***
1132 : :
1133 : : !-------------------------------------------------------------------
1134 : :
1135 : : !!****f* m_hexc/hexc_interp_free
1136 : : !! NAME
1137 : : !! hexc_interp_free
1138 : : !!
1139 : : !! FUNCTION
1140 : : !! Free dynamic memory.
1141 : : !!
1142 : : !! SOURCE
1143 : :
1144 : 23 : subroutine hexc_interp_free(hexc_i)
1145 : :
1146 : : !Arguments ---------------------------
1147 : : class(hexc_interp_t),intent(inout) :: hexc_i
1148 : : !*****************************************************************************
1149 : :
1150 [ + + ]: 23 : ABI_SFREE(hexc_i%kdense2div)
1151 [ + + ]: 23 : ABI_SFREE(hexc_i%div2kdense)
1152 [ + + ]: 23 : ABI_SFREE(hexc_i%diag_dense)
1153 [ + + ]: 23 : ABI_SFREE(hexc_i%hinterp)
1154 [ + + ]: 23 : ABI_SFREE(hexc_i%all_hmat)
1155 [ - + ]: 23 : ABI_SFREE(hexc_i%all_acoeffs)
1156 [ - + ]: 23 : ABI_SFREE(hexc_i%all_bcoeffs)
1157 [ - + ]: 23 : ABI_SFREE(hexc_i%all_ccoeffs)
1158 : :
1159 [ + + ]: 23 : if (associated(hexc_i%kmesh_dense)) then
1160 : 4 : nullify(hexc_i%kmesh_dense)
1161 : : end if
1162 : :
1163 [ + + ]: 23 : if (associated(hexc_i%vcp_dense)) then
1164 : 4 : nullify(hexc_i%vcp_dense)
1165 : : end if
1166 : :
1167 : 23 : call hexc_i%interpolator%free()
1168 : :
1169 : 23 : end subroutine hexc_interp_free
1170 : : !!***
1171 : :
1172 : : !----------------------------------------------------------------------
1173 : :
1174 : : !!****f* m_hexc/hexc_interp_matmul
1175 : : !! NAME
1176 : : !! hexc_interp_matmul
1177 : : !!
1178 : : !! FUNCTION
1179 : : !! Compute matrix-vector product Hmat * phi by interpolating coarse Hmat
1180 : : !!
1181 : : !! INPUTS
1182 : : !! BSp<type(excparam)>=Parameters defining the BS calculation
1183 : : !! hsize_coarse = Size of the coarse hamiltonian
1184 : : !! hsize_dense = Size of the dense hamiltonian
1185 : : !! hmat(hsize_coarse,hsize_coarse,8) = coarse hamiltonian
1186 : : !! phi(hsize_dense) = ket on which apply the matrix
1187 : : !! grid <double_grid_t> = Correspondence between coarse and dense k-mesh.
1188 : : !! nbnd_coarse = Total number of bands
1189 : : !! interpolator<interpolator_t> = Interpolator
1190 : : !! div2kdense = Mapping from coarse and division -> dense mesh
1191 : : !! kdense2div = Mapping from dense mesh -> division
1192 : : !!
1193 : : !! OUTPUT
1194 : : !! hphi(hsize_dense) = Interp(hmat)*phi
1195 : : !!
1196 : : !! SOURCE
1197 : :
1198 : 3263 : subroutine hexc_interp_matmul(BSp,hsize_coarse,hsize_dense,hmat,phi,hphi,grid,&
1199 : 3263 : nbnd_coarse,interpolator,div2kdense,kdense2div)
1200 : :
1201 : : !Arguments ------------------------------------
1202 : : !scalars
1203 : : integer,intent(in) :: hsize_coarse,hsize_dense,nbnd_coarse !,ntrans
1204 : : type(excparam),intent(in) :: BSp
1205 : : type(double_grid_t),intent(in) :: grid
1206 : : type(interpolator_t),intent(in) :: interpolator
1207 : : !arrays
1208 : : integer,intent(in) :: div2kdense(grid%nbz_coarse,grid%ndiv), kdense2div(grid%nbz_dense)
1209 : : complex(dp),intent(in) :: phi(hsize_dense)
1210 : : complex(dp),intent(in) :: hmat(hsize_coarse,hsize_coarse)
1211 : : complex(dp),intent(inout) :: hphi(hsize_dense)
1212 : :
1213 : : !Local variables ------------------------------
1214 : : !scalars
1215 : : integer :: itt,ik_dense,ik_coarse,it_coarse, ic,iv,iv1,ic1, ibnd_coarse
1216 : : integer :: ibnd_coarse1, ineighbour,idense,ikpt
1217 : : integer :: my_k1,my_k2,ind_with_nb,is, is1
1218 : : real(dp) :: factor
1219 : : complex(dp) :: tmp
1220 : : logical,parameter :: use_blas=.True.
1221 : : !arrays
1222 : 6526 : integer :: allindices(nbnd_coarse)
1223 : 6526 : complex(dp) :: allp(hsize_coarse,interpolator%nvert), test(hsize_coarse)
1224 : 6526 : complex(dp) :: ophi(grid%nbz_dense,interpolator%nvert,nbnd_coarse)
1225 : 3263 : complex(dp),allocatable :: b(:), c(:),A(:,:), tmp_array(:), tmp_array2(:,:)
1226 : : !************************************************************************
1227 : :
1228 : 3263 : factor = one/grid%ndiv
1229 : : !hphi = czero
1230 : :
1231 : : ! Outer index : k point in the dense zone
1232 : : ! Sum over vc
1233 : : ! Index of result : k point in the dense zone, v2,c2,neighbour
1234 : :
1235 : : ! Parallelization on nbz in the coarse mesh !
1236 : 3263 : my_k1 = 1
1237 : 3263 : my_k2 = grid%nbz_coarse
1238 : :
1239 [ - + + - : 13052 : ABI_MALLOC(A,(interpolator%nvert*nbnd_coarse,nbnd_coarse))
- + + - +
- + - - +
- + ]
1240 [ - + + - : 9789 : ABI_MALLOC(b,(nbnd_coarse))
+ - + - -
+ - + ]
1241 [ - + + - : 9789 : ABI_MALLOC(c,(interpolator%nvert*nbnd_coarse))
+ - + - -
+ + - ]
1242 : :
1243 [ + + + + : 14582842 : c = czero; ophi = czero
+ + + + ]
1244 : :
1245 [ + + ]: 212095 : do ik_dense = 1,grid%nbz_dense
1246 : : ! if( ik_dense is not in my set of k-points)
1247 : : ! ! continue
1248 : : !
1249 [ + + ]: 417664 : do is1 = 1, BSp%nsppol
1250 [ + + ]: 1044160 : do iv1 = BSp%lomo_spin(is1),Bsp%homo_spin(is1)
1251 [ + + ]: 3341312 : do ic1 = BSp%lumo_spin(is1),Bsp%humo_spin(is1)
1252 : 2505984 : ibnd_coarse = (iv1-BSp%lomo_spin(is1))*BSp%maxnbndc+(ic1-BSp%lumo_spin(is1)+1)
1253 : 2505984 : itt = BSp%vcks2t_interp(iv1,ic1,ik_dense,is1)
1254 : 3132480 : allindices(ibnd_coarse) = itt
1255 : : end do !ic1
1256 : : end do !iv1
1257 : : end do !is1
1258 : :
1259 [ + + ]: 2714816 : b(:) = phi(allindices(:))
1260 : :
1261 [ + + ]: 417664 : do is = 1, BSp%nsppol
1262 [ + + ]: 1044160 : do iv = BSp%lomo_spin(is),Bsp%homo_spin(is)
1263 [ + + ]: 3341312 : do ic = BSp%lumo_spin(is),Bsp%humo_spin(is)
1264 : 2505984 : ibnd_coarse = (iv-BSp%lomo_spin(is))*BSp%maxnbndc+(ic-BSp%lumo_spin(is)+1)
1265 : 2505984 : idense = Bsp%vcks2t_interp(iv,ic,ik_dense,is)
1266 : :
1267 [ + + ]: 17229120 : do ineighbour = 1,interpolator%nvert
1268 : 14096640 : ind_with_nb = (ineighbour-1)*(nbnd_coarse)+ibnd_coarse
1269 : :
1270 : : !A(ind_with_nb,:) = overlaps(allindices(:),ibnd_coarse,ineighbour)
1271 : :
1272 : : ! Should be optimized !!!
1273 [ + + ]: 58892544 : do iv1 = BSp%lomo_spin(is),Bsp%homo_spin(is)
1274 [ + + ]: 225546240 : do ic1 = BSp%lumo_spin(is),Bsp%humo_spin(is)
1275 : 169159680 : ibnd_coarse1 = (iv1-BSp%lomo_spin(is))*BSp%maxnbndc+(ic1-BSp%lumo_spin(is)+1)
1276 : : A(ind_with_nb,ibnd_coarse1) = GWPC_CONJG(interpolator%overlaps(iv,iv1,ineighbour,ik_dense,is)) &
1277 : 211449600 : & *interpolator%overlaps(ic,ic1,ineighbour,ik_dense,is)
1278 : : end do !ic1
1279 : : end do !iv1
1280 : : end do !ineighbour
1281 : : end do !ic
1282 : : end do !iv
1283 : : end do !is
1284 : :
1285 : : if(use_blas) then
1286 : 208832 : call xgemv('N',interpolator%nvert*nbnd_coarse,nbnd_coarse,cone,A,interpolator%nvert*nbnd_coarse,b,1,czero,c,1)
1287 : : else
1288 : : c = MATMUL(A,b)
1289 : : end if
1290 : :
1291 [ + + ]: 420927 : do is = 1, BSp%nsppol
1292 [ + + ]: 1044160 : do iv = BSp%lomo_spin(is),BSp%homo_spin(is)
1293 [ + + ]: 3341312 : do ic = BSp%lumo_spin(is),BSp%humo_spin(is)
1294 : 2505984 : ibnd_coarse = (iv-BSp%lomo_spin(is))*BSp%maxnbndc+(ic-BSp%lumo_spin(is)+1)
1295 [ + + ]: 17229120 : do ineighbour = 1,interpolator%nvert
1296 : 14096640 : ind_with_nb = (ineighbour-1)*(nbnd_coarse)+ibnd_coarse
1297 : 16602624 : ophi(ik_dense,ineighbour,ibnd_coarse) = c(ind_with_nb)
1298 : : end do !ineighbour
1299 : : end do !ic
1300 : : end do !iv
1301 : : end do !is
1302 : :
1303 : : end do !ik_dense
1304 : :
1305 [ - + ]: 3263 : ABI_FREE(A)
1306 [ - + ]: 3263 : ABI_FREE(b)
1307 [ - + ]: 3263 : ABI_FREE(c)
1308 : :
1309 : : !call xmpi_sum_(ophi,comm,ierr)
1310 : :
1311 : : ! Outer index : k,v,c in the coarse zone, ineighbour
1312 : : ! Sum over all k-dense relative to one coarse point
1313 : : ! Index of result : k,v,c in the coarse zone, ineighbour
1314 : :
1315 [ - + + - : 9789 : ABI_MALLOC(b,(grid%ndiv))
+ - + - -
+ - + ]
1316 [ - + + - : 6526 : ABI_MALLOC(c,(grid%ndiv))
- + + - ]
1317 : :
1318 [ + + + + ]: 1783698 : allp = czero
1319 : :
1320 [ + + ]: 6526 : do is = 1, BSp%nsppol
1321 [ + + ]: 24881 : do ineighbour = 1,interpolator%nvert
1322 : :
1323 [ + + ]: 1783698 : do it_coarse = 1, BSp%nreh(is)
1324 : : ibnd_coarse = (Bsp%trans(it_coarse,is)%v-BSp%lomo_spin(is))*BSp%maxnbndc+&
1325 : 1762080 : (BSp%Trans(it_coarse,is)%c-BSp%lumo_spin(is)+1)
1326 : 1762080 : ik_coarse = BSp%trans(it_coarse,is)%k
1327 : : !b(:) = interp_factors(it_coarse,ineighbour,:)
1328 [ + + ]: 15858720 : b(:) = interpolator%interp_factors(ineighbour,:)
1329 : : !c(:) = ophi(indices(it_coarse,:),ineighbour,ibnd_coarse)
1330 [ + + ]: 15858720 : c(:) = ophi(div2kdense(ik_coarse,:),ineighbour,ibnd_coarse)
1331 [ + + ]: 15858720 : tmp = DOT_PRODUCT(b,c)
1332 : 1780435 : allp(it_coarse,ineighbour) = tmp
1333 : : end do
1334 : :
1335 : : end do
1336 : : end do
1337 : :
1338 : : !call xmpi_sum_(allp,comm,ierr)
1339 : :
1340 [ - + ]: 3263 : ABI_FREE(b)
1341 [ - + ]: 3263 : ABI_FREE(c)
1342 : :
1343 [ - + + - : 9789 : ABI_MALLOC(tmp_array,(hsize_coarse))
+ - + - -
+ - + ]
1344 [ - + + - : 13052 : ABI_MALLOC(tmp_array2,(hsize_coarse,hsize_coarse))
- + + - +
- + - - +
+ - ]
1345 [ + + ]: 316511 : tmp_array(:) = czero
1346 [ + + + + ]: 30388319 : tmp_array2(:,:) = czero
1347 : :
1348 [ + + ]: 316511 : test = czero
1349 : :
1350 : : ! Second step : Multiplication by hmat
1351 [ + + ]: 21618 : do ineighbour = 1,interpolator%nvert
1352 [ + - + + : 172720550 : tmp_array2 = factor*hmat(:,interpolator%corresp(:,ineighbour,1)) ! 1 is for spin
+ - - + -
- + + +
+ ]
1353 : 3263 : if(use_blas) then
1354 : : !call xgemv('N',hsize_coarse,hsize_coarse,cone,factor*(hmat(:,:,ineighbour)),hsize_coarse,allp(:,ineighbour),1,czero,tmp_array,1)
1355 : : !tmp_array2 = hmat(:,:,ineighbour)
1356 : : !tmp_array2 = factor*tmp_array2
1357 : 18355 : call xgemv('N',hsize_coarse,hsize_coarse,cone,tmp_array2,hsize_coarse,allp(:,ineighbour),1,czero,tmp_array,1)
1358 [ + + ]: 1780435 : test = test + tmp_array
1359 : : else
1360 : : !test = test+MATMUL(factor*(hmat(:,interpolator%corresp(:,ineighbour,1))),allp(:,ineighbour))
1361 : : test = test+MATMUL(tmp_array2,allp(:,ineighbour))
1362 : : end if
1363 : : end do
1364 : :
1365 [ - + ]: 3263 : ABI_FREE(tmp_array)
1366 [ - + ]: 3263 : ABI_FREE(tmp_array2)
1367 : :
1368 : : ! Outer index : ineighbour
1369 : : ! Sum over all v c
1370 : : ! Index of result : ineighbour, k_dense, v,c
1371 [ - + + - : 13052 : ABI_MALLOC(A,(nbnd_coarse,nbnd_coarse))
- + + - +
- + - - +
- + ]
1372 [ - + + - : 6526 : ABI_MALLOC(b,(nbnd_coarse))
- + - + ]
1373 [ - + + - : 6526 : ABI_MALLOC(c,(nbnd_coarse))
- + + - ]
1374 [ + + ]: 42419 : c = czero
1375 : :
1376 [ + + ]: 21618 : do ineighbour = 1,interpolator%nvert
1377 [ + + ]: 1196338 : do ik_dense = 1,grid%nbz_dense
1378 : :
1379 [ + + ]: 2349440 : do is1 = 1, Bsp%nsppol
1380 [ + + ]: 5873600 : do iv1 = Bsp%lomo_spin(is1),Bsp%homo_spin(is1)
1381 [ + + ]: 18795520 : do ic1 = BSp%lumo_spin(is1), Bsp%humo_spin(is1)
1382 : 14096640 : ibnd_coarse = (iv1-BSp%lomo_spin(is1))*BSp%maxnbndc+(ic1-BSp%lumo_spin(is1)+1)
1383 : :
1384 : 14096640 : ik_coarse = grid%dense_to_coarse(ik_dense)
1385 : 14096640 : itt = BSp%vcks2t(iv1,ic1,ik_coarse,is1)
1386 : 17620800 : b(ibnd_coarse) = test(interpolator%corresp(itt,ineighbour,is1))
1387 : : end do ! ic1
1388 : : end do ! iv1
1389 : : end do ! is1
1390 : :
1391 [ + + ]: 2349440 : do is = 1, BSp%nsppol
1392 [ + + ]: 5873600 : do iv = BSp%lomo_spin(is),Bsp%homo_spin(is)
1393 [ + + ]: 18795520 : do ic = BSp%lumo_spin(is),BSp%humo_spin(is)
1394 : 14096640 : ibnd_coarse = (iv-BSp%lomo_spin(is))*Bsp%maxnbndc+(ic-BSp%lumo_spin(is)+1)
1395 : 14096640 : idense = BSp%vcks2t_interp(iv,ic,ik_dense,is)
1396 : :
1397 : : !A(ibnd_coarse,:) = CONJG(overlaps(idense,:,ineighbour))
1398 : :
1399 : : ! Should be optimized !!!
1400 [ + + ]: 59910720 : do iv1 = BSp%lomo_spin(is),Bsp%homo_spin(is)
1401 [ + + ]: 225546240 : do ic1 = BSp%lumo_spin(is),Bsp%humo_spin(is)
1402 : 169159680 : ibnd_coarse1 = (iv1-BSp%lomo_spin(is))*BSp%maxnbndc+(ic1-BSp%lumo_spin(is)+1)
1403 : : A(ibnd_coarse,ibnd_coarse1) = (interpolator%overlaps(iv1,iv,ineighbour,ik_dense,is)) &
1404 : 211449600 : *GWPC_CONJG(interpolator%overlaps(ic1,ic,ineighbour,ik_dense,is))
1405 : : end do !ic1
1406 : : end do !iv1
1407 : : end do ! ic
1408 : : end do !iv
1409 : : end do !is
1410 : :
1411 : : if(use_blas) then
1412 : 1174720 : call xgemv('N',nbnd_coarse,nbnd_coarse,cone,A,nbnd_coarse,b,1,czero,c,1)
1413 : : else
1414 : : c = MATMUL(A,b)
1415 : : end if
1416 : :
1417 [ + + ]: 2367795 : do is = 1, BSp%nsppol
1418 [ + + ]: 5873600 : do iv = BSp%lomo_spin(is),Bsp%homo_spin(is)
1419 [ + + ]: 18795520 : do ic = BSp%lumo_spin(is),BSp%humo_spin(is)
1420 : 14096640 : ibnd_coarse = (iv-BSp%lomo_spin(is))*BSp%maxnbndc+(ic-BSp%lumo_spin(is)+1)
1421 : 14096640 : idense = Bsp%vcks2t_interp(iv,ic,ik_dense,is)
1422 : : !ophi(ik_dense,ineighbour,ibnd_coarse) = c(idense)
1423 : 17620800 : ophi(ik_dense,ineighbour,ibnd_coarse) = c(ibnd_coarse)
1424 : : end do
1425 : : end do
1426 : : end do
1427 : :
1428 : : end do ! ik_dense
1429 : : end do ! ineighbour
1430 : :
1431 : : !call xmpi_sum_(ophi,comm,ierr)
1432 : :
1433 [ - + ]: 3263 : ABI_FREE(A)
1434 [ - + ]: 3263 : ABI_FREE(b)
1435 [ - + ]: 3263 : ABI_FREE(c)
1436 : :
1437 : : ! Outer indices : it_dense
1438 : : ! Sum over neighbours
1439 : : ! Index of result : it_dense (ik,ic,iv dense)
1440 : :
1441 [ - + + - : 9789 : ABI_MALLOC(b,(interpolator%nvert))
+ - + - -
+ - + ]
1442 [ - + + - : 6526 : ABI_MALLOC(c,(interpolator%nvert))
- + + - ]
1443 : :
1444 [ + + ]: 6526 : do is = 1, BSp%nsppol
1445 [ + + ]: 2512510 : do itt = 1,BSp%nreh_interp(is)
1446 : : ! From itt -> ik_ibz,ic,iv
1447 : 2505984 : ik_dense = BSp%Trans_interp(itt,is)%k
1448 : 2505984 : ic = BSp%Trans_interp(itt,is)%c
1449 : 2505984 : iv = BSp%Trans_interp(itt,is)%v
1450 : :
1451 : : ! From ik_ibz in the dense mesh -> indices_dense
1452 : 2505984 : ik_coarse = grid%dense_to_coarse(ik_dense)
1453 : 2505984 : it_coarse = BSp%vcks2t(iv,ic,ik_coarse,is)
1454 : :
1455 : 2505984 : ibnd_coarse = (iv-BSp%lomo_spin(is))*BSp%maxnbndc+(ic-BSp%lumo_spin(is)+1)
1456 : :
1457 : 2505984 : ikpt = kdense2div(ik_dense)
1458 : : !ikpt = -1
1459 : : !do ix = 1,grid%ndiv
1460 : : ! if (indices(it_coarse,ix) == ik_dense) then
1461 : : ! ikpt = ix
1462 : : ! exit
1463 : : ! end if
1464 : : !end do
1465 : : !ABI_CHECK(ikpt/=-1,"Cannot find ik_dense")
1466 : :
1467 : : !b = interp_factors(it_coarse,:,ikpt)
1468 [ + - - + : 19108608 : b = interpolator%interp_factors(:,ikpt)
- - + + ]
1469 [ + - - + : 19108608 : c = ophi(ik_dense,:,ibnd_coarse)
- - + + ]
1470 : :
1471 : 2509247 : hphi(itt) = hphi(itt) + xdotc(interpolator%nvert, b, 1, c, 1)
1472 : : end do
1473 : : end do
1474 : :
1475 [ - + ]: 3263 : ABI_FREE(b)
1476 [ - + ]: 3263 : ABI_FREE(c)
1477 : :
1478 : 3263 : end subroutine hexc_interp_matmul
1479 : : !!***
1480 : :
1481 : : !-------------------------------------------------------------------
1482 : :
1483 : : !!****f* m_hexc/hexc_matmul_tda
1484 : : !! NAME
1485 : : !! hexc_matmul_tda
1486 : : !!
1487 : : !! FUNCTION
1488 : : !! Compute H_exc |\psi>
1489 : : !!
1490 : : !! INPUTS
1491 : : !! hexc_i<hexc_interp_t> = Interpolated excitonic hamiltonian
1492 : : !! phi = Input ket
1493 : : !!
1494 : : !! OUTPUT
1495 : : !! hphi = hreso * phi
1496 : : !!
1497 : : !! SOURCE
1498 : :
1499 : 9809 : subroutine hexc_matmul_tda(hexc, hexc_i, phi, hphi)
1500 : :
1501 : : !Arguments ---------------------------
1502 : : class(hexc_t),intent(in) :: hexc
1503 : : type(hexc_interp_t),intent(in) :: hexc_i
1504 : : complex(dp),intent(in) :: phi(hexc%hsize)
1505 : : complex(dp),intent(out) :: hphi(hexc%hsize)
1506 : :
1507 : : !Local variables ---------------------
1508 : : integer :: ierr
1509 : : real(dp) :: tsec(2)
1510 : : !*****************************************************************************
1511 : :
1512 : 9809 : call timab(697,1,tsec)
1513 : :
1514 [ + + ]: 9809 : if (hexc%BSp%use_interp) then
1515 [ + + ]: 3361299 : hphi = hexc_i%diag_dense * phi
1516 : :
1517 [ + + ]: 4371 : if (any(hexc%BSp%interp_mode == [2,3,4])) then
1518 : : ! hphi = hphi + MATMUL(hinterp,phi)
1519 : 2213 : call xgemv('N',hexc_i%hsize_dense,hexc_i%hsize_dense,cone,hexc_i%hinterp,hexc_i%hsize_dense,phi,1,cone,hphi,1)
1520 [ + + ]: 2213 : if (any(hexc%BSp%interp_mode == [2,4])) then
1521 : 1108 : call timab(697,2,tsec)
1522 : 1108 : return ! We are done
1523 : : end if
1524 : : end if
1525 : :
1526 : : call hexc_interp_matmul(hexc%bsp, hexc%hsize_coarse, hexc_i%hsize_dense, hexc_i%all_hmat, phi, hphi, &
1527 : 3263 : hexc_i%interpolator%double_grid,hexc%nbnd_coarse, hexc_i%interpolator, hexc_i%div2kdense, hexc_i%kdense2div)
1528 : :
1529 : : else ! No interpolation
1530 : 5438 : call xgemv('N',hexc%hsize,hexc%my_nt,cone,hexc%hreso,hexc%hsize,phi,1,czero,hphi,1)
1531 : 5438 : call xmpi_sum(hphi,hexc%comm,ierr)
1532 : : end if
1533 : :
1534 : 8701 : call timab(697,2,tsec)
1535 : :
1536 : : end subroutine hexc_matmul_tda
1537 : : !!***
1538 : :
1539 : : !-------------------------------------------------------------------
1540 : :
1541 : : !!****f* m_hexc/hexc_matmul_elphon
1542 : : !! NAME
1543 : : !! hexc_matmul_elphon
1544 : : !!
1545 : : !! FUNCTION
1546 : : !! Compute H | \psi > + E_{elphon} | \psi >
1547 : : !!
1548 : : !! INPUTS
1549 : : !! hexc<hexc_t> = Excitonic hamiltonian
1550 : : !! phi = Input ket
1551 : : !! ep_renorm = vector with electron-phonon renorms
1552 : : !! op = 'N' for H | psi >, 'C' for H^\dagger | psi >
1553 : : !!
1554 : : !! OUTPUT
1555 : : !! hphi = hreso * phi + ep_renorm * phi
1556 : : !!
1557 : : !! SOURCE
1558 : :
1559 : 6466 : subroutine hexc_matmul_elphon(hexc, phi, hphi, op, ep_renorm)
1560 : :
1561 : : !Arguments ---------------------------
1562 : : class(hexc_t),intent(in) :: hexc
1563 : : character,intent(in) :: op
1564 : : complex(dp),intent(in) :: phi(hexc%my_nt)
1565 : : complex(dp),intent(out) :: hphi(hexc%hsize)
1566 : : complex(dp),intent(in) :: ep_renorm(hexc%hsize)
1567 : :
1568 : : !Local variables ---------------------
1569 : : integer :: ierr
1570 : : real(dp) :: tsec(2)
1571 : : !*****************************************************************************
1572 : :
1573 : 6466 : call timab(697,1,tsec)
1574 : :
1575 [ - + ]: 6466 : if(hexc%BSp%use_interp) then
1576 : 0 : ABI_ERROR('Not yet implemented with interpolation !')
1577 : : else ! No interpolation
1578 : : ! As our matrix is hermitian (hreso), we should always use 'N' here (it is stored column-wise !)
1579 : 6466 : call xgemv('N',hexc%hsize,hexc%my_nt,cone,hexc%hreso,hexc%hsize,phi,1,czero,hphi,1)
1580 : :
1581 : : !!! ep_renorm is stored on each cpu
1582 [ + + ]: 6466 : if (op == 'N') then
1583 [ + + ]: 313601 : hphi(hexc%my_t1:hexc%my_t2) = hphi(hexc%my_t1:hexc%my_t2) + ep_renorm(hexc%my_t1:hexc%my_t2) * phi
1584 [ + - ]: 3233 : else if(op == 'C') then
1585 [ + + ]: 313601 : hphi(hexc%my_t1:hexc%my_t2) = hphi(hexc%my_t1:hexc%my_t2) + CONJG(ep_renorm(hexc%my_t1:hexc%my_t2)) * phi
1586 : : end if
1587 : 6466 : call xmpi_sum(hphi,hexc%comm,ierr)
1588 : : end if
1589 : :
1590 : 6466 : call timab(697,2,tsec)
1591 : :
1592 : 6466 : end subroutine hexc_matmul_elphon
1593 : : !!***
1594 : :
1595 : : !-------------------------------------------------------------------
1596 : :
1597 : : !!****f* m_hexc/hexc_matmul_full
1598 : : !! NAME
1599 : : !! hexc_matmul_full
1600 : : !!
1601 : : !! FUNCTION
1602 : : !! Compute H | \psi >
1603 : : !!
1604 : : !! INPUTS
1605 : : !! hexc<hexc_t> = Excitonic hamiltonian
1606 : : !! hexc_i<hexc_interp_t> = Interpolated hamiltonian
1607 : : !! phi = Input ket
1608 : : !! parity = -1 or +1 parameter
1609 : : !!
1610 : : !! OUTPUT
1611 : : !! hphi = hreso * phi + parity * hcoup * CONJ(phi)
1612 : : !!
1613 : : !! SOURCE
1614 : :
1615 : 0 : subroutine hexc_matmul_full(hexc, hexc_i, phi, hphi, parity)
1616 : :
1617 : : !Arguments ---------------------------
1618 : : class(hexc_t),intent(in) :: hexc
1619 : : type(hexc_interp_t),intent(in) :: hexc_i
1620 : : complex(dp),intent(in) :: phi(hexc%hsize)
1621 : : complex(dp),intent(out) :: hphi(hexc%hsize)
1622 : : integer,intent(in) :: parity
1623 : :
1624 : : !Local variables ---------------------
1625 : : real(dp) :: tsec(2)
1626 : : !*****************************************************************************
1627 : :
1628 : 0 : call timab(697,1,tsec)
1629 : :
1630 : : ABI_UNUSED(hexc_i%hsize_dense)
1631 : :
1632 [ # # ]: 0 : if(hexc%BSp%use_interp) then
1633 : 0 : ABI_ERROR("Coupling is not yet implemented with interpolation")
1634 : : else
1635 : : ! No interpolation
1636 [ # # # # : 0 : hphi = MATMUL(hexc%hreso,phi) + parity * MATMUL(hexc%hcoup,CONJG(phi))
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # #
# ]
1637 : : end if
1638 : :
1639 : 0 : call timab(697,2,tsec)
1640 : :
1641 : 0 : end subroutine hexc_matmul_full
1642 : : !!***
1643 : :
1644 : : !-------------------------------------------------------------------
1645 : :
1646 [ # # # # : 0 : end module m_hexc
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # ]
1647 : : !!***
|
