TABLE OF CONTENTS
ABINIT/m_integrals [ 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 ]
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