Branch data Line data Source code
1 : : !!****m* ABINIT/m_bseinterp
2 : : !! NAME
3 : : !! m_bseinterp
4 : : !!
5 : : !! FUNCTION
6 : : !!
7 : : !! COPYRIGHT
8 : : !! Copyright (C) 2014-2026 ABINIT group (Y. Gillet, M.Giantomassi)
9 : : !! This file is distributed under the terms of the
10 : : !! GNU General Public License, see ~abinit/COPYING
11 : : !! or http://www.gnu.org/copyleft/gpl.txt .
12 : : !!
13 : : !! SOURCE
14 : :
15 : : #if defined HAVE_CONFIG_H
16 : : #include "config.h"
17 : : #endif
18 : :
19 : : #include "abi_common.h"
20 : :
21 : : MODULE m_bseinterp
22 : :
23 : : use defs_basis
24 : : use m_abicore
25 : : use m_bs_defs
26 : : use m_xmpi
27 : : use m_errors
28 : : use m_nctk
29 : : use m_haydock_io
30 : : use m_linalg_interfaces
31 : : use netcdf
32 : :
33 : : use m_fstrings, only : indent, strcat, sjoin, itoa
34 : : use defs_datatypes, only : pseudopotential_type
35 : : use m_hide_blas, only : xdotc
36 : : use m_fft_mesh, only : calc_ceigr
37 : : use m_crystal, only : crystal_t
38 : : use m_bz_mesh, only : kmesh_t
39 : : use m_double_grid, only : double_grid_t, get_kpt_from_indices_coarse
40 : : use m_wfd, only : wfdgw_t
41 : : use m_pawtab, only : pawtab_type
42 : :
43 : : implicit none
44 : :
45 : : private
46 : : !!***
47 : :
48 : : !----------------------------------------------------------------------
49 : :
50 : : !!****t* m_haydock/interpolator_t
51 : : !! NAME
52 : : !! interpolator_t
53 : : !!
54 : : !! FUNCTION
55 : : !! Store the overlap matrix elements needed for the interpolation of the BSE Hamiltonian
56 : : !!
57 : : !! TODO
58 : : !! Decide if we want to make the number of bands k-dependent.
59 : : !!
60 : : !! SOURCE
61 : :
62 : : type,public :: interpolator_t
63 : :
64 : : integer :: nvert=8
65 : : ! Number of vertices for interpolation
66 : :
67 : : integer :: method
68 : : ! Interpolation method (YG or Rohlfing & Louie or ...)
69 : :
70 : : integer :: mband_dense, mband_coarse
71 : : ! Max number of bands dense and coarse
72 : :
73 : : integer :: nsppol
74 : : ! Number of spin channels
75 : :
76 : : integer, allocatable :: corresp(:,:,:)
77 : : ! corresp(max_nreh,nvert,spin)
78 : : ! it_coarse, idiv -> it_coarse (idiv-th neighbour)
79 : :
80 : : real(dp),allocatable :: interp_factors(:,:)
81 : : ! interp_factors(nvert,ndiv)
82 : : ! index_in_fine_box -> k-point in Trans_interp
83 : :
84 : : complex(gwp),allocatable :: overlaps(:,:,:,:,:)
85 : : ! Overlaps between dense and coarse mesh
86 : : ! overlaps(mband_coarse,mband_dense,ivertex_coarse,double_grid%nkpt_dense,spin)
87 : :
88 : : complex(dp),allocatable :: btemp(:), ctemp(:)
89 : : ! Temporary arrays for work
90 : :
91 : : ! Pointers to datatypes that are already in memory
92 : : type(double_grid_t),pointer :: double_grid => null()
93 : : ! Mapping between coarse and dense mesh
94 : :
95 : : contains
96 : : procedure :: init => interpolator_init ! Construct the object
97 : : procedure :: free => interpolator_free ! Free memory
98 : : procedure :: normalize => interpolator_normalize ! Normalize the overlaps
99 : : procedure :: alloc_work => int_alloc_work ! Alloc temp memory
100 : : procedure :: int_free => int_free ! Free temp memory
101 : :
102 : : end type interpolator_t
103 : : !!***
104 : :
105 : : !----------------------------------------------------------------------
106 : :
107 : : CONTAINS !=======================================================================
108 : : !!***
109 : :
110 : : !!****f* m_bseinterp/interpolator_init
111 : : !! NAME
112 : : !! interpolator_init
113 : : !!
114 : : !! FUNCTION
115 : : !! Construct the interpolator object
116 : : !!
117 : : !! INPUTS
118 : : !!
119 : : !! OUTPUT
120 : : !!
121 : : !! SOURCE
122 : :
123 : 4 : subroutine interpolator_init(interpolator, double_grid, Wfd_dense, Wfd_coarse, &
124 : 4 : Kmesh_dense, Kmesh_coarse, BSp, Cryst, Psps, Pawtab, method)
125 : :
126 : : !Arguments ---------------------------
127 : : !scalars
128 : : class(interpolator_t),intent(inout) :: interpolator
129 : : integer,intent(in) :: method
130 : : type(double_grid_t),intent(in),target :: double_grid
131 : : type(wfdgw_t),intent(inout) :: Wfd_dense, Wfd_coarse
132 : : type(kmesh_t),intent(in) :: Kmesh_dense, Kmesh_coarse
133 : : type(excparam),intent(in) :: BSp
134 : : type(crystal_t),intent(in) :: Cryst
135 : : type(pseudopotential_type),intent(in) :: Psps
136 : : !arrays
137 : : type(pawtab_type),intent(in) :: Pawtab(Cryst%ntypat*Wfd_coarse%usepaw)
138 : :
139 : : !Local variables ---------------------
140 : : !scalars
141 : : integer :: nsppol, nvert
142 : : integer :: maxnreh, nreh1, nreh2
143 : : integer :: mbandc, mbandd, nbzd
144 : : real(dp),parameter :: threshold = 0.1_dp
145 : : !arrays
146 : : character(len=500) :: msg
147 : :
148 : : !*****************************************************************************
149 : :
150 [ - + ]: 4 : ABI_CHECK(Wfd_coarse%usepaw==0, "PAW not yet supported")
151 [ - + ]: 4 : ABI_CHECK(BSp%nsppol==1, "nsppol != 1 not yet implemented")
152 [ - + ]: 4 : ABI_CHECK(Wfd_coarse%nspinor==1, "nspinor != 1 not supported")
153 : :
154 : : ABI_UNUSED(Pawtab(1)%basis_size)
155 : : !paw_overlap(cprj1,cprj2,typat,pawtab,spinor_comm) result(onsite)
156 : :
157 : 4 : interpolator%double_grid => double_grid
158 : 4 : interpolator%mband_dense = Wfd_dense%mband
159 : 4 : interpolator%mband_coarse = Wfd_coarse%mband
160 : 4 : interpolator%method = method
161 [ + - ]: 4 : interpolator%nsppol = BSp%nsppol
162 : :
163 : : SELECT CASE(method)
164 : : CASE (BSE_INTERP_YG)
165 : : nvert = 8
166 : : CASE (BSE_INTERP_RL2)
167 : : nvert = 2
168 : : CASE (BSE_INTERP_RL)
169 : 0 : nvert = 1
170 : : CASE DEFAULT
171 : 0 : write(msg,'(a,i0)') "Wrong interpolation method: ",method
172 : 0 : ABI_ERROR(msg)
173 : : END SELECT
174 : :
175 : 4 : interpolator%nvert = nvert
176 : :
177 : 4 : mbandc = interpolator%mband_coarse
178 : 4 : mbandd = interpolator%mband_dense
179 : 4 : nbzd = double_grid%nbz_dense
180 : 4 : nsppol = interpolator%nsppol
181 [ - + + - : 28 : ABI_MALLOC(interpolator%overlaps,(mbandc,mbandd,nvert,nbzd,nsppol))
- + + - -
+ + - - +
+ - - + +
- + - + -
- + - + -
+ ]
182 : :
183 : : call int_compute_overlaps(interpolator,double_grid, Wfd_dense, Wfd_coarse, Kmesh_dense, &
184 : 4 : & Kmesh_coarse, BSp, Cryst, Psps, Pawtab)
185 : :
186 [ - + + - : 16 : ABI_MALLOC(interpolator%interp_factors,(nvert,double_grid%ndiv))
- + + - +
- + - - +
- + - + ]
187 : :
188 : 4 : call int_preprocess_tables(interpolator,double_grid)
189 : :
190 : 4 : nreh1 = BSp%nreh(1)
191 [ - + ]: 4 : nreh2 = nreh1; if(BSp%nsppol == 2) nreh2 = BSp%nreh(2)
192 : 4 : maxnreh = MAX(nreh1,nreh2)
193 : :
194 [ - + + - : 20 : ABI_MALLOC(interpolator%corresp,(maxnreh,interpolator%nvert,interpolator%nsppol))
- + + - -
+ + - + -
+ - - + -
+ - + ]
195 : :
196 : 4 : call int_compute_corresp(interpolator,BSp,double_grid)
197 : :
198 : 4 : end subroutine interpolator_init
199 : : !!***
200 : :
201 : : !-------------------------------------------------------------------
202 : :
203 : : !!****f* m_bseinterp/int_alloc_work
204 : : !! NAME
205 : : !! int_alloc_work
206 : : !!
207 : : !! FUNCTION
208 : : !! Allocate temporary arrays
209 : : !!
210 : : !! INPUTS
211 : : !!
212 : : !! SOURCE
213 : :
214 : 2 : subroutine int_alloc_work(interpolator, work_size)
215 : :
216 : : !Arguments ---------------------------
217 : : !scalars
218 : : class(interpolator_t),intent(inout) :: interpolator
219 : : integer,intent(in) :: work_size
220 : : !*****************************************************************************
221 : :
222 [ - + + - : 6 : ABI_MALLOC(interpolator%btemp,(work_size))
+ - + - -
+ - + -
+ ]
223 [ - + + - : 4 : ABI_MALLOC(interpolator%ctemp,(work_size))
- + - + -
+ ]
224 : :
225 : 2 : end subroutine int_alloc_work
226 : : !!***
227 : :
228 : : !-------------------------------------------------------------------
229 : :
230 : : !!****f* m_bseinterp/int_free
231 : : !! NAME
232 : : !! int_free
233 : : !!
234 : : !! FUNCTION
235 : : !! Deallocate temporary arrays
236 : : !!
237 : : !! SOURCE
238 : :
239 : 2 : subroutine int_free(interpolator)
240 : :
241 : : !Arguments ---------------------------
242 : : class(interpolator_t),intent(inout) :: interpolator
243 : : !*****************************************************************************
244 : :
245 [ + - ]: 2 : ABI_SFREE(interpolator%btemp)
246 [ + - ]: 2 : ABI_SFREE(interpolator%ctemp)
247 : :
248 : 2 : end subroutine int_free
249 : : !!***
250 : :
251 : : !-------------------------------------------------------------------
252 : :
253 : : !!****f* m_bseinterp/int_compute_overlaps
254 : : !! NAME
255 : : !! int_compute_overlaps
256 : : !!
257 : : !! FUNCTION
258 : : !! Compute the overlaps prefactors
259 : : !!
260 : : !! INPUTS
261 : : !!
262 : : !! OUTPUT
263 : : !!
264 : : !! SOURCE
265 : :
266 : 4 : subroutine int_compute_overlaps(interpolator, double_grid, Wfd_dense, Wfd_coarse, &
267 : 4 : Kmesh_dense, Kmesh_coarse, BSp, Cryst, Psps, Pawtab)
268 : :
269 : : !Arguments ---------------------------
270 : : !scalars
271 : : class(interpolator_t),intent(inout) :: interpolator
272 : : type(double_grid_t),intent(in),target :: double_grid
273 : : type(wfdgw_t),intent(inout) :: Wfd_dense, Wfd_coarse
274 : : type(kmesh_t),intent(in) :: Kmesh_dense, Kmesh_coarse
275 : : type(excparam),intent(in) :: BSp
276 : : type(crystal_t),intent(in) :: Cryst
277 : : type(pseudopotential_type),intent(in) :: Psps
278 : : !arrays
279 : : type(pawtab_type),intent(in) :: Pawtab(Cryst%ntypat*Wfd_coarse%usepaw)
280 : :
281 : : !Local variables ---------------------
282 : : !scalars
283 : : integer :: nprocs, my_rank, ierr
284 : : integer :: nfft, nspinor, nsppol, nvert
285 : : integer :: ib_coarse, ib_dense, ik_coarse, ik_dense
286 : : integer :: spin, iorder, ivertex, ix, iy, iz, bstart, bstop
287 : : real(dp),parameter :: threshold = 0.1_dp
288 : : complex(gwp) :: ovlp
289 : : !arrays
290 : : integer :: curindices_dense(6), curindices_coarse(3)
291 : : integer :: neighbour(3)
292 : : integer :: g0(3),g01(3),diffg0(3)
293 : 4 : complex(gwp),allocatable :: ur_coarse(:),ur_dense(:)
294 : 4 : complex(gwp),allocatable :: ceigr(:)
295 : : !*****************************************************************************
296 : :
297 : 4 : nprocs = xmpi_comm_size(Wfd_coarse%comm)
298 : 4 : my_rank = xmpi_comm_rank(Wfd_coarse%comm)
299 : :
300 : : ABI_UNUSED(Pawtab(1)%basis_size)
301 : :
302 : : ! Ensure Wfd and Wfd_coarse use the same FFT mesh.
303 : 4 : call wfd_dense%change_ngfft(Cryst,Psps,Wfd_coarse%ngfft)
304 : 4 : nfft = Wfd_coarse%nfft
305 : 4 : nspinor = Wfd_coarse%nspinor
306 : 4 : nsppol = Bsp%nsppol
307 : 4 : nvert = interpolator%nvert
308 : :
309 : : ! Allocate workspace for wavefunctions in real space.
310 [ - + + - : 12 : ABI_MALLOC(ur_coarse,(nfft*nspinor))
+ - + - -
+ - + ]
311 [ - + + - : 8 : ABI_MALLOC(ur_dense,(nfft*nspinor))
- + - + ]
312 [ - + + - : 8 : ABI_MALLOC(ceigr,(nfft*nspinor))
- + - + ]
313 : :
314 [ + + + + : 117064 : interpolator%overlaps = czero
+ + + + +
+ ]
315 : :
316 : : ! TODO
317 : : ! 1) Choose whether we want to compute only dvv, dcc or all dbb
318 : : ! 2) Check the ordering of the loops
319 : : ! 3) Improve vertex -> neighbour (in double_grid ?)
320 [ + + ]: 8 : do spin = 1,nsppol
321 [ + + ]: 264 : do ik_dense = 1,double_grid%nbz_dense
322 : :
323 : : ! MPI parallelization
324 : : ! We assume that each node owns in memory the full set of wavefunctions
325 : : ! both coarse and dense k-mesh and both spins.
326 [ + - ]: 256 : if (mod(ik_dense, nprocs) /= my_rank) cycle
327 : :
328 : : ! From ik_dense -> indices_dense
329 : 256 : iorder = double_grid%iktoint_dense(ik_dense)
330 [ + + ]: 1024 : g01 = double_grid%g0_dense(:,iorder)
331 [ + + ]: 1792 : curindices_dense = double_grid%indices_dense(:,iorder)
332 : :
333 [ + + ]: 1860 : do ivertex = 1,nvert
334 : :
335 : : ! From vertex to neighbour
336 : : ! TODO improve this part + permit to choose other neighbour (e.g. nearest neighbour for RL)
337 [ + + ]: 1600 : if(nvert > 1) then
338 : 1536 : ix = (ivertex-1)/4
339 : 1536 : iy = (ivertex-ix*4-1)/2
340 : 1536 : iz = (ivertex-ix*4-iy*2-1)
341 : : else
342 : 64 : ix = (BSp%rl_nb-1)/4
343 : 64 : iy = (BSp%rl_nb-ix*4-1)/2
344 : 64 : iz = (BSp%rl_nb-ix*4-iy*2-1)
345 : : end if
346 : :
347 [ + + ]: 6400 : neighbour = [ix,iy,iz]
348 : :
349 : : ! From indices_dense -> indices_coarse
350 [ + + ]: 6400 : curindices_coarse = curindices_dense(1:3) + neighbour(:)
351 : :
352 : : ! From indices_coarse -> ik_ibz in the coarse mesh
353 : : call get_kpt_from_indices_coarse(curindices_coarse,double_grid%maxcomp_coarse,&
354 : 1600 : & double_grid%inttoik_coarse,double_grid%g0_coarse,double_grid%nbz_closedcoarse,ik_coarse,g0)
355 : :
356 : : ! Take into account a possible umklapp between k_dense and k_coarse
357 [ + + ]: 6400 : diffg0 = g0 - g01
358 : :
359 [ + + + + : 4456 : if (ANY(diffg0/=0)) then
+ + ]
360 : : ! WARNING works only with nspinor = 1 !!!
361 : 888 : call calc_ceigr(diffg0,nfft,nspinor,Wfd_coarse%ngfft,ceigr)
362 : : end if
363 : :
364 [ + + ]: 13056 : do ib_dense = BSp%lomo_spin(spin), BSp%humo_spin(spin)
365 : : ! ur(ib_dense, ik_dense)
366 : 11200 : call wfd_dense%sym_ur(Cryst,Kmesh_dense,ib_dense,ik_dense,spin,ur_dense)
367 : :
368 [ + + + + : 31192 : if (ANY(diffg0/=0)) then
+ + ]
369 : : !ur_kbz = ur_kbz*e(ig0r)
370 [ + + ]: 25466952 : ur_dense(:) = ur_dense(:)*ceigr(:)
371 : : end if
372 : :
373 : : ! Uncomment for the complete overlap
374 : : !bstart = BSp%lomo_spin(spin); bstop = BSp%humo_spin(spin)
375 : :
376 : : ! Compute only dvv or dcc
377 [ + + ]: 11200 : if (ib_dense <= BSp%homo_spin(spin)) then
378 : : ! if ib_dense is a valence band => loop on valence bands
379 : 4800 : bstart = BSp%lomo_spin(spin); bstop = BSp%homo_spin(spin)
380 : : else
381 : : ! if ib_dense is a conduction band => loop on conduction bands
382 : 6400 : bstart = BSp%lumo_spin(spin); bstop = BSp%humo_spin(spin)
383 : : end if
384 : :
385 [ + + ]: 52800 : do ib_coarse = bstart, bstop
386 : : ! ur(ib_coarse, ik_coarse)
387 : 40000 : call wfd_coarse%sym_ur(Cryst,Kmesh_coarse,ib_coarse,ik_coarse,spin,ur_coarse)
388 : :
389 : : ! ovlp = < u_{ib_coarse,ik_coarse} | u_{ib_dense,ik_dense} >
390 : 40000 : ovlp = xdotc(nfft,ur_coarse,1,ur_dense,1)/nfft
391 : :
392 : : ! Filter too low values
393 [ + + ]: 40000 : if (ABS(ovlp) < threshold) ovlp = czero
394 : :
395 : 51200 : interpolator%overlaps(ib_coarse,ib_dense,ivertex,ik_dense,spin) = ovlp
396 : : end do ! ib_coarse
397 : :
398 : : !DBYG
399 : : ! write(std_out,*) "nb = ",neighbour
400 : : ! write(std_out,*) "(i1,i2,i3,j1,j2,j3) = ",curindices_dense
401 : : ! write(std_out,*) "ib = ",ib_dense
402 : : ! write(std_out,*) "Sum of dbb = ",REAL(SUM(GWPC_CONJG(interpolator%overlaps(:,ib_dense,ivertex,ik_dense,spin))*interpolator%overlaps(:,ib_dense,ivertex,ik_dense,spin))); call flush(std_out)
403 : : !ENDDBYG
404 : :
405 : : end do ! ib_dense
406 : : end do ! ivertex
407 : : end do ! ik_dense
408 : : end do ! spin
409 : :
410 [ - + ]: 4 : ABI_FREE(ur_coarse)
411 [ - + ]: 4 : ABI_FREE(ur_dense)
412 [ - + ]: 4 : ABI_FREE(ceigr)
413 : :
414 : : ! Gather results on each node.
415 : 4 : call xmpi_sum(interpolator%overlaps,Wfd_coarse%comm,ierr)
416 : :
417 : 4 : end subroutine int_compute_overlaps
418 : : !!***
419 : :
420 : : !----------------------------------------------------------------------
421 : :
422 : : !!****f* m_bseinterp/int_preprocess_tables
423 : : !! NAME
424 : : !! int_preprocess_tables
425 : : !!
426 : : !! FUNCTION
427 : : !! Pre-process tables to improve interpolation technique
428 : : !!
429 : : !! INPUTS
430 : : !!
431 : : !! OUTPUT
432 : : !!
433 : : !! SOURCE
434 : :
435 : 4 : subroutine int_preprocess_tables(interpolator,double_grid)
436 : :
437 : : !Argument ------------------------------------
438 : : !scalars
439 : : class(interpolator_t),intent(inout) :: interpolator
440 : : type(double_grid_t),intent(in) :: double_grid
441 : :
442 : : !Local variables -----------------------------
443 : : !scalars
444 : : integer :: iorder,ik_dense,ik_coarse
445 : : integer :: ix,iy,iz,ineighbour,curdim, curj
446 : : real(dp) :: interp_factor
447 : : !arrays
448 : : integer :: allxyz(3),curindices_dense(6)
449 : 4 : integer,allocatable :: curindex(:)
450 : : !*********************************************
451 : :
452 [ - + + - : 12 : ABI_MALLOC(curindex,(double_grid%nbz_coarse))
+ - + - -
+ - + ]
453 [ + + ]: 36 : curindex = 1
454 : :
455 [ + + + + ]: 236 : interpolator%interp_factors = zero
456 : :
457 [ + + ]: 260 : do ik_dense = 1,double_grid%nbz_dense
458 : :
459 : : ! From ik_ibz in the dense mesh -> indices_dense
460 : 256 : iorder = double_grid%iktoint_dense(ik_dense)
461 : : !g01 = double_grid%g0_dense(:,iorder)
462 : :
463 : : ! From indices_dense -> indices_coarse
464 [ + + ]: 1792 : curindices_dense = double_grid%indices_dense(:,iorder)
465 : :
466 : 256 : ik_coarse = double_grid%dense_to_coarse(ik_dense)
467 : :
468 : : ! Compute multi-linear interpolation factors
469 : : ! Loop over the neighbours
470 [ + + ]: 1856 : do ineighbour = 1,interpolator%nvert
471 : : !TODO helper function from [ix,iy,iz] -> ineighbour and vice versa
472 : 1600 : ix = (ineighbour-1)/4
473 : 1600 : iy = (ineighbour-ix*4-1)/2
474 : 1600 : iz = (ineighbour-ix*4-iy*2-1)
475 [ + + ]: 6400 : allxyz = [ix,iy,iz]
476 : 1600 : interp_factor = one
477 [ + + ]: 6400 : do curdim = 1,3
478 [ + + ]: 4800 : if (interpolator%method == BSE_INTERP_RL) then
479 : : cycle
480 [ - + ]: 4608 : else if(interpolator%method == BSE_INTERP_RL2) then
481 [ # # ]: 0 : if (curdim /= 3) cycle
482 : : end if
483 : 4608 : curj = curindices_dense(3+curdim)
484 : : interp_factor = interp_factor*((allxyz(curdim)*(curj*1.0/double_grid%kmult(curdim)))&
485 : 6400 : & +((1-allxyz(curdim))*(1-(curj*1.0/double_grid%kmult(curdim)))))
486 : : end do
487 : 1856 : interpolator%interp_factors(ineighbour,curindex(ik_coarse)) = interp_factor
488 : : end do
489 : :
490 : 260 : curindex(ik_coarse) = curindex(ik_coarse) + 1
491 : : end do
492 : :
493 [ - + ]: 4 : ABI_FREE(curindex)
494 : :
495 : 4 : end subroutine int_preprocess_tables
496 : : !!***
497 : :
498 : : !-------------------------------------------------------------------
499 : :
500 : : !!****f* m_haydock/int_compute_corresp
501 : : !! NAME
502 : : !! int_compute_corresp
503 : : !!
504 : : !! FUNCTION
505 : : !!
506 : : !! INPUTS
507 : : !! BSp<type(excparam)=The parameter for the Bethe-Salpeter run.
508 : : !! grid <double_grid_t>=Correspondence between coarse and fine k-grid
509 : : !! spin=Spin index.
510 : : !!
511 : : !! OUTPUT
512 : : !! corresp(Bsp%nreh(spin),8)= Correspondence between a transition on the
513 : : !! coarse mesh and its i-th neighbour for i in [1,2,..,8].
514 : : !!
515 : : !! TODO:
516 : : !! Some operations are faster if we allocate with shape (8,nreh(spin))
517 : : !!
518 : : !! SOURCE
519 : :
520 : 4 : subroutine int_compute_corresp(interpolator,BSp,double_grid)
521 : :
522 : : !Arguments ------------------------------------
523 : : class(interpolator_t),intent(inout) :: interpolator
524 : : type(excparam),intent(in) :: BSp
525 : : type(double_grid_t),intent(in) :: double_grid
526 : :
527 : : !Local variables ------------------------------
528 : : !scalars
529 : : integer :: spin
530 : : integer :: itt,ik_dense,ik_coarse,iorder,it_coarse
531 : : integer :: ic,iv,ik_coarse0,it_coarse0,iovlp,ix,iy,iz
532 : : !arrays
533 : : integer :: curindices_dense(6),curindices_coarse(3),g0(3),g01(3),neighbour(3)
534 : : !************************************************************************
535 : :
536 [ + + ]: 8 : do spin=1,interpolator%nsppol
537 [ + + ]: 3080 : do itt=1,BSp%nreh_interp(spin)
538 : : ! From dense itt -> ik_dense, ic, iv
539 : 3072 : ik_dense = BSp%Trans_interp(itt,spin)%k
540 : 3072 : ic = BSp%Trans_interp(itt,spin)%c
541 : 3072 : iv = BSp%Trans_interp(itt,spin)%v
542 : :
543 : : ! From ik_dense -> indices_dense
544 : 3072 : iorder = double_grid%iktoint_dense(ik_dense)
545 : 3072 : g01 = double_grid%g0_dense(:,iorder)
546 : :
547 : : ! Index of the k-point in the coarse mesh.
548 : 3072 : ik_coarse0 = double_grid%dense_to_coarse(ik_dense)
549 : 3072 : it_coarse0 = BSp%vcks2t(iv,ic,ik_coarse0,spin)
550 : :
551 : : ! From indices_dense -> indices_coarse
552 [ + + ]: 21504 : curindices_dense = double_grid%indices_dense(:,iorder)
553 : :
554 : : ! Loop over the 8 neighbors.
555 [ + + ]: 22276 : do iovlp = 1,interpolator%nvert
556 : :
557 : : !TODO : helper function from [ix,iy,iz] -> iovlp and vice versa
558 [ + + ]: 19200 : if(interpolator%nvert > 1) then
559 : 18432 : ix = (iovlp-1)/4
560 : 18432 : iy = (iovlp-ix*4-1)/2
561 : 18432 : iz = (iovlp-ix*4-iy*2-1)
562 : : else
563 : 768 : ix = (BSp%rl_nb-1)/4
564 : 768 : iy = (BSp%rl_nb-ix*4-1)/2
565 : 768 : iz = (BSp%rl_nb-ix*4-iy*2-1)
566 : : end if
567 [ + + ]: 76800 : neighbour = [ix,iy,iz]
568 : :
569 [ + + ]: 76800 : curindices_coarse = curindices_dense(1:3) + neighbour(:)
570 : :
571 : : ! From indices_coarse -> ik_ibz in the coarse mesh
572 : : call get_kpt_from_indices_coarse(curindices_coarse,double_grid%maxcomp_coarse,&
573 : 19200 : & double_grid%inttoik_coarse,double_grid%g0_coarse,double_grid%nbz_closedcoarse,ik_coarse,g0)
574 : :
575 : : ! From ik_coarse, ic, iv to it_coarse
576 : 19200 : it_coarse = BSp%vcks2t(iv,ic,ik_coarse,spin)
577 : :
578 : 22272 : interpolator%corresp(it_coarse0,iovlp,spin) = it_coarse
579 : : end do
580 : : end do ! itt
581 : : end do
582 : :
583 : 4 : end subroutine int_compute_corresp
584 : : !!***
585 : :
586 : : !----------------------------------------------------------------------
587 : :
588 : : !!****f* m_bseinterp/interpolator_normalize
589 : : !! NAME
590 : : !! interpolator_normalize
591 : : !!
592 : : !! FUNCTION
593 : : !! Normalize the overlaps so that \sum_{ib} | d_{kk'}^{b,ib} | ^2 = 1
594 : : !!
595 : : !! INPUTS
596 : : !!
597 : : !! OUTPUT
598 : : !!
599 : : !! SOURCE
600 : :
601 : 4 : subroutine interpolator_normalize(interpolator)
602 : :
603 : : !Arguments ---------------------------
604 : : class(interpolator_t),intent(inout) :: interpolator
605 : :
606 : : !Local variables ---------------------
607 : : !scalars
608 : : integer :: spin, ivertex, ib_dense, ik_dense
609 : : complex(gwp) :: sum_ovlp
610 : : !arrays
611 : 4 : complex(gwp),allocatable :: overlaps(:)
612 : : !*****************************************************************************
613 : :
614 [ - + + - : 12 : ABI_MALLOC(overlaps,(interpolator%mband_coarse))
+ - + - -
+ - + ]
615 [ + + ]: 8 : do spin = 1, interpolator%nsppol
616 [ + + ]: 33 : do ivertex = 1, interpolator%nvert
617 [ + + ]: 229 : do ib_dense = 1, interpolator%mband_dense
618 [ + + ]: 13025 : do ik_dense = 1, interpolator%double_grid%nbz_dense
619 [ + + ]: 115200 : overlaps(:) = interpolator%overlaps(:,ib_dense,ivertex,ik_dense,spin)
620 [ + + ]: 115200 : sum_ovlp = SQRT(REAL(SUM(GWPC_CONJG(overlaps(:))*overlaps(:))))
621 [ + + ]: 12800 : if (ABS(sum_ovlp) > tol6) then
622 [ + + ]: 99936 : overlaps(:) = overlaps(:)/sum_ovlp
623 : : else
624 [ + + ]: 15264 : overlaps(:) = czero
625 : : end if
626 [ + + ]: 115400 : interpolator%overlaps(:,ib_dense,ivertex,ik_dense,spin) = overlaps(:)
627 : : end do
628 : : end do
629 : : end do
630 : : end do
631 [ - + ]: 4 : ABI_FREE(overlaps)
632 : :
633 : 4 : end subroutine interpolator_normalize
634 : : !!***
635 : :
636 : : !-------------------------------------------------------------------
637 : :
638 : : !!****f* m_bseinterp/interpolator_free
639 : : !! NAME
640 : : !! interpolator_free
641 : : !!
642 : : !! FUNCTION
643 : : !! Destroy the interpolator object in memory
644 : : !!
645 : : !! SOURCE
646 : :
647 : 23 : subroutine interpolator_free(interpolator)
648 : :
649 : : !Arguments ---------------------------
650 : : class(interpolator_t),intent(inout) :: interpolator
651 : : !*****************************************************************************
652 : :
653 [ + + ]: 23 : ABI_SFREE(interpolator%overlaps)
654 [ + + ]: 23 : ABI_SFREE(interpolator%corresp)
655 [ + + ]: 23 : ABI_SFREE(interpolator%interp_factors)
656 [ + + ]: 23 : if( associated(interpolator%double_grid) ) nullify(interpolator%double_grid)
657 : :
658 : 23 : end subroutine interpolator_free
659 : : !!***
660 : :
661 : : !-------------------------------------------------------------------
662 : :
663 [ # # # # : 0 : end module m_bseinterp
# # # # #
# # # # #
# # # # #
# # # ]
664 : : !!***
|
