TABLE OF CONTENTS

ABINIT/m_integrals [ Modules ]

[ Top ] [ Modules ]

NAME

  m_integrals

FUNCTION

  Helper functions to compute integrals

COPYRIGHT

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

20 #if defined HAVE_CONFIG_H
21 #include "config.h"
22 #endif
23 
24 #include "abi_common.h"
25 
26 module m_integrals
27 
28  use defs_basis
29  use m_errors
30  use m_abicore
31 
32  use m_numeric_tools,   only : simpson_int
33 
34  implicit none
35 
36  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

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