TABLE OF CONTENTS


ABINIT/m_paw_init [ Modules ]

[ Top ] [ Modules ]

NAME

  m_paw_init

FUNCTION

  This module contains routines related tp PAW calculations initialization.

COPYRIGHT

 Copyright (C) 2018-2018 ABINIT group (FJ, MT)
 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

17 #if defined HAVE_CONFIG_H
18 #include "config.h"
19 #endif
20 
21 #include "abi_common.h"
22 
23 MODULE m_paw_init
24 
25  use defs_basis
26  use m_errors
27  use m_abicore
28  use m_splines
29 
30  use defs_abitypes,  only : dataset_type
31 
32  use m_time,    only : timab
33  use m_pawpsp,  only : pawpsp_nl
34  use m_paw_atom,only : atompaw_shpfun
35  use m_pawang,  only : pawang_type, pawang_init, pawang_free
36  use m_pawrad,  only : pawrad_type, simp_gen, nderiv_gen, poisson
37  use m_pawtab,  only : pawtab_type
38  use m_paw_numeric, only : paw_derfc
39 
40  implicit none
41 
42  private
43 
44 !public procedures.
45  public :: pawinit     ! Initialize some tabulated data for PAW calculations
46  public :: paw_gencond ! Test whether we have to call pawinit to regenerate tabulated data.
47 
48 CONTAINS  !========================================================================================

m_paw_init/paw_gencond [ Functions ]

[ Top ] [ m_paw_init ] [ Functions ]

NAME

   paw_gencond

FUNCTION

   This routine tests whether we have to call pawinit in one of the optdriver
   routines since important values have changed wrt to the previous dataset.
   The function uses an internal array to store of the previous values

   Usage example:

   call paw_gencond(Dtset,gnt_option,"test",call_pawinit)

   if (psp_gencond == 1 .or. call_pawinit) then
       call pawinit(...)
       call paw_gencond(Dtset, gnt_option, "save", call_pawinit)
   end if

  where psp_gencond is the value returned by pspini.

 INPUT
   Dtset<type(dataset_type)>=all input variables for this dataset
   gnt_option=flag activated if pawang%gntselect and pawang%realgnt have to be allocated
              also determine the size of these pointers
   mode= "test" to test if pawinit must be called
         "save" to update the internal variables.
         "reset" to reset the internal variables

OUTPUT

  call_pawinit=True if pawinit must be called. Meaninfull only if mode=="test"

SIDE EFFECTS

  mode=="save" updates the internal variables.
        "reset" reset the internal variables to -1

PARENTS

      bethe_salpeter,gstate,respfn,screening,sigma,wfk_analyze

CHILDREN

SOURCE

682 subroutine paw_gencond(Dtset,gnt_option,mode,call_pawinit)
683 
684 
685 !This section has been created automatically by the script Abilint (TD).
686 !Do not modify the following lines by hand.
687 #undef ABI_FUNC
688 #define ABI_FUNC 'paw_gencond'
689 !End of the abilint section
690 
691  implicit none
692 
693 !Arguments ------------------------------------
694  integer,intent(in) :: gnt_option
695  logical,intent(out) :: call_pawinit
696  character(len=*),intent(in) :: mode
697  type(dataset_type),intent(in) :: Dtset
698 
699 !Local variables-------------------------------
700 !scalars
701  integer,save :: gencond(9)=(/-1,-1,-1,-1,-1,-1,-1,-1,-1/)
702 
703 ! *********************************************************************
704 
705  call_pawinit = .False.
706  select case (mode)
707  case ("test")
708 
709    if (gencond(1)/=Dtset%pawlcutd .or.gencond(2)/=Dtset%pawlmix  .or.&
710 &   gencond(3)/=Dtset%pawnphi  .or.gencond(4)/=Dtset%pawntheta.or.&
711 &   gencond(5)/=Dtset%pawspnorb.or.gencond(6)/=Dtset%pawxcdev.or.&
712 &   gencond(7)/=Dtset%nsym     .or.gencond(8)/=gnt_option.or.&
713 &   gencond(9)/=Dtset%usepotzero) call_pawinit = .True.
714 
715  case ("save")
716     ! Update internal values
717    gencond(1)=Dtset%pawlcutd ; gencond(2)=Dtset%pawlmix
718    gencond(3)=Dtset%pawnphi  ; gencond(4)=Dtset%pawntheta
719    gencond(5)=Dtset%pawspnorb; gencond(6)=Dtset%pawxcdev
720    gencond(7)=Dtset%nsym     ; gencond(8)=gnt_option
721    gencond(9)=Dtset%usepotzero
722 
723  case ("reset")
724    gencond = -1
725 
726  case default
727    MSG_BUG("Wrong value for mode: "//trim(mode))
728  end select
729 
730 end subroutine paw_gencond

m_paw_init/pawinit [ Functions ]

[ Top ] [ m_paw_init ] [ Functions ]

NAME

 pawinit

FUNCTION

 Initialize some starting values of several arrays used in PAW calculations.

 1-Initialize data related to angular mesh
 2-Tabulate normalized shape function g(r)
 3-Compute indklmn indexes giving some l,m,n,lp,mp,np info
                           from klmn=[(l,m,n),(lp,mp,np)]
 4-Compute various factors/sizes (depending on (l,m,n))
 5-Compute $q_ijL=\displaystyle
                  \int_{0}^{r_c}{(\phi_i\phi_j-\widetilde{\phi_i}\widetilde{\phi_j}) r^l\,dr}
                   Gaunt(l_i m_i,l_j m_j,l m))$
           $S_ij=\displaystyle \sqrt{4 \pi} q_ij0$
 6-Compute $e_ijkl= vh1_ijkl - Vhatijkl - Bijkl - Cijkl$
     With:
       $vh1_ijkl =\sum_{L,m} {vh1*Gaunt(i,j,Lm)*Gaunt(k,l,Lm)}$
       $Vhat_ijkl=\sum_{L,m} {vhatijL*Gaunt(i,j,Lm)*q_klL}$
       $B_ijkl   =\sum_{L,m} {vhatijL*Gaunt(k,l,Lm)*q_ijL}$
       $C_ijkl   =\sum_{L,m} {intvhatL*q_ijL*q_klL}$
     and:
       vh1 according to eq. (A17) in Holzwarth et al., PRB 55, 2005 (1997) [[cite:Holzwarth1997]]
 7-Compute Ex-correlation energy for the core density

COPYRIGHT

 Copyright (C) 1998-2018 ABINIT group (FJ, MT)
 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

  gnt_option=flag activated if pawang%gntselect and pawang%realgnt have to be allocated
             also determine the size of these pointers
  gsqcut_shp=effective cut-off to determine shape functions in reciprocal space
  hyb_range_fock=range coefficient for screened hybrid XC functionals
  lcutdens=max. l for densities/potentials moments computations
  lmix=max. l for which spherical terms will be mixed durinf SCF cycle
  mpsang=1+maximum angular momentum
  nphi="phi" dimension of paw angular mesh
  nsym=Number of symmetry elements in space group
  ntheta="theta" dimension of paw angular mesh
  pawang <type(pawang_type)>=paw angular mesh and related data
     %lmax=Maximum value of angular momentum l+1
     %gntselect((2*l_max-1)**2,l_max**2,l_max**2)=
                     selection rules for Gaunt coefficients
  pawrad(ntypat) <type(pawrad_type)>=paw radial mesh and related data:
     %mesh_size=Dimension of radial mesh
     %rad(mesh_size)=The coordinates of all the points of the radial mesh
     %radfact(mesh_size)=Factor used to compute radial integrals
  pawspnorb=flag: 1 if spin-orbit coupling is activated
  pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data:
     %basis_size=Number of elements for the PAW nl basis
     %l_size=Maximum value of l+1 leading to non zero Gaunt coeffs
     %lmn_size=Number of (l,m,n) elements for the PAW basis
     %lmn2_size=lmn_size*(lmn_size+1)/2
     %dltij(lmn2_size)=factors used to compute sums over (ilmn,jlmn)
     %phi(mesh_size,basis_size)=PAW all electron wavefunctions
     %rshp=shape function radius (radius for compensation charge)
     %shape_type=Radial shape function type
     %shape_alpha=Alpha parameters in Bessel shape function
     %shape_lambda=Lambda parameter in gaussian shape function
     %shape_q=Q parameters in Bessel shape function
     %shape_sigma=Sigma parameter in gaussian shape function
     %tphi(mesh_size,basis_size)=PAW atomic pseudowavefunctions
  pawxcdev=Choice of XC development (0=no dev. (use of angular mesh) ; 1=dev. on moments)
  xclevel=XC functional level (1=LDA, 2=GGA)

OUTPUT

  pawang
     %gntselect(l_size_max**2,l_max**2*(l_max**2+1)/2)=selection rules for Gaunt coefficients
     %l_max=maximum value of angular momentum l+1
     %l_size_max=maximum value of angular momentum l_size=2*l_max-1
     %nsym=number of symmetry elements in space group
     %ngnt=number of non-zero Gaunt coefficients
     %realgnt(pawang%ngnt)=non-zero real Gaunt coefficients
     === only if pawxcdev==1 ==
       %anginit(3,angl_size)=for each point of the angular mesh, gives the coordinates
                             of the corresponding point on an unitary sphere
       %angl_size=dimension of paw angular mesh (angl_size=ntheta*nphi)
       %angwgth(angl_size)=for each point of the angular mesh, gives the weight
                           of the corresponding point on an unitary sphere
       %ntheta, nphi=dimensions of paw angular mesh
       %ylmr(l_size_max**2,angl_size)=real Ylm calculated in real space
       %ylmrgr(1:3,l_size_max**2,angl_size)=first gradients of real Ylm calculated in real space
     === only if pawspnorb==1 ==
     %ls_ylm(2,l_max**2,l_max**2,4)=LS operator in the real spherical harmonics basis
     %use_ls_ylm=flag activated if ls_ylm is allocated
     %usespnorb=flag activated for spin-orbit coupling
  pawtab(ntypat) <type(pawtab_type)>=paw tabulated data read at start:
     %lcut_size_=max. value of l+1 leading to non zero Gaunt coeffs modified by lcutdens
     %lmnmix_sz=number of (lmn,lmn_prime) verifying l<=lmix and l_prime<=lmix
     %mqgrid_shp=number of points in reciprocal space for shape function
     %indklmn(8,lmn2_size)=array giving klm, kln, abs(il-jl) and (il+jl), ilmn and jlmn for each klmn=(ilmn,jlmn)
     %dshpfunc(mesh_size,l_size,4)=derivatives of shape function (used only for numerical shape functions)
     %eijkl(lmn2_size,lmn2_size)=part of the Dij that depends only from the projected occupation coeffs
     %exccore=Exchange-correlation energy for the core density
     %gnorm(l_size)=normalization factor of radial shape function
     %phiphj(:,:)=useful product Phi(:,i)*Phi(:,j)
     %qgrid_shp(mqgrid_shp)=points in reciprocal space for shape function
     %qijl(l_size**2,lmn2_size)=moments of the difference charge density between AE and PS partial wave
     %rad_for_spline(mesh_size)=radial grid used for spline (copy of pawrad%rad)
     %shapefunc(mesh_size,l_size)=normalized radial shape function
     %shapefncg(mqgrid_shp,l_size)=normalized radial shape function in reciprocal space
     %sij(lmn2_size)=nonlocal part of the overlap operator
     %tphitphj(:,:)=useful product tPhi(:,i)*tPhi(:,j)

PARENTS

      bethe_salpeter,gstate,respfn,screening,sigma,wfk_analyze

CHILDREN

SOURCE

169 subroutine pawinit(gnt_option,gsqcut_eff,hyb_range_fock,lcutdens,lmix,mpsang,nphi,nsym,ntheta,&
170 &                  pawang,pawrad,pawspnorb,pawtab,pawxcdev,xclevel,usepotzero)
171 
172 
173 !This section has been created automatically by the script Abilint (TD).
174 !Do not modify the following lines by hand.
175 #undef ABI_FUNC
176 #define ABI_FUNC 'pawinit'
177 !End of the abilint section
178 
179  implicit none
180 
181 !Arguments ---------------------------------------------
182 !scalars
183  integer,intent(in) :: gnt_option,lcutdens,lmix,mpsang,nphi,nsym,ntheta
184  integer,intent(in) :: pawspnorb,pawxcdev,xclevel,usepotzero
185  real(dp),intent(in) :: gsqcut_eff,hyb_range_fock
186  type(pawang_type),intent(inout) :: pawang
187 !arrays
188  type(pawrad_type),intent(in) :: pawrad(:)
189  type(pawtab_type),target,intent(inout) :: pawtab(:)
190 
191 !Local variables ------------------------------
192 !scalars
193  integer,parameter :: mqgrid_shp_default=300
194  integer :: basis_size,i0lm,i0ln,ij_size,il,ilm,ilmn,iln,iloop,iq,isel,isel1
195  integer :: itypat,j0lm,j0lmn,j0ln,jl,jlm,jlmn,jln,klm,klm1
196  integer :: klmn,klmn1,kln,kln1,l_size,ll,lm0,lmax,lmax1,lmin,lmin1,lmn2_size
197  integer :: lmn_size,lmnmix,mesh_size,meshsz,mm,ntypat,usexcnhat,use_ls_ylm,use_ylm
198  real(dp) :: dq,gnrm,intg,ql,ql1,rg,rg1,vh1,yp1,ypn
199  character(len=500) :: message
200 !arrays
201  integer,allocatable :: indl(:,:),klm_diag(:),kmix_tmp(:)
202  integer, ABI_CONTIGUOUS pointer :: indlmn(:,:)
203  real(dp) :: tsec(2)
204  real(dp),allocatable :: ff(:),gg(:),hh(:),indklmn_(:,:),intvhatl(:)
205  real(dp),allocatable :: rad(:),rgl(:,:),vhatijl(:,:),vhatl(:),work(:)
206  real(dp),pointer :: eijkl(:,:)
207 
208 !************************************************************************
209 
210  DBG_ENTER("COLL")
211 
212  call timab(553,1,tsec)
213 
214  ntypat=size(pawtab)
215  if (size(pawrad)/=ntypat) then
216    MSG_BUG('pawrad and pawtab should have the same size!')
217  end if
218 
219  ! Immediately set the value of usepotzero
220  ! it will be used later on in this subroutine
221  pawtab%usepotzero=usepotzero
222 
223 !==================================================
224 !1- INITIALIZE DATA RELATED TO ANGULAR MESH
225 !* ANGULAR GRID
226 !* REAL SPHERICAL HARMONICS
227 !* REAL GAUNT COEFFICIENTS
228 
229  use_ylm=0;if (pawxcdev==0) use_ylm=1
230  use_ls_ylm=0;if (pawspnorb>0) use_ls_ylm=1
231  call pawang_free(pawang)
232  call pawang_init(pawang,gnt_option,mpsang-1,nphi,nsym,ntheta,pawxcdev,use_ls_ylm,use_ylm,xclevel)
233 
234  usexcnhat=maxval(pawtab(1:ntypat)%usexcnhat)
235 
236 !*******************
237 !Loop on atom types
238 !*******************
239  do itypat=1,ntypat
240    mesh_size=pawtab(itypat)%mesh_size
241    l_size=pawtab(itypat)%l_size
242    lmn_size=pawtab(itypat)%lmn_size
243    lmn2_size=pawtab(itypat)%lmn2_size
244    basis_size=pawtab(itypat)%basis_size
245    ij_size=pawtab(itypat)%ij_size
246    indlmn => pawtab(itypat)%indlmn(:,:)
247    ABI_ALLOCATE(indklmn_,(8,lmn2_size))
248    ABI_ALLOCATE(klm_diag,(lmn2_size))
249    ABI_ALLOCATE(ff,(mesh_size))
250    ABI_ALLOCATE(gg,(mesh_size))
251    ABI_ALLOCATE(hh,(mesh_size))
252    ABI_ALLOCATE(rad,(mesh_size))
253    rad(1:mesh_size)=pawrad(itypat)%rad(1:mesh_size)
254 
255    if (pawtab(itypat)%usexcnhat/=usexcnhat) then
256      write(message, '(7a)' )&
257 &     'You cannot simultaneously use atomic data with different',ch10,&
258 &     'formulation of XC [using compensation charge in XC or not] !',ch10,&
259 &     'Action: change at least one of your atomic data (psp) file',ch10,&
260 &     '        or use usexcnhat keyword in input file.'
261      MSG_ERROR(message)
262    end if
263 
264 !  ==================================================
265 !  2- TABULATE SHAPE FUNCTION
266 
267 !  Allocated shape function
268    if (pawtab(itypat)%shape_type/=-1) then
269      if (allocated(pawtab(itypat)%shapefunc))  then
270        ABI_DEALLOCATE(pawtab(itypat)%shapefunc)
271      end if
272      ABI_ALLOCATE(pawtab(itypat)%shapefunc,(mesh_size,l_size))
273    else if (.not.allocated(pawtab(itypat)%shapefunc))  then
274      message='shapefunc should be allocated with shape_type=-1'
275      MSG_ERROR(message)
276    end if
277    if (allocated(pawtab(itypat)%gnorm))  then
278      ABI_DEALLOCATE(pawtab(itypat)%gnorm)
279    end if
280    ABI_ALLOCATE(pawtab(itypat)%gnorm,(l_size))
281 
282 !  Compute shape function
283    do il=1,l_size
284      ll=il-1
285      call atompaw_shpfun(ll,pawrad(itypat),gnrm,pawtab(itypat),ff)
286      pawtab(itypat)%shapefunc(1:mesh_size,il)=ff(1:mesh_size)
287      pawtab(itypat)%gnorm(il)=gnrm
288    end do
289 !  In case of numerical shape function, compute some derivatives
290    if (pawtab(itypat)%shape_type==-1) then
291      if (allocated(pawtab(itypat)%dshpfunc))  then
292        ABI_DEALLOCATE(pawtab(itypat)%dshpfunc)
293      end if
294      ABI_ALLOCATE(pawtab(itypat)%dshpfunc,(mesh_size,l_size,4))
295      ABI_ALLOCATE(work,(mesh_size))
296      do il=1,l_size
297        call nderiv_gen(pawtab(itypat)%dshpfunc(:,il,1),pawtab(itypat)%shapefunc(:,il),pawrad(itypat))
298        yp1=pawtab(itypat)%dshpfunc(1,il,1);ypn=pawtab(itypat)%dshpfunc(mesh_size,il,1)
299        call spline(rad,pawtab(itypat)%shapefunc(:,il),mesh_size,yp1,ypn,pawtab(itypat)%dshpfunc(:,il,2))
300        yp1=pawtab(itypat)%dshpfunc(1,il,2);ypn=pawtab(itypat)%dshpfunc(mesh_size,il,2)
301        call spline(rad,pawtab(itypat)%dshpfunc(:,il,1),mesh_size,yp1,ypn,pawtab(itypat)%dshpfunc(:,il,3))
302        yp1=pawtab(itypat)%dshpfunc(1,il,3);ypn=pawtab(itypat)%dshpfunc(mesh_size,il,3)
303        call spline(rad,pawtab(itypat)%dshpfunc(:,il,2),mesh_size,yp1,ypn,pawtab(itypat)%dshpfunc(:,il,4))
304      end do
305      ABI_DEALLOCATE(work)
306    end if
307 
308 !  In some cases, has to store radial mesh for shape function in pawtab variable
309    if (pawtab(itypat)%shape_type==-1) then
310      if (allocated(pawtab(itypat)%rad_for_spline))  then
311        ABI_DEALLOCATE(pawtab(itypat)%rad_for_spline)
312      end if
313      ABI_ALLOCATE(pawtab(itypat)%rad_for_spline,(mesh_size))
314      pawtab(itypat)%rad_for_spline(1:mesh_size)=pawrad(itypat)%rad(1:mesh_size)
315    end if
316 
317 !  In some cases, has to store shape function in reciprocal space
318    if (pawtab(itypat)%has_shapefncg>0) then
319      if (gsqcut_eff<tol8) then
320        message='Computation of shapefncg only possible when gsqcut>0!'
321        MSG_BUG(message)
322      end if
323      pawtab(itypat)%mqgrid_shp=mqgrid_shp_default
324      if (allocated(pawtab(itypat)%shapefncg))  then
325        ABI_DEALLOCATE(pawtab(itypat)%shapefncg)
326      end if
327      if (allocated(pawtab(itypat)%qgrid_shp))  then
328        ABI_DEALLOCATE(pawtab(itypat)%qgrid_shp)
329      end if
330      ABI_ALLOCATE(pawtab(itypat)%shapefncg,(pawtab(itypat)%mqgrid_shp,2,l_size))
331      ABI_ALLOCATE(pawtab(itypat)%qgrid_shp,(pawtab(itypat)%mqgrid_shp))
332      dq=1.1_dp*sqrt(gsqcut_eff)/dble(pawtab(itypat)%mqgrid_shp-1)
333      do iq=1,pawtab(itypat)%mqgrid_shp
334        pawtab(itypat)%qgrid_shp(iq)=dble(iq-1)*dq
335      end do
336      ABI_ALLOCATE(indl,(6,l_size))
337      ABI_ALLOCATE(rgl,(mesh_size,il))
338      do il=1,l_size
339        indl(:,il)=0;indl(1,il)=il-1;indl(5,il)=il
340        rgl(1:mesh_size,il)=rad(1:mesh_size)*pawtab(itypat)%shapefunc(1:mesh_size,il)
341      end do
342      call pawpsp_nl(pawtab(itypat)%shapefncg,indl,l_size,l_size,&
343 &     pawtab(itypat)%mqgrid_shp,pawtab(itypat)%qgrid_shp,pawrad(itypat),rgl)
344      pawtab(itypat)%shapefncg=four_pi*pawtab(itypat)%shapefncg
345      ABI_DEALLOCATE(indl)
346      ABI_DEALLOCATE(rgl)
347    else
348      pawtab(itypat)%mqgrid_shp=0
349    end if
350 
351 !  ==================================================
352 !  3- COMPUTE indklmn INDEXES GIVING klm, kln, abs(il-jl) and (il+jl), ilmn and jlmn
353 !  for each klmn=(ilmn,jlmn)
354 
355    if (allocated(pawtab(itypat)%indklmn))  then
356      ABI_DEALLOCATE(pawtab(itypat)%indklmn)
357    end if
358    ABI_ALLOCATE(pawtab(itypat)%indklmn,(8,lmn2_size))
359 
360    klm_diag=0
361    do jlmn=1,lmn_size
362      jl= indlmn(1,jlmn);jlm=indlmn(4,jlmn);jln=indlmn(5,jlmn)
363      j0lmn=jlmn*(jlmn-1)/2
364      j0lm =jlm *(jlm -1)/2
365      j0ln =jln *(jln -1)/2
366      do ilmn=1,jlmn
367        il= indlmn(1,ilmn);ilm=indlmn(4,ilmn);iln=indlmn(5,ilmn)
368        klmn=j0lmn+ilmn
369        if (ilm<=jlm) then
370          indklmn_(1,klmn)=j0lm+ilm
371        else
372          i0lm=ilm*(ilm-1)/2
373          indklmn_(1,klmn)=i0lm+jlm
374        end if
375        if (iln<=jln) then
376          indklmn_(2,klmn)=j0ln+iln
377        else
378          i0ln=iln*(iln-1)/2
379          indklmn_(2,klmn)=i0ln+jln
380        end if
381        indklmn_(3,klmn)=min(abs(il-jl),lcutdens)
382        indklmn_(4,klmn)=min(il+jl,lcutdens)
383        indklmn_(5,klmn)=ilm
384        indklmn_(6,klmn)=jlm
385        indklmn_(7,klmn)=ilmn
386        indklmn_(8,klmn)=jlmn
387        pawtab(itypat)%indklmn(:,klmn)=indklmn_(:,klmn)
388        if (ilm==jlm) klm_diag(klmn)=1
389      end do
390    end do
391 
392 !  ==================================================
393 !  4- COMPUTE various FACTORS/SIZES (depending on (l,m,n))
394 
395    pawtab(itypat)%usespnorb=pawspnorb
396    pawtab(itypat)%lcut_size=min(l_size,lcutdens+1)
397 
398    if (allocated(pawtab(itypat)%dltij))  then
399      ABI_DEALLOCATE(pawtab(itypat)%dltij)
400    end if
401    ABI_ALLOCATE(pawtab(itypat)%dltij,(lmn2_size))
402    pawtab(itypat)%dltij(:)=two
403    do ilmn=1,lmn_size
404      pawtab(itypat)%dltij(ilmn*(ilmn+1)/2)=one
405    end do
406 
407    lmnmix=zero
408    ABI_ALLOCATE(kmix_tmp,(lmn2_size))
409    do jlmn=1,lmn_size
410      jl=indlmn(1,jlmn)
411      if (jl<=lmix) then
412        j0lmn=jlmn*(jlmn-1)/2
413        do ilmn=1,jlmn
414          il=indlmn(1,ilmn)
415          if (il<=lmix) then
416            lmnmix=lmnmix+1
417            kmix_tmp(lmnmix)=j0lmn+ilmn
418          end if
419        end do
420      end if
421    end do
422    if (allocated(pawtab(itypat)%kmix))  then
423      ABI_DEALLOCATE(pawtab(itypat)%kmix)
424    end if
425    ABI_ALLOCATE(pawtab(itypat)%kmix,(lmnmix))
426    pawtab(itypat)%lmnmix_sz=lmnmix
427    pawtab(itypat)%kmix(1:lmnmix)=kmix_tmp(1:lmnmix)
428    ABI_DEALLOCATE(kmix_tmp)
429 
430 !  ==================================================
431 !  5- COMPUTE Qijl TERMS AND Sij MATRIX
432 
433 !  Store some usefull quantities
434    if (allocated(pawtab(itypat)%phiphj))  then
435      ABI_DEALLOCATE(pawtab(itypat)%phiphj)
436    end if
437    if (allocated(pawtab(itypat)%tphitphj))  then
438      ABI_DEALLOCATE(pawtab(itypat)%tphitphj)
439    end if
440    ABI_ALLOCATE(pawtab(itypat)%phiphj,(mesh_size,ij_size))
441    ABI_ALLOCATE(pawtab(itypat)%tphitphj,(mesh_size,ij_size))
442    do jln=1,basis_size
443      j0ln=jln*(jln-1)/2
444      do iln=1,jln
445        kln=j0ln+iln
446        pawtab(itypat)%phiphj  (1:mesh_size,kln)=pawtab(itypat)%phi (1:mesh_size,iln)&
447 &       *pawtab(itypat)%phi (1:mesh_size,jln)
448        pawtab(itypat)%tphitphj(1:mesh_size,kln)=pawtab(itypat)%tphi(1:mesh_size,iln)&
449 &       *pawtab(itypat)%tphi(1:mesh_size,jln)
450      end do
451    end do
452 
453 !  Compute q_ijL and S_ij=q_ij0
454    if (allocated(pawtab(itypat)%qijl))  then
455      ABI_DEALLOCATE(pawtab(itypat)%qijl)
456    end if
457    if (allocated(pawtab(itypat)%sij))  then
458      ABI_DEALLOCATE(pawtab(itypat)%sij)
459    end if
460    ABI_ALLOCATE(pawtab(itypat)%qijl,(l_size*l_size,lmn2_size))
461    ABI_ALLOCATE(pawtab(itypat)%sij,(lmn2_size))
462    pawtab(itypat)%qijl=zero
463    pawtab(itypat)%sij=zero
464    do klmn=1,lmn2_size
465      klm=indklmn_(1,klmn);kln=indklmn_(2,klmn)
466      lmin=indklmn_(3,klmn);lmax=indklmn_(4,klmn)
467      do ll=lmin,lmax,2
468        lm0=ll*ll+ll+1;ff(1)=zero
469        ff(2:mesh_size)=(pawtab(itypat)%phiphj  (2:mesh_size,kln)&
470 &       -pawtab(itypat)%tphitphj(2:mesh_size,kln))&
471 &       *rad(2:mesh_size)**ll
472        call simp_gen(intg,ff,pawrad(itypat))
473        do mm=-ll,ll
474          isel=pawang%gntselect(lm0+mm,klm)
475          if (isel>0) pawtab(itypat)%qijl(lm0+mm,klmn)=intg*pawang%realgnt(isel)
476        end do
477      end do
478      if (klm_diag(klmn)==1) pawtab(itypat)%sij(klmn)= &
479 &     pawtab(itypat)%qijl(1,klmn)*sqrt(four_pi)
480    end do
481 
482 !  ==================================================
483 !  6- COMPUTE Eijkl TERMS (Hartree)
484 !     Compute eventually short-range screened version of Eijkl (Fock)
485 
486    if (allocated(pawtab(itypat)%eijkl))  then
487      ABI_DEALLOCATE(pawtab(itypat)%eijkl)
488    end if
489    ABI_ALLOCATE(pawtab(itypat)%eijkl,(lmn2_size,lmn2_size))
490    if (abs(hyb_range_fock)>tol8) then
491      if (allocated(pawtab(itypat)%eijkl_sr))  then
492        ABI_DEALLOCATE(pawtab(itypat)%eijkl_sr)
493      end if
494      ABI_ALLOCATE(pawtab(itypat)%eijkl_sr,(lmn2_size,lmn2_size))
495    end if
496 
497 !  First loop is for eijkl (Hartree)
498 !  2nd loop is for eijkl_sr (short-range screened Fock exchange)
499    do iloop=1,2
500      if (iloop==2.and.abs(hyb_range_fock)<=tol8) cycle
501      if (iloop==1) eijkl => pawtab(itypat)%eijkl
502      if (iloop==2) eijkl => pawtab(itypat)%eijkl_sr
503 
504 !    Compute:
505 !    vhatL(r) according to eq. (A14) in Holzwarth et al., PRB 55, 2005 (1997) [[cite:Holzwarth1997]]
506 !    intvhatL=$\int_{0}^{r_c}{vhatL(r) shapefunc_L(r) r^2\,dr}$
507 !    vhatijL =$\int_{0}^{r_c}{vhatL(r) \tilde{\phi}_i \tilde{\phi}_j \,dr}$
508 !    -----------------------------------------------------------------
509      ABI_ALLOCATE(vhatl,(mesh_size))
510      ABI_ALLOCATE(vhatijl,(lmn2_size,l_size))
511      ABI_ALLOCATE(intvhatl,(l_size))
512      intvhatl(:)=zero;vhatl(:)=zero;vhatijl(:,:)=zero
513      do il=1,l_size
514        vhatl(1)=zero;ff(1)=zero
515        ff(2:mesh_size)=pawtab(itypat)%shapefunc(2:mesh_size,il)*rad(2:mesh_size)**2
516        if (iloop==1) call poisson(ff,il-1,pawrad(itypat),vhatl)
517        if (iloop==2) call poisson(ff,il-1,pawrad(itypat),vhatl,screened_sr_separation=hyb_range_fock)
518        vhatl(2:mesh_size)=two*vhatl(2:mesh_size)/rad(2:mesh_size)
519        gg(1:mesh_size)=vhatl(1:mesh_size)*ff(1:mesh_size)
520        call simp_gen(intvhatl(il),gg,pawrad(itypat))
521        do klmn=1,lmn2_size
522          kln=indklmn_(2,klmn)
523          hh(1:mesh_size)=vhatl(1:mesh_size)*pawtab(itypat)%tphitphj(1:mesh_size,kln)
524          call simp_gen(vhatijl(klmn,il),hh,pawrad(itypat))
525        end do
526      end do
527      ABI_DEALLOCATE(vhatl)
528 
529 !    Compute:
530 !    eijkl=$ vh1_ijkl - Vhatijkl - Bijkl - Cijkl$
531 !    With:
532 !          $vh1_ijkl =\sum_{L,m} {vh1*Gaunt(i,j,Lm)*Gaunt(k,l,Lm)}$
533 !          $Vhat_ijkl=\sum_{L,m} {vhatijL*Gaunt(i,j,Lm)*q_klL}$
534 !          $B_ijkl   =\sum_{L,m} {vhatijL*Gaunt(k,l,Lm)*q_ijL}$
535 !          $C_ijkl   =\sum_{L,m} {intvhatL*q_ijL*q_klL}$
536 !    and:
537 !      vh1 according to eq. (A17) in Holzwarth et al., PRB 55, 2005 (1997) [[cite:Holzwarth1997]]
538 !    Warning: compute only eijkl for (i,j)<=(k,l)
539 !    -----------------------------------------------------------------
540      eijkl(:,:)=zero
541      meshsz=pawrad(itypat)%int_meshsz;if (mesh_size>meshsz) ff(meshsz+1:mesh_size)=zero
542      do klmn=1,lmn2_size
543        klm=indklmn_(1,klmn);kln=indklmn_(2,klmn)
544        lmin=indklmn_(3,klmn);lmax=indklmn_(4,klmn)
545        do ll=lmin,lmax,2
546          lm0=ll*ll+ll+1
547          ff(1:meshsz)=pawtab(itypat)%phiphj  (1:meshsz,kln)
548          if (iloop==1) call poisson(ff,ll,pawrad(itypat),gg)
549          if (iloop==2) call poisson(ff,ll,pawrad(itypat),gg,screened_sr_separation=hyb_range_fock)
550          ff(1:meshsz)=pawtab(itypat)%tphitphj(1:meshsz,kln)
551          if (iloop==1) call poisson(ff,ll,pawrad(itypat),hh)
552          if (iloop==2) call poisson(ff,ll,pawrad(itypat),hh,screened_sr_separation=hyb_range_fock)
553          do klmn1=klmn,lmn2_size
554            klm1=indklmn_(1,klmn1);kln1=indklmn_(2,klmn1)
555            lmin1=indklmn_(3,klmn1);lmax1=indklmn_(4,klmn1)
556            vh1=zero
557            if ((ll.ge.lmin1).and.(ll.le.lmax1)) then
558              ff(1)=zero
559              ff(2:meshsz)=(pawtab(itypat)%phiphj  (2:meshsz,kln1)*gg(2:meshsz)&
560 &             -pawtab(itypat)%tphitphj(2:meshsz,kln1)*hh(2:meshsz))&
561 &             *two/rad(2:meshsz)
562              call simp_gen(vh1,ff,pawrad(itypat))
563            end if
564            do mm=-ll,ll
565              isel =pawang%gntselect(lm0+mm,klm)
566              isel1=pawang%gntselect(lm0+mm,klm1)
567              if (isel>0.and.isel1>0) then
568                rg =pawang%realgnt(isel)
569                rg1=pawang%realgnt(isel1)
570                ql =pawtab(itypat)%qijl(lm0+mm,klmn)
571                ql1=pawtab(itypat)%qijl(lm0+mm,klmn1)
572                eijkl(klmn,klmn1)=eijkl(klmn,klmn1)&
573 &               +(   vh1                *rg *rg1&      ! vh1_ijkl
574 &              -    vhatijl(klmn ,ll+1)*rg *ql1&     ! Vhat_ijkl
575 &              -    vhatijl(klmn1,ll+1)*rg1*ql &     ! B_ijkl
576 &              -    intvhatl(ll+1)     *ql *ql1&     ! C_ijkl
577 &              )*two_pi
578              end if
579            end do
580          end do
581        end do
582      end do
583      ABI_DEALLOCATE(vhatijl)
584      ABI_DEALLOCATE(intvhatl)
585    end do ! iloop
586 
587 !  ==================================================
588 !  7- COMPUTE gamma_ij TERMS
589 !  corrections to get the background right
590 
591    if (pawtab(itypat)%usepotzero==1) then
592      if (allocated(pawtab(itypat)%gammaij))  then
593        ABI_DEALLOCATE(pawtab(itypat)%gammaij)
594      end if
595      ABI_ALLOCATE(pawtab(itypat)%gammaij,(lmn2_size))
596      ABI_ALLOCATE(work,(mesh_size))
597      do klmn=1,lmn2_size
598        if (klm_diag(klmn)==1) then
599          kln=indklmn_(2,klmn)
600          ff(1)=zero
601          ff(2:mesh_size)=pawtab(itypat)%phiphj(2:mesh_size,kln)-pawtab(itypat)%tphitphj(2:mesh_size,kln)
602          ! First, compute q_ij^00
603          call simp_gen(intg,ff,pawrad(itypat))
604          ! Second, compute phi^2 - tphi^2 - 4pi*r^2 q_ij^00 g_0(r)
605          ff(2:mesh_size)= ff(2:mesh_size) - intg*pawtab(itypat)%shapefunc(2:mesh_size,1)*rad(2:mesh_size)**2
606          call poisson(ff,0,pawrad(itypat),work)
607          ! work is r*vh; should be then multiplied by r to prepare the integration in the sphere
608          work(1)=zero ; work(2:mesh_size)=work(2:mesh_size)*rad(2:mesh_size)
609          ! Third, average it over the sphere
610          call simp_gen(intg,work,pawrad(itypat))
611          ! Finally, store it in pawtab%gammaij
612          pawtab(itypat)%gammaij(klmn)=intg*four_pi
613        else
614          pawtab(itypat)%gammaij(klmn)=zero
615        end if
616      end do
617      ABI_DEALLOCATE(work)
618    end if
619 
620 !  ***********************
621 !  End Loop on atom types
622 !  ***********************
623    ABI_DEALLOCATE(ff)
624    ABI_DEALLOCATE(gg)
625    ABI_DEALLOCATE(hh)
626    ABI_DEALLOCATE(indklmn_)
627    ABI_DEALLOCATE(klm_diag)
628    ABI_DEALLOCATE(rad)
629  end do
630 
631  call timab(553,2,tsec)
632 
633  DBG_EXIT("COLL")
634 
635 end subroutine pawinit