LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2025-04-25 17:54:44 +02:00
Code Releases Activity

cleaning in ao_one_e_ints

Browse Source
This commit is contained in:
eginer 2025-04-16 12:21:12 +02:00
parent b4b6ef4266
commit 25e7b33bae
21 changed files with 223 additions and 5674 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

65
src/ao_basis/spherical_to_cartesian.irp.f
View File

@ -846,6 +846,71 @@ END_PROVIDER
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_cart_to_sphe_coef, (ao_cart_num,ao_sphe_num)]
implicit none
BEGIN_DOC
! Coefficients to go from cartesian to spherical coordinates in the current
! basis set
!
! S_cart^-1 <cart|sphe>
END_DOC
integer :: i
integer, external :: ao_power_index
integer :: ibegin,j,k
integer :: prev, ao_sphe_count
prev = 0
ao_cart_to_sphe_coef(:,:) = 0.d0
if (ao_cartesian) then
! Identity matrix
do i=1,ao_sphe_num
ao_cart_to_sphe_coef(i,i) = 1.d0
enddo
else
! Assume order provided by ao_power_index
i = 1
ao_sphe_count = 0
do while (i <= ao_num)
select case ( ao_l(i) )
case (0)
ao_sphe_count += 1
ao_cart_to_sphe_coef(i,ao_sphe_count) = 1.d0
i += 1
BEGIN_TEMPLATE
case ($SHELL)
if (ao_power(i,1) == $SHELL) then
do k=1,size(cart_to_sphe_$SHELL,2)
do j=1,size(cart_to_sphe_$SHELL,1)
ao_cart_to_sphe_coef(i+j-1,ao_sphe_count+k) = cart_to_sphe_$SHELL(j,k)
enddo
enddo
i += size(cart_to_sphe_$SHELL,1)
ao_sphe_count += size(cart_to_sphe_$SHELL,2)
endif
SUBST [ SHELL ]
1;;
2;;
3;;
4;;
5;;
6;;
7;;
8;;
9;;
END_TEMPLATE
case default
stop 'Error in ao_cart_to_sphe : angular momentum too high'
end select
enddo
endif
if (ao_sphe_count /= ao_sphe_num) then
call qp_bug(irp_here, ao_sphe_count, "ao_sphe_count /= ao_sphe_num")
endif
END_PROVIDER

6
src/ao_cart_one_e_ints/one_e_coul_integrals_cgtos.irp.f
View File

@ -26,7 +26,7 @@ BEGIN_PROVIDER [double precision, ao_cart_integrals_n_e_cgtos, (ao_cart_num, ao_
ao_cart_integrals_n_e_cgtos = 0.d0
ao_cart_coul_n_e_cgtos = 0.d0
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
@ -37,7 +37,7 @@ BEGIN_PROVIDER [double precision, ao_cart_integrals_n_e_cgtos, (ao_cart_num, ao_
!$OMP ao_cart_power, nucl_num, nucl_charge, n_pt_max_integrals, &
!$OMP ao_cart_expo_cgtos_ord_transp, ao_cart_coef_cgtos_norm_ord_transp, &
!$OMP ao_cart_expo_pw_ord_transp, ao_cart_expo_phase_ord_transp, &
!$OMP ao_cart_integrals_n_e_cgtos)
!$OMP ao_cart_coul_n_e_cgtos)
!$OMP DO SCHEDULE (dynamic)
do j = 1, ao_cart_num
@ -91,7 +91,7 @@ BEGIN_PROVIDER [double precision, ao_cart_integrals_n_e_cgtos, (ao_cart_num, ao_
c = c - Z * 2.d0 * real(C1 * I1 + C2 * I2)
enddo
ao_cart_integrals_n_e_cgtos(i,j) += c * ao_cart_coef_cgtos_norm_ord_transp(n,j) &
ao_cart_coul_n_e_cgtos(i,j) += c * ao_cart_coef_cgtos_norm_ord_transp(n,j) &
* ao_cart_coef_cgtos_norm_ord_transp(l,i)
enddo
enddo

8
src/ao_cart_one_e_ints/pot_ao_erf_ints.irp.f
View File

@ -1,7 +1,7 @@
! ---
subroutine give_all_erf_kl_ao(integrals_ao,mu_in,C_center)
subroutine give_all_erf_kl_ao_cart(integrals_ao,mu_in,C_center)
implicit none
BEGIN_DOC
! Subroutine that returns all integrals over $r$ of type
@ -9,18 +9,18 @@ subroutine give_all_erf_kl_ao(integrals_ao,mu_in,C_center)
END_DOC
double precision, intent(in) :: mu_in,C_center(3)
double precision, intent(out) :: integrals_ao(ao_cart_num,ao_cart_num)
double precision :: NAI_pol_mult_erf_ao
double precision :: NAI_pol_mult_erf_ao_cart
integer :: i,j,l,k,m
do k = 1, ao_cart_num
do m = 1, ao_cart_num
integrals_ao(m,k) = NAI_pol_mult_erf_ao(m,k,mu_in,C_center)
integrals_ao(m,k) = NAI_pol_mult_erf_ao_cart(m,k,mu_in,C_center)
enddo
enddo
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_cart(i_ao, j_ao, mu_in, C_center)
BEGIN_DOC
!

40
src/ao_cart_one_e_ints/pot_ao_ints.irp.f
View File

@ -1,7 +1,7 @@
! ---
BEGIN_PROVIDER [ double precision, ao_cart_integrals_n_e, (ao_cart_num,ao_cart_num)]
BEGIN_PROVIDER [ double precision, ao_cart_coul_n_e, (ao_cart_num,ao_cart_num)]
BEGIN_DOC
! Nucleus-electron interaction, in the |AO| basis set.
@ -18,11 +18,11 @@ BEGIN_PROVIDER [ double precision, ao_cart_integrals_n_e, (ao_cart_num,ao_cart_n
double precision :: A_center(3),B_center(3),C_center(3)
double precision :: overlap_x,overlap_y,overlap_z,overlap,dx,NAI_pol_mult
ao_cart_integrals_n_e = 0.d0
ao_cart_coul_n_e = 0.d0
if (read_ao_cart_integrals_n_e) then
if (read_ao_cart_coul_n_e) then
call ezfio_get_ao_cart_one_e_ints_ao_cart_integrals_n_e(ao_cart_integrals_n_e)
call ezfio_get_ao_cart_one_e_ints_ao_cart_coul_n_e(ao_cart_coul_n_e)
print *, 'AO N-e integrals read from disk'
else
@ -31,7 +31,7 @@ BEGIN_PROVIDER [ double precision, ao_cart_integrals_n_e, (ao_cart_num,ao_cart_n
do j = 1, ao_cart_num
do i = 1, ao_cart_num
ao_cart_integrals_n_e(i,j) = ao_cart_integrals_n_e_cgtos(i,j)
ao_cart_coul_n_e(i,j) = ao_cart_coul_n_e_cgtos(i,j)
enddo
enddo
@ -42,7 +42,7 @@ BEGIN_PROVIDER [ double precision, ao_cart_integrals_n_e, (ao_cart_num,ao_cart_n
!$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_cart_num,ao_cart_prim_num,ao_cart_expo_ordered_transp,ao_cart_power,ao_cart_nucl,nucl_coord,ao_cart_coef_normalized_ordered_transp,&
!$OMP n_pt_max_integrals,ao_cart_integrals_n_e,nucl_num,nucl_charge)
!$OMP n_pt_max_integrals,ao_cart_coul_n_e,nucl_num,nucl_charge)
n_pt_in = n_pt_max_integrals
@ -86,7 +86,7 @@ BEGIN_PROVIDER [ double precision, ao_cart_integrals_n_e, (ao_cart_num,ao_cart_n
c = c - Z * c1
enddo
ao_cart_integrals_n_e(i,j) = ao_cart_integrals_n_e(i,j) &
ao_cart_coul_n_e(i,j) = ao_cart_coul_n_e(i,j) &
+ ao_cart_coef_normalized_ordered_transp(l,j) &
* ao_cart_coef_normalized_ordered_transp(m,i) * c
enddo
@ -99,25 +99,17 @@ BEGIN_PROVIDER [ double precision, ao_cart_integrals_n_e, (ao_cart_num,ao_cart_n
endif
IF(do_pseudo) THEN
ao_cart_integrals_n_e += ao_cart_pseudo_integrals
ENDIF
IF(point_charges) THEN
ao_cart_integrals_n_e += ao_cart_integrals_pt_chrg
ENDIF
endif
if (write_ao_cart_integrals_n_e) then
call ezfio_set_ao_cart_one_e_ints_ao_cart_integrals_n_e(ao_cart_integrals_n_e)
if (write_ao_cart_coul_n_e) then
call ezfio_set_ao_cart_one_e_ints_ao_cart_coul_n_e(ao_cart_coul_n_e)
print *, 'AO N-e integrals written to disk'
endif
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_cart_integrals_n_e_imag, (ao_cart_num,ao_cart_num)]
BEGIN_PROVIDER [ double precision, ao_cart_coul_n_e_imag, (ao_cart_num,ao_cart_num)]
BEGIN_DOC
! Nucleus-electron interaction, in the |AO| basis set.
!
@ -131,8 +123,8 @@ BEGIN_PROVIDER [ double precision, ao_cart_integrals_n_e_imag, (ao_cart_num,ao_c
integer :: i,j,k,l,n_pt_in,m
double precision :: overlap_x,overlap_y,overlap_z,overlap,dx,NAI_pol_mult
if (read_ao_cart_integrals_n_e) then
call ezfio_get_ao_cart_one_e_ints_ao_cart_integrals_n_e_imag(ao_cart_integrals_n_e_imag)
if (read_ao_cart_coul_n_e) then
call ezfio_get_ao_cart_one_e_ints_ao_cart_coul_n_e_imag(ao_cart_coul_n_e_imag)
print *, 'AO N-e integrals read from disk'
else
print *, irp_here, ': Not yet implemented'
@ -140,7 +132,7 @@ BEGIN_PROVIDER [ double precision, ao_cart_integrals_n_e_imag, (ao_cart_num,ao_c
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_cart_integrals_n_e_per_atom, (ao_cart_num,ao_cart_num,nucl_num)]
BEGIN_PROVIDER [ double precision, ao_cart_coul_n_e_per_atom, (ao_cart_num,ao_cart_num,nucl_num)]
BEGIN_DOC
! Nucleus-electron interaction in the |AO| basis set, per atom A.
!
@ -154,14 +146,14 @@ BEGIN_PROVIDER [ double precision, ao_cart_integrals_n_e_per_atom, (ao_cart_num,
integer :: i,j,k,l,n_pt_in,m
double precision :: overlap_x,overlap_y,overlap_z,overlap,dx,NAI_pol_mult
ao_cart_integrals_n_e_per_atom = 0.d0
ao_cart_coul_n_e_per_atom = 0.d0
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,l,m,alpha,beta,A_center,B_center,power_A,power_B,&
!$OMP num_A,num_B,c,n_pt_in,C_center) &
!$OMP SHARED (ao_cart_num,ao_cart_prim_num,ao_cart_expo_ordered_transp,ao_cart_power,ao_cart_nucl,nucl_coord,ao_cart_coef_normalized_ordered_transp,&
!$OMP n_pt_max_integrals,ao_cart_integrals_n_e_per_atom,nucl_num)
!$OMP n_pt_max_integrals,ao_cart_coul_n_e_per_atom,nucl_num)
n_pt_in = n_pt_max_integrals
!$OMP DO SCHEDULE (dynamic)
@ -197,7 +189,7 @@ BEGIN_PROVIDER [ double precision, ao_cart_integrals_n_e_per_atom, (ao_cart_num,
* ao_cart_coef_normalized_ordered_transp(m,i)
enddo
enddo
ao_cart_integrals_n_e_per_atom(i,j,k) = -c
ao_cart_coul_n_e_per_atom(i,j,k) = -c
enddo
enddo
enddo

1
src/ao_one_e_ints/EZFIO.cfg
View File

@ -16,7 +16,6 @@ doc: Read/Write |AO| nucleus-electron attraction integrals from/to disk [ Write
interface: ezfio,provider,ocaml
default: None
[ao_integrals_kinetic]
type: double precision
doc: Kinetic energy integrals in |AO| basis set

1
src/ao_one_e_ints/NEED
View File

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

16
src/ao_one_e_ints/ao_one_e_ints.irp.f
View File

@ -45,19 +45,3 @@ BEGIN_PROVIDER [ double precision, ao_one_e_integrals_imag,(ao_num,ao_num)]
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_sphe_one_e_integrals,(ao_sphe_num,ao_sphe_num)]
&BEGIN_PROVIDER [ double precision, ao_sphe_one_e_integrals_diag,(ao_sphe_num)]
implicit none
integer :: i,j,n,l
BEGIN_DOC
! One-electron Hamiltonian in the spherical |AO| basis.
END_DOC
ao_sphe_one_e_integrals = ao_sphe_integrals_n_e + ao_sphe_kinetic_integrals
do j = 1, ao_num
ao_sphe_one_e_integrals_diag(j) = ao_sphe_one_e_integrals(j,j)
enddo
END_PROVIDER

132
src/ao_one_e_ints/ao_ortho_canonical.irp.f
View File

@ -1,74 +1,3 @@
BEGIN_PROVIDER [ double precision, ao_cart_to_sphe_coef, (ao_num,ao_num)]
&BEGIN_PROVIDER [ double precision, ao_cart_to_sphe_normalization, (ao_num)]
implicit none
BEGIN_DOC
! Coefficients to go from cartesian to spherical coordinates in the current
! basis set
!
! S_cart^-1 <cart|sphe>
END_DOC
integer :: i
integer, external :: ao_power_index
integer :: ibegin,j,k
integer :: prev, ao_sphe_count
prev = 0
ao_cart_to_sphe_coef(:,:) = 0.d0
ao_cart_to_sphe_normalization(:) = 1.d0
if (ao_cartesian) then
! Identity matrix
do i=1,ao_sphe_num
ao_cart_to_sphe_coef(i,i) = 1.d0
enddo
else
! Assume order provided by ao_power_index
i = 1
ao_sphe_count = 0
do while (i <= ao_num)
select case ( ao_l(i) )
case (0)
ao_sphe_count += 1
ao_cart_to_sphe_coef(i,ao_sphe_count) = 1.d0
ao_cart_to_sphe_normalization(i) = 1.d0
i += 1
BEGIN_TEMPLATE
case ($SHELL)
if (ao_power(i,1) == $SHELL) then
do k=1,size(cart_to_sphe_$SHELL,2)
do j=1,size(cart_to_sphe_$SHELL,1)
ao_cart_to_sphe_coef(i+j-1,ao_sphe_count+k) = cart_to_sphe_$SHELL(j,k)
enddo
enddo
do j=1,size(cart_to_sphe_$SHELL,1)
ao_cart_to_sphe_normalization(i+j-1) = cart_to_sphe_norm_$SHELL(j)
enddo
i += size(cart_to_sphe_$SHELL,1)
ao_sphe_count += size(cart_to_sphe_$SHELL,2)
endif
SUBST [ SHELL ]
1;;
2;;
3;;
4;;
5;;
6;;
7;;
8;;
9;;
END_TEMPLATE
case default
stop 'Error in ao_cart_to_sphe : angular momentum too high'
end select
enddo
endif
if (ao_sphe_count /= ao_sphe_num) then
call qp_bug(irp_here, ao_sphe_count, "ao_sphe_count /= ao_sphe_num")
endif
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_cart_to_sphe_overlap, (ao_sphe_num,ao_sphe_num) ]
implicit none
BEGIN_DOC
@ -146,63 +75,10 @@ END_PROVIDER
enddo
call write_double(6, lin_dep_cutoff, "Linear dependencies cut-off")
if (ao_cartesian) then
ao_ortho_canonical_num = ao_num
call ortho_canonical(ao_overlap,size(ao_overlap,1), &
ao_num,ao_ortho_canonical_coef,size(ao_ortho_canonical_coef,1), &
ao_ortho_canonical_num, lin_dep_cutoff)
else
double precision, allocatable :: S(:,:)
allocate(S(ao_sphe_num,ao_sphe_num))
S = 0.d0
do i=1,ao_sphe_num
S(i,i) = 1.d0
enddo
ao_ortho_canonical_num = ao_sphe_num
call ortho_canonical(ao_cart_to_sphe_overlap, size(ao_cart_to_sphe_overlap,1), &
ao_sphe_num, S, size(S,1), ao_ortho_canonical_num, lin_dep_cutoff)
call dgemm('N','N', ao_num, ao_ortho_canonical_num, ao_sphe_num, 1.d0, &
ao_cart_to_sphe_coef, size(ao_cart_to_sphe_coef,1), &
S, size(S,1), &
0.d0, ao_ortho_canonical_coef, size(ao_ortho_canonical_coef,1))
deallocate(S)
endif
FREE ao_overlap
ao_ortho_canonical_num = ao_num
call ortho_canonical(ao_overlap,size(ao_overlap,1), &
ao_num,ao_ortho_canonical_coef,size(ao_ortho_canonical_coef,1), &
ao_ortho_canonical_num, lin_dep_cutoff)
END_PROVIDER
BEGIN_PROVIDER [double precision, ao_ortho_canonical_overlap, (ao_ortho_canonical_num,ao_ortho_canonical_num)]
implicit none
BEGIN_DOC
! overlap matrix of the ao_ortho_canonical.
! Expected to be the Identity
END_DOC
integer :: i,j,k,l
double precision :: c
do j=1, ao_ortho_canonical_num
do i=1, ao_ortho_canonical_num
ao_ortho_canonical_overlap(i,j) = 0.d0
enddo
enddo
do j=1, ao_ortho_canonical_num
do k=1, ao_num
c = 0.d0
do l=1, ao_num
c += ao_ortho_canonical_coef(l,j) * ao_overlap(l,k)
enddo
do i=1, ao_ortho_canonical_num
ao_ortho_canonical_overlap(i,j) += ao_ortho_canonical_coef(k,i) * c
enddo
enddo
enddo
END_PROVIDER

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

@ -2,9 +2,6 @@
! ---
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_DOC
! Overlap between atomic basis functions:
@ -19,78 +16,12 @@
double precision :: A_center(3), B_center(3)
ao_overlap = 0.d0
ao_overlap_x = 0.d0
ao_overlap_y = 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))
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
if(use_cgtos) then
do j = 1, ao_num
do i = 1, ao_num
ao_overlap (i,j) = ao_overlap_cgtos (i,j)
ao_overlap_x(i,j) = ao_overlap_cgtos_x(i,j)
ao_overlap_y(i,j) = ao_overlap_cgtos_y(i,j)
ao_overlap_z(i,j) = ao_overlap_cgtos_z(i,j)
enddo
enddo
else
dim1=100
!$OMP PARALLEL DO SCHEDULE(GUIDED) &
!$OMP DEFAULT(NONE) &
!$OMP PRIVATE(A_center,B_center,power_A,power_B,&
!$OMP overlap_x,overlap_y, overlap_z, overlap, &
!$OMP alpha, beta,i,j,n,l,c) &
!$OMP SHARED(nucl_coord,ao_power,ao_prim_num, &
!$OMP ao_overlap_x,ao_overlap_y,ao_overlap_z,ao_overlap,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
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_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
!$OMP END PARALLEL DO
endif
call ao_cart_to_ao_basis(ao_cart_overlap, ao_cart_num, ao_overlap, ao_num)
endif
if (write_ao_integrals_overlap) then
@ -151,45 +82,8 @@ BEGIN_PROVIDER [ double precision, ao_overlap_abs, (ao_num, ao_num) ]
enddo
else
dim1=100
lower_exp_val = 40.d0
!$OMP PARALLEL DO SCHEDULE(GUIDED) &
!$OMP DEFAULT(NONE) &
!$OMP PRIVATE(A_center,B_center,power_A,power_B, &
!$OMP overlap_x,overlap_y, overlap_z, &
!$OMP alpha, beta,i,j,dx) &
!$OMP SHARED(nucl_coord,ao_power,ao_prim_num, &
!$OMP ao_overlap_abs,ao_num,ao_coef_normalized_ordered_transp,ao_nucl,&
!$OMP ao_expo_ordered_transp,dim1,lower_exp_val)
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_overlap_abs(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_x_abs(A_center(1),B_center(1),alpha,beta,power_A(1),power_B(1),overlap_x,lower_exp_val,dx,dim1)
call overlap_x_abs(A_center(2),B_center(2),alpha,beta,power_A(2),power_B(2),overlap_y,lower_exp_val,dx,dim1)
call overlap_x_abs(A_center(3),B_center(3),alpha,beta,power_A(3),power_B(3),overlap_z,lower_exp_val,dx,dim1)
ao_overlap_abs(i,j) += abs(ao_coef_normalized_ordered_transp(n,j) * ao_coef_normalized_ordered_transp(l,i)) * overlap_x * overlap_y * overlap_z
enddo
enddo
enddo
enddo
!$OMP END PARALLEL DO
print*,'todo ! numerical integration on DFT grid !'
stop
endif
@ -308,26 +202,3 @@ BEGIN_PROVIDER [ double precision, S_half, (ao_num,ao_num) ]
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_sphe_overlap, (ao_sphe_num,ao_sphe_num) ]
implicit none
BEGIN_DOC
! |AO| overlap matrix in the spherical basis set
END_DOC
double precision, allocatable :: tmp(:,:)
allocate (tmp(ao_sphe_num,ao_num))
call dgemm('N','N',ao_sphe_num,ao_num,ao_num, 1.d0, &
ao_cart_to_sphe_inv,size(ao_cart_to_sphe_inv,1), &
ao_overlap,size(ao_overlap,1), 0.d0, &
tmp, size(tmp,1))
call dgemm('N','T',ao_sphe_num,ao_sphe_num,ao_num, 1.d0, &
tmp, size(tmp,1), &
ao_cart_to_sphe_inv,size(ao_cart_to_sphe_inv,1), 0.d0, &
ao_sphe_overlap,size(ao_sphe_overlap,1))
deallocate(tmp)
END_PROVIDER

124
src/ao_one_e_ints/aos_cgtos.irp.f
View File

@ -1,124 +0,0 @@
! ---
BEGIN_PROVIDER [double precision, ao_overlap_cgtos, (ao_num, ao_num)]
&BEGIN_PROVIDER [double precision, ao_overlap_cgtos_x, (ao_num, ao_num)]
&BEGIN_PROVIDER [double precision, ao_overlap_cgtos_y, (ao_num, ao_num)]
&BEGIN_PROVIDER [double precision, ao_overlap_cgtos_z, (ao_num, ao_num)]
implicit none
integer :: i, j, m, n, l, ii, jj, dim1, power_A(3), power_B(3)
double precision :: c, overlap, overlap_x, overlap_y, overlap_z
double precision :: KA2(3), phiA(3)
double precision :: KB2(3), phiB(3)
complex*16 :: alpha, alpha_inv, Ae_center(3), Ap_center(3)
complex*16 :: beta, beta_inv, Be_center(3), Bp_center(3)
complex*16 :: C1(1:4), C2(1:4)
complex*16 :: overlap1, overlap_x1, overlap_y1, overlap_z1
complex*16 :: overlap2, overlap_x2, overlap_y2, overlap_z2
ao_overlap_cgtos = 0.d0
ao_overlap_cgtos_x = 0.d0
ao_overlap_cgtos_y = 0.d0
ao_overlap_cgtos_z = 0.d0
dim1 = 100
!$OMP PARALLEL DO SCHEDULE(GUIDED) &
!$OMP DEFAULT(NONE) &
!$OMP PRIVATE(i, j, m, n, l, ii, jj, c, C1, C2, &
!$OMP alpha, alpha_inv, Ae_center, Ap_center, power_A, KA2, phiA, &
!$OMP beta, beta_inv, Be_center, Bp_center, power_B, KB2, phiB, &
!$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_coef_cgtos_norm_ord_transp, ao_expo_cgtos_ord_transp, &
!$OMP ao_expo_pw_ord_transp, ao_expo_phase_ord_transp, &
!$OMP ao_overlap_cgtos_x, ao_overlap_cgtos_y, ao_overlap_cgtos_z, &
!$OMP ao_overlap_cgtos)
do j = 1, ao_num
jj = ao_nucl(j)
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
ii = ao_nucl(i)
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_cgtos_ord_transp(n,j)
alpha_inv = (1.d0, 0.d0) / alpha
do m = 1, 3
phiA(m) = ao_expo_phase_ord_transp(m,n,j)
Ap_center(m) = nucl_coord(jj,m)
Ae_center(m) = nucl_coord(jj,m) - (0.d0, 0.5d0) * alpha_inv * ao_expo_pw_ord_transp(m,n,j)
KA2(m) = ao_expo_pw_ord_transp(m,n,j) * ao_expo_pw_ord_transp(m,n,j)
enddo
do l = 1, ao_prim_num(i)
beta = ao_expo_cgtos_ord_transp(l,i)
beta_inv = (1.d0, 0.d0) / beta
do m = 1, 3
phiB(m) = ao_expo_phase_ord_transp(m,l,i)
Bp_center(m) = nucl_coord(ii,m)
Be_center(m) = nucl_coord(ii,m) - (0.d0, 0.5d0) * beta_inv * ao_expo_pw_ord_transp(m,l,i)
KB2(m) = ao_expo_pw_ord_transp(m,l,i) * ao_expo_pw_ord_transp(m,l,i)
enddo
c = ao_coef_cgtos_norm_ord_transp(n,j) * ao_coef_cgtos_norm_ord_transp(l,i)
C1(1) = zexp((0.d0, 1.d0) * (-phiA(1) - phiB(1)) - 0.25d0 * (alpha_inv * KA2(1) + beta_inv * KB2(1)))
C1(2) = zexp((0.d0, 1.d0) * (-phiA(2) - phiB(2)) - 0.25d0 * (alpha_inv * KA2(2) + beta_inv * KB2(2)))
C1(3) = zexp((0.d0, 1.d0) * (-phiA(3) - phiB(3)) - 0.25d0 * (alpha_inv * KA2(3) + beta_inv * KB2(3)))
C1(4) = C1(1) * C1(2) * C1(3)
C2(1) = zexp((0.d0, 1.d0) * (phiA(1) - phiB(1)) - 0.25d0 * (conjg(alpha_inv) * KA2(1) + beta_inv * KB2(1)))
C2(2) = zexp((0.d0, 1.d0) * (phiA(2) - phiB(2)) - 0.25d0 * (conjg(alpha_inv) * KA2(2) + beta_inv * KB2(2)))
C2(3) = zexp((0.d0, 1.d0) * (phiA(3) - phiB(3)) - 0.25d0 * (conjg(alpha_inv) * KA2(3) + beta_inv * KB2(3)))
C2(4) = C2(1) * C2(2) * C2(3)
call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, &
Ap_center, Bp_center, overlap_x1, overlap_y1, overlap_z1, overlap1, dim1)
call overlap_cgaussian_xyz(conjg(Ae_center), Be_center, conjg(alpha), beta, power_A, power_B, &
conjg(Ap_center), Bp_center, overlap_x2, overlap_y2, overlap_z2, overlap2, dim1)
overlap_x = 2.d0 * real(C1(1) * overlap_x1 + C2(1) * overlap_x2)
overlap_y = 2.d0 * real(C1(2) * overlap_y1 + C2(2) * overlap_y2)
overlap_z = 2.d0 * real(C1(3) * overlap_z1 + C2(3) * overlap_z2)
overlap = 2.d0 * real(C1(4) * overlap1 + C2(4) * overlap2 )
ao_overlap_cgtos(i,j) = ao_overlap_cgtos(i,j) + c * overlap
if(isnan(ao_overlap_cgtos(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_cgtos_x(i,j) = ao_overlap_cgtos_x(i,j) + c * overlap_x
ao_overlap_cgtos_y(i,j) = ao_overlap_cgtos_y(i,j) + c * overlap_y
ao_overlap_cgtos_z(i,j) = ao_overlap_cgtos_z(i,j) + c * overlap_z
enddo
enddo
enddo
enddo
!$OMP END PARALLEL DO
END_PROVIDER
! ---

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

@ -16,123 +16,11 @@
END_DOC
implicit none
integer :: i, j, n, l, dim1, power_A(3), power_B(3)
double precision :: overlap, overlap_y, overlap_z
double precision :: overlap_x0, overlap_y0, overlap_z0
double precision :: alpha, beta, c
double precision :: A_center(3), B_center(3)
double precision :: d_a_2,d_2
integer :: i, j
if(use_cgtos) then
do j = 1, ao_num
do i = 1, ao_num
ao_deriv2_x(i,j) = ao_deriv2_cgtos_x(i,j)
ao_deriv2_y(i,j) = ao_deriv2_cgtos_y(i,j)
ao_deriv2_z(i,j) = ao_deriv2_cgtos_z(i,j)
enddo
enddo
else
dim1=100
! -- Dummy call to provide everything
A_center(:) = 0.d0
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) &
!$OMP DEFAULT(NONE) &
!$OMP PRIVATE(A_center,B_center,power_A,power_B,&
!$OMP overlap_y, overlap_z, overlap, &
!$OMP alpha, beta, n, l, 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(1) = power_A(1)-2
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
deriv_tmp = (-2.d0 * alpha * (2.d0 * dble(power_A(1)) +1.d0) * overlap_x0 &
+dble(power_A(1)) * (dble(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
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 * dble(power_A(2)) +1.d0 ) * overlap_y0 &
+dble(power_A(2)) * (dble(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 * dble(power_A(3)) +1.d0 ) * overlap_z0 &
+dble(power_A(3)) * (dble(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
!$OMP END PARALLEL DO
endif
call ao_cart_to_ao_basis(ao_cart_deriv2_x, ao_cart_num, ao_deriv2_x, ao_num)
call ao_cart_to_ao_basis(ao_cart_deriv2_y, ao_cart_num, ao_deriv2_y, ao_num)
call ao_cart_to_ao_basis(ao_cart_deriv2_z, ao_cart_num, ao_deriv2_z, ao_num)
END_PROVIDER
@ -190,25 +78,3 @@ BEGIN_PROVIDER [double precision, ao_kinetic_integrals_imag, (ao_num,ao_num)]
endif
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_sphe_kinetic_integrals, (ao_sphe_num,ao_sphe_num) ]
implicit none
BEGIN_DOC
! |AO| kinetic inntegrals matrix in the spherical basis set
END_DOC
double precision, allocatable :: tmp(:,:)
allocate (tmp(ao_sphe_num,ao_num))
call dgemm('N','N',ao_sphe_num,ao_num,ao_num, 1.d0, &
ao_cart_to_sphe_inv,size(ao_cart_to_sphe_inv,1), &
ao_kinetic_integrals,size(ao_kinetic_integrals,1), 0.d0, &
tmp, size(tmp,1))
call dgemm('N','T',ao_sphe_num,ao_sphe_num,ao_num, 1.d0, &
tmp, size(tmp,1), &
ao_cart_to_sphe_inv,size(ao_cart_to_sphe_inv,1), 0.d0, &
ao_sphe_kinetic_integrals,size(ao_sphe_kinetic_integrals,1))
deallocate(tmp)
END_PROVIDER

558
src/ao_one_e_ints/one_e_coul_integrals_cgtos.irp.f
View File

@ -1,558 +0,0 @@
! ---
BEGIN_PROVIDER [double precision, ao_integrals_n_e_cgtos, (ao_num, ao_num)]
BEGIN_DOC
!
! Nucleus-electron interaction, in the cgtos |AO| basis set.
!
! :math:`\langle \chi_i | -\sum_A \frac{1}{|r-R_A|} | \chi_j \rangle`
!
END_DOC
implicit none
integer :: power_A(3), power_B(3)
integer :: i, j, k, l, m, n, ii, jj
double precision :: c, Z, C_center(3)
double precision :: phiA, KA2
double precision :: phiB, KB2
complex*16 :: alpha, alpha_inv, Ae_center(3), Ap_center(3)
complex*16 :: beta, beta_inv, Be_center(3), Bp_center(3)
complex*16 :: C1, C2, I1, I2
complex*16, external :: NAI_pol_mult_cgtos
ao_integrals_n_e_cgtos = 0.d0
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i, j, k, l, m, n, ii, jj, C_center, Z, c, C1, C2, I1, I2, &
!$OMP alpha, alpha_inv, Ae_center, Ap_center, phiA, KA2, power_A, &
!$OMP beta, beta_inv, Be_center, Bp_center, phiB, KB2, power_B) &
!$OMP SHARED (ao_num, ao_prim_num, ao_nucl, nucl_coord, &
!$OMP ao_power, nucl_num, nucl_charge, n_pt_max_integrals, &
!$OMP ao_expo_cgtos_ord_transp, ao_coef_cgtos_norm_ord_transp, &
!$OMP ao_expo_pw_ord_transp, ao_expo_phase_ord_transp, &
!$OMP ao_integrals_n_e_cgtos)
!$OMP DO SCHEDULE (dynamic)
do j = 1, ao_num
jj = ao_nucl(j)
power_A(1:3) = ao_power(j,1:3)
do i = 1, ao_num
ii = ao_nucl(i)
power_B(1:3) = ao_power(i,1:3)
do n = 1, ao_prim_num(j)
alpha = ao_expo_cgtos_ord_transp(n,j)
alpha_inv = (1.d0, 0.d0) / alpha
do m = 1, 3
Ap_center(m) = nucl_coord(jj,m)
Ae_center(m) = nucl_coord(jj,m) - (0.d0, 0.5d0) * alpha_inv * ao_expo_pw_ord_transp(m,n,j)
enddo
phiA = ao_expo_phase_ord_transp(4,n,j)
KA2 = ao_expo_pw_ord_transp(4,n,j)
do l = 1, ao_prim_num(i)
beta = ao_expo_cgtos_ord_transp(l,i)
beta_inv = (1.d0, 0.d0) / beta
do m = 1, 3
Bp_center(m) = nucl_coord(ii,m)
Be_center(m) = nucl_coord(ii,m) - (0.d0, 0.5d0) * beta_inv * ao_expo_pw_ord_transp(m,l,i)
enddo
phiB = ao_expo_phase_ord_transp(4,l,i)
KB2 = ao_expo_pw_ord_transp(4,l,i)
C1 = zexp((0.d0, 1.d0) * (-phiA - phiB) - 0.25d0 * (alpha_inv * KA2 + beta_inv * KB2))
C2 = zexp((0.d0, 1.d0) * ( phiA - phiB) - 0.25d0 * (conjg(alpha_inv) * KA2 + beta_inv * KB2))
c = 0.d0
do k = 1, nucl_num
Z = nucl_charge(k)
C_center(1:3) = nucl_coord(k,1:3)
I1 = NAI_pol_mult_cgtos(Ae_center, Be_center, power_A, power_B, alpha, beta, &
Ap_center, Bp_center, C_center, n_pt_max_integrals)
I2 = NAI_pol_mult_cgtos(conjg(Ae_center), Be_center, power_A, power_B, conjg(alpha), beta, &
conjg(Ap_center), Bp_center, C_center, n_pt_max_integrals)
c = c - Z * 2.d0 * real(C1 * I1 + C2 * I2)
enddo
ao_integrals_n_e_cgtos(i,j) += c * ao_coef_cgtos_norm_ord_transp(n,j) &
* ao_coef_cgtos_norm_ord_transp(l,i)
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
END_PROVIDER
! ---
complex*16 function NAI_pol_mult_cgtos(Ae_center, Be_center, power_A, power_B, alpha, beta, &
Ap_center, Bp_center, C_center, n_pt_in)
BEGIN_DOC
!
! Computes the electron-nucleus attraction with two primitves cgtos.
!
! :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)
complex*16, intent(in) :: alpha, Ae_center(3), Ap_center(3)
complex*16, intent(in) :: beta, Be_center(3), Bp_center(3)
integer :: i, n_pt, n_pt_out
double precision :: dist_AB, dist_AC
complex*16 :: p, p_inv, rho, dist, dist_integral, const, const_factor, coeff, factor
complex*16 :: P_center(3)
complex*16 :: d(0:n_pt_in)
complex*16, external :: V_n_e_cgtos
complex*16, external :: crint_sum
complex*16, external :: crint_1
dist_AB = 0.d0
dist_AC = 0.d0
do i = 1, 3
dist_AB += abs(Ae_center(i) - Be_center(i))
dist_AC += abs(Ae_center(i) - C_center(i) * (1.d0, 0.d0))
enddo
if((dist_AB .gt. 1d-13) .or. (dist_AC .gt. 1d-13) .or. use_pw) then
continue
else
NAI_pol_mult_cgtos = V_n_e_cgtos(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 = (Ae_center(1) - Be_center(1)) * (Ae_center(1) - Be_center(1)) &
+ (Ae_center(2) - Be_center(2)) * (Ae_center(2) - Be_center(2)) &
+ (Ae_center(3) - Be_center(3)) * (Ae_center(3) - Be_center(3))
const_factor = dist * rho
if(real(const_factor) > 80.d0) then
NAI_pol_mult_cgtos = (0.d0, 0.d0)
return
endif
P_center(1) = (alpha * Ae_center(1) + beta * Be_center(1)) * p_inv
P_center(2) = (alpha * Ae_center(2) + beta * Be_center(2)) * p_inv
P_center(3) = (alpha * Ae_center(3) + beta * Be_center(3)) * p_inv
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))
const = p * dist_integral
factor = zexp(-const_factor)
coeff = dtwo_pi * factor * p_inv
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_cgtos = coeff * crint_1(0, const)
return
endif
d(0:n_pt_in) = (0.d0, 0.d0)
call give_cpolynomial_mult_center_one_e(Ap_center, Bp_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_cgtos = (0.d0, 0.d0)
return
endif
NAI_pol_mult_cgtos = coeff * crint_sum(n_pt_out, const, d)
return
end
! ---
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) :: C_center(3)
complex*16, intent(in) :: alpha, beta, A_center(3), B_center(3)
integer, intent(out) :: n_pt_out
complex*16, intent(inout) :: 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
P_center(1) = (alpha * A_center(1) + beta * B_center(1)) * p_inv
P_center(2) = (alpha * A_center(2) + beta * B_center(2)) * p_inv
P_center(3) = (alpha * A_center(3) + beta * B_center(3)) * p_inv
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_cgtos(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_cgtos(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_cgtos(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)
n_pt_out = 0
d1(0:n_pt_tmp) = (0.d0, 0.d0)
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
! ---
recursive subroutine I_x1_pol_mult_one_e_cgtos(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_cgtos(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_cgtos(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_cgtos(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_cgtos(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_cgtos(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_cgtos(a-1, c, R1x, R1xp, R2x, Y, ny, n_pt_in)
call multiply_cpoly(Y, ny, R1x, 2, d, nd)
endif
end
! ---
recursive subroutine I_x2_pol_mult_one_e_cgtos(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_cgtos(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_cgtos(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
! ---
complex*16 function V_n_e_cgtos(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_cgtos
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_cgtos = (0.d0, 0.d0)
else
V_n_e_cgtos = V_r_cgtos(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
! ---
complex*16 function V_r_cgtos(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_cgtos = 0.5d0 * fact(shiftr(n, 1)) / (alpha**(shiftr(n, 1) + 1))
else
V_r_cgtos = sqpi * fact(n) / fact(shiftr(n, 1)) * (0.5d0/zsqrt(alpha))**(n+1)
endif
end
! ---

318
src/ao_one_e_ints/one_e_kin_integrals_cgtos.irp.f
View File

@ -1,318 +0,0 @@
! ---
BEGIN_PROVIDER [double precision, ao_deriv2_cgtos_x, (ao_num, ao_num)]
&BEGIN_PROVIDER [double precision, ao_deriv2_cgtos_y, (ao_num, ao_num)]
&BEGIN_PROVIDER [double precision, ao_deriv2_cgtos_z, (ao_num, ao_num)]
implicit none
integer :: i, j, m, n, l, ii, jj, dim1, power_A(3), power_B(3)
double precision :: c, deriv_tmp
double precision :: KA2, phiA
double precision :: KB2, phiB
double precision :: aa
complex*16 :: alpha, alpha_inv, Ae_center(3), Ap_center(3), C1
complex*16 :: beta, beta_inv, Be_center(3), Bp_center(3), C2
complex*16 :: xa
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_m1_1, overlap_p1_1, overlap_p2_1
complex*16 :: overlap_m2_2, overlap_m1_2, overlap_p1_2, overlap_p2_2
complex*16 :: deriv_tmp_1, deriv_tmp_2
dim1 = 100
! -- Dummy call to provide everything
Ae_center(:) = (0.0d0, 0.d0)
Be_center(:) = (1.0d0, 0.d0)
Ap_center(:) = (0.0d0, 0.d0)
Bp_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(Ae_center, Be_center, alpha, beta, power_A, power_B, &
Ap_center, Bp_center, overlap_x0_1, overlap_y0_1, overlap_z0_1, overlap, dim1)
! ---
!$OMP PARALLEL &
!$OMP DEFAULT(NONE) &
!$OMP PRIVATE(i, j, m, n, l, ii, jj, c, aa, xa, C1, C2, &
!$OMP Ae_center, power_A, alpha, alpha_inv, KA2, phiA, Ap_center, &
!$OMP Be_center, power_B, beta, beta_inv, KB2, phiB, Bp_center, &
!$OMP deriv_tmp, deriv_tmp_1, deriv_tmp_2, &
!$OMP overlap_x, overlap_y, overlap_z, overlap, &
!$OMP overlap_m2_1, overlap_m1_1, overlap_p1_1, overlap_p2_1, &
!$OMP overlap_m2_2, overlap_m1_2, overlap_p1_2, overlap_p2_2, &
!$OMP overlap_x0_1, overlap_y0_1, overlap_z0_1, overlap_x0_2, &
!$OMP overlap_y0_2, overlap_z0_2) &
!$OMP SHARED(nucl_coord, ao_power, ao_prim_num, ao_num, ao_nucl, dim1, &
!$OMP ao_coef_cgtos_norm_ord_transp, ao_expo_cgtos_ord_transp, &
!$OMP ao_expo_pw_ord_transp, ao_expo_phase_ord_transp, &
!$OMP ao_deriv2_cgtos_x, ao_deriv2_cgtos_y, ao_deriv2_cgtos_z)
!$OMP DO SCHEDULE(GUIDED)
do j = 1, ao_num
jj = ao_nucl(j)
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
ii = ao_nucl(i)
power_B(1) = ao_power(i,1)
power_B(2) = ao_power(i,2)
power_B(3) = ao_power(i,3)
ao_deriv2_cgtos_x(i,j) = 0.d0
ao_deriv2_cgtos_y(i,j) = 0.d0
ao_deriv2_cgtos_z(i,j) = 0.d0
do n = 1, ao_prim_num(j)
alpha = ao_expo_cgtos_ord_transp(n,j)
alpha_inv = (1.d0, 0.d0) / alpha
do m = 1, 3
Ap_center(m) = nucl_coord(jj,m)
Ae_center(m) = nucl_coord(jj,m) - (0.d0, 0.5d0) * alpha_inv * ao_expo_pw_ord_transp(m,n,j)
enddo
phiA = ao_expo_phase_ord_transp(4,n,j)
KA2 = ao_expo_pw_ord_transp(4,n,j)
do l = 1, ao_prim_num(i)
beta = ao_expo_cgtos_ord_transp(l,i)
beta_inv = (1.d0, 0.d0) / beta
do m = 1, 3
Bp_center(m) = nucl_coord(ii,m)
Be_center(m) = nucl_coord(ii,m) - (0.d0, 0.5d0) * beta_inv * ao_expo_pw_ord_transp(m,l,i)
enddo
phiB = ao_expo_phase_ord_transp(4,l,i)
KB2 = ao_expo_pw_ord_transp(4,l,i)
c = ao_coef_cgtos_norm_ord_transp(n,j) * ao_coef_cgtos_norm_ord_transp(l,i)
C1 = zexp((0.d0, 1.d0) * (-phiA - phiB) - 0.25d0 * (alpha_inv * KA2 + beta_inv * KB2))
C2 = zexp((0.d0, 1.d0) * (-phiA + phiB) - 0.25d0 * (alpha_inv * KA2 + conjg(beta_inv) * KB2))
call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, &
Ap_center, Bp_center, overlap_x0_1, overlap_y0_1, overlap_z0_1, overlap, dim1)
call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, &
Ap_center, conjg(Bp_center), overlap_x0_2, overlap_y0_2, overlap_z0_2, overlap, dim1)
! ---
power_A(1) = power_A(1) - 1
if(power_A(1) > -1) then
call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, &
Ap_center, Bp_center, overlap_m1_1, overlap_y, overlap_z, overlap, dim1)
call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, &
Ap_center, conjg(Bp_center), overlap_m1_2, overlap_y, overlap_z, overlap, dim1)
power_A(1) = power_A(1) - 1
if(power_A(1) > -1) then
call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, &
Ap_center, Bp_center, overlap_m2_1, overlap_y, overlap_z, overlap, dim1)
call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, &
Ap_center, conjg(Bp_center), 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) + 1
else
overlap_m1_1 = (0.d0, 0.d0)
overlap_m1_2 = (0.d0, 0.d0)
overlap_m2_1 = (0.d0, 0.d0)
overlap_m2_2 = (0.d0, 0.d0)
endif
power_A(1) = power_A(1) + 1
power_A(1) = power_A(1) + 1
call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, &
Ap_center, Bp_center, overlap_p1_1, overlap_y, overlap_z, overlap, dim1)
call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, &
Ap_center, conjg(Bp_center), overlap_p1_2, overlap_y, overlap_z, overlap, dim1)
power_A(1) = power_A(1) + 1
call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, &
Ap_center, Bp_center, overlap_p2_1, overlap_y, overlap_z, overlap, dim1)
call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, &
Ap_center, conjg(Bp_center), overlap_p2_2, overlap_y, overlap_z, overlap, dim1)
power_A(1) = power_A(1) - 2
aa = dble(power_A(1))
xa = Ap_center(1) - Ae_center(1)
deriv_tmp_1 = aa * (aa - 1.d0) * overlap_m2_1 - 4.d0 * alpha * aa * xa * overlap_m1_1 &
+ 4.d0 * alpha * (alpha * xa * xa - aa - 0.5d0) * overlap_x0_1 &
+ 8.d0 * alpha * alpha * (xa * overlap_p1_1 + 0.5d0 * overlap_p2_1)
deriv_tmp_1 = deriv_tmp_1 * overlap_y0_1 * overlap_z0_1
deriv_tmp_2 = aa * (aa - 1.d0) * overlap_m2_2 - 4.d0 * alpha * aa * xa * overlap_m1_2 &
+ 4.d0 * alpha * (alpha * xa * xa - aa - 0.5d0) * overlap_x0_2 &
+ 8.d0 * alpha * alpha * (xa * overlap_p1_2 + 0.5d0 * overlap_p2_2)
deriv_tmp_2 = deriv_tmp_2 * overlap_y0_2 * overlap_z0_2
deriv_tmp = 2.d0 * real(C1 * deriv_tmp_1 + C2 * deriv_tmp_2)
ao_deriv2_cgtos_x(i,j) += c * deriv_tmp
! ---
power_A(2) = power_A(2) - 1
if(power_A(2) > -1) then
call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, &
Ap_center, Bp_center, overlap_x, overlap_m1_1, overlap_z, overlap, dim1)
call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, &
Ap_center, conjg(Bp_center), overlap_x, overlap_m1_2, overlap_z, overlap, dim1)
power_A(2) = power_A(2) - 1
if(power_A(2) > -1) then
call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, &
Ap_center, Bp_center, overlap_x, overlap_m2_1, overlap_z, overlap, dim1)
call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, &
Ap_center, conjg(Bp_center), overlap_x, overlap_m2_2, overlap_z, overlap, dim1)
else
overlap_m2_1 = (0.d0, 0.d0)
overlap_m2_2 = (0.d0, 0.d0)
endif
power_A(2) = power_A(2) + 1
else
overlap_m1_1 = (0.d0, 0.d0)
overlap_m1_2 = (0.d0, 0.d0)
overlap_m2_1 = (0.d0, 0.d0)
overlap_m2_2 = (0.d0, 0.d0)
endif
power_A(2) = power_A(2) + 1
power_A(2) = power_A(2) + 1
call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, &
Ap_center, Bp_center, overlap_x, overlap_p1_1, overlap_z, overlap, dim1)
call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, &
Ap_center, conjg(Bp_center), overlap_x, overlap_p1_2, overlap_z, overlap, dim1)
power_A(2) = power_A(2) + 1
call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, &
Ap_center, Bp_center, overlap_x, overlap_p2_1, overlap_z, overlap, dim1)
call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, &
Ap_center, conjg(Bp_center), overlap_x, overlap_p2_2, overlap_z, overlap, dim1)
power_A(2) = power_A(2) - 2
aa = dble(power_A(2))
xa = Ap_center(2) - Ae_center(2)
deriv_tmp_1 = aa * (aa - 1.d0) * overlap_m2_1 - 4.d0 * alpha * aa * xa * overlap_m1_1 &
+ 4.d0 * alpha * (alpha * xa * xa - aa - 0.5d0) * overlap_y0_1 &
+ 8.d0 * alpha * alpha * (xa * overlap_p1_1 + 0.5d0 * overlap_p2_1)
deriv_tmp_1 = deriv_tmp_1 * overlap_x0_1 * overlap_z0_1
deriv_tmp_2 = aa * (aa - 1.d0) * overlap_m2_2 - 4.d0 * alpha * aa * xa * overlap_m1_2 &
+ 4.d0 * alpha * (alpha * xa * xa - aa - 0.5d0) * overlap_y0_2 &
+ 8.d0 * alpha * alpha * (xa * overlap_p1_2 + 0.5d0 * overlap_p2_2)
deriv_tmp_2 = deriv_tmp_2 * overlap_x0_2 * overlap_z0_2
deriv_tmp = 2.d0 * real(C1 * deriv_tmp_1 + C2 * deriv_tmp_2)
ao_deriv2_cgtos_y(i,j) += c * deriv_tmp
! ---
power_A(3) = power_A(3) - 1
if(power_A(3) > -1) then
call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, &
Ap_center, Bp_center, overlap_x, overlap_y, overlap_m1_1, overlap, dim1)
call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, &
Ap_center, conjg(Bp_center), overlap_x, overlap_y, overlap_m1_2, overlap, dim1)
power_A(3) = power_A(3) - 1
if(power_A(3) > -1) then
call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, &
Ap_center, Bp_center, overlap_x, overlap_y, overlap_m2_1, overlap, dim1)
call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, &
Ap_center, conjg(Bp_center), 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) + 1
else
overlap_m1_1 = (0.d0, 0.d0)
overlap_m1_2 = (0.d0, 0.d0)
overlap_m2_1 = (0.d0, 0.d0)
overlap_m2_2 = (0.d0, 0.d0)
endif
power_A(3) = power_A(3) + 1
power_A(3) = power_A(3) + 1
call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, &
Ap_center, Bp_center, overlap_x, overlap_y, overlap_p1_1, overlap, dim1)
call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, &
Ap_center, conjg(Bp_center), overlap_x, overlap_y, overlap_p1_2, overlap, dim1)
power_A(3) = power_A(3) + 1
call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, &
Ap_center, Bp_center, overlap_x, overlap_y, overlap_p2_1, overlap, dim1)
call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, &
Ap_center, conjg(Bp_center), overlap_x, overlap_y, overlap_p2_2, overlap, dim1)
power_A(3) = power_A(3) - 2
aa = dble(power_A(3))
xa = Ap_center(3) - Ae_center(3)
deriv_tmp_1 = aa * (aa - 1.d0) * overlap_m2_1 - 4.d0 * alpha * aa * xa * overlap_m1_1 &
+ 4.d0 * alpha * (alpha * xa * xa - aa - 0.5d0) * overlap_z0_1 &
+ 8.d0 * alpha * alpha * (xa * overlap_p1_1 + 0.5d0 * overlap_p2_1)
deriv_tmp_1 = deriv_tmp_1 * overlap_x0_1 * overlap_y0_1
deriv_tmp_2 = aa * (aa - 1.d0) * overlap_m2_2 - 4.d0 * alpha * aa * xa * overlap_m1_2 &
+ 4.d0 * alpha * (alpha * xa * xa - aa - 0.5d0) * overlap_z0_2 &
+ 8.d0 * alpha * alpha * (xa * overlap_p1_2 + 0.5d0 * overlap_p2_2)
deriv_tmp_2 = deriv_tmp_2 * overlap_x0_2 * overlap_y0_2
deriv_tmp = 2.d0 * real(C1 * deriv_tmp_1 + C2 * deriv_tmp_2)
ao_deriv2_cgtos_z(i,j) += c * deriv_tmp
! ---
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
END_PROVIDER
! ---
BEGIN_PROVIDER [double precision, ao_kinetic_integrals_cgtos, (ao_num, ao_num)]
BEGIN_DOC
!
! Kinetic energy integrals in the cgtos |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_cgtos, ao_deriv2_cgtos_x, ao_deriv2_cgtos_y, ao_deriv2_cgtos_z)
do j = 1, ao_num
do i = 1, ao_num
ao_kinetic_integrals_cgtos(i,j) = -0.5d0 * (ao_deriv2_cgtos_x(i,j) + &
ao_deriv2_cgtos_y(i,j) + &
ao_deriv2_cgtos_z(i,j))
enddo
enddo
!$OMP END PARALLEL DO
END_PROVIDER
! ---

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

@ -9,13 +9,11 @@ subroutine give_all_erf_kl_ao(integrals_ao,mu_in,C_center)
END_DOC
double precision, intent(in) :: mu_in,C_center(3)
double precision, intent(out) :: integrals_ao(ao_num,ao_num)
double precision :: NAI_pol_mult_erf_ao
integer :: i,j,l,k,m
do k = 1, ao_num
do m = 1, ao_num
integrals_ao(m,k) = NAI_pol_mult_erf_ao(m,k,mu_in,C_center)
enddo
enddo
double precision, allocatable :: integrals_ao_cart(:,:)
allocate(integrals_ao_cart(ao_cart_num,ao_cart_num))
call give_all_erf_kl_ao_cart(integrals_ao_cart,mu_in,C_center)
call ao_cart_to_ao_basis(integrals_ao_cart, ao_cart_num, integrals_ao, ao_num)
end
! ---
@ -33,35 +31,39 @@ double precision function NAI_pol_mult_erf_ao(i_ao, j_ao, mu_in, C_center)
integer, intent(in) :: i_ao, j_ao
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
double precision :: NAI_pol_mult_erf
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)
num_B = ao_nucl(j_ao)
power_B(1:3) = ao_power(j_ao,1:3)
B_center(1:3) = nucl_coord(num_B,1:3)
n_pt_in = n_pt_max_integrals
NAI_pol_mult_erf_ao = 0.d0
do i = 1, ao_prim_num(i_ao)
alpha = ao_expo_ordered_transp(i,i_ao)
do j = 1, ao_prim_num(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)
enddo
enddo
double precision, allocatable :: integrals_ao(:,:),integrals_ao_cart(:,:)
allocate(integrals_ao(ao_num,ao_num),integrals_ao_cart(ao_cart_num,ao_cart_num))
call give_all_erf_kl_ao_cart(integrals_ao_cart,mu_in,C_center)
call ao_cart_to_ao_basis(integrals_ao_cart, ao_cart_num, integrals_ao, ao_num)
NAI_pol_mult_erf_ao = integrals_ao(i_ao, j_ao)
end function NAI_pol_mult_erf_ao
! ---
subroutine all_NAI_pol_mult_erf_ao_with1s(beta, B_center, mu_in, C_center)
BEGIN_DOC
!
! Computes ALL 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
double precision, intent(in) :: beta, B_center(3)
double precision, intent(in) :: mu_in, C_center(3)
double precision, intent(out) :: integrals_ao(ao_num,ao_num)
double precision :: NAI_pol_mult_erf_ao_cart_with1s
integer :: i,j
double precision, allocatable :: integrals_ao_cart(:,:)
allocate(integrals_ao_cart(ao_cart_num,ao_cart_num))
do i = 1, ao_cart_num
do j = 1, ao_cart_num
integrals_ao_car(j,i) = NAI_pol_mult_erf_ao_cart_with1s(i_ao, j_ao, beta, B_center, mu_in, C_center)
enddo
enddo
call ao_cart_to_ao_basis(integrals_ao_cart, ao_cart_num, integrals_ao, ao_num)
end
double precision function NAI_pol_mult_erf_ao_with1s(i_ao, j_ao, beta, B_center, mu_in, C_center)
@ -76,719 +78,17 @@ double precision function NAI_pol_mult_erf_ao_with1s(i_ao, j_ao, beta, B_center,
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
double precision :: NAI_pol_mult_erf_ao_cart_with1s
integer :: i,j
double precision, allocatable :: integrals_ao(:,:),integrals_ao_cart(:,:)
allocate(integrals_ao(ao_num,ao_num),integrals_ao_cart(ao_cart_num,ao_cart_num))
do i = 1, ao_cart_num
do j = 1, ao_cart_num
integrals_ao_car(j,i) = NAI_pol_mult_erf_ao_cart_with1s(i_ao, j_ao, beta, B_center, mu_in, C_center)
enddo
enddo
call ao_cart_to_ao_basis(integrals_ao_cart, ao_cart_num, integrals_ao, ao_num)
NAI_pol_mult_erf_ao_with1s = integrals_ao(i_ao, j_ao)
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)
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
integer, intent(in) :: power_A(3), power_B(3)
double precision, intent(in) :: C_center(3), A_center(3), B_center(3), alpha, beta, mu_in
integer :: i, n_pt, n_pt_out
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
double precision :: rint
p = alpha + beta
p_inv = 1.d0 / p
p_inv_2 = 0.5d0 * p_inv
rho = alpha * beta * p_inv
dist = 0.d0
dist_integral = 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
if(const_factor > 80.d0) then
NAI_pol_mult_erf = 0.d0
return
endif
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_A(1) + power_B(1)) + (power_A(2) + power_B(2)) + (power_A(3) + power_B(3)) )
const = p * dist_integral * p_new * p_new
if(n_pt == 0) then
NAI_pol_mult_erf = coeff * rint(0, const)
return
endif
do i = 0, n_pt_in
d(i) = 0.d0
enddo
! call 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)
p_new = p_new * p_new
call give_polynomial_mult_center_one_e_erf_opt(A_center, B_center, power_A, power_B, 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 = 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 = accu * coeff
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)
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
! ---
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
! ---
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)
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, LD_B, LD_C, LD_resv, n_points
integer, intent(in) :: power_A1(3), power_A2(3)
double precision, intent(in) :: A1_center(3), A2_center(3)
double precision, intent(in) :: C_center(LD_C,3), B_center(LD_B,3)
double precision, intent(in) :: alpha1, alpha2, beta, mu_in
double precision, intent(out) :: res_v(LD_resv)
integer :: i, n_pt, n_pt_out, ipoint
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, p_new2, coef_tmp, cons_tmp
double precision :: rint
res_V(1:LD_resv) = 0.d0
! 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
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_new = mu_in / dsqrt(p + mu_in * mu_in)
p_new2 = p_new * p_new
coef_tmp = dtwo_pi * p_inv * p_new
cons_tmp = p * p_new2
n_pt = 2 * (power_A1(1) + power_A2(1) + power_A1(2) + power_A2(2) + power_A1(3) + power_A2(3) )
if(n_pt == 0) then
do ipoint = 1, n_points
dist = (A12_center(1) - B_center(ipoint,1)) * (A12_center(1) - B_center(ipoint,1))&
+ (A12_center(2) - B_center(ipoint,2)) * (A12_center(2) - B_center(ipoint,2))&
+ (A12_center(3) - B_center(ipoint,3)) * (A12_center(3) - B_center(ipoint,3))
const_factor = const_factor12 + dist * rho
if(const_factor > 80.d0) cycle
coeff = coef_tmp * dexp(-const_factor)
P_center(1) = (alpha12 * A12_center(1) + beta * B_center(ipoint,1)) * p_inv
P_center(2) = (alpha12 * A12_center(2) + beta * B_center(ipoint,2)) * p_inv
P_center(3) = (alpha12 * A12_center(3) + beta * B_center(ipoint,3)) * p_inv
dist_integral = (P_center(1) - C_center(ipoint,1)) * (P_center(1) - C_center(ipoint,1))&
+ (P_center(2) - C_center(ipoint,2)) * (P_center(2) - C_center(ipoint,2))&
+ (P_center(3) - C_center(ipoint,3)) * (P_center(3) - C_center(ipoint,3))
const = cons_tmp * dist_integral
res_v(ipoint) = coeff * rint(0, const)
enddo
else
do ipoint = 1, n_points
dist = (A12_center(1) - B_center(ipoint,1)) * (A12_center(1) - B_center(ipoint,1))&
+ (A12_center(2) - B_center(ipoint,2)) * (A12_center(2) - B_center(ipoint,2))&
+ (A12_center(3) - B_center(ipoint,3)) * (A12_center(3) - B_center(ipoint,3))
const_factor = const_factor12 + dist * rho
if(const_factor > 80.d0) cycle
coeff = coef_tmp * dexp(-const_factor)
P_center(1) = (alpha12 * A12_center(1) + beta * B_center(ipoint,1)) * p_inv
P_center(2) = (alpha12 * A12_center(2) + beta * B_center(ipoint,2)) * p_inv
P_center(3) = (alpha12 * A12_center(3) + beta * B_center(ipoint,3)) * p_inv
dist_integral = (P_center(1) - C_center(ipoint,1)) * (P_center(1) - C_center(ipoint,1))&
+ (P_center(2) - C_center(ipoint,2)) * (P_center(2) - C_center(ipoint,2))&
+ (P_center(3) - C_center(ipoint,3)) * (P_center(3) - C_center(ipoint,3))
const = cons_tmp * dist_integral
do i = 0, n_pt_in
d(i) = 0.d0
enddo
!TODO: VECTORIZE HERE
call give_polynomial_mult_center_one_e_erf_opt(A1_center, A2_center, power_A1, power_A2, 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
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_with1s_v
! ---
subroutine give_polynomial_mult_center_one_e_erf_opt(A_center, B_center, power_A, power_B, C_center, n_pt_in, d, n_pt_out, p_inv_2, p_new, P_center)
BEGIN_DOC
! Returns the explicit polynomial in terms of the $t$ variable of the
! following polynomial:
!
! $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
integer, intent(in) :: power_A(3), power_B(3)
double precision, intent(in) :: A_center(3), B_center(3), C_center(3), p_inv_2, p_new, P_center(3)
integer, intent(out) :: n_pt_out
double precision, 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
integer :: n_pt_tmp
double precision :: d1(0:n_pt_in)
double precision :: d2(0:n_pt_in)
double precision :: d3(0:n_pt_in)
double precision :: accu
double precision :: R1x(0:2), B01(0:2), R1xp(0:2), R2x(0:2)
accu = 0.d0
ASSERT (n_pt_in > 1)
R1x(0) = (P_center(1) - A_center(1))
R1x(1) = 0.d0
R1x(2) = -(P_center(1) - C_center(1))* p_new
! R1x = (P_x - A_x) - (P_x - C_x) ( t * mu/sqrt(p+mu^2) )^2
R1xp(0) = (P_center(1) - B_center(1))
R1xp(1) = 0.d0
R1xp(2) =-(P_center(1) - C_center(1))* p_new
!R1xp = (P_x - B_x) - (P_x - C_x) ( t * mu/sqrt(p+mu^2) )^2
R2x(0) = p_inv_2
R2x(1) = 0.d0
R2x(2) = -p_inv_2 * p_new
!R2x = 0.5 / p - 0.5/p ( t * mu/sqrt(p+mu^2) )^2
do i = 0, n_pt_in
d (i) = 0.d0
d1(i) = 0.d0
d2(i) = 0.d0
d3(i) = 0.d0
enddo
n_pt1 = n_pt_in
n_pt2 = n_pt_in
n_pt3 = n_pt_in
a_x = power_A(1)
b_x = power_B(1)
call I_x1_pol_mult_one_e(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
enddo
return
endif
R1x(0) = (P_center(2) - A_center(2))
R1x(1) = 0.d0
R1x(2) = -(P_center(2) - C_center(2))* p_new
! R1x = (P_x - A_x) - (P_x - C_x) ( t * mu/sqrt(p+mu^2) )^2
R1xp(0) = (P_center(2) - B_center(2))
R1xp(1) = 0.d0
R1xp(2) =-(P_center(2) - C_center(2))* p_new
!R1xp = (P_x - B_x) - (P_x - C_x) ( t * mu/sqrt(p+mu^2) )^2
a_y = power_A(2)
b_y = power_B(2)
call I_x1_pol_mult_one_e(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
enddo
return
endif
R1x(0) = (P_center(3) - A_center(3))
R1x(1) = 0.d0
R1x(2) = -(P_center(3) - C_center(3)) * p_new
! R1x = (P_x - A_x) - (P_x - C_x) ( t * mu/sqrt(p+mu^2) )^2
R1xp(0) = (P_center(3) - B_center(3))
R1xp(1) = 0.d0
R1xp(2) =-(P_center(3) - C_center(3)) * p_new
!R2x = 0.5 / p - 0.5/p ( t * mu/sqrt(p+mu^2) )^2
a_z = power_A(3)
b_z = power_B(3)
call I_x1_pol_mult_one_e(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
enddo
return
endif
n_pt_tmp = 0
call multiply_poly(d1, n_pt1, d2, n_pt2, d, n_pt_tmp)
do i = 0, n_pt_tmp
d1(i) = 0.d0
enddo
n_pt_out = 0
call multiply_poly(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_polynomial_mult_center_one_e_erf_opt
! ---
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)
BEGIN_DOC
! Returns the explicit polynomial in terms of the $t$ variable of the
! following polynomial:
!
! $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
integer,intent(out) :: n_pt_out
double precision, intent(in) :: A_center(3), B_center(3),C_center(3)
double precision, intent(in) :: alpha,beta,mu_in
integer, intent(in) :: power_A(3), power_B(3)
integer :: a_x,b_x,a_y,b_y,a_z,b_z
double precision :: d(0:n_pt_in)
double precision :: d1(0:n_pt_in)
double precision :: d2(0:n_pt_in)
double precision :: d3(0:n_pt_in)
double precision :: accu, pq_inv, p10_1, p10_2, p01_1, p01_2
double precision :: p,P_center(3),rho,p_inv,p_inv_2
accu = 0.d0
ASSERT (n_pt_in > 1)
p = alpha+beta
p_inv = 1.d0/p
p_inv_2 = 0.5d0/p
do i =1, 3
P_center(i) = (alpha * A_center(i) + beta * B_center(i)) * p_inv
enddo
double precision :: R1x(0:2), B01(0:2), R1xp(0:2),R2x(0:2)
R1x(0) = (P_center(1) - A_center(1))
R1x(1) = 0.d0
R1x(2) = -(P_center(1) - C_center(1))* mu_in**2 / (p+mu_in*mu_in)
! R1x = (P_x - A_x) - (P_x - C_x) ( t * mu/sqrt(p+mu^2) )^2
R1xp(0) = (P_center(1) - B_center(1))
R1xp(1) = 0.d0
R1xp(2) =-(P_center(1) - C_center(1))* mu_in**2 / (p+mu_in*mu_in)
!R1xp = (P_x - B_x) - (P_x - C_x) ( t * mu/sqrt(p+mu^2) )^2
R2x(0) = p_inv_2
R2x(1) = 0.d0
R2x(2) = -p_inv_2* mu_in**2 / (p+mu_in*mu_in)
!R2x = 0.5 / p - 0.5/p ( t * mu/sqrt(p+mu^2) )^2
do i = 0,n_pt_in
d(i) = 0.d0
enddo
do i = 0,n_pt_in
d1(i) = 0.d0
enddo
do i = 0,n_pt_in
d2(i) = 0.d0
enddo
do i = 0,n_pt_in
d3(i) = 0.d0
enddo
integer :: n_pt1,n_pt2,n_pt3,dim,i
n_pt1 = n_pt_in
n_pt2 = n_pt_in
n_pt3 = n_pt_in
a_x = power_A(1)
b_x = power_B(1)
call I_x1_pol_mult_one_e(a_x,b_x,R1x,R1xp,R2x,d1,n_pt1,n_pt_in)
! print*,'passed the first I_x1'
if(n_pt1<0)then
n_pt_out = -1
do i = 0,n_pt_in
d(i) = 0.d0
enddo
return
endif
R1x(0) = (P_center(2) - A_center(2))
R1x(1) = 0.d0
R1x(2) = -(P_center(2) - C_center(2))* mu_in**2 / (p+mu_in*mu_in)
! R1x = (P_x - A_x) - (P_x - C_x) ( t * mu/sqrt(p+mu^2) )^2
R1xp(0) = (P_center(2) - B_center(2))
R1xp(1) = 0.d0
R1xp(2) =-(P_center(2) - C_center(2))* mu_in**2 / (p+mu_in*mu_in)
!R1xp = (P_x - B_x) - (P_x - C_x) ( t * mu/sqrt(p+mu^2) )^2
a_y = power_A(2)
b_y = power_B(2)
call I_x1_pol_mult_one_e(a_y,b_y,R1x,R1xp,R2x,d2,n_pt2,n_pt_in)
! print*,'passed the second I_x1'
if(n_pt2<0)then
n_pt_out = -1
do i = 0,n_pt_in
d(i) = 0.d0
enddo
return
endif
R1x(0) = (P_center(3) - A_center(3))
R1x(1) = 0.d0
R1x(2) = -(P_center(3) - C_center(3))* mu_in**2 / (p+mu_in*mu_in)
! R1x = (P_x - A_x) - (P_x - C_x) ( t * mu/sqrt(p+mu^2) )^2
R1xp(0) = (P_center(3) - B_center(3))
R1xp(1) = 0.d0
R1xp(2) =-(P_center(3) - C_center(3))* mu_in**2 / (p+mu_in*mu_in)
!R2x = 0.5 / p - 0.5/p ( t * mu/sqrt(p+mu^2) )^2
a_z = power_A(3)
b_z = power_B(3)
! print*,'a_z = ',a_z
! print*,'b_z = ',b_z
call I_x1_pol_mult_one_e(a_z,b_z,R1x,R1xp,R2x,d3,n_pt3,n_pt_in)
! print*,'passed the third I_x1'
if(n_pt3<0)then
n_pt_out = -1
do i = 0,n_pt_in
d(i) = 0.d0
enddo
return
endif
integer :: n_pt_tmp
n_pt_tmp = 0
call multiply_poly(d1,n_pt1,d2,n_pt2,d,n_pt_tmp)
do i = 0,n_pt_tmp
d1(i) = 0.d0
enddo
n_pt_out = 0
call multiply_poly(d ,n_pt_tmp ,d3,n_pt3,d1,n_pt_out)
do i = 0, n_pt_out
d(i) = d1(i)
enddo
end

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

@ -27,78 +27,7 @@ BEGIN_PROVIDER [ double precision, ao_integrals_n_e, (ao_num,ao_num)]
else
if(use_cgtos) then
do j = 1, ao_num
do i = 1, ao_num
ao_integrals_n_e(i,j) = ao_integrals_n_e_cgtos(i,j)
enddo
enddo
else
!$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)
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_ordered_transp(l,j)
do m=1,ao_prim_num(i)
beta = ao_expo_ordered_transp(m,i)
double precision :: c, c1
c = 0.d0
do k = 1, nucl_num
double precision :: Z
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 *, 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
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
endif
call ao_cart_to_ao_basis(ao_cart_coul_n_e, ao_cart_num, ao_integrals_n_e, ao_num)
IF(do_pseudo) THEN
ao_integrals_n_e += ao_pseudo_integrals
@ -147,487 +76,9 @@ 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`
END_DOC
implicit none
double precision :: alpha, beta
integer :: i_c,num_A,num_B
double precision :: A_center(3),B_center(3),C_center(3)
integer :: power_A(3),power_B(3)
integer :: i,j,k,l,n_pt_in,m
double precision :: overlap_x,overlap_y,overlap_z,overlap,dx,NAI_pol_mult
ao_integrals_n_e_per_atom = 0.d0
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,l,m,alpha,beta,A_center,B_center,power_A,power_B,&
!$OMP num_A,num_B,c,n_pt_in,C_center) &
!$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_per_atom,nucl_num)
n_pt_in = n_pt_max_integrals
!$OMP DO SCHEDULE (dynamic)
double precision :: c
do j = 1, ao_num
power_A(1)= ao_power(j,1)
power_A(2)= ao_power(j,2)
power_A(3)= ao_power(j,3)
num_A = ao_nucl(j)
A_center(1) = nucl_coord(num_A,1)
A_center(2) = nucl_coord(num_A,2)
A_center(3) = nucl_coord(num_A,3)
do k = 1, nucl_num
C_center(1) = nucl_coord(k,1)
C_center(2) = nucl_coord(k,2)
C_center(3) = nucl_coord(k,3)
do i = 1, ao_num
power_B(1)= ao_power(i,1)
power_B(2)= ao_power(i,2)
power_B(3)= ao_power(i,3)
num_B = ao_nucl(i)
B_center(1) = nucl_coord(num_B,1)
B_center(2) = nucl_coord(num_B,2)
B_center(3) = nucl_coord(num_B,3)
c = 0.d0
do l=1,ao_prim_num(j)
alpha = ao_expo_ordered_transp(l,j)
do m=1,ao_prim_num(i)
beta = ao_expo_ordered_transp(m,i)
c = c + NAI_pol_mult(A_center,B_center,power_A,power_B, &
alpha,beta,C_center,n_pt_in) &
* ao_coef_normalized_ordered_transp(l,j) &
* ao_coef_normalized_ordered_transp(m,i)
enddo
enddo
ao_integrals_n_e_per_atom(i,j,k) = -c
enddo
enddo
integer :: i
do i = 1, nucl_num
call ao_cart_to_ao_basis(ao_cart_coul_n_e_per_atom(1,1,i), ao_cart_num,ao_integrals_n_e_per_atom(1,1,i), ao_num)
enddo
!$OMP END DO
!$OMP END PARALLEL
END_PROVIDER
double precision function NAI_pol_mult(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.
!
! :math:`\langle g_i | \frac{1}{|r-R_c|} | g_j \rangle`
END_DOC
implicit none
integer, intent(in) :: n_pt_in
double precision,intent(in) :: C_center(3),A_center(3),B_center(3),alpha,beta
integer :: power_A(3),power_B(3)
integer :: i,j,k,l,n_pt
double precision :: P_center(3)
double precision :: d(0:n_pt_in),pouet,coeff,rho,dist,const,pouet_2,p,p_inv,factor
double precision :: I_n_special_exact,integrate_bourrin,I_n_bibi
double precision :: V_n_e,const_factor,dist_integral,tmp
double precision :: accu,epsilo,rint
integer :: n_pt_out,lmax
include 'utils/constants.include.F'
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 = V_n_e(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/p
rho = alpha * beta * p_inv
dist = 0.d0
dist_integral = 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
if(const_factor > 80.d0)then
NAI_pol_mult = 0.d0
return
endif
factor = dexp(-const_factor)
coeff = dtwo_pi * factor * p_inv
lmax = 20
! print*, "b"
do i = 0, n_pt_in
d(i) = 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
epsilo = 1.d0
pouet = rint(0,const)
NAI_pol_mult = coeff * pouet
return
endif
call 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)
if(n_pt_out<0)then
NAI_pol_mult = 0.d0
return
endif
accu = 0.d0
! 1/r1 standard attraction integral
epsilo = 1.d0
! sum of integrals of type : int {t,[0,1]} exp-(rho.(P-Q)^2 * t^2) * t^i
do i =0 ,n_pt_out,2
accu += d(i) * rint(i/2,const)
enddo
NAI_pol_mult = accu * coeff
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)
implicit none
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
integer, intent(in) :: n_pt_in
integer,intent(out) :: n_pt_out
double precision, intent(in) :: A_center(3), B_center(3),C_center(3)
double precision, intent(in) :: alpha,beta
integer, intent(in) :: power_A(3), power_B(3)
integer :: a_x,b_x,a_y,b_y,a_z,b_z
double precision :: d(0:n_pt_in)
double precision :: d1(0:n_pt_in)
double precision :: d2(0:n_pt_in)
double precision :: d3(0:n_pt_in)
double precision :: accu, pq_inv, p10_1, p10_2, p01_1, p01_2
double precision :: p,P_center(3),rho,p_inv,p_inv_2
accu = 0.d0
ASSERT (n_pt_in > 1)
p = alpha+beta
p_inv = 1.d0/p
p_inv_2 = 0.5d0/p
do i =1, 3
P_center(i) = (alpha * A_center(i) + beta * B_center(i)) * p_inv
enddo
double precision :: R1x(0:2), B01(0:2), R1xp(0:2),R2x(0:2)
R1x(0) = (P_center(1) - A_center(1))
R1x(1) = 0.d0
R1x(2) = -(P_center(1) - C_center(1))
R1xp(0) = (P_center(1) - B_center(1))
R1xp(1) = 0.d0
R1xp(2) =-(P_center(1) - C_center(1))
R2x(0) = p_inv_2
R2x(1) = 0.d0
R2x(2) = -p_inv_2
do i = 0,n_pt_in
d(i) = 0.d0
enddo
do i = 0,n_pt_in
d1(i) = 0.d0
enddo
do i = 0,n_pt_in
d2(i) = 0.d0
enddo
do i = 0,n_pt_in
d3(i) = 0.d0
enddo
integer :: n_pt1,n_pt2,n_pt3,dim,i
n_pt1 = n_pt_in
n_pt2 = n_pt_in
n_pt3 = n_pt_in
a_x = power_A(1)
b_x = power_B(1)
call I_x1_pol_mult_one_e(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
enddo
return
endif
R1x(0) = (P_center(2) - A_center(2))
R1x(1) = 0.d0
R1x(2) = -(P_center(2) - C_center(2))
R1xp(0) = (P_center(2) - B_center(2))
R1xp(1) = 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(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
enddo
return
endif
R1x(0) = (P_center(3) - A_center(3))
R1x(1) = 0.d0
R1x(2) = -(P_center(3) - C_center(3))
R1xp(0) = (P_center(3) - B_center(3))
R1xp(1) = 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(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
enddo
return
endif
integer :: n_pt_tmp
n_pt_tmp = 0
call multiply_poly(d1,n_pt1,d2,n_pt2,d,n_pt_tmp)
do i = 0,n_pt_tmp
d1(i) = 0.d0
enddo
n_pt_out = 0
call multiply_poly(d ,n_pt_tmp ,d3,n_pt3,d1,n_pt_out)
do i = 0, n_pt_out
d(i) = d1(i)
enddo
end
recursive subroutine I_x1_pol_mult_one_e(a,c,R1x,R1xp,R2x,d,nd,n_pt_in)
implicit none
BEGIN_DOC
! Recursive routine involved in the electron-nucleus potential
END_DOC
integer , intent(in) :: n_pt_in
double precision,intent(inout) :: d(0:n_pt_in)
integer,intent(inout) :: nd
integer, intent(in) :: a,c
double precision, intent(in) :: R1x(0:2),R1xp(0:2),R2x(0:2)
include 'utils/constants.include.F'
double precision :: X(0:max_dim)
double precision :: Y(0:max_dim)
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: X, Y
integer :: nx, ix,dim,iy,ny
dim = n_pt_in
! print*,'a,c = ',a,c
! print*,'nd_in = ',nd
if( (a==0) .and. (c==0))then
nd = 0
d(0) = 1.d0
return
elseif( (c<0).or.(nd<0) )then
nd = -1
return
else if ((a==0).and.(c.ne.0)) then
call I_x2_pol_mult_one_e(c,R1x,R1xp,R2x,d,nd,n_pt_in)
else if (a==1) then
nx = nd
do ix=0,n_pt_in
X(ix) = 0.d0
Y(ix) = 0.d0
enddo
call I_x2_pol_mult_one_e(c-1,R1x,R1xp,R2x,X,nx,n_pt_in)
do ix=0,nx
X(ix) *= dble(c)
enddo
! call multiply_poly(X,nx,R2x,2,d,nd)
call multiply_poly_c2(X,nx,R2x,d,nd)
ny=0
call I_x2_pol_mult_one_e(c,R1x,R1xp,R2x,Y,ny,n_pt_in)
! call multiply_poly(Y,ny,R1x,2,d,nd)
call multiply_poly_c2(Y,ny,R1x,d,nd)
else
do ix=0,n_pt_in
X(ix) = 0.d0
Y(ix) = 0.d0
enddo
nx = 0
call I_x1_pol_mult_one_e(a-2,c,R1x,R1xp,R2x,X,nx,n_pt_in)
do ix=0,nx
X(ix) *= dble(a-1)
enddo
! call multiply_poly(X,nx,R2x,2,d,nd)
call multiply_poly_c2(X,nx,R2x,d,nd)
nx = nd
do ix=0,n_pt_in
X(ix) = 0.d0
enddo
call I_x1_pol_mult_one_e(a-1,c-1,R1x,R1xp,R2x,X,nx,n_pt_in)
do ix=0,nx
X(ix) *= dble(c)
enddo
! call multiply_poly(X,nx,R2x,2,d,nd)
call multiply_poly_c2(X,nx,R2x,d,nd)
ny=0
call I_x1_pol_mult_one_e(a-1,c,R1x,R1xp,R2x,Y,ny,n_pt_in)
! call multiply_poly(Y,ny,R1x,2,d,nd)
call multiply_poly_c2(Y,ny,R1x,d,nd)
endif
end
recursive subroutine I_x2_pol_mult_one_e(c,R1x,R1xp,R2x,d,nd,dim)
implicit none
BEGIN_DOC
! Recursive routine involved in the electron-nucleus potential
END_DOC
integer , intent(in) :: dim
include 'utils/constants.include.F'
double precision :: d(0:max_dim)
integer,intent(inout) :: nd
integer, intent(in) :: c
double precision, intent(in) :: R1x(0:2),R1xp(0:2),R2x(0:2)
integer :: i
if(c==0) then
nd = 0
d(0) = 1.d0
return
elseif ((nd<0).or.(c<0))then
nd = -1
return
else
integer :: nx, ix,ny
double precision :: X(0:max_dim),Y(0:max_dim)
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: X, Y
do ix=0,dim
X(ix) = 0.d0
Y(ix) = 0.d0
enddo
nx = 0
call I_x1_pol_mult_one_e(0,c-2,R1x,R1xp,R2x,X,nx,dim)
do ix=0,nx
X(ix) *= dble(c-1)
enddo
! call multiply_poly(X,nx,R2x,2,d,nd)
call multiply_poly_c2(X,nx,R2x,d,nd)
ny = 0
do ix=0,dim
Y(ix) = 0.d0
enddo
call I_x1_pol_mult_one_e(0,c-1,R1x,R1xp,R2x,Y,ny,dim)
if(ny.ge.0)then
! call multiply_poly(Y,ny,R1xp,2,d,nd)
call multiply_poly_c2(Y,ny,R1xp,d,nd)
endif
endif
end
double precision function V_n_e(a_x,a_y,a_z,b_x,b_y,b_z,alpha,beta)
implicit none
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
integer :: a_x,a_y,a_z,b_x,b_y,b_z
double precision :: alpha,beta
double precision :: V_r, V_phi, V_theta
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 = 0.d0
else
V_n_e = V_r(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
double precision function int_gaus_pol(alpha,n)
implicit none
BEGIN_DOC
! Computes the integral:
!
! $\int_{-\infty}^{\infty} x^n \exp(-\alpha x^2) dx$.
END_DOC
double precision :: alpha
integer :: n
double precision :: dble_fact
include 'utils/constants.include.F'
int_gaus_pol = 0.d0
if(iand(n,1).eq.0)then
int_gaus_pol = dsqrt(alpha/pi)
double precision :: two_alpha
two_alpha = alpha+alpha
integer :: i
do i=1,n,2
int_gaus_pol = int_gaus_pol * two_alpha
enddo
int_gaus_pol = dble_fact(n -1) / int_gaus_pol
endif
end
double precision function V_r(n,alpha)
implicit none
BEGIN_DOC
! Computes the radial part of the nuclear attraction integral:
!
! $\int_{0}^{\infty} r^n \exp(-\alpha r^2) dr$
!
END_DOC
double precision :: alpha, fact
integer :: n
include 'utils/constants.include.F'
if(iand(n,1).eq.1)then
V_r = 0.5d0 * fact(shiftr(n,1)) / (alpha ** (shiftr(n,1) + 1))
else
V_r = sqpi * fact(n) / fact(shiftr(n,1)) * (0.5d0/sqrt(alpha)) ** (n+1)
endif
end
BEGIN_PROVIDER [ double precision, ao_sphe_integrals_n_e, (ao_sphe_num,ao_sphe_num) ]
implicit none
BEGIN_DOC
! |AO| VneVne inntegrals matrix in the spherical basis set
END_DOC
double precision, allocatable :: tmp(:,:)
allocate (tmp(ao_sphe_num,ao_num))
call dgemm('N','N',ao_sphe_num,ao_num,ao_num, 1.d0, &
ao_cart_to_sphe_inv,size(ao_cart_to_sphe_inv,1), &
ao_integrals_n_e,size(ao_integrals_n_e,1), 0.d0, &
tmp, size(tmp,1))
call dgemm('N','T',ao_sphe_num,ao_sphe_num,ao_num, 1.d0, &
tmp, size(tmp,1), &
ao_cart_to_sphe_inv,size(ao_cart_to_sphe_inv,1), 0.d0, &
ao_sphe_integrals_n_e,size(ao_sphe_integrals_n_e,1))
deallocate(tmp)
END_PROVIDER

357
src/ao_one_e_ints/pot_ao_pseudo_ints.irp.f
View File

@ -8,16 +8,7 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integrals, (ao_num,ao_num)]
call ezfio_get_ao_one_e_ints_ao_integrals_pseudo(ao_pseudo_integrals)
print *, 'AO pseudopotential integrals read from disk'
else
ao_pseudo_integrals = 0.d0
if (do_pseudo) then
if (pseudo_klocmax > 0) then
ao_pseudo_integrals += ao_pseudo_integrals_local
endif
if (pseudo_kmax > 0) then
ao_pseudo_integrals += ao_pseudo_integrals_non_local
endif
endif
call ao_cart_to_ao_basis(ao_cart_pseudo_integrals, ao_cart_num, ao_pseudo_integrals, ao_num)
endif
if (write_ao_integrals_pseudo) then
@ -27,349 +18,3 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integrals, (ao_num,ao_num)]
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_pseudo_integrals_local, (ao_num,ao_num)]
use omp_lib
implicit none
BEGIN_DOC
! Local pseudo-potential
END_DOC
include 'utils/constants.include.F'
double precision :: alpha, beta, gama, delta
integer :: num_A,num_B
double precision :: A_center(3),B_center(3),C_center(3)
integer :: power_A(3),power_B(3)
integer :: i,j,k,l,m
double precision :: Vloc, Vpseudo
double precision :: wall_1, wall_2, wall_0
integer :: thread_num
double precision :: c
double precision :: Z
PROVIDE ao_coef_normalized_ordered_transp
PROVIDE pseudo_v_k_transp pseudo_n_k_transp pseudo_klocmax pseudo_dz_k_transp
ao_pseudo_integrals_local = 0.d0
print*, 'Providing the nuclear electron pseudo integrals (local)'
! Dummy iteration for OpenMP
j=1
i=1
l=1
m=1
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)
power_B(1:3)= ao_power(i,1:3)
B_center(1:3) = nucl_coord(num_B,1:3)
alpha = ao_expo_ordered_transp(l,j)
beta = ao_expo_ordered_transp(m,i)
c = 0.d0
do k = 1, nucl_num
Z = nucl_charge(k)
C_center(1:3) = nucl_coord(k,1:3)
c = c + Vloc(pseudo_klocmax, &
pseudo_v_k_transp (1,k), &
pseudo_n_k_transp (1,k), &
pseudo_dz_k_transp(1,k), &
A_center,power_A,alpha,B_center,power_B,beta,C_center)
enddo
ao_pseudo_integrals_local = 0.d0
call wall_time(wall_1)
thread_num = 0
!$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, &
!$OMP wall_0,wall_2,thread_num) &
!$OMP SHARED (ao_num,ao_prim_num,ao_expo_ordered_transp,ao_power,ao_nucl,nucl_coord,ao_coef_normalized_ordered_transp,&
!$OMP ao_pseudo_integrals_local,nucl_num,nucl_charge, &
!$OMP pseudo_klocmax,pseudo_lmax,pseudo_kmax,pseudo_v_k_transp,pseudo_n_k_transp, pseudo_dz_k_transp,&
!$OMP wall_1)
!$ thread_num = omp_get_thread_num()
wall_0 = wall_1
!$OMP DO
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_ordered_transp(l,j)
do m=1,ao_prim_num(i)
if (dabs(ao_coef_normalized_ordered_transp(l,j)*ao_coef_normalized_ordered_transp(m,i))&
< thresh) then
cycle
endif
beta = ao_expo_ordered_transp(m,i)
c = 0.d0
do k = 1, nucl_num
Z = nucl_charge(k)
C_center(1:3) = nucl_coord(k,1:3)
c = c + Vloc(pseudo_klocmax, &
pseudo_v_k_transp (1,k), &
pseudo_n_k_transp (1,k), &
pseudo_dz_k_transp(1,k), &
A_center,power_A,alpha,B_center,power_B,beta,C_center)
enddo
ao_pseudo_integrals_local(i,j) = ao_pseudo_integrals_local(i,j) +&
ao_coef_normalized_ordered_transp(l,j)*ao_coef_normalized_ordered_transp(m,i)*c
enddo
enddo
enddo
call wall_time(wall_2)
if (thread_num == 0) then
if (wall_2 - wall_0 > 1.d0) then
wall_0 = wall_2
print*, 100.*float(j)/float(ao_num), '% in ', &
wall_2-wall_1, 's'
endif
endif
enddo
!$OMP END DO
!$OMP END PARALLEL
do i=1,ao_num
do j=1,i
ao_pseudo_integrals_local(j,i) = 0.5d0*(ao_pseudo_integrals_local(i,j) + ao_pseudo_integrals_local(i,j))
ao_pseudo_integrals_local(i,j) = ao_pseudo_integrals_local(i,j)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_pseudo_integrals_non_local, (ao_num,ao_num)]
use omp_lib
implicit none
BEGIN_DOC
! Non-local pseudo-potential
END_DOC
include 'utils/constants.include.F'
double precision :: alpha, beta, gama, delta
integer :: num_A,num_B
double precision :: A_center(3),B_center(3),C_center(3)
integer :: power_A(3),power_B(3)
integer :: i,j,k,l,m
double precision :: Vloc, Vpseudo
double precision :: wall_1, wall_2, wall_0
integer :: thread_num
double precision :: c
double precision :: Z
PROVIDE ao_coef_normalized_ordered_transp
PROVIDE pseudo_lmax pseudo_kmax pseudo_v_kl_transp pseudo_n_kl_transp pseudo_dz_kl_transp
ao_pseudo_integrals_non_local = 0.d0
print*, 'Providing the nuclear electron pseudo integrals (non-local)'
call wall_time(wall_1)
thread_num = 0
!$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, &
!$OMP wall_0,wall_2,thread_num) &
!$OMP SHARED (ao_num,ao_prim_num,ao_expo_ordered_transp,ao_power,ao_nucl,nucl_coord,ao_coef_normalized_ordered_transp,&
!$OMP ao_pseudo_integrals_non_local,nucl_num,nucl_charge,&
!$OMP pseudo_klocmax,pseudo_lmax,pseudo_kmax,pseudo_n_kl_transp, pseudo_v_kl_transp, pseudo_dz_kl_transp,&
!$OMP wall_1)
!$ thread_num = omp_get_thread_num()
wall_0 = wall_1
!$OMP DO SCHEDULE (guided)
!
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_ordered_transp(l,j)
do m=1,ao_prim_num(i)
if (dabs(ao_coef_normalized_ordered_transp(l,j)*ao_coef_normalized_ordered_transp(m,i))&
< thresh) then
cycle
endif
beta = ao_expo_ordered_transp(m,i)
c = 0.d0
do k = 1, nucl_num
Z = nucl_charge(k)
C_center(1:3) = nucl_coord(k,1:3)
c = c + Vpseudo(pseudo_lmax,pseudo_kmax, &
pseudo_v_kl_transp(1,0,k), &
pseudo_n_kl_transp(1,0,k), &
pseudo_dz_kl_transp(1,0,k), &
A_center,power_A,alpha,B_center,power_B,beta,C_center)
enddo
ao_pseudo_integrals_non_local(i,j) = ao_pseudo_integrals_non_local(i,j) +&
ao_coef_normalized_ordered_transp(l,j)*ao_coef_normalized_ordered_transp(m,i)*c
enddo
enddo
enddo
call wall_time(wall_2)
if (thread_num == 0) then
if (wall_2 - wall_0 > 1.d0) then
wall_0 = wall_2
print*, 100.*float(j)/float(ao_num), '% in ', &
wall_2-wall_1, 's'
endif
endif
enddo
!$OMP END DO
!$OMP END PARALLEL
do i=1,ao_num
do j=1,i
ao_pseudo_integrals_non_local(j,i) = 0.5d0*(ao_pseudo_integrals_non_local(i,j) + ao_pseudo_integrals_non_local(i,j))
ao_pseudo_integrals_non_local(i,j) = ao_pseudo_integrals_non_local(i,j)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, pseudo_v_k_transp, (pseudo_klocmax,nucl_num) ]
&BEGIN_PROVIDER [ integer , pseudo_n_k_transp, (pseudo_klocmax,nucl_num) ]
&BEGIN_PROVIDER [ double precision, pseudo_dz_k_transp, (pseudo_klocmax,nucl_num)]
implicit none
BEGIN_DOC
! Transposed arrays for pseudopotentials
END_DOC
integer :: i,j
do j=1,nucl_num
do i=1,pseudo_klocmax
pseudo_v_k_transp (i,j) = pseudo_v_k (j,i)
pseudo_n_k_transp (i,j) = pseudo_n_k (j,i)
pseudo_dz_k_transp(i,j) = pseudo_dz_k(j,i)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, pseudo_v_kl_transp, (pseudo_kmax,0:pseudo_lmax,nucl_num) ]
&BEGIN_PROVIDER [ integer , pseudo_n_kl_transp, (pseudo_kmax,0:pseudo_lmax,nucl_num) ]
&BEGIN_PROVIDER [ double precision, pseudo_dz_kl_transp, (pseudo_kmax,0:pseudo_lmax,nucl_num)]
implicit none
BEGIN_DOC
! Transposed arrays for pseudopotentials
END_DOC
integer :: i,j,l
do j=1,nucl_num
do l=0,pseudo_lmax
do i=1,pseudo_kmax
pseudo_v_kl_transp (i,l,j) = pseudo_v_kl (j,i,l)
pseudo_n_kl_transp (i,l,j) = pseudo_n_kl (j,i,l)
pseudo_dz_kl_transp(i,l,j) = pseudo_dz_kl(j,i,l)
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_sphe_pseudo_integrals_local, (ao_sphe_num,ao_sphe_num) ]
implicit none
BEGIN_DOC
! |AO| pseudo_integrals_local matrix in the spherical basis set
END_DOC
double precision, allocatable :: tmp(:,:)
allocate (tmp(ao_sphe_num,ao_num))
call dgemm('T','N',ao_sphe_num,ao_num,ao_num, 1.d0, &
ao_cart_to_sphe_inv,size(ao_cart_to_sphe_inv,1), &
ao_pseudo_integrals_local,size(ao_pseudo_integrals_local,1), 0.d0, &
tmp, size(tmp,1))
call dgemm('N','N',ao_sphe_num,ao_sphe_num,ao_num, 1.d0, &
tmp, size(tmp,1), &
ao_cart_to_sphe_inv,size(ao_cart_to_sphe_inv,1), 0.d0, &
ao_sphe_pseudo_integrals_local,size(ao_sphe_pseudo_integrals_local,1))
deallocate(tmp)
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_sphe_pseudo_integrals_non_local, (ao_sphe_num,ao_sphe_num) ]
implicit none
BEGIN_DOC
! |AO| pseudo_integrals_non_local matrix in the spherical basis set
END_DOC
double precision, allocatable :: tmp(:,:)
allocate (tmp(ao_sphe_num,ao_num))
call dgemm('T','N',ao_sphe_num,ao_num,ao_num, 1.d0, &
ao_cart_to_sphe_inv,size(ao_cart_to_sphe_inv,1), &
ao_pseudo_integrals_non_local,size(ao_pseudo_integrals_non_local,1), 0.d0, &
tmp, size(tmp,1))
call dgemm('N','N',ao_sphe_num,ao_sphe_num,ao_num, 1.d0, &
tmp, size(tmp,1), &
ao_cart_to_sphe_inv,size(ao_cart_to_sphe_inv,1), 0.d0, &
ao_sphe_pseudo_integrals_non_local,size(ao_sphe_pseudo_integrals_non_local,1))
deallocate(tmp)
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_sphe_pseudo_integrals, (ao_sphe_num,ao_sphe_num)]
implicit none
BEGIN_DOC
! Pseudo-potential integrals in the |AO| basis set.
END_DOC
ao_sphe_pseudo_integrals = 0.d0
if (do_pseudo) then
if (pseudo_klocmax > 0) then
ao_sphe_pseudo_integrals += ao_sphe_pseudo_integrals_local
endif
if (pseudo_kmax > 0) then
ao_sphe_pseudo_integrals += ao_sphe_pseudo_integrals_non_local
endif
endif
END_PROVIDER

95
src/ao_one_e_ints/pot_pt_charges.irp.f
View File

@ -10,99 +10,6 @@ BEGIN_PROVIDER [ double precision, ao_integrals_pt_chrg, (ao_num,ao_num)]
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 :: alpha, beta
double precision :: A_center(3),B_center(3),C_center(3)
double precision :: overlap_x,overlap_y,overlap_z,overlap,dx,NAI_pol_mult
ao_integrals_pt_chrg = 0.d0
! if (read_ao_integrals_pt_chrg) then
!
! call ezfio_get_ao_one_e_ints_ao_integrals_pt_chrg(ao_integrals_pt_chrg)
! print *, 'AO N-e integrals read from disk'
!
! else
! if(use_cgtos) then
! !print *, " use_cgtos for ao_integrals_pt_chrg ?", use_cgtos
!
! do j = 1, ao_num
! do i = 1, ao_num
! ao_integrals_pt_chrg(i,j) = ao_integrals_pt_chrg_cgtos(i,j)
! enddo
! enddo
!
! else
!$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,pts_charge_coord,ao_coef_normalized_ordered_transp,nucl_coord,&
!$OMP n_pt_max_integrals,ao_integrals_pt_chrg,n_pts_charge,pts_charge_z)
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_ordered_transp(l,j)
do m=1,ao_prim_num(i)
beta = ao_expo_ordered_transp(m,i)
double precision :: c, c1
c = 0.d0
do k = 1, n_pts_charge
double precision :: Z
Z = pts_charge_z(k)
C_center(1:3) = pts_charge_coord(k,1:3)
c1 = NAI_pol_mult( A_center, B_center, power_A, power_B &
, alpha, beta, C_center, n_pt_in )
c = c - Z * c1
enddo
ao_integrals_pt_chrg(i,j) = ao_integrals_pt_chrg(i,j) &
+ ao_coef_normalized_ordered_transp(l,j) &
* ao_coef_normalized_ordered_transp(m,i) * c
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
! endif
! IF(do_pseudo) THEN
! ao_integrals_pt_chrg += ao_pseudo_integrals
! ENDIF
! endif
! if (write_ao_integrals_pt_chrg) then
! call ezfio_set_ao_one_e_ints_ao_integrals_pt_chrg(ao_integrals_pt_chrg)
! print *, 'AO N-e integrals written to disk'
! endif
call ao_cart_to_ao_basis(ao_cart_integrals_pt_chrg, ao_cart_num, ao_integrals_pt_chrg, ao_num)
END_PROVIDER

2127
src/ao_one_e_ints/pseudopot.f90
View File

File diff suppressed because it is too large Load Diff

2
src/ao_one_e_ints/screening.irp.f
View File

@ -3,7 +3,7 @@ logical function ao_one_e_integral_zero(i,k)
integer, intent(in) :: i,k
ao_one_e_integral_zero = .False.
if (.not.((io_ao_integrals_overlap/='None').or.is_periodic.or.use_cgtos)) then
if((io_ao_integrals_overlap=='None').and.(.not.use_cgtos)) then
if (ao_overlap_abs(i,k) < ao_one_e_integrals_threshold) then
ao_one_e_integral_zero = .True.
return

394
src/ao_one_e_ints/spread_dipole_ao.irp.f
View File

@ -1,376 +1,68 @@
BEGIN_PROVIDER [ double precision, ao_spread_x, (ao_num,ao_num)]
&BEGIN_PROVIDER [ double precision, ao_spread_y, (ao_num,ao_num)]
&BEGIN_PROVIDER [ double precision, ao_spread_z, (ao_num,ao_num)]
BEGIN_PROVIDER [ double precision, ao_spread_x, (ao_num,ao_num)]
implicit none
BEGIN_DOC
! * array of the integrals of AO_i * x^2 AO_j
!
END_DOC
call ao_cart_to_ao_basis(ao_cart_spread_x, ao_cart_num, ao_spread_x, ao_num)
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_spread_y, (ao_num,ao_num)]
implicit none
BEGIN_DOC
! * array of the integrals of AO_i * y^2 AO_j
!
END_DOC
call ao_cart_to_ao_basis(ao_cart_spread_y, ao_cart_num, ao_spread_y, ao_num)
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_spread_z, (ao_num,ao_num)]
implicit none
BEGIN_DOC
! * array of the integrals of AO_i * z^2 AO_j
END_DOC
implicit none
integer :: i,j,n,l
double precision :: f, tmp
integer :: dim1
double precision :: overlap, overlap_x, overlap_y, overlap_z
double precision :: alpha, beta
double precision :: A_center(3), B_center(3)
integer :: power_A(3), power_B(3)
double precision :: lower_exp_val, dx, c,accu_x,accu_y,accu_z
dim1=500
lower_exp_val = 40.d0
ao_spread_x= 0.d0
ao_spread_y= 0.d0
ao_spread_z= 0.d0
!$OMP PARALLEL DO SCHEDULE(GUIDED) &
!$OMP DEFAULT(NONE) &
!$OMP PRIVATE(A_center,B_center,power_A,power_B,&
!$OMP overlap_x,overlap_y, overlap_z, overlap, &
!$OMP alpha, beta,i,j,dx,tmp,c,accu_x,accu_y,accu_z) &
!$OMP SHARED(nucl_coord,ao_power,ao_prim_num, &
!$OMP ao_spread_x,ao_spread_y,ao_spread_z,ao_num,ao_coef_normalized_ordered_transp,ao_nucl, &
!$OMP ao_expo_ordered_transp,dim1,lower_exp_val)
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
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 )
accu_x = 0.d0
accu_y = 0.d0
accu_z = 0.d0
do n = 1,ao_prim_num(j)
alpha = ao_expo_ordered_transp(n,j)
do l = 1, ao_prim_num(i)
c = ao_coef_normalized_ordered_transp(n,j)*ao_coef_normalized_ordered_transp(l,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)
call overlap_bourrin_spread(A_center(1),B_center(1),alpha,beta,power_A(1),power_B(1),tmp,lower_exp_val,dx,dim1)
accu_x += c*tmp*overlap_y*overlap_z
call overlap_bourrin_spread(A_center(2),B_center(2),alpha,beta,power_A(2),power_B(2),tmp,lower_exp_val,dx,dim1)
accu_y += c*tmp*overlap_x*overlap_z
call overlap_bourrin_spread(A_center(3),B_center(3),alpha,beta,power_A(3),power_B(3),tmp,lower_exp_val,dx,dim1)
accu_z += c*tmp*overlap_y*overlap_x
enddo
enddo
ao_spread_x(i,j) = accu_x
ao_spread_y(i,j) = accu_y
ao_spread_z(i,j) = accu_z
enddo
enddo
!$OMP END PARALLEL DO
call ao_cart_to_ao_basis(ao_cart_spread_z, ao_cart_num, ao_spread_z, ao_num)
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_dipole_x, (ao_num,ao_num)]
&BEGIN_PROVIDER [ double precision, ao_dipole_y, (ao_num,ao_num)]
&BEGIN_PROVIDER [ double precision, ao_dipole_z, (ao_num,ao_num)]
BEGIN_PROVIDER [ double precision, ao_dipole_x, (ao_num,ao_num)]
BEGIN_DOC
! * array of the integrals of AO_i * x AO_j
!
END_DOC
implicit none
call ao_cart_to_ao_basis(ao_cart_dipole_x, ao_cart_num, ao_dipole_x, ao_num)
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_dipole_y, (ao_num,ao_num)]
BEGIN_DOC
! * array of the integrals of AO_i * y AO_j
!
END_DOC
implicit none
call ao_cart_to_ao_basis(ao_cart_dipole_y, ao_cart_num, ao_dipole_y, ao_num)
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_dipole_z, (ao_num,ao_num)]
BEGIN_DOC
! * array of the integrals of AO_i * z AO_j
END_DOC
implicit none
integer :: i,j,n,l
double precision :: f, tmp
integer :: dim1
double precision :: overlap, overlap_x, overlap_y, overlap_z,accu_x,accu_y,accu_z
double precision :: alpha, beta
double precision :: A_center(3), B_center(3)
integer :: power_A(3), power_B(3)
double precision :: lower_exp_val, dx, c
dim1=500
lower_exp_val = 40.d0
ao_dipole_x= 0.d0
ao_dipole_y= 0.d0
ao_dipole_z= 0.d0
!$OMP PARALLEL DO SCHEDULE(GUIDED) &
!$OMP DEFAULT(NONE) &
!$OMP PRIVATE(A_center,B_center,power_A,power_B,&
!$OMP overlap_x,overlap_y, overlap_z, overlap, &
!$OMP alpha, beta,i,j,dx,tmp,c,accu_x,accu_y,accu_z) &
!$OMP SHARED(nucl_coord,ao_power,ao_prim_num, &
!$OMP ao_dipole_x,ao_dipole_y,ao_dipole_z,ao_num,ao_coef_normalized_ordered_transp,ao_nucl, &
!$OMP ao_expo_ordered_transp,dim1,lower_exp_val)
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
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 )
accu_x = 0.d0
accu_y = 0.d0
accu_z = 0.d0
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)
c = ao_coef_normalized_ordered_transp(l,i)*ao_coef_normalized_ordered_transp(n,j)
call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,overlap_y,overlap_z,overlap,dim1)
call overlap_bourrin_dipole(A_center(1),B_center(1),alpha,beta,power_A(1),power_B(1),tmp,lower_exp_val,dx,dim1)
accu_x = accu_x + c*tmp*overlap_y*overlap_z
call overlap_bourrin_dipole(A_center(2),B_center(2),alpha,beta,power_A(2),power_B(2),tmp,lower_exp_val,dx,dim1)
accu_y = accu_y + c*tmp*overlap_x*overlap_z
call overlap_bourrin_dipole(A_center(3),B_center(3),alpha,beta,power_A(3),power_B(3),tmp,lower_exp_val,dx,dim1)
accu_z = accu_z + c*tmp*overlap_y*overlap_x
enddo
enddo
ao_dipole_x(i,j) = accu_x
ao_dipole_y(i,j) = accu_y
ao_dipole_z(i,j) = accu_z
enddo
enddo
!$OMP END PARALLEL DO
call ao_cart_to_ao_basis(ao_cart_dipole_z, ao_cart_num, ao_dipole_z, ao_num)
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_deriv_1_x, (ao_num,ao_num)]
&BEGIN_PROVIDER [ double precision, ao_deriv_1_y, (ao_num,ao_num)]
&BEGIN_PROVIDER [ double precision, ao_deriv_1_z, (ao_num,ao_num)]
BEGIN_DOC
! * array of the integrals of AO_i * d/dx AO_j
!
END_DOC
implicit none
call ao_cart_to_ao_basis(ao_cart_deriv_1_x, ao_cart_num, ao_deriv_1_x, ao_num)
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_deriv_1_y, (ao_num,ao_num)]
BEGIN_DOC
! * array of the integrals of AO_i * d/dy AO_j
!
END_DOC
implicit none
call ao_cart_to_ao_basis(ao_cart_deriv_1_y, ao_cart_num, ao_deriv_1_y, ao_num)
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_deriv_1_z, (ao_num,ao_num)]
BEGIN_DOC
! * array of the integrals of AO_i * d/dz AO_j
END_DOC
implicit none
integer :: i,j,n,l
double precision :: f, tmp
integer :: dim1
double precision :: overlap, overlap_x, overlap_y, overlap_z
double precision :: alpha, beta
double precision :: A_center(3), B_center(3)
integer :: power_A(3), power_B(3)
double precision :: lower_exp_val, dx, c,accu_x,accu_y,accu_z
integer :: i_component
dim1=500
lower_exp_val = 40.d0
ao_deriv_1_x= 0.d0
ao_deriv_1_y= 0.d0
ao_deriv_1_z= 0.d0
!$OMP PARALLEL DO SCHEDULE(GUIDED) &
!$OMP DEFAULT(NONE) &
!$OMP PRIVATE(A_center,B_center,power_A,power_B,&
!$OMP overlap_x,overlap_y, overlap_z, overlap, &
!$OMP alpha, beta,i,j,dx,tmp,c,i_component,accu_x,accu_y,accu_z) &
!$OMP SHARED(nucl_coord,ao_power,ao_prim_num, &
!$OMP ao_deriv_1_x,ao_deriv_1_y,ao_deriv_1_z,ao_num,ao_coef_normalized_ordered_transp,ao_nucl, &
!$OMP ao_expo_ordered_transp,dim1,lower_exp_val)
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
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 )
accu_x = 0.d0
accu_y = 0.d0
accu_z = 0.d0
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_x,overlap_y,overlap_z,overlap,dim1)
c = ao_coef_normalized_ordered_transp(l,i) * ao_coef_normalized_ordered_transp(n,j)
i_component = 1
call overlap_bourrin_deriv_x(i_component,A_center,B_center,alpha,beta,power_A,power_B,dx,lower_exp_val,tmp,dim1)
accu_x += c*(tmp*overlap_y*overlap_z)
i_component = 2
call overlap_bourrin_deriv_x(i_component,A_center,B_center,alpha,beta,power_A,power_B,dx,lower_exp_val,tmp,dim1)
accu_y += c*(tmp*overlap_x*overlap_z)
i_component = 3
call overlap_bourrin_deriv_x(i_component,A_center,B_center,alpha,beta,power_A,power_B,dx,lower_exp_val,tmp,dim1)
accu_z += c*(tmp*overlap_y*overlap_x)
enddo
enddo
ao_deriv_1_x(i,j) = accu_x
ao_deriv_1_y(i,j) = accu_y
ao_deriv_1_z(i,j) = accu_z
enddo
enddo
!$OMP END PARALLEL DO
call ao_cart_to_ao_basis(ao_cart_deriv_1_z, ao_cart_num, ao_deriv_1_z, ao_num)
END_PROVIDER
subroutine overlap_bourrin_spread(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,lower_exp_val,dx,nx)
BEGIN_DOC
! Computes the following integral :
! int [-infty ; +infty] of [(x-A_center)^(power_A) * (x-B_center)^power_B * exp(-alpha(x-A_center)^2) * exp(-beta(x-B_center)^2) * x^2 ]
! needed for the dipole and those things
END_DOC
implicit none
integer :: i,j,k,l
integer,intent(in) :: power_A,power_B
double precision, intent(in) :: lower_exp_val
double precision,intent(in) :: A_center, B_center,alpha,beta
double precision, intent(out) :: overlap_x,dx
integer, intent(in) :: nx
double precision :: x_min,x_max,domain,x,factor,dist,p,p_inv,rho
double precision :: P_center,pouet_timy
if(power_A.lt.0.or.power_B.lt.0)then
overlap_x = 0.d0
dx = 0.d0
return
endif
p = alpha + beta
p_inv= 1.d0/p
rho = alpha * beta * p_inv
dist = (A_center - B_center)*(A_center - B_center)
P_center = (alpha * A_center + beta * B_center) * p_inv
factor = dexp(-rho * dist)
if(factor.lt.0.000001d0)then
! print*,'factor = ',factor
dx = 0.d0
overlap_x = 0.d0
return
endif
pouet_timy = dsqrt(lower_exp_val/p)
x_min = P_center - pouet_timy
x_max = P_center + pouet_timy
domain = x_max-x_min
dx = domain/dble(nx)
overlap_x = 0.d0
x = x_min
do i = 1, nx
x += dx
overlap_x += (x-A_center)**(power_A) * (x-B_center)**(power_B) * dexp(-p * (x-P_center)*(x-P_center)) * x * x
enddo
overlap_x *= factor * dx
end
subroutine overlap_bourrin_dipole(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,lower_exp_val,dx,nx)
! compute the following integral :
! int [-infty ; +infty] of [(x-A_center)^(power_A) * (x-B_center)^power_B * exp(-alpha(x-A_center)^2) * exp(-beta(x-B_center)^2) * x ]
! needed for the dipole and those things
implicit none
integer :: i,j,k,l
integer,intent(in) :: power_A,power_B
double precision, intent(in) :: lower_exp_val
double precision,intent(in) :: A_center, B_center,alpha,beta
double precision, intent(out) :: overlap_x,dx
integer, intent(in) :: nx
double precision :: x_min,x_max,domain,x,factor,dist,p,p_inv,rho
double precision :: P_center
if(power_A.lt.0.or.power_B.lt.0)then
overlap_x = 0.d0
dx = 0.d0
return
endif
p = alpha + beta
p_inv= 1.d0/p
rho = alpha * beta * p_inv
dist = (A_center - B_center)*(A_center - B_center)
P_center = (alpha * A_center + beta * B_center) * p_inv
factor = dexp(-rho * dist)
if(power_B == 0 .and. power_A ==0)then
double precision :: F_integral
overlap_x = P_center * F_integral(0,p) * factor
dx = 0.d0
return
endif
double precision :: pouet_timy
pouet_timy = dsqrt(lower_exp_val/p)
x_min = P_center - pouet_timy
x_max = P_center + pouet_timy
domain = x_max-x_min
dx = domain/dble(nx)
overlap_x = 0.d0
x = x_min
do i = 1, nx
x += dx
overlap_x += (x-A_center)**(power_A) * (x-B_center)**(power_B) * dexp(-p * (x-P_center)*(x-P_center)) * x
enddo
overlap_x *= factor * dx
end
subroutine overlap_bourrin_deriv_x(i_component,A_center,B_center,alpha,beta,power_A,power_B,dx,lower_exp_val,overlap_x,nx)
implicit none
integer :: i,j,k,l
integer,intent(in) :: power_A(3),power_B(3),i_component
double precision,intent(in) :: A_center(3), B_center(3),alpha,beta,lower_exp_val
double precision, intent(out) :: overlap_x,dx
integer, intent(in) :: nx
double precision :: overlap_first, overlap_second
! computes : <phi_i|d/dx|phi_j> = (a_x_i <phi_i_x|phi_j_x(a_x_j-1)> - 2 alpha <phi_i_x|phi_j_w(a_x_j+1)>)
call overlap_bourrin_x(A_center(i_component),B_center(i_component),alpha,beta,power_A(i_component)-1,power_B(i_component),overlap_first,lower_exp_val,dx,nx)
call overlap_bourrin_x(A_center(i_component),B_center(i_component),alpha,beta,power_A(i_component)+1,power_B(i_component),overlap_second,lower_exp_val,dx,nx)
overlap_x = (power_A(i_component) * overlap_first - 2.d0 * alpha * overlap_second)
end
subroutine overlap_bourrin_x(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,lower_exp_val,dx,nx)
implicit none
! compute the following integral :
! int [-infty ; +infty] of [(x-A_center)^(power_A) * (x-B_center)^power_B * exp(-alpha(x-A_center)^2) * exp(-beta(x-B_center)^2) ]
integer :: i,j,k,l
integer,intent(in) :: power_A,power_B
double precision, intent(in) :: lower_exp_val
double precision,intent(in) :: A_center, B_center,alpha,beta
double precision, intent(out) :: overlap_x,dx
integer, intent(in) :: nx
double precision :: x_min,x_max,domain,x,factor,dist,p,p_inv,rho
double precision :: P_center,pouet_timy
if(power_A.lt.0.or.power_B.lt.0)then
overlap_x = 0.d0
dx = 0.d0
return
endif
p = alpha + beta
p_inv= 1.d0/p
rho = alpha * beta * p_inv
dist = (A_center - B_center)*(A_center - B_center)
P_center = (alpha * A_center + beta * B_center) * p_inv
factor = dexp(-rho * dist)
if(factor.lt.0.000001d0)then
dx = 0.d0
overlap_x = 0.d0
return
endif
pouet_timy = dsqrt(lower_exp_val/p)
x_min = P_center - pouet_timy
x_max = P_center + pouet_timy
domain = x_max-x_min
dx = domain/dble(nx)
overlap_x = 0.d0
x = x_min
do i = 1, nx
x += dx
overlap_x += (x-A_center)**(power_A) * (x-B_center)**(power_B) * dexp(-p * (x-P_center)*(x-P_center))
enddo
overlap_x *= factor * dx
end

27
src/ao_one_e_ints/test_ao_one_ints.irp.f Normal file
View File

@ -0,0 +1,27 @@
BEGIN_PROVIDER [double precision, ao_ortho_canonical_overlap, (ao_ortho_canonical_num,ao_ortho_canonical_num)]
implicit none
BEGIN_DOC
! overlap matrix of the ao_ortho_canonical.
! Expected to be the Identity
END_DOC
integer :: i,j,k,l
double precision :: c
do j=1, ao_ortho_canonical_num
do i=1, ao_ortho_canonical_num
ao_ortho_canonical_overlap(i,j) = 0.d0
enddo
enddo
do j=1, ao_ortho_canonical_num
do k=1, ao_num
c = 0.d0
do l=1, ao_num
c += ao_ortho_canonical_coef(l,j) * ao_overlap(l,k)
enddo
do i=1, ao_ortho_canonical_num
ao_ortho_canonical_overlap(i,j) += ao_ortho_canonical_coef(k,i) * c
enddo
enddo
enddo
END_PROVIDER