10
0
mirror of https://gitlab.com/scemama/eplf synced 2024-10-31 19:23:55 +01:00

Improved eplf function

This commit is contained in:
Anthony Scemama 2010-10-06 13:27:18 +02:00
parent 0e1c2f74a9
commit e9733253c3

View File

@ -21,7 +21,7 @@ BEGIN_PROVIDER [ double precision, ao_eplf_integral_matrix, (ao_num,ao_num) ]
double precision :: ao_eplf_integral double precision :: ao_eplf_integral
do i=1,ao_num do i=1,ao_num
do j=i,ao_num do j=i,ao_num
ao_eplf_integral_matrix(j,i) = ao_eplf_integral(j,i,eplf_gamma,point) ao_eplf_integral_matrix(j,i) = ao_eplf_integral(j,i,dble(eplf_gamma),point)
ao_eplf_integral_matrix(i,j) = ao_eplf_integral_matrix(j,i) ao_eplf_integral_matrix(i,j) = ao_eplf_integral_matrix(j,i)
enddo enddo
enddo enddo
@ -237,8 +237,8 @@ BEGIN_PROVIDER [ real, eplf_value_p ]
if ( (aa > 0.d0).and.(ab > 0.d0) ) then if ( (aa > 0.d0).and.(ab > 0.d0) ) then
aa = min(1.d0,aa) aa = min(1.d0,aa)
ab = min(1.d0,ab) ab = min(1.d0,ab)
aa = -dlog(aa)/eplf_gamma aa = -dlog(aa)
ab = -dlog(ab)/eplf_gamma ab = -dlog(ab)
aa = dsqrt(aa) aa = dsqrt(aa)
ab = dsqrt(ab) ab = dsqrt(ab)
eplf_value_p = (aa-ab)/(aa+ab+eps) eplf_value_p = (aa-ab)/(aa+ab+eps)
@ -249,93 +249,94 @@ BEGIN_PROVIDER [ real, eplf_value_p ]
END_PROVIDER END_PROVIDER
double precision function ao_eplf_integral_primitive_oneD_numeric(a,xa,i,b,xb,j,gmma,xr) !double precision function ao_eplf_integral_primitive_oneD_numeric(a,xa,i,b,xb,j,gmma,xr)
implicit none ! implicit none
include 'constants.F' ! include 'constants.F'
!
real, intent(in) :: a,b,gmma ! Exponents ! real, intent(in) :: a,b,gmma ! Exponents
real, intent(in) :: xa,xb,xr ! Centers ! real, intent(in) :: xa,xb,xr ! Centers
integer, intent(in) :: i,j ! Powers of xa and xb ! integer, intent(in) :: i,j ! Powers of xa and xb
integer,parameter :: Npoints=10000 ! integer,parameter :: Npoints=10000
real :: x, xmin, xmax, dx ! real :: x, xmin, xmax, dx
!
ASSERT (a>0.) ! ASSERT (a>0.)
ASSERT (b>0.) ! ASSERT (b>0.)
ASSERT (i>=0) ! ASSERT (i>=0)
ASSERT (j>=0) ! ASSERT (j>=0)
!
xmin = min(xa,xb) ! xmin = min(xa,xb)
xmax = max(xa,xb) ! xmax = max(xa,xb)
xmin = min(xmin,xr) - 10. ! xmin = min(xmin,xr) - 10.
xmax = max(xmax,xr) + 10. ! xmax = max(xmax,xr) + 10.
dx = (xmax-xmin)/real(Npoints) ! dx = (xmax-xmin)/real(Npoints)
!
real :: dtemp ! real :: dtemp
dtemp = 0. ! dtemp = 0.
x = xmin ! x = xmin
integer :: k ! integer :: k
do k=1,Npoints ! do k=1,Npoints
dtemp += & ! dtemp += &
(x-xa)**i * (x-xb)**j * exp(-(a*(x-xa)**2+b*(x-xb)**2+gmma*(x-xr)**2)) ! (x-xa)**i * (x-xb)**j * exp(-(a*(x-xa)**2+b*(x-xb)**2+gmma*(x-xr)**2))
x = x+dx ! x = x+dx
enddo ! enddo
ao_eplf_integral_primitive_oneD_numeric = dtemp*dx ! ao_eplf_integral_primitive_oneD_numeric = dtemp*dx
!
end function !end function
!
double precision function ao_eplf_integral_numeric(i,j,gmma,center) !double precision function ao_eplf_integral_numeric(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
double precision :: ao_eplf_integral_primitive_oneD_numeric ! double precision :: ao_eplf_integral_primitive_oneD_numeric
real :: gmma, center(3), c ! real :: gmma, center(3), c
!
!
ao_eplf_integral_numeric = 0.d0 ! ao_eplf_integral_numeric = 0.d0
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)
c = ao_coef(i,p)*ao_coef(j,q) ! c = ao_coef(i,p)*ao_coef(j,q)
integral = & ! integral = &
ao_eplf_integral_primitive_oneD_numeric( & ! ao_eplf_integral_primitive_oneD_numeric( &
ao_expo(i,p), & ! ao_expo(i,p), &
nucl_coord(ao_nucl(i),1), & ! nucl_coord(ao_nucl(i),1), &
ao_power(i,1), & ! ao_power(i,1), &
ao_expo(j,q), & ! ao_expo(j,q), &
nucl_coord(ao_nucl(j),1), & ! nucl_coord(ao_nucl(j),1), &
ao_power(j,1), & ! ao_power(j,1), &
gmma, & ! gmma, &
center(1)) * & ! center(1)) * &
ao_eplf_integral_primitive_oneD_numeric( & ! ao_eplf_integral_primitive_oneD_numeric( &
ao_expo(i,p), & ! ao_expo(i,p), &
nucl_coord(ao_nucl(i),2), & ! nucl_coord(ao_nucl(i),2), &
ao_power(i,2), & ! ao_power(i,2), &
ao_expo(j,q), & ! ao_expo(j,q), &
nucl_coord(ao_nucl(j),2), & ! nucl_coord(ao_nucl(j),2), &
ao_power(j,2), & ! ao_power(j,2), &
gmma, & ! gmma, &
center(2)) * & ! center(2)) * &
ao_eplf_integral_primitive_oneD_numeric( & ! ao_eplf_integral_primitive_oneD_numeric( &
ao_expo(i,p), & ! ao_expo(i,p), &
nucl_coord(ao_nucl(i),3), & ! nucl_coord(ao_nucl(i),3), &
ao_power(i,3), & ! ao_power(i,3), &
ao_expo(j,q), & ! ao_expo(j,q), &
nucl_coord(ao_nucl(j),3), & ! nucl_coord(ao_nucl(j),3), &
ao_power(j,3), & ! ao_power(j,3), &
gmma, & ! gmma, &
center(3)) ! center(3))
ao_eplf_integral_numeric = ao_eplf_integral_numeric + c*integral ! ao_eplf_integral_numeric = ao_eplf_integral_numeric + c*integral
enddo ! enddo
enddo ! enddo
!
end function !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'
real, intent(in) :: a,b,gmma ! Exponents real, intent(in) :: a,b ! Exponents
double precision , intent(in) :: gmma ! eplf_gamma
real, intent(in) :: xa,xb,xr ! Centers real, intent(in) :: xa,xb,xr ! Centers
integer, intent(in) :: i,j ! Powers of xa and xb integer, intent(in) :: i,j ! Powers of xa and xb
integer :: ii, jj, kk, ll integer :: ii, jj, kk, ll
@ -350,7 +351,6 @@ double precision function ao_eplf_integral_primitive_oneD(a,xa,i,b,xb,j,gmma,xr)
ASSERT (j>=0) ASSERT (j>=0)
! Gaussian product ! Gaussian product
! Inlined Gaussian products (same as call gaussian_product)
real :: t(2), xab(2), ab(2) real :: t(2), xab(2), ab(2)
inv_p(1) = 1.d0/(a+b) inv_p(1) = 1.d0/(a+b)
p1 = a+b p1 = a+b
@ -364,17 +364,23 @@ double precision function ao_eplf_integral_primitive_oneD(a,xa,i,b,xb,j,gmma,xr)
t(2) = (p1*xp1+gmma*xr) t(2) = (p1*xp1+gmma*xr)
xab(2) = xp1-xr xab(2) = xp1-xr
xp = t(2)*inv_p(2) xp = t(2)*inv_p(2)
c = real(ab(1)*inv_p(1)*xab(1)**2 + & c = dble(ab(1)*inv_p(1)*xab(1)**2 + &
ab(2)*inv_p(2)*xab(2)**2) ab(2)*inv_p(2)*xab(2)**2)
if ( c > 32.d0 ) then
! double precision, save :: c_accu(2)
! c_accu(1) += c
! c_accu(2) += 1.d0
! print *, c_accu(1)/c_accu(2)
if ( c > 32.d0 ) then ! Cut-off on exp(-32)
ao_eplf_integral_primitive_oneD = 0.d0 ao_eplf_integral_primitive_oneD = 0.d0
return return
endif endif
c = exp(-c) c = exp(-c)
!S(0,0) = dsqrt(pi*inv_p(2))*c
S(0,0) = 1.d0 ! Factor is applied at the end
! Obara-Saika recursion ! Obara-Saika recursion
S(0,0) = 1.d0
do ii=1,max(i,j) do ii=1,max(i,j)
di(ii) = 0.5d0*inv_p(2)*dble(ii) di(ii) = 0.5d0*inv_p(2)*dble(ii)
@ -382,14 +388,14 @@ double precision function ao_eplf_integral_primitive_oneD(a,xa,i,b,xb,j,gmma,xr)
xab(1) = xp-xa xab(1) = xp-xa
xab(2) = xp-xb xab(2) = xp-xb
S(1,0) = xab(1) * S(0,0) S(1,0) = xab(1) ! * S(0,0)
if (i>1) then if (i>1) then
do ii=1,i-1 do ii=1,i-1
S(ii+1,0) = xab(1) * S(ii,0) + di(ii)*S(ii-1,0) S(ii+1,0) = xab(1) * S(ii,0) + di(ii)*S(ii-1,0)
enddo enddo
endif endif
S(0,1) = xab(2) * S(0,0) S(0,1) = xab(2) ! * S(0,0)
if (j>1) then if (j>1) then
do jj=1,j-1 do jj=1,j-1
S(0,jj+1) = xab(2) * S(0,jj) + di(jj)*S(0,jj-1) S(0,jj+1) = xab(2) * S(0,jj) + di(jj)*S(0,jj-1)
@ -403,7 +409,7 @@ double precision function ao_eplf_integral_primitive_oneD(a,xa,i,b,xb,j,gmma,xr)
enddo enddo
enddo enddo
ao_eplf_integral_primitive_oneD = dsqrt(pi*inv_p(2))*S(i,j)*c ! Application of the factor of S(0,0) ao_eplf_integral_primitive_oneD = dsqrt(pi*inv_p(2))*c*S(i,j)
end function end function
@ -412,6 +418,7 @@ double precision function ao_eplf_integral_primitive_oneD(a,xa,i,b,xb,j,gmma,xr)
! include 'constants.F' ! include 'constants.F'
!! !!
! real, intent(in) :: a,b,gmma ! Exponents ! real, intent(in) :: a,b,gmma ! Exponents
! double precision, intent(in) :: gmma
! real, intent(in) :: xa,xb,xr ! Centers ! real, intent(in) :: xa,xb,xr ! Centers
! integer, intent(in) :: i,j ! Powers of xa and xb ! integer, intent(in) :: i,j ! Powers of xa and xb
! integer :: ii, jj, kk, ll ! integer :: ii, jj, kk, ll
@ -485,16 +492,18 @@ double precision function ao_eplf_integral_primitive_oneD(a,xa,i,b,xb,j,gmma,xr)
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
real, intent(in) :: center(3)
double precision, intent(in) :: gmma
!DEC$ ATTRIBUTES FORCEINLINE
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)
double precision :: buffer(100) double precision :: buffer(100)
ASSERT(i>0) ASSERT(i>0)
ASSERT(j>0) ASSERT(j>0)
ASSERT(ao_prim_num_max < 100)
ASSERT(i<=ao_num) ASSERT(i<=ao_num)
ASSERT(j<=ao_num) ASSERT(j<=ao_num)