LCPQ/eplf
10
0
Fork 0
mirror of https://gitlab.com/scemama/eplf synced 2024-12-22 04:14:17 +01:00
Code Issues Releases Wiki Activity

Optimized eplf integrals

Browse Source
This commit is contained in:
Anthony Scemama 2012-05-31 23:26:30 +02:00
parent 8a1c7a03a8
commit 40ad59efa8
2 changed files with 66 additions and 20 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

24
src/Util.irp.f
View File

@ -23,7 +23,18 @@ double precision function binom(k,n)
!DEC$ ATTRIBUTES ALIGN : 32 :: memo
integer, save :: ifirst
double precision :: fact
if (ifirst == 0) then
if (ifirst /= 0) then
if (n<20) then
binom = memo(k,n)
return
else
binom=fact(n)/(fact(k)*fact(n-k))
return
endif
else
ifirst = 1
integer :: i,j
do j=0,20
@ -31,12 +42,13 @@ double precision function binom(k,n)
memo(i,j) = fact(j)/(fact(i)*fact(j-i))
enddo
enddo
endif
if (n<20) then
binom = memo(k,n)
else
binom=fact(n)/(fact(k)*fact(n-k))
if (n<20) then
binom = memo(k,n)
return
else
binom=fact(n)/(fact(k)*fact(n-k))
endif
endif
end

62
src/eplf_function.irp.f
View File

@ -243,6 +243,16 @@ double precision function ao_eplf_integral_primitive_oneD(a,xa,i,b,xb,j,gmma,xr)
real, intent(in) :: xa,xb,xr ! Centers
integer, intent(in) :: i,j ! Powers of xa and xb
integer :: ii, jj, kk
integer, save :: ifirst=0
double precision, save :: f2(32)
if (ifirst == 0) then
ifirst = 1
do ii=1,32
f2(ii) = 2.d0**(-ii/2) * fact2(ii-1)
enddo
endif
ASSERT (a>0.)
ASSERT (b>0.)
@ -271,16 +281,43 @@ double precision function ao_eplf_integral_primitive_oneD(a,xa,i,b,xb,j,gmma,xr)
v=db_inv-xb
double precision :: fact2, binom, f, ac
double precision :: ui(-1:i), vj(-1:j), wij(-1:i+j+3)
!DEC$ ATTRIBUTES ALIGN : 32 :: ui, vj, wij
ui(i) = 1.d0
ui(i-1) = u
do ii=i-2,0,-1
ui(ii) = ui(ii+1) * u
enddo
vj(j) = 1.d0
vj(j-1) = v
do jj=j-2,0,-1
vj(jj) = vj(jj+1) * v
enddo
wij(0) = 1.d0
wij(1) = w
wij(2) = w*w
do ii=3,i+j
wij(ii) = wij(ii-1) * w
enddo
do ii=2,i+j
wij(ii) *= f2(ii)
enddo
ac = 0.d0
integer :: istart(0:1)
istart(0) = 0
istart(1) = 1
do ii=0,i
do jj=istart(mod(ii,2)),j,2
kk= -(ii+jj)/2
f=binom(ii,i)*binom(jj,j)*fact2(ii+jj-1)
ac += u**(i-ii) * v**(j-jj) * w**(ii+jj) * f * 2.d0**kk
enddo
if (istart(mod(ii,2))<=j) then
f=ui(ii)*binom(ii,i)
do jj=istart(mod(ii,2)),j,2
ac += vj(jj) * wij(ii+jj) * f *binom(jj,j)
enddo
endif
enddo
ao_eplf_integral_primitive_oneD = dsqrt_pi*w*c*ac
@ -305,12 +342,12 @@ double precision function ao_eplf_integral(i,j,gmma,center)
ASSERT(j<=ao_num)
ao_eplf_integral = 0.d0
do p=1,ao_prim_num(i)
buffer(p) = 0.d0
enddo
do q=1,ao_prim_num(j)
do p=1,ao_prim_num(i)
if (abs(ao_coef(q,j)*ao_coef(p,i)) < 1.e-9) then
cycle
endif
integral = &
ao_eplf_integral_primitive_oneD( &
ao_expo(p,i), &
@ -321,7 +358,7 @@ double precision function ao_eplf_integral(i,j,gmma,center)
ao_power(j,1), &
gmma, &
center(1))
if (integral /= 0.d0) then
if (dabs(integral) > 1.d-15) then
integral *= &
ao_eplf_integral_primitive_oneD( &
ao_expo(p,i), &
@ -332,7 +369,7 @@ double precision function ao_eplf_integral(i,j,gmma,center)
ao_power(j,2), &
gmma, &
center(2))
if (integral /= 0.d0) then
if (dabs(integral) > 1.d-15) then
integral *= &
ao_eplf_integral_primitive_oneD( &
ao_expo(p,i), &
@ -343,14 +380,11 @@ double precision function ao_eplf_integral(i,j,gmma,center)
ao_power(j,3), &
gmma, &
center(3))
buffer(p) += integral*ao_coef(p,i)*ao_coef(q,j)
ao_eplf_integral += integral*ao_coef(p,i)*ao_coef(q,j)
endif
endif
enddo
enddo
do p=1,ao_prim_num(i)
ao_eplf_integral += buffer(p)
enddo
end function