TABLE OF CONTENTS


ABINIT/exc_build_ham [ Functions ]

[ Top ] [ Functions ]

NAME

  exc_build_ham

FUNCTION

  Calculate and write the excitonic Hamiltonian on an external binary file (Fortran file open
  in random mode) for subsequent treatment in the Bethe-Salpeter code.

COPYRIGHT

 Copyright (C) 1992-2009 EXC group (L.Reining, V.Olevano, F.Sottile, S.Albrecht, G.Onida)
 Copyright (C) 2009-2018 ABINIT group (L.Reining, V.Olevano, F.Sottile, S.Albrecht, G.Onida, M.Giantomassi)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt .

INPUTS

  BSp<excparam>=The parameters for the Bethe-Salpeter calculation.
  BS_files<excfiles>=File names internally used in the BS code.
  Cryst<crystal_t>=Info on the crystalline structure.
  Kmesh<kmesh_t>=The list of k-points in the BZ, IBZ and symmetry tables.
  Qmesh<kmesh_t>=The list of q-points for epsilon^{-1} and related symmetry tables.
  ktabr(nfftot_osc,BSp%nkbz)=The FFT index of $(R^{-1}(r-\tau))$ where R is symmetry needed to obtains
    the k-points from the irreducible image.  Used to symmetrize u_Sk where S = \transpose R^{-1}
  Gsph_x<gsphere_t>=Info on the G-sphere used to describe wavefunctions and W (the largest one is actually stored).
  Gsph_c<gsphere_t>=Info on the G-sphere used to describe the correlation part.
  Vcp<vcoul_t>=The Coulomb interaction in reciprocal space. A cutoff can be used
  W<screen_t>=Data type gathering info and data for W.
  nfftot_osc=Total Number of FFT points used for the oscillator matrix elements.
  ngfft_osc(18)=Info on the FFT algorithm used to calculate the oscillator matrix elements.
  Psps<Pseudopotential_type>=Variables related to pseudopotentials
  Pawtab(Psps%ntypat)<pawtab_type>=PAW tabulated starting data.
  Pawang<pawang_type>=PAW angular mesh and related data.
  Paw_pwff(Cryst%ntypat*Wfd%usepaw)<pawpwff_t>=Form factor used to calculate the onsite mat. elements of a plane wave.
  Wfd<wfd_t>=Handler for the wavefunctions.

OUTPUT

  The excitonic Hamiltonian is saved on an external binary file (see below).

PARENTS

      bethe_salpeter

CHILDREN

      cwtime,get_bz_item,gsph_fft_tabs,paw_rho_tw_g,pawcprj_alloc
      pawcprj_free,pawpwij_free,pawpwij_init,rho_tw_g,timab,wfd_change_ngfft
      wfd_get_cprj,wfd_get_ur,wrtout,xmpi_distab,xmpi_sum

SOURCE

 51 #if defined HAVE_CONFIG_H
 52 #include "config.h"
 53 #endif
 54 
 55 #include "abi_common.h"
 56 
 57 subroutine exc_build_ham(BSp,BS_files,Cryst,Kmesh,Qmesh,ktabr,Gsph_x,Gsph_c,Vcp,&
 58 & Wfd,W,Hdr_bse,nfftot_osc,ngfft_osc,Psps,Pawtab,Pawang,Paw_pwff)
 59 
 60  use defs_basis
 61  use defs_datatypes
 62  use defs_abitypes
 63  use m_profiling_abi
 64  use m_bs_defs
 65  use m_xmpi
 66  use m_errors
 67 
 68  use m_gwdefs,       only : czero_gw
 69  use m_crystal,      only : crystal_t
 70  use m_gsphere,      only : gsphere_t
 71  use m_vcoul,        only : vcoul_t
 72  use m_bz_mesh,      only : kmesh_t
 73  use m_pawpwij,     only : pawpwff_t
 74  use m_pawang,       only : pawang_type
 75  use m_pawtab,       only : pawtab_type
 76  use m_pawcprj,      only : pawcprj_type
 77  use m_wfd,          only : wfd_t
 78  use m_screen,       only : screen_t
 79 
 80 !This section has been created automatically by the script Abilint (TD).
 81 !Do not modify the following lines by hand.
 82 #undef ABI_FUNC
 83 #define ABI_FUNC 'exc_build_ham'
 84  use interfaces_14_hidewrite
 85  use interfaces_18_timing
 86  use interfaces_71_bse, except_this_one => exc_build_ham
 87 !End of the abilint section
 88 
 89  implicit none
 90 
 91 !Arguments ------------------------------------
 92 !scalars
 93  integer,intent(in) :: nfftot_osc
 94  type(excparam),intent(in) :: BSp
 95  type(excfiles),intent(in) :: BS_files
 96  type(screen_t),intent(inout) :: W
 97  type(kmesh_t),intent(in) :: Kmesh,Qmesh
 98  type(crystal_t),intent(in) :: Cryst
 99  type(vcoul_t),intent(in) :: Vcp
100  type(gsphere_t),intent(in) :: Gsph_x,Gsph_c
101  type(Pseudopotential_type),intent(in) :: Psps
102  type(Hdr_type),intent(inout) :: Hdr_bse
103  type(pawang_type),intent(in) :: Pawang
104  type(wfd_t),target,intent(inout) :: Wfd
105 !arrays
106  integer,intent(in) :: ngfft_osc(18)
107  integer,intent(in) :: ktabr(nfftot_osc,Kmesh%nbz)
108  type(Pawtab_type),intent(in) :: Pawtab(Psps%ntypat*Wfd%usepaw)
109  type(pawpwff_t),intent(in) :: Paw_pwff(Psps%ntypat*Wfd%usepaw)
110 
111 !Local variables ------------------------------
112 !scalars
113  logical :: do_resonant,do_coupling
114  !character(len=500) :: msg
115 !arrays
116  real(dp) :: tsec(2)
117  complex(gwpc),allocatable :: all_mgq0(:,:,:,:,:)
118 
119 !************************************************************************
120 
121  call timab(670,1,tsec)
122 
123  ABI_CHECK(Wfd%nspinor==1,"nspinor==2 not coded")
124  ABI_CHECK(nfftot_osc==PRODUCT(ngfft_osc(1:3)),"mismatch in FFT size")
125 
126  if (BSp%have_complex_ene) then
127    MSG_ERROR("Complex energies are not supported yet")
128  end if
129 
130  ! Do we have to compute some block?
131  do_resonant = (BS_files%in_hreso == BSE_NOFILE)
132  do_coupling = (BS_files%in_hcoup == BSE_NOFILE)
133 
134  if (BSp%use_coupling == 0) then
135    if (.not.do_resonant) then
136      call wrtout(std_out,"Will skip the calculation of resonant block (will use BSR file)","COLL")
137      goto 100
138    end if
139  else
140    if (.not. do_resonant .and. .not. do_coupling) then
141      call wrtout(std_out,"Will skip the calculation of both resonant and coupling block (will use BSR and BSC files)","COLL")
142      goto 100
143    end if
144  end if
145 
146  ! Compute M_{k,q=0}^{b,b}(G) for all k-points in the IBZ and each pair b, b'
147  ! used for the exchange part and part of the Coulomb term.
148  call wrtout(std_out," Calculating all matrix elements for q=0 to save CPU time","COLL")
149 
150  call wfd_all_mgq0(Wfd,Cryst,Qmesh,Gsph_x,Vcp,Psps,Pawtab,Paw_pwff,&
151 &  Bsp%lomo_spin,Bsp%homo_spin,Bsp%humo_spin,nfftot_osc,ngfft_osc,Bsp%npweps,all_mgq0)
152 
153  ! ========================
154  ! ==== Resonant Block ====
155  ! ========================
156  if (do_resonant) then
157    call timab(672,1,tsec)
158    call exc_build_block(BSp,Cryst,Kmesh,Qmesh,ktabr,Gsph_x,Gsph_c,Vcp,&
159 &    Wfd,W,Hdr_bse,nfftot_osc,ngfft_osc,Psps,Pawtab,Pawang,Paw_pwff,all_mgq0,.TRUE.,BS_files%out_hreso)
160    call timab(672,2,tsec)
161  end if
162 
163  ! ========================
164  ! ==== Coupling Block ====
165  ! ========================
166  if (do_coupling.and.BSp%use_coupling>0) then
167    call timab(673,1,tsec)
168    call exc_build_block(BSp,Cryst,Kmesh,Qmesh,ktabr,Gsph_x,Gsph_c,Vcp,&
169 &    Wfd,W,Hdr_bse,nfftot_osc,ngfft_osc,Psps,Pawtab,Pawang,Paw_pwff,all_mgq0,.FALSE.,BS_files%out_hcoup)
170    call timab(673,2,tsec)
171  end if
172  !
173  ! * Free memory.
174  ABI_FREE(all_mgq0)
175 
176 100 call timab(670,2,tsec)
177 
178 end subroutine exc_build_ham

ABINIT/wfd_all_mgq0 [ Functions ]

[ Top ] [ Functions ]

NAME

  wfd_all_mgq0

FUNCTION

INPUTS

  Wfd<wfd_t>=Handler for the wavefunctions.
  Cryst<crystal_t>=Info on the crystalline structure.
  Qmesh<kmesh_t>=The list of q-points for epsilon^{-1} and related symmetry tables.
  Gsph_x<gsphere_t>=G-sphere with the G-vectors in mgq0.
  Vcp<vcoul_t>=The Coulomb interaction in reciprocal space. A cutoff can be used
  Psps<Pseudopotential_type>=Variables related to pseudopotentials
  Pawtab(Psps%ntypat)<pawtab_type>=PAW tabulated starting data.
  Paw_pwff(Cryst%ntypat*Wfd%usepaw)<pawpwff_t>=Form factor used to calculate the onsite mat. elements of a plane wave.
  lomo_spin(Wfd%nsppol)=Lowest occupied band for each spin
  homo_spin(Wfd%nsppol)=Highest occupied band for each spin
  humo_spin(Wfd%nsppol)=Highest unoccupied band for each spin
  nfftot_osc=Total Number of FFT points used for the oscillator matrix elements.
  ngfft_osc(18)=Info on the FFT algorithm used to calculate the oscillator matrix elements.
  npweps=Number of G-vectors in mgq0.

OUTPUT

   mgq0(npweps,lomo_min:humo_max,lomo_min:humo_max,Wfd%nkibz,Wfd%nsppol)
     Allocated here and filled with the matrix elements on each node.

PARENTS

      exc_build_ham

CHILDREN

      cwtime,get_bz_item,gsph_fft_tabs,paw_rho_tw_g,pawcprj_alloc
      pawcprj_free,pawpwij_free,pawpwij_init,rho_tw_g,timab,wfd_change_ngfft
      wfd_get_cprj,wfd_get_ur,wrtout,xmpi_distab,xmpi_sum

SOURCE

217 subroutine wfd_all_mgq0(Wfd,Cryst,Qmesh,Gsph_x,Vcp,&
218 & Psps,Pawtab,Paw_pwff,lomo_spin,homo_spin,humo_spin,nfftot_osc,ngfft_osc,npweps,mgq0)
219 
220  use defs_basis
221  use m_profiling_abi
222  use m_xmpi
223  use m_errors
224 
225  use defs_datatypes,  only : pseudopotential_type
226  use m_gwdefs,        only : czero_gw
227  use m_time,          only : cwtime
228  use m_crystal,       only : crystal_t
229  use m_gsphere,       only : gsphere_t, gsph_fft_tabs
230  use m_vcoul,         only : vcoul_t
231  use m_bz_mesh,       only : kmesh_t, get_BZ_item
232  use m_pawpwij,      only : pawpwff_t, pawpwij_t, pawpwij_init, pawpwij_free, paw_rho_tw_g
233  use m_pawtab,        only : pawtab_type
234  use m_pawcprj,       only : pawcprj_type, pawcprj_alloc, pawcprj_free
235  use m_wfd,           only : wfd_t, wfd_get_ur, wfd_get_cprj, wfd_change_ngfft, wfd_ihave_ur, wfd_distribute_bbp
236  use m_oscillators,   only : rho_tw_g
237 
238 !This section has been created automatically by the script Abilint (TD).
239 !Do not modify the following lines by hand.
240 #undef ABI_FUNC
241 #define ABI_FUNC 'wfd_all_mgq0'
242  use interfaces_14_hidewrite
243  use interfaces_18_timing
244 !End of the abilint section
245 
246  implicit none
247 
248 !Arguments ------------------------------------
249 !scalars
250  integer,intent(in) :: nfftot_osc,npweps
251  type(kmesh_t),intent(in) :: Qmesh
252  type(crystal_t),intent(in) :: Cryst
253  type(vcoul_t),intent(in) :: Vcp
254  type(gsphere_t),intent(in) :: Gsph_x
255  type(Pseudopotential_type),intent(in) :: Psps
256  type(wfd_t),target,intent(inout) :: Wfd
257 !arrays
258  integer,intent(in) :: lomo_spin(Wfd%nsppol),homo_spin(Wfd%nsppol),humo_spin(Wfd%nsppol)
259  integer,intent(in) :: ngfft_osc(18)
260  complex(gwpc),allocatable,intent(out) :: mgq0(:,:,:,:,:)
261  type(Pawtab_type),intent(in) :: Pawtab(Psps%ntypat)
262  type(pawpwff_t),intent(in) :: Paw_pwff(Psps%ntypat*Wfd%usepaw)
263 
264 !Local variables ------------------------------
265 !scalars
266  integer,parameter :: map2sphere1=1,dim_rtwg1=1,ndat1=1
267  integer :: use_padfft,mgfft_osc,fftalga_osc,ii
268  integer :: ik_ibz,itim_k,isym_k,iq_bz,iq_ibz,isym_q,itim_q,iqbz0
269  integer :: ierr,iv,ic,spin,lomo_min,humo_max !,inv_ipw,ipw
270  real(dp) :: cpu,wall,gflops !q0vol,fcc_const
271  complex(dpc) :: ph_mkt
272  character(len=500) :: msg
273 !arrays
274  integer,allocatable :: igfftg0(:),task_distrib(:,:,:,:)
275  integer,allocatable :: gbound(:,:),id_tab(:)
276  real(dp) :: qbz(3),spinrot_k(4),tsec(2)
277  complex(gwpc),allocatable :: rhotwg1(:)
278  complex(gwpc),target,allocatable :: ur1(:),ur2(:)
279  complex(gwpc),ABI_CONTIGUOUS pointer :: ptr_ur1(:),ptr_ur2(:)
280  type(pawcprj_type),allocatable :: Cp1(:,:),Cp2(:,:)
281  type(pawpwij_t),allocatable :: Pwij_q0(:)
282 
283 !************************************************************************
284 
285  call timab(671,1,tsec)
286 
287  ABI_CHECK(Wfd%nspinor==1,"nspinor==2 not coded")
288  ABI_CHECK(nfftot_osc==PRODUCT(ngfft_osc(1:3)),"mismatch in FFT size")
289 
290  lomo_min = MINVAL(lomo_spin); humo_max = MAXVAL(humo_spin)
291 
292  if ( ANY(ngfft_osc(1:3) /= Wfd%ngfft(1:3)) ) then
293    call wfd_change_ngfft(Wfd,Cryst,Psps,ngfft_osc)
294  end if
295 
296  mgfft_osc   = MAXVAL(ngfft_osc(1:3))
297  fftalga_osc = ngfft_osc(7)/100 !; fftalgc_osc=MOD(ngfft_osc(7),10)
298 
299  ! (temporary) Table used for the wavefunction in the IBZ.
300  ABI_MALLOC(id_tab, (Wfd%nfft))
301  id_tab = (/(ii, ii=1,Wfd%nfft)/)
302 
303  ! Analytic integration of 4pi/q^2 over the volume element:
304  ! $4pi/V \int_V d^3q 1/q^2 =4pi bz_geometric_factor V^(-2/3)$
305  ! i_sz=4*pi*bz_geometry_factor*q0_vol**(-two_thirds) where q0_vol= V_BZ/N_k
306  ! bz_geometry_factor: sphere=7.79, fcc=7.44, sc=6.188, bcc=6.946, wz=5.255
307  ! (see gwa.pdf, appendix A.4)
308 
309  ! If q=0 and C=V then set up rho-twiddle(G=0) to reflect an
310  ! analytic integration of q**-2 over the volume element:
311  ! <q**-2> = 7.44 V**(-2/3)   (for fcc cell)
312 
313  ! q0vol = (8.0*pi**3) / (Cryst%ucvol*Kmesh%nbz)
314  ! fcc_const = SQRT(7.44*q0vol**(-2.0/3.0))
315  ! rtw = (6.0*pi**2/(Cryst%ucvol*Kmesh%nkbz))**(1./3.)
316  ! Average of (q+q')**-2 integration for head of Coulomb matrix
317  ! INTRTW(QL) = (2*pi*rtw + pi*(rtw**2/QL-QL)*LOG((QL+rtw)/(QL-rtw)))
318  ! &              * (Cryst%ucvol*Kmesh%nbz)/(2*pi)**3. * QL*QL
319 
320  if (Wfd%usepaw==1) then
321    ABI_DT_MALLOC(Cp1,(Wfd%natom,Wfd%nspinor))
322    call pawcprj_alloc(Cp1,0,Wfd%nlmn_atm)
323    ABI_DT_MALLOC(Cp2,(Wfd%natom,Wfd%nspinor))
324    call pawcprj_alloc(Cp2,0,Wfd%nlmn_atm)
325  end if
326 
327  ABI_MALLOC(ur1,(nfftot_osc*Wfd%nspinor))
328  ABI_MALLOC(ur2,(nfftot_osc*Wfd%nspinor))
329 
330  ! Identify q==0
331  iqbz0=0
332  do iq_bz=1,Qmesh%nbz
333    if (ALL(ABS(Qmesh%bz(:,iq_bz))<tol3)) iqbz0 = iq_bz
334  end do
335  ABI_CHECK(iqbz0/=0,"q=0 not found in q-point list!")
336 
337  ! * Get iq_ibz, and symmetries from iqbz0.
338  call get_BZ_item(Qmesh,iqbz0,qbz,iq_ibz,isym_q,itim_q)
339 
340  if (Wfd%usepaw==1) then ! Prepare onsite contributions at q==0
341    ABI_DT_MALLOC(Pwij_q0,(Cryst%ntypat))
342    call pawpwij_init(Pwij_q0,npweps,Qmesh%bz(:,iqbz0),Gsph_x%gvec,Cryst%rprimd,Psps,Pawtab,Paw_pwff)
343  end if
344  !
345  ! Tables for the FFT of the oscillators.
346  !  a) FFT index of the G sphere (only vertical transitions, unlike cchi0, no need to shift the sphere).
347  !  b) gbound table for the zero-padded FFT performed in rhotwg.
348  ABI_MALLOC(igfftg0,(Gsph_x%ng))
349  ABI_MALLOC(gbound,(2*mgfft_osc+8,2))
350  call gsph_fft_tabs(Gsph_x,(/0,0,0/),mgfft_osc,ngfft_osc,use_padfft,gbound,igfftg0)
351  if ( ANY(fftalga_osc == (/2,4/)) ) use_padfft=0 ! Pad-FFT is not coded in rho_tw_g
352  if (use_padfft==0) then
353    ABI_FREE(gbound)
354    ABI_MALLOC(gbound,(2*mgfft_osc+8,2*use_padfft))
355  end if
356 
357  ABI_MALLOC(rhotwg1,(npweps))
358 
359  ABI_STAT_MALLOC(mgq0, (npweps,lomo_min:humo_max,lomo_min:humo_max,Wfd%nkibz,Wfd%nsppol), ierr)
360  ABI_CHECK(ierr==0, "out-of-memory in mgq0")
361  mgq0 = czero
362 
363  call cwtime(cpu,wall,gflops,"start")
364 
365  do spin=1,Wfd%nsppol
366    ! Distribute the calculation of the matrix elements.
367    ! processors have the entire set of wavefunctions hence we divide the workload
368    ! without checking if the pair of states is available. Last dimension is fake.
369    ABI_MALLOC(task_distrib,(lomo_spin(spin):humo_spin(spin),lomo_spin(spin):humo_spin(spin),Wfd%nkibz,1))
370    call xmpi_distab(Wfd%nproc,task_distrib)
371 
372    ! loop over the k-points in IBZ
373    do ik_ibz=1,Wfd%nkibz
374      if ( ALL(task_distrib(:,:,ik_ibz,1)/= Wfd%my_rank) ) CYCLE
375 
376      ! Don't need to symmetrize the wavefunctions.
377      itim_k=1; isym_k=1; ph_mkt=cone; spinrot_k=Cryst%spinrot(:,isym_k)
378 
379      do iv=lomo_spin(spin),humo_spin(spin) ! Loop over band V
380        if ( ALL(task_distrib(:,iv,ik_ibz,1)/=Wfd%my_rank) ) CYCLE
381 
382        if (wfd_ihave_ur(Wfd,iv,ik_ibz,spin,how="Stored")) then
383          ptr_ur1 =>  Wfd%Wave(iv,ik_ibz,spin)%ur
384        else
385          call wfd_get_ur(Wfd,iv,ik_ibz,spin,ur1)
386          ptr_ur1 =>  ur1
387        end if
388 
389        if (Wfd%usepaw==1) then
390          call wfd_get_cprj(Wfd,iv,ik_ibz,spin,Cryst,Cp1,sorted=.FALSE.)
391        end if
392 
393        ! Loop over band C
394        do ic=lomo_spin(spin),humo_spin(spin)
395          if ( task_distrib(ic,iv,ik_ibz,1)/=Wfd%my_rank ) CYCLE
396 
397          if (wfd_ihave_ur(Wfd,ic,ik_ibz,spin,how="Stored")) then
398            ptr_ur2 =>  Wfd%Wave(ic,ik_ibz,spin)%ur
399          else
400            call wfd_get_ur(Wfd,ic,ik_ibz,spin,ur2)
401            ptr_ur2 =>  ur2
402          end if
403 
404          if (Wfd%usepaw==1) then
405            call wfd_get_cprj(Wfd,ic,ik_ibz,spin,Cryst,Cp2,sorted=.FALSE.)
406          end if
407 
408          call rho_tw_g(Wfd%nspinor,npweps,nfftot_osc,ndat1,ngfft_osc,map2sphere1,use_padfft,igfftg0,gbound,&
409 &          ptr_ur1,1,id_tab,ph_mkt,spinrot_k,&
410 &          ptr_ur2,1,id_tab,ph_mkt,spinrot_k,&
411 &          dim_rtwg1,rhotwg1)
412 
413          if (Wfd%usepaw==1) then ! Add PAW onsite contribution.
414            call paw_rho_tw_g(npweps,dim_rtwg1,Wfd%nspinor,Cryst%natom,Cryst%ntypat,Cryst%typat,Cryst%xred,Gsph_x%gvec,&
415 &            Cp1,Cp2,Pwij_q0,rhotwg1)
416          end if
417 
418          ! If q=0 treat Exchange and Coulomb-term independently
419          if (iv <= homo_spin(spin) .and. ic <= homo_spin(spin) .or. &
420 &            iv >  homo_spin(spin) .and. ic >  homo_spin(spin)) then
421 
422            if (iv/=ic) then !COULOMB term: C/=V: ignore them
423              rhotwg1(1) = czero_gw
424            else
425              ! If q=0 and C=V then set up rho-twiddle(G=0) to reflect an
426              ! analytic integration of q**-2 over the volume element:
427              ! <q**-2> = 7.44 V**(-2/3)   (for fcc cell)
428              !rhotwg1(1) = fcc_const * qpg(1,iqbz0)
429              rhotwg1(1) = SQRT(GWPC_CMPLX(Vcp%i_sz,zero)) / Vcp%vcqlwl_sqrt(1,1)
430              !if (vcut) rhotwg1(1) = 1.0
431            end if
432 
433          else
434            ! At present this term is set to zero
435            ! EXCHANGE term: limit value.
436            ! Set up rho-twiddle(G=0) using small vector q instead of zero and k.p perturbation theory (see notes)
437            rhotwg1(1) = czero_gw
438          end if
439 
440          mgq0(:,iv,ic,ik_ibz,spin) = rhotwg1(:)
441        end do !ic
442      end do !iv
443    end do !ik_ibz
444 
445    ABI_FREE(task_distrib)
446  end do !spin
447 
448  ! TODO: One can speedup the calculation by computing the upper triangle of the
449  ! matrix in (b,b') space and then take advantage of the symmetry property:
450  !
451  !   M_{k,0}{{bb'}(G)^* = M{k,0}{b'b'}(-G)
452 
453 #if 0
454  !!!! $OMP PARALLEL DO COLLAPSE(3) PRIVATE(inv_ipw)
455  do spin=1,Wfd%nsppol
456    do ik_ibz=1,Wfd%nkibz
457      do iv=lomo_spin(spin),humo_spin(spin)
458        do ic=1,iv-1
459          do ipw=1,npweps
460            inv_ipw = gsph_x%g2mg(ipw)
461            mgq0(inv_ipw,ic,iv,ik_ibz,spin) = mgq0(ipw,iv,ic,ik_ibz,spin)
462          end do
463        end do
464      end do
465    end do
466  end do
467 #endif
468  !
469  ! Gather matrix elements on each node.
470  call xmpi_sum(mgq0,Wfd%comm,ierr)
471 
472  call cwtime(cpu,wall,gflops,"stop")
473  write(msg,'(2(a,f9.6))')"cpu_time = ",cpu,", wall_time = ",wall
474  call wrtout(std_out,msg,"PERS")
475 
476  ABI_FREE(rhotwg1)
477  ABI_FREE(igfftg0)
478  ABI_FREE(gbound)
479  ABI_FREE(ur1)
480  ABI_FREE(ur2)
481  ABI_FREE(id_tab)
482 
483  if (Wfd%usepaw==1) then
484    ! Deallocation for PAW.
485    call pawpwij_free(Pwij_q0)
486    ABI_DT_FREE(Pwij_q0)
487    call pawcprj_free(Cp1)
488    ABI_DT_FREE(Cp1)
489    call pawcprj_free(Cp2)
490    ABI_DT_FREE(Cp2)
491  end if
492 
493  call timab(671,2,tsec)
494 
495 end subroutine wfd_all_mgq0