TABLE OF CONTENTS


ABINIT/m_integrals [ Modules ]

[ Top ] [ Modules ]

NAME

  m_integrals

FUNCTION

  Helper functions to compute integrals

COPYRIGHT

  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 .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 module m_integrals
28 
29  use defs_basis
30  use m_errors
31  use m_abicore
32 
33  use m_numeric_tools,   only : simpson_int
34 
35  implicit none
36 
37  private

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)

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

 78 subroutine radsintr(funr,funq,mqgrid,mrgrid,qgrid,rgrid,yq1,yqn)
 79 
 80 
 81 !This section has been created automatically by the script Abilint (TD).
 82 !Do not modify the following lines by hand.
 83 #undef ABI_FUNC
 84 #define ABI_FUNC 'radsintr'
 85 !End of the abilint section
 86 
 87  implicit none
 88 
 89 !Arguments ------------------------------------
 90 !scalars
 91  integer , intent(in)  :: mqgrid,mrgrid
 92  real(dp), intent(out) :: yq1,yqn
 93 !arrays
 94  real(dp), intent(in)  :: funr(mrgrid),qgrid(mqgrid),rgrid(mrgrid)
 95  real(dp), intent(out) :: funq(mqgrid)
 96 !Local variables-------------------------------
 97 !scalars
 98  integer :: iq,ir,irmax
 99  real(dp) :: arg,r0tor1,r1torm,rmax,rmtoin,rstep
100  logical :: begin_r0
101 !arrays
102  real(dp),allocatable :: ff(:),intg(:),rzf(:)
103 
104 ! *************************************************************************
105 
106 !DEBUG
107 !write (std_out,*) ' radsintr : enter'
108 !ENDDEBUG
109 
110  rstep=rgrid(2)-rgrid(1);rmax=min(20._dp,rgrid(mrgrid))
111  irmax=int(tol8+rmax/rstep)+1
112 
113 !Particular case of a null fuction to transform
114  if (maxval(abs(funr(1:irmax)))<=1.e-20_dp) then
115    funq=zero;yq1=zero;yqn=zero
116    return
117  end if
118 
119  ABI_ALLOCATE(ff,(mrgrid))
120  ABI_ALLOCATE(intg,(mrgrid))
121  ABI_ALLOCATE(rzf,(mrgrid))
122  ff=zero;rzf=zero
123 
124 !Is mesh beginning with r=0 ?
125  begin_r0=(rgrid(1)<1.e-20_dp)
126 
127 !Store r.F(r)
128  do ir=1,irmax
129    ff(ir)=rgrid(ir)*funr(ir) !ff is part of the integrand
130  end do
131 
132 !===========================================
133 !=== Compute v(q) for q=0 separately
134 !===========================================
135 
136 !Integral from 0 to r1 (only if r1<>0)
137  r0tor1=zero;if (.not.begin_r0) r0tor1=(funr(1)*third)*rgrid(1)**three
138 
139 !Integral from r1 to rmax
140  do ir=1,irmax
141    if (abs(ff(ir))>1.e-20_dp) then
142      rzf(ir)=rgrid(ir)*ff(ir) !sin(2\pi*q*r)/(2*\pi*q*r)-->1 for q=0
143    end if                     !so the integrand is 4*\pi*r^2*F(r)
144  end do
145  call simpson_int(mrgrid,rstep,rzf,intg)
146  r1torm=intg(mrgrid)
147 
148 !Integral from rmax to infinity
149 !This part is neglected... might be improved.
150  rmtoin=zero
151 
152 !Sum of the three parts
153 
154  funq(1)=four_pi*(r0tor1+r1torm+rmtoin)
155 
156 !===========================================
157 !=== Compute v(q) for other q''s
158 !===========================================
159 
160 !Loop over q values
161  do iq=2,mqgrid
162    arg=two_pi*qgrid(iq)
163 
164 !  Integral from 0 to r1 (only if r1<>0)
165    r0tor1=zero
166    if (.not.begin_r0) r0tor1=(funr(1)/arg) &
167 &   *(sin(arg*rgrid(1))/arg-rgrid(1)*cos(arg*rgrid(1)))
168 
169 !  Integral from r1 to rmax
170    rzf=zero
171    do ir=1,irmax
172      if (abs(ff(ir))>1.e-20_dp) rzf(ir)=sin(arg*rgrid(ir))*ff(ir)
173    end do
174    call simpson_int(mrgrid,rstep,rzf,intg)
175    r1torm=intg(mrgrid)
176 
177 !  Integral from rmax to infinity
178 !  This part is neglected... might be improved.
179    rmtoin=zero
180 
181 !  Store F(q)
182    funq(iq) = two/qgrid(iq)*(r0tor1+r1torm+rmtoin)
183  end do
184 
185 
186 !===========================================
187 !=== Compute derivatives of F(q)
188 !=== at ends of interval
189 !===========================================
190 
191 !yq(0)=zero
192  yq1=zero
193 
194 !yp(qmax)=$
195  arg=two_pi*qgrid(mqgrid)
196 
197 !Integral from 0 to r1 (only if r1<>0)
198  r0tor1=zero
199  if (.not.begin_r0) r0tor1=(funr(1)/arg) &
200 & *(sin(arg*rgrid(1))*(rgrid(1)**2-three/arg**2) &
201 & +three*rgrid(1)*cos(arg*rgrid(1))/arg)
202 
203 !Integral from r1 to rmax
204  do ir=1,irmax
205    if (abs(ff(ir))>1.e-20_dp) rzf(ir)=(rgrid(ir)*cos(arg*rgrid(ir)) &
206 &   -sin(arg*rgrid(ir))/arg)*ff(ir)
207  end do
208  call simpson_int(mrgrid,rstep,rzf,intg)
209  r1torm=intg(mrgrid)
210 
211 !Integral from rmax to infinity
212 !This part is neglected... might be improved.
213  rmtoin=zero
214 
215 !Sum of the three parts
216  yqn=(four*pi/qgrid(mqgrid))*(r0tor1+r1torm+rmtoin)
217 
218  ABI_DEALLOCATE(ff)
219  ABI_DEALLOCATE(intg)
220  ABI_DEALLOCATE(rzf)
221 
222 end subroutine radsintr