## ABINIT/calc_coh [ Functions ]

[ Top ] [ Functions ]

NAME

calc_coh

FUNCTION

Calculates the partial contribution to the COH part of the COHSEX self-energy for a given q-point.

INPUTS

iqibz=index of the irreducible q-point in the array qibz, point which is
related by a symmetry operation to the point q summed over (see csigme).
This index is also used to treat the integrable coulombian singularity at q=0
ngfft(18)=contain all needed information about 3D FFT for GW wavefuntions,
see ~abinit/doc/variables/vargs.htm#ngfft
nsig_ab=Number of components in the self-energy operator (1 for collinear magnetism)
npwc=number of plane waves in $\tilde epsilon^{-1}$
nspinor=Number of spinorial components.
nfftot=number of points in real space
i_sz=contribution arising from the integrable coulombian singularity at q==0
(see csigme for the method used), note that in case of 3-D systems the factor
4pi in the coulombian potential is included in the definition of i_sz
gvec(3,npwc)=G vectors in reduced coordinates
vc_sqrt(npwc)= square root of the coulombian matrix elements for this q-point
epsm1q_o(npwc,npwc)= contains $\tilde epsilon^{-1}(q,w=0) - \delta_{G Gp}$ for
the particular q-point considered in the sum
wfg2_jk(nsig_ab*nfftot)= Fourier Transform of $\u_{jb k}^*(r) u_{kb k}$
jb,kb=left and righ band indeces definining the left and right states where the
partial contribution to the matrix element of $\Sigma_{COH}$ is evaluated

OUTPUT

sigcohme=partial contribution to the matrix element of $<jb k \sigma|\Sigma_{COH} | kb k \sigma>$
coming from this single q-point

PARENTS

cohsex_me

CHILDREN

wrtout

SOURCE

897 subroutine calc_coh(nspinor,nsig_ab,nfftot,ngfft,npwc,gvec,wfg2_jk,epsm1q_o,vc_sqrt,i_sz,iqibz,same_band,sigcohme)
898
899
900 !This section has been created automatically by the script Abilint (TD).
901 !Do not modify the following lines by hand.
902 #undef ABI_FUNC
903 #define ABI_FUNC 'calc_coh'
904 !End of the abilint section
905
906  implicit none
907
908 !Arguments ------------------------------------
909 !scalars
910  integer,intent(in) :: iqibz,nfftot,npwc,nsig_ab,nspinor
911  real(dp),intent(in) :: i_sz
912  logical,intent(in) :: same_band
913 !arrays
914  integer,intent(in) :: gvec(3,npwc),ngfft(18)
915  complex(gwpc),intent(in) :: epsm1q_o(npwc,npwc),vc_sqrt(npwc)
916  complex(gwpc),intent(in) :: wfg2_jk(nfftot*nsig_ab)
917  complex(gwpc),intent(out) :: sigcohme(nsig_ab)
918
919 !Local variables-------------------------------
920 !scalars
921  integer,save :: enough=0
923 !arrays
924  integer :: g2mg1(3)
925
926 ! *************************************************************************
927
928  DBG_ENTER("COLL")
929
930  ! === Partial contribution to the matrix element of Sigma_c ===
931  ! * For nspinor==2, the closure relation reads:
932  !  $\sum_s \psi_a^*(1)\psi_b(2) = \delta_{ab} \delta(1-2)$
933  !  where a,b are the spinor components. As a consequence, Sigma_{COH} is always
934  !  diagonal in spin-space and only diagonal matrix elements have to be calculated.
935  ! MG  TODO wfg2_jk should be calculated on an augmented FFT box to avoid spurious wrapping of G1-G2.
936  ! MG: One has to make sure G1-G2 is still in the FFT mesh for each G1 and G2 in chi0 (not always true)
937  ! MODULO wraps G1-G2 in the FFT box but the Fourier components are not periodic!
938
939  ! * Treat the case q --> 0 adequately.
940  ! TODO Better treatment of wings, check cutoff in the coulombian interaction.
941  igmin=1; if (iqibz==1) igmin=2
942
943  sigcohme(:)=czero_gw
944
945  do ispinor=1,nspinor
947    outofbox=0
948
949    do igp=igmin,npwc
950      do ig=igmin,npwc
951
952       g2mg1 = gvec(:,igp)-gvec(:,ig)
953       if (ANY(g2mg1(:)>ngfft(1:3)/2) .or. ANY(g2mg1(:)<-(ngfft(1:3)-1)/2)) then
954         outofbox = outofbox+1; CYCLE
955       end if
956
957       ig4x=MODULO(g2mg1(1),ngfft(1))
958       ig4y=MODULO(g2mg1(2),ngfft(2))
959       ig4z=MODULO(g2mg1(3),ngfft(3))
960       ig4= 1+ig4x+ig4y*ngfft(1)+ig4z*ngfft(1)*ngfft(2)
961
962       sigcohme(ispinor) = sigcohme(ispinor) + &
964      end do !ig
965    end do !igp
966
967    if (iqibz==1.and.same_band) then
968      sigcohme(ispinor) = sigcohme(ispinor) + half*wfg2_jk(spad+ig4)*epsm1q_o(1,1)*i_sz
969    end if
970  end do !ispinor
971
972  if (outofbox/=0) then
973    enough=enough+1
974    if (enough<=50) then
975      MSG_WARNING(sjoin(' Number of G1-G2 pairs outside the G-sphere for Wfns:', itoa(outofbox)))
976      if (enough==50) then
977        call wrtout(std_out,' ========== Stop writing Warnings ==========')
978      end if
979    end if
980  end if
981
982  DBG_EXIT("COLL")
983
984 end subroutine calc_coh

## ABINIT/cohsex_me [ Functions ]

[ Top ] [ Functions ]

NAME

cohsex_me

FUNCTION

Calculate diagonal and off-diagonal matrix elements of the SEX or COHSEX self-energy operator.

Copyright (C) 1999-2018 ABINIT group (FB, GMR, VO, LR, RWG, MG, RShaltaf)
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

sigmak_ibz=Index of the k-point in the IBZ.
minbnd, maxbnd= min and Max band index for GW correction (for this k-point)
iomode=Option defining the file format of the SCR file (Fortran, NETCDF)
Er <Epsilonm1_results> (see the definition of this structured datatype)
%mqmem=if 0 use out-of-core method in which a single q-slice of espilon is read inside the loop over k
%nomega_i=Number of imaginary frequencies.
%nomega_r=Number of real frequencies.
%nomega=Total number of frequencies.
Gsph_c<gsphere_t>= info on the G-sphere for Sigma_x
%nsym=number of symmetry operations
%rottb(ng,timrev,nsym)=index of (IS) G where I is the identity or the inversion
operation and G is one of the ng vectors in reciprocal space
%timrev=2 if time-reversal symmetry is used, 1 otherwise
%gvec(3,Sigp%npwc)=integer coordinates of each plane wave in reciprocal space
ikcalc=index in the array Sigp%kptgw2bz of the k-point where GW corrections are calculated
Ltg_k datatype containing information on the little group
Kmesh <kmesh_t>
%nbz=Number of points in the BZ
%nibz=Number of points in IBZ
%kibz(3,nibz)=k-point coordinates, irreducible Brillouin zone
%kbz(3,nbz)=k-point coordinates, full Brillouin zone
%ktab(nbz)= table giving for each k-point in the BZ (kBZ), the corresponding
%ktabi(nbz)= for each k-point in the BZ defines whether inversion has to be considered
%ktabp(nbz)= phase factor associated to tnons
gwc_ngfft(18)=Information about 3D FFT for the oscillator strengths used for the correlation part,
Vcp <vcoul_t datatype> containing information on the cutoff technique
%vc_sqrt(npwc,nqibz)= square-root of the coulombian potential for q-points in the IBZ
Pawtab(Psps%ntypat) <type(pawtab_type)>=paw tabulated starting data
Pawang <type(pawang_type)>=paw angular mesh and related data
Psps <type(pseudopotential_type)>=variables related to pseudopotentials
%usepaw=1 for PAW, 0 for NC pseudopotentials.
Qmesh <kmesh_t> : datatype gathering information of the q-mesh used
%ibz=q points where $\tilde\epsilon^{-1}$ has been computed
%bz(3,nqbz)=coordinates of all q-points in BZ
Sigp <sigparams_t> (see the definition of this structured datatype)
Cryst<crystal_t>=Info on unit cell and symmetries
%natom=number of atoms in unit cell
%ucvol=unit cell volume
%nsym=number of symmetry operations
%typat(natom)=type of each atom
much slower but it requires less memory
QP_BSt<ebands_t>=Datatype gathering info on the QP energies (KS if one shot)
eig(Sigp%nbnds,Kmesh%nibz,%nsppol)=KS or QP energies for k-points, bands and spin
occ(Sigp%nbnds,Kmesh%nibz,%nsppol)=occupation numbers, for each k point in IBZ, each band and spin
Paw_pwff<pawpwff_t>=Form factor used to calculate the onsite mat. elements of a plane wave.
allQP_sym(%nkibz,%nsppol)<esymm_t>=Datatype collecting data on the irreducible representaions of the
little group of kcalc in the KS representation as well as the symmetry of the bdgw_k states.
Sr=sigma_t (see the definition of this structured datatype)

OUTPUT

NOTES

1) The treatment of the divergence of Gygi+Baldereschi (PRB 1986) [[cite:Gigy1986]] is included.
2) The calculation of energy derivative is based on finite elements.
3) On the symmetrization of Sigma matrix elements ***/
If  Sk = k+G0 then  M_G(k, Sq)= e^{-i( Sq+G).t} M_{ S^-1(G}   (k,q)
If -Sk = k+G0 then  M_G(k,-Sq)= e^{-i(-Sq+G).t} M_{-S^-1(G)}^*(k,q)

Notice the absence of G0 in the expression. Moreover, when we sum over the little group, it turns out
that there is a cancellation of the phase factor associated to the non-symmorphic operations due to a
similar term coming from the symmetrization of \epsilon^{-1}. Mind however that the nonsymmorphic phase
has to be considered when epsilon^-1 is reconstructed starting from the q-points in the IBZ.

4) The unitary transformation relating wavefunctions
at symmetric k-points should be taken into account during the symmetrization
of the oscillator matrix elements. In case of G_oW_o and GW_o calculations, however,
it is possible to make an invariant by just including all the degenerate states and
averaging the final results over the degenerate subset.

PARENTS

sigma

CHILDREN

calc_coh,calc_wfwfg,epsm1_symmetrizer,esymm_symmetrize_mels,findqg0
get_bz_item,get_epsm1,get_gftt,gsph_fft_tabs,hermitianize
littlegroup_print,paw_rho_tw_g,paw_symcprj,pawcprj_alloc,pawcprj_copy
pawcprj_free,pawpwij_free,pawpwij_init,rho_tw_g,rotate_fft_mesh
sigma_distribute_bks,timab,wfd_change_ngfft,wfd_get_cprj
wfd_get_many_ur,wfd_get_ur,wrtout,xgemv,xmpi_sum

SOURCE

166 subroutine cohsex_me(sigmak_ibz,ikcalc,nomega_sigc,minbnd,maxbnd,Cryst,QP_BSt,Sigp,Sr,Er,Gsph_c,Vcp,&
167 & Kmesh,Qmesh,Ltg_k,Pawtab,Pawang,Paw_pwff,Psps,Wfd,allQP_sym,gwc_ngfft,iomode,prtvol,sigcme_tmp)
168
169
170 !This section has been created automatically by the script Abilint (TD).
171 !Do not modify the following lines by hand.
172 #undef ABI_FUNC
173 #define ABI_FUNC 'cohsex_me'
174 !End of the abilint section
175
176  implicit none
177
178 !Arguments ------------------------------------
179 !scalars
180  integer,intent(in) :: sigmak_ibz,ikcalc,prtvol,iomode,nomega_sigc,minbnd,maxbnd
181  type(crystal_t),intent(in) :: Cryst
182  type(ebands_t),target,intent(in) :: QP_BSt
183  type(kmesh_t),intent(in) :: Kmesh,Qmesh
184  type(vcoul_t),intent(in) :: Vcp
185  type(Epsilonm1_results),intent(inout) :: Er
186  type(gsphere_t),intent(in) :: Gsph_c
187  type(littlegroup_t),intent(in) :: Ltg_k
188  type(Pseudopotential_type),intent(in) :: Psps
189  type(pawang_type),intent(in) :: pawang
190  type(sigparams_t),target,intent(in) :: Sigp
191  type(sigma_t),intent(in) :: Sr
192  type(wfd_t),target,intent(inout) :: Wfd
193 !arrays
194  integer,intent(in) :: gwc_ngfft(18)
195  complex(dpc),intent(out) :: sigcme_tmp(nomega_sigc,minbnd:maxbnd,minbnd:maxbnd,Wfd%nsppol*Sigp%nsig_ab)
196  type(Pawtab_type),intent(in) :: Pawtab(Psps%ntypat)
197  type(pawpwff_t),intent(in) :: Paw_pwff(Psps%ntypat*Psps%usepaw)
198  type(esymm_t),target,intent(in) :: allQP_sym(Wfd%nkibz,Wfd%nsppol)
199
200 !Local variables ------------------------------
201 !scalars
202  integer,parameter :: tim_fourdp=2,ndat1=1
203  integer :: iab,ib,ib1,ib2,ierr,ig,ii,iik,itim_q,i1,i2,npwc
204  integer :: ik_bz,ik_ibz,io,isym_q,iq_bz,iq_ibz,spin,isym,jb,is_idx
205  integer :: band,band1,band2,idle,rank
206  integer :: jik,jk_bz,jk_ibz,kb,nspinor,nsppol
207  integer :: nomega_tot,nq_summed,ispinor,ibsp,dimcprj_gw
209  integer :: ndegs,wtqm,wtqp,mod10
211  real(dp) :: fact_sp,theta_mu_minus_e0i,tol_empty,gw_gsq
212  complex(dpc) :: ctmp,ph_mkgwt,ph_mkt
213  logical :: iscompatibleFFT,q_is_gamma
214  character(len=500) :: msg
215 !arrays
217  integer,allocatable :: proc_distrb(:,:,:),coh_distrb(:,:,:,:),degtab(:,:,:)
218  integer,allocatable :: igfftcg0(:),gw_gfft(:,:),gw_gbound(:,:),irottb(:,:),ktabr(:,:)
219  integer :: got(Wfd%nproc)
220  real(dp) :: ksum(3),kgw(3),kgw_m_ksum(3),q0(3),tsec(2),qbz(3),spinrot_kbz(4),spinrot_kgw(4)
221  real(dp),ABI_CONTIGUOUS pointer :: qp_ene(:,:,:),qp_occ(:,:,:)
222  complex(gwpc) :: sigcohme(Sigp%nsig_ab)
223  complex(gwpc),allocatable :: vc_sqrt_qbz(:),rhotwg(:),rhotwgp(:),sigsex(:)
224  complex(gwpc),allocatable :: epsm1_qbz(:,:,:)
225  complex(gwpc),allocatable :: sigc_ket(:,:)
226  complex(gwpc),allocatable :: rhotwg_ki(:,:)
227  complex(gwpc),allocatable :: sigctmp(:,:)
228  complex(gwpc),allocatable :: wfr_bdgw(:,:),ur_sum(:),wf1swf2_g(:)
229  complex(gwpc),ABI_CONTIGUOUS pointer :: cg_jb(:),cg_sum(:)
230  complex(dpc),allocatable :: sym_cme(:,:,:,:),sigc(:,:,:,:,:)
233  type(sigijtab_t),pointer :: Sigcij_tab(:)
234  type(pawcprj_type),allocatable :: Cprj_kgw(:,:),Cprj_ksum(:,:)
235  type(pawpwij_t),allocatable :: Pwij_qg(:),Pwij_fft(:)
236  type(esymm_t),pointer :: QP_sym(:)
237
238 !************************************************************************
239
240  DBG_ENTER("COLL")
241
242  call timab(423,1,tsec) ! cohsex_me
243
244  ! Initial check
245  ABI_CHECK(Sr%nomega_r == Sigp%nomegasr,"")
246  ABI_CHECK(Sr%nomega4sd == Sigp%nomegasrd,"")
247  !ABI_CHECK(Sigp%npwc==Gsph_c%ng,"")
248
249  ! Initialize some values
250  nspinor = Wfd%nspinor; nsppol = Wfd%nsppol
251  npwc = sigp%npwc
252  spinor_padc = RESHAPE([0, 0, npwc, npwc, 0, npwc, npwc,0], [2, 4])
253
254  qp_ene => QP_BSt%eig; qp_occ => QP_BSt%occ
255
256  ! Extract the symmetries of the bands for this k-point
257  QP_sym => allQP_sym(sigmak_ibz,1:nsppol)
258
259  ! Index of the GW point in the BZ array, its image in IBZ and time-reversal ===
260  jk_bz=Sigp%kptgw2bz(ikcalc)
261  call get_BZ_item(Kmesh,jk_bz,kgw,jk_ibz,isym_kgw,jik,ph_mkgwt)
262  !%call get_IBZ_item(Kmesh,jk_ibz,kibz,wtk)
263  spinrot_kgw=Cryst%spinrot(:,isym_kgw)
264  ib1 = minbnd; ib2 = maxbnd
265
266  write(msg,'(2a,3f8.3,2a,2(i3,a))')ch10,&
267 &  ' Calculating <nk|Sigma_c(omega)|nk> at k = ',kgw(:),ch10,&
268 &  ' bands n = from ',ib1,' to ',ib2,ch10
269  call wrtout(std_out,msg,'COLL')
270
271  if (ANY(gwc_ngfft(1:3) /= Wfd%ngfft(1:3)) ) call wfd_change_ngfft(Wfd,Cryst,Psps,gwc_ngfft)
272  gwc_mgfft = MAXVAL(gwc_ngfft(1:3))
273  gwc_fftalga = gwc_ngfft(7)/100 !; gwc_fftalgc=MOD(gwc_ngfft(7),10)
274
275  can_symmetrize = .FALSE.
276  if (Sigp%symsigma>0) then
277    can_symmetrize = .TRUE.
278    if (Sigp%gwcalctyp >= 20) then
279     do spin=1,Wfd%nsppol
280       can_symmetrize(spin) = .not. esymm_failed(QP_sym(spin))
281       if (.not.can_symmetrize(spin)) then
282         write(msg,'(a,i0,4a)')" Symmetrization cannot be performed for spin: ",spin,ch10,&
283 &         " band classification encountered the following problem: ",ch10,TRIM(QP_sym(spin)%err_msg)
284         MSG_WARNING(msg)
285       end if
286     end do
287    end if
288    if (nspinor == 2) MSG_WARNING('Symmetrization with nspinor=2 not implemented')
289  end if
290
291  mod10=MOD(Sigp%gwcalctyp, 10)
292
293  call timab(491,1,tsec) ! csigme(tot) Overall clock. TODO check this
294  call timab(495,1,tsec) ! csigme (SigC)
295
296  ! Normalization of theta_mu_minus_e0i
297  ! If nsppol==2, qp_occ $\in [0,1]$
298  SELECT CASE (nsppol)
299  CASE (1)
300    fact_sp=half; tol_empty=0.01   ! below this value the state is assumed empty
301    if (nspinor==2) then
302     fact_sp=one; tol_empty=0.005  ! below this value the state is assumed empty
303    end if
304  CASE (2)
305    fact_sp=one; tol_empty=0.005   ! to be consistent and obtain similar results if a metallic
306  CASE DEFAULT                     ! spin unpolarized system is treated using nsppol==2
307    MSG_BUG('Wrong nsppol')
308  END SELECT
309
310  call timab(442,1,tsec) ! csigme(init0)
311
312  ! Precalculate the FFT index of $(R^{-1}(r-\tau))$ ===
313  ! S=\transpose R^{-1} and k_BZ = S k_IBZ
314  ! irottb is the FFT index of $R^{-1} (r-\tau)$ used to symmetrize u_Sk.
315  gwc_nfftot = PRODUCT(gwc_ngfft(1:3))
316  ABI_MALLOC(irottb,(gwc_nfftot,Cryst%nsym))
317  call rotate_FFT_mesh(Cryst%nsym,Cryst%symrel,Cryst%tnons,gwc_ngfft,irottb,iscompatibleFFT)
318  if (.not.iscompatibleFFT) then
319    MSG_WARNING("FFT mesh is not compatible with symmetries. Results might be affected by large errors!")
320  end if
321
322  ABI_MALLOC(ktabr,(gwc_nfftot, Kmesh%nbz))
323  do ik_bz=1,Kmesh%nbz
324    isym=Kmesh%tabo(ik_bz)
325    do ifft=1,gwc_nfftot
326      ktabr(ifft,ik_bz)=irottb(ifft,isym)
327    end do
328  end do
329  ABI_FREE(irottb)
330
331  ! The number of occupied states for each point in the IBZ and spin.
332  ! nbv_ks(:,:) = COUNT(qp_occ>=tol_empty,DIM=1)  MG: g95 returns random numbers, likely a bug in the compiler
333  do spin=1,nsppol
334    do ik_ibz=1,Kmesh%nibz
335      nbv_ks(ik_ibz,spin) = COUNT(qp_occ(:,ik_ibz,spin)>=tol_empty)
336    end do
337  end do
338
339  ! (b,k,s) mask for MPI distribution of the sum over occupied states in the BZ.
342  do spin=1,nsppol
343    do ik_bz=1,Kmesh%nbz
344       ik_ibz = Kmesh%tab(ik_bz)
346    end do
347  end do
348
349  ! Distribute the individual terms of the sum over the BZ taking into account symmetries and MPI memory distribution.
350  ! got is used to optimize the distribution if more than one node can calculate the same (b,k,s) element.
351  got=0
352  ABI_MALLOC(proc_distrb,(Wfd%mband,Kmesh%nbz,nsppol))
353  call sigma_distribute_bks(Wfd,Kmesh,Ltg_k,Qmesh,nsppol,can_symmetrize,kgw,Sigp%mg0,my_nbks,&
355
357
358  write(msg,'(a,i0,a)')" Will sum ",my_nbks," (b,k,s) occupied states in (COHSEX|SEX)."
359  call wrtout(std_out,msg,'PERS')
360
361  Sigcij_tab => Sigp%Sigcij_tab(ikcalc,1:nsppol)
362
363  if (mod10==SIG_COHSEX) then
364    ! Distribute the COHSEX terms, taking into account the symmetries of the Sigma_ij matrix.
365    ABI_MALLOC(coh_distrb,(ib1:ib2,ib1:ib2,Kmesh%nbz,nsppol))
366
367    coh_distrb = xmpi_undefined_rank
368    do spin=1,nsppol
369      do ik_bz=1,Kmesh%nbz
370         if (ANY(proc_distrb(:,ik_bz,spin) /= xmpi_undefined_rank) ) then ! This BZ point will be calculated.
371            rank_mask = .FALSE. ! To select only those nodes that will treat (k,s).
372            do band=1,Wfd%mband
373              rank = proc_distrb(band,ik_bz,spin)
374              if (rank /= xmpi_undefined_rank) rank_mask(rank+1)=.TRUE.
375            end do
376            do band2=ib1,ib2
377              do irow=1,Sigcij_tab(spin)%col(band2)%size1   ! Looping over the upper triangle of sigma_ij with non-zero elements.
378                band1 = Sigcij_tab(spin)%col(band2)%bidx(irow)
380                got(idle) = got(idle)+1
381                coh_distrb(band1,band2,ik_bz,spin) = idle-1
382              end do
383            end do
384         end if
385      end do
386    end do
387
388    write(msg,'(a,i0,a)')" will treat ",COUNT(coh_distrb==Wfd%my_rank)," COH terms."
389    call wrtout(std_out,msg,'PERS')
390  end if
391
392  ABI_MALLOC(rhotwg_ki, (npwc * nspinor, minbnd:maxbnd))
393  rhotwg_ki=czero_gw
394  ABI_MALLOC(rhotwg, (npwc * nspinor))
395  ABI_MALLOC(rhotwgp  ,(npwc * nspinor))
396  ABI_MALLOC(vc_sqrt_qbz, (npwc))
397
398  ! Additional allocations for PAW
399  if (Psps%usepaw==1) then
400    ABI_DT_MALLOC(Cprj_ksum,(Cryst%natom,nspinor))
401    call pawcprj_alloc(Cprj_ksum,0,Wfd%nlmn_atm)
402
403    ! For COHSEX we need the onsite terms of the PW on the FFT mesh.
404    ! gw_gfft is the set of plane waves in the FFT Box for the oscillators.
405    if (mod10==SIG_COHSEX) then
406      ABI_MALLOC(gw_gfft,(3,gwc_nfftot))
407      q0=zero
408      call get_gftt(gwc_ngfft,q0,Cryst%gmet,gw_gsq,gw_gfft)
409      ABI_DT_MALLOC(Pwij_fft,(Psps%ntypat))
410      call pawpwij_init(Pwij_fft,gwc_nfftot,(/zero,zero,zero/),gw_gfft,Cryst%rprimd,Psps,Pawtab,Paw_pwff)
411    end if
412  end if ! usepaw==1
413
414  ! === Calculate total number of frequencies and allocate related arrays ===
415  ! sigcme2 is used to accumulate the diagonal matrix elements over k-points and
416  ! GW bands, used only in case of ppmodel 3 and 4 (TODO save memory)
417  nomega_tot=Sr%nomega_r+Sr%nomega4sd
418
419  ABI_MALLOC(sigctmp, (nomega_sigc,Sigp%nsig_ab))
420  sigctmp = czero_gw
421  ABI_MALLOC(sigc_ket, (npwc*nspinor, nomega_sigc))
422
423  if (mod10==SIG_COHSEX)  then
424    ABI_MALLOC(wf1swf2_g,(gwc_nfftot*nspinor))
425  end if
426
427  ! Arrays storing the contribution given by the Hermitian/anti-Hermitian part of \Sigma_c
428  !allocate(aherm_sigc_ket(npwc*nspinor,nomega_sigc))
429  !allocate( herm_sigc_ket(npwc*nspinor,nomega_sigc))
430  ABI_MALLOC(sigsex,(npwc))
431  sigcme_tmp=czero
432
433  ABI_MALLOC(sigc,(2,nomega_sigc,ib1:ib2,ib1:ib2,nsppol*Sigp%nsig_ab))
434  sigc=czero
435
436  ! Here we divide the states where the QP energies are required into complexes. Note however that this approach is not
437  ! based on group theory, and it might lead to spurious results in case of accidental degeneracies.
438  nq_summed=Kmesh%nbz
439  if (Sigp%symsigma>0) then
440    call littlegroup_print(Ltg_k,std_out,prtvol,'COLL')
441    nq_summed=SUM(Ltg_k%ibzq(:))
442    !
443    ! Find number of degenerate states and number of bands in each subspace
444    ! The tolerance is a little bit arbitrary (0.001 eV)
445    ! It could be reduced, in particular in case of nearly accidental degeneracies
446    ABI_MALLOC(degtab,(ib1:ib2,ib1:ib2,nsppol))
447    degtab=0
448    do spin=1,nsppol
449      do ib=ib1,ib2
450        do jb=ib1,ib2
451         if (ABS(qp_ene(ib,jk_ibz,spin)-qp_ene(jb,jk_ibz,spin))<0.001/Ha_ev) then
452           degtab(ib,jb,spin)=1
453         end if
454        end do
455      end do
456    end do
457  end if !symsigma
458
459  write(msg,'(2a,i6,a)')ch10,' calculation status ( ',nq_summed,' to be completed):'
460  call wrtout(std_out,msg,'COLL')
461
462  ! TODO if single q (ex molecule) dont allocate epsm1q, avoid waste of memory
463  ABI_STAT_MALLOC(epsm1_qbz, (npwc, npwc, 1), ierr)
464  ABI_CHECK(ierr==0, "out-of-memory in epsm1_qbz")
465  ABI_MALLOC(igfftcg0,(Gsph_c%ng))
466
467  ! Out-of-core solution for epsilon.
468  if (Er%mqmem==0) then
469    MSG_COMMENT('Reading q-slices from file. Slower but less memory.')
470  end if
471
472  call timab(442,2,tsec)
473
474  ! ==========================================
475  ! ==== Fat loop over k_i in the full BZ ====
476  ! ==========================================
477  ABI_MALLOC(ur_sum,(gwc_nfftot*nspinor))
478
479  do spin=1,nsppol
480    if (ALL(proc_distrb(:,:,spin)/=Wfd%my_rank)) CYCLE
481
482    ABI_MALLOC(wfr_bdgw,(gwc_nfftot*nspinor,ib1:ib2))
483    call wfd_get_many_ur(Wfd, [(jb, jb=ib1,ib2)], jk_ibz, spin, wfr_bdgw)
484
485    if (Wfd%usepaw==1) then
486      ! Load cprj for GW states, note the indexing.
487      dimcprj_gw=nspinor*(ib2-ib1+1)
488      ABI_DT_MALLOC(Cprj_kgw,(Cryst%natom,ib1:ib1+dimcprj_gw-1))
489      call pawcprj_alloc(Cprj_kgw,0,Wfd%nlmn_atm)
490      ibsp=ib1
491      do jb=ib1,ib2
492        call wfd_get_cprj(Wfd,jb,jk_ibz,spin,Cryst,Cprj_ksum,sorted=.FALSE.)
493        call paw_symcprj(jk_bz,nspinor,1,Cryst,Kmesh,Pawtab,Pawang,Cprj_ksum)
494        call pawcprj_copy(Cprj_ksum,Cprj_kgw(:,ibsp:ibsp+(nspinor-1)))
495        ibsp=ibsp+nspinor
496      end do
497    end if
498
499    do ik_bz=1,Kmesh%nbz
500      ! Parallelization over k-points and spin
501      ! For the spin there is another check in the inner loop.
502      if (ALL(proc_distrb(:,ik_bz,spin)/=Wfd%my_rank)) CYCLE
503
504      call timab(443,1,tsec) ! csigme (initq)
505
506      ! Find the corresponding irreducible k-point
507      call get_BZ_item(Kmesh,ik_bz,ksum,ik_ibz,isym_ki,iik,ph_mkt)
508      spinrot_kbz(:)=Cryst%spinrot(:,isym_ki)
509
510      ! Identify q and G0 where q+G0=k_GW-k_i
511      kgw_m_ksum=kgw-ksum
512      call findqg0(iq_bz,g0,kgw_m_ksum,Qmesh%nbz,Qmesh%bz,Sigp%mG0)
513
514      ! Symmetrize the matrix elements.
515      ! Sum only q"s in IBZ_k. In this case elements are weighted
516      ! according to wtqp and wtqm. wtqm is for time-reversal.
517      wtqp=1; wtqm=0
518      if (can_symmetrize(spin)) then
519        if (Ltg_k%ibzq(iq_bz)/=1) CYCLE
520        wtqp=0; wtqm=0
521        do isym=1,Ltg_k%nsym_sg
522          wtqp=wtqp+Ltg_k%wtksym(1,isym,iq_bz)
523          wtqm=wtqm+Ltg_k%wtksym(2,isym,iq_bz)
524        end do
525      end if
526
527      !%write(msg,'(2(a,i4),a,i3)')' csigme : ik_bz ',ik_bz,'/',Kmesh%nbz,' done by processor ',Wfd%my_rank
528      !%call wrtout(std_out,msg,'PERS')
529
530      ! Find the corresponding irred q-point.
531      call get_BZ_item(Qmesh,iq_bz,qbz,iq_ibz,isym_q,itim_q)
532      q_is_gamma = (normv(qbz,Cryst%gmet,"G") < GW_TOL_W0)
533
534      ! Tables for the FFT of the oscillators.
535      !  a) FFT index of the G-G0.
536      !  b) gw_gbound table for the zero-padded FFT performed in rhotwg.
537      ABI_MALLOC(gw_gbound,(2*gwc_mgfft+8,2))
539      if ( ANY(gwc_fftalga == [2, 4]) ) use_padfft=0 ! Pad-FFT is not coded in rho_tw_g
541        ABI_FREE(gw_gbound)
543      end if
544
545      if (Psps%usepaw==1) then
546        ! Get PAW oscillator matrix elements $<phj/r|e^{-i(q+G)}|phi/r> - <tphj/r|e^{-i(q+G)}|tphi/r>$ in packed form.
547        ABI_DT_MALLOC(Pwij_qg,(Psps%ntypat))
548        q0 = qbz !;if (q_is_gamma) q0 = (/0.00001_dp,0.00001_dp,0.00001_dp/) ! GW_Q0_DEFAULT
549        call pawpwij_init(Pwij_qg, npwc, q0, Gsph_c%gvec, Cryst%rprimd, Psps, Pawtab, Paw_pwff)
550      end if
551
552      if (Er%mqmem==0) then
553        ! Read q-slice of epsilon^{-1}|chi0 in Er%epsm1(:,:,:,1) (much slower but less memory).
554        call get_epsm1(Er,Vcp,0,0,iomode,xmpi_comm_self,iqibzA=iq_ibz)
555      end if
556
557      ! Only omega==0 for SEX or COHSEX
558      call Epsm1_symmetrizer(iq_bz, 1, npwc, Er, Gsph_c, Qmesh, .True., epsm1_qbz)
559
560      ! Get Fourier components of the Coulomb interaction in the BZ.
561      ! In 3D systems, neglecting umklapp,  vc(Sq,sG)=vc(q,G)=4pi/|q+G|
562      ! The same relation holds for 0-D systems, but not in 1-D or 2D systems. It depends on S.
563      do ig=1,npwc
564        vc_sqrt_qbz(Gsph_c%rottb(ig,itim_q,isym_q)) = Vcp%vc_sqrt(ig,iq_ibz)
565      end do
566
567      call timab(443,2,tsec) ! csigme (initq)
568
569      ! Sum over bands.
570      do ib=1,Sigp%nbnds
571        ! Parallelism over spin
572        ! This processor has this k-point but what about spin?
573        if (proc_distrb(ib,ik_bz,spin)/=Wfd%my_rank) CYCLE
574
575        ! Skip empty state ib for HF, SEX, and COHSEX.
576        if (qp_occ(ib,ik_ibz,spin)<tol_empty) CYCLE
577
578        theta_mu_minus_e0i=fact_sp*qp_occ(ib,ik_ibz,spin)
579
580        call wfd_get_ur(Wfd,ib,ik_ibz,spin,ur_sum)
581
582        if (Psps%usepaw==1) then
583          ! Load cprj for point ksum, this spin or spinor and *THIS* band.
584          ! TODO MG I could avoid doing this but I have to exchange spin and bands ???
585          ! For sure there is a better way to do this!
586          call wfd_get_cprj(Wfd,ib,ik_ibz,spin,Cryst,Cprj_ksum,sorted=.FALSE.)
587          call paw_symcprj(ik_bz,nspinor,1,Cryst,Kmesh,Pawtab,Pawang,Cprj_ksum)
588        end if
589
590        do jb=ib1,ib2
591          ! Get all <k-q,ib,s|e^{-i(q+G).r}|s,jb,k>, at once.
593 &          ur_sum       ,iik,ktabr(:,ik_bz),ph_mkt  ,spinrot_kbz,  &
594 &          wfr_bdgw(:,jb),jik,ktabr(:,jk_bz),ph_mkgwt,spinrot_kgw,&
595 &          nspinor,rhotwg_ki(:,jb))
596
597          if (Psps%usepaw==1) then
598            ! Add on-site contribution, projectors are already in BZ !TODO Recheck this!
599            i2=jb; if (nspinor==2) i2=(2*jb-1)
601            call paw_rho_tw_g(npwc,nspinor,nspinor,Cryst%natom,Cryst%ntypat,Cryst%typat,Cryst%xred,Gsph_c%gvec,&
603          end if
604
605          ! Multiply by the square root of the Coulomb term.
606          ! In 3-D systems, the factor sqrt(4pi) is included)
607          do ii=1,nspinor
608            spad = (ii-1) * npwc
610          end do
611
612          ! === Treat analytically the case q --> 0 ===
613          ! * The oscillator is evaluated at q=O as it is considered constant in the small cube around Gamma
614          !   while the Colulomb term is integrated out
615          ! * In the scalar case we have nonzero contribution only if ib==jb
616          ! * For nspinor==2 evalute <ib,up|jb,up> and <ib,dwn|jb,dwn>,
617          !   impose orthonormalization since npwwfn might be < npwvec.
618          if (ik_bz == jk_bz) then
619            if (nspinor == 1) then
620              rhotwg_ki(1, jb) = czero_gw
621              if (ib==jb) rhotwg_ki(1, jb)=CMPLX(SQRT(Vcp%i_sz),0.0_gwp)
622            else
623              npw_k = Wfd%npwarr(ik_ibz)
624              rhotwg_ki(1, jb) = zero; rhotwg_ki(npwc+1, jb) = zero
625              if (ib == jb) then
626                cg_sum => Wfd%Wave(ib, ik_ibz, spin)%ug
627                cg_jb  => Wfd%Wave(jb, jk_ibz, spin)%ug
628                ctmp = xdotc(npw_k, cg_sum(1:), 1, cg_jb(1:), 1)
629                rhotwg_ki(1,jb) = CMPLX(SQRT(Vcp%i_sz),0.0_gwp) * real(ctmp)
630                ctmp = xdotc(npw_k, cg_sum(npw_k+1:), 1, cg_jb(npw_k+1:), 1)
631                rhotwg_ki(npwc+1,jb) = CMPLX(SQRT(Vcp%i_sz),0.0_gwp) * real(ctmp)
632                ! PAW is missing
633
634                !rhotwg_ki(1,jb)=CMPLX(SQRT(Vcp%i_sz),0.0_gwp) * sqrt(half)
635                !rhotwg_ki(npwc+1,jb) = CMPLX(SQRT(Vcp%i_sz),0.0_gwp) * sqrt(half)
636              end if
637            end if
638          end if
639        end do !jb  Got all matrix elements from minbnd up to maxbnd.
640
641        do kb=ib1,ib2
642          ! Get the ket \Sigma|\phi_{k,kb}> according to the method.
643          rhotwgp(:) = rhotwg_ki(:,kb)
644          sigc_ket = czero_gw
645
646          ! SEX part. TODO add check on theta_mu_minus_e0i
647          do ispinor=1,nspinor
648            spadc = (ispinor-1) * npwc
650
651            sigsex(:)= -theta_mu_minus_e0i*sigsex(:)
652
653            do io=1,nomega_tot ! nomega==1 as SEX is energy independent.
655            end do
656          end do
657
658          ! Loop over the non-zero row elements of this column.
659          ! 1) If gwcalctyp<20 : only diagonal elements since QP==KS.
660          ! 2) If gwcalctyp>=20:
661          !     * Only off-diagonal elements connecting states with same character.
662          !     * Only the upper triangle if HF, SEX, or COHSEX.
663          do irow=1,Sigcij_tab(spin)%col(kb)%size1
664            jb = Sigcij_tab(spin)%col(kb)%bidx(irow)
665            rhotwg = rhotwg_ki(:,jb)
666
667            ! Calculate <\phi_j|\Sigma_c|\phi_k>
668            ! Different freqs according to method (AC or Perturbative), see nomega_sigc.
669            do iab=1,Sigp%nsig_ab
671              do io=1,nomega_sigc
673              end do
674            end do
675
676            ! TODO: save wf1swf2_g to avoid having to recalculate it at each q-point.
677            if (mod10==SIG_COHSEX) then
678              ! Evaluate Static COH. TODO add spinor.
679              if (coh_distrb(jb,kb,ik_bz,spin) == Wfd%my_rank) then
680                ! COH term is done only once for each k-point.
681                ! It does not depend on the index ib summed over.
682                coh_distrb(jb,kb,ik_bz,spin) = xmpi_undefined_rank
683
684 #if 1
685                call calc_wfwfg(ktabr(:,jk_ibz),jik, spinrot_kgw, & ! TODO why jk_ibz?
686 &                gwc_nfftot,nspinor,gwc_ngfft,wfr_bdgw(:,jb),wfr_bdgw(:,kb),wf1swf2_g)
687 #else
688                ABI_CHECK(jik==1,"jik")
689                call calc_wfwfg(ktabr(:,jk_bz),jik, spinrot_kgw, &
690                  gwc_nfftot,nspinor,gwc_ngfft,wfr_bdgw(:,jb),wfr_bdgw(:,kb),wf1swf2_g)
691 #endif
692
693                if (Psps%usepaw==1) then
694                  i1=jb; i2=kb
695                  if (nspinor==2) then
696                    i1=(2*jb-1); i2=(2*kb-1)
697                  end if
699                  call paw_rho_tw_g(gwc_nfftot,Sigp%nsig_ab,nspinor,Cryst%natom,Cryst%ntypat,Cryst%typat,Cryst%xred,&
701                end if
702
703                call calc_coh(nspinor,Sigp%nsig_ab,gwc_nfftot,gwc_ngfft,npwc,Gsph_c%gvec,wf1swf2_g,epsm1_qbz(:,:,1),&
704 &                vc_sqrt_qbz,Vcp%i_sz,iq_ibz,(jb==kb),sigcohme)
705
706                do io=1,nomega_sigc ! Should be 1
707                  sigctmp(io,:) = sigctmp(io,:)+sigcohme(:)
708                end do
709
710              end if
711            end if ! COHSEX
712
713            ! Accumulate and, in case, symmetrize matrix elements of Sigma_c.
714            do iab=1,Sigp%nsig_ab
715              is_idx = spin; if (nspinor==2) is_idx=iab
716
717              sigcme_tmp(:,jb,kb,is_idx)=sigcme_tmp(:,jb,kb,is_idx) + &
718 &              (wtqp+wtqm)*DBLE(sigctmp(:,iab)) + (wtqp-wtqm)*j_gw*AIMAG(sigctmp(:,iab))
719
720              sigc(1,:,jb,kb,is_idx)=sigc(1,:,jb,kb,is_idx) + wtqp*      sigctmp(:,iab)
721              sigc(2,:,jb,kb,is_idx)=sigc(2,:,jb,kb,is_idx) + wtqm*CONJG(sigctmp(:,iab))
722              ! TODO this should be the contribution coming from the anti-hermitian part.
723            end do
724          end do !jb used to calculate matrix elements of $\Sigma$
725
726        end do !kb to calculate matrix elements of $\Sigma$
727      end do !ib
728
729      ! Deallocate k-dependent quantities.
730      ABI_FREE(gw_gbound)
731      if (Psps%usepaw==1) then
732        call pawpwij_free(Pwij_qg)
733        ABI_DT_FREE(Pwij_qg)
734      end if
735    end do !ik_bz
736
737    ABI_FREE(wfr_bdgw)
738    if (Wfd%usepaw==1) then
739      call pawcprj_free(Cprj_kgw)
740      ABI_DT_FREE(Cprj_kgw)
741    end if
742  end do !spin
743
744  ABI_FREE(igfftcg0)
745
746  ! Gather contributions from all the CPUs.
747  call xmpi_sum(sigcme_tmp, wfd%comm, ierr)
748  call xmpi_sum(sigc, wfd%comm, ierr)
749
750  ! Multiply by constants
751  ! For 3D systems sqrt(4pi) is included in vc_sqrt_qbz.
752  sigcme_tmp = sigcme_tmp /(Cryst%ucvol*Kmesh%nbz)
753  sigc       = sigc       /(Cryst%ucvol*Kmesh%nbz)
754
755  ! If we have summed over the IBZ_q now we have to average over degenerate states.
756  ! Presently only diagonal terms are considered
757  ! TODO it does not work if nspinor==2.
758  do spin=1,nsppol
759    if (can_symmetrize(spin)) then
760      ABI_MALLOC(sym_cme, (nomega_tot, ib1:ib2, ib1:ib2, sigp%nsig_ab))
761      sym_cme=czero
762
763      ! Average over degenerate diagonal elements.
764      ! NOTE: frequencies for \Sigma_c(\omega) should be equal to avoid spurious results.
765      ! another good reason to use a strict criterion for the tolerance on eigenvalues.
766      do ib=ib1,ib2
767        ndegs=0
768        do jb=ib1,ib2
769          if (degtab(ib,jb,spin)==1) then
770            if (nspinor == 1) then
771              sym_cme(:, ib, ib, 1) = sym_cme(:, ib, ib, 1) + SUM(sigc(:, :, jb, jb, spin), dim=1)
772            else
773              do ii=1,sigp%nsig_ab
774                sym_cme(:, ib, ib, ii) = sym_cme(:, ib, ib, ii) + SUM(sigc(:, :, jb, jb, ii), dim=1)
775              end do
776            end if
777
778          end if
779          ndegs = ndegs + degtab(ib,jb,spin)
780        end do
781        sym_cme(:,ib,ib,:) = sym_cme(:,ib,ib,:) / ndegs
782      end do
783
784      if (Sigp%gwcalctyp >= 20) then
785        call esymm_symmetrize_mels(QP_sym(spin),ib1,ib2,sigc(:,1,:,:,spin),sym_cme(1,:,:,1))
786      end if
787
788      ! Copy symmetrized values.
789      do ib=ib1,ib2
790        do jb=ib1,ib2
791          if (nspinor == 1) then
792            sigcme_tmp(:,ib,jb,spin) = sym_cme(:,ib,jb,1)
793          else
794            sigcme_tmp(:,ib,jb,:) = sym_cme(:,ib,jb,:)
795          end if
796        end do
797      end do
798      ABI_FREE(sym_cme)
799    end if
800  end do
801
802  ! Reconstruct the full sigma matrix from the upper triangle (only for HF, SEX and COHSEX)
803  if (Sigp%gwcalctyp>=20 .and. sigma_is_herm(Sigp) ) then
804    ABI_CHECK(nspinor==1,"cannot hermitianize non-collinear sigma!")
805    do spin=1,nsppol
806      do io=1,nomega_sigc
807        call hermitianize(sigcme_tmp(io,:,:,spin),"Upper")
808      end do
809    end do
810  end if
811
812  ! ===========================
813  ! ==== Deallocate memory ====
814  ! ===========================
815  if (Psps%usepaw==1) then
816    if (allocated(gw_gfft)) then
817      ABI_FREE(gw_gfft)
818    end if
819    call pawcprj_free(Cprj_ksum)
820    ABI_DT_FREE(Cprj_ksum)
821    if (allocated(Pwij_fft)) then
822      call pawpwij_free(Pwij_fft)
823      ABI_DT_FREE(Pwij_fft)
824    end if
825  end if
826
827  ABI_FREE(ktabr)
828  ABI_FREE(ur_sum)
829  ABI_FREE(rhotwg_ki)
830  ABI_FREE(rhotwg)
831  ABI_FREE(rhotwgp)
832  ABI_FREE(vc_sqrt_qbz)
833  ABI_FREE(sigc_ket)
834  ABI_FREE(epsm1_qbz)
835  ABI_FREE(sigctmp)
836  ABI_FREE(sigc)
837  ABI_FREE(sigsex)
838  ABI_FREE(proc_distrb)
839  if (mod10==SIG_COHSEX) then
840    ABI_FREE(wf1swf2_g)
841    ABI_FREE(coh_distrb)
842  end if
843
844  if (allocated(degtab)) then
845    ABI_FREE(degtab)
846  end if
847
848  call timab(495,2,tsec) ! csigme(SigC)
849  call timab(491,2,tsec)
850  call timab(423,2,tsec) ! cohsex_me
851
852  DBG_EXIT("COLL")
853
854 end subroutine cohsex_me

## ABINIT/m_cohsex [ Modules ]

[ Top ] [ Modules ]

NAME

m_cohsex

FUNCTION

Calculate diagonal and off-diagonal matrix elements of the SEX or COHSEX self-energy operator.

Copyright (C) 1999-2018 ABINIT group (FB, GMR, VO, LR, RWG, MG, RShaltaf)
GNU General Public License, see ~abinit/COPYING
or http://www.gnu.org/copyleft/gpl.txt .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24
25 #include "abi_common.h"
26
27 module m_cohsex
28
29  use defs_basis
30  use defs_datatypes
31  use m_defs_ptgroups
32  use m_gwdefs !,        only : czero_gw, cone_gw, j_gw, sigparams_t, sigma_type_from_key, sigma_is_herm
33  use m_xmpi
34  use m_errors
35  use m_abicore
36
37  use m_time,          only : timab
38  use m_fstrings,      only : sjoin, itoa
39  use m_hide_blas,     only : xdotc, xgemv
40  use m_numeric_tools, only : hermitianize, imin_loc
41  use m_geometry,      only : normv
42  use m_crystal,       only : crystal_t
43  use m_bz_mesh,       only : kmesh_t, get_BZ_item, findqg0, littlegroup_t, littlegroup_print
44  use m_gsphere,       only : gsphere_t, gsph_fft_tabs
45  use m_fft_mesh,      only : get_gftt, rotate_fft_mesh, cigfft
46  use m_vcoul,         only : vcoul_t
47  use m_pawpwij,       only : pawpwff_t, pawpwij_t, pawpwij_init, pawpwij_free, paw_rho_tw_g
48  use m_wfd,           only : wfd_get_ur, wfd_t, wfd_get_cprj, wfd_change_ngfft, wfd_get_many_ur, wfd_sym_ur
49  use m_oscillators,   only : rho_tw_g, calc_wfwfg
50  use m_screening,     only : epsm1_symmetrizer, get_epsm1, epsilonm1_results
51  use m_esymm,         only : esymm_t, esymm_symmetrize_mels, esymm_failed
52  use m_sigma,         only : sigma_t, sigma_distribute_bks
53  use m_pawang,        only : pawang_type
54  use m_pawtab,        only : pawtab_type
55  use m_pawcprj,       only : pawcprj_type, pawcprj_alloc, pawcprj_free, pawcprj_copy, paw_overlap
56  use m_paw_sym,       only : paw_symcprj
57
58  implicit none
59
60  private