TABLE OF CONTENTS


ABINIT/m_paw_onsite [ Modules ]

[ Top ] [ Modules ]

NAME

  m_paw_onsite

FUNCTION

  This module contains a set of routines to compute various PAW on-site quantities
  i.e. quantities expressed with <Phi_i|...|Phi_j> and/or <tild_Phi_i|...|tild_Phi_j>.

COPYRIGHT

 Copyright (C) 2013-2022 ABINIT group (MT,FJ)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .

NOTES

  FOR DEVELOPPERS: in order to preserve the portability of libPAW library,
  please consult ~abinit/src/??_libpaw/libpaw-coding-rules.txt

SOURCE

21 #include "libpaw.h"
22 
23 MODULE m_paw_onsite
24 
25  USE_DEFS
26  USE_MSG_HANDLING
27  USE_MEMORY_PROFILING
28 
29  use m_pawrad,      only : pawrad_type, pawrad_deducer0, simp_gen, nderiv_gen
30  use m_pawtab,      only : pawtab_type
31  use m_paw_sphharm, only : setnabla_ylm
32 
33  implicit none
34 
35  private
36 
37 !public procedures.
38  public ::  pawnabla_init      ! Evaluate valence-valence on-site contribs of the nabla operator in cart. coord.
39  public ::  pawnabla_core_init ! Evaluate core-valence on-site contribs of the nabla operator in cart. coord.

m_paw_onsite/pawnabla_core_init [ Functions ]

[ Top ] [ m_paw_onsite ] [ Functions ]

NAME

 pawnabla_core_init

FUNCTION

 Evaluate core-valence onsite contributions of the nabla operator in cartesian coordinates,
  i.e. <Phi_i|Nabla|Phi_core_j>-<tPhi_i|Nabla|tPhi_core_j>.
 Core wave-functions are only given for one atom type.

INPUTS

  mpsang=1+maximum angular momentum
  ntypat=Number of types of atoms in cell
  Pawrad(ntypat)<Pawrad_type>=PAW radial mesh and related data:
    %rad(mesh_size)=The coordinates of all the points of the radial mesh
  Pawtab(ntypat) <type(pawtab_type>=PAW tabulated starting data:
    %mesh_size=Dimension of radial mesh
    %lmn_size=Number of (l,m,n) elements for the PAW basis
  phi_cor(mesh_size,nphicor)=core wave-functions for the current type of atoms.
  indlmn_cor(6,nlmn_core)= array giving l,m,n,lm,ln,s for i=lmn, for the core wave functions.

OUTPUT

  See side effects

SIDE EFFECTS

  Pawtab(ntypat) <type(pawtab_type>=PAW tabulated starting data:
    %has_nabla=set to 1 in matrix elements are calculated and stored
    %nabla_ij(3,lmn_size,lmn_size)= <phi_i|nabla|phi_core_j>-<tphi_i|nabla|tphi_core_j>

NOTES

  MG extracted this piece of code from optics_paw.F90 in order to have something more
  reusable! Note however the storage mode of nabla_ij differs from optics_paw
  (here Cartesian coordinates run faster). Besides nabla_ij contains the matrix
  elements of \nabla instead of the elements of the momentum operator p.

SOURCE

261 subroutine pawnabla_core_init(mpsang,ntypat,pawrad,pawtab,phi_cor,indlmn_cor,diracflag)
262 
263 !Arguments ------------------------------------
264 !scalars
265  integer,intent(in) :: mpsang,ntypat
266  integer,intent(in),optional :: diracflag
267 !arrays
268  integer,intent(in) :: indlmn_cor(:,:)
269  real(dp),intent(in) :: phi_cor(:,:)
270  type(pawtab_type),target,intent(inout) :: pawtab(ntypat)
271  type(pawrad_type),intent(in) :: pawrad(ntypat)
272 
273 !Local variables-------------------------------
274 !scalars
275  integer :: nln,nln_cor,ilm,ilmn,iln,itypat,sgnkappa
276  integer :: jl,jm,jlm,jlmn,jln,js,jlm_re,jlm_im,jm_re,jm_im
277  integer :: lmn_size,lmncmax,lcmax,ltmax,mesh_size,mesh_size_cor
278  real(dp) :: intg,jmj,cgc
279  logical :: dirac
280  character(len=500) :: msg
281 !arrays
282  integer, LIBPAW_CONTIGUOUS pointer :: indlmn(:,:)
283  real(dp) , allocatable:: ang_phipphj(:,:,:)
284  real(dp),allocatable :: dphi(:),ff(:),int1(:,:),int2(:,:),rad(:)
285 
286 ! *************************************************************************
287 
288  mesh_size_cor=size(phi_cor,1)
289  nln_cor=size(phi_cor,2)
290 
291  dirac=.false.
292  if(present(diracflag)) then
293     dirac=(diracflag==1)
294     if(diracflag==1.and.size(indlmn_cor,1)<8) then
295       msg='Wrong first dimension of indlmn_cor in pawnabla_core_init for diracrelativism!'
296       LIBPAW_BUG(msg)
297     endif
298  endif
299 
300 !To be checked
301  lmncmax=size(indlmn_cor,2) !Includes spinors if diracrel core wf are used
302  lcmax=maxval(indlmn_cor(1,:))
303  ltmax=max(lcmax+1,mpsang)
304  LIBPAW_ALLOCATE(ang_phipphj,(ltmax**2,ltmax**2,8))
305 
306  if (ltmax>4)then
307    write(msg,'(3a)')&
308 &   'Not designed for angular momentum greater than 3!',ch10,&
309 &   'Modification in the table defined in routine setnabla_ylm is required.'
310    LIBPAW_BUG(msg)
311  end if
312 
313 !if (mesh_size_cor/=pawrad(1)%mesh_size) then
314 !  write(msg,'(a)') 'Wrong mesh_size_cor value (1)!'
315 !  LIBPAW_BUG(msg)
316 !end if
317 !if (any(mesh_size_cor/=pawtab(:)%mesh_size)) then
318 !  write(msg,'(3a)') 'Wrong mesh_size_cor value (2)!',ch10,&
319 !&                    'Should have only one type of atom.'
320 !  LIBPAW_ERROR(msg)
321 !end if
322 
323 !Integration of the angular part: all angular integrals have been computed
324 !outside Abinit and tabulated for each (l,m) value up to l=3
325  call setnabla_ylm(ang_phipphj,ltmax)
326 
327  do itypat=1,ntypat
328 
329 !  COMPUTE nabla_ij := <phi_i|nabla|phi_cor> for this type
330    mesh_size=min(pawtab(itypat)%partialwave_mesh_size,pawrad(itypat)%mesh_size)
331    mesh_size=min(mesh_size_cor,mesh_size)
332    lmn_size=pawtab(itypat)%lmn_size
333    nln=pawtab(itypat)%basis_size
334 
335    if (mesh_size_cor<mesh_size) then
336      msg='mesh_size and mesh_sier_cor not compatible!'
337      LIBPAW_BUG(msg)
338    endif
339 
340    if (allocated(pawtab(itypat)%nabla_ij)) then
341      LIBPAW_DEALLOCATE(pawtab(itypat)%nabla_ij)
342    end if
343    LIBPAW_ALLOCATE(pawtab(itypat)%nabla_ij,(3,lmn_size,lmncmax))
344 
345    if (dirac) then
346      if (allocated(pawtab(itypat)%nabla_im_ij)) then
347        LIBPAW_DEALLOCATE(pawtab(itypat)%nabla_im_ij)
348      end if
349      LIBPAW_ALLOCATE(pawtab(itypat)%nabla_im_ij,(3,lmn_size,lmncmax))
350    end if
351 
352    pawtab(itypat)%has_nabla=1
353 
354    LIBPAW_ALLOCATE(ff,(mesh_size))
355    LIBPAW_ALLOCATE(rad,(mesh_size))
356    LIBPAW_ALLOCATE(int1,(lmn_size,lmncmax))
357    LIBPAW_ALLOCATE(int2,(lmn_size,lmncmax))
358    LIBPAW_ALLOCATE(dphi,(mesh_size))
359 
360    indlmn => pawtab(itypat)%indlmn
361    rad(1:mesh_size)=pawrad(itypat)%rad(1:mesh_size)
362 
363 !  int1= \int  Phi d/dr(Phi_core) r^2 dr
364 !      = \int (phi d/dr(phi_core) - phi phj_core/r) dr
365 !    with Phi=phi/r and Phi_core=phi_core/r
366    do jln=1,nln_cor
367      ff(1:mesh_size)=phi_cor(1:mesh_size,jln)
368      call nderiv_gen(dphi,ff,pawrad(itypat))
369      do iln=1,nln
370        ff(2:mesh_size)=pawtab(itypat)%phi(2:mesh_size,iln)*dphi(2:mesh_size) &
371 &       -pawtab(itypat)%phi(2:mesh_size,iln)*phi_cor(2:mesh_size,jln)/rad(2:mesh_size)
372        call pawrad_deducer0(ff,mesh_size,pawrad(itypat))
373        call simp_gen(intg,ff,pawrad(itypat))
374        int1(iln,jln)=intg
375      end do
376    end do
377 
378 !  int2= \int Phi Phi_core /r r^2 dr = \int phi phi_core /r dr
379 !    with Phi=phi/r and Phi_core=phi_core/r
380    do jln=1,nln_cor
381      do iln=1,nln
382        ff(2:mesh_size)=(pawtab(itypat)%phi(2:mesh_size,iln)*phi_cor(2:mesh_size,jln))/rad(2:mesh_size)
383        call pawrad_deducer0(ff,mesh_size,pawrad(itypat))
384        call simp_gen(intg,ff,pawrad(itypat))
385        int2(iln,jln)=intg
386      end do
387    end do
388 
389 !  ===== FULLY-RELATIVISTIC =====
390    if(dirac) then
391 
392 !    Integration of the radial part, Note unpacked loop
393      do jlmn=1,lmncmax
394        jl=indlmn_cor(1,jlmn)
395        jm=indlmn_cor(2,jlmn)
396 
397        sgnkappa=indlmn_cor(3,jlmn)
398        jmj=half*indlmn_cor(8,jlmn) ! 2mj is stored in indlmn_cor
399        js=indlmn_cor(6,jlmn)       ! 1 is up, 2 is down
400 
401 !      Calculate spinor dependend coeffs
402 !        (Clebsch-Gordan, I guess)
403        cgc=one ! so nothing changes without core spinors
404        if (sgnkappa==1) then
405          if(js==1) then
406            cgc= sqrt((dble(jl)-jmj+half)/dble(2*jl+1))
407          else
408            cgc=-sqrt((dble(jl)+jmj+half)/dble(2*jl+1))
409          endif
410        else
411          if(js==1) then
412            cgc= sqrt((dble(jl)+jmj+half)/dble(2*jl+1))
413          else
414            cgc= sqrt((dble(jl)-jmj+half)/dble(2*jl+1))
415          endif
416        endif
417 
418        jlm=indlmn_cor(4,jlmn)
419        jln=indlmn_cor(5,jlmn)
420 
421        do ilmn=1,lmn_size
422          ilm=indlmn(4,ilmn)
423          iln=indlmn(5,ilmn)
424 
425 !        jl was set as a flag for invalid combinations
426 !          i.e. m=-(l+1) or m=(l+1)
427 !        In these cases, cgc=0 ; so nabla_ij=0 
428          if(jl==-1) then
429            pawtab(itypat)%nabla_ij(1:3,ilmn,jlmn)= zero
430            pawtab(itypat)%nabla_im_ij(1:3,ilmn,jlmn) = zero
431 
432          else
433 
434 !          if jm<>0, need to convert from complex
435 !            to real spherical harmonics
436            if(jm<0) then
437              jm_re=abs(jm)
438              jm_im=-abs(jm)
439              jlm_re=jl*(jl+1)+jm_re+1
440              jlm_im=jl*(jl+1)+jm_im+1
441              pawtab(itypat)%nabla_ij(1,ilmn,jlmn)=half_sqrt2*cgc*( &
442 &              int1(iln,jln)* ang_phipphj(ilm,jlm_re,1) &
443 &             +int2(iln,jln)*(ang_phipphj(ilm,jlm_re,2)+ang_phipphj(ilm,jlm_re,3)))
444              pawtab(itypat)%nabla_ij(2,ilmn,jlmn)=half_sqrt2*cgc*( &
445 &              int1(iln,jln)* ang_phipphj(ilm,jlm_re,4) &
446 &             +int2(iln,jln)*(ang_phipphj(ilm,jlm_re,5)+ang_phipphj(ilm,jlm_re,6)))
447              pawtab(itypat)%nabla_ij(3,ilmn,jlmn)=half_sqrt2*cgc*( &
448 &              int1(iln,jln)* ang_phipphj(ilm,jlm_re,7) &
449 &             +int2(iln,jln)* ang_phipphj(ilm,jlm_re,8))
450              pawtab(itypat)%nabla_im_ij(1,ilmn,jlmn)=-half_sqrt2*cgc*( &
451 &              int1(iln,jln)* ang_phipphj(ilm,jlm_im,1) &
452 &             +int2(iln,jln)*(ang_phipphj(ilm,jlm_im,2)+ang_phipphj(ilm,jlm_im,3)))
453              pawtab(itypat)%nabla_im_ij(2,ilmn,jlmn)=-half_sqrt2*cgc*( &
454 &              int1(iln,jln)* ang_phipphj(ilm,jlm_im,4) &
455 &             +int2(iln,jln)*(ang_phipphj(ilm,jlm_im,5)+ang_phipphj(ilm,jlm_im,6)))
456              pawtab(itypat)%nabla_im_ij(3,ilmn,jlmn)=-half_sqrt2*cgc*( &
457 &              int1(iln,jln)* ang_phipphj(ilm,jlm_im,7) &
458 &             +int2(iln,jln)* ang_phipphj(ilm,jlm_im,8))
459 
460            else if (jm>0) then
461              jm_re=abs(jm)
462              jm_im=-abs(jm)
463              jlm_re=jl*(jl+1)+jm_re+1
464              jlm_im=jl*(jl+1)+jm_im+1
465              pawtab(itypat)%nabla_ij(1,ilmn,jlmn)=((-1)**jm)*half_sqrt2*cgc*( &
466 &              int1(iln,jln)* ang_phipphj(ilm,jlm_re,1) &
467 &             +int2(iln,jln)*(ang_phipphj(ilm,jlm_re,2)+ang_phipphj(ilm,jlm_re,3)))
468              pawtab(itypat)%nabla_ij(2,ilmn,jlmn)=((-1)**jm)*half_sqrt2*cgc*( &
469 &              int1(iln,jln)* ang_phipphj(ilm,jlm_re,4) &
470 &             +int2(iln,jln)*(ang_phipphj(ilm,jlm_re,5)+ang_phipphj(ilm,jlm_re,6)))
471              pawtab(itypat)%nabla_ij(3,ilmn,jlmn)=((-1)**jm)*half_sqrt2*cgc*( &
472 &              int1(iln,jln)* ang_phipphj(ilm,jlm_re,7) &
473 &             +int2(iln,jln)* ang_phipphj(ilm,jlm_re,8))
474              pawtab(itypat)%nabla_im_ij(1,ilmn,jlmn)=((-1)**jm)*half_sqrt2*cgc*( &
475 &              int1(iln,jln)* ang_phipphj(ilm,jlm_im,1) &
476 &             +int2(iln,jln)*(ang_phipphj(ilm,jlm_im,2)+ang_phipphj(ilm,jlm_im,3)))
477              pawtab(itypat)%nabla_im_ij(2,ilmn,jlmn)=((-1)**jm)*half_sqrt2*cgc*( &
478 &              int1(iln,jln)* ang_phipphj(ilm,jlm_im,4) &
479 &             +int2(iln,jln)*(ang_phipphj(ilm,jlm_im,5)+ang_phipphj(ilm,jlm_im,6)))
480              pawtab(itypat)%nabla_im_ij(3,ilmn,jlmn)=((-1)**jm)*half_sqrt2*cgc*( &
481 &              int1(iln,jln)* ang_phipphj(ilm,jlm_im,7) &
482 &             +int2(iln,jln)* ang_phipphj(ilm,jlm_im,8))
483 
484            else ! jm=0 : no conversion necessary if m=0
485              pawtab(itypat)%nabla_ij(1,ilmn,jlmn)=cgc*( &
486 &              int1(iln,jln)* ang_phipphj(ilm,jlm,1) &
487 &             +int2(iln,jln)*(ang_phipphj(ilm,jlm,2)+ang_phipphj(ilm,jlm,3)))
488              pawtab(itypat)%nabla_ij(2,ilmn,jlmn)=cgc*( &
489 &              int1(iln,jln)* ang_phipphj(ilm,jlm,4) &
490 &             +int2(iln,jln)*(ang_phipphj(ilm,jlm,5)+ang_phipphj(ilm,jlm,6)))
491              pawtab(itypat)%nabla_ij(3,ilmn,jlmn)=cgc*( &
492 &              int1(iln,jln)* ang_phipphj(ilm,jlm,7) &
493 &             +int2(iln,jln)* ang_phipphj(ilm,jlm,8))
494              pawtab(itypat)%nabla_im_ij(1,ilmn,jlmn)= zero
495              pawtab(itypat)%nabla_im_ij(2,ilmn,jlmn)= zero
496              pawtab(itypat)%nabla_im_ij(3,ilmn,jlmn)= zero
497            end if
498 
499          endif ! jl==-1?
500 
501        end do !ilmn
502      end do !jlmn
503 
504      pawtab(itypat)%has_nabla=4
505 
506 !  ===== NON-RELATIVISTIC OR SCALAR-RELATICISTIC =====
507    else
508 
509 !    Integration of the radial part, Note unpacked loop
510      do jlmn=1,lmncmax
511        jl=indlmn_cor(1,jlmn)
512        jlm=indlmn_cor(4,jlmn)
513        jln =indlmn_cor(5,jlmn)
514        do ilmn=1,lmn_size
515          ilm=indlmn(4,ilmn)
516          iln =indlmn(5,ilmn)
517          pawtab(itypat)%nabla_ij(1,ilmn,jlmn)=( &
518 &          int1(iln,jln)* ang_phipphj(ilm,jlm,1) &
519 &         +int2(iln,jln)*(ang_phipphj(ilm,jlm,2)+ang_phipphj(ilm,jlm,3)))
520 
521          pawtab(itypat)%nabla_ij(2,ilmn,jlmn)=( &
522 &          int1(iln,jln)* ang_phipphj(ilm,jlm,4) &
523 &         +int2(iln,jln)*(ang_phipphj(ilm,jlm,5)+ang_phipphj(ilm,jlm,6)))
524 
525          pawtab(itypat)%nabla_ij(3,ilmn,jlmn)=( &
526 &          int1(iln,jln)* ang_phipphj(ilm,jlm,7) &
527 &         +int2(iln,jln)* ang_phipphj(ilm,jlm,8))
528        end do !ilmn
529      end do !jlmn
530 
531      pawtab(itypat)%has_nabla=3
532 
533    end if ! Relativistic?
534 
535    LIBPAW_DEALLOCATE(ff)
536    LIBPAW_DEALLOCATE(rad)
537    LIBPAW_DEALLOCATE(int1)
538    LIBPAW_DEALLOCATE(int2)
539    LIBPAW_DEALLOCATE(dphi)
540 
541  end do !itypat
542 
543  LIBPAW_DEALLOCATE(ang_phipphj)
544 
545 end subroutine pawnabla_core_init

m_paw_onsite/pawnabla_init [ Functions ]

[ Top ] [ m_paw_onsite ] [ Functions ]

NAME

 pawnabla_init

FUNCTION

 Evaluate all valence-valence onsite contributions of the nabla operator in cartesian coordinates,
  i.e. <Phi_i|Nabla|Phi_j>-<tPhi_i|Nabla|tPhi_j>.

INPUTS

  mpsang=1+maximum angular momentum
  ntypat=Number of types of atoms in cell
  Pawrad(ntypat)<Pawrad_type>=PAW radial mesh and related data:
    %rad(mesh_size)=The coordinates of all the points of the radial mesh
  Pawtab(ntypat) <type(pawtab_type>=PAW tabulated starting data:
    %mesh_size=Dimension of radial mesh
    %lmn_size=Number of (l,m,n) elements for the PAW basis

OUTPUT

  See side effects

SIDE EFFECTS

  Pawtab(ntypat) <type(pawtab_type>=PAW tabulated starting data:
    %has_nabla=set to 1 in matrix elements are calculated and stored
    %nabla_ij(3,lmn_size,lmn_size)= <phi_i|nabla|phi_j>-<tphi_i|nabla|tphi_j>

NOTES

  MG extracted this piece of code from optics_paw.F90 in order to have something more
  reusable! Note however the storage mode of nabla_ij differs from optics_paw
  (here Cartesian coordinates run faster). Besides nabla_ij contains the matrix
  elements of \nabla instead of the elements of the momentum operator p.

SOURCE

 83 subroutine pawnabla_init(mpsang,ntypat,pawrad,pawtab)
 84 
 85 !Arguments ------------------------------------
 86 !scalars
 87  integer,intent(in) :: mpsang,ntypat
 88 !arrays
 89  type(pawtab_type),target,intent(inout) :: pawtab(ntypat)
 90  type(pawrad_type),intent(in) :: pawrad(ntypat)
 91 
 92 !Local variables-------------------------------
 93 !scalars
 94  integer :: ii,nln,il,ilm,ilmn,iln,itypat
 95  integer :: jl,jlm,jlmn,jln,lmn_size,mesh_size 
 96  real(dp) :: avg,intg
 97  character(len=500) :: msg
 98 !arrays
 99  integer, LIBPAW_CONTIGUOUS pointer :: indlmn(:,:)
100  real(dp) :: ang_phipphj(mpsang**2,mpsang**2,8)
101  real(dp),allocatable :: dphi(:),dtphi(:),ff(:),int1(:,:),int2(:,:),rad(:)
102  
103 ! *************************************************************************
104 
105  if (mpsang>4)then
106    write(msg,'(3a)')&
107 &   'Not designed for angular momentum greater than 3 ',ch10,&
108 &   'Modification in the table defined in routine setnabla_ylm is required.'
109    LIBPAW_BUG(msg)
110  end if
111 
112 !Integration of the angular part: all angular integrals have been computed 
113 !outside Abinit and tabulated for each (l,m) value up to l=3
114  call setnabla_ylm(ang_phipphj,mpsang)
115 
116  do itypat=1,ntypat
117 
118 !  COMPUTE nabla_ij := <phi_i|nabla|phi_j>-<tphi_i|nabla|tphi_j> for this type
119    mesh_size=pawtab(itypat)%mesh_size
120    lmn_size=pawtab(itypat)%lmn_size
121    nln=pawtab(itypat)%basis_size
122 
123    if (allocated(pawtab(itypat)%nabla_ij)) then
124      LIBPAW_DEALLOCATE(pawtab(itypat)%nabla_ij)
125    end if
126    LIBPAW_ALLOCATE(pawtab(itypat)%nabla_ij,(3,lmn_size,lmn_size))
127    pawtab(itypat)%has_nabla=1
128 
129    LIBPAW_ALLOCATE(ff,(mesh_size))
130    LIBPAW_ALLOCATE(rad,(mesh_size))
131    LIBPAW_ALLOCATE(int1,(lmn_size,lmn_size))
132    LIBPAW_ALLOCATE(int2,(lmn_size,lmn_size))
133    LIBPAW_ALLOCATE(dphi,(mesh_size))
134    LIBPAW_ALLOCATE(dtphi,(mesh_size))
135 
136    indlmn => pawtab(itypat)%indlmn
137    rad(1:mesh_size)=pawrad(itypat)%rad(1:mesh_size)
138 
139 !  int1= \int [ Phi d/dr(Phj) - tPhi d/dr(tPhj) ] r^2 dr
140 !      = \int [ (phi d/dr(phj) - phi phj/r) - (tphi d/dr(tphj) - tphi tphj/r) ] dr
141 !    with Phi=phi/r and tPhi=phi/r
142    do jln=1,nln
143      ff(1:mesh_size)=pawtab(itypat)%phi(1:mesh_size,jln)
144      call nderiv_gen(dphi,ff,pawrad(itypat))
145      ff(1:mesh_size)=pawtab(itypat)%tphi(1:mesh_size,jln)
146      call nderiv_gen(dtphi,ff,pawrad(itypat))
147      do iln=1,nln
148        ff(2:mesh_size)= &
149 &       pawtab(itypat)%phi (2:mesh_size,iln)*dphi (2:mesh_size) &
150 &       -pawtab(itypat)%phi (2:mesh_size,iln)*pawtab(itypat)%phi (2:mesh_size,jln)/rad(2:mesh_size) &
151 &       -( pawtab(itypat)%tphi(2:mesh_size,iln)*dtphi(2:mesh_size) &
152 &       -pawtab(itypat)%tphi(2:mesh_size,iln)*pawtab(itypat)%tphi(2:mesh_size,jln)/rad(2:mesh_size) )
153        call pawrad_deducer0(ff,mesh_size,pawrad(itypat))
154        call simp_gen(intg,ff,pawrad(itypat))
155        int1(iln,jln)=intg
156      end do
157    end do
158 
159 !  int2= \int [ Phi Phj /r - \int tPhi tPhj /r ] r^2 dr
160 !      = \int [ phi phj /r - \int tphi tphj /r ] dr
161 !    with Phi=phi/r and tPhi=phi/r
162    do jln=1,nln
163      do iln=1,nln
164        ff(2:mesh_size)= ( &
165 &       pawtab(itypat)%phi (2:mesh_size,iln)*pawtab(itypat)%phi (2:mesh_size,jln) &
166 &       -pawtab(itypat)%tphi(2:mesh_size,iln)*pawtab(itypat)%tphi(2:mesh_size,jln) ) /rad(2:mesh_size)
167        call pawrad_deducer0(ff,mesh_size,pawrad(itypat))
168        call simp_gen(intg,ff,pawrad(itypat))
169        int2(iln,jln)=intg
170      end do
171    end do
172 
173 !  Integration of the radial part, Note unpacked loop
174    do jlmn=1,lmn_size
175      jlm=indlmn(4,jlmn)
176      jl =indlmn(5,jlmn)
177      do ilmn=1,lmn_size
178        ilm=indlmn(4,ilmn)
179        il =indlmn(5,ilmn)
180 
181        pawtab(itypat)%nabla_ij(1,ilmn,jlmn)= &
182 &        int1(il,jl)* ang_phipphj(ilm,jlm,1) &
183 &       +int2(il,jl)*(ang_phipphj(ilm,jlm,2)+ang_phipphj(ilm,jlm,3))
184 
185        pawtab(itypat)%nabla_ij(2,ilmn,jlmn)= &
186 &        int1(il,jl)* ang_phipphj(ilm,jlm,4) &
187 &       +int2(il,jl)*(ang_phipphj(ilm,jlm,5)+ang_phipphj(ilm,jlm,6))
188 
189        pawtab(itypat)%nabla_ij(3,ilmn,jlmn)= &
190 &        int1(il,jl)* ang_phipphj(ilm,jlm,7) &
191 &       +int2(il,jl)* ang_phipphj(ilm,jlm,8)
192 
193      end do !ilmn
194    end do !jlmn
195 
196 !  Symetrization
197    if (lmn_size>1) then
198      do jlmn=2,lmn_size
199        do ilmn=1,jlmn-1
200          do ii=1,3
201            avg=half*(pawtab(itypat)%nabla_ij(ii,ilmn,jlmn)-pawtab(itypat)%nabla_ij(ii,jlmn,ilmn))
202            pawtab(itypat)%nabla_ij(ii,ilmn,jlmn)= avg
203            pawtab(itypat)%nabla_ij(ii,jlmn,ilmn)=-avg
204          end do           
205        end do
206      end do
207    end if
208 
209 !  End
210    pawtab(itypat)%has_nabla=2
211    LIBPAW_DEALLOCATE(ff)
212    LIBPAW_DEALLOCATE(rad)
213    LIBPAW_DEALLOCATE(int2)
214    LIBPAW_DEALLOCATE(int1)
215    LIBPAW_DEALLOCATE(dphi)
216    LIBPAW_DEALLOCATE(dtphi)
217 
218  end do !itypat
219 
220 end subroutine pawnabla_init