75% acceleration

This commit is contained in:
Anthony Scemama 2009-12-22 00:07:34 +01:00
parent 4a6ce8c3d8
commit ec34d430dc
1 changed files with 23 additions and 8 deletions

View File

@ -282,8 +282,8 @@ double precision function ao_eplf_integral_primitive_oneD(a,xa,i,b,xb,j,gmma,xr)
integer :: ii, jj, kk, ll
real :: xp1,xp
real :: p1,p
double precision :: S00, xpa, xpb
double precision :: inv_p,c,c1
double precision :: xpa, xpb
double precision :: inv_p(2),S00, c
double precision :: ObaraS
ASSERT (a>0.)
@ -291,15 +291,30 @@ double precision function ao_eplf_integral_primitive_oneD(a,xa,i,b,xb,j,gmma,xr)
ASSERT (i>=0)
ASSERT (j>=0)
! Gaussian product
call gaussian_product(a,xa,b,xb,c1,p1,xp1)
call gaussian_product(p1,xp1,gmma,xr,c,p,xp)
! Gaussian products
!double precision :: t(2), xab(2), ab(2)
real :: t(2), xab(2), ab(2)
!call gaussian_product(a,xa,b,xb,c1,p1,xp1)
inv_p(1) = 1.d0/(a+b)
p1 = a+b
ab(1) = a*b
inv_p(2) = 1.d0/(p1+gmma)
t(1) = (a*xa+b*xb)
xab(1) = xa-xb
xp1 = t(1)*inv_p(1)
p = p1+gmma
ab(2) = p1*gmma
t(2) = (p1*xp1+gmma*xr)
xab(2) = xp1-xr
xp = t(2)*inv_p(2)
c = exp(- real(ab(1)*inv_p(1)*xab(1)**2 + &
ab(2)*inv_p(2)*xab(2)**2) )
inv_p = 1.d0/p
S00 = dsqrt(pi*inv_p)*c*c1
xpa = xp-xa
xpb = xp-xb
ao_eplf_integral_primitive_oneD = ObaraS(i,j,xpa,xpb,inv_p,S00)
S00 = dsqrt(pi*inv_p(2))*c
ao_eplf_integral_primitive_oneD = ObaraS(i,j,xpa,xpb,inv_p(2),S00)
end function