./32_util/radsintr.F90

ABINIT/radsintr

ABINIT/radsintr [ 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.

  Copyright (C) 2010-2018 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 .

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)

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

PARENTS

      m_xc_vdw,test_radsintr

CHILDREN

      simpson_int

SOURCE

 40 #if defined HAVE_CONFIG_H
 41 #include "config.h"
 42 #endif
 43 
 44 #include "abi_common.h"
 45 
 46 subroutine radsintr(funr,funq,mqgrid,mrgrid,qgrid,rgrid,yq1,yqn)
 47 
 48  use defs_basis
 49  use m_profiling_abi
 50 
 51  use m_numeric_tools,   only : simpson_int
 52 
 53 !This section has been created automatically by the script Abilint (TD).
 54 !Do not modify the following lines by hand.
 55 #undef ABI_FUNC
 56 #define ABI_FUNC 'radsintr'
 57 !End of the abilint section
 58 
 59  implicit none
 60 
 61 !Arguments ------------------------------------
 62 !scalars
 63  integer , intent(in)  :: mqgrid,mrgrid
 64  real(dp), intent(out) :: yq1,yqn
 65 !arrays
 66  real(dp), intent(in)  :: funr(mrgrid),qgrid(mqgrid),rgrid(mrgrid)
 67  real(dp), intent(out) :: funq(mqgrid)
 68 !Local variables-------------------------------
 69 !scalars
 70  integer :: iq,ir,irmax
 71  real(dp) :: arg,r0tor1,r1torm,rmax,rmtoin,rstep
 72  logical :: begin_r0
 73 !arrays
 74  real(dp),allocatable :: ff(:),intg(:),rzf(:)
 75 
 76 ! *************************************************************************
 77 
 78 !DEBUG
 79 !write (std_out,*) ' radsintr : enter'
 80 !ENDDEBUG
 81 
 82  rstep=rgrid(2)-rgrid(1);rmax=min(20._dp,rgrid(mrgrid))
 83  irmax=int(tol8+rmax/rstep)+1
 84 
 85 !Particular case of a null fuction to transform
 86  if (maxval(abs(funr(1:irmax)))<=1.e-20_dp) then
 87    funq=zero;yq1=zero;yqn=zero
 88    return
 89  end if
 90 
 91  ABI_ALLOCATE(ff,(mrgrid))
 92  ABI_ALLOCATE(intg,(mrgrid))
 93  ABI_ALLOCATE(rzf,(mrgrid))
 94  ff=zero;rzf=zero
 95 
 96 !Is mesh beginning with r=0 ?
 97  begin_r0=(rgrid(1)<1.e-20_dp)
 98 
 99 !Store r.F(r)
100  do ir=1,irmax
101    ff(ir)=rgrid(ir)*funr(ir) !ff is part of the integrand
102  end do
103 
104 !===========================================
105 !=== Compute v(q) for q=0 separately
106 !===========================================
107 
108 !Integral from 0 to r1 (only if r1<>0)
109  r0tor1=zero;if (.not.begin_r0) r0tor1=(funr(1)*third)*rgrid(1)**three
110 
111 !Integral from r1 to rmax
112  do ir=1,irmax
113    if (abs(ff(ir))>1.e-20_dp) then
114      rzf(ir)=rgrid(ir)*ff(ir) !sin(2\pi*q*r)/(2*\pi*q*r)-->1 for q=0
115    end if                     !so the integrand is 4*\pi*r^2*F(r)
116  end do
117  call simpson_int(mrgrid,rstep,rzf,intg)
118  r1torm=intg(mrgrid)
119 
120 !Integral from rmax to infinity
121 !This part is neglected... might be improved.
122  rmtoin=zero
123 
124 !Sum of the three parts
125 
126  funq(1)=four_pi*(r0tor1+r1torm+rmtoin)
127 
128 !===========================================
129 !=== Compute v(q) for other q''s
130 !===========================================
131 
132 !Loop over q values
133  do iq=2,mqgrid
134    arg=two_pi*qgrid(iq)
135 
136 !  Integral from 0 to r1 (only if r1<>0)
137    r0tor1=zero
138    if (.not.begin_r0) r0tor1=(funr(1)/arg) &
139 &   *(sin(arg*rgrid(1))/arg-rgrid(1)*cos(arg*rgrid(1)))
140 
141 !  Integral from r1 to rmax
142    rzf=zero
143    do ir=1,irmax
144      if (abs(ff(ir))>1.e-20_dp) rzf(ir)=sin(arg*rgrid(ir))*ff(ir)
145    end do
146    call simpson_int(mrgrid,rstep,rzf,intg)
147    r1torm=intg(mrgrid)
148 
149 !  Integral from rmax to infinity
150 !  This part is neglected... might be improved.
151    rmtoin=zero
152 
153 !  Store F(q)
154    funq(iq) = two/qgrid(iq)*(r0tor1+r1torm+rmtoin)
155  end do
156 
157 
158 !===========================================
159 !=== Compute derivatives of F(q)
160 !=== at ends of interval
161 !===========================================
162 
163 !yq(0)=zero
164  yq1=zero
165 
166 !yp(qmax)=$ 
167  arg=two_pi*qgrid(mqgrid)
168 
169 !Integral from 0 to r1 (only if r1<>0)
170  r0tor1=zero
171  if (.not.begin_r0) r0tor1=(funr(1)/arg) &
172 & *(sin(arg*rgrid(1))*(rgrid(1)**2-three/arg**2) &
173 & +three*rgrid(1)*cos(arg*rgrid(1))/arg)
174 
175 !Integral from r1 to rmax
176  do ir=1,irmax
177    if (abs(ff(ir))>1.e-20_dp) rzf(ir)=(rgrid(ir)*cos(arg*rgrid(ir)) &
178 &   -sin(arg*rgrid(ir))/arg)*ff(ir)
179  end do
180  call simpson_int(mrgrid,rstep,rzf,intg)
181  r1torm=intg(mrgrid)
182 
183 !Integral from rmax to infinity
184 !This part is neglected... might be improved.
185  rmtoin=zero
186 
187 !Sum of the three parts
188  yqn=(four*pi/qgrid(mqgrid))*(r0tor1+r1torm+rmtoin)
189 
190  ABI_DEALLOCATE(ff)
191  ABI_DEALLOCATE(intg)
192  ABI_DEALLOCATE(rzf)
193 
194 end subroutine radsintr

TABLE OF CONTENTS

ABINIT/radsintr [ Functions ]