From e9733253c3dc46dacbe5c278a6f3d77b4850d097 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 6 Oct 2010 13:27:18 +0200 Subject: [PATCH] Improved eplf function --- src/eplf_function.irp.f | 233 +++++++++++++++++++++------------------- 1 file changed, 121 insertions(+), 112 deletions(-) diff --git a/src/eplf_function.irp.f b/src/eplf_function.irp.f index 67976d1..eb63902 100644 --- a/src/eplf_function.irp.f +++ b/src/eplf_function.irp.f @@ -21,7 +21,7 @@ BEGIN_PROVIDER [ double precision, ao_eplf_integral_matrix, (ao_num,ao_num) ] double precision :: ao_eplf_integral do i=1,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) enddo enddo @@ -237,8 +237,8 @@ BEGIN_PROVIDER [ real, eplf_value_p ] if ( (aa > 0.d0).and.(ab > 0.d0) ) then aa = min(1.d0,aa) ab = min(1.d0,ab) - aa = -dlog(aa)/eplf_gamma - ab = -dlog(ab)/eplf_gamma + aa = -dlog(aa) + ab = -dlog(ab) aa = dsqrt(aa) ab = dsqrt(ab) eplf_value_p = (aa-ab)/(aa+ab+eps) @@ -249,93 +249,94 @@ BEGIN_PROVIDER [ real, eplf_value_p ] END_PROVIDER -double precision function ao_eplf_integral_primitive_oneD_numeric(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,parameter :: Npoints=10000 - real :: x, xmin, xmax, dx - - ASSERT (a>0.) - ASSERT (b>0.) - ASSERT (i>=0) - ASSERT (j>=0) - - xmin = min(xa,xb) - xmax = max(xa,xb) - xmin = min(xmin,xr) - 10. - xmax = max(xmax,xr) + 10. - dx = (xmax-xmin)/real(Npoints) - - real :: dtemp - dtemp = 0. - x = xmin - integer :: k - do k=1,Npoints - dtemp += & - (x-xa)**i * (x-xb)**j * exp(-(a*(x-xa)**2+b*(x-xb)**2+gmma*(x-xr)**2)) - x = x+dx - enddo - ao_eplf_integral_primitive_oneD_numeric = dtemp*dx - -end function - -double precision function ao_eplf_integral_numeric(i,j,gmma,center) - implicit none - integer, intent(in) :: i, j - integer :: p,q,k - double precision :: integral - double precision :: ao_eplf_integral_primitive_oneD_numeric - real :: gmma, center(3), c - - - ao_eplf_integral_numeric = 0.d0 - do q=1,ao_prim_num(j) - do p=1,ao_prim_num(i) - c = ao_coef(i,p)*ao_coef(j,q) - integral = & - ao_eplf_integral_primitive_oneD_numeric( & - ao_expo(i,p), & - nucl_coord(ao_nucl(i),1), & - ao_power(i,1), & - ao_expo(j,q), & - nucl_coord(ao_nucl(j),1), & - ao_power(j,1), & - gmma, & - center(1)) * & - ao_eplf_integral_primitive_oneD_numeric( & - ao_expo(i,p), & - nucl_coord(ao_nucl(i),2), & - ao_power(i,2), & - ao_expo(j,q), & - nucl_coord(ao_nucl(j),2), & - ao_power(j,2), & - gmma, & - center(2)) * & - ao_eplf_integral_primitive_oneD_numeric( & - ao_expo(i,p), & - nucl_coord(ao_nucl(i),3), & - ao_power(i,3), & - ao_expo(j,q), & - nucl_coord(ao_nucl(j),3), & - ao_power(j,3), & - gmma, & - center(3)) - ao_eplf_integral_numeric = ao_eplf_integral_numeric + c*integral - enddo - enddo - -end function +!double precision function ao_eplf_integral_primitive_oneD_numeric(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,parameter :: Npoints=10000 +! real :: x, xmin, xmax, dx +! +! ASSERT (a>0.) +! ASSERT (b>0.) +! ASSERT (i>=0) +! ASSERT (j>=0) +! +! xmin = min(xa,xb) +! xmax = max(xa,xb) +! xmin = min(xmin,xr) - 10. +! xmax = max(xmax,xr) + 10. +! dx = (xmax-xmin)/real(Npoints) +! +! real :: dtemp +! dtemp = 0. +! x = xmin +! integer :: k +! do k=1,Npoints +! dtemp += & +! (x-xa)**i * (x-xb)**j * exp(-(a*(x-xa)**2+b*(x-xb)**2+gmma*(x-xr)**2)) +! x = x+dx +! enddo +! ao_eplf_integral_primitive_oneD_numeric = dtemp*dx +! +!end function +! +!double precision function ao_eplf_integral_numeric(i,j,gmma,center) +! implicit none +! integer, intent(in) :: i, j +! integer :: p,q,k +! double precision :: integral +! double precision :: ao_eplf_integral_primitive_oneD_numeric +! real :: gmma, center(3), c +! +! +! ao_eplf_integral_numeric = 0.d0 +! do q=1,ao_prim_num(j) +! do p=1,ao_prim_num(i) +! c = ao_coef(i,p)*ao_coef(j,q) +! integral = & +! ao_eplf_integral_primitive_oneD_numeric( & +! ao_expo(i,p), & +! nucl_coord(ao_nucl(i),1), & +! ao_power(i,1), & +! ao_expo(j,q), & +! nucl_coord(ao_nucl(j),1), & +! ao_power(j,1), & +! gmma, & +! center(1)) * & +! ao_eplf_integral_primitive_oneD_numeric( & +! ao_expo(i,p), & +! nucl_coord(ao_nucl(i),2), & +! ao_power(i,2), & +! ao_expo(j,q), & +! nucl_coord(ao_nucl(j),2), & +! ao_power(j,2), & +! gmma, & +! center(2)) * & +! ao_eplf_integral_primitive_oneD_numeric( & +! ao_expo(i,p), & +! nucl_coord(ao_nucl(i),3), & +! ao_power(i,3), & +! ao_expo(j,q), & +! nucl_coord(ao_nucl(j),3), & +! ao_power(j,3), & +! gmma, & +! center(3)) +! ao_eplf_integral_numeric = ao_eplf_integral_numeric + c*integral +! enddo +! enddo +! +!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) :: a,b ! Exponents + 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 @@ -350,31 +351,36 @@ double precision function ao_eplf_integral_primitive_oneD(a,xa,i,b,xb,j,gmma,xr) ASSERT (j>=0) ! Gaussian product - ! 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 = real(ab(1)*inv_p(1)*xab(1)**2 + & - ab(2)*inv_p(2)*xab(2)**2) - if ( c > 32.d0 ) then - ao_eplf_integral_primitive_oneD = 0.d0 - return - endif - c = exp(-c) - !S(0,0) = dsqrt(pi*inv_p(2))*c - S(0,0) = 1.d0 ! Factor is applied at the end + 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) + + if ( c > 32.d0 ) then ! Cut-off on exp(-32) + ao_eplf_integral_primitive_oneD = 0.d0 + return + endif + + c = exp(-c) ! Obara-Saika recursion + S(0,0) = 1.d0 do ii=1,max(i,j) 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(2) = xp-xb - S(1,0) = xab(1) * S(0,0) + 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) + 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) @@ -403,7 +409,7 @@ double precision function ao_eplf_integral_primitive_oneD(a,xa,i,b,xb,j,gmma,xr) 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 @@ -412,6 +418,7 @@ double precision function ao_eplf_integral_primitive_oneD(a,xa,i,b,xb,j,gmma,xr) ! 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 @@ -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) implicit none integer, intent(in) :: i, j + real, intent(in) :: center(3) + double precision, intent(in) :: gmma +!DEC$ ATTRIBUTES FORCEINLINE integer :: p,q,k double precision :: integral -!DEC$ ATTRIBUTES FORCEINLINE double precision :: ao_eplf_integral_primitive_oneD - real :: gmma, center(3) double precision :: buffer(100) ASSERT(i>0) ASSERT(j>0) + ASSERT(ao_prim_num_max < 100) ASSERT(i<=ao_num) ASSERT(j<=ao_num)