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-2022 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

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

SOURCE

735 subroutine paw_gencond(Dtset,gnt_option,mode,call_pawinit)
736 
737 !Arguments ------------------------------------
738  integer,intent(in) :: gnt_option
739  logical,intent(out) :: call_pawinit
740  character(len=*),intent(in) :: mode
741  type(dataset_type),intent(in) :: Dtset
742 
743 !Local variables-------------------------------
744 !scalars
745  integer,save :: gencond(10)=(/-1,-1,-1,-1,-1,-1,-1,-1,-1,-1/)
746 
747 ! *********************************************************************
748 
749  call_pawinit = .False.
750  select case (mode)
751  case ("test")
752 
753    if (gencond(1)/=Dtset%pawlcutd  .or.gencond(2) /=Dtset%pawlmix  .or.&
754 &      gencond(3)/=Dtset%pawnphi   .or.gencond(4) /=Dtset%pawntheta.or.&
755 &      gencond(5)/=Dtset%pawspnorb .or.gencond(6) /=Dtset%pawxcdev .or.&
756 &      gencond(7)/=Dtset%nsym      .or.gencond(8) /=gnt_option     .or.&
757 &      gencond(9)/=Dtset%usepotzero.or.gencond(10)/=Dtset%usekden) call_pawinit = .True.
758 
759  case ("save")
760     ! Update internal values
761    gencond(1)=Dtset%pawlcutd  ; gencond(2) =Dtset%pawlmix
762    gencond(3)=Dtset%pawnphi   ; gencond(4) =Dtset%pawntheta
763    gencond(5)=Dtset%pawspnorb ; gencond(6) =Dtset%pawxcdev
764    gencond(7)=Dtset%nsym      ; gencond(8) =gnt_option
765    gencond(9)=Dtset%usepotzero; gencond(10)=Dtset%usekden
766 
767  case ("reset")
768    gencond = -1
769 
770  case default
771    ABI_BUG("Wrong value for mode: "//trim(mode))
772  end select
773 
774 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-2022 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

  effmass_free=effective mass for electrons (1. in common case)
  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
  ixc=choice of exchange-correlation functional
  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)
  usekden= 1 is kinetic energy density has to be computed, 0 otherwise

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)

SOURCE

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