!{\src2tex{textfont=tt}}
!!****f* ABINIT/make_efg_ion
!! NAME
!! make_efg_ion
!!
!! FUNCTION
!! compute the electric field gradient due to ionic cores
!!
!! 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
!! natom, number of atoms in the unit cell
!! nsym=number of symmetries in space group
!! ntypat, the number of types of atoms in the unit cell
!! rprimd(3,3), the matrix giving the transformation from crystal to cartesian coordinates
!! symrel(3,3,nsym)=symmetry operators in terms of action on primitive translations
!! tnons(3,nsym) = nonsymmorphic translations
!! typat(natom), the type of each atom in the unit cell
!! ucvol, the volume of the unit cell in atomic units
!! xred(3,natom) the location of each atom in the cell in crystallographic coordinates
!! zion(ntypat) the net charge on each type of atom 
!!
!! OUTPUT
!! efg(3,3,natom), the 3x3 efg tensors at each atomic site
!!
!! SIDE EFFECTS
!!
!! 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 ionic cores, 
!! at each atomic site in the unit cell. 
!! Key references: 
!! 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); 
!! A. Honma, ``Dipolar lattice-sums with applications to the exciton bands of anthracene crystal and 
!! the crystal field due to point charges'', J. Phys. Soc. Jpn. 42, 1129--1135 (1977);
!! and Kresse and Joubert, ``From ultrasoft pseudopotentials to the projector augmented wave method'', 
!! Phys. Rev. B. 59, 1758--1775 (1999). In Kresse and Joubert's notation, the ionic cores are $n_{Zc}$; 
!! these charges are given by the net core charges on the pseudoatoms. Due to otherwise slow convergence, 
!! the sum over atoms is carried out by an Ewald method as detailed in the Honma reference, specifically 
!! his Eq. 4.8.
!!
!! PARENTS
!!      calc_efg
!!
!! CHILDREN
!!      matr3inv,xred2xcart
!!
!! SOURCE

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

#include "abi_common.h"

subroutine make_efg_ion(efg,natom,nsym,ntypat,rprimd,symrel,tnons,typat,ucvol,xred,zion)

 use defs_basis
 use m_profiling

 use m_special_funcs,  only : abi_derfc

!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_ion'
 use interfaces_32_util
 use interfaces_41_geometry
 use interfaces_45_geomoptim
!End of the abilint section

 implicit none

!Arguments ------------------------------------
!scalars
 integer,intent(in) :: natom,nsym,ntypat
 real(dp) :: ucvol
!arrays
 integer,intent(in) :: symrel(3,3,nsym),typat(natom)
 real(dp),intent(in) :: rprimd(3,3),tnons(3,nsym)
 real(dp),intent(in) :: zion(ntypat)
 real(dp),intent(inout) :: xred(3,natom)
 real(dp),intent(out) :: efg(3,3,natom)
!Local variables-------------------------------
!scalars
 integer :: iatom,ishell,ii,jatom,jj,nshell,sx,sy,sz
 real(dp) :: cph,dampfac,derfc_karg,derivs,gsq,karg
 real(dp) :: lenrho,phase,qk,rlkcut,trace,xi0
 real(dp) :: glkcut
!arrays
 real(dp) :: cvec(3),gvec(3),gpl(3),gprimd(3,3)
 real(dp) :: rhok(3),rhored(3),rpl(3)
 real(dp),allocatable :: efg_g(:,:,:),efg_r(:,:,:)
 real(dp),allocatable :: xcart(:,:)

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

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

 ABI_ALLOCATE(efg_g,(3,3,natom))
 ABI_ALLOCATE(efg_r,(3,3,natom))
 ABI_ALLOCATE(xcart,(3,natom))
 efg(:,:,:) = zero ! final efg tensor
 efg_g(:,:,:) = zero ! part of tensor accumulated in G space
 efg_r(:,:,:) = zero ! part of tensor accumulated in R space

 call xred2xcart(natom,rprimd,xcart,xred) ! get atomic locations in cartesian coords

 do ii = 1, 3 ! generate the lengths of the unit cell edges in atomic units
   rpl(ii) = sqrt(rprimd(1,ii)**2+rprimd(2,ii)**2+rprimd(3,ii)**2)
 end do
 xi0 = sqrt(pi/(maxval(rpl)*minval(rpl))) ! this estimate for xi0 is from Honma's paper

 call matr3inv(rprimd,gprimd) ! gprimd holds the inverse transpose of rprimd
!remember ordering: rprimd( (x_comp,y_comp,z_comp), (edge 1, edge 2, edge 3) )
!while gprimd( (edge 1, edge 2, edge 3),(x_comp, y_comp, z_comp) )
 do ii = 1, 3 ! generate the lengths of the reciprocal cell edges 
   gpl(ii) = sqrt(gprimd(ii,1)**2+gprimd(ii,2)**2+gprimd(ii,3)**2)
 end do

!go out enough shells such that g**2/4*xi0**2 is of order 30
 nshell = int(anint(sqrt(30.0)*xi0/(pi*minval(gpl))))
 glkcut = (0.95*nshell*two*pi*minval(gpl))**2

 do ishell = 0, nshell ! loop over shells
   do sx = -ishell, ishell
     do sy = -ishell, ishell
       do sz = -ishell, ishell
         if ( .not. (sx==0 .and. sy==0 .and. sz==0) ) then ! avoid origin
!          constrain to be on shell surface, not interior
           if ( abs(sx)==ishell .or. abs(sy)==ishell .or. abs(sz)==ishell ) then
             cvec(1)=sx;cvec(2)=sy;cvec(3)=sz
!            make the g vector in cartesian coords
             gvec(:) = zero
             do ii = 1, 3
               do jj = 1, 3
                 gvec(ii) = gvec(ii) + gprimd(ii,jj)*cvec(jj)*two*pi
               end do
             end do
             gsq = dot_product(gvec,gvec)
             if(gsq < glkcut) then
               dampfac = exp(-gsq/(4.0*xi0*xi0)) ! see Honma eq. 4.8
               do iatom = 1, natom
                 do jatom = 1, natom
                   qk = zion(typat(jatom)) ! charge on neighbor atom
                   rhok = xcart(:,jatom)-xcart(:,iatom)
                   phase = dot_product(gvec,rhok)
                   cph = cos(phase)
                   do ii = 1, 3
                     do jj = 1, 3
                       derivs = -3.0*gvec(ii)*gvec(jj)/gsq
                       if (ii == jj) derivs = 1.0 + derivs
                       efg_g(ii,jj,iatom) = efg_g(ii,jj,iatom) + &
&                       qk*cph*derivs*dampfac 
                     end do ! end loop over jj
                   end do ! end loop over ii
                 end do ! end loop over jatom
               end do ! end loop over iatom
             end if ! constrain to gsq < glkcut 
           end if ! end selection on shell edge
         end if ! end avoidance of origin
       end do ! end loop over sz
     end do ! end loop over sy
   end do ! end loop over sx
 end do ! end loop over ishell

!sum in real space begins here

!go out enough shells such that (r*xi0)**2 is of order 30
 nshell = int(anint(sqrt(30.)/(minval(rpl)*xi0)))
 rlkcut = nshell*minval(rpl)*0.95
!
!go out enough shells so that rlkcut is of order 30 bohr
!nshell=int(anint(30.0/minval(rpl)))
!rlkcut = 0.95*nshell*minval(rpl)

 do ishell = 0, nshell ! total set of cells to loop over
   do sx = -ishell, ishell ! loop over all cells in each dimension
     do sy = -ishell, ishell
       do sz = -ishell, ishell
!        constrain to shell surface, not interior
         if ( abs(sx)==ishell .or. abs(sy)==ishell .or. abs(sz)==ishell ) then
           do jatom = 1, natom ! loop over atoms in shell cell
             do iatom = 1, natom ! loop over atoms in central unit cell
               if (.NOT. (jatom == iatom .AND. sx == 0 .AND. sy == 0 .AND. sz == 0)) then ! avoid self term
                 qk = zion(typat(jatom)) ! charge on each neighbor atom
!                ! rhored is the vector in crystal coords from neighbor to target
                 rhored(1) = xred(1,jatom) + sx - xred(1,iatom) 
                 rhored(2) = xred(2,jatom) + sy - xred(2,iatom)
                 rhored(3) = xred(3,jatom) + sz - xred(3,iatom)
!                !  rhok is rhored in cartesian coords
                 rhok(1) = rprimd(1,1)*rhored(1)+rprimd(1,2)*rhored(2)+rprimd(1,3)*rhored(3) 
                 rhok(2) = rprimd(2,1)*rhored(1)+rprimd(2,2)*rhored(2)+rprimd(2,3)*rhored(3)
                 rhok(3) = rprimd(3,1)*rhored(1)+rprimd(3,2)*rhored(2)+rprimd(3,3)*rhored(3)
                 trace = dot_product(rhok,rhok)
                 lenrho = sqrt(trace)
                 if (lenrho < rlkcut) then ! this restriction is critical as it ensures 
!                  ! that we sum over a sphere of atoms in real space
!                  ! no matter what shape the unit cell has
                   karg = xi0*lenrho
                   derfc_karg = abi_derfc(karg)
!                  see Honma eq. 2.10 for derivation of the following damping factor
                   dampfac = (1.0+3.0/(2.0*karg*karg))*exp(-karg*karg)+3.0*sqrt(pi)*derfc_karg/(4.0*karg**3) 
                   do ii = 1, 3 ! loop over tensor elements
                     do jj = 1, 3 ! loop over tensor elements
                       derivs = -3.0*rhok(ii)*rhok(jj)/trace 
                       if(ii == jj) derivs = derivs + 1.0 ! see Honma eq 4.8 re: sign
!                      accumulate real space tensor element, 
!                      weighted by charge of neighbor and Ewald damping factor
                       efg_r(ii,jj,iatom) = efg_r(ii,jj,iatom) + qk*derivs*dampfac 
                     end do ! end loop over jj in efg(ii,jj,iatom)
                   end do ! end loop over ii in efg(ii,jj,iatom)
                 end if ! end if statement restricting to a sphere of radius rlkcut
               end if ! end if statement avoiding the self atom term
             end do ! end loop over i atoms in cell
           end do ! end loop over j atoms in cell
         end if ! end selection on outer shell of cells only
       end do ! end loop over sz cells
     end do ! end loop over sy cells
   end do ! end loop over sx cells
 end do ! end loop over shells

!now combine the g-space and r-space parts, properly weighted (see Honma)
 do iatom = 1, natom
   do ii = 1, 3
     do jj = 1, 3
       efg(ii,jj,iatom) = four_pi*efg_g(ii,jj,iatom)/(three*ucvol)-&
&       four*xi0**3*efg_r(ii,jj,iatom)/(three*sqrt(pi)) 
!      note extra factor of two: compare Honma eq. 4.6
     end do
   end do
 end do   

! 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(efg_g)
 ABI_DEALLOCATE(efg_r)
 ABI_DEALLOCATE(xcart)

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

 end subroutine make_efg_ion
!!***
