./67_common/hartre1.F90

ABINIT/hartre1

ABINIT/hartre1 [ Functions ]

NAME

 hartre1

FUNCTION

 Compute the Hartree energy from the density in reciprocal space.

 Copyright (C) 1998-2017 ABINIT group (DCA, MF, 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

 The Hartree energy is given for the Coulomb interaction with a
 real space cutoff, $\theta(r-R)/r$, as
  ehvalues(3)=$\displaystyle \frac{1}{2\pi}\sum_{\vec G}
          \frac{|n(\vec G)|^2}{G^2}\left( 1-cos(GR) \right)$,
 and for the true Coulomb interaction, $1/r$, as
  ehvalues(2)=$\displaystyle \frac{1}{2\pi}\sum_{\vec G \neq 0}\frac{|n(\vec G)|^2}{G^2}$.
 Also give
  ehvalues(1)=$\displaystyle \frac{1}{2\pi}\sum_{\vec G \neq 0}
          \frac{|n(\vec G)|^2}{G^2}\left( 1-cos(GR) \right)$.

INPUTS

  cplex= if 1, vhartr is REAL, if 2, vhartr is COMPLEX
  gmet(3,3)=metrix tensor in G space in Bohr**-2.
  gsqcut=cutoff value on G**2 for sphere inside fft box.
     (gsqcut=(boxcut**2)*ecut/(2._dp*(Pi**2))
  nfft=(effective) number of FFT grid points (for this processor)
  ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/input_variables/vargs.htm#ngfft
  qphon(3)=reduced coordinates for the phonon wavelength (needed if cplex==2).
  rcut=cutoff radius for Coulomb interaction in bohr
  rhog(2,nfft)=electron density in G space
  ucvol=unit cell volume

OUTPUT

  ehvalues(1)=Hartree energy with cutoff interaction without $G=0$ term
  ehvalues(2)=Hartree energy with full interaction without $G=0$ term
  ehvalues(2)=Hartree energy with cutoff interaction including $G=0$ term
  vhart(cplex*nfft)=Hartree potential in real space, either REAL or COMPLEX

WARNINGS

 Case cplex=2 not implemented.
 Hartree potential is not computed.

TODO

 Implement case cplex=2, clean up and calculate Hartree potential.

PARENTS

CHILDREN

      fourdp,timab

SOURCE

 59 #if defined HAVE_CONFIG_H
 60 #include "config.h"
 61 #endif
 62 
 63 #include "abi_common.h"
 64 
 65 subroutine hartre1(cplex,gmet,gsqcut,nfft,ngfft,paral_kgb,qphon,rhog,vhartr,ehvalues,rcut,ucvol)
 66 
 67  use defs_basis
 68  use defs_abitypes
 69  use m_errors
 70  use m_profiling_abi
 71 
 72 !This section has been created automatically by the script Abilint (TD).
 73 !Do not modify the following lines by hand.
 74 #undef ABI_FUNC
 75 #define ABI_FUNC 'hartre1'
 76  use interfaces_18_timing
 77  use interfaces_53_ffts
 78 !End of the abilint section
 79 
 80  implicit none
 81 
 82 !Arguments ------------------------------------
 83 !scalars
 84  integer,intent(in) :: cplex,nfft,paral_kgb
 85  real(dp),intent(in) :: gsqcut,rcut,ucvol
 86 !arrays
 87  integer,intent(in) :: ngfft(18)
 88  real(dp),intent(in) :: gmet(3,3),qphon(3),rhog(2,nfft)
 89  real(dp),intent(out) :: ehvalues(3),vhartr(cplex*nfft)
 90 
 91 !Local variables-------------------------------
 92 !scalars
 93  integer,parameter :: im=2,re=1
 94  integer :: i1,i2,i23,i3,ig,ii,ii1,ing,qeq0
 95  real(dp),parameter :: tolfix=1.000000001_dp
 96  real(dp) :: cutoff,den,gqg2p3,gqgm12,gqgm13,gqgm23,gs,gs2,gs3
 97  character(len=500) :: message
 98  type(MPI_type) :: mpi_enreg
 99 !arrays
100  integer :: id(3)
101  real(dp) :: tsec(2)
102  real(dp),allocatable :: gq(:,:),work1(:,:)
103 
104 ! *************************************************************************
105 
106 !Keep track of total time spent in hartre
107  call timab(10,1,tsec)
108 
109 !Check that cplex has an allowed value
110  if(cplex/=1 .and. cplex/=2)then
111    write(message, '(a,i0,a,a)' )&
112 &   'From the calling routine, cplex=',cplex,ch10,&
113 &   'but the only value allowed are 1 and 2.'
114    MSG_BUG(message)
115  end if
116 
117 !Initialize a few quantities
118  cutoff=gsqcut*tolfix
119 !This is to allow q=0
120  qeq0=0
121  if(qphon(1)**2+qphon(2)**2+qphon(3)**2<1.d-15) qeq0=1
122 
123 !If cplex=1 then qphon should be 0 0 0
124  if (cplex==1.and. qeq0/=1) then
125    write(message,'(a,3e12.4,a,a)')&
126 &   'cplex=1 but qphon=',qphon,ch10,&
127 &   'qphon should be 0 0 0.'
128    MSG_BUG(message)
129  end if
130 
131 !In order to speed the routine, precompute the components of g+q
132 !Also check if the booked space was large enough...
133  ABI_ALLOCATE(gq,(3,max(ngfft(1),ngfft(2),ngfft(3))))
134  do ii=1,3
135    id(ii)=ngfft(ii)/2+2
136    do ing=1,ngfft(ii)
137      ig=ing-(ing/id(ii))*ngfft(ii)-1
138      gq(ii,ing)=ig+qphon(ii)
139    end do
140  end do
141 
142  ABI_ALLOCATE(work1,(2,nfft))
143 
144 !Triple loop on each dimension
145  ehvalues(:)=zero
146  do i3=1,ngfft(3)
147 
148 !  Precompute some products that do not depend on i2 and i1
149    gs3=gq(3,i3)*gq(3,i3)*gmet(3,3)
150    gqgm23=gq(3,i3)*gmet(2,3)*2
151    gqgm13=gq(3,i3)*gmet(1,3)*2
152 
153    do i2=1,ngfft(2)
154      gs2=gs3+ gq(2,i2)*(gq(2,i2)*gmet(2,2)+gqgm23)
155      gqgm12=gq(2,i2)*gmet(1,2)*2
156      gqg2p3=gqgm13+gqgm12
157 
158      i23=ngfft(1)*((i2-1)+ngfft(2)*(i3-1))
159 !    Do the test that eliminates the Gamma point outside
160 !    of the inner loop
161      ii1=1
162      if(i23==0 .and. qeq0==1)then
163        ii1=2
164        work1(re,1+i23)=0.0_dp
165        work1(im,1+i23)=0.0_dp
166      end if
167 
168 !    Final inner loop on the first dimension
169 !    (note the lower limit)
170      do i1=ii1,ngfft(1)
171        gs=gs2+ gq(1,i1)*(gq(1,i1)*gmet(1,1)+gqg2p3)
172        ii=i1+i23
173        if(gs<=cutoff)then
174          den=piinv/gs
175          work1(re,ii)=rhog(re,ii)*den
176          work1(im,ii)=rhog(im,ii)*den
177 !        MF evaluate reciprocal space hartree energy
178          den=(rhog(re,ii)**2+rhog(im,ii)**2)/gs
179 !        MF use cutoff Coulomb interaction
180          ehvalues(re)=ehvalues(re)+den*(1._dp-cos(2._dp*pi*sqrt(gs)*rcut))
181 !        MF use full Coulomb interaction
182          ehvalues(im)=ehvalues(im)+den
183        else
184          work1(re,ii)=0.0_dp
185          work1(im,ii)=0.0_dp
186        end if
187 !      End loop on i1
188      end do
189 
190 !    End loop on i2
191    end do
192 
193 !  End loop on i3
194  end do
195  ehvalues(3)=ehvalues(re)+(rhog(re,1)**2+rhog(im,1)**2)*0.5_dp*(2*pi*rcut)**2
196 
197  ABI_DEALLOCATE(gq)
198 
199 !MF
200 !these are the correct pi factors,
201 !numerator:   $1/2$ [double counting] $\times  4\pi$ [Poisson eq] $= 2\pi$
202 !denominator: $2\pi$ [reciprocal lattice vectors] squared         $= (2\pi)^2$
203 !gives the same result as real space evaluation via
204 !$\frac{1}{2}\int n(\vec r)*V_{H}(\vec r) d^3r$
205  ehvalues(:)=ehvalues(:)*ucvol*2*pi/(2*pi)**2
206 !MF
207 
208 !DEBUG
209 !write(std_out,*)' hartre : before fourdp'
210 !write(std_out,*)' cplex,nfft,ngfft',cplex,nfft,ngfft
211 !write(std_out,*)' maxval work1=',maxval(abs(work1))
212 !ENDDEBUG
213 
214 !Fourier Transform Vhartree.
215 !Vh in reciprocal space was stored in work1
216 !MF no need to calculate potential
217  if(.false.)then
218    mpi_enreg%me_fft=0
219    mpi_enreg%nproc_fft=1
220    call fourdp(cplex,work1,vhartr,1,mpi_enreg,nfft,ngfft,paral_kgb,0)
221  end if
222 
223  ABI_DEALLOCATE(work1)
224 
225  call timab(10,2,tsec)
226 
227 end subroutine hartre1

TABLE OF CONTENTS

ABINIT/hartre1 [ Functions ]