LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-21 11:03:29 +01:00
Code Releases Activity

cos x GTOs integ added

Browse Source
This commit is contained in:
Abdallah Ammar 2023-03-04 17:49:48 +01:00
parent 4a23b63a42
commit 3c161384cc
22 changed files with 4591 additions and 674 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

26
src/ao_basis/aos_in_r.irp.f
View File

@ -12,21 +12,21 @@ double precision function ao_value(i,r)
integer :: power_ao(3) integer :: power_ao(3)
double precision :: accu,dx,dy,dz,r2 double precision :: accu,dx,dy,dz,r2
num_ao = ao_nucl(i) num_ao = ao_nucl(i)
! power_ao(1:3)= ao_power(i,1:3) power_ao(1:3)= ao_power(i,1:3)
! center_ao(1:3) = nucl_coord(num_ao,1:3) center_ao(1:3) = nucl_coord(num_ao,1:3)
! dx = (r(1) - center_ao(1)) dx = (r(1) - center_ao(1))
! dy = (r(2) - center_ao(2)) dy = (r(2) - center_ao(2))
! dz = (r(3) - center_ao(3)) dz = (r(3) - center_ao(3))
! r2 = dx*dx + dy*dy + dz*dz r2 = dx*dx + dy*dy + dz*dz
! dx = dx**power_ao(1) dx = dx**power_ao(1)
! dy = dy**power_ao(2) dy = dy**power_ao(2)
! dz = dz**power_ao(3) dz = dz**power_ao(3)
accu = 0.d0 accu = 0.d0
! do m=1,ao_prim_num(i) do m=1,ao_prim_num(i)
! beta = ao_expo_ordered_transp(m,i) beta = ao_expo_ordered_transp(m,i)
! accu += ao_coef_normalized_ordered_transp(m,i) * dexp(-beta*r2) accu += ao_coef_normalized_ordered_transp(m,i) * dexp(-beta*r2)
! enddo enddo
ao_value = accu * dx * dy * dz ao_value = accu * dx * dy * dz
end end

1
src/ao_one_e_ints/NEED
View File

@ -1,2 +1,3 @@
ao_basis ao_basis
pseudo pseudo
cosgtos_ao_int

200
src/ao_one_e_ints/ao_overlap.irp.f
View File

@ -1,75 +1,99 @@
BEGIN_PROVIDER [ double precision, ao_overlap,(ao_num,ao_num) ]
&BEGIN_PROVIDER [ double precision, ao_overlap_x,(ao_num,ao_num) ] ! ---
&BEGIN_PROVIDER [ double precision, ao_overlap_y,(ao_num,ao_num) ]
&BEGIN_PROVIDER [ double precision, ao_overlap_z,(ao_num,ao_num) ] BEGIN_PROVIDER [ double precision, ao_overlap , (ao_num, ao_num) ]
implicit none &BEGIN_PROVIDER [ double precision, ao_overlap_x, (ao_num, ao_num) ]
&BEGIN_PROVIDER [ double precision, ao_overlap_y, (ao_num, ao_num) ]
&BEGIN_PROVIDER [ double precision, ao_overlap_z, (ao_num, ao_num) ]
BEGIN_DOC BEGIN_DOC
! Overlap between atomic basis functions: ! Overlap between atomic basis functions:
! !
! :math:`\int \chi_i(r) \chi_j(r) dr` ! :math:`\int \chi_i(r) \chi_j(r) dr`
END_DOC END_DOC
integer :: i,j,n,l
double precision :: f implicit none
integer :: dim1 integer :: i, j, n, l, dim1, power_A(3), power_B(3)
double precision :: overlap, overlap_x, overlap_y, overlap_z double precision :: overlap, overlap_x, overlap_y, overlap_z
double precision :: alpha, beta, c double precision :: alpha, beta, c
double precision :: A_center(3), B_center(3) double precision :: A_center(3), B_center(3)
integer :: power_A(3), power_B(3)
ao_overlap = 0.d0 ao_overlap = 0.d0
ao_overlap_x = 0.d0 ao_overlap_x = 0.d0
ao_overlap_y = 0.d0 ao_overlap_y = 0.d0
ao_overlap_z = 0.d0 ao_overlap_z = 0.d0
if (read_ao_integrals_overlap) then
call ezfio_get_ao_one_e_ints_ao_integrals_overlap(ao_overlap(1:ao_num, 1:ao_num)) if(read_ao_integrals_overlap) then
print *, 'AO overlap integrals read from disk'
call ezfio_get_ao_one_e_ints_ao_integrals_overlap(ao_overlap(1:ao_num, 1:ao_num))
print *, 'AO overlap integrals read from disk'
else else
dim1=100 if(use_cosgtos) then
!$OMP PARALLEL DO SCHEDULE(GUIDED) & !print*, ' use_cosgtos for ao_overlap ?', use_cosgtos
!$OMP DEFAULT(NONE) &
!$OMP PRIVATE(A_center,B_center,power_A,power_B,& do j = 1, ao_num
!$OMP overlap_x,overlap_y, overlap_z, overlap, & do i = 1, ao_num
!$OMP alpha, beta,i,j,c) & ao_overlap (i,j) = ao_overlap_cosgtos (i,j)
!$OMP SHARED(nucl_coord,ao_power,ao_prim_num, & ao_overlap_x(i,j) = ao_overlap_cosgtos_x(i,j)
!$OMP ao_overlap_x,ao_overlap_y,ao_overlap_z,ao_overlap,ao_num,ao_coef_normalized_ordered_transp,ao_nucl, & ao_overlap_y(i,j) = ao_overlap_cosgtos_y(i,j)
!$OMP ao_expo_ordered_transp,dim1) ao_overlap_z(i,j) = ao_overlap_cosgtos_z(i,j)
do j=1,ao_num enddo
A_center(1) = nucl_coord( ao_nucl(j), 1 ) enddo
A_center(2) = nucl_coord( ao_nucl(j), 2 )
A_center(3) = nucl_coord( ao_nucl(j), 3 ) else
power_A(1) = ao_power( j, 1 )
power_A(2) = ao_power( j, 2 ) dim1=100
power_A(3) = ao_power( j, 3 ) !$OMP PARALLEL DO SCHEDULE(GUIDED) &
do i= 1,ao_num !$OMP DEFAULT(NONE) &
B_center(1) = nucl_coord( ao_nucl(i), 1 ) !$OMP PRIVATE(A_center,B_center,power_A,power_B,&
B_center(2) = nucl_coord( ao_nucl(i), 2 ) !$OMP overlap_x,overlap_y, overlap_z, overlap, &
B_center(3) = nucl_coord( ao_nucl(i), 3 ) !$OMP alpha, beta,i,j,c) &
power_B(1) = ao_power( i, 1 ) !$OMP SHARED(nucl_coord,ao_power,ao_prim_num, &
power_B(2) = ao_power( i, 2 ) !$OMP ao_overlap_x,ao_overlap_y,ao_overlap_z,ao_overlap,ao_num,ao_coef_normalized_ordered_transp,ao_nucl, &
power_B(3) = ao_power( i, 3 ) !$OMP ao_expo_ordered_transp,dim1)
do n = 1,ao_prim_num(j) do j=1,ao_num
alpha = ao_expo_ordered_transp(n,j) A_center(1) = nucl_coord( ao_nucl(j), 1 )
do l = 1, ao_prim_num(i) A_center(2) = nucl_coord( ao_nucl(j), 2 )
beta = ao_expo_ordered_transp(l,i) A_center(3) = nucl_coord( ao_nucl(j), 3 )
call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,overlap_y,overlap_z,overlap,dim1) power_A(1) = ao_power( j, 1 )
c = ao_coef_normalized_ordered_transp(n,j) * ao_coef_normalized_ordered_transp(l,i) power_A(2) = ao_power( j, 2 )
ao_overlap(i,j) += c * overlap power_A(3) = ao_power( j, 3 )
if(isnan(ao_overlap(i,j)))then do i= 1,ao_num
print*,'i,j',i,j B_center(1) = nucl_coord( ao_nucl(i), 1 )
print*,'l,n',l,n B_center(2) = nucl_coord( ao_nucl(i), 2 )
print*,'c,overlap',c,overlap B_center(3) = nucl_coord( ao_nucl(i), 3 )
print*,overlap_x,overlap_y,overlap_z power_B(1) = ao_power( i, 1 )
stop power_B(2) = ao_power( i, 2 )
endif power_B(3) = ao_power( i, 3 )
ao_overlap_x(i,j) += c * overlap_x do n = 1,ao_prim_num(j)
ao_overlap_y(i,j) += c * overlap_y alpha = ao_expo_ordered_transp(n,j)
ao_overlap_z(i,j) += c * overlap_z do l = 1, ao_prim_num(i)
beta = ao_expo_ordered_transp(l,i)
call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,overlap_y,overlap_z,overlap,dim1)
c = ao_coef_normalized_ordered_transp(n,j) * ao_coef_normalized_ordered_transp(l,i)
ao_overlap(i,j) += c * overlap
if(isnan(ao_overlap(i,j)))then
print*,'i,j',i,j
print*,'l,n',l,n
print*,'c,overlap',c,overlap
print*,overlap_x,overlap_y,overlap_z
stop
endif
ao_overlap_x(i,j) += c * overlap_x
ao_overlap_y(i,j) += c * overlap_y
ao_overlap_z(i,j) += c * overlap_z
enddo
enddo
enddo enddo
enddo enddo
enddo !$OMP END PARALLEL DO
enddo
!$OMP END PARALLEL DO endif
endif endif
if (write_ao_integrals_overlap) then if (write_ao_integrals_overlap) then
call ezfio_set_ao_one_e_ints_ao_integrals_overlap(ao_overlap(1:ao_num, 1:ao_num)) call ezfio_set_ao_one_e_ints_ao_integrals_overlap(ao_overlap(1:ao_num, 1:ao_num))
print *, 'AO overlap integrals written to disk' print *, 'AO overlap integrals written to disk'
@ -77,6 +101,8 @@
END_PROVIDER END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, ao_overlap_imag, (ao_num, ao_num) ] BEGIN_PROVIDER [ double precision, ao_overlap_imag, (ao_num, ao_num) ]
implicit none implicit none
BEGIN_DOC BEGIN_DOC
@ -85,6 +111,8 @@ BEGIN_PROVIDER [ double precision, ao_overlap_imag, (ao_num, ao_num) ]
ao_overlap_imag = 0.d0 ao_overlap_imag = 0.d0
END_PROVIDER END_PROVIDER
! ---
BEGIN_PROVIDER [ complex*16, ao_overlap_complex, (ao_num, ao_num) ] BEGIN_PROVIDER [ complex*16, ao_overlap_complex, (ao_num, ao_num) ]
implicit none implicit none
BEGIN_DOC BEGIN_DOC
@ -98,41 +126,43 @@ BEGIN_PROVIDER [ complex*16, ao_overlap_complex, (ao_num, ao_num) ]
enddo enddo
END_PROVIDER END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, ao_overlap_abs, (ao_num, ao_num) ]
BEGIN_PROVIDER [ double precision, ao_overlap_abs,(ao_num,ao_num) ]
implicit none
BEGIN_DOC BEGIN_DOC
! Overlap between absolute values of atomic basis functions: ! Overlap between absolute values of atomic basis functions:
! !
! :math:`\int |\chi_i(r)| |\chi_j(r)| dr` ! :math:`\int |\chi_i(r)| |\chi_j(r)| dr`
END_DOC END_DOC
integer :: i,j,n,l
double precision :: f implicit none
integer :: dim1 integer :: i, j, n, l, dim1, power_A(3), power_B(3)
double precision :: overlap, overlap_x, overlap_y, overlap_z double precision :: overlap_x, overlap_y, overlap_z
double precision :: alpha, beta double precision :: alpha, beta
double precision :: A_center(3), B_center(3) double precision :: A_center(3), B_center(3)
integer :: power_A(3), power_B(3)
double precision :: lower_exp_val, dx double precision :: lower_exp_val, dx
if (is_periodic) then
do j=1,ao_num if(is_periodic) then
do i= 1,ao_num
ao_overlap_abs(i,j)= cdabs(ao_overlap_complex(i,j)) do j = 1, ao_num
do i = 1, ao_num
ao_overlap_abs(i,j) = cdabs(ao_overlap_complex(i,j))
enddo enddo
enddo enddo
else else
dim1=100 dim1=100
lower_exp_val = 40.d0 lower_exp_val = 40.d0
!$OMP PARALLEL DO SCHEDULE(GUIDED) & !$OMP PARALLEL DO SCHEDULE(GUIDED) &
!$OMP DEFAULT(NONE) & !$OMP DEFAULT(NONE) &
!$OMP PRIVATE(A_center,B_center,power_A,power_B, & !$OMP PRIVATE(A_center,B_center,power_A,power_B, &
!$OMP overlap_x,overlap_y, overlap_z, overlap, & !$OMP overlap_x,overlap_y, overlap_z, &
!$OMP alpha, beta,i,j,dx) & !$OMP alpha, beta,i,j,dx) &
!$OMP SHARED(nucl_coord,ao_power,ao_prim_num, & !$OMP SHARED(nucl_coord,ao_power,ao_prim_num, &
!$OMP ao_overlap_abs,ao_num,ao_coef_normalized_ordered_transp,ao_nucl,& !$OMP ao_overlap_abs,ao_num,ao_coef_normalized_ordered_transp,ao_nucl,&
!$OMP ao_expo_ordered_transp,dim1,lower_exp_val) !$OMP ao_expo_ordered_transp,dim1,lower_exp_val)
do j=1,ao_num do j=1,ao_num
A_center(1) = nucl_coord( ao_nucl(j), 1 ) A_center(1) = nucl_coord( ao_nucl(j), 1 )
A_center(2) = nucl_coord( ao_nucl(j), 2 ) A_center(2) = nucl_coord( ao_nucl(j), 2 )
@ -160,10 +190,14 @@ BEGIN_PROVIDER [ double precision, ao_overlap_abs,(ao_num,ao_num) ]
enddo enddo
enddo enddo
enddo enddo
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
endif endif
END_PROVIDER END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, S_inv,(ao_num,ao_num) ] BEGIN_PROVIDER [ double precision, S_inv,(ao_num,ao_num) ]
implicit none implicit none
BEGIN_DOC BEGIN_DOC

206
src/ao_one_e_ints/kin_ao_ints.irp.f
View File

@ -1,7 +1,10 @@
BEGIN_PROVIDER [ double precision, ao_deriv2_x,(ao_num,ao_num) ]
&BEGIN_PROVIDER [ double precision, ao_deriv2_y,(ao_num,ao_num) ] ! ---
&BEGIN_PROVIDER [ double precision, ao_deriv2_z,(ao_num,ao_num) ]
implicit none BEGIN_PROVIDER [ double precision, ao_deriv2_x, (ao_num, ao_num) ]
&BEGIN_PROVIDER [ double precision, ao_deriv2_y, (ao_num, ao_num) ]
&BEGIN_PROVIDER [ double precision, ao_deriv2_z, (ao_num, ao_num) ]
BEGIN_DOC BEGIN_DOC
! Second derivative matrix elements in the |AO| basis. ! Second derivative matrix elements in the |AO| basis.
! !
@ -11,114 +14,131 @@
! \langle \chi_i(x,y,z) | \frac{\partial^2}{\partial x^2} |\chi_j (x,y,z) \rangle ! \langle \chi_i(x,y,z) | \frac{\partial^2}{\partial x^2} |\chi_j (x,y,z) \rangle
! !
END_DOC END_DOC
integer :: i,j,n,l
double precision :: f implicit none
integer :: dim1 integer :: i, j, n, l, dim1, power_A(3), power_B(3)
double precision :: overlap, overlap_y, overlap_z double precision :: overlap, overlap_y, overlap_z
double precision :: overlap_x0, overlap_y0, overlap_z0 double precision :: overlap_x0, overlap_y0, overlap_z0
double precision :: alpha, beta, c double precision :: alpha, beta, c
double precision :: A_center(3), B_center(3) double precision :: A_center(3), B_center(3)
integer :: power_A(3), power_B(3)
double precision :: d_a_2,d_2 double precision :: d_a_2,d_2
dim1=100
! -- Dummy call to provide everything if(use_cosgtos) then
A_center(:) = 0.d0 !print*, 'use_cosgtos for ao_kinetic_integrals ?', use_cosgtos
B_center(:) = 1.d0
alpha = 1.d0
beta = .1d0
power_A = 1
power_B = 0
call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_y,d_a_2,overlap_z,overlap,dim1)
! --
!$OMP PARALLEL DO SCHEDULE(GUIDED) & do j = 1, ao_num
!$OMP DEFAULT(NONE) & do i = 1, ao_num
!$OMP PRIVATE(A_center,B_center,power_A,power_B,& ao_deriv2_x(i,j) = ao_deriv2_cosgtos_x(i,j)
!$OMP overlap_y, overlap_z, overlap, & ao_deriv2_y(i,j) = ao_deriv2_cosgtos_y(i,j)
!$OMP alpha, beta,i,j,c,d_a_2,d_2,deriv_tmp, & ao_deriv2_z(i,j) = ao_deriv2_cosgtos_z(i,j)
!$OMP overlap_x0,overlap_y0,overlap_z0) & enddo
!$OMP SHARED(nucl_coord,ao_power,ao_prim_num, & enddo
!$OMP ao_deriv2_x,ao_deriv2_y,ao_deriv2_z,ao_num,ao_coef_normalized_ordered_transp,ao_nucl, &
!$OMP ao_expo_ordered_transp,dim1)
do j=1,ao_num
A_center(1) = nucl_coord( ao_nucl(j), 1 )
A_center(2) = nucl_coord( ao_nucl(j), 2 )
A_center(3) = nucl_coord( ao_nucl(j), 3 )
power_A(1) = ao_power( j, 1 )
power_A(2) = ao_power( j, 2 )
power_A(3) = ao_power( j, 3 )
do i= 1,ao_num
ao_deriv2_x(i,j)= 0.d0
ao_deriv2_y(i,j)= 0.d0
ao_deriv2_z(i,j)= 0.d0
B_center(1) = nucl_coord( ao_nucl(i), 1 )
B_center(2) = nucl_coord( ao_nucl(i), 2 )
B_center(3) = nucl_coord( ao_nucl(i), 3 )
power_B(1) = ao_power( i, 1 )
power_B(2) = ao_power( i, 2 )
power_B(3) = ao_power( i, 3 )
do n = 1,ao_prim_num(j)
alpha = ao_expo_ordered_transp(n,j)
do l = 1, ao_prim_num(i)
beta = ao_expo_ordered_transp(l,i)
call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x0,overlap_y0,overlap_z0,overlap,dim1)
c = ao_coef_normalized_ordered_transp(n,j) * ao_coef_normalized_ordered_transp(l,i)
power_A(1) = power_A(1)-2 else
if (power_A(1)>-1) then
call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,d_a_2,overlap_y,overlap_z,overlap,dim1)
else
d_a_2 = 0.d0
endif
power_A(1) = power_A(1)+4
call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,d_2,overlap_y,overlap_z,overlap,dim1)
power_A(1) = power_A(1)-2
double precision :: deriv_tmp dim1=100
deriv_tmp = (-2.d0 * alpha * (2.d0 * power_A(1) +1.d0) * overlap_x0 &
+power_A(1) * (power_A(1)-1.d0) * d_a_2 &
+4.d0 * alpha * alpha * d_2 )*overlap_y0*overlap_z0
ao_deriv2_x(i,j) += c*deriv_tmp ! -- Dummy call to provide everything
power_A(2) = power_A(2)-2 A_center(:) = 0.d0
if (power_A(2)>-1) then B_center(:) = 1.d0
call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_y,d_a_2,overlap_z,overlap,dim1) alpha = 1.d0
else beta = .1d0
d_a_2 = 0.d0 power_A = 1
endif power_B = 0
power_A(2) = power_A(2)+4 call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_y,d_a_2,overlap_z,overlap,dim1)
call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_y,d_2,overlap_z,overlap,dim1) ! --
power_A(2) = power_A(2)-2
deriv_tmp = (-2.d0 * alpha * (2.d0 * power_A(2) +1.d0 ) * overlap_y0 & !$OMP PARALLEL DO SCHEDULE(GUIDED) &
+power_A(2) * (power_A(2)-1.d0) * d_a_2 & !$OMP DEFAULT(NONE) &
+4.d0 * alpha * alpha * d_2 )*overlap_x0*overlap_z0 !$OMP PRIVATE(A_center,B_center,power_A,power_B,&
ao_deriv2_y(i,j) += c*deriv_tmp !$OMP overlap_y, overlap_z, overlap, &
!$OMP alpha, beta,i,j,c,d_a_2,d_2,deriv_tmp, &
!$OMP overlap_x0,overlap_y0,overlap_z0) &
!$OMP SHARED(nucl_coord,ao_power,ao_prim_num, &
!$OMP ao_deriv2_x,ao_deriv2_y,ao_deriv2_z,ao_num,ao_coef_normalized_ordered_transp,ao_nucl, &
!$OMP ao_expo_ordered_transp,dim1)
do j=1,ao_num
A_center(1) = nucl_coord( ao_nucl(j), 1 )
A_center(2) = nucl_coord( ao_nucl(j), 2 )
A_center(3) = nucl_coord( ao_nucl(j), 3 )
power_A(1) = ao_power( j, 1 )
power_A(2) = ao_power( j, 2 )
power_A(3) = ao_power( j, 3 )
do i= 1,ao_num
ao_deriv2_x(i,j)= 0.d0
ao_deriv2_y(i,j)= 0.d0
ao_deriv2_z(i,j)= 0.d0
B_center(1) = nucl_coord( ao_nucl(i), 1 )
B_center(2) = nucl_coord( ao_nucl(i), 2 )
B_center(3) = nucl_coord( ao_nucl(i), 3 )
power_B(1) = ao_power( i, 1 )
power_B(2) = ao_power( i, 2 )
power_B(3) = ao_power( i, 3 )
do n = 1,ao_prim_num(j)
alpha = ao_expo_ordered_transp(n,j)
do l = 1, ao_prim_num(i)
beta = ao_expo_ordered_transp(l,i)
call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x0,overlap_y0,overlap_z0,overlap,dim1)
c = ao_coef_normalized_ordered_transp(n,j) * ao_coef_normalized_ordered_transp(l,i)
power_A(3) = power_A(3)-2 power_A(1) = power_A(1)-2
if (power_A(3)>-1) then if (power_A(1)>-1) then
call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_y,overlap_z,d_a_2,overlap,dim1) call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,d_a_2,overlap_y,overlap_z,overlap,dim1)
else else
d_a_2 = 0.d0 d_a_2 = 0.d0
endif endif
power_A(3) = power_A(3)+4 power_A(1) = power_A(1)+4
call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_y,overlap_z,d_2,overlap,dim1) call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,d_2,overlap_y,overlap_z,overlap,dim1)
power_A(3) = power_A(3)-2 power_A(1) = power_A(1)-2
deriv_tmp = (-2.d0 * alpha * (2.d0 * power_A(3) +1.d0 ) * overlap_z0 & double precision :: deriv_tmp
+power_A(3) * (power_A(3)-1.d0) * d_a_2 & deriv_tmp = (-2.d0 * alpha * (2.d0 * power_A(1) +1.d0) * overlap_x0 &
+4.d0 * alpha * alpha * d_2 )*overlap_x0*overlap_y0 +power_A(1) * (power_A(1)-1.d0) * d_a_2 &
ao_deriv2_z(i,j) += c*deriv_tmp +4.d0 * alpha * alpha * d_2 )*overlap_y0*overlap_z0
ao_deriv2_x(i,j) += c*deriv_tmp
power_A(2) = power_A(2)-2
if (power_A(2)>-1) then
call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_y,d_a_2,overlap_z,overlap,dim1)
else
d_a_2 = 0.d0
endif
power_A(2) = power_A(2)+4
call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_y,d_2,overlap_z,overlap,dim1)
power_A(2) = power_A(2)-2
deriv_tmp = (-2.d0 * alpha * (2.d0 * power_A(2) +1.d0 ) * overlap_y0 &
+power_A(2) * (power_A(2)-1.d0) * d_a_2 &
+4.d0 * alpha * alpha * d_2 )*overlap_x0*overlap_z0
ao_deriv2_y(i,j) += c*deriv_tmp
power_A(3) = power_A(3)-2
if (power_A(3)>-1) then
call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_y,overlap_z,d_a_2,overlap,dim1)
else
d_a_2 = 0.d0
endif
power_A(3) = power_A(3)+4
call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_y,overlap_z,d_2,overlap,dim1)
power_A(3) = power_A(3)-2
deriv_tmp = (-2.d0 * alpha * (2.d0 * power_A(3) +1.d0 ) * overlap_z0 &
+power_A(3) * (power_A(3)-1.d0) * d_a_2 &
+4.d0 * alpha * alpha * d_2 )*overlap_x0*overlap_y0
ao_deriv2_z(i,j) += c*deriv_tmp
enddo
enddo
enddo enddo
enddo enddo
enddo !$OMP END PARALLEL DO
enddo
!$OMP END PARALLEL DO endif
END_PROVIDER END_PROVIDER
! ---
BEGIN_PROVIDER [double precision, ao_kinetic_integrals, (ao_num,ao_num)] BEGIN_PROVIDER [double precision, ao_kinetic_integrals, (ao_num,ao_num)]
implicit none implicit none
BEGIN_DOC BEGIN_DOC

531
src/ao_one_e_ints/pot_ao_erf_ints.irp.f
View File

@ -1,3 +1,6 @@
! ---
subroutine give_all_erf_kl_ao(integrals_ao,mu_in,C_center) subroutine give_all_erf_kl_ao(integrals_ao,mu_in,C_center)
implicit none implicit none
BEGIN_DOC BEGIN_DOC
@ -15,36 +18,104 @@ subroutine give_all_erf_kl_ao(integrals_ao,mu_in,C_center)
enddo enddo
end end
! ---
double precision function NAI_pol_mult_erf_ao(i_ao, j_ao, mu_in, C_center)
double precision function NAI_pol_mult_erf_ao(i_ao,j_ao,mu_in,C_center)
implicit none
BEGIN_DOC BEGIN_DOC
!
! Computes the following integral : ! Computes the following integral :
! $\int_{-\infty}^{infty} dr \chi_i(r) \chi_j(r) \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$. ! $\int_{-\infty}^{infty} dr \chi_i(r) \chi_j(r) \frac{\erf(\mu |r - R_C|)}{|r - R_C|}$.
!
END_DOC END_DOC
integer, intent(in) :: i_ao,j_ao
implicit none
integer, intent(in) :: i_ao, j_ao
double precision, intent(in) :: mu_in, C_center(3) double precision, intent(in) :: mu_in, C_center(3)
integer :: i,j,num_A,num_B, power_A(3), power_B(3), n_pt_in
double precision :: A_center(3), B_center(3),integral, alpha,beta integer :: i, j, num_A, num_B, power_A(3), power_B(3), n_pt_in
double precision :: A_center(3), B_center(3), integral, alpha, beta
double precision :: NAI_pol_mult_erf double precision :: NAI_pol_mult_erf
num_A = ao_nucl(i_ao)
power_A(1:3)= ao_power(i_ao,1:3) num_A = ao_nucl(i_ao)
power_A(1:3) = ao_power(i_ao,1:3)
A_center(1:3) = nucl_coord(num_A,1:3) A_center(1:3) = nucl_coord(num_A,1:3)
num_B = ao_nucl(j_ao) num_B = ao_nucl(j_ao)
power_B(1:3)= ao_power(j_ao,1:3) power_B(1:3) = ao_power(j_ao,1:3)
B_center(1:3) = nucl_coord(num_B,1:3) B_center(1:3) = nucl_coord(num_B,1:3)
n_pt_in = n_pt_max_integrals n_pt_in = n_pt_max_integrals
NAI_pol_mult_erf_ao = 0.d0 NAI_pol_mult_erf_ao = 0.d0
do i = 1, ao_prim_num(i_ao) do i = 1, ao_prim_num(i_ao)
alpha = ao_expo_ordered_transp(i,i_ao) alpha = ao_expo_ordered_transp(i,i_ao)
do j = 1, ao_prim_num(j_ao) do j = 1, ao_prim_num(j_ao)
beta = ao_expo_ordered_transp(j,j_ao) beta = ao_expo_ordered_transp(j,j_ao)
integral = NAI_pol_mult_erf(A_center,B_center,power_A,power_B,alpha,beta,C_center,n_pt_in,mu_in)
NAI_pol_mult_erf_ao += integral * ao_coef_normalized_ordered_transp(j,j_ao)*ao_coef_normalized_ordered_transp(i,i_ao) integral = NAI_pol_mult_erf(A_center, B_center, power_A, power_B, alpha, beta, C_center, n_pt_in,mu_in)
NAI_pol_mult_erf_ao += integral * ao_coef_normalized_ordered_transp(j,j_ao) * ao_coef_normalized_ordered_transp(i,i_ao)
enddo enddo
enddo enddo
end end function NAI_pol_mult_erf_ao
! ---
double precision function NAI_pol_mult_erf_ao_with1s(i_ao, j_ao, beta, B_center, mu_in, C_center)
BEGIN_DOC
!
! Computes the following integral :
! $\int_{-\infty}^{infty} dr \chi_i(r) \chi_j(r) e^{-\beta (r - B_center)^2} \frac{\erf(\mu |r - R_C|)}{|r - R_C|}$.
!
END_DOC
implicit none
integer, intent(in) :: i_ao, j_ao
double precision, intent(in) :: beta, B_center(3)
double precision, intent(in) :: mu_in, C_center(3)
integer :: i, j, power_A1(3), power_A2(3), n_pt_in
double precision :: A1_center(3), A2_center(3), alpha1, alpha2, coef12, coef1, integral
double precision, external :: NAI_pol_mult_erf_with1s, NAI_pol_mult_erf_ao
ASSERT(beta .ge. 0.d0)
if(beta .lt. 1d-10) then
NAI_pol_mult_erf_ao_with1s = NAI_pol_mult_erf_ao(i_ao, j_ao, mu_in, C_center)
return
endif
power_A1(1:3) = ao_power(i_ao,1:3)
power_A2(1:3) = ao_power(j_ao,1:3)
A1_center(1:3) = nucl_coord(ao_nucl(i_ao),1:3)
A2_center(1:3) = nucl_coord(ao_nucl(j_ao),1:3)
n_pt_in = n_pt_max_integrals
NAI_pol_mult_erf_ao_with1s = 0.d0
do i = 1, ao_prim_num(i_ao)
alpha1 = ao_expo_ordered_transp (i,i_ao)
coef1 = ao_coef_normalized_ordered_transp(i,i_ao)
do j = 1, ao_prim_num(j_ao)
alpha2 = ao_expo_ordered_transp(j,j_ao)
coef12 = coef1 * ao_coef_normalized_ordered_transp(j,j_ao)
if(dabs(coef12) .lt. 1d-14) cycle
integral = NAI_pol_mult_erf_with1s( A1_center, A2_center, power_A1, power_A2, alpha1, alpha2 &
, beta, B_center, C_center, n_pt_in, mu_in )
NAI_pol_mult_erf_ao_with1s += integral * coef12
enddo
enddo
end function NAI_pol_mult_erf_ao_with1s
! ---
double precision function NAI_pol_mult_erf(A_center, B_center, power_A, power_B, alpha, beta, C_center, n_pt_in, mu_in) double precision function NAI_pol_mult_erf(A_center, B_center, power_A, power_B, alpha, beta, C_center, n_pt_in, mu_in)
@ -127,58 +198,221 @@ end function NAI_pol_mult_erf
! --- ! ---
subroutine NAI_pol_mult_erf_v(A_center, B_center, power_A, power_B, alpha, beta, C_center, LD_C, n_pt_in, mu_in, res_v, LD_resv, n_points)
double precision function NAI_pol_mult_erf_ao_with1s(i_ao, j_ao, beta, B_center, mu_in, C_center)
BEGIN_DOC BEGIN_DOC
! !
! Computes the following integral : ! Computes the following integral :
! $\int_{-\infty}^{infty} dr \chi_i(r) \chi_j(r) e^{-\beta (r - B_center)^2} \frac{\erf(\mu |r - R_C|)}{|r - R_C|}$. !
! .. math::
!
! \int dr (x-A_x)^a (x-B_x)^b \exp(-\alpha (x-A_x)^2 - \beta (x-B_x)^2 )
! \frac{\erf(\mu |r - R_C |)}{| r - R_C |}$.
! !
END_DOC END_DOC
include 'utils/constants.include.F'
implicit none implicit none
integer, intent(in) :: i_ao, j_ao
double precision, intent(in) :: beta, B_center(3)
double precision, intent(in) :: mu_in, C_center(3)
integer :: i, j, power_A1(3), power_A2(3), n_pt_in integer, intent(in) :: n_pt_in, n_points, LD_C, LD_resv
double precision :: A1_center(3), A2_center(3), alpha1, alpha2, coef12, coef1, integral integer, intent(in) :: power_A(3), power_B(3)
double precision, intent(in) :: A_center(3), B_center(3), alpha, beta, mu_in
double precision, intent(in) :: C_center(LD_C,3)
double precision, intent(out) :: res_v(LD_resv)
double precision, external :: NAI_pol_mult_erf_with1s, NAI_pol_mult_erf_ao integer :: i, n_pt, n_pt_out, ipoint
double precision :: P_center(3)
double precision :: d(0:n_pt_in), coeff, dist, const, factor
double precision :: const_factor, dist_integral
double precision :: accu, p_inv, p, rho, p_inv_2
double precision :: p_new, p_new2, coef_tmp
ASSERT(beta .ge. 0.d0) double precision :: rint
if(beta .lt. 1d-10) then
NAI_pol_mult_erf_ao_with1s = NAI_pol_mult_erf_ao(i_ao, j_ao, mu_in, C_center) res_V(1:LD_resv) = 0.d0
p = alpha + beta
p_inv = 1.d0 / p
p_inv_2 = 0.5d0 * p_inv
rho = alpha * beta * p_inv
p_new = mu_in / dsqrt(p + mu_in * mu_in)
p_new2 = p_new * p_new
coef_tmp = p * p_new2
dist = 0.d0
do i = 1, 3
P_center(i) = (alpha * A_center(i) + beta * B_center(i)) * p_inv
dist += (A_center(i) - B_center(i)) * (A_center(i) - B_center(i))
enddo
const_factor = dist * rho
if(const_factor > 80.d0) then
return
endif
factor = dexp(-const_factor)
coeff = dtwo_pi * factor * p_inv * p_new
n_pt = 2 * ( power_A(1) + power_B(1) + power_A(2) + power_B(2) + power_A(3) + power_B(3) )
if(n_pt == 0) then
do ipoint = 1, n_points
dist_integral = 0.d0
do i = 1, 3
dist_integral += (P_center(i) - C_center(ipoint,i)) * (P_center(i) - C_center(ipoint,i))
enddo
const = coef_tmp * dist_integral
res_v(ipoint) = coeff * rint(0, const)
enddo
else
do ipoint = 1, n_points
dist_integral = 0.d0
do i = 1, 3
dist_integral += (P_center(i) - C_center(ipoint,i)) * (P_center(i) - C_center(ipoint,i))
enddo
const = coef_tmp * dist_integral
do i = 0, n_pt_in
d(i) = 0.d0
enddo
call give_polynomial_mult_center_one_e_erf_opt(A_center, B_center, power_A, power_B, C_center(ipoint,1:3), n_pt_in, d, n_pt_out, p_inv_2, p_new2, P_center)
if(n_pt_out < 0) then
res_v(ipoint) = 0.d0
cycle
endif
! sum of integrals of type : int {t,[0,1]} exp-(rho.(P-Q)^2 * t^2) * t^i
accu = 0.d0
do i = 0, n_pt_out, 2
accu += d(i) * rint(i/2, const)
enddo
res_v(ipoint) = accu * coeff
enddo
endif
end subroutine NAI_pol_mult_erf_v
! ---
double precision function NAI_pol_mult_erf_with1s( A1_center, A2_center, power_A1, power_A2, alpha1, alpha2 &
, beta, B_center, C_center, n_pt_in, mu_in )
BEGIN_DOC
!
! Computes the following integral :
!
! .. math::
!
! \int dx (x - A1_x)^a_1 (x - B1_x)^a_2 \exp(-\alpha_1 (x - A1_x)^2 - \alpha_2 (x - A2_x)^2)
! \int dy (y - A1_y)^b_1 (y - B1_y)^b_2 \exp(-\alpha_1 (y - A1_y)^2 - \alpha_2 (y - A2_y)^2)
! \int dz (x - A1_z)^c_1 (z - B1_z)^c_2 \exp(-\alpha_1 (z - A1_z)^2 - \alpha_2 (z - A2_z)^2)
! \exp(-\beta (r - B)^2)
! \frac{\erf(\mu |r - R_C|)}{|r - R_C|}$.
!
END_DOC
include 'utils/constants.include.F'
implicit none
integer, intent(in) :: n_pt_in
integer, intent(in) :: power_A1(3), power_A2(3)
double precision, intent(in) :: C_center(3), A1_center(3), A2_center(3), B_center(3)
double precision, intent(in) :: alpha1, alpha2, beta, mu_in
integer :: i, n_pt, n_pt_out
double precision :: alpha12, alpha12_inv, alpha12_inv_2, rho12, A12_center(3), dist12, const_factor12
double precision :: p, p_inv, p_inv_2, rho, P_center(3), dist, const_factor
double precision :: dist_integral
double precision :: d(0:n_pt_in), coeff, const, factor
double precision :: accu
double precision :: p_new
double precision :: rint
! e^{-alpha1 (r - A1)^2} e^{-alpha2 (r - A2)^2} = e^{-K12} e^{-alpha12 (r - A12)^2}
alpha12 = alpha1 + alpha2
alpha12_inv = 1.d0 / alpha12
alpha12_inv_2 = 0.5d0 * alpha12_inv
rho12 = alpha1 * alpha2 * alpha12_inv
A12_center(1) = (alpha1 * A1_center(1) + alpha2 * A2_center(1)) * alpha12_inv
A12_center(2) = (alpha1 * A1_center(2) + alpha2 * A2_center(2)) * alpha12_inv
A12_center(3) = (alpha1 * A1_center(3) + alpha2 * A2_center(3)) * alpha12_inv
dist12 = (A1_center(1) - A2_center(1)) * (A1_center(1) - A2_center(1)) &
+ (A1_center(2) - A2_center(2)) * (A1_center(2) - A2_center(2)) &
+ (A1_center(3) - A2_center(3)) * (A1_center(3) - A2_center(3))
const_factor12 = dist12 * rho12
if(const_factor12 > 80.d0) then
NAI_pol_mult_erf_with1s = 0.d0
return return
endif endif
power_A1(1:3) = ao_power(i_ao,1:3) ! ---
power_A2(1:3) = ao_power(j_ao,1:3)
A1_center(1:3) = nucl_coord(ao_nucl(i_ao),1:3) ! e^{-K12} e^{-alpha12 (r - A12)^2} e^{-beta (r - B)^2} = e^{-K} e^{-p (r - P)^2}
A2_center(1:3) = nucl_coord(ao_nucl(j_ao),1:3) p = alpha12 + beta
p_inv = 1.d0 / p
p_inv_2 = 0.5d0 * p_inv
rho = alpha12 * beta * p_inv
P_center(1) = (alpha12 * A12_center(1) + beta * B_center(1)) * p_inv
P_center(2) = (alpha12 * A12_center(2) + beta * B_center(2)) * p_inv
P_center(3) = (alpha12 * A12_center(3) + beta * B_center(3)) * p_inv
dist = (A12_center(1) - B_center(1)) * (A12_center(1) - B_center(1)) &
+ (A12_center(2) - B_center(2)) * (A12_center(2) - B_center(2)) &
+ (A12_center(3) - B_center(3)) * (A12_center(3) - B_center(3))
n_pt_in = n_pt_max_integrals const_factor = const_factor12 + dist * rho
if(const_factor > 80.d0) then
NAI_pol_mult_erf_with1s = 0.d0
return
endif
NAI_pol_mult_erf_ao_with1s = 0.d0 dist_integral = (P_center(1) - C_center(1)) * (P_center(1) - C_center(1)) &
do i = 1, ao_prim_num(i_ao) + (P_center(2) - C_center(2)) * (P_center(2) - C_center(2)) &
alpha1 = ao_expo_ordered_transp (i,i_ao) + (P_center(3) - C_center(3)) * (P_center(3) - C_center(3))
coef1 = ao_coef_normalized_ordered_transp(i,i_ao)
do j = 1, ao_prim_num(j_ao) ! ---
alpha2 = ao_expo_ordered_transp(j,j_ao)
coef12 = coef1 * ao_coef_normalized_ordered_transp(j,j_ao)
if(dabs(coef12) .lt. 1d-14) cycle
integral = NAI_pol_mult_erf_with1s( A1_center, A2_center, power_A1, power_A2, alpha1, alpha2 & p_new = mu_in / dsqrt(p + mu_in * mu_in)
, beta, B_center, C_center, n_pt_in, mu_in ) factor = dexp(-const_factor)
coeff = dtwo_pi * factor * p_inv * p_new
NAI_pol_mult_erf_ao_with1s += integral * coef12 n_pt = 2 * ( (power_A1(1) + power_A2(1)) + (power_A1(2) + power_A2(2)) + (power_A1(3) + power_A2(3)) )
enddo const = p * dist_integral * p_new * p_new
if(n_pt == 0) then
NAI_pol_mult_erf_with1s = coeff * rint(0, const)
return
endif
do i = 0, n_pt_in
d(i) = 0.d0
enddo enddo
p_new = p_new * p_new
call give_polynomial_mult_center_one_e_erf_opt( A1_center, A2_center, power_A1, power_A2, C_center, n_pt_in, d, n_pt_out, p_inv_2, p_new, P_center)
end function NAI_pol_mult_erf_ao_with1s if(n_pt_out < 0) then
NAI_pol_mult_erf_with1s = 0.d0
return
endif
! sum of integrals of type : int {t,[0,1]} exp-(rho.(P-Q)^2 * t^2) * t^i
accu = 0.d0
do i = 0, n_pt_out, 2
accu += d(i) * rint(i/2, const)
enddo
NAI_pol_mult_erf_with1s = accu * coeff
end function NAI_pol_mult_erf_with1s
! ---
subroutine NAI_pol_mult_erf_with1s_v(A1_center, A2_center, power_A1, power_A2, alpha1, alpha2, beta, B_center, LD_B, C_center, LD_C, n_pt_in, mu_in, res_v, LD_resv, n_points) subroutine NAI_pol_mult_erf_with1s_v(A1_center, A2_center, power_A1, power_A2, alpha1, alpha2, beta, B_center, LD_B, C_center, LD_C, n_pt_in, mu_in, res_v, LD_resv, n_points)
@ -428,107 +662,6 @@ subroutine give_polynomial_mult_center_one_e_erf_opt(A_center, B_center, power_A
end subroutine give_polynomial_mult_center_one_e_erf_opt end subroutine give_polynomial_mult_center_one_e_erf_opt
! --- ! ---
subroutine NAI_pol_mult_erf_v(A_center, B_center, power_A, power_B, alpha, beta, C_center, LD_C, n_pt_in, mu_in, res_v, LD_resv, n_points)
BEGIN_DOC
!
! Computes the following integral :
!
! .. math::
!
! \int dr (x-A_x)^a (x-B_x)^b \exp(-\alpha (x-A_x)^2 - \beta (x-B_x)^2 )
! \frac{\erf(\mu |r - R_C |)}{| r - R_C |}$.
!
END_DOC
include 'utils/constants.include.F'
implicit none
integer, intent(in) :: n_pt_in, n_points, LD_C, LD_resv
integer, intent(in) :: power_A(3), power_B(3)
double precision, intent(in) :: A_center(3), B_center(3), alpha, beta, mu_in
double precision, intent(in) :: C_center(LD_C,3)
double precision, intent(out) :: res_v(LD_resv)
integer :: i, n_pt, n_pt_out, ipoint
double precision :: P_center(3)
double precision :: d(0:n_pt_in), coeff, dist, const, factor
double precision :: const_factor, dist_integral
double precision :: accu, p_inv, p, rho, p_inv_2
double precision :: p_new, p_new2, coef_tmp
double precision :: rint
res_V(1:LD_resv) = 0.d0
p = alpha + beta
p_inv = 1.d0 / p
p_inv_2 = 0.5d0 * p_inv
rho = alpha * beta * p_inv
p_new = mu_in / dsqrt(p + mu_in * mu_in)
p_new2 = p_new * p_new
coef_tmp = p * p_new2
dist = 0.d0
do i = 1, 3
P_center(i) = (alpha * A_center(i) + beta * B_center(i)) * p_inv
dist += (A_center(i) - B_center(i)) * (A_center(i) - B_center(i))
enddo
const_factor = dist * rho
if(const_factor > 80.d0) then
return
endif
factor = dexp(-const_factor)
coeff = dtwo_pi * factor * p_inv * p_new
n_pt = 2 * ( power_A(1) + power_B(1) + power_A(2) + power_B(2) + power_A(3) + power_B(3) )
if(n_pt == 0) then
do ipoint = 1, n_points
dist_integral = 0.d0
do i = 1, 3
dist_integral += (P_center(i) - C_center(ipoint,i)) * (P_center(i) - C_center(ipoint,i))
enddo
const = coef_tmp * dist_integral
res_v(ipoint) = coeff * rint(0, const)
enddo
else
do ipoint = 1, n_points
dist_integral = 0.d0
do i = 1, 3
dist_integral += (P_center(i) - C_center(ipoint,i)) * (P_center(i) - C_center(ipoint,i))
enddo
const = coef_tmp * dist_integral
do i = 0, n_pt_in
d(i) = 0.d0
enddo
call give_polynomial_mult_center_one_e_erf_opt(A_center, B_center, power_A, power_B, C_center(ipoint,1:3), n_pt_in, d, n_pt_out, p_inv_2, p_new2, P_center)
if(n_pt_out < 0) then
res_v(ipoint) = 0.d0
cycle
endif
! sum of integrals of type : int {t,[0,1]} exp-(rho.(P-Q)^2 * t^2) * t^i
accu = 0.d0
do i = 0, n_pt_out, 2
accu += d(i) * rint(i/2, const)
enddo
res_v(ipoint) = accu * coeff
enddo
endif
end subroutine NAI_pol_mult_erf_v
subroutine give_polynomial_mult_center_one_e_erf(A_center,B_center,alpha,beta,power_A,power_B,C_center,n_pt_in,d,n_pt_out,mu_in) subroutine give_polynomial_mult_center_one_e_erf(A_center,B_center,alpha,beta,power_A,power_B,C_center,n_pt_in,d,n_pt_out,mu_in)
@ -659,113 +792,3 @@ subroutine give_polynomial_mult_center_one_e_erf(A_center,B_center,alpha,beta,po
end end
double precision function NAI_pol_mult_erf_with1s( A1_center, A2_center, power_A1, power_A2, alpha1, alpha2 &
, beta, B_center, C_center, n_pt_in, mu_in )
BEGIN_DOC
!
! Computes the following integral :
!
! .. math::
!
! \int dx (x - A1_x)^a_1 (x - B1_x)^a_2 \exp(-\alpha_1 (x - A1_x)^2 - \alpha_2 (x - A2_x)^2)
! \int dy (y - A1_y)^b_1 (y - B1_y)^b_2 \exp(-\alpha_1 (y - A1_y)^2 - \alpha_2 (y - A2_y)^2)
! \int dz (x - A1_z)^c_1 (z - B1_z)^c_2 \exp(-\alpha_1 (z - A1_z)^2 - \alpha_2 (z - A2_z)^2)
! \exp(-\beta (r - B)^2)
! \frac{\erf(\mu |r - R_C|)}{|r - R_C|}$.
!
END_DOC
include 'utils/constants.include.F'
implicit none
integer, intent(in) :: n_pt_in
integer, intent(in) :: power_A1(3), power_A2(3)
double precision, intent(in) :: C_center(3), A1_center(3), A2_center(3), B_center(3)
double precision, intent(in) :: alpha1, alpha2, beta, mu_in
integer :: i, n_pt, n_pt_out
double precision :: alpha12, alpha12_inv, alpha12_inv_2, rho12, A12_center(3), dist12, const_factor12
double precision :: p, p_inv, p_inv_2, rho, P_center(3), dist, const_factor
double precision :: dist_integral
double precision :: d(0:n_pt_in), coeff, const, factor
double precision :: accu
double precision :: p_new
double precision :: rint
! e^{-alpha1 (r - A1)^2} e^{-alpha2 (r - A2)^2} = e^{-K12} e^{-alpha12 (r - A12)^2}
alpha12 = alpha1 + alpha2
alpha12_inv = 1.d0 / alpha12
alpha12_inv_2 = 0.5d0 * alpha12_inv
rho12 = alpha1 * alpha2 * alpha12_inv
A12_center(1) = (alpha1 * A1_center(1) + alpha2 * A2_center(1)) * alpha12_inv
A12_center(2) = (alpha1 * A1_center(2) + alpha2 * A2_center(2)) * alpha12_inv
A12_center(3) = (alpha1 * A1_center(3) + alpha2 * A2_center(3)) * alpha12_inv
dist12 = (A1_center(1) - A2_center(1)) * (A1_center(1) - A2_center(1)) &
+ (A1_center(2) - A2_center(2)) * (A1_center(2) - A2_center(2)) &
+ (A1_center(3) - A2_center(3)) * (A1_center(3) - A2_center(3))
const_factor12 = dist12 * rho12
if(const_factor12 > 80.d0) then
NAI_pol_mult_erf_with1s = 0.d0
return
endif
! ---
! e^{-K12} e^{-alpha12 (r - A12)^2} e^{-beta (r - B)^2} = e^{-K} e^{-p (r - P)^2}
p = alpha12 + beta
p_inv = 1.d0 / p
p_inv_2 = 0.5d0 * p_inv
rho = alpha12 * beta * p_inv
P_center(1) = (alpha12 * A12_center(1) + beta * B_center(1)) * p_inv
P_center(2) = (alpha12 * A12_center(2) + beta * B_center(2)) * p_inv
P_center(3) = (alpha12 * A12_center(3) + beta * B_center(3)) * p_inv
dist = (A12_center(1) - B_center(1)) * (A12_center(1) - B_center(1)) &
+ (A12_center(2) - B_center(2)) * (A12_center(2) - B_center(2)) &
+ (A12_center(3) - B_center(3)) * (A12_center(3) - B_center(3))
const_factor = const_factor12 + dist * rho
if(const_factor > 80.d0) then
NAI_pol_mult_erf_with1s = 0.d0
return
endif
dist_integral = (P_center(1) - C_center(1)) * (P_center(1) - C_center(1)) &
+ (P_center(2) - C_center(2)) * (P_center(2) - C_center(2)) &
+ (P_center(3) - C_center(3)) * (P_center(3) - C_center(3))
! ---
p_new = mu_in / dsqrt(p + mu_in * mu_in)
factor = dexp(-const_factor)
coeff = dtwo_pi * factor * p_inv * p_new
n_pt = 2 * ( (power_A1(1) + power_A2(1)) + (power_A1(2) + power_A2(2)) + (power_A1(3) + power_A2(3)) )
const = p * dist_integral * p_new * p_new
if(n_pt == 0) then
NAI_pol_mult_erf_with1s = coeff * rint(0, const)
return
endif
do i = 0, n_pt_in
d(i) = 0.d0
enddo
p_new = p_new * p_new
call give_polynomial_mult_center_one_e_erf_opt( A1_center, A2_center, power_A1, power_A2, C_center, n_pt_in, d, n_pt_out, p_inv_2, p_new, P_center)
if(n_pt_out < 0) then
NAI_pol_mult_erf_with1s = 0.d0
return
endif
! sum of integrals of type : int {t,[0,1]} exp-(rho.(P-Q)^2 * t^2) * t^i
accu = 0.d0
do i = 0, n_pt_out, 2
accu += d(i) * rint(i/2, const)
enddo
NAI_pol_mult_erf_with1s = accu * coeff
end function NAI_pol_mult_erf_with1s

179
src/ao_one_e_ints/pot_ao_ints.irp.f
View File

@ -1,4 +1,8 @@
! ---
BEGIN_PROVIDER [ double precision, ao_integrals_n_e, (ao_num,ao_num)] BEGIN_PROVIDER [ double precision, ao_integrals_n_e, (ao_num,ao_num)]
BEGIN_DOC BEGIN_DOC
! Nucleus-electron interaction, in the |AO| basis set. ! Nucleus-electron interaction, in the |AO| basis set.
! !
@ -6,84 +10,100 @@ BEGIN_PROVIDER [ double precision, ao_integrals_n_e, (ao_num,ao_num)]
! !
! These integrals also contain the pseudopotential integrals. ! These integrals also contain the pseudopotential integrals.
END_DOC END_DOC
implicit none implicit none
double precision :: alpha, beta, gama, delta integer :: num_A, num_B, power_A(3), power_B(3)
integer :: num_A,num_B integer :: i, j, k, l, n_pt_in, m
double precision :: A_center(3),B_center(3),C_center(3) double precision :: alpha, beta
integer :: power_A(3),power_B(3) double precision :: A_center(3),B_center(3),C_center(3)
integer :: i,j,k,l,n_pt_in,m double precision :: overlap_x,overlap_y,overlap_z,overlap,dx,NAI_pol_mult
double precision :: overlap_x,overlap_y,overlap_z,overlap,dx,NAI_pol_mult
ao_integrals_n_e = 0.d0
if (read_ao_integrals_n_e) then if (read_ao_integrals_n_e) then
call ezfio_get_ao_one_e_ints_ao_integrals_n_e(ao_integrals_n_e) call ezfio_get_ao_one_e_ints_ao_integrals_n_e(ao_integrals_n_e)
print *, 'AO N-e integrals read from disk' print *, 'AO N-e integrals read from disk'
else else
ao_integrals_n_e = 0.d0 if(use_cosgtos) then
!print *, " use_cosgtos for ao_integrals_n_e ?", use_cosgtos
! _ do j = 1, ao_num
! /| / |_) do i = 1, ao_num
! | / | \ ao_integrals_n_e(i,j) = ao_integrals_n_e_cosgtos(i,j)
! enddo
enddo
!$OMP PARALLEL & else
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,l,m,alpha,beta,A_center,B_center,C_center,power_A,power_B,&
!$OMP num_A,num_B,Z,c,n_pt_in) &
!$OMP SHARED (ao_num,ao_prim_num,ao_expo_ordered_transp,ao_power,ao_nucl,nucl_coord,ao_coef_normalized_ordered_transp,&
!$OMP n_pt_max_integrals,ao_integrals_n_e,nucl_num,nucl_charge)
n_pt_in = n_pt_max_integrals !$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,l,m,alpha,beta,A_center,B_center,C_center,power_A,power_B,&
!$OMP num_A,num_B,Z,c,c1,n_pt_in) &
!$OMP SHARED (ao_num,ao_prim_num,ao_expo_ordered_transp,ao_power,ao_nucl,nucl_coord,ao_coef_normalized_ordered_transp,&
!$OMP n_pt_max_integrals,ao_integrals_n_e,nucl_num,nucl_charge)
!$OMP DO SCHEDULE (dynamic) n_pt_in = n_pt_max_integrals
do j = 1, ao_num !$OMP DO SCHEDULE (dynamic)
num_A = ao_nucl(j)
power_A(1:3)= ao_power(j,1:3)
A_center(1:3) = nucl_coord(num_A,1:3)
do i = 1, ao_num do j = 1, ao_num
num_A = ao_nucl(j)
power_A(1:3)= ao_power(j,1:3)
A_center(1:3) = nucl_coord(num_A,1:3)
num_B = ao_nucl(i) do i = 1, ao_num
power_B(1:3)= ao_power(i,1:3)
B_center(1:3) = nucl_coord(num_B,1:3)
do l=1,ao_prim_num(j) num_B = ao_nucl(i)
alpha = ao_expo_ordered_transp(l,j) power_B(1:3)= ao_power(i,1:3)
B_center(1:3) = nucl_coord(num_B,1:3)
do m=1,ao_prim_num(i) do l=1,ao_prim_num(j)
beta = ao_expo_ordered_transp(m,i) alpha = ao_expo_ordered_transp(l,j)
double precision :: c do m=1,ao_prim_num(i)
c = 0.d0 beta = ao_expo_ordered_transp(m,i)
do k = 1, nucl_num double precision :: c, c1
double precision :: Z c = 0.d0
Z = nucl_charge(k)
C_center(1:3) = nucl_coord(k,1:3) do k = 1, nucl_num
double precision :: Z
Z = nucl_charge(k)
c = c - Z * NAI_pol_mult(A_center,B_center, & C_center(1:3) = nucl_coord(k,1:3)
power_A,power_B,alpha,beta,C_center,n_pt_in)
!print *, ' '
!print *, A_center, B_center, C_center, power_A, power_B
!print *, alpha, beta
c1 = NAI_pol_mult( A_center, B_center, power_A, power_B &
, alpha, beta, C_center, n_pt_in )
!print *, ' c1 = ', c1
c = c - Z * c1
enddo
ao_integrals_n_e(i,j) = ao_integrals_n_e(i,j) &
+ ao_coef_normalized_ordered_transp(l,j) &
* ao_coef_normalized_ordered_transp(m,i) * c
enddo enddo
ao_integrals_n_e(i,j) = ao_integrals_n_e(i,j) &
+ ao_coef_normalized_ordered_transp(l,j) &
* ao_coef_normalized_ordered_transp(m,i) * c
enddo enddo
enddo enddo
enddo enddo
enddo
!$OMP END DO !$OMP END DO
!$OMP END PARALLEL !$OMP END PARALLEL
IF (DO_PSEUDO) THEN
endif
IF(do_pseudo) THEN
ao_integrals_n_e += ao_pseudo_integrals ao_integrals_n_e += ao_pseudo_integrals
ENDIF ENDIF
IF(point_charges) THEN
ao_integrals_n_e += ao_integrals_pt_chrg
ENDIF
endif endif
@ -102,7 +122,7 @@ BEGIN_PROVIDER [ double precision, ao_integrals_n_e_imag, (ao_num,ao_num)]
! :math:`\langle \chi_i | -\sum_A \frac{1}{|r-R_A|} | \chi_j \rangle` ! :math:`\langle \chi_i | -\sum_A \frac{1}{|r-R_A|} | \chi_j \rangle`
END_DOC END_DOC
implicit none implicit none
double precision :: alpha, beta, gama, delta double precision :: alpha, beta
integer :: num_A,num_B integer :: num_A,num_B
double precision :: A_center(3),B_center(3),C_center(3) double precision :: A_center(3),B_center(3),C_center(3)
integer :: power_A(3),power_B(3) integer :: power_A(3),power_B(3)
@ -125,7 +145,7 @@ BEGIN_PROVIDER [ double precision, ao_integrals_n_e_per_atom, (ao_num,ao_num,nuc
! :math:`\langle \chi_i | -\frac{1}{|r-R_A|} | \chi_j \rangle` ! :math:`\langle \chi_i | -\frac{1}{|r-R_A|} | \chi_j \rangle`
END_DOC END_DOC
implicit none implicit none
double precision :: alpha, beta, gama, delta double precision :: alpha, beta
integer :: i_c,num_A,num_B integer :: i_c,num_A,num_B
double precision :: A_center(3),B_center(3),C_center(3) double precision :: A_center(3),B_center(3),C_center(3)
integer :: power_A(3),power_B(3) integer :: power_A(3),power_B(3)
@ -268,6 +288,7 @@ double precision function NAI_pol_mult(A_center,B_center,power_A,power_B,alpha,b
end end
! ---
subroutine give_polynomial_mult_center_one_e(A_center,B_center,alpha,beta,power_A,power_B,C_center,n_pt_in,d,n_pt_out) subroutine give_polynomial_mult_center_one_e(A_center,B_center,alpha,beta,power_A,power_B,C_center,n_pt_in,d,n_pt_out)
implicit none implicit none
@ -579,61 +600,3 @@ double precision function V_r(n,alpha)
end end
double precision function V_phi(n,m)
implicit none
BEGIN_DOC
! Computes the angular $\phi$ part of the nuclear attraction integral:
!
! $\int_{0}^{2 \pi} \cos(\phi)^n \sin(\phi)^m d\phi$.
END_DOC
integer :: n,m, i
double precision :: prod, Wallis
prod = 1.d0
do i = 0,shiftr(n,1)-1
prod = prod/ (1.d0 + dfloat(m+1)/dfloat(n-i-i-1))
enddo
V_phi = 4.d0 * prod * Wallis(m)
end
double precision function V_theta(n,m)
implicit none
BEGIN_DOC
! Computes the angular $\theta$ part of the nuclear attraction integral:
!
! $\int_{0}^{\pi} \cos(\theta)^n \sin(\theta)^m d\theta$
END_DOC
integer :: n,m,i
double precision :: Wallis, prod
include 'utils/constants.include.F'
V_theta = 0.d0
prod = 1.d0
do i = 0,shiftr(n,1)-1
prod = prod / (1.d0 + dfloat(m+1)/dfloat(n-i-i-1))
enddo
V_theta = (prod+prod) * Wallis(m)
end
double precision function Wallis(n)
implicit none
BEGIN_DOC
! Wallis integral:
!
! $\int_{0}^{\pi} \cos(\theta)^n d\theta$.
END_DOC
double precision :: fact
integer :: n,p
include 'utils/constants.include.F'
if(iand(n,1).eq.0)then
Wallis = fact(shiftr(n,1))
Wallis = pi * fact(n) / (dble(ibset(0_8,n)) * (Wallis+Wallis)*Wallis)
else
p = shiftr(n,1)
Wallis = fact(p)
Wallis = dble(ibset(0_8,p+p)) * Wallis*Wallis / fact(p+p+1)
endif
end

7
src/ao_two_e_ints/EZFIO.cfg
View File

@ -4,13 +4,6 @@ doc: Read/Write |AO| integrals from/to disk [ Write | Read | None ]
interface: ezfio,provider,ocaml interface: ezfio,provider,ocaml
default: None default: None
[ao_integrals_threshold]
type: Threshold
doc: If | (pq|rs) | < `ao_integrals_threshold` then (pq|rs) is zero
interface: ezfio,provider,ocaml
default: 1.e-15
ezfio_name: threshold_ao
[do_direct_integrals] [do_direct_integrals]
type: logical type: logical
doc: Compute integrals on the fly (very slow, only for debugging) doc: Compute integrals on the fly (very slow, only for debugging)

195
src/ao_two_e_ints/two_e_integrals.irp.f
View File

@ -1,102 +1,123 @@
double precision function ao_two_e_integral(i,j,k,l)
implicit none ! ---
double precision function ao_two_e_integral(i, j, k, l)
BEGIN_DOC BEGIN_DOC
! integral of the AO basis <ik|jl> or (ij|kl) ! integral of the AO basis <ik|jl> or (ij|kl)
! i(r1) j(r1) 1/r12 k(r2) l(r2) ! i(r1) j(r1) 1/r12 k(r2) l(r2)
END_DOC END_DOC
integer,intent(in) :: i,j,k,l implicit none
integer :: p,q,r,s
double precision :: I_center(3),J_center(3),K_center(3),L_center(3)
integer :: num_i,num_j,num_k,num_l,dim1,I_power(3),J_power(3),K_power(3),L_power(3)
double precision :: integral
include 'utils/constants.include.F' include 'utils/constants.include.F'
integer, intent(in) :: i, j, k, l
integer :: p, q, r, s
integer :: num_i,num_j,num_k,num_l,dim1,I_power(3),J_power(3),K_power(3),L_power(3)
integer :: iorder_p(3), iorder_q(3)
double precision :: I_center(3), J_center(3), K_center(3), L_center(3)
double precision :: integral
double precision :: P_new(0:max_dim,3),P_center(3),fact_p,pp double precision :: P_new(0:max_dim,3),P_center(3),fact_p,pp
double precision :: Q_new(0:max_dim,3),Q_center(3),fact_q,qq double precision :: Q_new(0:max_dim,3),Q_center(3),fact_q,qq
integer :: iorder_p(3), iorder_q(3)
double precision :: ao_two_e_integral_schwartz_accel double precision :: ao_two_e_integral_schwartz_accel
if (ao_prim_num(i) * ao_prim_num(j) * ao_prim_num(k) * ao_prim_num(l) > 1024 ) then double precision :: ao_two_e_integral_cosgtos
ao_two_e_integral = ao_two_e_integral_schwartz_accel(i,j,k,l)
else
dim1 = n_pt_max_integrals
num_i = ao_nucl(i) if(use_cosgtos) then
num_j = ao_nucl(j) !print *, ' use_cosgtos for ao_two_e_integral ?', use_cosgtos
num_k = ao_nucl(k)
num_l = ao_nucl(l)
ao_two_e_integral = 0.d0
if (num_i /= num_j .or. num_k /= num_l .or. num_j /= num_k)then ao_two_e_integral = ao_two_e_integral_cosgtos(i, j, k, l)
do p = 1, 3
I_power(p) = ao_power(i,p)
J_power(p) = ao_power(j,p)
K_power(p) = ao_power(k,p)
L_power(p) = ao_power(l,p)
I_center(p) = nucl_coord(num_i,p)
J_center(p) = nucl_coord(num_j,p)
K_center(p) = nucl_coord(num_k,p)
L_center(p) = nucl_coord(num_l,p)
enddo
double precision :: coef1, coef2, coef3, coef4 else
double precision :: p_inv,q_inv
double precision :: general_primitive_integral
do p = 1, ao_prim_num(i) if (ao_prim_num(i) * ao_prim_num(j) * ao_prim_num(k) * ao_prim_num(l) > 1024 ) then
coef1 = ao_coef_normalized_ordered_transp(p,i)
do q = 1, ao_prim_num(j) ao_two_e_integral = ao_two_e_integral_schwartz_accel(i,j,k,l)
coef2 = coef1*ao_coef_normalized_ordered_transp(q,j)
call give_explicit_poly_and_gaussian(P_new,P_center,pp,fact_p,iorder_p,&
ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j), &
I_power,J_power,I_center,J_center,dim1)
p_inv = 1.d0/pp
do r = 1, ao_prim_num(k)
coef3 = coef2*ao_coef_normalized_ordered_transp(r,k)
do s = 1, ao_prim_num(l)
coef4 = coef3*ao_coef_normalized_ordered_transp(s,l)
call give_explicit_poly_and_gaussian(Q_new,Q_center,qq,fact_q,iorder_q,&
ao_expo_ordered_transp(r,k),ao_expo_ordered_transp(s,l), &
K_power,L_power,K_center,L_center,dim1)
q_inv = 1.d0/qq
integral = general_primitive_integral(dim1, &
P_new,P_center,fact_p,pp,p_inv,iorder_p, &
Q_new,Q_center,fact_q,qq,q_inv,iorder_q)
ao_two_e_integral = ao_two_e_integral + coef4 * integral
enddo ! s
enddo ! r
enddo ! q
enddo ! p
else else
do p = 1, 3 dim1 = n_pt_max_integrals
I_power(p) = ao_power(i,p)
J_power(p) = ao_power(j,p)
K_power(p) = ao_power(k,p)
L_power(p) = ao_power(l,p)
enddo
double precision :: ERI
do p = 1, ao_prim_num(i) num_i = ao_nucl(i)
coef1 = ao_coef_normalized_ordered_transp(p,i) num_j = ao_nucl(j)
do q = 1, ao_prim_num(j) num_k = ao_nucl(k)
coef2 = coef1*ao_coef_normalized_ordered_transp(q,j) num_l = ao_nucl(l)
do r = 1, ao_prim_num(k) ao_two_e_integral = 0.d0
coef3 = coef2*ao_coef_normalized_ordered_transp(r,k)
do s = 1, ao_prim_num(l) if (num_i /= num_j .or. num_k /= num_l .or. num_j /= num_k)then
coef4 = coef3*ao_coef_normalized_ordered_transp(s,l) do p = 1, 3
integral = ERI( & I_power(p) = ao_power(i,p)
ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j),ao_expo_ordered_transp(r,k),ao_expo_ordered_transp(s,l),& J_power(p) = ao_power(j,p)
I_power(1),J_power(1),K_power(1),L_power(1), & K_power(p) = ao_power(k,p)
I_power(2),J_power(2),K_power(2),L_power(2), & L_power(p) = ao_power(l,p)
I_power(3),J_power(3),K_power(3),L_power(3)) I_center(p) = nucl_coord(num_i,p)
ao_two_e_integral = ao_two_e_integral + coef4 * integral J_center(p) = nucl_coord(num_j,p)
enddo ! s K_center(p) = nucl_coord(num_k,p)
enddo ! r L_center(p) = nucl_coord(num_l,p)
enddo ! q enddo
enddo ! p
double precision :: coef1, coef2, coef3, coef4
double precision :: p_inv,q_inv
double precision :: general_primitive_integral
do p = 1, ao_prim_num(i)
coef1 = ao_coef_normalized_ordered_transp(p,i)
do q = 1, ao_prim_num(j)
coef2 = coef1*ao_coef_normalized_ordered_transp(q,j)
call give_explicit_poly_and_gaussian(P_new,P_center,pp,fact_p,iorder_p,&
ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j), &
I_power,J_power,I_center,J_center,dim1)
p_inv = 1.d0/pp
do r = 1, ao_prim_num(k)
coef3 = coef2*ao_coef_normalized_ordered_transp(r,k)
do s = 1, ao_prim_num(l)
coef4 = coef3*ao_coef_normalized_ordered_transp(s,l)
call give_explicit_poly_and_gaussian(Q_new,Q_center,qq,fact_q,iorder_q,&
ao_expo_ordered_transp(r,k),ao_expo_ordered_transp(s,l), &
K_power,L_power,K_center,L_center,dim1)
q_inv = 1.d0/qq
integral = general_primitive_integral(dim1, &
P_new,P_center,fact_p,pp,p_inv,iorder_p, &
Q_new,Q_center,fact_q,qq,q_inv,iorder_q)
ao_two_e_integral = ao_two_e_integral + coef4 * integral
enddo ! s
enddo ! r
enddo ! q
enddo ! p
else
do p = 1, 3
I_power(p) = ao_power(i,p)
J_power(p) = ao_power(j,p)
K_power(p) = ao_power(k,p)
L_power(p) = ao_power(l,p)
enddo
double precision :: ERI
do p = 1, ao_prim_num(i)
coef1 = ao_coef_normalized_ordered_transp(p,i)
do q = 1, ao_prim_num(j)
coef2 = coef1*ao_coef_normalized_ordered_transp(q,j)
do r = 1, ao_prim_num(k)
coef3 = coef2*ao_coef_normalized_ordered_transp(r,k)
do s = 1, ao_prim_num(l)
coef4 = coef3*ao_coef_normalized_ordered_transp(s,l)
integral = ERI( &
ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j),ao_expo_ordered_transp(r,k),ao_expo_ordered_transp(s,l),&
I_power(1),J_power(1),K_power(1),L_power(1), &
I_power(2),J_power(2),K_power(2),L_power(2), &
I_power(3),J_power(3),K_power(3),L_power(3))
ao_two_e_integral = ao_two_e_integral + coef4 * integral
enddo ! s
enddo ! r
enddo ! q
enddo ! p
endif
endif endif
@ -104,6 +125,8 @@ double precision function ao_two_e_integral(i,j,k,l)
end end
! ---
double precision function ao_two_e_integral_schwartz_accel(i,j,k,l) double precision function ao_two_e_integral_schwartz_accel(i,j,k,l)
implicit none implicit none
BEGIN_DOC BEGIN_DOC
@ -421,14 +444,17 @@ BEGIN_PROVIDER [ logical, ao_two_e_integrals_in_map ]
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_two_e_integral_schwartz,(ao_num,ao_num) ] ! ---
implicit none
BEGIN_PROVIDER [ double precision, ao_two_e_integral_schwartz, (ao_num, ao_num) ]
BEGIN_DOC BEGIN_DOC
! Needed to compute Schwartz inequalities ! Needed to compute Schwartz inequalities
END_DOC END_DOC
integer :: i,k implicit none
double precision :: ao_two_e_integral,cpu_1,cpu_2, wall_1, wall_2 integer :: i, k
double precision :: ao_two_e_integral,cpu_1,cpu_2, wall_1, wall_2
ao_two_e_integral_schwartz(1,1) = ao_two_e_integral(1,1,1,1) ao_two_e_integral_schwartz(1,1) = ao_two_e_integral(1,1,1,1)
!$OMP PARALLEL DO PRIVATE(i,k) & !$OMP PARALLEL DO PRIVATE(i,k) &
@ -445,6 +471,7 @@ BEGIN_PROVIDER [ double precision, ao_two_e_integral_schwartz,(ao_num,ao_num) ]
END_PROVIDER END_PROVIDER
! ---
double precision function general_primitive_integral(dim, & double precision function general_primitive_integral(dim, &
P_new,P_center,fact_p,p,p_inv,iorder_p, & P_new,P_center,fact_p,p,p_inv,iorder_p, &

19
src/cosgtos_ao_int/EZFIO.cfg Normal file
View File

@ -0,0 +1,19 @@
[ao_expoim_cosgtos]
type: double precision
doc: imag part for Exponents for each primitive of each cosGTOs |AO|
size: (ao_basis.ao_num,ao_basis.ao_prim_num_max)
interface: ezfio, provider
[use_cosgtos]
type: logical
doc: If true, use cosgtos for AO integrals
interface: ezfio,provider,ocaml
default: False
[ao_integrals_threshold]
type: Threshold
doc: If | (pq|rs) | < `ao_integrals_threshold` then (pq|rs) is zero
interface: ezfio,provider,ocaml
default: 1.e-15
ezfio_name: threshold_ao

4
src/cosgtos_ao_int/README.rst Normal file
View File

@ -0,0 +1,4 @@
==============
cosgtos_ao_int
==============

210
src/cosgtos_ao_int/aos_cosgtos.irp.f Normal file
View File

@ -0,0 +1,210 @@
! ---
BEGIN_PROVIDER [ double precision, ao_coef_norm_ord_transp_cosgtos, (ao_prim_num_max, ao_num) ]
implicit none
integer :: i, j
do j = 1, ao_num
do i = 1, ao_prim_num_max
ao_coef_norm_ord_transp_cosgtos(i,j) = ao_coef_norm_ord_cosgtos(j,i)
enddo
enddo
END_PROVIDER
! ---
BEGIN_PROVIDER [ complex*16, ao_expo_ord_transp_cosgtos, (ao_prim_num_max, ao_num) ]
implicit none
integer :: i, j
do j = 1, ao_num
do i = 1, ao_prim_num_max
ao_expo_ord_transp_cosgtos(i,j) = ao_expo_ord_cosgtos(j,i)
enddo
enddo
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, ao_coef_norm_cosgtos, (ao_num, ao_prim_num_max) ]
implicit none
integer :: i, j, powA(3), nz
double precision :: norm
complex*16 :: overlap_x, overlap_y, overlap_z, C_A(3)
complex*16 :: integ1, integ2, expo
nz = 100
C_A(1) = (0.d0, 0.d0)
C_A(2) = (0.d0, 0.d0)
C_A(3) = (0.d0, 0.d0)
ao_coef_norm_cosgtos = 0.d0
do i = 1, ao_num
powA(1) = ao_power(i,1)
powA(2) = ao_power(i,2)
powA(3) = ao_power(i,3)
! Normalization of the primitives
if(primitives_normalized) then
do j = 1, ao_prim_num(i)
expo = ao_expo(i,j) + (0.d0, 1.d0) * ao_expoim_cosgtos(i,j)
call overlap_cgaussian_xyz(C_A, C_A, expo, expo, powA, powA, overlap_x, overlap_y, overlap_z, integ1, nz)
call overlap_cgaussian_xyz(C_A, C_A, conjg(expo), expo, powA, powA, overlap_x, overlap_y, overlap_z, integ2, nz)
norm = 2.d0 * real( integ1 + integ2 )
ao_coef_norm_cosgtos(i,j) = ao_coef(i,j) / dsqrt(norm)
enddo
else
do j = 1, ao_prim_num(i)
ao_coef_norm_cosgtos(i,j) = ao_coef(i,j)
enddo
endif
enddo
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, ao_coef_norm_ord_cosgtos, (ao_num, ao_prim_num_max) ]
&BEGIN_PROVIDER [ complex*16 , ao_expo_ord_cosgtos, (ao_num, ao_prim_num_max) ]
implicit none
integer :: i, j
integer :: iorder(ao_prim_num_max)
double precision :: d(ao_prim_num_max,3)
d = 0.d0
do i = 1, ao_num
do j = 1, ao_prim_num(i)
iorder(j) = j
d(j,1) = ao_expo(i,j)
d(j,2) = ao_coef_norm_cosgtos(i,j)
d(j,3) = ao_expoim_cosgtos(i,j)
enddo
call dsort (d(1,1), iorder, ao_prim_num(i))
call dset_order(d(1,2), iorder, ao_prim_num(i))
call dset_order(d(1,3), iorder, ao_prim_num(i))
do j = 1, ao_prim_num(i)
ao_expo_ord_cosgtos (i,j) = d(j,1) + (0.d0, 1.d0) * d(j,3)
ao_coef_norm_ord_cosgtos(i,j) = d(j,2)
enddo
enddo
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, ao_overlap_cosgtos, (ao_num, ao_num) ]
&BEGIN_PROVIDER [ double precision, ao_overlap_cosgtos_x, (ao_num, ao_num) ]
&BEGIN_PROVIDER [ double precision, ao_overlap_cosgtos_y, (ao_num, ao_num) ]
&BEGIN_PROVIDER [ double precision, ao_overlap_cosgtos_z, (ao_num, ao_num) ]
implicit none
integer :: i, j, n, l, dim1, power_A(3), power_B(3)
double precision :: c, overlap, overlap_x, overlap_y, overlap_z
complex*16 :: alpha, beta, A_center(3), B_center(3)
complex*16 :: overlap1, overlap_x1, overlap_y1, overlap_z1
complex*16 :: overlap2, overlap_x2, overlap_y2, overlap_z2
ao_overlap_cosgtos = 0.d0
ao_overlap_cosgtos_x = 0.d0
ao_overlap_cosgtos_y = 0.d0
ao_overlap_cosgtos_z = 0.d0
dim1 = 100
!$OMP PARALLEL DO SCHEDULE(GUIDED) &
!$OMP DEFAULT(NONE) &
!$OMP PRIVATE( A_center, B_center, power_A, power_B, alpha, beta, i, j, n, l, c &
!$OMP , overlap_x , overlap_y , overlap_z , overlap &
!$OMP , overlap_x1, overlap_y1, overlap_z1, overlap1 &
!$OMP , overlap_x2, overlap_y2, overlap_z2, overlap2 ) &
!$OMP SHARED( nucl_coord, ao_power, ao_prim_num, ao_num, ao_nucl, dim1 &
!$OMP , ao_overlap_cosgtos_x, ao_overlap_cosgtos_y, ao_overlap_cosgtos_z, ao_overlap_cosgtos &
!$OMP , ao_coef_norm_ord_transp_cosgtos, ao_expo_ord_transp_cosgtos )
do j = 1, ao_num
A_center(1) = nucl_coord(ao_nucl(j),1) * (1.d0, 0.d0)
A_center(2) = nucl_coord(ao_nucl(j),2) * (1.d0, 0.d0)
A_center(3) = nucl_coord(ao_nucl(j),3) * (1.d0, 0.d0)
power_A(1) = ao_power(j,1)
power_A(2) = ao_power(j,2)
power_A(3) = ao_power(j,3)
do i = 1, ao_num
B_center(1) = nucl_coord(ao_nucl(i),1) * (1.d0, 0.d0)
B_center(2) = nucl_coord(ao_nucl(i),2) * (1.d0, 0.d0)
B_center(3) = nucl_coord(ao_nucl(i),3) * (1.d0, 0.d0)
power_B(1) = ao_power(i,1)
power_B(2) = ao_power(i,2)
power_B(3) = ao_power(i,3)
do n = 1, ao_prim_num(j)
alpha = ao_expo_ord_transp_cosgtos(n,j)
do l = 1, ao_prim_num(i)
c = ao_coef_norm_ord_transp_cosgtos(n,j) * ao_coef_norm_ord_transp_cosgtos(l,i)
beta = ao_expo_ord_transp_cosgtos(l,i)
call overlap_cgaussian_xyz( A_center, B_center, alpha, beta, power_A, power_B &
, overlap_x1, overlap_y1, overlap_z1, overlap1, dim1 )
call overlap_cgaussian_xyz( A_center, B_center, conjg(alpha), beta, power_A, power_B &
, overlap_x2, overlap_y2, overlap_z2, overlap2, dim1 )
overlap_x = 2.d0 * real( overlap_x1 + overlap_x2 )
overlap_y = 2.d0 * real( overlap_y1 + overlap_y2 )
overlap_z = 2.d0 * real( overlap_z1 + overlap_z2 )
overlap = 2.d0 * real( overlap1 + overlap2 )
ao_overlap_cosgtos(i,j) = ao_overlap_cosgtos(i,j) + c * overlap
if( isnan(ao_overlap_cosgtos(i,j)) ) then
print*,'i, j', i, j
print*,'l, n', l, n
print*,'c, overlap', c, overlap
print*, overlap_x, overlap_y, overlap_z
stop
endif
ao_overlap_cosgtos_x(i,j) = ao_overlap_cosgtos_x(i,j) + c * overlap_x
ao_overlap_cosgtos_y(i,j) = ao_overlap_cosgtos_y(i,j) + c * overlap_y
ao_overlap_cosgtos_z(i,j) = ao_overlap_cosgtos_z(i,j) + c * overlap_z
enddo
enddo
enddo
enddo
!$OMP END PARALLEL DO
END_PROVIDER
! ---

7
src/cosgtos_ao_int/cosgtos_ao_int.irp.f Normal file
View File

@ -0,0 +1,7 @@
program cosgtos_ao_int
implicit none
BEGIN_DOC
! TODO : Put the documentation of the program here
END_DOC
print *, 'Hello world'
end

0
src/ao_two_e_ints/gauss_legendre.irp.f → src/cosgtos_ao_int/gauss_legendre.irp.f
View File

535
src/cosgtos_ao_int/one_e_Coul_integrals.irp.f Normal file
View File

@ -0,0 +1,535 @@
! ---
BEGIN_PROVIDER [ double precision, ao_integrals_n_e_cosgtos, (ao_num, ao_num)]
BEGIN_DOC
!
! Nucleus-electron interaction, in the cosgtos |AO| basis set.
!
! :math:`\langle \chi_i | -\sum_A \frac{1}{|r-R_A|} | \chi_j \rangle`
!
END_DOC
implicit none
integer :: num_A, num_B, power_A(3), power_B(3)
integer :: i, j, k, l, n_pt_in, m
double precision :: c, Z, A_center(3), B_center(3), C_center(3)
complex*16 :: alpha, beta, c1, c2
complex*16 :: NAI_pol_mult_cosgtos
ao_integrals_n_e_cosgtos = 0.d0
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE ( i, j, k, l, m, alpha, beta, A_center, B_center, C_center &
!$OMP , power_A, power_B, num_A, num_B, Z, c, c1, c2, n_pt_in ) &
!$OMP SHARED ( ao_num, ao_prim_num, ao_nucl, nucl_coord, ao_power, nucl_num, nucl_charge &
!$OMP , ao_expo_ord_transp_cosgtos, ao_coef_norm_ord_transp_cosgtos &
!$OMP , n_pt_max_integrals, ao_integrals_n_e_cosgtos )
n_pt_in = n_pt_max_integrals
!$OMP DO SCHEDULE (dynamic)
do j = 1, ao_num
num_A = ao_nucl(j)
power_A(1:3) = ao_power(j,1:3)
A_center(1:3) = nucl_coord(num_A,1:3)
do i = 1, ao_num
num_B = ao_nucl(i)
power_B(1:3) = ao_power(i,1:3)
B_center(1:3) = nucl_coord(num_B,1:3)
do l = 1, ao_prim_num(j)
alpha = ao_expo_ord_transp_cosgtos(l,j)
do m = 1, ao_prim_num(i)
beta = ao_expo_ord_transp_cosgtos(m,i)
c = 0.d0
do k = 1, nucl_num
Z = nucl_charge(k)
C_center(1:3) = nucl_coord(k,1:3)
!print *, ' '
!print *, A_center, B_center, C_center, power_A, power_B
!print *, real(alpha), real(beta)
c1 = NAI_pol_mult_cosgtos( A_center, B_center, power_A, power_B &
, alpha, beta, C_center, n_pt_in )
!c2 = c1
c2 = NAI_pol_mult_cosgtos( A_center, B_center, power_A, power_B &
, conjg(alpha), beta, C_center, n_pt_in )
!print *, ' c1 = ', real(c1)
!print *, ' c2 = ', real(c2)
c = c - Z * 2.d0 * real(c1 + c2)
enddo
ao_integrals_n_e_cosgtos(i,j) = ao_integrals_n_e_cosgtos(i,j) &
+ ao_coef_norm_ord_transp_cosgtos(l,j) &
* ao_coef_norm_ord_transp_cosgtos(m,i) * c
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
END_PROVIDER
! ---
complex*16 function NAI_pol_mult_cosgtos(A_center, B_center, power_A, power_B, alpha, beta, C_center, n_pt_in)
BEGIN_DOC
!
! Computes the electron-nucleus attraction with two primitves cosgtos.
!
! :math:`\langle g_i | \frac{1}{|r-R_c|} | g_j \rangle`
!
END_DOC
implicit none
include 'utils/constants.include.F'
integer, intent(in) :: n_pt_in, power_A(3), power_B(3)
double precision, intent(in) :: C_center(3), A_center(3), B_center(3)
complex*16, intent(in) :: alpha, beta
integer :: i, n_pt, n_pt_out
double precision :: dist, const_mod
complex*16 :: p, p_inv, rho, dist_integral, const, const_factor, coeff, factor
complex*16 :: accu, P_center(3)
complex*16 :: d(0:n_pt_in)
complex*16 :: V_n_e_cosgtos
complex*16 :: crint
if ( (A_center(1)/=B_center(1)) .or. (A_center(2)/=B_center(2)) .or. (A_center(3)/=B_center(3)) .or. &
(A_center(1)/=C_center(1)) .or. (A_center(2)/=C_center(2)) .or. (A_center(3)/=C_center(3)) ) then
continue
else
NAI_pol_mult_cosgtos = V_n_e_cosgtos( power_A(1), power_A(2), power_A(3) &
, power_B(1), power_B(2), power_B(3) &
, alpha, beta )
return
endif
p = alpha + beta
p_inv = (1.d0, 0.d0) / p
rho = alpha * beta * p_inv
dist = 0.d0
dist_integral = (0.d0, 0.d0)
do i = 1, 3
P_center(i) = (alpha * A_center(i) + beta * B_center(i)) * p_inv
dist += (A_center(i) - B_center(i)) * (A_center(i) - B_center(i))
dist_integral += (P_center(i) - C_center(i)) * (P_center(i) - C_center(i))
enddo
const_factor = dist * rho
const = p * dist_integral
const_mod = dsqrt(real(const_factor)*real(const_factor) + aimag(const_factor)*aimag(const_factor))
if(const_mod > 80.d0) then
NAI_pol_mult_cosgtos = (0.d0, 0.d0)
return
endif
factor = zexp(-const_factor)
coeff = dtwo_pi * factor * p_inv
do i = 0, n_pt_in
d(i) = (0.d0, 0.d0)
enddo
n_pt = 2 * ( (power_A(1) + power_B(1)) + (power_A(2) + power_B(2)) + (power_A(3) + power_B(3)) )
if(n_pt == 0) then
NAI_pol_mult_cosgtos = coeff * crint(0, const)
return
endif
call give_cpolynomial_mult_center_one_e( A_center, B_center, alpha, beta &
, power_A, power_B, C_center, n_pt_in, d, n_pt_out)
if(n_pt_out < 0) then
NAI_pol_mult_cosgtos = (0.d0, 0.d0)
return
endif
accu = (0.d0, 0.d0)
do i = 0, n_pt_out, 2
accu += crint(shiftr(i, 1), const) * d(i)
! print *, shiftr(i, 1), real(const), real(d(i)), real(crint(shiftr(i, 1), const))
enddo
NAI_pol_mult_cosgtos = accu * coeff
end function NAI_pol_mult_cosgtos
! ---
subroutine give_cpolynomial_mult_center_one_e( A_center, B_center, alpha, beta &
, power_A, power_B, C_center, n_pt_in, d, n_pt_out)
BEGIN_DOC
! Returns the explicit polynomial in terms of the "t" variable of the following
!
! $I_{x1}(a_x, d_x,p,q) \times I_{x1}(a_y, d_y,p,q) \times I_{x1}(a_z, d_z,p,q)$.
END_DOC
implicit none
integer, intent(in) :: n_pt_in, power_A(3), power_B(3)
double precision, intent(in) :: A_center(3), B_center(3), C_center(3)
complex*16, intent(in) :: alpha, beta
integer, intent(out) :: n_pt_out
complex*16, intent(out) :: d(0:n_pt_in)
integer :: a_x, b_x, a_y, b_y, a_z, b_z
integer :: n_pt1, n_pt2, n_pt3, dim, i, n_pt_tmp
complex*16 :: p, P_center(3), rho, p_inv, p_inv_2
complex*16 :: R1x(0:2), B01(0:2), R1xp(0:2),R2x(0:2)
complex*16 :: d1(0:n_pt_in), d2(0:n_pt_in), d3(0:n_pt_in)
ASSERT (n_pt_in > 1)
p = alpha + beta
p_inv = (1.d0, 0.d0) / p
p_inv_2 = 0.5d0 * p_inv
do i = 1, 3
P_center(i) = (alpha * A_center(i) + beta * B_center(i)) * p_inv
enddo
do i = 0, n_pt_in
d(i) = (0.d0, 0.d0)
d1(i) = (0.d0, 0.d0)
d2(i) = (0.d0, 0.d0)
d3(i) = (0.d0, 0.d0)
enddo
! ---
n_pt1 = n_pt_in
R1x(0) = (P_center(1) - A_center(1))
R1x(1) = (0.d0, 0.d0)
R1x(2) = -(P_center(1) - C_center(1))
R1xp(0) = (P_center(1) - B_center(1))
R1xp(1) = (0.d0, 0.d0)
R1xp(2) = -(P_center(1) - C_center(1))
R2x(0) = p_inv_2
R2x(1) = (0.d0, 0.d0)
R2x(2) = -p_inv_2
a_x = power_A(1)
b_x = power_B(1)
call I_x1_pol_mult_one_e_cosgtos(a_x, b_x, R1x, R1xp, R2x, d1, n_pt1, n_pt_in)
if(n_pt1 < 0) then
n_pt_out = -1
do i = 0, n_pt_in
d(i) = (0.d0, 0.d0)
enddo
return
endif
! ---
n_pt2 = n_pt_in
R1x(0) = (P_center(2) - A_center(2))
R1x(1) = (0.d0, 0.d0)
R1x(2) = -(P_center(2) - C_center(2))
R1xp(0) = (P_center(2) - B_center(2))
R1xp(1) = (0.d0, 0.d0)
R1xp(2) = -(P_center(2) - C_center(2))
a_y = power_A(2)
b_y = power_B(2)
call I_x1_pol_mult_one_e_cosgtos(a_y, b_y, R1x, R1xp, R2x, d2, n_pt2, n_pt_in)
if(n_pt2 < 0) then
n_pt_out = -1
do i = 0, n_pt_in
d(i) = (0.d0, 0.d0)
enddo
return
endif
! ---
n_pt3 = n_pt_in
R1x(0) = (P_center(3) - A_center(3))
R1x(1) = (0.d0, 0.d0)
R1x(2) = -(P_center(3) - C_center(3))
R1xp(0) = (P_center(3) - B_center(3))
R1xp(1) = (0.d0, 0.d0)
R1xp(2) = -(P_center(3) - C_center(3))
a_z = power_A(3)
b_z = power_B(3)
call I_x1_pol_mult_one_e_cosgtos(a_z, b_z, R1x, R1xp, R2x, d3, n_pt3, n_pt_in)
if(n_pt3 < 0) then
n_pt_out = -1
do i = 0, n_pt_in
d(i) = (0.d0, 0.d0)
enddo
return
endif
! ---
n_pt_tmp = 0
call multiply_cpoly(d1, n_pt1, d2, n_pt2, d, n_pt_tmp)
do i = 0, n_pt_tmp
d1(i) = (0.d0, 0.d0)
enddo
n_pt_out = 0
call multiply_cpoly(d, n_pt_tmp, d3, n_pt3, d1, n_pt_out)
do i = 0, n_pt_out
d(i) = d1(i)
enddo
end subroutine give_cpolynomial_mult_center_one_e
! ---
recursive subroutine I_x1_pol_mult_one_e_cosgtos(a, c, R1x, R1xp, R2x, d, nd, n_pt_in)
BEGIN_DOC
! Recursive routine involved in the electron-nucleus potential
END_DOC
implicit none
include 'utils/constants.include.F'
integer, intent(in) :: a, c, n_pt_in
complex*16, intent(in) :: R1x(0:2), R1xp(0:2), R2x(0:2)
integer, intent(inout) :: nd
complex*16, intent(inout) :: d(0:n_pt_in)
integer :: nx, ix, dim, iy, ny
complex*16 :: X(0:max_dim)
complex*16 :: Y(0:max_dim)
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: X, Y
dim = n_pt_in
if( (a==0) .and. (c==0)) then
nd = 0
d(0) = (1.d0, 0.d0)
return
elseif( (c < 0) .or. (nd < 0) ) then
nd = -1
return
elseif((a == 0) .and. (c .ne. 0)) then
call I_x2_pol_mult_one_e_cosgtos(c, R1x, R1xp, R2x, d, nd, n_pt_in)
elseif(a == 1) then
nx = nd
do ix = 0, n_pt_in
X(ix) = (0.d0, 0.d0)
Y(ix) = (0.d0, 0.d0)
enddo
call I_x2_pol_mult_one_e_cosgtos(c-1, R1x, R1xp, R2x, X, nx, n_pt_in)
do ix = 0, nx
X(ix) *= dble(c)
enddo
call multiply_cpoly(X, nx, R2x, 2, d, nd)
ny = 0
call I_x2_pol_mult_one_e_cosgtos(c, R1x, R1xp, R2x, Y, ny, n_pt_in)
call multiply_cpoly(Y, ny, R1x, 2, d, nd)
else
nx = 0
do ix = 0, n_pt_in
X(ix) = (0.d0, 0.d0)
Y(ix) = (0.d0, 0.d0)
enddo
call I_x1_pol_mult_one_e_cosgtos(a-2, c, R1x, R1xp, R2x, X, nx, n_pt_in)
do ix = 0, nx
X(ix) *= dble(a-1)
enddo
call multiply_cpoly(X, nx, R2x, 2, d, nd)
nx = nd
do ix = 0, n_pt_in
X(ix) = (0.d0, 0.d0)
enddo
call I_x1_pol_mult_one_e_cosgtos(a-1, c-1, R1x, R1xp, R2x, X, nx, n_pt_in)
do ix = 0, nx
X(ix) *= dble(c)
enddo
call multiply_cpoly(X, nx, R2x, 2, d, nd)
ny = 0
call I_x1_pol_mult_one_e_cosgtos(a-1, c, R1x, R1xp, R2x, Y, ny, n_pt_in)
call multiply_cpoly(Y, ny, R1x, 2, d, nd)
endif
end subroutine I_x1_pol_mult_one_e_cosgtos
! ---
recursive subroutine I_x2_pol_mult_one_e_cosgtos(c, R1x, R1xp, R2x, d, nd, dim)
BEGIN_DOC
! Recursive routine involved in the electron-nucleus potential
END_DOC
implicit none
include 'utils/constants.include.F'
integer, intent(in) :: dim, c
complex*16, intent(in) :: R1x(0:2), R1xp(0:2), R2x(0:2)
integer, intent(inout) :: nd
complex*16, intent(out) :: d(0:max_dim)
integer :: i, nx, ix, ny
complex*16 :: X(0:max_dim), Y(0:max_dim)
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: X, Y
if(c == 0) then
nd = 0
d(0) = (1.d0, 0.d0)
return
elseif((nd < 0) .or. (c < 0)) then
nd = -1
return
else
nx = 0
do ix = 0, dim
X(ix) = (0.d0, 0.d0)
Y(ix) = (0.d0, 0.d0)
enddo
call I_x1_pol_mult_one_e_cosgtos(0, c-2, R1x, R1xp, R2x, X, nx, dim)
do ix = 0, nx
X(ix) *= dble(c-1)
enddo
call multiply_cpoly(X, nx, R2x, 2, d, nd)
ny = 0
do ix = 0, dim
Y(ix) = (0.d0, 0.d0)
enddo
call I_x1_pol_mult_one_e_cosgtos(0, c-1, R1x, R1xp, R2x, Y, ny, dim)
if(ny .ge. 0) then
call multiply_cpoly(Y, ny, R1xp, 2, d, nd)
endif
endif
end subroutine I_x2_pol_mult_one_e_cosgtos
! ---
complex*16 function V_n_e_cosgtos(a_x, a_y, a_z, b_x, b_y, b_z, alpha, beta)
BEGIN_DOC
! Primitve nuclear attraction between the two primitves centered on the same atom.
!
! $p_1 = x^{a_x} y^{a_y} z^{a_z} \exp(-\alpha r^2)$
!
! $p_2 = x^{b_x} y^{b_y} z^{b_z} \exp(-\beta r^2)$
END_DOC
implicit none
integer, intent(in) :: a_x, a_y, a_z, b_x, b_y, b_z
complex*16, intent(in) :: alpha, beta
double precision :: V_phi, V_theta
complex*16 :: V_r_cosgtos
if( (iand(a_x + b_x, 1) == 1) .or. &
(iand(a_y + b_y, 1) == 1) .or. &
(iand(a_z + b_z, 1) == 1) ) then
V_n_e_cosgtos = (0.d0, 0.d0)
else
V_n_e_cosgtos = V_r_cosgtos(a_x + b_x + a_y + b_y + a_z + b_z + 1, alpha + beta) &
* V_phi(a_x + b_x, a_y + b_y) &
* V_theta(a_z + b_z, a_x + b_x + a_y + b_y + 1)
endif
end function V_n_e_cosgtos
! ---
complex*16 function V_r_cosgtos(n, alpha)
BEGIN_DOC
! Computes the radial part of the nuclear attraction integral:
!
! $\int_{0}^{\infty} r^n \exp(-\alpha r^2) dr$
!
END_DOC
implicit none
include 'utils/constants.include.F'
integer , intent(in) :: n
complex*16, intent(in) :: alpha
double precision :: fact
if(iand(n, 1) .eq. 1) then
V_r_cosgtos = 0.5d0 * fact(shiftr(n, 1)) / (alpha**(shiftr(n, 1) + 1))
else
V_r_cosgtos = sqpi * fact(n) / fact(shiftr(n, 1)) * (0.5d0/zsqrt(alpha))**(n+1)
endif
end function V_r_cosgtos
! ---

223
src/cosgtos_ao_int/one_e_kin_integrals.irp.f Normal file
View File

@ -0,0 +1,223 @@
! ---
BEGIN_PROVIDER [ double precision, ao_deriv2_cosgtos_x, (ao_num, ao_num) ]
&BEGIN_PROVIDER [ double precision, ao_deriv2_cosgtos_y, (ao_num, ao_num) ]
&BEGIN_PROVIDER [ double precision, ao_deriv2_cosgtos_z, (ao_num, ao_num) ]
implicit none
integer :: i, j, n, l, dim1, power_A(3), power_B(3)
double precision :: c, deriv_tmp
complex*16 :: alpha, beta, A_center(3), B_center(3)
complex*16 :: overlap_x, overlap_y, overlap_z, overlap
complex*16 :: overlap_x0_1, overlap_y0_1, overlap_z0_1
complex*16 :: overlap_x0_2, overlap_y0_2, overlap_z0_2
complex*16 :: overlap_m2_1, overlap_p2_1
complex*16 :: overlap_m2_2, overlap_p2_2
complex*16 :: deriv_tmp_1, deriv_tmp_2
dim1 = 100
! -- Dummy call to provide everything
A_center(:) = (0.0d0, 0.d0)
B_center(:) = (1.0d0, 0.d0)
alpha = (1.0d0, 0.d0)
beta = (0.1d0, 0.d0)
power_A = 1
power_B = 0
call overlap_cgaussian_xyz( A_center, B_center, alpha, beta, power_A, power_B &
, overlap_x0_1, overlap_y0_1, overlap_z0_1, overlap, dim1 )
! ---
!$OMP PARALLEL DO SCHEDULE(GUIDED) &
!$OMP DEFAULT(NONE) &
!$OMP PRIVATE( A_center, B_center, power_A, power_B, alpha, beta, i, j, l, n, c &
!$OMP , deriv_tmp, deriv_tmp_1, deriv_tmp_2 &
!$OMP , overlap_x, overlap_y, overlap_z, overlap &
!$OMP , overlap_m2_1, overlap_p2_1, overlap_m2_2, overlap_p2_2 &
!$OMP , overlap_x0_1, overlap_y0_1, overlap_z0_1, overlap_x0_2, overlap_y0_2, overlap_z0_2 ) &
!$OMP SHARED( nucl_coord, ao_power, ao_prim_num, ao_num, ao_nucl, dim1 &
!$OMP , ao_coef_norm_ord_transp_cosgtos, ao_expo_ord_transp_cosgtos &
!$OMP , ao_deriv2_cosgtos_x, ao_deriv2_cosgtos_y, ao_deriv2_cosgtos_z )
do j = 1, ao_num
A_center(1) = nucl_coord(ao_nucl(j),1) * (1.d0, 0.d0)
A_center(2) = nucl_coord(ao_nucl(j),2) * (1.d0, 0.d0)
A_center(3) = nucl_coord(ao_nucl(j),3) * (1.d0, 0.d0)
power_A(1) = ao_power(j,1)
power_A(2) = ao_power(j,2)
power_A(3) = ao_power(j,3)
do i = 1, ao_num
B_center(1) = nucl_coord(ao_nucl(i),1) * (1.d0, 0.d0)
B_center(2) = nucl_coord(ao_nucl(i),2) * (1.d0, 0.d0)
B_center(3) = nucl_coord(ao_nucl(i),3) * (1.d0, 0.d0)
power_B(1) = ao_power(i,1)
power_B(2) = ao_power(i,2)
power_B(3) = ao_power(i,3)
ao_deriv2_cosgtos_x(i,j) = 0.d0
ao_deriv2_cosgtos_y(i,j) = 0.d0
ao_deriv2_cosgtos_z(i,j) = 0.d0
do n = 1, ao_prim_num(j)
alpha = ao_expo_ord_transp_cosgtos(n,j)
do l = 1, ao_prim_num(i)
c = ao_coef_norm_ord_transp_cosgtos(n,j) * ao_coef_norm_ord_transp_cosgtos(l,i)
beta = ao_expo_ord_transp_cosgtos(l,i)
call overlap_cgaussian_xyz( A_center, B_center, alpha, beta, power_A, power_B &
, overlap_x0_1, overlap_y0_1, overlap_z0_1, overlap, dim1 )
call overlap_cgaussian_xyz( A_center, B_center, alpha, conjg(beta), power_A, power_B &
, overlap_x0_2, overlap_y0_2, overlap_z0_2, overlap, dim1 )
! ---
power_A(1) = power_A(1) - 2
if(power_A(1) > -1) then
call overlap_cgaussian_xyz( A_center, B_center, alpha, beta, power_A, power_B &
, overlap_m2_1, overlap_y, overlap_z, overlap, dim1 )
call overlap_cgaussian_xyz( A_center, B_center, alpha, conjg(beta), power_A, power_B &
, overlap_m2_2, overlap_y, overlap_z, overlap, dim1 )
else
overlap_m2_1 = (0.d0, 0.d0)
overlap_m2_2 = (0.d0, 0.d0)
endif
power_A(1) = power_A(1) + 4
call overlap_cgaussian_xyz( A_center, B_center, alpha, beta, power_A, power_B &
, overlap_p2_1, overlap_y, overlap_z, overlap, dim1 )
call overlap_cgaussian_xyz( A_center, B_center, alpha, conjg(beta), power_A, power_B &
, overlap_p2_2, overlap_y, overlap_z, overlap, dim1 )
power_A(1) = power_A(1) - 2
deriv_tmp_1 = ( -2.d0 * alpha * (2.d0 * power_A(1) + 1.d0) * overlap_x0_1 &
+ power_A(1) * (power_A(1) - 1.d0) * overlap_m2_1 &
+ 4.d0 * alpha * alpha * overlap_p2_1 ) * overlap_y0_1 * overlap_z0_1
deriv_tmp_2 = ( -2.d0 * alpha * (2.d0 * power_A(1) + 1.d0) * overlap_x0_2 &
+ power_A(1) * (power_A(1) - 1.d0) * overlap_m2_2 &
+ 4.d0 * alpha * alpha * overlap_p2_2 ) * overlap_y0_2 * overlap_z0_2
deriv_tmp = 2.d0 * real(deriv_tmp_1 + deriv_tmp_2)
ao_deriv2_cosgtos_x(i,j) += c * deriv_tmp
! ---
power_A(2) = power_A(2) - 2
if(power_A(2) > -1) then
call overlap_cgaussian_xyz( A_center, B_center, alpha, beta, power_A, power_B &
, overlap_x, overlap_m2_1, overlap_y, overlap, dim1 )
call overlap_cgaussian_xyz( A_center, B_center, alpha, conjg(beta), power_A, power_B &
, overlap_x, overlap_m2_2, overlap_y, overlap, dim1 )
else
overlap_m2_1 = (0.d0, 0.d0)
overlap_m2_2 = (0.d0, 0.d0)
endif
power_A(2) = power_A(2) + 4
call overlap_cgaussian_xyz( A_center, B_center, alpha, beta, power_A, power_B &
, overlap_x, overlap_p2_1, overlap_y, overlap, dim1 )
call overlap_cgaussian_xyz( A_center, B_center, alpha, conjg(beta), power_A, power_B &
, overlap_x, overlap_p2_2, overlap_y, overlap, dim1 )
power_A(2) = power_A(2) - 2
deriv_tmp_1 = ( -2.d0 * alpha * (2.d0 * power_A(2) + 1.d0) * overlap_y0_1 &
+ power_A(2) * (power_A(2) - 1.d0) * overlap_m2_1 &
+ 4.d0 * alpha * alpha * overlap_p2_1 ) * overlap_x0_1 * overlap_z0_1
deriv_tmp_2 = ( -2.d0 * alpha * (2.d0 * power_A(2) + 1.d0) * overlap_y0_2 &
+ power_A(2) * (power_A(2) - 1.d0) * overlap_m2_2 &
+ 4.d0 * alpha * alpha * overlap_p2_2 ) * overlap_x0_2 * overlap_z0_2
deriv_tmp = 2.d0 * real(deriv_tmp_1 + deriv_tmp_2)
ao_deriv2_cosgtos_y(i,j) += c * deriv_tmp
! ---
power_A(3) = power_A(3) - 2
if(power_A(3) > -1) then
call overlap_cgaussian_xyz( A_center, B_center, alpha, beta, power_A, power_B &
, overlap_x, overlap_y, overlap_m2_1, overlap, dim1 )
call overlap_cgaussian_xyz( A_center, B_center, alpha, conjg(beta), power_A, power_B &
, overlap_x, overlap_y, overlap_m2_2, overlap, dim1 )
else
overlap_m2_1 = (0.d0, 0.d0)
overlap_m2_2 = (0.d0, 0.d0)
endif
power_A(3) = power_A(3) + 4
call overlap_cgaussian_xyz( A_center, B_center, alpha, beta, power_A, power_B &
, overlap_x, overlap_y, overlap_p2_1, overlap, dim1 )
call overlap_cgaussian_xyz( A_center, B_center, alpha, conjg(beta), power_A, power_B &
, overlap_x, overlap_y, overlap_p2_2, overlap, dim1 )
power_A(3) = power_A(3) - 2
deriv_tmp_1 = ( -2.d0 * alpha * (2.d0 * power_A(3) + 1.d0) * overlap_z0_1 &
+ power_A(3) * (power_A(3) - 1.d0) * overlap_m2_1 &
+ 4.d0 * alpha * alpha * overlap_p2_1 ) * overlap_x0_1 * overlap_y0_1
deriv_tmp_2 = ( -2.d0 * alpha * (2.d0 * power_A(3) + 1.d0) * overlap_z0_2 &
+ power_A(3) * (power_A(3) - 1.d0) * overlap_m2_2 &
+ 4.d0 * alpha * alpha * overlap_p2_2 ) * overlap_x0_2 * overlap_y0_2
deriv_tmp = 2.d0 * real(deriv_tmp_1 + deriv_tmp_2)
ao_deriv2_cosgtos_z(i,j) += c * deriv_tmp
! ---
enddo
enddo
enddo
enddo
!$OMP END PARALLEL DO
END_PROVIDER
! ---
BEGIN_PROVIDER [double precision, ao_kinetic_integrals_cosgtos, (ao_num, ao_num)]
BEGIN_DOC
!
! Kinetic energy integrals in the cosgtos |AO| basis.
!
! $\langle \chi_i |\hat{T}| \chi_j \rangle$
!
END_DOC
implicit none
integer :: i, j
!$OMP PARALLEL DO DEFAULT(NONE) &
!$OMP PRIVATE(i, j) &
!$OMP SHARED(ao_num, ao_kinetic_integrals_cosgtos, ao_deriv2_cosgtos_x, ao_deriv2_cosgtos_y, ao_deriv2_cosgtos_z)
do j = 1, ao_num
do i = 1, ao_num
ao_kinetic_integrals_cosgtos(i,j) = -0.5d0 * ( ao_deriv2_cosgtos_x(i,j) &
+ ao_deriv2_cosgtos_y(i,j) &
+ ao_deriv2_cosgtos_z(i,j) )
enddo
enddo
!$OMP END PARALLEL DO
END_PROVIDER
! ---

1584
src/cosgtos_ao_int/two_e_Coul_integrals.irp.f Normal file
View File

File diff suppressed because it is too large Load Diff

120
src/utils/cgtos_one_e.irp.f Normal file
View File

@ -0,0 +1,120 @@
! ---
complex*16 function overlap_cgaussian_x(A_center, B_center, alpha, beta, power_A, power_B, dim)
BEGIN_DOC
!
! \int_{-infty}^{+infty} (x-A_x)^ax (x-B_x)^bx exp(-alpha (x-A_x)^2) exp(- beta(x-B_X)^2) dx
! with complex arguments
!
END_DOC
implicit none
include 'constants.include.F'
integer, intent(in) :: dim, power_A, power_B
complex*16, intent(in) :: A_center, B_center, alpha, beta
integer :: i, iorder_p
double precision :: fact_p_mod
complex*16 :: P_new(0:max_dim), P_center, fact_p, p, inv_sq_p
complex*16 :: Fc_integral
call give_explicit_cpoly_and_cgaussian_x( P_new, P_center, p, fact_p, iorder_p &
, alpha, beta, power_A, power_B, A_center, B_center, dim)
fact_p_mod = dsqrt(real(fact_p)*real(fact_p) + aimag(fact_p)*aimag(fact_p))
if(fact_p_mod .lt. 1.d-14) then
overlap_cgaussian_x = (0.d0, 0.d0)
return
endif
inv_sq_p = (1.d0, 0.d0) / zsqrt(p)
overlap_cgaussian_x = (0.d0, 0.d0)
do i = 0, iorder_p
overlap_cgaussian_x += P_new(i) * Fc_integral(i, inv_sq_p)
enddo
overlap_cgaussian_x *= fact_p
end function overlap_cgaussian_x
! ---
subroutine overlap_cgaussian_xyz( A_center, B_center, alpha, beta, power_A, power_B &
, overlap_x, overlap_y, overlap_z, overlap, dim )
BEGIN_DOC
!
! S_x = \int (x-A_x)^{a_x} exp(-\alpha(x-A_x)^2) (x-B_x)^{b_x} exp(-beta(x-B_x)^2) dx
! S = S_x S_y S_z
! for complex arguments
!
END_DOC
implicit none
include 'constants.include.F'
integer, intent(in) :: dim, power_A(3), power_B(3)
complex*16, intent(in) :: A_center(3), B_center(3), alpha, beta
complex*16, intent(out) :: overlap_x, overlap_y, overlap_z, overlap
integer :: i, nmax, iorder_p(3)
double precision :: fact_p_mod
complex*16 :: P_new(0:max_dim,3), P_center(3), fact_p, p, inv_sq_p
complex*16 :: F_integral_tab(0:max_dim)
complex*16 :: Fc_integral
call give_explicit_cpoly_and_cgaussian(P_new, P_center, p, fact_p, iorder_p, alpha, beta, power_A, power_B, A_center, B_center, dim)
fact_p_mod = dsqrt(real(fact_p)*real(fact_p) + aimag(fact_p)*aimag(fact_p))
if(fact_p_mod .lt. 1.d-14) then
overlap_x = (1.d-10, 0.d0)
overlap_y = (1.d-10, 0.d0)
overlap_z = (1.d-10, 0.d0)
overlap = (1.d-10, 0.d0)
return
endif
nmax = maxval(iorder_p)
inv_sq_p = (1.d0, 0.d0) / zsqrt(p)
do i = 0, nmax
F_integral_tab(i) = Fc_integral(i, inv_sq_p)
enddo
overlap_x = P_new(0,1) * F_integral_tab(0)
overlap_y = P_new(0,2) * F_integral_tab(0)
overlap_z = P_new(0,3) * F_integral_tab(0)
do i = 1, iorder_p(1)
overlap_x = overlap_x + P_new(i,1) * F_integral_tab(i)
enddo
call cgaussian_product_x(alpha, A_center(1), beta, B_center(1), fact_p, p, P_center(1))
overlap_x *= fact_p
do i = 1, iorder_p(2)
overlap_y = overlap_y + P_new(i,2) * F_integral_tab(i)
enddo
call cgaussian_product_x(alpha, A_center(2), beta, B_center(2), fact_p, p, P_center(2))
overlap_y *= fact_p
do i = 1, iorder_p(3)
overlap_z = overlap_z + P_new(i,3) * F_integral_tab(i)
enddo
call cgaussian_product_x(alpha, A_center(3), beta, B_center(3), fact_p, p, P_center(3))
overlap_z *= fact_p
overlap = overlap_x * overlap_y * overlap_z
end subroutine overlap_cgaussian_xyz
! ---

780
src/utils/cgtos_utils.irp.f Normal file
View File

@ -0,0 +1,780 @@
! ---
subroutine give_explicit_cpoly_and_cgaussian_x(P_new, P_center, p, fact_k, iorder, alpha, beta, a, b, A_center, B_center, dim)
BEGIN_DOC
! Transform the product of
! (x-x_A)^a (x-x_B)^b exp(-(r-A)^2 alpha) exp(-(r-B)^2 beta)
! into
! fact_k \sum_{i=0}^{iorder} (x-x_P)^i exp(-p(r-P)^2)
END_DOC
implicit none
include 'constants.include.F'
integer, intent(in) :: dim
integer, intent(in) :: a, b
complex*16, intent(in) :: alpha, beta, A_center, B_center
integer, intent(out) :: iorder
complex*16, intent(out) :: p, P_center, fact_k
complex*16, intent(out) :: P_new(0:max_dim)
integer :: n_new, i, j
double precision :: tmp_mod
complex*16 :: P_a(0:max_dim), P_b(0:max_dim)
complex*16 :: p_inv, ab, d_AB, tmp
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: P_a, P_b
P_new = (0.d0, 0.d0)
! new exponent
p = alpha + beta
! new center
p_inv = (1.d0, 0.d0) / p
ab = alpha * beta
P_center = (alpha * A_center + beta * B_center) * p_inv
! get the factor
d_AB = (A_center - B_center) * (A_center - B_center)
tmp = ab * p_inv * d_AB
tmp_mod = dsqrt(REAL(tmp)*REAL(tmp) + AIMAG(tmp)*AIMAG(tmp))
if(tmp_mod .lt. 50.d0) then
fact_k = zexp(-tmp)
else
fact_k = (0.d0, 0.d0)
endif
! Recenter the polynomials P_a and P_b on P_center
!DIR$ FORCEINLINE
call recentered_cpoly2(P_a(0), A_center, P_center, a, P_b(0), B_center, P_center, b)
n_new = 0
!DIR$ FORCEINLINE
call multiply_cpoly(P_a(0), a, P_b(0), b, P_new(0), n_new)
iorder = a + b
end subroutine give_explicit_cpoly_and_cgaussian_x
! ---
subroutine give_explicit_cpoly_and_cgaussian(P_new, P_center, p, fact_k, iorder, alpha, beta, a, b, A_center, B_center, dim)
BEGIN_DOC
! Transforms the product of
! (x-x_A)^a(1) (x-x_B)^b(1) (y-y_A)^a(2) (y-y_B)^b(2) (z-z_A)^a(3) (z-z_B)^b(3) exp(-(r-A)^2 alpha) exp(-(r-B)^2 beta)
! into
! fact_k * [ sum (l_x = 0,i_order(1)) P_new(l_x,1) * (x-P_center(1))^l_x ] exp (- p (x-P_center(1))^2 )
! * [ sum (l_y = 0,i_order(2)) P_new(l_y,2) * (y-P_center(2))^l_y ] exp (- p (y-P_center(2))^2 )
! * [ sum (l_z = 0,i_order(3)) P_new(l_z,3) * (z-P_center(3))^l_z ] exp (- p (z-P_center(3))^2 )
!
! WARNING ::: IF fact_k is too smal then:
! returns a "s" function centered in zero
! with an inifinite exponent and a zero polynom coef
END_DOC
implicit none
include 'constants.include.F'
integer, intent(in) :: dim, a(3), b(3)
complex*16, intent(in) :: alpha, beta, A_center(3), B_center(3)
integer, intent(out) :: iorder(3)
complex*16, intent(out) :: p, P_center(3), fact_k, P_new(0:max_dim,3)
integer :: n_new, i, j
double precision :: tmp_mod
complex*16 :: P_a(0:max_dim,3), P_b(0:max_dim,3)
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: P_a, P_b
iorder(1) = 0
iorder(2) = 0
iorder(3) = 0
P_new(0,1) = (0.d0, 0.d0)
P_new(0,2) = (0.d0, 0.d0)
P_new(0,3) = (0.d0, 0.d0)
!DIR$ FORCEINLINE
call cgaussian_product(alpha, A_center, beta, B_center, fact_k, p, P_center)
! IF fact_k is too smal then: returns a "s" function centered in zero
! with an inifinite exponent and a zero polynom coef
tmp_mod = dsqrt(REAL(fact_k)*REAL(fact_k) + AIMAG(fact_k)*AIMAG(fact_k))
if(tmp_mod < 1d-14) then
iorder = 0
p = (1.d+14, 0.d0)
fact_k = (0.d0 , 0.d0)
P_new(0:max_dim,1:3) = (0.d0 , 0.d0)
P_center(1:3) = (0.d0 , 0.d0)
return
endif
!DIR$ FORCEINLINE
call recentered_cpoly2(P_a(0,1), A_center(1), P_center(1), a(1), P_b(0,1), B_center(1), P_center(1), b(1))
iorder(1) = a(1) + b(1)
do i = 0, iorder(1)
P_new(i,1) = 0.d0
enddo
n_new = 0
!DIR$ FORCEINLINE
call multiply_cpoly(P_a(0,1), a(1), P_b(0,1), b(1), P_new(0,1), n_new)
!DIR$ FORCEINLINE
call recentered_cpoly2(P_a(0,2), A_center(2), P_center(2), a(2), P_b(0,2), B_center(2), P_center(2), b(2))
iorder(2) = a(2) + b(2)
do i = 0, iorder(2)
P_new(i,2) = 0.d0
enddo
n_new = 0
!DIR$ FORCEINLINE
call multiply_cpoly(P_a(0,2), a(2), P_b(0,2), b(2), P_new(0,2), n_new)
!DIR$ FORCEINLINE
call recentered_cpoly2(P_a(0,3), A_center(3), P_center(3), a(3), P_b(0,3), B_center(3), P_center(3), b(3))
iorder(3) = a(3) + b(3)
do i = 0, iorder(3)
P_new(i,3) = 0.d0
enddo
n_new = 0
!DIR$ FORCEINLINE
call multiply_cpoly(P_a(0,3), a(3), P_b(0,3), b(3), P_new(0,3), n_new)
end subroutine give_explicit_cpoly_and_cgaussian
! ---
!subroutine give_explicit_poly_and_gaussian_double(P_new,P_center,p,fact_k,iorder,alpha,beta,gama,a,b,A_center,B_center,Nucl_center,dim)
! BEGIN_DOC
! ! Transforms the product of
! ! (x-x_A)^a(1) (x-x_B)^b(1) (x-x_A)^a(2) (y-y_B)^b(2) (z-z_A)^a(3) (z-z_B)^b(3)
! ! exp(-(r-A)^2 alpha) exp(-(r-B)^2 beta) exp(-(r-Nucl_center)^2 gama
! !
! ! into
! ! fact_k * [ sum (l_x = 0,i_order(1)) P_new(l_x,1) * (x-P_center(1))^l_x ] exp (- p (x-P_center(1))^2 )
! ! * [ sum (l_y = 0,i_order(2)) P_new(l_y,2) * (y-P_center(2))^l_y ] exp (- p (y-P_center(2))^2 )
! ! * [ sum (l_z = 0,i_order(3)) P_new(l_z,3) * (z-P_center(3))^l_z ] exp (- p (z-P_center(3))^2 )
! END_DOC
! implicit none
! include 'constants.include.F'
! integer, intent(in) :: dim
! integer, intent(in) :: a(3),b(3) ! powers : (x-xa)**a_x = (x-A(1))**a(1)
! double precision, intent(in) :: alpha, beta, gama ! exponents
! double precision, intent(in) :: A_center(3) ! A center
! double precision, intent(in) :: B_center (3) ! B center
! double precision, intent(in) :: Nucl_center(3) ! B center
! double precision, intent(out) :: P_center(3) ! new center
! double precision, intent(out) :: p ! new exponent
! double precision, intent(out) :: fact_k ! constant factor
! double precision, intent(out) :: P_new(0:max_dim,3)! polynomial
! integer , intent(out) :: iorder(3) ! i_order(i) = order of the polynomials
!
! double precision :: P_center_tmp(3) ! new center
! double precision :: p_tmp ! new exponent
! double precision :: fact_k_tmp,fact_k_bis ! constant factor
! double precision :: P_new_tmp(0:max_dim,3)! polynomial
! integer :: i,j
! double precision :: binom_func
!
! ! First you transform the two primitives into a sum of primitive with the same center P_center_tmp and gaussian exponent p_tmp
! call give_explicit_cpoly_and_cgaussian(P_new_tmp,P_center_tmp,p_tmp,fact_k_tmp,iorder,alpha,beta,a,b,A_center,B_center,dim)
! ! Then you create the new gaussian from the product of the new one per the Nuclei one
! call cgaussian_product(p_tmp,P_center_tmp,gama,Nucl_center,fact_k_bis,p,P_center)
! fact_k = fact_k_bis * fact_k_tmp
!
! ! Then you build the coefficient of the new polynom
! do i = 0, iorder(1)
! P_new(i,1) = 0.d0
! do j = i,iorder(1)
! P_new(i,1) = P_new(i,1) + P_new_tmp(j,1) * binom_func(j,j-i) * (P_center(1) - P_center_tmp(1))**(j-i)
! enddo
! enddo
! do i = 0, iorder(2)
! P_new(i,2) = 0.d0
! do j = i,iorder(2)
! P_new(i,2) = P_new(i,2) + P_new_tmp(j,2) * binom_func(j,j-i) * (P_center(2) - P_center_tmp(2))**(j-i)
! enddo
! enddo
! do i = 0, iorder(3)
! P_new(i,3) = 0.d0
! do j = i,iorder(3)
! P_new(i,3) = P_new(i,3) + P_new_tmp(j,3) * binom_func(j,j-i) * (P_center(3) - P_center_tmp(3))**(j-i)
! enddo
! enddo
!
!end
! ---
subroutine cgaussian_product(a, xa, b, xb, k, p, xp)
BEGIN_DOC
! complex Gaussian product
! e^{-a (r-r_A)^2} e^{-b (r-r_B)^2} = k e^{-p (r-r_P)^2}
END_DOC
implicit none
complex*16, intent(in) :: a, b, xa(3), xb(3)
complex*16, intent(out) :: p, k, xp(3)
double precision :: tmp_mod
complex*16 :: p_inv, xab(3), ab
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: xab
ASSERT (REAL(a) > 0.)
ASSERT (REAL(b) > 0.)
! new exponent
p = a + b
xab(1) = xa(1) - xb(1)
xab(2) = xa(2) - xb(2)
xab(3) = xa(3) - xb(3)
p_inv = (1.d0, 0.d0) / p
ab = a * b * p_inv
k = ab * (xab(1)*xab(1) + xab(2)*xab(2) + xab(3)*xab(3))
tmp_mod = dsqrt(REAL(k)*REAL(k) + AIMAG(k)*AIMAG(k))
if(tmp_mod .gt. 40.d0) then
k = (0.d0, 0.d0)
xp(1:3) = (0.d0, 0.d0)
return
endif
k = zexp(-k)
xp(1) = ( a * xa(1) + b * xb(1) ) * p_inv
xp(2) = ( a * xa(2) + b * xb(2) ) * p_inv
xp(3) = ( a * xa(3) + b * xb(3) ) * p_inv
end subroutine cgaussian_product
! ---
subroutine cgaussian_product_x(a, xa, b, xb, k, p, xp)
BEGIN_DOC
! complex Gaussian product in 1D.
! e^{-a (x-x_A)^2} e^{-b (x-x_B)^2} = K e^{-p (x-x_P)^2}
END_DOC
implicit none
complex*16, intent(in) :: a, b, xa, xb
complex*16, intent(out) :: p, k, xp
double precision :: tmp_mod
complex*16 :: p_inv
complex*16 :: xab, ab
ASSERT (REAL(a) > 0.)
ASSERT (REAL(b) > 0.)
! new center
p = a + b
xab = xa - xb
p_inv = (1.d0, 0.d0) / p
ab = a * b * p_inv
k = ab * xab*xab
tmp_mod = dsqrt(REAL(k)*REAL(k) + AIMAG(k)*AIMAG(k))
if(tmp_mod > 40.d0) then
k = (0.d0, 0.d0)
xp = (0.d0, 0.d0)
return
endif
k = zexp(-k)
xp = (a*xa + b*xb) * p_inv
end subroutine cgaussian_product_x
! ---
subroutine multiply_cpoly(b, nb, c, nc, d, nd)
BEGIN_DOC
! Multiply two complex polynomials
! D(t) += B(t) * C(t)
END_DOC
implicit none
integer, intent(in) :: nb, nc
complex*16, intent(in) :: b(0:nb), c(0:nc)
complex*16, intent(inout) :: d(0:nb+nc)
integer, intent(out) :: nd
integer :: ndtmp, ib, ic
double precision :: tmp_mod
complex*16 :: tmp
if(ior(nc, nb) >= 0) then ! True if nc>=0 and nb>=0
continue
else
return
endif
ndtmp = nb + nc
do ic = 0, nc
d(ic) = d(ic) + c(ic) * b(0)
enddo
do ib = 1, nb
d(ib) = d(ib) + c(0) * b(ib)
do ic = 1, nc
d(ib+ic) = d(ib+ic) + c(ic) * b(ib)
enddo
enddo
do nd = ndtmp, 0, -1
tmp = d(nd)
tmp_mod = dsqrt(REAL(tmp)*REAL(tmp) + AIMAG(tmp)*AIMAG(tmp))
if(tmp_mod .lt. 1.d-15) cycle
exit
enddo
end subroutine multiply_cpoly
! ---
subroutine add_cpoly(b, nb, c, nc, d, nd)
BEGIN_DOC
! Add two complex polynomials
! D(t) += B(t) + C(t)
END_DOC
implicit none
complex*16, intent(in) :: b(0:nb), c(0:nc)
integer, intent(inout) :: nb, nc
integer, intent(out) :: nd
complex*16, intent(out) :: d(0:nb+nc)
integer :: ib
double precision :: tmp_mod
complex*16 :: tmp
nd = nb + nc
do ib = 0, max(nb, nc)
d(ib) = d(ib) + c(ib) + b(ib)
enddo
tmp = d(nd)
tmp_mod = dsqrt(REAL(tmp)*REAL(tmp) + AIMAG(tmp)*AIMAG(tmp))
do while( (tmp_mod .lt. 1.d-15) .and. (nd >= 0) )
nd -= 1
tmp = d(nd)
tmp_mod = dsqrt(REAL(tmp)*REAL(tmp) + AIMAG(tmp)*AIMAG(tmp))
if(nd < 0) exit
enddo
end subroutine add_cpoly
! ---
subroutine add_cpoly_multiply(b, nb, cst, d, nd)
BEGIN_DOC
! Add a complex polynomial multiplied by a complex constant
! D(t) += cst * B(t)
END_DOC
implicit none
integer, intent(in) :: nb
complex*16, intent(in) :: b(0:nb), cst
integer, intent(inout) :: nd
complex*16, intent(inout) :: d(0:max(nb, nd))
integer :: ib
double precision :: tmp_mod
complex*16 :: tmp
nd = max(nd, nb)
if(nd /= -1) then
do ib = 0, nb
d(ib) = d(ib) + cst * b(ib)
enddo
tmp = d(nd)
tmp_mod = dsqrt(REAL(tmp)*REAL(tmp) + AIMAG(tmp)*AIMAG(tmp))
do while(tmp_mod .lt. 1.d-15)
nd -= 1
if(nd < 0) exit
tmp = d(nd)
tmp_mod = dsqrt(REAL(tmp)*REAL(tmp) + AIMAG(tmp)*AIMAG(tmp))
enddo
endif
end subroutine add_cpoly_multiply
! ---
subroutine recentered_cpoly2(P_A, x_A, x_P, a, P_B, x_B, x_Q, b)
BEGIN_DOC
!
! write two complex polynomials (x-x_A)^a (x-x_B)^b
! as P_A(x-x_P) and P_B(x-x_Q)
!
END_DOC
implicit none
integer, intent(in) :: a, b
complex*16, intent(in) :: x_A, x_P, x_B, x_Q
complex*16, intent(out) :: P_A(0:a), P_B(0:b)
integer :: i, minab, maxab
complex*16 :: pows_a(-2:a+b+4), pows_b(-2:a+b+4)
double precision :: binom_func
if((a<0) .or. (b<0)) return
maxab = max(a, b)
minab = max(min(a, b), 0)
pows_a(0) = (1.d0, 0.d0)
pows_a(1) = x_P - x_A
pows_b(0) = (1.d0, 0.d0)
pows_b(1) = x_Q - x_B
do i = 2, maxab
pows_a(i) = pows_a(i-1) * pows_a(1)
pows_b(i) = pows_b(i-1) * pows_b(1)
enddo
P_A(0) = pows_a(a)
P_B(0) = pows_b(b)
do i = 1, min(minab, 20)
P_A(i) = binom_transp(a-i,a) * pows_a(a-i)
P_B(i) = binom_transp(b-i,b) * pows_b(b-i)
enddo
do i = minab+1, min(a, 20)
P_A(i) = binom_transp(a-i,a) * pows_a(a-i)
enddo
do i = minab+1, min(b, 20)
P_B(i) = binom_transp(b-i,b) * pows_b(b-i)
enddo
do i = 101, a
P_A(i) = binom_func(a,a-i) * pows_a(a-i)
enddo
do i = 101, b
P_B(i) = binom_func(b,b-i) * pows_b(b-i)
enddo
end subroutine recentered_cpoly2
! ---
complex*16 function Fc_integral(n, inv_sq_p)
BEGIN_DOC
! function that calculates the following integral
! \int_{\-infty}^{+\infty} x^n \exp(-p x^2) dx
! for complex valued p
END_DOC
implicit none
include 'constants.include.F'
integer, intent(in) :: n
complex*16, intent(in) :: inv_sq_p
! (n)!
double precision :: fact
if(n < 0) then
Fc_integral = (0.d0, 0.d0)
return
endif
! odd n
if(iand(n, 1) .ne. 0) then
Fc_integral = (0.d0, 0.d0)
return
endif
if(n == 0) then
Fc_integral = sqpi * inv_sq_p
return
endif
Fc_integral = sqpi * 0.5d0**n * inv_sq_p**dble(n+1) * fact(n) / fact(shiftr(n, 1))
end function Fc_integral
! ---
complex*16 function crint(n, rho)
implicit none
include 'constants.include.F'
integer, intent(in) :: n
complex*16, intent(in) :: rho
integer :: i, mmax
double precision :: rho_mod, rho_re, rho_im
double precision :: sq_rho_re, sq_rho_im
double precision :: n_tmp
complex*16 :: sq_rho, rho_inv, rho_exp
complex*16 :: crint_smallz, cpx_erf
rho_re = REAL (rho)
rho_im = AIMAG(rho)
rho_mod = dsqrt(rho_re*rho_re + rho_im*rho_im)
if(rho_mod < 10.d0) then
! small z
if(rho_mod .lt. 1.d-10) then
crint = 1.d0 / dble(n + n + 1)
else
crint = crint_smallz(n, rho)
endif
else
! large z
if(rho_mod .gt. 40.d0) then
n_tmp = dble(n) + 0.5d0
crint = 0.5d0 * gamma(n_tmp) / (rho**n_tmp)
else
! get \sqrt(rho)
sq_rho_re = sq_op5 * dsqrt(rho_re + rho_mod)
sq_rho_im = 0.5d0 * rho_im / sq_rho_re
sq_rho = sq_rho_re + (0.d0, 1.d0) * sq_rho_im
rho_exp = 0.5d0 * zexp(-rho)
rho_inv = (1.d0, 0.d0) / rho
crint = 0.5d0 * sqpi * cpx_erf(sq_rho_re, sq_rho_im) / sq_rho
mmax = n
if(mmax .gt. 0) then
do i = 0, mmax-1
crint = ((dble(i) + 0.5d0) * crint - rho_exp) * rho_inv
enddo
endif
! ***
endif
endif
! print *, n, real(rho), real(crint)
end function crint
! ---
complex*16 function crint_sum(n_pt_out, rho, d1)
implicit none
include 'constants.include.F'
integer, intent(in) :: n_pt_out
complex*16, intent(in) :: rho, d1(0:n_pt_out)
integer :: n, i, mmax
double precision :: rho_mod, rho_re, rho_im
double precision :: sq_rho_re, sq_rho_im
complex*16 :: sq_rho, F0
complex*16 :: rho_tmp, rho_inv, rho_exp
complex*16, allocatable :: Fm(:)
complex*16 :: crint_smallz, cpx_erf
rho_re = REAL (rho)
rho_im = AIMAG(rho)
rho_mod = dsqrt(rho_re*rho_re + rho_im*rho_im)
if(rho_mod < 10.d0) then
! small z
if(rho_mod .lt. 1.d-10) then
! print *, ' 111'
! print *, ' rho = ', rho
crint_sum = d1(0)
! print *, 0, 1
do i = 2, n_pt_out, 2
n = shiftr(i, 1)
crint_sum = crint_sum + d1(i) / dble(n+n+1)
! print *, n, 1.d0 / dble(n+n+1)
enddo
! ***
else
! print *, ' 222'
! print *, ' rho = ', real(rho)
! if(abs(aimag(rho)) .gt. 1d-15) then
! print *, ' complex rho', rho
! stop
! endif
crint_sum = d1(0) * crint_smallz(0, rho)
! print *, 0, real(d1(0)), real(crint_smallz(0, rho))
! if(abs(aimag(d1(0))) .gt. 1d-15) then
! print *, ' complex d1(0)', d1(0)
! stop
! endif
do i = 2, n_pt_out, 2
n = shiftr(i, 1)
crint_sum = crint_sum + d1(i) * crint_smallz(n, rho)
! print *, n, real(d1(i)), real(crint_smallz(n, rho))
! if(abs(aimag(d1(i))) .gt. 1d-15) then
! print *, ' complex d1(i)', i, d1(i)
! stop
! endif
enddo
! print *, 'sum = ', real(crint_sum)
! if(abs(aimag(crint_sum)) .gt. 1d-15) then
! print *, ' complex crint_sum', crint_sum
! stop
! endif
! ***
endif
else
! large z
if(rho_mod .gt. 40.d0) then
! print *, ' 333'
! print *, ' rho = ', rho
rho_inv = (1.d0, 0.d0) / rho
rho_tmp = 0.5d0 * sqpi * zsqrt(rho_inv)
crint_sum = rho_tmp * d1(0)
! print *, 0, rho_tmp
do i = 2, n_pt_out, 2
n = shiftr(i, 1)
rho_tmp = rho_tmp * (dble(n) + 0.5d0) * rho_inv
crint_sum = crint_sum + rho_tmp * d1(i)
! print *, n, rho_tmp
enddo
! ***
else
! print *, ' 444'
! print *, ' rho = ', rho
! get \sqrt(rho)
sq_rho_re = sq_op5 * dsqrt(rho_re + rho_mod)
sq_rho_im = 0.5d0 * rho_im / sq_rho_re
sq_rho = sq_rho_re + (0.d0, 1.d0) * sq_rho_im
!sq_rho = zsqrt(rho)
F0 = 0.5d0 * sqpi * cpx_erf(sq_rho_re, sq_rho_im) / sq_rho
crint_sum = F0 * d1(0)
! print *, 0, F0
rho_exp = 0.5d0 * zexp(-rho)
rho_inv = (1.d0, 0.d0) / rho
mmax = shiftr(n_pt_out, 1)
if(mmax .gt. 0) then
allocate( Fm(mmax) )
Fm(1:mmax) = (0.d0, 0.d0)
do n = 0, mmax-1
F0 = ((dble(n) + 0.5d0) * F0 - rho_exp) * rho_inv
Fm(n+1) = F0
! print *, n, F0
enddo
do i = 2, n_pt_out, 2
n = shiftr(i, 1)
crint_sum = crint_sum + Fm(n) * d1(i)
enddo
deallocate(Fm)
endif
! ***
endif
endif
end function crint_sum
! ---
complex*16 function crint_smallz(n, rho)
BEGIN_DOC
! Standard version of rint
END_DOC
implicit none
integer, intent(in) :: n
complex*16, intent(in) :: rho
integer, parameter :: kmax = 40
double precision, parameter :: eps = 1.d-13
integer :: k
double precision :: delta_mod
complex*16 :: rho_k, ct, delta_k
ct = 0.5d0 * zexp(-rho) * gamma(dble(n) + 0.5d0)
rho_k = (1.d0, 0.d0)
crint_smallz = ct * rho_k / gamma(dble(n) + 1.5d0)
do k = 1, kmax
rho_k = rho_k * rho
delta_k = ct * rho_k / gamma(dble(n+k) + 1.5d0)
crint_smallz = crint_smallz + delta_k
delta_mod = dsqrt(REAL(delta_k)*REAL(delta_k) + AIMAG(delta_k)*AIMAG(delta_k))
if(delta_mod .lt. eps) return
enddo
if(delta_mod > eps) then
write(*,*) ' pb in crint_smallz !'
write(*,*) ' n, rho = ', n, rho
write(*,*) ' delta_mod = ', delta_mod
stop 1
endif
end function crint_smallz
! ---

204
src/utils/cpx_erf.irp.f Normal file
View File

@ -0,0 +1,204 @@
! ---
complex*16 function cpx_erf(x, y)
BEGIN_DOC
!
! compute erf(z) for z = x + i y
!
! REF: Abramowitz and Stegun
!
END_DOC
implicit none
double precision, intent(in) :: x, y
double precision :: yabs
complex*16 :: erf_tmp1, erf_tmp2, erf_tmp3, erf_tot
double precision :: erf_F
complex*16 :: erf_E, erf_G, erf_H
yabs = dabs(y)
if(yabs .lt. 1.d-15) then
cpx_erf = (1.d0, 0.d0) * derf(x)
return
else
erf_tmp1 = (1.d0, 0.d0) * derf(x)
erf_tmp2 = erf_E(x, yabs) + erf_F(x, yabs)
erf_tmp3 = zexp(-(0.d0, 2.d0) * x * yabs) * ( erf_G(x, yabs) + erf_H(x, yabs) )
erf_tot = erf_tmp1 + erf_tmp2 - erf_tmp3
endif
if(y .gt. 0.d0) then
cpx_erf = erf_tot
else
cpx_erf = CONJG(erf_tot)
endif
end function cpx_erf
! ---
complex*16 function erf_E(x, yabs)
implicit none
include 'constants.include.F'
double precision, intent(in) :: x, yabs
if( (dabs(x).gt.6.d0) .or. (x==0.d0) ) then
erf_E = (0.d0, 0.d0)
return
endif
if(dabs(x) .lt. 1.d-7) then
erf_E = -inv_pi * (0.d0, 1.d0) * yabs
else
erf_E = 0.5d0 * inv_pi * dexp(-x*x) &
* ((1.d0, 0.d0) - zexp(-(2.d0, 0.d0) * x * yabs)) / x
endif
end function erf_E
! ---
double precision function erf_F(x, yabs)
implicit none
include 'constants.include.F'
double precision, intent(in) :: x, yabs
integer, parameter :: Nmax = 13
integer :: i
double precision :: tmp1, tmp2, x2, ct
if(dabs(x) .gt. 5.8d0) then
erf_F = 0.d0
else
x2 = x * x
ct = x * inv_pi
erf_F = 0.d0
do i = 1, Nmax
tmp1 = 0.25d0 * dble(i) * dble(i) + x2
tmp2 = dexp(-tmp1) / tmp1
erf_F = erf_F + tmp2
if(dabs(tmp2) .lt. 1d-15) exit
enddo
erf_F = ct * erf_F
endif
end function erf_F
! ---
complex*16 function erf_G(x, yabs)
implicit none
include 'constants.include.F'
double precision, intent(in) :: x, yabs
integer, parameter :: Nmax = 13
integer :: i, tmpi, imin, imax
double precision :: tmp0, tmp1, x2, idble
complex*16 :: tmp2
if(x .eq. 0.d0) then
erf_G = (0.d0, 0.d0)
return
endif
tmpi = int(2.d0 * yabs)
imin = max(1, tmpi-Nmax)
imax = tmpi + Nmax
x2 = x * x
erf_G = 0.d0
do i = imin, imax
idble = dble(i)
tmp0 = 0.5d0 * idble
tmp1 = tmp0 * tmp0 + x2
tmp2 = dexp( idble * yabs - tmp1 - dlog(tmp1) - dlog_2pi) * (x - (0.d0, 1.d0)*tmp0)
erf_G = erf_G + tmp2
enddo
end function erf_G
! ---
complex*16 function erf_H(x, yabs)
implicit none
include 'constants.include.F'
double precision, intent(in) :: x, yabs
integer, parameter :: Nmax = 13
integer :: i
double precision :: tmp0, tmp1, tmp_mod, x2, ct, idble
complex*16 :: tmp2
if(x .eq. 0.d0) then
erf_H = (0.d0, 0.d0)
return
endif
if( (dabs(x) .lt. 10d0) .and. (yabs .lt. 6.1d0) ) then
x2 = x * x
ct = 0.5d0 * inv_pi
erf_H = 0.d0
do i = 1, Nmax
idble = dble(i)
tmp0 = 0.5d0 * idble
tmp1 = tmp0 * tmp0 + x2
tmp2 = dexp(-tmp1-idble*yabs) * (x + (0.d0, 1.d0)*tmp0) / tmp1
erf_H = erf_H + tmp2
tmp_mod = dsqrt(REAL(tmp2)*REAL(tmp2) + AIMAG(tmp2)*AIMAG(tmp2))
if(tmp_mod .lt. 1d-15) exit
enddo
erf_H = ct * erf_H
else
erf_H = (0.d0, 0.d0)
endif
end function erf_H
! ---

135
src/utils/integration.irp.f
View File

@ -133,7 +133,7 @@ subroutine give_explicit_poly_and_gaussian_v(P_new, ldp, P_center, p, fact_k, io
BEGIN_DOC BEGIN_DOC
! Transforms the product of ! Transforms the product of
! (x-x_A)^a(1) (x-x_B)^b(1) (x-x_A)^a(2) (y-y_B)^b(2) (z-z_A)^a(3) (z-z_B)^b(3) exp(-(r-A)^2 alpha) exp(-(r-B)^2 beta) ! (x-x_A)^a(1) (x-x_B)^b(1) (y-y_A)^a(2) (y-y_B)^b(2) (z-z_A)^a(3) (z-z_B)^b(3) exp(-(r-A)^2 alpha) exp(-(r-B)^2 beta)
! into ! into
! fact_k * [ sum (l_x = 0,i_order(1)) P_new(l_x,1) * (x-P_center(1))^l_x ] exp (- p (x-P_center(1))^2 ) ! fact_k * [ sum (l_x = 0,i_order(1)) P_new(l_x,1) * (x-P_center(1))^l_x ] exp (- p (x-P_center(1))^2 )
! * [ sum (l_y = 0,i_order(2)) P_new(l_y,2) * (y-P_center(2))^l_y ] exp (- p (y-P_center(2))^2 ) ! * [ sum (l_y = 0,i_order(2)) P_new(l_y,2) * (y-P_center(2))^l_y ] exp (- p (y-P_center(2))^2 )
@ -427,6 +427,46 @@ subroutine gaussian_product_x(a,xa,b,xb,k,p,xp)
end subroutine end subroutine
!-
subroutine gaussian_product_x_v(a,xa,b,xb,k,p,xp,n_points)
implicit none
BEGIN_DOC
! Gaussian product in 1D with multiple xa
! e^{-a (x-x_A)^2} e^{-b (x-x_B)^2} = K_{ab}^x e^{-p (x-x_P)^2}
END_DOC
integer, intent(in) :: n_points
double precision , intent(in) :: a,b ! Exponents
double precision , intent(in) :: xa(n_points),xb ! Centers
double precision , intent(out) :: p(n_points) ! New exponent
double precision , intent(out) :: xp(n_points) ! New center
double precision , intent(out) :: k(n_points) ! Constant
double precision :: p_inv
integer :: ipoint
ASSERT (a>0.)
ASSERT (b>0.)
double precision :: xab, ab
p = a+b
p_inv = 1.d0/(a+b)
ab = a*b*p_inv
do ipoint = 1, n_points
xab = xa(ipoint)-xb
k(ipoint) = ab*xab*xab
if (k(ipoint) > 40.d0) then
k(ipoint)=0.d0
cycle
endif
k(ipoint) = exp(-k(ipoint))
xp(ipoint) = (a*xa(ipoint)+b*xb)*p_inv
enddo
end subroutine
@ -506,8 +546,10 @@ subroutine multiply_poly_v(b,nb,c,nc,d,nd,n_points)
enddo enddo
enddo enddo
enddo enddo
end end
subroutine add_poly(b,nb,c,nc,d,nd) subroutine add_poly(b,nb,c,nc,d,nd)
implicit none implicit none
BEGIN_DOC BEGIN_DOC
@ -1041,3 +1083,94 @@ double precision function rint1(n,rho)
write(*,*)'pb in rint1 k too large!' write(*,*)'pb in rint1 k too large!'
stop 1 stop 1
end end
! ---
double precision function V_phi(n, m)
BEGIN_DOC
! Computes the angular $\phi$ part of the nuclear attraction integral:
!
! $\int_{0}^{2 \pi} \cos(\phi)^n \sin(\phi)^m d\phi$.
END_DOC
implicit none
integer, intent(in) :: n, m
integer :: i
double precision :: prod
double precision :: Wallis
prod = 1.d0
do i = 0, shiftr(n, 1)-1
prod = prod/ (1.d0 + dfloat(m+1)/dfloat(n-i-i-1))
enddo
V_phi = 4.d0 * prod * Wallis(m)
end function V_phi
! ---
double precision function V_theta(n, m)
BEGIN_DOC
! Computes the angular $\theta$ part of the nuclear attraction integral:
!
! $\int_{0}^{\pi} \cos(\theta)^n \sin(\theta)^m d\theta$
END_DOC
implicit none
include 'utils/constants.include.F'
integer, intent(in) :: n, m
integer :: i
double precision :: prod
double precision :: Wallis
V_theta = 0.d0
prod = 1.d0
do i = 0, shiftr(n, 1)-1
prod = prod / (1.d0 + dfloat(m+1)/dfloat(n-i-i-1))
enddo
V_theta = (prod + prod) * Wallis(m)
end function V_theta
! ---
double precision function Wallis(n)
BEGIN_DOC
! Wallis integral:
!
! $\int_{0}^{\pi} \cos(\theta)^n d\theta$.
END_DOC
implicit none
include 'utils/constants.include.F'
integer, intent(in) :: n
integer :: p
double precision :: fact
if(iand(n, 1) .eq. 0) then
Wallis = fact(shiftr(n, 1))
Wallis = pi * fact(n) / (dble(ibset(0_8, n)) * (Wallis + Wallis) * Wallis)
else
p = shiftr(n, 1)
Wallis = fact(p)
Wallis = dble(ibset(0_8, p+p)) * Wallis * Wallis / fact(p+p+1)
endif
end function Wallis
! ---

72
src/utils/one_e_integration.irp.f
View File

@ -32,9 +32,8 @@ double precision function overlap_gaussian_x(A_center,B_center,alpha,beta,power_
end end
subroutine overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,& subroutine overlap_gaussian_xyz(A_center, B_center, alpha, beta, power_A, power_B, overlap_x, overlap_y, overlap_z, overlap, dim)
power_B,overlap_x,overlap_y,overlap_z,overlap,dim)
implicit none
BEGIN_DOC BEGIN_DOC
!.. math:: !.. math::
! !
@ -42,7 +41,10 @@ subroutine overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,&
! S = S_x S_y S_z ! S = S_x S_y S_z
! !
END_DOC END_DOC
include 'constants.include.F' include 'constants.include.F'
implicit none
integer,intent(in) :: dim ! dimension maximum for the arrays representing the polynomials integer,intent(in) :: dim ! dimension maximum for the arrays representing the polynomials
double precision,intent(in) :: A_center(3),B_center(3) ! center of the x1 functions double precision,intent(in) :: A_center(3),B_center(3) ! center of the x1 functions
double precision, intent(in) :: alpha,beta double precision, intent(in) :: alpha,beta
@ -51,17 +53,18 @@ subroutine overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,&
double precision :: P_new(0:max_dim,3),P_center(3),fact_p,p double precision :: P_new(0:max_dim,3),P_center(3),fact_p,p
double precision :: F_integral_tab(0:max_dim) double precision :: F_integral_tab(0:max_dim)
integer :: iorder_p(3) integer :: iorder_p(3)
call give_explicit_poly_and_gaussian(P_new,P_center,p,fact_p,iorder_p,alpha,beta,power_A,power_B,A_center,B_center,dim)
if(fact_p.lt.1d-20)then
overlap_x = 1.d-10
overlap_y = 1.d-10
overlap_z = 1.d-10
overlap = 1.d-10
return
endif
integer :: nmax integer :: nmax
double precision :: F_integral double precision :: F_integral
call give_explicit_poly_and_gaussian(P_new, P_center, p, fact_p, iorder_p, alpha, beta, power_A, power_B, A_center, B_center, dim)
if(fact_p.lt.1d-20)then
overlap_x = 1.d-10
overlap_y = 1.d-10
overlap_z = 1.d-10
overlap = 1.d-10
return
endif
nmax = maxval(iorder_p) nmax = maxval(iorder_p)
do i = 0,nmax do i = 0,nmax
F_integral_tab(i) = F_integral(i,p) F_integral_tab(i) = F_integral(i,p)
@ -93,40 +96,47 @@ subroutine overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,&
end end
! ---
subroutine overlap_x_abs(A_center, B_center, alpha, beta, power_A, power_B, overlap_x, lower_exp_val, dx, nx)
subroutine overlap_x_abs(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,lower_exp_val,dx,nx)
implicit none
BEGIN_DOC BEGIN_DOC
! .. math :: ! .. math ::
! !
! \int_{-infty}^{+infty} (x-A_center)^(power_A) * (x-B_center)^power_B * exp(-alpha(x-A_center)^2) * exp(-beta(x-B_center)^2) dx ! \int_{-infty}^{+infty} (x-A_center)^(power_A) * (x-B_center)^power_B * exp(-alpha(x-A_center)^2) * exp(-beta(x-B_center)^2) dx
! !
END_DOC END_DOC
integer :: i,j,k,l
integer,intent(in) :: power_A,power_B implicit none
double precision, intent(in) :: lower_exp_val
double precision,intent(in) :: A_center, B_center,alpha,beta integer, intent(in) :: power_A, power_B, nx
double precision, intent(out) :: overlap_x,dx double precision, intent(in) :: lower_exp_val, A_center, B_center, alpha, beta
integer, intent(in) :: nx double precision, intent(out) :: overlap_x, dx
double precision :: x_min,x_max,domain,x,factor,dist,p,p_inv,rho
double precision :: P_center integer :: i, j, k, l
if(power_A.lt.0.or.power_B.lt.0)then double precision :: x_min, x_max, domain, x, factor, dist, p, p_inv, rho
double precision :: P_center
double precision :: tmp
if(power_A.lt.0 .or. power_B.lt.0) then
overlap_x = 0.d0 overlap_x = 0.d0
dx = 0.d0 dx = 0.d0
return return
endif endif
p = alpha + beta
p_inv= 1.d0/p p = alpha + beta
rho = alpha * beta * p_inv p_inv = 1.d0/p
dist = (A_center - B_center)*(A_center - B_center) rho = alpha * beta * p_inv
dist = (A_center - B_center)*(A_center - B_center)
P_center = (alpha * A_center + beta * B_center) * p_inv P_center = (alpha * A_center + beta * B_center) * p_inv
if(rho*dist.gt.80.d0)then
if(rho*dist.gt.80.d0) then
overlap_x= 0.d0 overlap_x= 0.d0
return return
endif endif
factor = dexp(-rho * dist) factor = dexp(-rho * dist)
double precision :: tmp
tmp = dsqrt(lower_exp_val/p) tmp = dsqrt(lower_exp_val/p)
x_min = P_center - tmp x_min = P_center - tmp
@ -143,7 +153,7 @@ subroutine overlap_x_abs(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,
overlap_x = factor * dx * overlap_x overlap_x = factor * dx * overlap_x
end end
! ---
subroutine overlap_gaussian_xyz_v(A_center, B_center, alpha, beta, power_A, power_B, overlap, n_points) subroutine overlap_gaussian_xyz_v(A_center, B_center, alpha, beta, power_A, power_B, overlap, n_points)
@ -173,7 +183,7 @@ subroutine overlap_gaussian_xyz_v(A_center, B_center, alpha, beta, power_A, powe
double precision :: F_integral double precision :: F_integral
double precision, allocatable :: P_new(:,:,:), P_center(:,:), fact_p(:) double precision, allocatable :: P_new(:,:,:), P_center(:,:), fact_p(:)
ldp = maxval(power_A(1:3) + power_B(1:3)) ldp = maxval( power_A(1:3) + power_B(1:3) )
allocate(P_new(n_points,0:ldp,3), P_center(n_points,3), fact_p(n_points)) allocate(P_new(n_points,0:ldp,3), P_center(n_points,3), fact_p(n_points))

27
src/utils/util.irp.f
View File

@ -460,6 +460,33 @@ subroutine v2_over_x(v,x,res)
end end
! ---
subroutine check_sym(A, n)
implicit none
integer, intent(in) :: n
double precision, intent(in) :: A(n,n)
integer :: i, j
double precision :: dev_sym, norm, tmp
dev_sym = 0.d0
norm = 0.d0
do i = 1, n
do j = i+1, n
tmp = A(j,i) - A(i,j)
dev_sym += tmp * tmp
norm += A(j,i) * A(j,i)
enddo
enddo
print*, ' deviation from sym = ', dev_sym
print*, ' norm = ', norm
end subroutine check_sym
! ---
subroutine sum_A_At(A, N) subroutine sum_A_At(A, N)
!BEGIN_DOC !BEGIN_DOC