TABLE OF CONTENTS


ABINIT/m_integrals [ Modules ]

[ Top ] [ Modules ]

NAME

  m_integrals

FUNCTION

  Helper functions to compute integrals

COPYRIGHT

  Copyright (C) 2010-2022 ABINIT group (Camilo Espejo)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

SOURCE

16 #if defined HAVE_CONFIG_H
17 #include "config.h"
18 #endif
19 
20 #include "abi_common.h"
21 
22 module m_integrals
23 
24  use defs_basis
25  use m_errors
26  use m_abicore
27 
28  use m_numeric_tools,   only : simpson_int
29  use m_macroave,        only : four1
30 
31  implicit none
32 
33 ! private
34 
35 
36  public :: radsintr

ABINIT/radsintr [ Functions ]

[ Top ] [ Functions ]

NAME

  radsintr

FUNCTION

 Computes the sine transformation of a radial function F(r) to V(q).
 Computes integrals using corrected Simpson integration on a linear grid.

INPUTS

  funr(mrgrid)=F(r) on radial grid
  mqgrid=number of grid points in q from 0 to qmax
  mrgrid=number of grid points in r from 0 to rmax
  qgrid(mqgrid)=q grid values (bohr**-1).
  rgrid(mrgrid)=r grid values (bohr)
  rfttype=Type of radial sin Fourier transform   

OUTPUT

  funq(mqgrid)=\int_0^inf 4\pi\frac{\sin(2\pi r)}{2\pi r}r^2F(r)dr
  yq1, yqn: d/dq (F(q)) at the ends of the interval

SIDE EFFECTS

NOTES

  This routine is a modified version of /src/65_psp/psp7lo.F90

SOURCE

 69 subroutine radsintr(funr,funq,mqgrid,mrgrid,qgrid,rgrid,yq1,yqn,rfttype)
 70 
 71 !Arguments ------------------------------------
 72 !scalars
 73  integer , intent(in)  :: mqgrid,mrgrid,rfttype
 74  real(dp), intent(out) :: yq1,yqn
 75 !arrays
 76  real(dp), intent(in)  :: funr(0:mrgrid),qgrid(0:mqgrid),rgrid(0:mrgrid)
 77  real(dp), intent(out) :: funq(0:mqgrid)
 78 !Local variables-------------------------------
 79 character(len=512) :: msg
 80 !scalars
 81  integer :: iq,ir,irmax,jr,nq
 82  real(dp) :: arg,cc,dq,dr,fr,qmax,qq,r0tor1,r1torm,rmax,rmtoin,rn,rr,rstep
 83  logical :: begin_r0
 84 !arrays
 85  real(dp),allocatable :: ff(:),intg(:),rzf(:)
 86  real(dp) :: fn(2,0:2*mrgrid),pp(2)
 87 
 88 ! *************************************************************************
 89 
 90  select case(rfttype)
 91 
 92  case(1)
 93    rstep=rgrid(2)-rgrid(1)
 94    rmax = rgrid(mrgrid)
 95   ! rmax=min(20._dp,rgrid(mrgrid))
 96    irmax=int(tol8+rmax/rstep)+1
 97   
 98   !Particular case of a null fuction to transform
 99    if (maxval(abs(funr(1:irmax)))<=1.e-20_dp) then
100      funq=zero;yq1=zero;yqn=zero
101      return
102    end if
103   
104    ABI_MALLOC(ff,(mrgrid))
105    ABI_MALLOC(intg,(mrgrid))
106    ABI_MALLOC(rzf,(mrgrid))
107    ff=zero;rzf=zero
108   
109   !Is mesh beginning with r=0 ?
110    begin_r0=(rgrid(1)<1.e-20_dp)
111   
112   !Store r.F(r)
113    do ir=1,irmax
114      ff(ir)=rgrid(ir)*funr(ir) !ff is part of the integrand
115    end do
116   
117   !===========================================
118   !=== Compute v(q) for q=0 separately
119   !===========================================
120   
121   !Integral from 0 to r1 (only if r1<>0)
122    r0tor1=zero;if (.not.begin_r0) r0tor1=(funr(1)*third)*rgrid(1)**three
123   
124   !Integral from r1 to rmax
125    r1torm = zero
126    do ir=1,irmax
127      if (abs(ff(ir))>1.e-20_dp) then
128        rzf(ir)=rgrid(ir)*ff(ir) !sin(2\pi*q*r)/(2*\pi*q*r)-->1 for q=0
129      end if                     !so the integrand is 4*\pi*r^2*F(r)
130      r1torm = r1torm + rzf(ir)
131    end do
132   ! call simpson_int(mrgrid,rstep,rzf,intg)
133   ! r1torm=intg(mrgrid)
134   ! r1torm = r1torm*rstep 
135   !Integral from rmax to infinity
136   !This part is neglected... might be improved.
137    rmtoin=zero
138   
139   !Sum of the three parts
140   
141    funq(1)=four_pi*(r0tor1+r1torm+rmtoin)
142   
143   !===========================================
144   !=== Compute v(q) for other q''s
145   !===========================================
146   
147   !Loop over q values
148    do iq=2,mqgrid
149      arg=two_pi*qgrid(iq)
150   
151   !  Integral from 0 to r1 (only if r1<>0)
152      r0tor1=zero
153      if (.not.begin_r0) r0tor1=(funr(1)/arg) &
154   &   *(sin(arg*rgrid(1))/arg-rgrid(1)*cos(arg*rgrid(1)))
155   
156   !  Integral from r1 to rmax
157      rzf=zero
158      r1torm = zero
159      do ir=1,irmax
160        if (abs(ff(ir))>1.e-20_dp) then 
161          rzf(ir)=sin(arg*rgrid(ir))*ff(ir)
162        end if
163        r1torm = r1torm + rzf(ir) 
164      end do
165   !   r1torm = r1torm*rstep 
166   !   call simpson_int(mrgrid,rstep,rzf,intg)
167   !   r1torm=intg(mrgrid)
168      
169   !  Integral from rmax to infinity
170   !  This part is neglected... might be improved.
171      rmtoin=zero
172   
173   !  Store F(q)
174      funq(iq) = two/qgrid(iq)*(r0tor1+r1torm+rmtoin)
175    end do
176   
177   
178   !===========================================
179   !=== Compute derivatives of F(q)
180   !=== at ends of interval
181   !===========================================
182   
183   !yq(0)=zero
184    yq1=zero
185   
186   !yp(qmax)=$
187    arg=two_pi*qgrid(mqgrid)
188   
189   !Integral from 0 to r1 (only if r1<>0)
190    r0tor1=zero
191    if (.not.begin_r0) r0tor1=(funr(1)/arg) &
192   & *(sin(arg*rgrid(1))*(rgrid(1)**2-three/arg**2) &
193   & +three*rgrid(1)*cos(arg*rgrid(1))/arg)
194   
195   !Integral from r1 to rmax
196    do ir=1,irmax
197      if (abs(ff(ir))>1.e-20_dp) rzf(ir)=(rgrid(ir)*cos(arg*rgrid(ir)) &
198   &   -sin(arg*rgrid(ir))/arg)*ff(ir)
199    end do
200    call simpson_int(mrgrid,rstep,rzf,intg)
201    r1torm=intg(mrgrid)
202   
203   !Integral from rmax to infinity
204   !This part is neglected... might be improved.
205    rmtoin=zero
206   
207   !Sum of the three parts
208    yqn=(four*pi/qgrid(mqgrid))*(r0tor1+r1torm+rmtoin)
209   
210    ABI_FREE(ff)
211    ABI_FREE(intg)
212    ABI_FREE(rzf)
213 
214  case(2)  
215 
216        rmax = rgrid(mrgrid)
217        nq = mrgrid  !NQ = NR
218        dr = rgrid(mrgrid) / mrgrid !RMAX / NR
219        dq = pi / rmax
220        qmax = nq * dq
221        cc = dr / sqrt(2 * pi)
222        pp(1) =  0.0_dp
223        pp(2) = -1.0_dp
224 ! C Initialize accumulation array -------------------------------------
225        do iq = 0,nq
226          funq(iq) = 0.0_dp !gg(iq) = 0.0_dp
227        end do
228 ! C -------------------------------------------------------------------
229 ! 
230 ! C Iterate on terms of the j_l(q*r) polynomial -----------------------
231 !       DO N = 0,L
232 ! 
233 ! C       Set up function to be fast fourier transformed
234          fn(1,0) = 0.0_dp
235          fn(2,0) = 0.0_dp
236          do jr = 1, 2*mrgrid-1
237  
238            if (jr .lt. mrgrid) then
239              ir = jr
240              rr = ir * dr
241              fr = funr(ir)
242            elseif (jr .eq. mrgrid) then
243              ir = jr
244              rr = ir * dr
245              fr = 0.0_dp
246            else
247              ir = 2*mrgrid - jr
248              rr = - (ir * dr)
249              fr = funr(ir) !* (-1.D0)**L
250            end if
251  
252 ! C         Find  r**2 * r**n / r**(l+1)
253            rn = rr !**(N-L+1)
254  
255            fn(1,jr) = cc * fr * rn * pp(1)
256            fn(2,jr) = cc * fr * rn * pp(2)
257          end do
258  
259 ! C       Perform one-dimensional complex FFT
260 ! !
261 ! !       Only the elements from 0 to 2*NR-1 of FN are used.
262 ! !       (a total of 2*NR). Four1 will receive a one-dimensional
263 ! !       array of size 2*NR.
264 ! !
265          call four1( fn, 2*mrgrid, +1 )
266 ! 
267 ! C       Accumulate contribution
268          do iq = 1,nq
269            qq = iq * dq
270            funq(iq) = ( funq(iq) + fn(1,iq) ) / qq
271          end do
272 ! 
273 !       ENDDO
274 ! C -------------------------------------------------------------------
275 ! 
276 ! C Special case for Q=0 ---------------------------------------------
277          funq(0) = 0.0_dp
278         !gg(0) = 0.0_dp
279      !  if ( L .EQ. 0 ) THEN
280          do ir = 1,mrgrid
281            rr = ir * dr
282            funq(0) = funq(0) + rr*rr * funr(ir)
283          end do
284          funq(0) = funq(0) * 2.0_dp * cc
285      !  ENDIF
286  case default 
287    
288    write(msg,'(2x,a,1x,i2)') " Unknown radial sine FFT type", rfttype
289    ABI_ERROR(msg)
290 
291  end select
292 
293 end subroutine radsintr
294 
295 end module m_integrals