mirror of https://gitlab.com/scemama/eplf
Added recursive fortran func
This commit is contained in:
parent
5069a720a0
commit
5de7b72028
5
Makefile
5
Makefile
|
|
@ -1,6 +1,7 @@
|
|||
IRPF90 = irpf90 -DMPI#-a -d
|
||||
FC = mpif90
|
||||
FCFLAGS= -O3 -xT
|
||||
IRPF90 = irpf90 -DMPI#-a -d
|
||||
FC = mpif90 -xT
|
||||
FCFLAGS= -O3 -g
|
||||
|
||||
#FC = gfortran -g -ffree-line-length-none
|
||||
#FCFLAGS=
|
||||
|
|
|
|||
101
eplf.irp.f
101
eplf.irp.f
|
|
@ -33,13 +33,15 @@ BEGIN_PROVIDER [ double precision, mo_eplf_integral_matrix, (mo_num,mo_num) ]
|
|||
do j=i,mo_num
|
||||
mo_eplf_integral_matrix(j,i) = 0.
|
||||
enddo
|
||||
enddo
|
||||
|
||||
do k=1,ao_num
|
||||
do i=1,mo_num
|
||||
do k=1,ao_num
|
||||
if (mo_coef(k,i) /= 0.) then
|
||||
do j=i,mo_num
|
||||
do l=1,ao_num
|
||||
mo_eplf_integral_matrix(j,i) = mo_eplf_integral_matrix(j,i) + &
|
||||
mo_coef(k,i)*mo_coef(l,j)*ao_eplf_integral_matrix(k,l)
|
||||
mo_coef(k,i)*mo_coef(l,j)*ao_eplf_integral_matrix(l,k)
|
||||
enddo
|
||||
enddo
|
||||
endif
|
||||
|
|
@ -203,7 +205,7 @@ double precision function ao_eplf_integral_numeric(i,j,gmma,center)
|
|||
|
||||
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_oneD2(a,xa,i,b,xb,j,gmma,xr)
|
||||
implicit none
|
||||
include 'constants.F'
|
||||
|
||||
|
|
@ -213,7 +215,7 @@ 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 :: S(0:i+1,0:j)
|
||||
double precision :: S(0:i+1,0:j+1)
|
||||
double precision :: inv_p, di(max(i,j)), dj(j), c, c1
|
||||
|
||||
ASSERT (a>0.)
|
||||
|
|
@ -229,23 +231,18 @@ double precision function ao_eplf_integral_primitive_oneD(a,xa,i,b,xb,j,gmma,xr)
|
|||
|
||||
! Obara-Saika recursion
|
||||
|
||||
if (i>0) then
|
||||
S(1,0) = (xp-xa) * S(0,0)
|
||||
endif
|
||||
if (j>0) then
|
||||
S(0,1) = (xp-xb) * S(0,0)
|
||||
endif
|
||||
|
||||
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)
|
||||
|
|
@ -259,7 +256,69 @@ double precision function ao_eplf_integral_primitive_oneD(a,xa,i,b,xb,j,gmma,xr)
|
|||
enddo
|
||||
enddo
|
||||
|
||||
ao_eplf_integral_primitive_oneD = S(i,j)
|
||||
ao_eplf_integral_primitive_oneD2 = S(i,j)
|
||||
|
||||
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 :: S00, xpa, xpb
|
||||
double precision :: inv_p,c,c1
|
||||
double precision :: ObaraS
|
||||
|
||||
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
|
||||
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)
|
||||
|
||||
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 = 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 = 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 = res + 0.5d0*dble(i-1)*inv_p*ObaraS(i-2,j,xpa,xpb,inv_p,S00)
|
||||
endif
|
||||
res = 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
|
||||
|
||||
|
|
@ -313,3 +372,21 @@ double precision function ao_eplf_integral(i,j,gmma,center)
|
|||
|
||||
end function
|
||||
|
||||
double precision function mo_eplf_integral(i,j)
|
||||
implicit none
|
||||
integer :: i, j, k, l
|
||||
PROVIDE ao_eplf_integral_matrix
|
||||
PROVIDE mo_coef
|
||||
mo_eplf_integral = 0.d0
|
||||
|
||||
do k=1,ao_num
|
||||
if (mo_coef(k,i) /= 0.) then
|
||||
do l=1,ao_num
|
||||
mo_eplf_integral = mo_eplf_integral + &
|
||||
mo_coef(k,i)*mo_coef(l,j)*ao_eplf_integral_matrix(k,l)
|
||||
enddo
|
||||
endif
|
||||
enddo
|
||||
|
||||
end function
|
||||
|
||||
|
|
|
|||
Loading…
Reference in New Issue