mirror of
https://gitlab.com/scemama/eplf
synced 2024-12-22 12:23:50 +01:00
Improved EPLF function for large gammas
This commit is contained in:
Anthony Scemama
parent
b638782709
commit
42d8ba45ec
@ -15,6 +15,12 @@
|
||||
!
|
||||
!end function
|
||||
|
||||
|
||||
double precision function binom(k,n)
|
||||
implicit double precision(a-h,o-z)
|
||||
binom=fact(n)/(fact(k)*fact(n-k))
|
||||
end
|
||||
|
||||
double precision function fact2(n)
|
||||
implicit none
|
||||
integer :: n
|
||||
@ -80,7 +86,7 @@ double precision function rintgauss(n)
|
||||
|
||||
integer :: n
|
||||
|
||||
rintgauss = sqrt(pi)
|
||||
rintgauss = dsqrt_pi
|
||||
if ( n == 0 ) then
|
||||
return
|
||||
else if ( n == 1 ) then
|
||||
|
||||
@ -1 +1,2 @@
|
||||
double precision, parameter :: pi=3.14159265359d0
|
||||
double precision, parameter :: dsqrt_pi=1.77245385091d0
|
||||
|
||||
@ -5,9 +5,12 @@ BEGIN_PROVIDER [ real, eplf_gamma ]
|
||||
END_DOC
|
||||
include 'constants.F'
|
||||
real :: eps, N
|
||||
N = 0.1
|
||||
eps = -real(dlog(tiny(1.d0)))
|
||||
eplf_gamma = (4./3.*pi*density_value_p/N)**(2./3.) * eps
|
||||
! N = 0.1
|
||||
! eps = -real(dlog(tiny(1.d0)))
|
||||
! eplf_gamma = (4./3.*pi*density_value_p/N)**(2./3.) * eps
|
||||
N=0.001
|
||||
eps = 50.
|
||||
eplf_gamma = eps**2 *0.5d0*(4./9.*pi*density_value_p * N**(-1./3.))**(2./3.)
|
||||
!eplf_gamma = 1.e10
|
||||
!eplf_gamma = 1.e5
|
||||
END_PROVIDER
|
||||
@ -339,155 +342,67 @@ double precision function ao_eplf_integral_primitive_oneD(a,xa,i,b,xb,j,gmma,xr)
|
||||
double precision , intent(in) :: gmma ! eplf_gamma
|
||||
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(2), di(max(i,j)), dj(j), c
|
||||
integer :: ii, jj, kk
|
||||
|
||||
ASSERT (a>0.)
|
||||
ASSERT (b>0.)
|
||||
ASSERT (i>=0)
|
||||
ASSERT (j>=0)
|
||||
|
||||
! 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 = dble(ab(1)*inv_p(1)*xab(1)**2 + &
|
||||
ab(2)*inv_p(2)*xab(2)**2)
|
||||
|
||||
! double precision, save :: c_accu(2)
|
||||
! c_accu(1) += c
|
||||
! c_accu(2) += 1.d0
|
||||
! print *, c_accu(1)/c_accu(2)
|
||||
double precision :: alpha, beta, delta, c
|
||||
alpha = a*xa*xa + b*xb*xb + gmma*xr*xr
|
||||
delta = a*xa + b*xb + gmma*xr
|
||||
beta = a + b + gmma
|
||||
c = alpha-delta**2/beta
|
||||
|
||||
if ( c > 32.d0 ) then ! Cut-off on exp(-32)
|
||||
ao_eplf_integral_primitive_oneD = 0.d0
|
||||
return
|
||||
endif
|
||||
|
||||
double precision :: u,v,w
|
||||
c = exp(-c)
|
||||
w = 1.d0/sqrt(beta)
|
||||
u=delta/beta-xa
|
||||
v=delta/beta-xb
|
||||
|
||||
! Obara-Saika recursion
|
||||
S(0,0) = 1.d0
|
||||
|
||||
do ii=1,max(i,j)
|
||||
di(ii) = 0.5d0*inv_p(2)*dble(ii)
|
||||
enddo
|
||||
|
||||
xab(1) = xp-xa
|
||||
xab(2) = xp-xb
|
||||
S(1,0) = xab(1) ! * S(0,0)
|
||||
if (i>1) then
|
||||
do ii=1,i-1
|
||||
S(ii+1,0) = xab(1) * S(ii,0) + di(ii)*S(ii-1,0)
|
||||
enddo
|
||||
endif
|
||||
|
||||
S(0,1) = xab(2) ! * S(0,0)
|
||||
if (j>1) then
|
||||
do jj=1,j-1
|
||||
S(0,jj+1) = xab(2) * S(0,jj) + di(jj)*S(0,jj-1)
|
||||
enddo
|
||||
endif
|
||||
|
||||
do jj=1,j
|
||||
S(1,jj) = xab(1) * S(0,jj) + di(jj) * S(0,jj-1)
|
||||
do ii=2,i
|
||||
S(ii,jj) = xab(1) * S(ii-1,jj) + di(ii-1) * S(ii-2,jj) + di(jj) * S(ii-1,jj-1)
|
||||
double precision :: fact2, binom, f, ac
|
||||
ac = 0.d0
|
||||
integer :: istart(0:1)
|
||||
istart = (/ 0, 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/dble(2**kk)
|
||||
enddo
|
||||
enddo
|
||||
ao_eplf_integral_primitive_oneD = dsqrt_pi*w*c*ac
|
||||
|
||||
ao_eplf_integral_primitive_oneD = dsqrt(pi*inv_p(2))*c*S(i,j)
|
||||
end function
|
||||
|
||||
end function
|
||||
double precision function rint(i,xa,a,j,xb,b,xr,gmma)
|
||||
implicit double precision(a-h,o-z)
|
||||
include 'constants.F'
|
||||
beta=a+b+gmma
|
||||
w=1.d0/dsqrt(beta)
|
||||
alpha=a*xa**2+b*xb**2+gmma*xr**2
|
||||
delta=xa*a+xb*b+xr*gmma
|
||||
u=delta/beta-xa
|
||||
v=delta/beta-xb
|
||||
ac=0.d0
|
||||
do n=0,i
|
||||
do m=0,j
|
||||
if(mod(n+m,2).eq.0)then
|
||||
k=(n+m)/2
|
||||
f=binom(n,i)*binom(m,j)*fact2(n+m-1)
|
||||
ac=ac+u**(i-n)*v**(j-m)*w**(n+m)*f/2.d0**k
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
rint=ac*dsqrt(pi)*w*dexp(-alpha+delta**2/beta)
|
||||
end
|
||||
|
||||
!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
|
||||
! double precision, intent(in) :: gmma
|
||||
! 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)
|
||||
implicit none
|
||||
|
||||
Loading…
Reference in New Issue
Block a user