!{\src2tex{textfont=tt}}
!!****f* ABINIT/make_efg_el
!! NAME
!! make_efg_el
!!
!! FUNCTION
!! compute the electric field gradient due to electron density
!!
!! COPYRIGHT
!! Copyright (C) 2005-2014 ABINIT group (JJ)
!! 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
!! gcart(ngfft(1),ngfft(2),ngfft(3),3), G vectors in cartesian space
!! natom, number of atoms in unit cell
!! nfft,ngfft(18), number of FFT points and details of FFT
!! nspden, number of spin components
!! nsym=number of symmetries in space group
!! rhor(nfft,nspden), valence electron density, here $\tilde{n} + \hat{n}$
!! 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 atomic site due to rhor
!!
!! 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). This 
!! routine computes the second derivatives of the potential generated by $\tilde{n}$ (see Kresse and Joubert
!! for notation, Fourier-transforming the density, doing the sum in G space, and then transforming back at
!! each atomic site. The final formula is
!! \begin{displaymath}
!! \frac{\partial^2 V}{\partial x_\alpha\partial x_\beta} = -4\pi^2\sum_G (G_\alpha G_\beta - \delta_{\alpha,\beta}G^2/3) 
!! \left(\frac{\tilde{n}(G)}{\pi G^2}\right)e^{2\pi i G\cdot R}
!! \end{displaymath}
!! 
!!
!! PARENTS
!!      calc_efg
!!
!! CHILDREN
!!      destroy_mpi_enreg,fourdp,init_distribfft_seq,initmpi_seq,xred2xcart
!!
!! SOURCE

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

subroutine make_efg_el(efg,gcart,natom,nfft,ngfft,nspden,nsym,paral_kgb,rhor,rprimd,symrel,tnons,xred)

 use defs_basis
 use defs_abitypes
 use m_profiling

 use m_mpinfo,            only : destroy_mpi_enreg

 use m_mpinfo,            only : destroy_mpi_enreg

!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_el'
 use interfaces_41_geometry
 use interfaces_45_geomoptim
 use interfaces_51_manage_mpi
 use interfaces_53_ffts
!End of the abilint section

 implicit none

!Arguments ------------------------------------
!scalars
 integer,intent(in) :: natom,nfft,nspden,nsym,paral_kgb
!arrays
 integer,intent(in) :: ngfft(18),symrel(3,3,nsym)
 real(dp),intent(in) :: gcart(ngfft(1),ngfft(2),ngfft(3),3)
 real(dp),intent(in) :: rhor(nfft,nspden),rprimd(3,3),tnons(3,nsym),xred(3,natom)
 real(dp),intent(out) :: efg(3,3,natom)

!Local variables-------------------------------
!scalars
 integer :: cplex,iatom,igfft1,igfft2,igfft3,ii,index,isign,jj
 integer :: tim_fourdp
 real(dp) :: cph,derivs,phase,sph,trace
 type(MPI_type) :: mpi_enreg_seq
!arrays
 real(dp) :: gvec(3),ratom(3)
 real(dp),allocatable :: efg_ij_c(:,:,:,:),fofg(:,:),fofr(:),xcart(:,:)

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

!DEBUG
!write(std_out,*)' make_efg_el : enter'
!ENDDEBUG

 ABI_ALLOCATE(fofg,(2,nfft))
 ABI_ALLOCATE(fofr,(nfft))
 ABI_ALLOCATE(efg_ij_c,(2,3,3,natom))
 ABI_ALLOCATE(xcart,(3,natom))

 efg(:,:,:) = zero
 efg_ij_c(:,:,:,:) = zero
 call xred2xcart(natom,rprimd,xcart,xred) ! get atomic locations in cartesian coords

 tim_fourdp = 0 ! timing code, not using
 isign = -1 ! FT from R to G
 cplex = 1 ! fofr is real
!here we are only interested in the total charge density including nhat, which is rhor(:,1)
!regardless of the value of nspden. This may change in the future depending on 
!developments with noncollinear magnetization and so forth. Such a change will
!require an additional loop over nspden.
!Multiply by -1 to convert the electron density to the charge density
 fofr(:) = -rhor(:,1)
 call initmpi_seq(mpi_enreg_seq) !fake mpi type
 call init_distribfft_seq(mpi_enreg_seq%distribfft,'c',ngfft(2),ngfft(3),'fourdp')
 call fourdp(cplex,fofg,fofr,isign,mpi_enreg_seq,nfft,ngfft,paral_kgb,tim_fourdp) ! construct charge density in G space
 call destroy_mpi_enreg(mpi_enreg_seq)

 do igfft1 = 1, ngfft(1) ! summing over G space vector components...
   do igfft2 = 1, ngfft(2)
     do igfft3 = 1, ngfft(3)
       index = (igfft3-1)*ngfft(2)*ngfft(1) + (igfft2-1)*ngfft(1) + igfft1
       gvec(:) = gcart(igfft1,igfft2,igfft3,:) ! gvec is the current vector in G space
       trace = dot_product(gvec,gvec)
       if (trace > zero) then ! avoid G = 0 point
         do ii = 1, 3 ! sum over components of efg tensor
           do jj = 1, 3 ! sum over components of efg tensor
             derivs = gvec(ii)*gvec(jj) ! This term is $G_\alpha G_\beta$
            if (ii == jj) derivs = derivs - trace/three
             do iatom = 1, natom ! sum over atoms in unit cell
               ratom(:) = xcart(:,iatom) ! extract location of atom iatom
               phase = two_pi*dot_product(gvec,ratom) ! argument of $e^{2\pi i G\cdot R}$
               cph = cos(phase)
               sph = sin(phase)
               efg_ij_c(1,ii,jj,iatom) = efg_ij_c(1,ii,jj,iatom) - &
&               four_pi*derivs*(fofg(1,index)*cph-fofg(2,index)*sph)/trace ! real part of efg tensor
               efg_ij_c(2,ii,jj,iatom) = efg_ij_c(2,ii,jj,iatom) - &
&               four_pi*derivs*(fofg(1,index)*sph+fofg(2,index)*cph)/trace ! imaginary part, should be zero if G grid dense enough
             end do ! end loop over atoms in cell
           end do ! end loop over jj in V_ij
         end do ! end loop over ii in V_ij
       end if ! end statement avoiding G = 0 point
     end do ! end loop over ngfft(3)
   end do ! end loop over ngfft(2)
 end do ! end loop over ngfft(1)

 efg(:,:,:) = efg_ij_c(1,:,:,:) ! extract real part only of final tensor

! 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

 ABI_DEALLOCATE(fofg)
 ABI_DEALLOCATE(fofr)
 ABI_DEALLOCATE(efg_ij_c)
 ABI_DEALLOCATE(xcart)

!DEBUG
!write(std_out,*)' make_efg_el : exit '
!stop
!ENDDEBUG

end subroutine make_efg_el
!!***
