2009-05-14 17:48:27 +02:00
|
|
|
BEGIN_PROVIDER [ real, eplf_gamma ]
|
|
|
|
implicit none
|
|
|
|
BEGIN_DOC
|
|
|
|
! Value of the gaussian for the EPLF
|
|
|
|
END_DOC
|
2009-10-29 18:57:46 +01:00
|
|
|
include 'constants.F'
|
2010-06-09 15:10:14 +02:00
|
|
|
real :: eps, N
|
2010-10-27 17:11:38 +02:00
|
|
|
! 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.)
|
2010-05-05 13:08:36 +02:00
|
|
|
!eplf_gamma = 1.e10
|
2010-05-28 18:23:27 +02:00
|
|
|
!eplf_gamma = 1.e5
|
2009-05-14 17:48:27 +02:00
|
|
|
END_PROVIDER
|
|
|
|
|
2009-05-15 17:18:39 +02:00
|
|
|
BEGIN_PROVIDER [ double precision, ao_eplf_integral_matrix, (ao_num,ao_num) ]
|
2009-05-14 17:48:27 +02:00
|
|
|
implicit none
|
|
|
|
BEGIN_DOC
|
|
|
|
! Array of all the <chi_i chi_j | exp(-gamma r^2)> for EPLF
|
|
|
|
END_DOC
|
|
|
|
integer :: i, j
|
2009-05-15 17:18:39 +02:00
|
|
|
double precision :: ao_eplf_integral
|
2009-05-14 17:48:27 +02:00
|
|
|
do i=1,ao_num
|
|
|
|
do j=i,ao_num
|
2010-10-06 13:27:18 +02:00
|
|
|
ao_eplf_integral_matrix(j,i) = ao_eplf_integral(j,i,dble(eplf_gamma),point)
|
2009-05-15 17:18:39 +02:00
|
|
|
ao_eplf_integral_matrix(i,j) = ao_eplf_integral_matrix(j,i)
|
|
|
|
enddo
|
|
|
|
enddo
|
|
|
|
END_PROVIDER
|
|
|
|
|
2009-05-18 15:00:54 +02:00
|
|
|
BEGIN_PROVIDER [ double precision, mo_eplf_integral_matrix, (mo_num,mo_num) ]
|
2009-05-15 17:18:39 +02:00
|
|
|
implicit none
|
|
|
|
BEGIN_DOC
|
2010-06-04 15:24:54 +02:00
|
|
|
! Array of all the <phi_i phi_j | exp(-gamma r^2)> for EPLF
|
2009-05-15 17:18:39 +02:00
|
|
|
END_DOC
|
|
|
|
integer :: i, j, k, l
|
2009-06-15 11:33:16 +02:00
|
|
|
double precision :: t
|
2009-05-15 17:43:02 +02:00
|
|
|
PROVIDE ao_eplf_integral_matrix
|
2010-12-17 11:31:05 +01:00
|
|
|
PROVIDE mo_coef mo_coef_transp
|
2009-05-18 15:00:54 +02:00
|
|
|
do i=1,mo_num
|
|
|
|
do j=i,mo_num
|
2009-10-12 17:37:07 +02:00
|
|
|
mo_eplf_integral_matrix(j,i) = 0.d0
|
2009-05-15 17:43:02 +02:00
|
|
|
enddo
|
|
|
|
|
2009-06-15 11:33:16 +02:00
|
|
|
|
|
|
|
do k=1,ao_num
|
2010-04-29 17:15:07 +02:00
|
|
|
if (mo_coef(k,i) /= 0.) then
|
2009-06-15 11:33:16 +02:00
|
|
|
do l=1,ao_num
|
2010-06-04 15:24:54 +02:00
|
|
|
if (abs(ao_eplf_integral_matrix(l,k))>1.d-16) then
|
|
|
|
t = mo_coef(k,i)*ao_eplf_integral_matrix(l,k)
|
|
|
|
do j=i,mo_num
|
|
|
|
mo_eplf_integral_matrix(j,i) += t*mo_coef_transp(j,l)
|
|
|
|
enddo
|
|
|
|
endif
|
2009-06-15 11:33:16 +02:00
|
|
|
enddo
|
|
|
|
endif
|
2010-06-04 15:24:54 +02:00
|
|
|
|
2009-05-15 17:43:02 +02:00
|
|
|
enddo
|
|
|
|
|
2009-06-15 11:33:16 +02:00
|
|
|
enddo
|
|
|
|
|
|
|
|
do i=1,mo_num
|
|
|
|
do j=i+1,mo_num
|
2009-05-15 17:18:39 +02:00
|
|
|
mo_eplf_integral_matrix(i,j) = mo_eplf_integral_matrix(j,i)
|
2009-05-14 17:48:27 +02:00
|
|
|
enddo
|
|
|
|
enddo
|
|
|
|
END_PROVIDER
|
|
|
|
|
2011-02-11 09:11:15 +01:00
|
|
|
|
2009-05-14 17:48:27 +02:00
|
|
|
BEGIN_PROVIDER [ double precision, eplf_up_up ]
|
|
|
|
&BEGIN_PROVIDER [ double precision, eplf_up_dn ]
|
|
|
|
implicit none
|
|
|
|
BEGIN_DOC
|
|
|
|
! Value of the d_upup and d_updn quantities needed for EPLF
|
|
|
|
END_DOC
|
|
|
|
|
2010-04-29 12:59:30 +02:00
|
|
|
integer :: i, j
|
2009-05-14 17:48:27 +02:00
|
|
|
|
2009-10-12 17:37:07 +02:00
|
|
|
eplf_up_up = 0.d0
|
|
|
|
eplf_up_dn = 0.d0
|
2009-05-14 17:48:27 +02:00
|
|
|
|
|
|
|
PROVIDE mo_coef_transp
|
|
|
|
|
2010-06-04 15:24:54 +02:00
|
|
|
do i=1,mo_closed_num
|
|
|
|
do j=1,mo_closed_num
|
2011-02-11 09:11:15 +01:00
|
|
|
double precision :: temp
|
|
|
|
temp = mo_value_prod_p(i,i)*mo_eplf_integral_matrix(j,j)
|
|
|
|
eplf_up_up += temp - mo_value_prod_p(j,i)*mo_eplf_integral_matrix(j,i)
|
|
|
|
eplf_up_dn += temp
|
2009-05-15 17:18:39 +02:00
|
|
|
enddo
|
|
|
|
enddo
|
2010-06-04 15:24:54 +02:00
|
|
|
eplf_up_up *= 2.d0
|
|
|
|
eplf_up_dn *= 2.d0
|
2009-05-14 17:48:27 +02:00
|
|
|
|
2011-02-11 09:11:15 +01:00
|
|
|
integer :: k,l,m
|
|
|
|
|
2011-02-14 12:04:15 +01:00
|
|
|
do m=1,two_e_density_num
|
|
|
|
i=two_e_density_indice(1,m)
|
|
|
|
j=two_e_density_indice(2,m)
|
|
|
|
k=two_e_density_indice(3,m)
|
|
|
|
l=two_e_density_indice(4,m)
|
2011-02-11 09:11:15 +01:00
|
|
|
temp = mo_value_prod_p(i,j)*mo_eplf_integral_matrix(k,l)
|
2011-02-14 12:04:15 +01:00
|
|
|
eplf_up_up += two_e_density_value(1,m)*temp
|
|
|
|
eplf_up_dn += two_e_density_value(2,m)*temp
|
2011-02-11 09:11:15 +01:00
|
|
|
enddo
|
|
|
|
|
|
|
|
END_PROVIDER
|
|
|
|
|
2009-05-14 17:48:27 +02:00
|
|
|
|
2009-11-06 00:27:24 +01:00
|
|
|
BEGIN_PROVIDER [ real, eplf_value_p ]
|
2009-05-14 17:48:27 +02:00
|
|
|
implicit none
|
|
|
|
BEGIN_DOC
|
|
|
|
! Value of the EPLF at the current point.
|
|
|
|
END_DOC
|
|
|
|
double precision :: aa, ab
|
2009-09-11 17:35:23 +02:00
|
|
|
double precision, parameter :: eps = tiny(1.d0)
|
2009-05-14 17:48:27 +02:00
|
|
|
|
|
|
|
aa = eplf_up_up
|
|
|
|
ab = eplf_up_dn
|
2009-09-11 17:35:23 +02:00
|
|
|
if ( (aa > 0.d0).and.(ab > 0.d0) ) then
|
|
|
|
aa = min(1.d0,aa)
|
|
|
|
ab = min(1.d0,ab)
|
2010-10-06 13:27:18 +02:00
|
|
|
aa = -dlog(aa)
|
2009-09-11 17:35:23 +02:00
|
|
|
aa = dsqrt(aa)
|
2010-10-27 18:16:57 +02:00
|
|
|
ab = -dlog(ab)
|
2009-09-11 17:35:23 +02:00
|
|
|
ab = dsqrt(ab)
|
2009-11-06 00:27:24 +01:00
|
|
|
eplf_value_p = (aa-ab)/(aa+ab+eps)
|
2009-09-11 17:35:23 +02:00
|
|
|
else
|
2009-11-06 00:27:24 +01:00
|
|
|
eplf_value_p = 0.d0
|
2009-09-11 17:35:23 +02:00
|
|
|
endif
|
2009-05-14 17:48:27 +02:00
|
|
|
|
|
|
|
END_PROVIDER
|
|
|
|
|
|
|
|
|
2010-10-06 13:27:18 +02:00
|
|
|
!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
|
2009-05-14 17:48:27 +02:00
|
|
|
|
2010-10-05 19:07:42 +02:00
|
|
|
|
2010-06-04 15:24:54 +02:00
|
|
|
double precision function ao_eplf_integral_primitive_oneD(a,xa,i,b,xb,j,gmma,xr)
|
2009-12-22 00:35:27 +01:00
|
|
|
implicit none
|
|
|
|
include 'constants.F'
|
|
|
|
|
2010-10-06 13:27:18 +02:00
|
|
|
real, intent(in) :: a,b ! Exponents
|
|
|
|
double precision , intent(in) :: gmma ! eplf_gamma
|
2009-12-22 00:35:27 +01:00
|
|
|
real, intent(in) :: xa,xb,xr ! Centers
|
|
|
|
integer, intent(in) :: i,j ! Powers of xa and xb
|
2010-10-27 17:11:38 +02:00
|
|
|
integer :: ii, jj, kk
|
2009-12-22 00:35:27 +01:00
|
|
|
|
|
|
|
ASSERT (a>0.)
|
|
|
|
ASSERT (b>0.)
|
|
|
|
ASSERT (i>=0)
|
|
|
|
ASSERT (j>=0)
|
|
|
|
|
2010-10-27 17:11:38 +02:00
|
|
|
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
|
2010-10-06 13:27:18 +02:00
|
|
|
|
|
|
|
if ( c > 32.d0 ) then ! Cut-off on exp(-32)
|
|
|
|
ao_eplf_integral_primitive_oneD = 0.d0
|
|
|
|
return
|
|
|
|
endif
|
|
|
|
|
2010-10-27 17:11:38 +02:00
|
|
|
double precision :: u,v,w
|
2010-10-06 13:27:18 +02:00
|
|
|
c = exp(-c)
|
2010-10-27 17:11:38 +02:00
|
|
|
w = 1.d0/sqrt(beta)
|
|
|
|
u=delta/beta-xa
|
|
|
|
v=delta/beta-xb
|
|
|
|
|
|
|
|
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)
|
2009-12-22 00:35:27 +01:00
|
|
|
enddo
|
|
|
|
enddo
|
2010-10-27 17:11:38 +02:00
|
|
|
ao_eplf_integral_primitive_oneD = dsqrt_pi*w*c*ac
|
2009-12-22 00:35:27 +01:00
|
|
|
|
2010-10-27 17:11:38 +02:00
|
|
|
end function
|
|
|
|
|
2009-05-14 17:48:27 +02:00
|
|
|
|
2009-05-15 17:18:39 +02:00
|
|
|
double precision function ao_eplf_integral(i,j,gmma,center)
|
2009-05-14 17:48:27 +02:00
|
|
|
implicit none
|
|
|
|
integer, intent(in) :: i, j
|
2010-10-06 13:27:18 +02:00
|
|
|
real, intent(in) :: center(3)
|
|
|
|
double precision, intent(in) :: gmma
|
|
|
|
!DEC$ ATTRIBUTES FORCEINLINE
|
2009-05-14 17:48:27 +02:00
|
|
|
integer :: p,q,k
|
2009-05-18 15:00:54 +02:00
|
|
|
double precision :: integral
|
2009-05-15 17:18:39 +02:00
|
|
|
double precision :: ao_eplf_integral_primitive_oneD
|
2009-12-22 00:35:27 +01:00
|
|
|
double precision :: buffer(100)
|
|
|
|
|
2009-05-14 17:48:27 +02:00
|
|
|
|
|
|
|
ASSERT(i>0)
|
|
|
|
ASSERT(j>0)
|
2010-10-06 13:27:18 +02:00
|
|
|
ASSERT(ao_prim_num_max < 100)
|
2009-05-14 17:48:27 +02:00
|
|
|
ASSERT(i<=ao_num)
|
|
|
|
ASSERT(j<=ao_num)
|
|
|
|
|
2009-12-22 00:35:27 +01:00
|
|
|
ao_eplf_integral = 0.d0
|
|
|
|
do p=1,ao_prim_num_max
|
|
|
|
buffer(p) = 0.d0
|
|
|
|
enddo
|
|
|
|
|
2009-05-14 17:48:27 +02:00
|
|
|
do q=1,ao_prim_num(j)
|
|
|
|
do p=1,ao_prim_num(i)
|
2009-05-18 15:00:54 +02:00
|
|
|
integral = &
|
|
|
|
ao_eplf_integral_primitive_oneD( &
|
2009-10-12 17:37:07 +02:00
|
|
|
ao_expo(i,p), &
|
2009-05-14 17:48:27 +02:00
|
|
|
nucl_coord(ao_nucl(i),1), &
|
|
|
|
ao_power(i,1), &
|
2009-10-12 17:37:07 +02:00
|
|
|
ao_expo(j,q), &
|
2009-05-14 17:48:27 +02:00
|
|
|
nucl_coord(ao_nucl(j),1), &
|
|
|
|
ao_power(j,1), &
|
|
|
|
gmma, &
|
2010-10-05 19:07:42 +02:00
|
|
|
center(1))
|
|
|
|
if (integral /= 0.d0) then
|
|
|
|
integral *= &
|
|
|
|
ao_eplf_integral_primitive_oneD( &
|
|
|
|
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))
|
|
|
|
if (integral /= 0.d0) then
|
|
|
|
integral *= &
|
|
|
|
ao_eplf_integral_primitive_oneD( &
|
|
|
|
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))
|
|
|
|
buffer(p) += integral*ao_coef(i,p)*ao_coef(j,q)
|
|
|
|
endif
|
|
|
|
endif
|
2009-05-14 17:48:27 +02:00
|
|
|
enddo
|
|
|
|
enddo
|
2009-12-22 00:35:27 +01:00
|
|
|
do p=1,ao_prim_num_max
|
|
|
|
ao_eplf_integral += buffer(p)
|
|
|
|
enddo
|
2009-05-14 17:48:27 +02:00
|
|
|
|
|
|
|
end function
|
|
|
|
|
2009-10-12 17:37:07 +02:00
|
|
|
double precision function mo_eplf_integral(i,j)
|
2009-05-20 17:57:35 +02:00
|
|
|
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
|
2009-12-21 23:39:00 +01:00
|
|
|
mo_eplf_integral += &
|
2009-05-20 17:57:35 +02:00
|
|
|
mo_coef(k,i)*mo_coef(l,j)*ao_eplf_integral_matrix(k,l)
|
|
|
|
enddo
|
|
|
|
endif
|
|
|
|
enddo
|
|
|
|
|
2009-10-12 17:37:07 +02:00
|
|
|
end function
|
2009-05-20 17:57:35 +02:00
|
|
|
|