TABLE OF CONTENTS


ABINIT/hartrestr [ Functions ]

[ Top ] [ Functions ]

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-2018 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

  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/variables/vargs.htm#ngfft
  rhog(2,nfft)=array for Fourier transform of GS electron density
  rprimd(3,3)=dimensional primitive translations in real space (bohr)

OUTPUT

  vhartr1(nfft)=Inhomogeneous term in strain-perturbation-induced Hartree
   potential in real space,

PARENTS

      dfpt_nselt,dfpt_nstpaw,dfpt_rhotov

CHILDREN

      fourdp,metric,ptabs_fourdp

SOURCE

 55 #if defined HAVE_CONFIG_H
 56 #include "config.h"
 57 #endif
 58 
 59 #include "abi_common.h"
 60 
 61 
 62 subroutine hartrestr(gsqcut,idir,ipert,mpi_enreg,natom,nfft,ngfft,&
 63 &  paral_kgb,rhog,rprimd,vhartr1)
 64 
 65  use defs_basis
 66  use defs_abitypes
 67  use m_errors
 68  use m_profiling_abi
 69 
 70  use m_mpinfo,     only : ptabs_fourdp
 71  use m_geometry,     only : metric
 72 
 73 !This section has been created automatically by the script Abilint (TD).
 74 !Do not modify the following lines by hand.
 75 #undef ABI_FUNC
 76 #define ABI_FUNC 'hartrestr'
 77  use interfaces_53_ffts
 78 !End of the abilint section
 79 
 80  implicit none
 81 
 82 !Arguments ------------------------------------
 83 !scalars
 84  integer,intent(in) :: idir,ipert,natom,nfft,paral_kgb
 85  real(dp),intent(in) :: gsqcut
 86  type(MPI_type),intent(in) :: mpi_enreg
 87 !arrays
 88  integer,intent(in) :: ngfft(18)
 89  real(dp),intent(in) :: rhog(2,nfft),rprimd(3,3)
 90  real(dp),intent(out) :: vhartr1(nfft)
 91 
 92 !Local variables-------------------------------
 93 !scalars
 94  integer,parameter :: im=2,re=1
 95  integer :: i1,i2,i23,i3,id2,id3,ig,ig2,ig3,ii,ii1,ing,istr,ka,kb,n1,n2,n3
 96  real(dp),parameter :: tolfix=1.000000001_dp
 97  real(dp) :: cutoff,ddends,den,dgsds,gqg2p3,gqgm12,gqgm13,gqgm23,gs,gs2,gs3
 98  real(dp) :: term,ucvol
 99  character(len=500) :: message
100 !arrays
101  integer,save :: idx(12)=(/1,1,2,2,3,3,3,2,3,1,2,1/)
102  integer :: id(3)
103  integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:)
104  integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:)
105  real(dp) :: dgmetds(3,3),gmet(3,3),gprimd(3,3),gqr(3),rmet(3,3)
106  real(dp),allocatable :: gq(:,:),work1(:,:)
107 
108 ! *************************************************************************
109 
110  if( .not. (ipert==natom+3 .or. ipert==natom+4))then
111    write(message, '(a,i0,a,a)' )&
112 &   'From the calling routine, ipert=',ipert,ch10,&
113 &   'so this routine for the strain perturbation should not be called.'
114    MSG_BUG(message)
115  end if
116 
117  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
118 
119  n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
120 
121 !Get the distrib associated with this fft_grid
122  call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
123 
124 !Initialize a few quantities
125  cutoff=gsqcut*tolfix
126 
127  istr=idir + 3*(ipert-natom-3)
128 
129  if(istr<1 .or. istr>6)then
130    write(message, '(a,i10,a,a,a)' )&
131 &   'Input dir gives istr=',istr,' not allowed.',ch10,&
132 &   'Possible values are 1,2,3,4,5,6 only.'
133    MSG_BUG(message)
134  end if
135 
136  ka=idx(2*istr-1);kb=idx(2*istr)
137  do ii = 1,3
138    dgmetds(:,ii)=-(gprimd(ka,:)*gprimd(kb,ii)+gprimd(kb,:)*gprimd(ka,ii))
139  end do
140 !For historical reasons:
141  dgmetds(:,:)=0.5_dp*dgmetds(:,:)
142 
143 !In order to speed the routine, precompute the components of g+q
144 !Also check if the booked space was large enough...
145  ABI_ALLOCATE(gq,(3,max(n1,n2,n3)))
146  do ii=1,3
147    id(ii)=ngfft(ii)/2+2
148    do ing=1,ngfft(ii)
149      ig=ing-(ing/id(ii))*ngfft(ii)-1
150      gq(ii,ing)=ig
151    end do
152  end do
153 
154  ABI_ALLOCATE(work1,(2,nfft))
155  id2=n2/2+2
156  id3=n3/2+2
157 !Triple loop on each dimension
158  do i3=1,n3
159    ig3=i3-(i3/id3)*n3-1
160 !  Precompute some products that do not depend on i2 and i1
161    gqr(3)=gq(3,i3)
162    gs3=gq(3,i3)*gq(3,i3)*gmet(3,3)
163    gqgm23=gq(3,i3)*gmet(2,3)*2
164    gqgm13=gq(3,i3)*gmet(1,3)*2
165 
166    do i2=1,n2
167      if (fftn2_distrib(i2)==mpi_enreg%me_fft) then
168        gqr(2)=gq(2,i2)
169        gs2=gs3+ gq(2,i2)*(gq(2,i2)*gmet(2,2)+gqgm23)
170        gqgm12=gq(2,i2)*gmet(1,2)*2
171        gqg2p3=gqgm13+gqgm12
172        ig2=i2-(i2/id2)*n2-1
173 !      i23=n1*((i2-1)+n2*(i3-1))
174        i23=n1*((ffti2_local(i2)-1)+(n2/mpi_enreg%nproc_fft)*(i3-1))
175 !      Do the test that eliminates the Gamma point outside
176 !      of the inner loop
177        ii1=1
178        if(i23==0  .and. ig2==0 .and. ig3==0)then
179          ii1=2
180          work1(re,1+i23)=0.0_dp
181          work1(im,1+i23)=0.0_dp
182        end if
183 
184 !      Final inner loop on the first dimension
185 !      (note the lower limit)
186        do i1=ii1,n1
187          gs=gs2+ gq(1,i1)*(gq(1,i1)*gmet(1,1)+gqg2p3)
188          ii=i1+i23
189          if(gs<=cutoff)then
190            den=piinv/gs
191            gqr(1)=gq(1,i1)
192            dgsds=&
193 &           (gqr(1)*(dgmetds(1,1)*gqr(1)+dgmetds(1,2)*gqr(2)+dgmetds(1,3)*gqr(3))+  &
194 &           gqr(2)*(dgmetds(2,1)*gqr(1)+dgmetds(2,2)*gqr(2)+dgmetds(2,3)*gqr(3))+  &
195 &           gqr(3)*(dgmetds(3,1)*gqr(1)+dgmetds(3,2)*gqr(2)+dgmetds(3,3)*gqr(3)) )
196            ddends=-piinv*dgsds/gs**2
197            if(istr<=3)then
198              term=2.0_dp*ddends-den
199            else
200              term=2.0_dp*ddends
201            end if
202            work1(re,ii)=rhog(re,ii)*term
203            work1(im,ii)=rhog(im,ii)*term
204          else
205            work1(re,ii)=0.0_dp
206            work1(im,ii)=0.0_dp
207          end if
208 
209        end do ! End loop on i1
210      end if
211    end do ! End loop on i2
212  end do !  End loop on i3
213 
214  ABI_DEALLOCATE(gq)
215 
216 !Fourier Transform Vhartree.
217 !Vh in reciprocal space was stored in work1
218  call fourdp(1,work1,vhartr1,1,mpi_enreg,nfft,ngfft,paral_kgb,0)
219 
220  ABI_DEALLOCATE(work1)
221 
222 end subroutine hartrestr