!{\src2tex{textfont=tt}}
!!****f* ABINIT/make_efg_onsite
!! NAME
!! make_efg_onsite
!!
!! FUNCTION
!! Compute the electric field gradient due to onsite densities
!!
!! COPYRIGHT
!! Copyright (C) 2005-2014 ABINIT group (JZ,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 .
!!
!! INPUTS
!!  mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
!!  mpi_comm_atom=--optional-- MPI communicator over atoms
!!  my_natom=number of atoms treated by current processor
!!  natom=number of atoms in cell.
!!  nsym=number of symmetries in space group
!!  ntypat=number of atom types
!!  paw_an(my_natom) <type(paw_an_type)>=paw arrays given on angular mesh
!!  pawang <type(pawang_type)>=paw angular mesh and related data
!!  pawrad(ntypat) <type(pawrad_type)>=paw radial mesh and related data
!!  pawrhoij(my_natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data
!!  pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
!!  rprimd(3,3), conversion from crystal coordinates to cartesian coordinates
!!  symrel(3,3,nsym)=symmetry operators in terms of action on primitive translations
!!  tnons(3,nsym) = nonsymmorphic translations
!!  xred(3,natom), location of atoms in crystal coordinates.
!!
!! OUTPUT
!!  efg(3,3,natom), the 3x3 efg tensor at each site due to nhat

!! NOTES
!! This routine computes the electric field gradient, specifically the components
!! $\partial^2 V/\partial x_\alpha \partial x_\beta$ of the potential generated by the valence
!! electrons, at each atomic site in the unit cell. Key references: Kresse and Joubert, ``From
!! ultrasoft pseudopotentials to the projector augmented wave method'', Phys. Rev. B. 59, 1758--1775 (1999),
!! and Profeta, Mauri, and Pickard, ``Accurate first principles prediction of $^{17}$O NMR parameters in
!! SiO$_2$: Assignment of the zeolite ferrierite spectrum'', J. Am. Chem. Soc. 125, 541--548 (2003). See in particular
!! Eq. 11 and 12 of Profeta et al., but note that their sum over occupied states times 2 for occupation number is
!! replaced in the Kresse and Joubert formulation by the sum over $\rho_{ij}$ occupations for each basis element pair.
!!
!! PARENTS
!!      calc_efg
!!
!! CHILDREN
!!      free_my_atmtab,get_my_atmtab,pawdensities,pawrad_deducer0,simp_gen
!!      xmpi_sum
!!
!! SOURCE

#if defined HAVE_CONFIG_H
#include "config.h"
#endif

#include "abi_common.h"

subroutine make_efg_onsite(efg,my_natom,natom,nsym,ntypat,paw_an,pawang,pawrhoij,pawrad,pawtab, &
&                          rprimd,symrel,tnons,xred,&
&                          mpi_atmtab,mpi_comm_atom) ! optional arguments (parallelism)

 use defs_basis
 use m_errors
 use m_profiling
 use m_xmpi, only : xmpi_sum

 use m_pawang, only : pawang_type
 use m_pawtab, only : pawtab_type
 use m_pawrad, only : pawrad_type, pawrad_deducer0, simp_gen
 use m_pawtab, only : pawtab_type
 use m_paw_an, only : paw_an_type
 use m_pawrhoij, only : pawrhoij_type
 use m_paral_atom, only : get_my_atmtab, free_my_atmtab

!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'make_efg_onsite'
 use interfaces_45_geomoptim
 use interfaces_66_paw, except_this_one => make_efg_onsite
!End of the abilint section

 implicit none

!Arguments ------------------------------------
!scalars
 integer,intent(in) :: my_natom,natom,nsym,ntypat
 integer,optional,intent(in) :: mpi_comm_atom
 type(pawang_type),intent(in) :: pawang
!arrays
 integer,intent(in) :: symrel(3,3,nsym)
 integer,optional,target,intent(in) :: mpi_atmtab(:)
 real(dp),intent(in) :: rprimd(3,3),tnons(3,nsym),xred(3,natom)
 real(dp),intent(out) :: efg(3,3,natom)
 type(paw_an_type),intent(in) :: paw_an(my_natom)
 type(pawrad_type),intent(in) :: pawrad(ntypat)
 type(pawrhoij_type),intent(in) :: pawrhoij(my_natom)
 type(pawtab_type),intent(in) :: pawtab(ntypat)

!Local variables-------------------------------
!scalars
 integer :: cplex,iatom,iatom_tot,ictr,ierr,imesh_size,ispden,itypat
 integer :: lm,lm_size,lnspden,mesh_size,nzlmopt,nspden
 integer :: opt_compch,opt_dens,opt_l,opt_print
 logical :: my_atmtab_allocated,paral_atom
 real(dp) :: c1,c2,c3,compch_sph,intg
!arrays
 integer,pointer :: my_atmtab(:)
 logical,allocatable :: lmselectin(:),lmselectout(:)
 real(dp),allocatable :: ff(:),nhat1(:,:,:),rho1(:,:,:),trho1(:,:,:)

! ************************************************************************

 DBG_ENTER("COLL")

 efg(:,:,:) = zero

!Set up parallelism over atoms
 paral_atom=(present(mpi_comm_atom).and.(my_natom/=natom))
 nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
 my_atmtab_allocated = .False.
 if (paral_atom) then
   call get_my_atmtab(mpi_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom)
 end if

!the following factors arise in expanding the angular dependence of the electric field gradient tensor in
!terms of real spherical harmonics. The real spherical harmonics are as in the routine initylmr.F90; see
!in particular also http://www.unioviedo.es/qcg/art/Theochem419-19-ov-BF97-rotation-matrices.pdf
 c1 = sqrt(16.0*pi/5.0)
 c2 = sqrt(4.0*pi/5.0)
 c3 = sqrt(12.0*pi/5.0)

!loop over atoms in cell
 do iatom = 1, my_natom
   iatom_tot=iatom;if (paral_atom) iatom_tot=my_atmtab(iatom)
   itypat=pawrhoij(iatom)%itypat

   lm_size = paw_an(iatom)%lm_size
   if (lm_size < 5) cycle ! if lm_size < 5 then the on-site densities for this atom have no L=2 component
!  and therefore nothing to contribute to the on-site electric field gradient

   mesh_size=pawrad(itypat)%mesh_size
   ABI_ALLOCATE(ff,(mesh_size))

   cplex = pawrhoij(iatom)%cplex
   nspden = pawrhoij(iatom)%nspden
   ABI_ALLOCATE(lmselectin,(lm_size))
   ABI_ALLOCATE(lmselectout,(lm_size))
   lmselectin = .true. ! compute all moments of densities
   nzlmopt = -1
   opt_compch = 0
   compch_sph = zero
   opt_dens = 0 ! compute all densities
   opt_l = -1 ! all moments contribute
   opt_print = 0 ! do not print out moments

   ABI_ALLOCATE(nhat1,(cplex*mesh_size,lm_size,nspden))
   ABI_ALLOCATE(rho1,(cplex*mesh_size,lm_size,nspden))
   ABI_ALLOCATE(trho1,(cplex*mesh_size,lm_size,nspden))

!  loop over spin components
!  nspden = 1: just a single component
!  nspden = 2: two components, loop over them
!  nspden = 4: total is in component 1, only one of interest
   if ( nspden == 2 ) then
     lnspden = 2
   else
     lnspden = 1
   end if
   do ispden=1,lnspden

!    construct multipole expansion of on-site charge densities for this atom
     call pawdensities(compch_sph,cplex,iatom_tot,lmselectin,lmselectout,lm_size,&
&     nhat1,nspden,nzlmopt,opt_compch,opt_dens,opt_l,opt_print,&
&     pawang,0,pawrad(itypat),pawrhoij(iatom),pawtab(itypat),&
&     rho1,trho1)

     do lm = 5, 9 ! loop on L=2 components of multipole expansion

       if(.not. lmselectout(lm)) cycle ! skip moments that contributes zero

!      the following is r^2*(n1-tn1-nhat)/r^3 for this multipole moment
!      use only the real part of the density in case of cplex==2
       do imesh_size = 2, mesh_size
         ictr = cplex*(imesh_size - 1) + 1
         ff(imesh_size)=(rho1(ictr,lm,ispden)-trho1(ictr,lm,ispden)-&
&         nhat1(ictr,lm,ispden))/&
&         pawrad(itypat)%rad(imesh_size)
       end do
       call pawrad_deducer0(ff,mesh_size,pawrad(itypat))
       call simp_gen(intg,ff,pawrad(itypat))
       select case (lm)
         case (5) ! S_{2,-2}
           efg(1,2,iatom_tot) = efg(1,2,iatom_tot) - c3*intg ! xy case
         case (6) ! S_{2,-1}
           efg(2,3,iatom_tot) = efg(2,3,iatom_tot) - c3*intg ! yz case
         case (7) ! S_{2, 0}
           efg(1,1,iatom_tot) = efg(1,1,iatom_tot) + c2*intg ! xx case
           efg(2,2,iatom_tot) = efg(2,2,iatom_tot) + c2*intg ! yy case
           efg(3,3,iatom_tot) = efg(3,3,iatom_tot) - c1*intg ! zz case
         case (8) ! S_{2,+1}
           efg(1,3,iatom_tot) = efg(1,3,iatom_tot) - c3*intg ! xz case
         case (9) ! S_{2,+2}
           efg(1,1,iatom_tot) = efg(1,1,iatom_tot) - c3*intg ! xx case
           efg(2,2,iatom_tot) = efg(2,2,iatom_tot) + c3*intg ! yy case
       end select

     end do  ! end loop over LM components with L=2

   end do    ! Loop on spin components

!  Symmetrization of EFG
   efg(2,1,iatom_tot) = efg(1,2,iatom_tot)
   efg(3,1,iatom_tot) = efg(1,3,iatom_tot)
   efg(3,2,iatom_tot) = efg(2,3,iatom_tot)

   ABI_DEALLOCATE(lmselectin)
   ABI_DEALLOCATE(lmselectout)
   ABI_DEALLOCATE(ff)
   ABI_DEALLOCATE(nhat1)
   ABI_DEALLOCATE(rho1)
   ABI_DEALLOCATE(trho1)

 end do     ! Loop on atoms

!Reduction in case of parallelisation over atoms
 if (paral_atom) then
   call xmpi_sum(efg,mpi_comm_atom,ierr)
 end if

! symmetrize tensor at each atomic site using point symmetry operations
 do iatom = 1, natom
   call matpointsym(iatom,efg(:,:,iatom),natom,nsym,rprimd,symrel,tnons,xred)
 end do

!Destroy atom table used for parallelism
 call free_my_atmtab(my_atmtab,my_atmtab_allocated)

 DBG_EXIT("COLL")

 end subroutine make_efg_onsite
!!***
