mirror of https://gitlab.com/scemama/eplf
d_up_up is negative for open shell => bug
This commit is contained in:
parent
1284eab958
commit
2a64971195
|
|
@ -0,0 +1,107 @@
|
|||
BEGIN_PROVIDER [ real, ao_axis_power_p, (-2:ao_power_max,3,nucl_num) ]
|
||||
implicit none
|
||||
|
||||
BEGIN_DOC
|
||||
! Evaluation of power of x, y, z at the current point for each
|
||||
! nucleus. Negative power -> 0.
|
||||
END_DOC
|
||||
|
||||
integer :: i,k,l
|
||||
do i=1,nucl_num
|
||||
do l=1,3
|
||||
ao_axis_power_p(-2,l,i) = 0.
|
||||
ao_axis_power_p(-1,l,i) = 0.
|
||||
ao_axis_power_p(0,l,i) = 0.
|
||||
ao_axis_power_p(0,l,i) = 1.
|
||||
do k=1,ao_power_max_nucl(i,l)
|
||||
ao_axis_power_p(k,l,i) = point_nucl_dist_vec(i,l)*ao_axis_power_p(k-1,l,i)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ real, ao_axis_p, (ao_num) ]
|
||||
implicit none
|
||||
include 'types.F'
|
||||
|
||||
BEGIN_DOC
|
||||
! Cartesian polynomial part of the atomic orbitals.
|
||||
END_DOC
|
||||
integer :: i
|
||||
|
||||
do i=1,ao_num
|
||||
ao_axis_p(i) &
|
||||
= ao_axis_power_p( ao_power(i,1) , 1 , ao_nucl(i) ) &
|
||||
* ao_axis_power_p( ao_power(i,2) , 2 , ao_nucl(i) ) &
|
||||
* ao_axis_power_p( ao_power(i,3) , 3 , ao_nucl(i) )
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ real, ao_axis_grad_p, (ao_num,3) ]
|
||||
implicit none
|
||||
include 'types.F'
|
||||
|
||||
BEGIN_DOC
|
||||
! Gradients of the cartesian polynomial part of the atomic orbitals.
|
||||
END_DOC
|
||||
integer :: i, l
|
||||
real:: real_of_int(-1:10)
|
||||
data real_of_int /0.,0.,1.,2.,3.,4.,5.,6.,7.,8.,9.,10./
|
||||
|
||||
do i=1,ao_num
|
||||
ao_axis_grad_p(i,1) = real_of_int(ao_power(i,1)) &
|
||||
* ao_axis_power_p( ao_power(i,1)-1, 1 , ao_nucl(i) ) &
|
||||
* ao_axis_power_p( ao_power(i,2) , 2 , ao_nucl(i) ) &
|
||||
* ao_axis_power_p( ao_power(i,3) , 3 , ao_nucl(i) )
|
||||
enddo
|
||||
|
||||
do i=1,ao_num
|
||||
ao_axis_grad_p(i,2) = real_of_int(ao_power(i,2)) &
|
||||
* ao_axis_power_p( ao_power(i,1) , 1 , ao_nucl(i) ) &
|
||||
* ao_axis_power_p( ao_power(i,2)-1, 2 , ao_nucl(i) ) &
|
||||
* ao_axis_power_p( ao_power(i,3) , 3 , ao_nucl(i) )
|
||||
enddo
|
||||
|
||||
do i=1,ao_num
|
||||
ao_axis_grad_p(i,3) = real_of_int(ao_power(i,3)) &
|
||||
* ao_axis_power_p( ao_power(i,1) , 1 , ao_nucl(i) ) &
|
||||
* ao_axis_power_p( ao_power(i,2) , 2 , ao_nucl(i) ) &
|
||||
* ao_axis_power_p( ao_power(i,3)-1, 3 , ao_nucl(i) )
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ real, ao_axis_lapl_p, (ao_num) ]
|
||||
implicit none
|
||||
include 'types.F'
|
||||
BEGIN_DOC
|
||||
! Laplacian of the cartesian atomic orbitals
|
||||
END_DOC
|
||||
integer :: i, j, l
|
||||
|
||||
do i=1,ao_num
|
||||
real:: real_of_int(-2:10)
|
||||
data real_of_int /0.,0.,0.,1.,2.,3.,4.,5.,6.,7.,8.,9.,10./
|
||||
|
||||
ao_axis_lapl_p(i) &
|
||||
= real_of_int(ao_power(i,1)) &
|
||||
* real_of_int(ao_power(i,1)-1) &
|
||||
* ao_axis_power_p( ao_power(i,1)-2, 1 , ao_nucl(i) ) &
|
||||
* ao_axis_power_p( ao_power(i,2) , 2 , ao_nucl(i) ) &
|
||||
* ao_axis_power_p( ao_power(i,3) , 3 , ao_nucl(i) ) &
|
||||
+ real_of_int(ao_power(i,2)) &
|
||||
* real_of_int(ao_power(i,2)-1) &
|
||||
* ao_axis_power_p( ao_power(i,1) , 1 , ao_nucl(i) ) &
|
||||
* ao_axis_power_p( ao_power(i,2)-2, 2 , ao_nucl(i) ) &
|
||||
* ao_axis_power_p( ao_power(i,3) , 3 , ao_nucl(i) ) &
|
||||
+ real_of_int(ao_power(i,3)) &
|
||||
* real_of_int(ao_power(i,3)-1) &
|
||||
* ao_axis_power_p( ao_power(i,1) , 1 , ao_nucl(i) ) &
|
||||
* ao_axis_power_p( ao_power(i,2) , 2 , ao_nucl(i) ) &
|
||||
* ao_axis_power_p( ao_power(i,3)-2, 3 , ao_nucl(i) )
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
|
|
@ -0,0 +1,51 @@
|
|||
BEGIN_PROVIDER [ real, ao_value_p, (ao_num) ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Values of the atomic orbitals
|
||||
END_DOC
|
||||
integer :: i
|
||||
do i=1,ao_num
|
||||
ao_value_p(i) = ao_oneD_p(i) * ao_axis_p(i)
|
||||
enddo
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ real, ao_grad_p, (ao_num,3) ]
|
||||
implicit none
|
||||
include 'types.F'
|
||||
BEGIN_DOC
|
||||
! Gradients of the atomic orbitals
|
||||
END_DOC
|
||||
|
||||
integer :: i,l
|
||||
do l=1,3
|
||||
do i=1,ao_num
|
||||
ao_grad_p(i,l) = ao_oneD_p(i) * ao_axis_grad_p(i,l)
|
||||
enddo
|
||||
do i=1,ao_num
|
||||
ao_grad_p(i,l) = ao_grad_p(i,l) + ao_oneD_grad_p(i,l) * ao_axis_p(i)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ real, ao_lapl_p, (ao_num) ]
|
||||
implicit none
|
||||
include 'types.F'
|
||||
BEGIN_DOC
|
||||
! Laplacian of the atomic orbitals
|
||||
END_DOC
|
||||
integer :: i,l
|
||||
do i=1,ao_num
|
||||
ao_lapl_p(i) = ao_oneD_p(i) * ao_axis_lapl_p(i)
|
||||
enddo
|
||||
do i=1,ao_num
|
||||
ao_lapl_p(i) = ao_lapl_p(i) + ao_oneD_lapl_p(i) * ao_axis_p(i)
|
||||
enddo
|
||||
do l=1,3
|
||||
do i=1,ao_num
|
||||
ao_lapl_p(i) = ao_lapl_p(i) + 2.*ao_oneD_grad_p(i,l) * ao_axis_grad_p(i,l)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
|
@ -0,0 +1,128 @@
|
|||
BEGIN_PROVIDER [ real, ao_oneD_prim_p, (ao_prim_num_max,ao_num) ]
|
||||
implicit none
|
||||
include 'types.F'
|
||||
|
||||
BEGIN_DOC
|
||||
! Exponentials of the primitive AOs
|
||||
END_DOC
|
||||
integer :: i, k
|
||||
real:: r2, rtemp
|
||||
|
||||
! Compute alpha*r or alpha*r^2
|
||||
do i=1,ao_num
|
||||
r2 = point_nucl_dist_2(ao_nucl(i))
|
||||
do k=1,ao_prim_num(i)
|
||||
ao_oneD_prim_p(k,i) = r2
|
||||
enddo
|
||||
enddo
|
||||
|
||||
! Compute exp(-alpha*r) or exp(-alpha*r^2)
|
||||
do i=1,ao_num
|
||||
do k=1,ao_prim_num(i)
|
||||
ao_oneD_prim_p(k,i) = exp(-ao_oneD_prim_p(k,i)*ao_expo(k,i))
|
||||
enddo
|
||||
! Cut below 1.d-12
|
||||
do k=1,ao_prim_num(i)
|
||||
if ( abs(ao_oneD_prim_p(k,i)) < 1.e-12 ) then
|
||||
ao_oneD_prim_p(k,i) = 0.
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ real, ao_oneD_p, (ao_num) ]
|
||||
implicit none
|
||||
include 'types.F'
|
||||
|
||||
BEGIN_DOC
|
||||
! One-dimensional component of the AOs
|
||||
END_DOC
|
||||
|
||||
integer :: i, k
|
||||
|
||||
do i=1,ao_num
|
||||
ao_oneD_p(i) = 0.
|
||||
do k=1,ao_prim_num(i)
|
||||
ao_oneD_p(i) = ao_oneD_p(i) + ao_coef(k,i)*ao_oneD_prim_p(k,i)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
BEGIN_PROVIDER [ real, ao_oneD_prim_grad_p, (ao_prim_num_max,ao_num,3) ]
|
||||
implicit none
|
||||
include 'types.F'
|
||||
|
||||
BEGIN_DOC
|
||||
! Gradients of the one-dimensional component of the primitive AOs
|
||||
END_DOC
|
||||
integer :: i, k, l
|
||||
real:: factor
|
||||
do l=1,3
|
||||
do i=1,ao_num
|
||||
factor = -2.*point_nucl_dist_vec(ao_nucl(i),l)
|
||||
do k=1,ao_prim_num(i)
|
||||
ao_oneD_prim_grad_p(k,i,l) = factor*ao_expo(k,i)*ao_oneD_prim_p(k,i)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ real, ao_oneD_grad_p, (ao_num,3) ]
|
||||
implicit none
|
||||
include 'types.F'
|
||||
|
||||
BEGIN_DOC
|
||||
! Gradients of the one-dimensional component of the AOs
|
||||
END_DOC
|
||||
integer :: i, k, l
|
||||
do l=1,3
|
||||
do i=1,ao_num
|
||||
ao_oneD_grad_p(i,l) = 0.
|
||||
do k=1,ao_prim_num(i)
|
||||
ao_oneD_grad_p(i,l) = ao_oneD_grad_p(i,l) + ao_coef(k,i)*ao_oneD_prim_grad_p(k,i,l)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ real, ao_oneD_prim_lapl_p, (ao_prim_num_max,ao_num) ]
|
||||
implicit none
|
||||
include 'types.F'
|
||||
|
||||
BEGIN_DOC
|
||||
! Laplacian of the one-dimensional component of the primitive AOs
|
||||
END_DOC
|
||||
integer :: i, k, l
|
||||
do i=1,ao_num
|
||||
do k=1,ao_prim_num(i)
|
||||
ao_oneD_prim_lapl_p(k,i) = ao_oneD_prim_p(k,i) * ao_expo(k,i) * &
|
||||
( 4.*ao_expo(k,i)*point_nucl_dist_2(ao_nucl(i)) - 6. )
|
||||
enddo
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ real, ao_oneD_lapl_p, (ao_num) ]
|
||||
implicit none
|
||||
include 'types.F'
|
||||
|
||||
BEGIN_DOC
|
||||
! Laplacian of the one-dimensional component of the AOs
|
||||
END_DOC
|
||||
|
||||
integer :: i, k
|
||||
|
||||
do i=1,ao_num
|
||||
ao_oneD_lapl_p(i) = 0.
|
||||
do k=1,ao_prim_num(i)
|
||||
ao_oneD_lapl_p(i) = ao_oneD_lapl_p(i) + ao_coef(k,i)*ao_oneD_prim_lapl_p(k,i)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
|
@ -0,0 +1,54 @@
|
|||
BEGIN_PROVIDER [ integer, elec_alpha_num ]
|
||||
|
||||
BEGIN_DOC
|
||||
! Number of alpha electrons
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
elec_alpha_num = -1
|
||||
!$OMP CRITICAL (qcio_critical)
|
||||
call qcio_get_system_num_alpha(elec_alpha_num)
|
||||
!$OMP END CRITICAL (qcio_critical)
|
||||
ASSERT (elec_alpha_num > 0)
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ integer, elec_beta_num ]
|
||||
|
||||
BEGIN_DOC
|
||||
! Number of beta electrons
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
elec_beta_num = -1
|
||||
!$OMP CRITICAL (qcio_critical)
|
||||
call qcio_get_system_num_beta(elec_beta_num)
|
||||
!$OMP END CRITICAL (qcio_critical)
|
||||
ASSERT (elec_beta_num >= 0)
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ integer, elec_num ]
|
||||
|
||||
BEGIN_DOC
|
||||
! Number of electrons
|
||||
END_DOC
|
||||
implicit none
|
||||
|
||||
elec_num = elec_alpha_num + elec_beta_num
|
||||
|
||||
ASSERT ( elec_num > 0 )
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ integer, elec_num_2, (2) ]
|
||||
|
||||
BEGIN_DOC
|
||||
! Number of alpha and beta electrons in an array
|
||||
END_DOC
|
||||
|
||||
elec_num_2(1) = elec_alpha_num
|
||||
elec_num_2(2) = elec_beta_num
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
|
@ -0,0 +1,337 @@
|
|||
BEGIN_PROVIDER [ real, eplf_gamma ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Value of the gaussian for the EPLF
|
||||
END_DOC
|
||||
eplf_gamma = 10.
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ double precision, eplf_integral_matrix, (ao_num,ao_num) ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Array of all the <chi_i chi_j | exp(-gamma r^2)> for EPLF
|
||||
END_DOC
|
||||
integer :: i, j
|
||||
double precision :: eplf_integral
|
||||
do i=1,ao_num
|
||||
do j=i,ao_num
|
||||
eplf_integral_matrix(j,i) = eplf_integral(j,i,eplf_gamma,point)
|
||||
eplf_integral_matrix(i,j) = eplf_integral_matrix(j,i)
|
||||
enddo
|
||||
enddo
|
||||
END_PROVIDER
|
||||
|
||||
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
|
||||
|
||||
integer :: i, j, k, l, m, n
|
||||
|
||||
eplf_up_up = 0.
|
||||
eplf_up_dn = 0.
|
||||
|
||||
PROVIDE mo_coef_transp
|
||||
|
||||
do m=1,ao_num
|
||||
do n=1,ao_num
|
||||
double precision :: ao_mn
|
||||
ao_mn = eplf_integral_matrix(m,n)
|
||||
|
||||
if (ao_mn /= 0.d0) then
|
||||
|
||||
do k=1,ao_num
|
||||
do l=1,ao_num
|
||||
double precision :: ao_klmn
|
||||
ao_klmn = ao_mn*ao_value_p(k)*ao_value_p(l)
|
||||
|
||||
if (ao_klmn /= 0.d0) then
|
||||
|
||||
do j=1,elec_beta_num
|
||||
if (mo_coef_transp(j,n) /= 0.d0) then
|
||||
do i=1,elec_beta_num
|
||||
eplf_up_up = eplf_up_up + 2.d0*ao_klmn* &
|
||||
mo_coef_transp(i,k)*mo_coef_transp(j,n) * &
|
||||
(mo_coef_transp(i,l)*mo_coef_transp(j,m) - &
|
||||
mo_coef_transp(i,m)*mo_coef_transp(j,l) )
|
||||
enddo
|
||||
endif
|
||||
enddo
|
||||
|
||||
do j=1+elec_beta_num, elec_alpha_num
|
||||
if (mo_coef_transp(j,n) /= 0.d0) then
|
||||
do i=1,elec_beta_num
|
||||
eplf_up_up = eplf_up_up + 2.d0*ao_klmn* &
|
||||
mo_coef_transp(i,k)*mo_coef_transp(j,n) * &
|
||||
(mo_coef_transp(i,l)*mo_coef_transp(j,m) - &
|
||||
mo_coef_transp(i,m)*mo_coef_transp(j,l) )
|
||||
enddo
|
||||
do i=1+elec_beta_num,elec_alpha_num
|
||||
eplf_up_up = eplf_up_up + ao_klmn* &
|
||||
mo_coef_transp(i,k)*mo_coef_transp(j,n) * &
|
||||
(mo_coef_transp(i,l)*mo_coef_transp(j,m) - &
|
||||
mo_coef_transp(i,m)*mo_coef_transp(j,l) )
|
||||
enddo
|
||||
endif
|
||||
enddo
|
||||
|
||||
do j=1,elec_beta_num
|
||||
if ( (mo_coef_transp(j,n) /= 0.d0).and. &
|
||||
(mo_coef_transp(j,m) /= 0.d0) ) then
|
||||
do i=1,elec_alpha_num
|
||||
eplf_up_dn = eplf_up_dn + 2.d0*ao_klmn* &
|
||||
mo_coef_transp(i,k)*mo_coef_transp(j,n) * &
|
||||
mo_coef_transp(i,l)*mo_coef_transp(j,m)
|
||||
enddo
|
||||
endif
|
||||
enddo
|
||||
|
||||
endif
|
||||
|
||||
enddo
|
||||
enddo
|
||||
|
||||
endif
|
||||
|
||||
enddo
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
BEGIN_PROVIDER [ real, eplf_value ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Value of the EPLF at the current point.
|
||||
END_DOC
|
||||
double precision :: aa, ab
|
||||
|
||||
aa = eplf_up_up
|
||||
ab = eplf_up_dn
|
||||
aa = min(1.d0,aa)
|
||||
ab = min(1.d0,ab)
|
||||
aa = max(1.d-30,aa)
|
||||
ab = max(1.d-30,ab)
|
||||
aa = -dlog(aa)/eplf_gamma
|
||||
ab = -dlog(ab)/eplf_gamma
|
||||
aa = dsqrt(aa)
|
||||
ab = dsqrt(ab)
|
||||
|
||||
eplf_value = (aa-ab)/(aa+ab)
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
double precision function 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=1000
|
||||
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 = dtemp + &
|
||||
(x-xa)**i * (x-xb)**j * exp(-(a*(x-xa)**2+b*(x-xb)**2+gmma*(x-xr)**2))
|
||||
x = x+dx
|
||||
enddo
|
||||
eplf_integral_primitive_oneD_numeric = dtemp*dx
|
||||
|
||||
end function
|
||||
|
||||
double precision function eplf_integral_numeric(i,j,gmma,center)
|
||||
implicit none
|
||||
integer, intent(in) :: i, j
|
||||
integer :: p,q,k
|
||||
double precision :: integral(ao_prim_num_max,ao_prim_num_max)
|
||||
double precision :: eplf_integral_primitive_oneD_numeric
|
||||
real :: gmma, center(3)
|
||||
|
||||
do q=1,ao_prim_num(j)
|
||||
do p=1,ao_prim_num(i)
|
||||
integral(p,q) = &
|
||||
eplf_integral_primitive_oneD_numeric( &
|
||||
ao_expo(p,i), &
|
||||
nucl_coord(ao_nucl(i),1), &
|
||||
ao_power(i,1), &
|
||||
ao_expo(q,j), &
|
||||
nucl_coord(ao_nucl(j),1), &
|
||||
ao_power(j,1), &
|
||||
gmma, &
|
||||
center(1)) * &
|
||||
eplf_integral_primitive_oneD_numeric( &
|
||||
ao_expo(p,i), &
|
||||
nucl_coord(ao_nucl(i),2), &
|
||||
ao_power(i,2), &
|
||||
ao_expo(q,j), &
|
||||
nucl_coord(ao_nucl(j),2), &
|
||||
ao_power(j,2), &
|
||||
gmma, &
|
||||
center(2)) * &
|
||||
eplf_integral_primitive_oneD_numeric( &
|
||||
ao_expo(p,i), &
|
||||
nucl_coord(ao_nucl(i),3), &
|
||||
ao_power(i,3), &
|
||||
ao_expo(q,j), &
|
||||
nucl_coord(ao_nucl(j),3), &
|
||||
ao_power(j,3), &
|
||||
gmma, &
|
||||
center(3))
|
||||
enddo
|
||||
enddo
|
||||
|
||||
do q=1,ao_prim_num(j)
|
||||
do p=1,ao_prim_num(i)
|
||||
integral(p,q) = integral(p,q)*ao_coef(p,i)*ao_coef(q,j)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
eplf_integral_numeric = 0.
|
||||
do q=1,ao_prim_num(j)
|
||||
do p=1,ao_prim_num(i)
|
||||
eplf_integral_numeric = eplf_integral_numeric + integral(p,q)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
end function
|
||||
|
||||
double precision function 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 :: S(0:i+1,0:j)
|
||||
double precision :: inv_p, di(max(i,j)), dj(j), c, c1
|
||||
|
||||
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
|
||||
S(0,0) = dsqrt(pi*inv_p)*c*c1
|
||||
|
||||
! 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
|
||||
|
||||
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
|
||||
|
||||
if (j>1) then
|
||||
do jj=1,j-1
|
||||
S(0,jj+1) = (xp-xb) * S(0,jj) + di(jj)*S(0,jj-1)
|
||||
enddo
|
||||
endif
|
||||
|
||||
do jj=1,j
|
||||
S(1,jj) = (xp-xa) * S(0,jj) + di(jj) * S(0,jj-1)
|
||||
do ii=2,i
|
||||
S(ii,jj) = (xp-xa) * S(ii-1,jj) + di(ii-1) * S(ii-2,jj) + di(jj) * S(ii-1,jj-1)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
eplf_integral_primitive_oneD = S(i,j)
|
||||
|
||||
end function
|
||||
|
||||
double precision function eplf_integral(i,j,gmma,center)
|
||||
implicit none
|
||||
integer, intent(in) :: i, j
|
||||
integer :: p,q,k
|
||||
double precision :: integral(ao_prim_num_max,ao_prim_num_max)
|
||||
double precision :: eplf_integral_primitive_oneD
|
||||
real :: gmma, center(3)
|
||||
|
||||
ASSERT(i>0)
|
||||
ASSERT(j>0)
|
||||
ASSERT(i<=ao_num)
|
||||
ASSERT(j<=ao_num)
|
||||
|
||||
do q=1,ao_prim_num(j)
|
||||
do p=1,ao_prim_num(i)
|
||||
integral(p,q) = &
|
||||
eplf_integral_primitive_oneD( &
|
||||
ao_expo(p,i), &
|
||||
nucl_coord(ao_nucl(i),1), &
|
||||
ao_power(i,1), &
|
||||
ao_expo(q,j), &
|
||||
nucl_coord(ao_nucl(j),1), &
|
||||
ao_power(j,1), &
|
||||
gmma, &
|
||||
center(1)) * &
|
||||
eplf_integral_primitive_oneD( &
|
||||
ao_expo(p,i), &
|
||||
nucl_coord(ao_nucl(i),2), &
|
||||
ao_power(i,2), &
|
||||
ao_expo(q,j), &
|
||||
nucl_coord(ao_nucl(j),2), &
|
||||
ao_power(j,2), &
|
||||
gmma, &
|
||||
center(2)) * &
|
||||
eplf_integral_primitive_oneD( &
|
||||
ao_expo(p,i), &
|
||||
nucl_coord(ao_nucl(i),3), &
|
||||
ao_power(i,3), &
|
||||
ao_expo(q,j), &
|
||||
nucl_coord(ao_nucl(j),3), &
|
||||
ao_power(j,3), &
|
||||
gmma, &
|
||||
center(3))
|
||||
enddo
|
||||
enddo
|
||||
|
||||
do q=1,ao_prim_num(j)
|
||||
do p=1,ao_prim_num(i)
|
||||
integral(p,q) = integral(p,q)*ao_coef(p,i)*ao_coef(q,j)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
eplf_integral = 0.
|
||||
do q=1,ao_prim_num(j)
|
||||
do p=1,ao_prim_num(i)
|
||||
eplf_integral = eplf_integral + integral(p,q)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
end function
|
||||
|
||||
|
|
@ -0,0 +1,136 @@
|
|||
BEGIN_PROVIDER [ integer, mo_closed_num ]
|
||||
implicit none
|
||||
|
||||
BEGIN_DOC
|
||||
! Number of closed shell Molecular orbitals
|
||||
END_DOC
|
||||
|
||||
!$OMP CRITICAL (qcio_critical)
|
||||
call qcio_get_mo_num_closed(mo_closed_num)
|
||||
!$OMP END CRITICAL (qcio_critical)
|
||||
ASSERT (mo_closed_num >= 0)
|
||||
ASSERT (mo_closed_num <= elec_alpha_num)
|
||||
ASSERT (mo_closed_num <= elec_beta_num)
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ integer, mo_active_num ]
|
||||
implicit none
|
||||
|
||||
BEGIN_DOC
|
||||
! Number of active Molecular orbitals
|
||||
END_DOC
|
||||
|
||||
!$OMP CRITICAL (qcio_critical)
|
||||
call qcio_get_mo_num_active(mo_active_num)
|
||||
!$OMP END CRITICAL (qcio_critical)
|
||||
ASSERT (mo_active_num >= 0)
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ integer, mo_num ]
|
||||
implicit none
|
||||
|
||||
BEGIN_DOC
|
||||
! Number of Molecular orbitals
|
||||
END_DOC
|
||||
|
||||
mo_num = mo_closed_num + mo_active_num
|
||||
ASSERT(mo_num > 0)
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
BEGIN_PROVIDER [ real, mo_coef, (ao_num,mo_num) ]
|
||||
implicit none
|
||||
|
||||
BEGIN_DOC
|
||||
! Molecular orbital coefficients
|
||||
END_DOC
|
||||
|
||||
integer :: i, j
|
||||
double precision :: buffer(ao_num,mo_tot_num)
|
||||
|
||||
!$OMP CRITICAL (qcio_critical)
|
||||
call qcio_get_mo_matrix(buffer)
|
||||
!$OMP END CRITICAL (qcio_critical)
|
||||
do j=1,mo_num
|
||||
do i=1,ao_num
|
||||
mo_coef(i,j) = buffer(i,j)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ logical, mo_non_zero, (mo_num,ao_num) ]
|
||||
implicit none
|
||||
|
||||
BEGIN_DOC
|
||||
! Values where the mo coefficients are /= 0.
|
||||
END_DOC
|
||||
|
||||
integer :: i, j, k
|
||||
|
||||
do j=1,ao_num
|
||||
do k=1,mo_num
|
||||
mo_non_zero(k,j) = mo_coef_transp(j,k) /= 0.
|
||||
enddo
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ real, mo_coef_transp, (mo_num,ao_num) ]
|
||||
implicit none
|
||||
|
||||
BEGIN_DOC
|
||||
! Transposed array of mo coefficients
|
||||
END_DOC
|
||||
|
||||
integer :: i, j, k
|
||||
|
||||
do j=1,ao_num
|
||||
do k=1,mo_num
|
||||
mo_coef_transp(k,j) = mo_coef(j,k)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ logical, mo_is_closed, (mo_num) ]
|
||||
&BEGIN_PROVIDER [ logical, mo_is_active, (mo_num) ]
|
||||
implicit none
|
||||
|
||||
BEGIN_DOC
|
||||
! mo_is_closed : True if mo(i) is a closed-shell
|
||||
! mo_is_active : True if mo(i) is an active orbital
|
||||
END_DOC
|
||||
|
||||
character :: buffer(mo_tot_num)
|
||||
|
||||
integer :: i
|
||||
do i=1,mo_num
|
||||
if ( buffer(i) == 'c' ) then
|
||||
mo_is_closed(i) = .True.
|
||||
mo_is_active(i) = .False.
|
||||
else if ( buffer(i) == 'a' ) then
|
||||
mo_is_closed(i) = .False.
|
||||
mo_is_active(i) = .True.
|
||||
else
|
||||
mo_is_closed(i) = .False.
|
||||
mo_is_active(i) = .False.
|
||||
endif
|
||||
enddo
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ integer, mo_tot_num ]
|
||||
|
||||
BEGIN_DOC
|
||||
! Total number of MOs in the QCIO file
|
||||
END_DOC
|
||||
!$OMP CRITICAL (qcio_critical)
|
||||
call qcio_get_mo_num_orb_tot(mo_tot_num)
|
||||
!$OMP END CRITICAL (qcio_critical)
|
||||
ASSERT (mo_tot_num > 0)
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
|
@ -0,0 +1,48 @@
|
|||
BEGIN_PROVIDER [ real, point, (3) ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Coordinates of the current point
|
||||
END_DOC
|
||||
point(1) = 0.
|
||||
point(2) = 0.
|
||||
point(3) = 0.
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ real, point_nucl_dist_vec, (nucl_num,3) ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Distance vector between the current point and the nuclei
|
||||
END_DOC
|
||||
integer :: k
|
||||
do k=1,nucl_num
|
||||
point_nucl_dist_vec(k,1) = point(1)-nucl_coord(k,1)
|
||||
point_nucl_dist_vec(k,2) = point(2)-nucl_coord(k,2)
|
||||
point_nucl_dist_vec(k,3) = point(3)-nucl_coord(k,3)
|
||||
enddo
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ real, point_nucl_dist_2, (nucl_num) ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Square of the distance between the current point and the nuclei
|
||||
END_DOC
|
||||
integer :: k,l
|
||||
do k=1,nucl_num
|
||||
point_nucl_dist_2(k) = 0.
|
||||
do l=1,3
|
||||
point_nucl_dist_2(k) = point_nucl_dist_2(k) + point_nucl_dist_vec(k,l)**2
|
||||
enddo
|
||||
enddo
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ real, point_nucl_dist, (nucl_num) ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Distance between the current point and the nuclei
|
||||
END_DOC
|
||||
integer :: k
|
||||
do k=1,nucl_num
|
||||
point_nucl_dist(k) = sqrt(point_nucl_dist_2(k))
|
||||
enddo
|
||||
END_PROVIDER
|
||||
|
||||
|
|
@ -0,0 +1,18 @@
|
|||
program debug
|
||||
PROVIDE ao_prim_num_max
|
||||
read(*,*) eplf_gamma
|
||||
TOUCH eplf_gamma
|
||||
call run()
|
||||
end
|
||||
|
||||
subroutine run
|
||||
implicit none
|
||||
point(1) = 0.
|
||||
point(2) = 0.
|
||||
integer :: i
|
||||
do i=- 60,40
|
||||
point(3) = real(i)/10.
|
||||
TOUCH point
|
||||
print *, point(3), eplf_value, eplf_up_up, eplf_up_dn
|
||||
enddo
|
||||
end
|
||||
Loading…
Reference in New Issue