## ABINIT/m_integrals [ Modules ]

[ Top ] [ Modules ]

NAME

  m_integrals


FUNCTION

  Helper functions to compute integrals


  Copyright (C) 2010-2018 ABINIT group (Camilo Espejo)
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


[ 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
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