TABLE OF CONTENTS

ABINIT/m_bseinterp [ Modules ]

[ Top ] [ Modules ]

NAME

 m_bseinterp

FUNCTION

COPYRIGHT

  Copyright (C) 2014-2024 ABINIT group (M.Giantomassi, Y. Gillet)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

SOURCE

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

m_bseinterp/int_alloc_work [ Functions ]

[ Top ] [ m_bseinterp ] [ Functions ]

NAME

 int_alloc_work

FUNCTION

 Allocate temporary arrays

INPUTS

OUTPUT

SOURCE

216 subroutine int_alloc_work(interpolator, work_size)
217 
218 !Arguments ---------------------------
219 !scalars
220  integer,intent(in) :: work_size
221  type(interpolator_t),intent(inout) :: interpolator
222 
223 !*****************************************************************************
224 
225  ABI_MALLOC(interpolator%btemp,(work_size))
226  ABI_MALLOC(interpolator%ctemp,(work_size))
227 
228 end subroutine int_alloc_work

m_bseinterp/int_compute_overlaps [ Functions ]

[ Top ] [ m_bseinterp ] [ Functions ]

NAME

 int_compute_overlaps

FUNCTION

 Compute the overlaps prefactors

INPUTS

OUTPUT

SOURCE

275 subroutine int_compute_overlaps(interpolator, double_grid, Wfd_dense, Wfd_coarse, &
276 &   Kmesh_dense, Kmesh_coarse, BSp, Cryst, Psps, Pawtab)
277 
278 !Arguments ---------------------------
279 !scalars
280  type(interpolator_t),intent(inout) :: interpolator
281  type(double_grid_t),intent(in),target :: double_grid
282  type(wfdgw_t),intent(inout) :: Wfd_dense, Wfd_coarse
283  type(kmesh_t),intent(in) :: Kmesh_dense, Kmesh_coarse
284  type(excparam),intent(in) :: BSp
285  type(crystal_t),intent(in) :: Cryst
286  type(pseudopotential_type),intent(in) :: Psps
287 !arrays
288  type(pawtab_type),intent(in) :: Pawtab(Cryst%ntypat*Wfd_coarse%usepaw)
289 
290 !Local variables ---------------------
291 !scalars
292  integer :: nprocs, my_rank
293  integer :: ierr
294  integer :: nfft, nspinor, nsppol, nvert
295  integer :: ib_coarse, ib_dense
296  integer :: ik_coarse, ik_dense
297  integer :: spin
298  integer :: iorder
299  integer :: ivertex, ix, iy, iz
300  integer :: bstart, bstop
301  real(dp),parameter :: threshold = 0.1_dp
302  complex(gwpc) :: ovlp
303 !arrays
304  integer :: curindices_dense(6), curindices_coarse(3)
305  integer :: neighbour(3)
306  integer :: g0(3),g01(3),diffg0(3)
307  complex(gwpc),allocatable :: ur_coarse(:),ur_dense(:)
308  complex(gwpc),allocatable :: ceigr(:)
309 !arrays
310 
311 !*****************************************************************************
312 
313  nprocs = xmpi_comm_size(Wfd_coarse%comm)
314  my_rank = xmpi_comm_rank(Wfd_coarse%comm)
315 
316  ABI_UNUSED(Pawtab(1)%basis_size)
317 
318  ! Ensure Wfd and Wfd_coarse use the same FFT mesh.
319  call wfd_dense%change_ngfft(Cryst,Psps,Wfd_coarse%ngfft)
320  nfft = Wfd_coarse%nfft
321  nspinor = Wfd_coarse%nspinor
322  nsppol = Bsp%nsppol
323  nvert = interpolator%nvert
324 
325  ! Allocate workspace for wavefunctions in real space.
326  ABI_MALLOC(ur_coarse,(nfft*nspinor))
327  ABI_MALLOC(ur_dense,(nfft*nspinor))
328  ABI_MALLOC(ceigr,(nfft*nspinor))
329 
330  interpolator%overlaps = czero
331 
332  ! TODO
333  ! 1) Choose whether we want to compute only dvv, dcc or all dbb
334  ! 2) Check the ordering of the loops
335  ! 3) Improve vertex -> neighbour (in double_grid ?)
336  do spin = 1,nsppol
337    do ik_dense = 1,double_grid%nbz_dense
338 
339      ! MPI parallelization
340      ! We assume that each node owns in memory the full set of wavefunctions
341      ! both coarse and dense k-mesh and both spins.
342      if (mod(ik_dense, nprocs) /= my_rank) cycle
343 
344      ! From ik_dense -> indices_dense
345      iorder = double_grid%iktoint_dense(ik_dense)
346      g01 = double_grid%g0_dense(:,iorder)
347      curindices_dense = double_grid%indices_dense(:,iorder)
348 
349      do ivertex = 1,nvert
350 
351        ! From vertex to neighbour
352        ! TODO improve this part + permit to choose other neighbour (e.g. nearest neighbour for RL)
353        if(nvert > 1) then
354          ix = (ivertex-1)/4
355          iy = (ivertex-ix*4-1)/2
356          iz = (ivertex-ix*4-iy*2-1)
357        else
358          ix = (BSp%rl_nb-1)/4
359          iy = (BSp%rl_nb-ix*4-1)/2
360          iz = (BSp%rl_nb-ix*4-iy*2-1)
361        end if
362 
363        neighbour = [ix,iy,iz]
364 
365        ! From indices_dense -> indices_coarse
366        curindices_coarse = curindices_dense(1:3) + neighbour(:)
367 
368        ! From indices_coarse -> ik_ibz in the coarse mesh
369        call get_kpt_from_indices_coarse(curindices_coarse,double_grid%maxcomp_coarse,&
370 &        double_grid%inttoik_coarse,double_grid%g0_coarse,double_grid%nbz_closedcoarse,ik_coarse,g0)
371 
372        ! Take into account a possible umklapp between k_dense and k_coarse
373        diffg0 = g0 - g01
374 
375        if (ANY(diffg0/=0)) then
376          ! WARNING works only with nspinor = 1 !!!
377          call calc_ceigr(diffg0,nfft,nspinor,Wfd_coarse%ngfft,ceigr)
378        end if
379 
380        do ib_dense = BSp%lomo_spin(spin), BSp%humo_spin(spin)
381          ! ur(ib_dense, ik_dense)
382          call wfd_dense%sym_ur(Cryst,Kmesh_dense,ib_dense,ik_dense,spin,ur_dense)
383 
384          if (ANY(diffg0/=0)) then
385            !ur_kbz = ur_kbz*e(ig0r)
386            ur_dense(:) = ur_dense(:)*ceigr(:)
387          end if
388 
389          ! Uncomment for the complete overlap
390          !bstart = BSp%lomo_spin(spin); bstop = BSp%humo_spin(spin)
391 
392          ! Compute only dvv or dcc
393          if (ib_dense <= BSp%homo_spin(spin)) then
394            ! if ib_dense is a valence band => loop on valence bands
395            bstart = BSp%lomo_spin(spin); bstop = BSp%homo_spin(spin)
396          else
397            ! if ib_dense is a conduction band => loop on conduction bands
398            bstart = BSp%lumo_spin(spin); bstop = BSp%humo_spin(spin)
399          end if
400 
401          do ib_coarse = bstart, bstop
402            ! ur(ib_coarse, ik_coarse)
403            call wfd_coarse%sym_ur(Cryst,Kmesh_coarse,ib_coarse,ik_coarse,spin,ur_coarse)
404 
405            ! ovlp = < u_{ib_coarse,ik_coarse} | u_{ib_dense,ik_dense} >
406            ovlp =  xdotc(nfft,ur_coarse,1,ur_dense,1)/nfft
407 
408            ! Filter too low values
409            if (ABS(ovlp) < threshold) ovlp = czero
410 
411            interpolator%overlaps(ib_coarse,ib_dense,ivertex,ik_dense,spin) = ovlp
412          end do ! ib_coarse
413 
414          !DBYG
415          ! write(std_out,*) "nb = ",neighbour
416          ! write(std_out,*) "(i1,i2,i3,j1,j2,j3) = ",curindices_dense
417          ! write(std_out,*) "ib = ",ib_dense
418          ! 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)
419          !ENDDBYG
420 
421        end do ! ib_dense
422      end do ! ivertex
423    end do ! ik_dense
424  end do ! spin
425 
426  ABI_FREE(ur_coarse)
427  ABI_FREE(ur_dense)
428  ABI_FREE(ceigr)
429 
430  ! Gather results on each node.
431  call xmpi_sum(interpolator%overlaps,Wfd_coarse%comm,ierr)
432 
433 end subroutine int_compute_overlaps

m_bseinterp/int_free_work [ Functions ]

[ Top ] [ m_bseinterp ] [ Functions ]

NAME

 int_free_work

FUNCTION

 Deallocate temporary arrays

INPUTS

OUTPUT

SOURCE

246 subroutine int_free_work(interpolator)
247 
248 !Arguments ---------------------------
249 !scalars
250  type(interpolator_t),intent(inout) :: interpolator
251 
252 !*****************************************************************************
253 
254  ABI_SFREE(interpolator%btemp)
255  ABI_SFREE(interpolator%ctemp)
256 
257 end subroutine int_free_work

m_bseinterp/int_preprocess_tables [ Functions ]

[ Top ] [ m_bseinterp ] [ Functions ]

NAME

 int_preprocess_tables

FUNCTION

 Pre-process tables to improve interpolation technique

INPUTS

OUTPUT

SOURCE

451 subroutine int_preprocess_tables(interpolator,double_grid)
452 
453 !Argument ------------------------------------
454 !scalars
455  type(interpolator_t),intent(inout) :: interpolator
456  type(double_grid_t),intent(in) :: double_grid
457 !arrays
458 
459 !Local variables -----------------------------
460 !scalars
461  integer :: iorder,ik_dense,ik_coarse
462  integer :: ix,iy,iz,ineighbour,curdim, curj
463  real(dp) :: interp_factor
464 !arrays
465  integer :: allxyz(3),curindices_dense(6)
466  integer,allocatable :: curindex(:)
467 
468 !*********************************************
469  ABI_MALLOC(curindex,(double_grid%nbz_coarse))
470  curindex = 1
471 
472  interpolator%interp_factors = zero
473 
474  do ik_dense = 1,double_grid%nbz_dense
475 
476    ! From ik_ibz in the dense mesh -> indices_dense
477    iorder = double_grid%iktoint_dense(ik_dense)
478    !g01 = double_grid%g0_dense(:,iorder)
479 
480    ! From indices_dense -> indices_coarse
481    curindices_dense = double_grid%indices_dense(:,iorder)
482 
483    ik_coarse = double_grid%dense_to_coarse(ik_dense)
484 
485    ! Compute multi-linear interpolation factors
486    ! Loop over the neighbours
487    do ineighbour = 1,interpolator%nvert
488      !TODO helper function from [ix,iy,iz] -> ineighbour and vice versa
489      ix = (ineighbour-1)/4
490      iy = (ineighbour-ix*4-1)/2
491      iz = (ineighbour-ix*4-iy*2-1)
492      allxyz = [ix,iy,iz]
493      interp_factor = one
494      do curdim = 1,3
495        if (interpolator%method == BSE_INTERP_RL) then
496          cycle
497        else if(interpolator%method == BSE_INTERP_RL2) then
498          if (curdim /= 3) cycle
499        end if
500        curj = curindices_dense(3+curdim)
501        interp_factor = interp_factor*((allxyz(curdim)*(curj*1.0/double_grid%kmult(curdim)))&
502 &                               +((1-allxyz(curdim))*(1-(curj*1.0/double_grid%kmult(curdim)))))
503      end do
504      interpolator%interp_factors(ineighbour,curindex(ik_coarse)) = interp_factor
505    end do
506 
507    curindex(ik_coarse) = curindex(ik_coarse) + 1
508  end do
509 
510  ABI_FREE(curindex)
511 
512 end subroutine int_preprocess_tables

m_bseinterp/interpolator_free [ Functions ]

[ Top ] [ m_bseinterp ] [ Functions ]

NAME

 interpolator_free

FUNCTION

 Destroy the interpolator object in memory

INPUTS

OUTPUT

SOURCE

674 subroutine interpolator_free(interpolator)
675 
676 !Arguments ---------------------------
677  type(interpolator_t),intent(inout) :: interpolator
678 
679 !*****************************************************************************
680 
681  ABI_SFREE(interpolator%overlaps)
682  ABI_SFREE(interpolator%corresp)
683  ABI_SFREE(interpolator%interp_factors)
684  if( associated(interpolator%double_grid) ) nullify(interpolator%double_grid)
685 
686 end subroutine interpolator_free

m_bseinterp/interpolator_init [ Functions ]

[ Top ] [ m_bseinterp ] [ Functions ]

NAME

 interpolator_init

FUNCTION

 Construct the interpolator object

INPUTS

OUTPUT

SOURCE

123 subroutine interpolator_init(interpolator, double_grid, Wfd_dense, Wfd_coarse, &
124 &    Kmesh_dense, Kmesh_coarse, BSp, Cryst, Psps, Pawtab, method)
125 
126 !Arguments ---------------------------
127 !scalars
128  integer,intent(in) :: method
129  type(interpolator_t),intent(inout) :: interpolator
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  ABI_CHECK(Wfd_coarse%usepaw==0, "PAW not yet supported")
151  ABI_CHECK(BSp%nsppol==1, "nsppol != 1 not yet implemented")
152  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  interpolator%double_grid => double_grid
158  interpolator%mband_dense = Wfd_dense%mband
159  interpolator%mband_coarse = Wfd_coarse%mband
160  interpolator%method = method
161  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    nvert = 1
170  CASE DEFAULT
171    write(msg,'(a,i0)') "Wrong interpolation method: ",method
172    ABI_ERROR(msg)
173  END SELECT
174 
175  interpolator%nvert = nvert
176 
177  mbandc = interpolator%mband_coarse
178  mbandd = interpolator%mband_dense
179  nbzd = double_grid%nbz_dense
180  nsppol = interpolator%nsppol
181  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 &   Kmesh_coarse, BSp, Cryst, Psps, Pawtab)
185 
186  ABI_MALLOC(interpolator%interp_factors,(nvert,double_grid%ndiv))
187 
188  call int_preprocess_tables(interpolator,double_grid)
189 
190  nreh1 = BSp%nreh(1)
191  nreh2 = nreh1; if(BSp%nsppol == 2) nreh2 = BSp%nreh(2)
192  maxnreh = MAX(nreh1,nreh2)
193 
194  ABI_MALLOC(interpolator%corresp,(maxnreh,interpolator%nvert,interpolator%nsppol))
195 
196  call int_compute_corresp(interpolator,BSp,double_grid)
197 
198 end subroutine interpolator_init

m_bseinterp/interpolator_normalize [ Functions ]

[ Top ] [ m_bseinterp ] [ Functions ]

NAME

 interpolator_normalize

FUNCTION

 Normalize the overlaps so that \sum_{ib} | d_{kk'}^{b,ib} | ^2 = 1

INPUTS

OUTPUT

SOURCE

620 subroutine interpolator_normalize(interpolator)
621 
622 !Arguments ---------------------------
623 !scalars
624  type(interpolator_t),intent(inout) :: interpolator
625 
626 !Local variables ---------------------
627 !scalars
628  integer :: spin, ivertex, ib_dense, ik_dense
629  complex(gwpc) :: sum_ovlp
630 !arrays
631  complex(gwpc),allocatable :: overlaps(:)
632 
633 !*****************************************************************************
634 
635  ABI_MALLOC(overlaps,(interpolator%mband_coarse))
636  do spin = 1, interpolator%nsppol
637    do ivertex = 1, interpolator%nvert
638      do ib_dense = 1, interpolator%mband_dense
639        do ik_dense = 1, interpolator%double_grid%nbz_dense
640          overlaps(:) = interpolator%overlaps(:,ib_dense,ivertex,ik_dense,spin)
641          sum_ovlp = SQRT(REAL(SUM(GWPC_CONJG(overlaps(:))*overlaps(:))))
642          if (ABS(sum_ovlp) > tol6) then
643            overlaps(:) = overlaps(:)/sum_ovlp
644          else
645            overlaps(:) = czero
646          end if
647          interpolator%overlaps(:,ib_dense,ivertex,ik_dense,spin) = overlaps(:)
648        end do
649      end do
650    end do
651  end do
652  ABI_FREE(overlaps)
653 
654 end subroutine interpolator_normalize

m_haydock/int_compute_corresp [ Functions ]

[ Top ] [ m_haydock ] [ Functions ]

NAME

 int_compute_corresp

FUNCTION

INPUTS

 BSp<type(excparam)=The parameter for the Bethe-Salpeter run.
 grid <double_grid_t>=Correspondence between coarse and fine k-grid
 spin=Spin index.

OUTPUT

 corresp(Bsp%nreh(spin),8)= Correspondence between a transition on the
   coarse mesh and its i-th neighbour for i in [1,2,..,8].

 TODO:
  Some operations are faster if we allocate with shape (8,nreh(spin))

SOURCE

537 subroutine int_compute_corresp(interpolator,BSp,double_grid)
538 
539 !Arguments ------------------------------------
540 !scalars
541  type(excparam),intent(in) :: BSp
542  type(double_grid_t),intent(in) :: double_grid
543  type(interpolator_t),intent(inout) :: interpolator
544 
545 !Local variables ------------------------------
546 !scalars
547  integer :: spin
548  integer :: itt,ik_dense,ik_coarse,iorder,it_coarse
549  integer :: ic,iv,ik_coarse0,it_coarse0,iovlp,ix,iy,iz
550 !arrays
551  integer :: curindices_dense(6),curindices_coarse(3),g0(3),g01(3),neighbour(3)
552 
553 !************************************************************************
554 
555  do spin=1,interpolator%nsppol
556    do itt=1,BSp%nreh_interp(spin)
557      ! From dense itt -> ik_dense, ic, iv
558      ik_dense = BSp%Trans_interp(itt,spin)%k
559      ic = BSp%Trans_interp(itt,spin)%c
560      iv = BSp%Trans_interp(itt,spin)%v
561 
562      ! From ik_dense -> indices_dense
563      iorder = double_grid%iktoint_dense(ik_dense)
564      g01 = double_grid%g0_dense(:,iorder)
565 
566      ! Index of the k-point in the coarse mesh.
567      ik_coarse0 = double_grid%dense_to_coarse(ik_dense)
568      it_coarse0 = BSp%vcks2t(iv,ic,ik_coarse0,spin)
569 
570      ! From indices_dense -> indices_coarse
571      curindices_dense = double_grid%indices_dense(:,iorder)
572 
573      ! Loop over the 8 neighbors.
574      do iovlp = 1,interpolator%nvert
575 
576        !TODO : helper function from [ix,iy,iz] -> iovlp and vice versa
577        if(interpolator%nvert > 1) then
578          ix = (iovlp-1)/4
579          iy = (iovlp-ix*4-1)/2
580          iz = (iovlp-ix*4-iy*2-1)
581        else
582          ix = (BSp%rl_nb-1)/4
583          iy = (BSp%rl_nb-ix*4-1)/2
584          iz = (BSp%rl_nb-ix*4-iy*2-1)
585        end if
586        neighbour = [ix,iy,iz]
587 
588        curindices_coarse = curindices_dense(1:3) + neighbour(:)
589 
590        ! From indices_coarse -> ik_ibz in the coarse mesh
591        call get_kpt_from_indices_coarse(curindices_coarse,double_grid%maxcomp_coarse,&
592 &        double_grid%inttoik_coarse,double_grid%g0_coarse,double_grid%nbz_closedcoarse,ik_coarse,g0)
593 
594        ! From ik_coarse, ic, iv to it_coarse
595        it_coarse = BSp%vcks2t(iv,ic,ik_coarse,spin)
596 
597        interpolator%corresp(it_coarse0,iovlp,spin) = it_coarse
598      end do
599    end do ! itt
600  end do
601 
602 end subroutine int_compute_corresp

m_haydock/interpolator_t [ Types ]

[ Top ] [ m_haydock ] [ Types ]

NAME

 interpolator_t

FUNCTION

  Store the overlap matrix elements needed for the interpolation of the BSE Hamiltonian

TODO

  Decide if we want to make the number of bands k-dependent.

SOURCE

 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(gwpc),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(dpc),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  end type interpolator_t
 96 
 97  public :: interpolator_init    ! Construct the object
 98  public :: interpolator_free    ! Free memory
 99  public :: interpolator_normalize ! Normalize the overlaps
100  public :: int_alloc_work       ! Alloc temp memory
101  public :: int_free_work        ! Free temp memory