TABLE OF CONTENTS


ABINIT/hartre [ Functions ]

[ Top ] [ Functions ]

NAME

 hartre

FUNCTION

 Given rho(G), compute Hartree potential (=FFT of rho(G)/pi/(G+q)**2)
 When cplex=1, assume q=(0 0 0), and vhartr will be REAL
 When cplex=2, q must be taken into account, and vhartr will be COMPLEX

COPYRIGHT

 Copyright (C) 1998-2018 ABINIT group (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

 *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

  cplex= if 1, vhartr is REAL, if 2, vhartr is COMPLEX
  gsqcut=cutoff value on G**2 for sphere inside fft box.
         (gsqcut=(boxcut**2)*ecut/(2.d0*(Pi**2))
  izero=if 1, unbalanced components of Vhartree(g) are set to zero
  mpi_enreg=information about MPI parallelization
  nfft=(effective) number of FFT grid points (for this processor)
  ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
  [qpt(3)=reduced coordinates for a wavevector to be combined with the G vectors (needed if cplex==2).]
  rhog(2,nfft)=electron density in G space
  rprimd(3,3)=dimensional primitive translations in real space (bohr)
  divgq0= [optional argument] value of the integration of the Coulomb singularity 4pi\int_BZ 1/q^2 dq

OUTPUT

  vhartr(cplex*nfft)=Hartree potential in real space, either REAL or COMPLEX

PARENTS

      dfpt_rhotov,dfptnl_loop,energy,fock_getghc,m_kxc,nonlinear,nres2vres
      odamix,prcref,prcref_PMA,respfn,rhotov,setup_positron,setvtr,tddft

CHILDREN

      fourdp,metric,ptabs_fourdp,timab,zerosym

SOURCE

 51 #if defined HAVE_CONFIG_H
 52 #include "config.h"
 53 #endif
 54 
 55 #include "abi_common.h"
 56 
 57 subroutine hartre(cplex,gsqcut,izero,mpi_enreg,nfft,ngfft,paral_kgb,rhog,rprimd,vhartr,&
 58 &  divgq0,qpt) ! Optional argument
 59 
 60  use defs_basis
 61  use defs_abitypes
 62  use m_errors
 63  use m_profiling_abi
 64 
 65  use m_geometry, only : metric
 66  use m_mpinfo,   only : ptabs_fourdp
 67  use m_fft,      only : zerosym
 68 
 69 !This section has been created automatically by the script Abilint (TD).
 70 !Do not modify the following lines by hand.
 71 #undef ABI_FUNC
 72 #define ABI_FUNC 'hartre'
 73  use interfaces_18_timing
 74  use interfaces_53_ffts
 75 !End of the abilint section
 76 
 77  implicit none
 78 
 79 !Arguments ------------------------------------
 80 !scalars
 81  integer,intent(in) :: cplex,izero,nfft,paral_kgb
 82  real(dp),intent(in) :: gsqcut
 83  type(MPI_type),intent(in) :: mpi_enreg
 84 !arrays
 85  integer,intent(in) :: ngfft(18)
 86  real(dp),intent(in) :: rprimd(3,3),rhog(2,nfft)
 87  real(dp),intent(in),optional :: divgq0
 88  real(dp),intent(in),optional :: qpt(3)
 89  real(dp),intent(out) :: vhartr(cplex*nfft)
 90 
 91 !Local variables-------------------------------
 92 !scalars
 93  integer,parameter :: im=2,re=1
 94  integer :: i1,i2,i23,i2_local,i3,id1,id2,id3
 95  integer :: ig,ig1min,ig1,ig1max,ig2,ig2min,ig2max,ig3,ig3min,ig3max
 96  integer :: ii,ii1,ing,n1,n2,n3,qeq0,qeq05,me_fft,nproc_fft
 97  real(dp),parameter :: tolfix=1.000000001e0_dp
 98  real(dp) :: cutoff,den,gqg2p3,gqgm12,gqgm13,gqgm23,gs,gs2,gs3,ucvol
 99  character(len=500) :: message
100 !arrays
101  integer :: id(3)
102  integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:)
103  integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:)
104  real(dp) :: gmet(3,3),gprimd(3,3),qpt_(3),rmet(3,3),tsec(2)
105  real(dp),allocatable :: gq(:,:),work1(:,:)
106 
107 ! *************************************************************************
108 
109  ! Keep track of total time spent in hartre
110  call timab(10,1,tsec)
111 
112  ! Check that cplex has an allowed value
113  if(cplex/=1 .and. cplex/=2)then
114    write(message, '(a,i0,a,a)' )&
115 &   'From the calling routine, cplex=',cplex,ch10,&
116 &   'but the only value allowed are 1 and 2.'
117    MSG_BUG(message)
118  end if
119 
120  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
121 
122  n1=ngfft(1); n2=ngfft(2); n3=ngfft(3)
123  nproc_fft = mpi_enreg%nproc_fft; me_fft = mpi_enreg%me_fft
124 
125  ! Get the distrib associated with this fft_grid
126  call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
127 
128  ! Initialize a few quantities
129  cutoff=gsqcut*tolfix
130  if(present(qpt))then
131    qpt_=qpt
132  else
133    qpt_=zero
134  end if
135  qeq0=0
136  if(qpt_(1)**2+qpt_(2)**2+qpt_(3)**2<1.d-15) qeq0=1
137  qeq05=0
138  if (qeq0==0) then
139    if (abs(abs(qpt_(1))-half)<tol12.or.abs(abs(qpt_(2))-half)<tol12.or. &
140 &   abs(abs(qpt_(3))-half)<tol12) qeq05=1
141  end if
142 
143  ! If cplex=1 then qpt_ should be 0 0 0
144  if (cplex==1.and. qeq0/=1) then
145    write(message,'(a,3e12.4,a,a)')&
146 &   'cplex=1 but qpt=',qpt_,ch10,&
147 &   'qpt should be 0 0 0.'
148    MSG_BUG(message)
149  end if
150 
151  ! If FFT parallelism then qpt should not be 1/2
152  if (nproc_fft>1.and.qeq05==1) then
153    write(message, '(a,3e12.4,a,a)' )&
154 &   'FFT parallelism selected but qpt',qpt_,ch10,&
155 &   'qpt(i) should not be 1/2...'
156    MSG_ERROR(message)
157  end if
158 
159  ! In order to speed the routine, precompute the components of g+q
160  ! Also check if the booked space was large enough...
161  ABI_ALLOCATE(gq,(3,max(n1,n2,n3)))
162  do ii=1,3
163    id(ii)=ngfft(ii)/2+2
164    do ing=1,ngfft(ii)
165      ig=ing-(ing/id(ii))*ngfft(ii)-1
166      gq(ii,ing)=ig+qpt_(ii)
167    end do
168  end do
169  ig1max=-1;ig2max=-1;ig3max=-1
170  ig1min=n1;ig2min=n2;ig3min=n3
171 
172  ABI_ALLOCATE(work1,(2,nfft))
173  id1=n1/2+2;id2=n2/2+2;id3=n3/2+2
174 
175  ! Triple loop on each dimension
176  do i3=1,n3
177    ig3=i3-(i3/id3)*n3-1
178    ! Precompute some products that do not depend on i2 and i1
179    gs3=gq(3,i3)*gq(3,i3)*gmet(3,3)
180    gqgm23=gq(3,i3)*gmet(2,3)*2
181    gqgm13=gq(3,i3)*gmet(1,3)*2
182 
183    do i2=1,n2
184      ig2=i2-(i2/id2)*n2-1
185      if (fftn2_distrib(i2) == me_fft) then
186        gs2=gs3+ gq(2,i2)*(gq(2,i2)*gmet(2,2)+gqgm23)
187        gqgm12=gq(2,i2)*gmet(1,2)*2
188        gqg2p3=gqgm13+gqgm12
189 
190        i2_local = ffti2_local(i2)
191        i23=n1*(i2_local-1 +(n2/nproc_fft)*(i3-1))
192        ! Do the test that eliminates the Gamma point outside of the inner loop
193        ii1=1
194        if(i23==0 .and. qeq0==1  .and. ig2==0 .and. ig3==0)then
195          ii1=2
196          work1(re,1+i23)=zero
197          work1(im,1+i23)=zero
198          ! If the value of the integration of the Coulomb singularity 4pi\int_BZ 1/q^2 dq is given, use it
199          if (PRESENT(divgq0)) then
200            work1(re,1+i23)=rhog(re,1+i23)*divgq0*piinv
201            work1(im,1+i23)=rhog(im,1+i23)*divgq0*piinv
202          end if
203        end if
204 
205        ! Final inner loop on the first dimension (note the lower limit)
206        do i1=ii1,n1
207          gs=gs2+ gq(1,i1)*(gq(1,i1)*gmet(1,1)+gqg2p3)
208          ii=i1+i23
209 
210          if(gs<=cutoff)then
211            ! Identify min/max indexes (to cancel unbalanced contributions later)
212            ! Count (q+g)-vectors with similar norm
213            if ((qeq05==1).and.(izero==1)) then
214              ig1=i1-(i1/id1)*n1-1
215              ig1max=max(ig1max,ig1); ig1min=min(ig1min,ig1)
216              ig2max=max(ig2max,ig2); ig2min=min(ig2min,ig2)
217              ig3max=max(ig3max,ig3); ig3min=min(ig3min,ig3)
218            end if
219 
220            den=piinv/gs
221            work1(re,ii)=rhog(re,ii)*den
222            work1(im,ii)=rhog(im,ii)*den
223          else
224            ! gs>cutoff
225            work1(re,ii)=zero
226            work1(im,ii)=zero
227          end if
228 
229        end do ! End loop on i1
230      end if
231    end do ! End loop on i2
232  end do ! End loop on i3
233 
234  ABI_DEALLOCATE(gq)
235 
236  if (izero==1) then
237    ! Set contribution of unbalanced components to zero
238 
239    if (qeq0==1) then !q=0
240      call zerosym(work1,2,n1,n2,n3,comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft)
241 
242    else if (qeq05==1) then
243      !q=1/2; this doesn't work in parallel
244      ig1=-1;if (mod(n1,2)==0) ig1=1+n1/2
245      ig2=-1;if (mod(n2,2)==0) ig2=1+n2/2
246      ig3=-1;if (mod(n3,2)==0) ig3=1+n3/2
247      if (abs(abs(qpt_(1))-half)<tol12) then
248        if (abs(ig1min)<abs(ig1max)) ig1=abs(ig1max)
249        if (abs(ig1min)>abs(ig1max)) ig1=n1-abs(ig1min)
250      end if
251      if (abs(abs(qpt_(2))-half)<tol12) then
252        if (abs(ig2min)<abs(ig2max)) ig2=abs(ig2max)
253        if (abs(ig2min)>abs(ig2max)) ig2=n2-abs(ig2min)
254      end if
255      if (abs(abs(qpt_(3))-half)<tol12) then
256        if (abs(ig3min)<abs(ig3max)) ig3=abs(ig3max)
257        if (abs(ig3min)>abs(ig3max)) ig3=n3-abs(ig3min)
258      end if
259      call zerosym(work1,2,n1,n2,n3,ig1=ig1,ig2=ig2,ig3=ig3,&
260 &     comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft)
261    end if
262  end if
263 
264  ! Fourier Transform Vhartree. Vh in reciprocal space was stored in work1
265  call fourdp(cplex,work1,vhartr,1,mpi_enreg,nfft,ngfft,paral_kgb,0)
266 
267  ABI_DEALLOCATE(work1)
268 
269  call timab(10,2,tsec)
270 
271 end subroutine hartre