!{\src2tex{textfont=tt}}
!!****f* ABINIT/hartrestr
!! NAME
!! hartrestr
!!
!! FUNCTION
!! To be called for strain perturbation only
!! Compute the inhomogenous terms generated by the strain derivative of
!! Hartree potential due to the ground state charge rho(G)
!!
!!  FFT of (rho(G)/pi)*[d(1/G**2)/d(strain)
!!
!!          - delta(diagonal strain)*(1/G**2)]
!!
!! COPYRIGHT
!! Copyright (C) 1998-2014 ABINIT group (DRH, DCA, XG, GMR).
!! 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 .
!!
!! NOTES
!! *based largely on hartre.f
!! *Modified code to avoid if statements inside loops to skip G=0.
!!  Replaced if statement on G^2>gsqcut to skip G s outside where
!!  rho(G) should be 0.  Effect is negligible but gsqcut should be
!!  used to be strictly consistent with usage elsewhere in code.
!! *The speed-up is provided by doing a few precomputations outside
!!  the inner loop. One variable size array is needed for this (gq).
!!
!! INPUTS
!!  gmet(3,3)=metrix tensor in G space in Bohr**-2.
!!  gprimd(3,3)=reciprocal space dimensional primitive translations
!!  gsqcut=cutoff value on G**2 for sphere inside fft box.
!!  idir=direction of the current perturbation
!!  ipert=type of the perturbation
!!  mpi_enreg=informations about MPI parallelization
!!  natom=number of atoms in cell.
!!  nfft=number of fft grid points (gsqcut=(boxcut**2)*ecut/(2._dp*(Pi**2))
!!  ngfft(18)=contain all needed information about 3D FFT,
!!     see ~abinit/doc/input_variables/vargs.htm#ngfft
!!  rhog(2,nfft)=array for Fourier transform of GS electron density
!!
!! OUTPUT
!!  vhartr1(nfft)=Inhomogeneous term in strain-perturbation-induced Hartree
!!   potential in real space,
!!
!! PARENTS
!!      nselt3,rhotov3
!!
!! CHILDREN
!!      fourdp
!!
!! SOURCE

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

#include "abi_common.h"


subroutine hartrestr(gmet,gprimd,gsqcut,idir,ipert,mpi_enreg,natom,nfft,ngfft,&
&  paral_kgb,rhog,vhartr1)

 use defs_basis
 use defs_abitypes
 use m_errors
 use m_profiling

!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 'hartrestr'
 use interfaces_53_ffts
!End of the abilint section

 implicit none

!Arguments ------------------------------------
!scalars
 integer,intent(in) :: idir,ipert,natom,nfft,paral_kgb
 real(dp),intent(in) :: gsqcut
 type(MPI_type),intent(in) :: mpi_enreg
!arrays
 integer,intent(in) :: ngfft(18)
 real(dp),intent(in) :: gmet(3,3),gprimd(3,3),rhog(2,nfft)
 real(dp),intent(out) :: vhartr1(nfft)

!Local variables-------------------------------
!scalars
 integer,parameter :: im=2,re=1
 integer :: i1,i2,i23,i3,id2,id3,ig,ig2,ig3,ii,ii1,ing,istr,ka,kb,n1,n2,n3
 real(dp),parameter :: tolfix=1.000000001_dp
 real(dp) :: cutoff,ddends,den,dgsds,gqg2p3,gqgm12,gqgm13,gqgm23,gs,gs2,gs3
 real(dp) :: term
 character(len=500) :: message
!arrays
 integer,save :: idx(12)=(/1,1,2,2,3,3,3,2,3,1,2,1/)
 integer :: id(3)
 integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:)
 real(dp) :: dgmetds(3,3),gqr(3)
 real(dp),allocatable :: gq(:,:),work1(:,:)

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

 if( .not. (ipert==natom+3 .or. ipert==natom+4))then
   write(message, '(a,i0,a,a)' )&
&   '  From the calling routine, ipert=',ipert,ch10,&
&   '  so this routine for the strain perturbation should not be called.'
   MSG_BUG(message)
 end if

 n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)

!Get the distrib associated with this fft_grid
 if (n2 == mpi_enreg%distribfft%n2_coarse) then
   fftn2_distrib => mpi_enreg%distribfft%tab_fftdp2_distrib
   ffti2_local => mpi_enreg%distribfft%tab_fftdp2_local
 else if (n2 == mpi_enreg%distribfft%n2_fine) then
   fftn2_distrib => mpi_enreg%distribfft%tab_fftdp2dg_distrib
   ffti2_local => mpi_enreg%distribfft%tab_fftdp2dg_local
 else
   MSG_BUG("Unable to find an allocated distrib for this fft grid")
 end if

!Initialize a few quantities
 cutoff=gsqcut*tolfix

 istr=idir + 3*(ipert-natom-3)

 if(istr<1 .or. istr>6)then
   write(message, '(a,i10,a,a,a)' )&
&   '  Input dir gives istr=',istr,' not allowed.',ch10,&
&   '  Possible values are 1,2,3,4,5,6 only.'
   MSG_BUG(message)
 end if

 ka=idx(2*istr-1);kb=idx(2*istr)
 do ii = 1,3
   dgmetds(:,ii)=-(gprimd(ka,:)*gprimd(kb,ii)+gprimd(kb,:)*gprimd(ka,ii))
 end do
!For historical reasons:
 dgmetds(:,:)=0.5_dp*dgmetds(:,:)

!In order to speed the routine, precompute the components of g+q
!Also check if the booked space was large enough...
 ABI_ALLOCATE(gq,(3,max(n1,n2,n3)))
 do ii=1,3
   id(ii)=ngfft(ii)/2+2
   do ing=1,ngfft(ii)
     ig=ing-(ing/id(ii))*ngfft(ii)-1
     gq(ii,ing)=ig
   end do
 end do

 ABI_ALLOCATE(work1,(2,nfft))
 id2=n2/2+2
 id3=n3/2+2
!Triple loop on each dimension
 do i3=1,n3
   ig3=i3-(i3/id3)*n3-1
!  Precompute some products that do not depend on i2 and i1
   gqr(3)=gq(3,i3)
   gs3=gq(3,i3)*gq(3,i3)*gmet(3,3)
   gqgm23=gq(3,i3)*gmet(2,3)*2
   gqgm13=gq(3,i3)*gmet(1,3)*2

   do i2=1,n2
     if (fftn2_distrib(i2)==mpi_enreg%me_fft) then
       gqr(2)=gq(2,i2)
       gs2=gs3+ gq(2,i2)*(gq(2,i2)*gmet(2,2)+gqgm23)
       gqgm12=gq(2,i2)*gmet(1,2)*2
       gqg2p3=gqgm13+gqgm12
       ig2=i2-(i2/id2)*n2-1
!      i23=n1*((i2-1)+n2*(i3-1))
       i23=n1*((ffti2_local(i2)-1)+(n2/mpi_enreg%nproc_fft)*(i3-1))
!      Do the test that eliminates the Gamma point outside
!      of the inner loop
       ii1=1
       if(i23==0  .and. ig2==0 .and. ig3==0)then
         ii1=2
         work1(re,1+i23)=0.0_dp
         work1(im,1+i23)=0.0_dp
       end if

!      Final inner loop on the first dimension
!      (note the lower limit)
       do i1=ii1,n1
         gs=gs2+ gq(1,i1)*(gq(1,i1)*gmet(1,1)+gqg2p3)
         ii=i1+i23
         if(gs<=cutoff)then
           den=piinv/gs
           gqr(1)=gq(1,i1)
           dgsds=&
&           (gqr(1)*(dgmetds(1,1)*gqr(1)+dgmetds(1,2)*gqr(2)+dgmetds(1,3)*gqr(3))+  &
&           gqr(2)*(dgmetds(2,1)*gqr(1)+dgmetds(2,2)*gqr(2)+dgmetds(2,3)*gqr(3))+  &
&           gqr(3)*(dgmetds(3,1)*gqr(1)+dgmetds(3,2)*gqr(2)+dgmetds(3,3)*gqr(3)) )
           ddends=-piinv*dgsds/gs**2
           if(istr<=3)then
             term=2.0_dp*ddends-den
           else
             term=2.0_dp*ddends
           end if
           work1(re,ii)=rhog(re,ii)*term
           work1(im,ii)=rhog(im,ii)*term
         else
           work1(re,ii)=0.0_dp
           work1(im,ii)=0.0_dp
         end if
!        End loop on i1
       end do
     end if
!    End loop on i2
   end do

!  End loop on i3
 end do

 ABI_DEALLOCATE(gq)

!Fourier Transform Vhartree.
!Vh in reciprocal space was stored in work1
 call fourdp(1,work1,vhartr1,1,mpi_enreg,nfft,ngfft,paral_kgb,0)

 ABI_DEALLOCATE(work1)

end subroutine hartrestr
!!***
