mirror of
https://gitlab.com/scemama/eplf
synced 2024-12-23 04:43:52 +01:00
EPLF acceleration x2
This commit is contained in:
Anthony Scemama
parent
ec34d430dc
commit
07b0dacf6c
@ -217,61 +217,6 @@ double precision function ao_eplf_integral_numeric(i,j,gmma,center)
|
|||||||
|
|
||||||
end function
|
end function
|
||||||
|
|
||||||
!double precision function ao_eplf_integral_primitive_oneD(a,xa,i,b,xb,j,gmma,xr)
|
|
||||||
! implicit none
|
|
||||||
! include 'constants.F'
|
|
||||||
!
|
|
||||||
! real, intent(in) :: a,b,gmma ! Exponents
|
|
||||||
! real, intent(in) :: xa,xb,xr ! Centers
|
|
||||||
! integer, intent(in) :: i,j ! Powers of xa and xb
|
|
||||||
! integer :: ii, jj, kk, ll
|
|
||||||
! real :: xp1,xp
|
|
||||||
! real :: p1,p
|
|
||||||
! double precision :: S(0:i+1,0:j+1)
|
|
||||||
! double precision :: inv_p, di(max(i,j)), dj(j), c, c1
|
|
||||||
!
|
|
||||||
! ASSERT (a>0.)
|
|
||||||
! ASSERT (b>0.)
|
|
||||||
! 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)
|
|
||||||
! inv_p = 1.d0/p
|
|
||||||
! S(0,0) = dsqrt(pi*inv_p)*c*c1
|
|
||||||
!
|
|
||||||
! ! Obara-Saika recursion
|
|
||||||
!
|
|
||||||
! do ii=1,max(i,j)
|
|
||||||
! di(ii) = 0.5d0*inv_p*dble(ii)
|
|
||||||
! enddo
|
|
||||||
!
|
|
||||||
! S(1,0) = (xp-xa) * S(0,0)
|
|
||||||
! if (i>1) then
|
|
||||||
! do ii=1,i-1
|
|
||||||
! S(ii+1,0) = (xp-xa) * S(ii,0) + di(ii)*S(ii-1,0)
|
|
||||||
! enddo
|
|
||||||
! endif
|
|
||||||
!
|
|
||||||
! S(0,1) = (xp-xb) * S(0,0)
|
|
||||||
! if (j>1) then
|
|
||||||
! do jj=1,j-1
|
|
||||||
! S(0,jj+1) = (xp-xb) * S(0,jj) + di(jj)*S(0,jj-1)
|
|
||||||
! enddo
|
|
||||||
! endif
|
|
||||||
!
|
|
||||||
! do jj=1,j
|
|
||||||
! S(1,jj) = (xp-xa) * S(0,jj) + di(jj) * S(0,jj-1)
|
|
||||||
! do ii=2,i
|
|
||||||
! S(ii,jj) = (xp-xa) * S(ii-1,jj) + di(ii-1) * S(ii-2,jj) + di(jj) * S(ii-1,jj-1)
|
|
||||||
! enddo
|
|
||||||
! enddo
|
|
||||||
!
|
|
||||||
! ao_eplf_integral_primitive_oneD = S(i,j)
|
|
||||||
!
|
|
||||||
!end function
|
|
||||||
|
|
||||||
double precision function ao_eplf_integral_primitive_oneD(a,xa,i,b,xb,j,gmma,xr)
|
double precision function ao_eplf_integral_primitive_oneD(a,xa,i,b,xb,j,gmma,xr)
|
||||||
implicit none
|
implicit none
|
||||||
include 'constants.F'
|
include 'constants.F'
|
||||||
@ -282,19 +227,17 @@ double precision function ao_eplf_integral_primitive_oneD(a,xa,i,b,xb,j,gmma,xr)
|
|||||||
integer :: ii, jj, kk, ll
|
integer :: ii, jj, kk, ll
|
||||||
real :: xp1,xp
|
real :: xp1,xp
|
||||||
real :: p1,p
|
real :: p1,p
|
||||||
double precision :: xpa, xpb
|
double precision :: S(0:i+1,0:j+1)
|
||||||
double precision :: inv_p(2),S00, c
|
double precision :: inv_p(2), di(max(i,j)), dj(j), c
|
||||||
double precision :: ObaraS
|
|
||||||
|
|
||||||
ASSERT (a>0.)
|
ASSERT (a>0.)
|
||||||
ASSERT (b>0.)
|
ASSERT (b>0.)
|
||||||
ASSERT (i>=0)
|
ASSERT (i>=0)
|
||||||
ASSERT (j>=0)
|
ASSERT (j>=0)
|
||||||
|
|
||||||
! Gaussian products
|
! Gaussian product
|
||||||
!double precision :: t(2), xab(2), ab(2)
|
! Inlined Gaussian products (same as call gaussian_product)
|
||||||
real :: 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)
|
inv_p(1) = 1.d0/(a+b)
|
||||||
p1 = a+b
|
p1 = a+b
|
||||||
ab(1) = a*b
|
ab(1) = a*b
|
||||||
@ -310,60 +253,136 @@ double precision function ao_eplf_integral_primitive_oneD(a,xa,i,b,xb,j,gmma,xr)
|
|||||||
c = exp(- real(ab(1)*inv_p(1)*xab(1)**2 + &
|
c = exp(- real(ab(1)*inv_p(1)*xab(1)**2 + &
|
||||||
ab(2)*inv_p(2)*xab(2)**2) )
|
ab(2)*inv_p(2)*xab(2)**2) )
|
||||||
|
|
||||||
|
! inv_p = 1.d0/p
|
||||||
|
S(0,0) = dsqrt(pi*inv_p(2))*c
|
||||||
|
|
||||||
xpa = xp-xa
|
! Obara-Saika recursion
|
||||||
xpb = xp-xb
|
|
||||||
S00 = dsqrt(pi*inv_p(2))*c
|
|
||||||
ao_eplf_integral_primitive_oneD = ObaraS(i,j,xpa,xpb,inv_p(2),S00)
|
|
||||||
|
|
||||||
end function
|
do ii=1,max(i,j)
|
||||||
|
di(ii) = 0.5d0*inv_p(2)*dble(ii)
|
||||||
|
enddo
|
||||||
|
|
||||||
recursive double precision function ObaraS(i,j,xpa,xpb,inv_p,S00) result(res)
|
S(1,0) = (xp-xa) * S(0,0)
|
||||||
implicit none
|
if (i>1) then
|
||||||
integer, intent(in) :: i, j
|
do ii=1,i-1
|
||||||
double precision, intent(in) :: xpa, xpb, inv_p
|
S(ii+1,0) = (xp-xa) * S(ii,0) + di(ii)*S(ii-1,0)
|
||||||
double precision,intent(in) :: S00
|
enddo
|
||||||
|
endif
|
||||||
|
|
||||||
if (i == 0) then
|
S(0,1) = (xp-xb) * S(0,0)
|
||||||
if (j == 0) then
|
|
||||||
res = S00
|
|
||||||
else ! (j>0)
|
|
||||||
res = xpb*ObaraS(0,j-1,xpa,xpb,inv_p,S00)
|
|
||||||
if (j>1) then
|
if (j>1) then
|
||||||
res += 0.5d0*dble(j-1)*inv_p*ObaraS(0,j-2,xpa,xpb,inv_p,S00)
|
do jj=1,j-1
|
||||||
|
S(0,jj+1) = (xp-xb) * S(0,jj) + di(jj)*S(0,jj-1)
|
||||||
|
enddo
|
||||||
endif
|
endif
|
||||||
endif ! (i==0).and.(j>0)
|
|
||||||
else ! (i>0)
|
do jj=1,j
|
||||||
if (j==0) then
|
S(1,jj) = (xp-xa) * S(0,jj) + di(jj) * S(0,jj-1)
|
||||||
res = xpa*ObaraS(i-1,0,xpa,xpb,inv_p,S00)
|
do ii=2,i
|
||||||
if (i>1) then
|
S(ii,jj) = (xp-xa) * S(ii-1,jj) + di(ii-1) * S(ii-2,jj) + di(jj) * S(ii-1,jj-1)
|
||||||
res += 0.5d0*dble(i-1)*inv_p*ObaraS(i-2,0,xpa,xpb,inv_p,S00)
|
enddo
|
||||||
endif
|
enddo
|
||||||
else ! (i>0).and.(j>0)
|
|
||||||
res = xpa * ObaraS(i-1,j,xpa,xpb,inv_p,S00)
|
ao_eplf_integral_primitive_oneD = S(i,j)
|
||||||
if (i>1) then
|
|
||||||
res += 0.5d0*dble(i-1)*inv_p*ObaraS(i-2,j,xpa,xpb,inv_p,S00)
|
|
||||||
endif
|
|
||||||
res += 0.5d0*dble(j)*inv_p*ObaraS(i-1,j-1,xpa,xpb,inv_p,S00)
|
|
||||||
endif ! (i>0).and.(j>0)
|
|
||||||
endif ! (i>0)
|
|
||||||
|
|
||||||
end function
|
end function
|
||||||
|
|
||||||
|
!double precision function ao_eplf_integral_primitive_oneD(a,xa,i,b,xb,j,gmma,xr)
|
||||||
|
! implicit none
|
||||||
|
! include 'constants.F'
|
||||||
|
!!
|
||||||
|
! real, intent(in) :: a,b,gmma ! Exponents
|
||||||
|
! real, intent(in) :: xa,xb,xr ! Centers
|
||||||
|
! integer, intent(in) :: i,j ! Powers of xa and xb
|
||||||
|
! integer :: ii, jj, kk, ll
|
||||||
|
! real :: xp1,xp
|
||||||
|
! real :: p1,p
|
||||||
|
! double precision :: xpa, xpb
|
||||||
|
! double precision :: inv_p(2),S00, c
|
||||||
|
! double precision :: ObaraS
|
||||||
|
!!
|
||||||
|
! ASSERT (a>0.)
|
||||||
|
! ASSERT (b>0.)
|
||||||
|
! ASSERT (i>=0)
|
||||||
|
! ASSERT (j>=0)
|
||||||
|
!!
|
||||||
|
! ! Inlined Gaussian products (same as call gaussian_product)
|
||||||
|
! real :: t(2), xab(2), ab(2)
|
||||||
|
! 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) )
|
||||||
|
!!
|
||||||
|
! xpa = xp-xa
|
||||||
|
! xpb = xp-xb
|
||||||
|
! S00 = sqrt(real(pi*inv_p(2)))*c
|
||||||
|
! ao_eplf_integral_primitive_oneD = ObaraS(i,j,xpa,xpb,inv_p(2),S00)
|
||||||
|
!!
|
||||||
|
!end function
|
||||||
|
!!
|
||||||
|
!recursive double precision function ObaraS(i,j,xpa,xpb,inv_p,S00) result(res)
|
||||||
|
! implicit none
|
||||||
|
! integer, intent(in) :: i, j
|
||||||
|
! double precision, intent(in) :: xpa, xpb, inv_p
|
||||||
|
! double precision,intent(in) :: S00
|
||||||
|
!!
|
||||||
|
! if (i == 0) then
|
||||||
|
! if (j == 0) then
|
||||||
|
! res = S00
|
||||||
|
! else ! (j>0)
|
||||||
|
! res = xpb*ObaraS(0,j-1,xpa,xpb,inv_p,S00)
|
||||||
|
! if (j>1) then
|
||||||
|
! res += 0.5d0*dble(j-1)*inv_p*ObaraS(0,j-2,xpa,xpb,inv_p,S00)
|
||||||
|
! endif
|
||||||
|
! endif ! (i==0).and.(j>0)
|
||||||
|
! else ! (i>0)
|
||||||
|
! if (j==0) then
|
||||||
|
! res = xpa*ObaraS(i-1,0,xpa,xpb,inv_p,S00)
|
||||||
|
! if (i>1) then
|
||||||
|
! res += 0.5d0*dble(i-1)*inv_p*ObaraS(i-2,0,xpa,xpb,inv_p,S00)
|
||||||
|
! endif
|
||||||
|
! else ! (i>0).and.(j>0)
|
||||||
|
! res = xpa * ObaraS(i-1,j,xpa,xpb,inv_p,S00)
|
||||||
|
! if (i>1) then
|
||||||
|
! res += 0.5d0*dble(i-1)*inv_p*ObaraS(i-2,j,xpa,xpb,inv_p,S00)
|
||||||
|
! endif
|
||||||
|
! res += 0.5d0*dble(j)*inv_p*ObaraS(i-1,j-1,xpa,xpb,inv_p,S00)
|
||||||
|
! endif ! (i>0).and.(j>0)
|
||||||
|
! endif ! (i>0)
|
||||||
|
!!
|
||||||
|
!end function
|
||||||
|
|
||||||
double precision function ao_eplf_integral(i,j,gmma,center)
|
double precision function ao_eplf_integral(i,j,gmma,center)
|
||||||
implicit none
|
implicit none
|
||||||
integer, intent(in) :: i, j
|
integer, intent(in) :: i, j
|
||||||
integer :: p,q,k
|
integer :: p,q,k
|
||||||
double precision :: integral
|
double precision :: integral
|
||||||
|
!DEC$ ATTRIBUTES FORCEINLINE
|
||||||
double precision :: ao_eplf_integral_primitive_oneD
|
double precision :: ao_eplf_integral_primitive_oneD
|
||||||
real :: gmma, center(3)
|
real :: gmma, center(3)
|
||||||
|
double precision :: buffer(100)
|
||||||
|
|
||||||
|
|
||||||
ASSERT(i>0)
|
ASSERT(i>0)
|
||||||
ASSERT(j>0)
|
ASSERT(j>0)
|
||||||
ASSERT(i<=ao_num)
|
ASSERT(i<=ao_num)
|
||||||
ASSERT(j<=ao_num)
|
ASSERT(j<=ao_num)
|
||||||
|
|
||||||
ao_eplf_integral = 0.
|
ao_eplf_integral = 0.d0
|
||||||
|
do p=1,ao_prim_num_max
|
||||||
|
buffer(p) = 0.d0
|
||||||
|
enddo
|
||||||
|
|
||||||
do q=1,ao_prim_num(j)
|
do q=1,ao_prim_num(j)
|
||||||
do p=1,ao_prim_num(i)
|
do p=1,ao_prim_num(i)
|
||||||
integral = &
|
integral = &
|
||||||
@ -394,9 +413,13 @@ double precision function ao_eplf_integral(i,j,gmma,center)
|
|||||||
ao_power(j,3), &
|
ao_power(j,3), &
|
||||||
gmma, &
|
gmma, &
|
||||||
center(3))
|
center(3))
|
||||||
ao_eplf_integral += integral*ao_coef(i,p)*ao_coef(j,q)
|
! ao_eplf_integral += integral*ao_coef(i,p)*ao_coef(j,q)
|
||||||
|
buffer(p) += integral*ao_coef(i,p)*ao_coef(j,q)
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
|
do p=1,ao_prim_num_max
|
||||||
|
ao_eplf_integral += buffer(p)
|
||||||
|
enddo
|
||||||
|
|
||||||
end function
|
end function
|
||||||
|
|
||||||
|
|||||||
Loading…
Reference in New Issue
Block a user