mirror of
https://github.com/QuantumPackage/qp2.git
synced 2024-09-02 01:03:45 +02:00
Merge pull request #3 from AbdAmmar/dev-stable-cosgtos
cos x GTOs integ added
This commit is contained in:
commit
b61125cf3a
@ -12,21 +12,21 @@ double precision function ao_value(i,r)
|
||||
integer :: power_ao(3)
|
||||
double precision :: accu,dx,dy,dz,r2
|
||||
num_ao = ao_nucl(i)
|
||||
! power_ao(1:3)= ao_power(i,1:3)
|
||||
! center_ao(1:3) = nucl_coord(num_ao,1:3)
|
||||
! dx = (r(1) - center_ao(1))
|
||||
! dy = (r(2) - center_ao(2))
|
||||
! dz = (r(3) - center_ao(3))
|
||||
! r2 = dx*dx + dy*dy + dz*dz
|
||||
! dx = dx**power_ao(1)
|
||||
! dy = dy**power_ao(2)
|
||||
! dz = dz**power_ao(3)
|
||||
power_ao(1:3)= ao_power(i,1:3)
|
||||
center_ao(1:3) = nucl_coord(num_ao,1:3)
|
||||
dx = (r(1) - center_ao(1))
|
||||
dy = (r(2) - center_ao(2))
|
||||
dz = (r(3) - center_ao(3))
|
||||
r2 = dx*dx + dy*dy + dz*dz
|
||||
dx = dx**power_ao(1)
|
||||
dy = dy**power_ao(2)
|
||||
dz = dz**power_ao(3)
|
||||
|
||||
accu = 0.d0
|
||||
! do m=1,ao_prim_num(i)
|
||||
! beta = ao_expo_ordered_transp(m,i)
|
||||
! accu += ao_coef_normalized_ordered_transp(m,i) * dexp(-beta*r2)
|
||||
! enddo
|
||||
do m=1,ao_prim_num(i)
|
||||
beta = ao_expo_ordered_transp(m,i)
|
||||
accu += ao_coef_normalized_ordered_transp(m,i) * dexp(-beta*r2)
|
||||
enddo
|
||||
ao_value = accu * dx * dy * dz
|
||||
|
||||
end
|
||||
|
@ -1,2 +1,3 @@
|
||||
ao_basis
|
||||
pseudo
|
||||
cosgtos_ao_int
|
||||
|
@ -1,75 +1,99 @@
|
||||
BEGIN_PROVIDER [ double precision, ao_overlap,(ao_num,ao_num) ]
|
||||
&BEGIN_PROVIDER [ double precision, ao_overlap_x,(ao_num,ao_num) ]
|
||||
&BEGIN_PROVIDER [ double precision, ao_overlap_y,(ao_num,ao_num) ]
|
||||
&BEGIN_PROVIDER [ double precision, ao_overlap_z,(ao_num,ao_num) ]
|
||||
implicit none
|
||||
|
||||
! ---
|
||||
|
||||
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:
|
||||
!
|
||||
! :math:`\int \chi_i(r) \chi_j(r) dr`
|
||||
! Overlap between atomic basis functions:
|
||||
!
|
||||
! :math:`\int \chi_i(r) \chi_j(r) dr`
|
||||
END_DOC
|
||||
integer :: i,j,n,l
|
||||
double precision :: f
|
||||
integer :: dim1
|
||||
|
||||
implicit none
|
||||
integer :: i, j, n, l, dim1, power_A(3), power_B(3)
|
||||
double precision :: overlap, overlap_x, overlap_y, overlap_z
|
||||
double precision :: alpha, beta, c
|
||||
double precision :: A_center(3), B_center(3)
|
||||
integer :: power_A(3), power_B(3)
|
||||
ao_overlap = 0.d0
|
||||
|
||||
ao_overlap = 0.d0
|
||||
ao_overlap_x = 0.d0
|
||||
ao_overlap_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'
|
||||
|
||||
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'
|
||||
|
||||
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,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
|
||||
if(use_cosgtos) then
|
||||
!print*, ' use_cosgtos for ao_overlap ?', use_cosgtos
|
||||
|
||||
do j = 1, ao_num
|
||||
do i = 1, ao_num
|
||||
ao_overlap (i,j) = ao_overlap_cosgtos (i,j)
|
||||
ao_overlap_x(i,j) = ao_overlap_cosgtos_x(i,j)
|
||||
ao_overlap_y(i,j) = ao_overlap_cosgtos_y(i,j)
|
||||
ao_overlap_z(i,j) = ao_overlap_cosgtos_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,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
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END PARALLEL DO
|
||||
!$OMP END PARALLEL DO
|
||||
|
||||
endif
|
||||
|
||||
endif
|
||||
|
||||
if (write_ao_integrals_overlap) then
|
||||
call ezfio_set_ao_one_e_ints_ao_integrals_overlap(ao_overlap(1:ao_num, 1:ao_num))
|
||||
print *, 'AO overlap integrals written to disk'
|
||||
@ -77,6 +101,8 @@
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [ double precision, ao_overlap_imag, (ao_num, ao_num) ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
@ -85,6 +111,8 @@ BEGIN_PROVIDER [ double precision, ao_overlap_imag, (ao_num, ao_num) ]
|
||||
ao_overlap_imag = 0.d0
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [ complex*16, ao_overlap_complex, (ao_num, ao_num) ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
@ -98,41 +126,43 @@ BEGIN_PROVIDER [ complex*16, ao_overlap_complex, (ao_num, ao_num) ]
|
||||
enddo
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [ double precision, ao_overlap_abs, (ao_num, ao_num) ]
|
||||
|
||||
|
||||
BEGIN_PROVIDER [ double precision, ao_overlap_abs,(ao_num,ao_num) ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Overlap between absolute values of atomic basis functions:
|
||||
!
|
||||
! :math:`\int |\chi_i(r)| |\chi_j(r)| dr`
|
||||
! Overlap between absolute values of atomic basis functions:
|
||||
!
|
||||
! :math:`\int |\chi_i(r)| |\chi_j(r)| dr`
|
||||
END_DOC
|
||||
integer :: i,j,n,l
|
||||
double precision :: f
|
||||
integer :: dim1
|
||||
double precision :: overlap, overlap_x, overlap_y, overlap_z
|
||||
|
||||
implicit none
|
||||
integer :: i, j, n, l, dim1, power_A(3), power_B(3)
|
||||
double precision :: 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
|
||||
if (is_periodic) then
|
||||
do j=1,ao_num
|
||||
do i= 1,ao_num
|
||||
ao_overlap_abs(i,j)= cdabs(ao_overlap_complex(i,j))
|
||||
|
||||
if(is_periodic) then
|
||||
|
||||
do j = 1, ao_num
|
||||
do i = 1, ao_num
|
||||
ao_overlap_abs(i,j) = cdabs(ao_overlap_complex(i,j))
|
||||
enddo
|
||||
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, overlap, &
|
||||
!$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)
|
||||
!$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 )
|
||||
@ -160,10 +190,14 @@ BEGIN_PROVIDER [ double precision, ao_overlap_abs,(ao_num,ao_num) ]
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END PARALLEL DO
|
||||
!$OMP END PARALLEL DO
|
||||
|
||||
endif
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [ double precision, S_inv,(ao_num,ao_num) ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
|
@ -1,7 +1,10 @@
|
||||
BEGIN_PROVIDER [ double precision, ao_deriv2_x,(ao_num,ao_num) ]
|
||||
&BEGIN_PROVIDER [ double precision, ao_deriv2_y,(ao_num,ao_num) ]
|
||||
&BEGIN_PROVIDER [ double precision, ao_deriv2_z,(ao_num,ao_num) ]
|
||||
implicit none
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [ double precision, ao_deriv2_x, (ao_num, ao_num) ]
|
||||
&BEGIN_PROVIDER [ double precision, ao_deriv2_y, (ao_num, ao_num) ]
|
||||
&BEGIN_PROVIDER [ double precision, ao_deriv2_z, (ao_num, ao_num) ]
|
||||
|
||||
BEGIN_DOC
|
||||
! Second derivative matrix elements in the |AO| basis.
|
||||
!
|
||||
@ -11,114 +14,131 @@
|
||||
! \langle \chi_i(x,y,z) | \frac{\partial^2}{\partial x^2} |\chi_j (x,y,z) \rangle
|
||||
!
|
||||
END_DOC
|
||||
integer :: i,j,n,l
|
||||
double precision :: f
|
||||
integer :: dim1
|
||||
|
||||
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)
|
||||
integer :: power_A(3), power_B(3)
|
||||
double precision :: d_a_2,d_2
|
||||
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)
|
||||
! --
|
||||
if(use_cosgtos) then
|
||||
!print*, 'use_cosgtos for ao_kinetic_integrals ?', use_cosgtos
|
||||
|
||||
!$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,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)
|
||||
do j = 1, ao_num
|
||||
do i = 1, ao_num
|
||||
ao_deriv2_x(i,j) = ao_deriv2_cosgtos_x(i,j)
|
||||
ao_deriv2_y(i,j) = ao_deriv2_cosgtos_y(i,j)
|
||||
ao_deriv2_z(i,j) = ao_deriv2_cosgtos_z(i,j)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
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
|
||||
else
|
||||
|
||||
double precision :: deriv_tmp
|
||||
deriv_tmp = (-2.d0 * alpha * (2.d0 * power_A(1) +1.d0) * overlap_x0 &
|
||||
+power_A(1) * (power_A(1)-1.d0) * d_a_2 &
|
||||
+4.d0 * alpha * alpha * d_2 )*overlap_y0*overlap_z0
|
||||
dim1=100
|
||||
|
||||
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
|
||||
! -- 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)
|
||||
! --
|
||||
|
||||
deriv_tmp = (-2.d0 * alpha * (2.d0 * power_A(2) +1.d0 ) * overlap_y0 &
|
||||
+power_A(2) * (power_A(2)-1.d0) * d_a_2 &
|
||||
+4.d0 * alpha * alpha * d_2 )*overlap_x0*overlap_z0
|
||||
ao_deriv2_y(i,j) += c*deriv_tmp
|
||||
!$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,i,j,c,d_a_2,d_2,deriv_tmp, &
|
||||
!$OMP overlap_x0,overlap_y0,overlap_z0) &
|
||||
!$OMP SHARED(nucl_coord,ao_power,ao_prim_num, &
|
||||
!$OMP ao_deriv2_x,ao_deriv2_y,ao_deriv2_z,ao_num,ao_coef_normalized_ordered_transp,ao_nucl, &
|
||||
!$OMP ao_expo_ordered_transp,dim1)
|
||||
do j=1,ao_num
|
||||
A_center(1) = nucl_coord( ao_nucl(j), 1 )
|
||||
A_center(2) = nucl_coord( ao_nucl(j), 2 )
|
||||
A_center(3) = nucl_coord( ao_nucl(j), 3 )
|
||||
power_A(1) = ao_power( j, 1 )
|
||||
power_A(2) = ao_power( j, 2 )
|
||||
power_A(3) = ao_power( j, 3 )
|
||||
do i= 1,ao_num
|
||||
ao_deriv2_x(i,j)= 0.d0
|
||||
ao_deriv2_y(i,j)= 0.d0
|
||||
ao_deriv2_z(i,j)= 0.d0
|
||||
B_center(1) = nucl_coord( ao_nucl(i), 1 )
|
||||
B_center(2) = nucl_coord( ao_nucl(i), 2 )
|
||||
B_center(3) = nucl_coord( ao_nucl(i), 3 )
|
||||
power_B(1) = ao_power( i, 1 )
|
||||
power_B(2) = ao_power( i, 2 )
|
||||
power_B(3) = ao_power( i, 3 )
|
||||
do n = 1,ao_prim_num(j)
|
||||
alpha = ao_expo_ordered_transp(n,j)
|
||||
do l = 1, ao_prim_num(i)
|
||||
beta = ao_expo_ordered_transp(l,i)
|
||||
call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x0,overlap_y0,overlap_z0,overlap,dim1)
|
||||
c = ao_coef_normalized_ordered_transp(n,j) * ao_coef_normalized_ordered_transp(l,i)
|
||||
|
||||
power_A(3) = power_A(3)-2
|
||||
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
|
||||
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
|
||||
|
||||
deriv_tmp = (-2.d0 * alpha * (2.d0 * power_A(3) +1.d0 ) * overlap_z0 &
|
||||
+power_A(3) * (power_A(3)-1.d0) * d_a_2 &
|
||||
+4.d0 * alpha * alpha * d_2 )*overlap_x0*overlap_y0
|
||||
ao_deriv2_z(i,j) += c*deriv_tmp
|
||||
double precision :: deriv_tmp
|
||||
deriv_tmp = (-2.d0 * alpha * (2.d0 * power_A(1) +1.d0) * overlap_x0 &
|
||||
+power_A(1) * (power_A(1)-1.d0) * d_a_2 &
|
||||
+4.d0 * alpha * alpha * d_2 )*overlap_y0*overlap_z0
|
||||
|
||||
ao_deriv2_x(i,j) += c*deriv_tmp
|
||||
power_A(2) = power_A(2)-2
|
||||
if (power_A(2)>-1) then
|
||||
call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_y,d_a_2,overlap_z,overlap,dim1)
|
||||
else
|
||||
d_a_2 = 0.d0
|
||||
endif
|
||||
power_A(2) = power_A(2)+4
|
||||
call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_y,d_2,overlap_z,overlap,dim1)
|
||||
power_A(2) = power_A(2)-2
|
||||
|
||||
deriv_tmp = (-2.d0 * alpha * (2.d0 * power_A(2) +1.d0 ) * overlap_y0 &
|
||||
+power_A(2) * (power_A(2)-1.d0) * d_a_2 &
|
||||
+4.d0 * alpha * alpha * d_2 )*overlap_x0*overlap_z0
|
||||
ao_deriv2_y(i,j) += c*deriv_tmp
|
||||
|
||||
power_A(3) = power_A(3)-2
|
||||
if (power_A(3)>-1) then
|
||||
call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_y,overlap_z,d_a_2,overlap,dim1)
|
||||
else
|
||||
d_a_2 = 0.d0
|
||||
endif
|
||||
power_A(3) = power_A(3)+4
|
||||
call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_y,overlap_z,d_2,overlap,dim1)
|
||||
power_A(3) = power_A(3)-2
|
||||
|
||||
deriv_tmp = (-2.d0 * alpha * (2.d0 * power_A(3) +1.d0 ) * overlap_z0 &
|
||||
+power_A(3) * (power_A(3)-1.d0) * d_a_2 &
|
||||
+4.d0 * alpha * alpha * d_2 )*overlap_x0*overlap_y0
|
||||
ao_deriv2_z(i,j) += c*deriv_tmp
|
||||
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END PARALLEL DO
|
||||
!$OMP END PARALLEL DO
|
||||
|
||||
endif
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [double precision, ao_kinetic_integrals, (ao_num,ao_num)]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
|
@ -1,3 +1,6 @@
|
||||
|
||||
! ---
|
||||
|
||||
subroutine give_all_erf_kl_ao(integrals_ao,mu_in,C_center)
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
@ -15,36 +18,104 @@ subroutine give_all_erf_kl_ao(integrals_ao,mu_in,C_center)
|
||||
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(i_ao,j_ao,mu_in,C_center)
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
!
|
||||
! Computes the following integral :
|
||||
! $\int_{-\infty}^{infty} dr \chi_i(r) \chi_j(r) \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$.
|
||||
! $\int_{-\infty}^{infty} dr \chi_i(r) \chi_j(r) \frac{\erf(\mu |r - R_C|)}{|r - R_C|}$.
|
||||
!
|
||||
END_DOC
|
||||
integer, intent(in) :: i_ao,j_ao
|
||||
|
||||
implicit none
|
||||
integer, intent(in) :: i_ao, j_ao
|
||||
double precision, intent(in) :: mu_in, C_center(3)
|
||||
integer :: i,j,num_A,num_B, power_A(3), power_B(3), n_pt_in
|
||||
double precision :: A_center(3), B_center(3),integral, alpha,beta
|
||||
|
||||
integer :: i, j, num_A, num_B, power_A(3), power_B(3), n_pt_in
|
||||
double precision :: A_center(3), B_center(3), integral, alpha, beta
|
||||
|
||||
double precision :: NAI_pol_mult_erf
|
||||
num_A = ao_nucl(i_ao)
|
||||
power_A(1:3)= ao_power(i_ao,1:3)
|
||||
|
||||
num_A = ao_nucl(i_ao)
|
||||
power_A(1:3) = ao_power(i_ao,1:3)
|
||||
A_center(1:3) = nucl_coord(num_A,1:3)
|
||||
num_B = ao_nucl(j_ao)
|
||||
power_B(1:3)= ao_power(j_ao,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)
|
||||
|
||||
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
|
||||
|
||||
end
|
||||
end function NAI_pol_mult_erf_ao
|
||||
|
||||
! ---
|
||||
|
||||
double precision function NAI_pol_mult_erf_ao_with1s(i_ao, j_ao, beta, B_center, mu_in, C_center)
|
||||
|
||||
BEGIN_DOC
|
||||
!
|
||||
! Computes the following integral :
|
||||
! $\int_{-\infty}^{infty} dr \chi_i(r) \chi_j(r) e^{-\beta (r - B_center)^2} \frac{\erf(\mu |r - R_C|)}{|r - R_C|}$.
|
||||
!
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
integer, intent(in) :: i_ao, j_ao
|
||||
double precision, intent(in) :: beta, B_center(3)
|
||||
double precision, intent(in) :: mu_in, C_center(3)
|
||||
|
||||
integer :: i, j, power_A1(3), power_A2(3), n_pt_in
|
||||
double precision :: A1_center(3), A2_center(3), alpha1, alpha2, coef12, coef1, integral
|
||||
|
||||
double precision, external :: NAI_pol_mult_erf_with1s, NAI_pol_mult_erf_ao
|
||||
|
||||
ASSERT(beta .ge. 0.d0)
|
||||
if(beta .lt. 1d-10) then
|
||||
NAI_pol_mult_erf_ao_with1s = NAI_pol_mult_erf_ao(i_ao, j_ao, mu_in, C_center)
|
||||
return
|
||||
endif
|
||||
|
||||
power_A1(1:3) = ao_power(i_ao,1:3)
|
||||
power_A2(1:3) = ao_power(j_ao,1:3)
|
||||
|
||||
A1_center(1:3) = nucl_coord(ao_nucl(i_ao),1:3)
|
||||
A2_center(1:3) = nucl_coord(ao_nucl(j_ao),1:3)
|
||||
|
||||
n_pt_in = n_pt_max_integrals
|
||||
|
||||
NAI_pol_mult_erf_ao_with1s = 0.d0
|
||||
do i = 1, ao_prim_num(i_ao)
|
||||
alpha1 = ao_expo_ordered_transp (i,i_ao)
|
||||
coef1 = ao_coef_normalized_ordered_transp(i,i_ao)
|
||||
|
||||
do j = 1, ao_prim_num(j_ao)
|
||||
alpha2 = ao_expo_ordered_transp(j,j_ao)
|
||||
coef12 = coef1 * ao_coef_normalized_ordered_transp(j,j_ao)
|
||||
if(dabs(coef12) .lt. 1d-14) cycle
|
||||
|
||||
integral = NAI_pol_mult_erf_with1s( A1_center, A2_center, power_A1, power_A2, alpha1, alpha2 &
|
||||
, beta, B_center, C_center, n_pt_in, mu_in )
|
||||
|
||||
NAI_pol_mult_erf_ao_with1s += integral * coef12
|
||||
enddo
|
||||
enddo
|
||||
|
||||
end function NAI_pol_mult_erf_ao_with1s
|
||||
|
||||
! ---
|
||||
|
||||
double precision function NAI_pol_mult_erf(A_center, B_center, power_A, power_B, alpha, beta, C_center, n_pt_in, mu_in)
|
||||
|
||||
@ -127,58 +198,221 @@ end function NAI_pol_mult_erf
|
||||
|
||||
! ---
|
||||
|
||||
|
||||
double precision function NAI_pol_mult_erf_ao_with1s(i_ao, j_ao, beta, B_center, mu_in, C_center)
|
||||
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 :
|
||||
! $\int_{-\infty}^{infty} dr \chi_i(r) \chi_j(r) e^{-\beta (r - B_center)^2} \frac{\erf(\mu |r - R_C|)}{|r - R_C|}$.
|
||||
!
|
||||
! .. math::
|
||||
!
|
||||
! \int dr (x-A_x)^a (x-B_x)^b \exp(-\alpha (x-A_x)^2 - \beta (x-B_x)^2 )
|
||||
! \frac{\erf(\mu |r - R_C |)}{| r - R_C |}$.
|
||||
!
|
||||
END_DOC
|
||||
|
||||
include 'utils/constants.include.F'
|
||||
|
||||
implicit none
|
||||
integer, intent(in) :: i_ao, j_ao
|
||||
double precision, intent(in) :: beta, B_center(3)
|
||||
double precision, intent(in) :: mu_in, C_center(3)
|
||||
|
||||
integer :: i, j, power_A1(3), power_A2(3), n_pt_in
|
||||
double precision :: A1_center(3), A2_center(3), alpha1, alpha2, coef12, coef1, integral
|
||||
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)
|
||||
|
||||
double precision, external :: NAI_pol_mult_erf_with1s, NAI_pol_mult_erf_ao
|
||||
integer :: i, n_pt, n_pt_out, ipoint
|
||||
double precision :: P_center(3)
|
||||
double precision :: d(0:n_pt_in), coeff, dist, const, factor
|
||||
double precision :: const_factor, dist_integral
|
||||
double precision :: accu, p_inv, p, rho, p_inv_2
|
||||
double precision :: p_new, p_new2, coef_tmp
|
||||
|
||||
ASSERT(beta .ge. 0.d0)
|
||||
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)
|
||||
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
|
||||
|
||||
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)
|
||||
! 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))
|
||||
|
||||
n_pt_in = n_pt_max_integrals
|
||||
const_factor = const_factor12 + dist * rho
|
||||
if(const_factor > 80.d0) then
|
||||
NAI_pol_mult_erf_with1s = 0.d0
|
||||
return
|
||||
endif
|
||||
|
||||
NAI_pol_mult_erf_ao_with1s = 0.d0
|
||||
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)
|
||||
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))
|
||||
|
||||
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 )
|
||||
p_new = mu_in / dsqrt(p + mu_in * mu_in)
|
||||
factor = dexp(-const_factor)
|
||||
coeff = dtwo_pi * factor * p_inv * p_new
|
||||
|
||||
NAI_pol_mult_erf_ao_with1s += integral * coef12
|
||||
enddo
|
||||
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)
|
||||
|
||||
end function NAI_pol_mult_erf_ao_with1s
|
||||
if(n_pt_out < 0) then
|
||||
NAI_pol_mult_erf_with1s = 0.d0
|
||||
return
|
||||
endif
|
||||
|
||||
! sum of integrals of type : int {t,[0,1]} exp-(rho.(P-Q)^2 * t^2) * t^i
|
||||
accu = 0.d0
|
||||
do i = 0, n_pt_out, 2
|
||||
accu += d(i) * rint(i/2, const)
|
||||
enddo
|
||||
NAI_pol_mult_erf_with1s = accu * coeff
|
||||
|
||||
end function NAI_pol_mult_erf_with1s
|
||||
|
||||
! ---
|
||||
|
||||
subroutine NAI_pol_mult_erf_with1s_v(A1_center, A2_center, power_A1, power_A2, alpha1, alpha2, beta, B_center, LD_B, C_center, LD_C, n_pt_in, mu_in, res_v, LD_resv, n_points)
|
||||
|
||||
@ -428,107 +662,6 @@ subroutine give_polynomial_mult_center_one_e_erf_opt(A_center, B_center, power_A
|
||||
end subroutine give_polynomial_mult_center_one_e_erf_opt
|
||||
|
||||
! ---
|
||||
subroutine NAI_pol_mult_erf_v(A_center, B_center, power_A, power_B, alpha, beta, C_center, LD_C, n_pt_in, mu_in, res_v, LD_resv, n_points)
|
||||
|
||||
BEGIN_DOC
|
||||
!
|
||||
! Computes the following integral :
|
||||
!
|
||||
! .. math::
|
||||
!
|
||||
! \int dr (x-A_x)^a (x-B_x)^b \exp(-\alpha (x-A_x)^2 - \beta (x-B_x)^2 )
|
||||
! \frac{\erf(\mu |r - R_C |)}{| r - R_C |}$.
|
||||
!
|
||||
END_DOC
|
||||
|
||||
include 'utils/constants.include.F'
|
||||
|
||||
implicit none
|
||||
|
||||
integer, intent(in) :: n_pt_in, n_points, LD_C, LD_resv
|
||||
integer, intent(in) :: power_A(3), power_B(3)
|
||||
double precision, intent(in) :: A_center(3), B_center(3), alpha, beta, mu_in
|
||||
double precision, intent(in) :: C_center(LD_C,3)
|
||||
double precision, intent(out) :: res_v(LD_resv)
|
||||
|
||||
integer :: i, n_pt, n_pt_out, ipoint
|
||||
double precision :: P_center(3)
|
||||
double precision :: d(0:n_pt_in), coeff, dist, const, factor
|
||||
double precision :: const_factor, dist_integral
|
||||
double precision :: accu, p_inv, p, rho, p_inv_2
|
||||
double precision :: p_new, p_new2, coef_tmp
|
||||
|
||||
double precision :: rint
|
||||
|
||||
res_V(1:LD_resv) = 0.d0
|
||||
|
||||
p = alpha + beta
|
||||
p_inv = 1.d0 / p
|
||||
p_inv_2 = 0.5d0 * p_inv
|
||||
rho = alpha * beta * p_inv
|
||||
p_new = mu_in / dsqrt(p + mu_in * mu_in)
|
||||
p_new2 = p_new * p_new
|
||||
coef_tmp = p * p_new2
|
||||
|
||||
dist = 0.d0
|
||||
do i = 1, 3
|
||||
P_center(i) = (alpha * A_center(i) + beta * B_center(i)) * p_inv
|
||||
dist += (A_center(i) - B_center(i)) * (A_center(i) - B_center(i))
|
||||
enddo
|
||||
|
||||
const_factor = dist * rho
|
||||
if(const_factor > 80.d0) then
|
||||
return
|
||||
endif
|
||||
factor = dexp(-const_factor)
|
||||
coeff = dtwo_pi * factor * p_inv * p_new
|
||||
|
||||
n_pt = 2 * ( power_A(1) + power_B(1) + power_A(2) + power_B(2) + power_A(3) + power_B(3) )
|
||||
|
||||
if(n_pt == 0) then
|
||||
|
||||
do ipoint = 1, n_points
|
||||
dist_integral = 0.d0
|
||||
do i = 1, 3
|
||||
dist_integral += (P_center(i) - C_center(ipoint,i)) * (P_center(i) - C_center(ipoint,i))
|
||||
enddo
|
||||
const = coef_tmp * dist_integral
|
||||
|
||||
res_v(ipoint) = coeff * rint(0, const)
|
||||
enddo
|
||||
|
||||
else
|
||||
|
||||
do ipoint = 1, n_points
|
||||
dist_integral = 0.d0
|
||||
do i = 1, 3
|
||||
dist_integral += (P_center(i) - C_center(ipoint,i)) * (P_center(i) - C_center(ipoint,i))
|
||||
enddo
|
||||
const = coef_tmp * dist_integral
|
||||
|
||||
do i = 0, n_pt_in
|
||||
d(i) = 0.d0
|
||||
enddo
|
||||
call give_polynomial_mult_center_one_e_erf_opt(A_center, B_center, power_A, power_B, C_center(ipoint,1:3), n_pt_in, d, n_pt_out, p_inv_2, p_new2, P_center)
|
||||
|
||||
if(n_pt_out < 0) then
|
||||
res_v(ipoint) = 0.d0
|
||||
cycle
|
||||
endif
|
||||
|
||||
! sum of integrals of type : int {t,[0,1]} exp-(rho.(P-Q)^2 * t^2) * t^i
|
||||
accu = 0.d0
|
||||
do i = 0, n_pt_out, 2
|
||||
accu += d(i) * rint(i/2, const)
|
||||
enddo
|
||||
|
||||
res_v(ipoint) = accu * coeff
|
||||
enddo
|
||||
|
||||
endif
|
||||
|
||||
end subroutine NAI_pol_mult_erf_v
|
||||
|
||||
|
||||
subroutine give_polynomial_mult_center_one_e_erf(A_center,B_center,alpha,beta,power_A,power_B,C_center,n_pt_in,d,n_pt_out,mu_in)
|
||||
|
||||
@ -659,113 +792,3 @@ subroutine give_polynomial_mult_center_one_e_erf(A_center,B_center,alpha,beta,po
|
||||
|
||||
end
|
||||
|
||||
double precision function NAI_pol_mult_erf_with1s( A1_center, A2_center, power_A1, power_A2, alpha1, alpha2 &
|
||||
, beta, B_center, C_center, n_pt_in, mu_in )
|
||||
|
||||
BEGIN_DOC
|
||||
!
|
||||
! Computes the following integral :
|
||||
!
|
||||
! .. math::
|
||||
!
|
||||
! \int dx (x - A1_x)^a_1 (x - B1_x)^a_2 \exp(-\alpha_1 (x - A1_x)^2 - \alpha_2 (x - A2_x)^2)
|
||||
! \int dy (y - A1_y)^b_1 (y - B1_y)^b_2 \exp(-\alpha_1 (y - A1_y)^2 - \alpha_2 (y - A2_y)^2)
|
||||
! \int dz (x - A1_z)^c_1 (z - B1_z)^c_2 \exp(-\alpha_1 (z - A1_z)^2 - \alpha_2 (z - A2_z)^2)
|
||||
! \exp(-\beta (r - B)^2)
|
||||
! \frac{\erf(\mu |r - R_C|)}{|r - R_C|}$.
|
||||
!
|
||||
END_DOC
|
||||
|
||||
include 'utils/constants.include.F'
|
||||
|
||||
implicit none
|
||||
integer, intent(in) :: n_pt_in
|
||||
integer, intent(in) :: power_A1(3), power_A2(3)
|
||||
double precision, intent(in) :: C_center(3), A1_center(3), A2_center(3), B_center(3)
|
||||
double precision, intent(in) :: alpha1, alpha2, beta, mu_in
|
||||
|
||||
integer :: i, n_pt, n_pt_out
|
||||
double precision :: alpha12, alpha12_inv, alpha12_inv_2, rho12, A12_center(3), dist12, const_factor12
|
||||
double precision :: p, p_inv, p_inv_2, rho, P_center(3), dist, const_factor
|
||||
double precision :: dist_integral
|
||||
double precision :: d(0:n_pt_in), coeff, const, factor
|
||||
double precision :: accu
|
||||
double precision :: p_new
|
||||
|
||||
double precision :: rint
|
||||
|
||||
|
||||
! e^{-alpha1 (r - A1)^2} e^{-alpha2 (r - A2)^2} = e^{-K12} e^{-alpha12 (r - A12)^2}
|
||||
alpha12 = alpha1 + alpha2
|
||||
alpha12_inv = 1.d0 / alpha12
|
||||
alpha12_inv_2 = 0.5d0 * alpha12_inv
|
||||
rho12 = alpha1 * alpha2 * alpha12_inv
|
||||
A12_center(1) = (alpha1 * A1_center(1) + alpha2 * A2_center(1)) * alpha12_inv
|
||||
A12_center(2) = (alpha1 * A1_center(2) + alpha2 * A2_center(2)) * alpha12_inv
|
||||
A12_center(3) = (alpha1 * A1_center(3) + alpha2 * A2_center(3)) * alpha12_inv
|
||||
dist12 = (A1_center(1) - A2_center(1)) * (A1_center(1) - A2_center(1)) &
|
||||
+ (A1_center(2) - A2_center(2)) * (A1_center(2) - A2_center(2)) &
|
||||
+ (A1_center(3) - A2_center(3)) * (A1_center(3) - A2_center(3))
|
||||
|
||||
const_factor12 = dist12 * rho12
|
||||
if(const_factor12 > 80.d0) then
|
||||
NAI_pol_mult_erf_with1s = 0.d0
|
||||
return
|
||||
endif
|
||||
|
||||
! ---
|
||||
|
||||
! e^{-K12} e^{-alpha12 (r - A12)^2} e^{-beta (r - B)^2} = e^{-K} e^{-p (r - P)^2}
|
||||
p = alpha12 + beta
|
||||
p_inv = 1.d0 / p
|
||||
p_inv_2 = 0.5d0 * p_inv
|
||||
rho = alpha12 * beta * p_inv
|
||||
P_center(1) = (alpha12 * A12_center(1) + beta * B_center(1)) * p_inv
|
||||
P_center(2) = (alpha12 * A12_center(2) + beta * B_center(2)) * p_inv
|
||||
P_center(3) = (alpha12 * A12_center(3) + beta * B_center(3)) * p_inv
|
||||
dist = (A12_center(1) - B_center(1)) * (A12_center(1) - B_center(1)) &
|
||||
+ (A12_center(2) - B_center(2)) * (A12_center(2) - B_center(2)) &
|
||||
+ (A12_center(3) - B_center(3)) * (A12_center(3) - B_center(3))
|
||||
|
||||
const_factor = const_factor12 + dist * rho
|
||||
if(const_factor > 80.d0) then
|
||||
NAI_pol_mult_erf_with1s = 0.d0
|
||||
return
|
||||
endif
|
||||
|
||||
dist_integral = (P_center(1) - C_center(1)) * (P_center(1) - C_center(1)) &
|
||||
+ (P_center(2) - C_center(2)) * (P_center(2) - C_center(2)) &
|
||||
+ (P_center(3) - C_center(3)) * (P_center(3) - C_center(3))
|
||||
|
||||
! ---
|
||||
|
||||
p_new = mu_in / dsqrt(p + mu_in * mu_in)
|
||||
factor = dexp(-const_factor)
|
||||
coeff = dtwo_pi * factor * p_inv * p_new
|
||||
|
||||
n_pt = 2 * ( (power_A1(1) + power_A2(1)) + (power_A1(2) + power_A2(2)) + (power_A1(3) + power_A2(3)) )
|
||||
const = p * dist_integral * p_new * p_new
|
||||
if(n_pt == 0) then
|
||||
NAI_pol_mult_erf_with1s = coeff * rint(0, const)
|
||||
return
|
||||
endif
|
||||
|
||||
do i = 0, n_pt_in
|
||||
d(i) = 0.d0
|
||||
enddo
|
||||
p_new = p_new * p_new
|
||||
call give_polynomial_mult_center_one_e_erf_opt( A1_center, A2_center, power_A1, power_A2, C_center, n_pt_in, d, n_pt_out, p_inv_2, p_new, P_center)
|
||||
|
||||
if(n_pt_out < 0) then
|
||||
NAI_pol_mult_erf_with1s = 0.d0
|
||||
return
|
||||
endif
|
||||
|
||||
! sum of integrals of type : int {t,[0,1]} exp-(rho.(P-Q)^2 * t^2) * t^i
|
||||
accu = 0.d0
|
||||
do i = 0, n_pt_out, 2
|
||||
accu += d(i) * rint(i/2, const)
|
||||
enddo
|
||||
NAI_pol_mult_erf_with1s = accu * coeff
|
||||
|
||||
end function NAI_pol_mult_erf_with1s
|
||||
|
@ -1,4 +1,8 @@
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [ double precision, ao_integrals_n_e, (ao_num,ao_num)]
|
||||
|
||||
BEGIN_DOC
|
||||
! Nucleus-electron interaction, in the |AO| basis set.
|
||||
!
|
||||
@ -6,84 +10,100 @@ BEGIN_PROVIDER [ double precision, ao_integrals_n_e, (ao_num,ao_num)]
|
||||
!
|
||||
! These integrals also contain the pseudopotential integrals.
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
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,n_pt_in,m
|
||||
double precision :: overlap_x,overlap_y,overlap_z,overlap,dx,NAI_pol_mult
|
||||
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_n_e = 0.d0
|
||||
|
||||
if (read_ao_integrals_n_e) then
|
||||
|
||||
call ezfio_get_ao_one_e_ints_ao_integrals_n_e(ao_integrals_n_e)
|
||||
print *, 'AO N-e integrals read from disk'
|
||||
|
||||
else
|
||||
|
||||
ao_integrals_n_e = 0.d0
|
||||
if(use_cosgtos) then
|
||||
!print *, " use_cosgtos for ao_integrals_n_e ?", use_cosgtos
|
||||
|
||||
! _
|
||||
! /| / |_)
|
||||
! | / | \
|
||||
!
|
||||
do j = 1, ao_num
|
||||
do i = 1, ao_num
|
||||
ao_integrals_n_e(i,j) = ao_integrals_n_e_cosgtos(i,j)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
!$OMP PARALLEL &
|
||||
!$OMP DEFAULT (NONE) &
|
||||
!$OMP PRIVATE (i,j,k,l,m,alpha,beta,A_center,B_center,C_center,power_A,power_B,&
|
||||
!$OMP num_A,num_B,Z,c,n_pt_in) &
|
||||
!$OMP SHARED (ao_num,ao_prim_num,ao_expo_ordered_transp,ao_power,ao_nucl,nucl_coord,ao_coef_normalized_ordered_transp,&
|
||||
!$OMP n_pt_max_integrals,ao_integrals_n_e,nucl_num,nucl_charge)
|
||||
else
|
||||
|
||||
n_pt_in = n_pt_max_integrals
|
||||
!$OMP PARALLEL &
|
||||
!$OMP DEFAULT (NONE) &
|
||||
!$OMP PRIVATE (i,j,k,l,m,alpha,beta,A_center,B_center,C_center,power_A,power_B,&
|
||||
!$OMP num_A,num_B,Z,c,c1,n_pt_in) &
|
||||
!$OMP SHARED (ao_num,ao_prim_num,ao_expo_ordered_transp,ao_power,ao_nucl,nucl_coord,ao_coef_normalized_ordered_transp,&
|
||||
!$OMP n_pt_max_integrals,ao_integrals_n_e,nucl_num,nucl_charge)
|
||||
|
||||
!$OMP DO SCHEDULE (dynamic)
|
||||
n_pt_in = n_pt_max_integrals
|
||||
|
||||
do j = 1, ao_num
|
||||
num_A = ao_nucl(j)
|
||||
power_A(1:3)= ao_power(j,1:3)
|
||||
A_center(1:3) = nucl_coord(num_A,1:3)
|
||||
!$OMP DO SCHEDULE (dynamic)
|
||||
|
||||
do i = 1, ao_num
|
||||
do j = 1, ao_num
|
||||
num_A = ao_nucl(j)
|
||||
power_A(1:3)= ao_power(j,1:3)
|
||||
A_center(1:3) = nucl_coord(num_A,1:3)
|
||||
|
||||
num_B = ao_nucl(i)
|
||||
power_B(1:3)= ao_power(i,1:3)
|
||||
B_center(1:3) = nucl_coord(num_B,1:3)
|
||||
do i = 1, ao_num
|
||||
|
||||
do l=1,ao_prim_num(j)
|
||||
alpha = ao_expo_ordered_transp(l,j)
|
||||
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 m=1,ao_prim_num(i)
|
||||
beta = ao_expo_ordered_transp(m,i)
|
||||
do l=1,ao_prim_num(j)
|
||||
alpha = ao_expo_ordered_transp(l,j)
|
||||
|
||||
double precision :: c
|
||||
c = 0.d0
|
||||
do m=1,ao_prim_num(i)
|
||||
beta = ao_expo_ordered_transp(m,i)
|
||||
|
||||
do k = 1, nucl_num
|
||||
double precision :: Z
|
||||
Z = nucl_charge(k)
|
||||
double precision :: c, c1
|
||||
c = 0.d0
|
||||
|
||||
C_center(1:3) = nucl_coord(k,1:3)
|
||||
do k = 1, nucl_num
|
||||
double precision :: Z
|
||||
Z = nucl_charge(k)
|
||||
|
||||
c = c - Z * NAI_pol_mult(A_center,B_center, &
|
||||
power_A,power_B,alpha,beta,C_center,n_pt_in)
|
||||
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
|
||||
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
|
||||
IF (DO_PSEUDO) THEN
|
||||
|
||||
endif
|
||||
|
||||
|
||||
IF(do_pseudo) THEN
|
||||
ao_integrals_n_e += ao_pseudo_integrals
|
||||
ENDIF
|
||||
IF(point_charges) THEN
|
||||
ao_integrals_n_e += ao_integrals_pt_chrg
|
||||
ENDIF
|
||||
|
||||
|
||||
endif
|
||||
|
||||
@ -102,7 +122,7 @@ BEGIN_PROVIDER [ double precision, ao_integrals_n_e_imag, (ao_num,ao_num)]
|
||||
! :math:`\langle \chi_i | -\sum_A \frac{1}{|r-R_A|} | \chi_j \rangle`
|
||||
END_DOC
|
||||
implicit none
|
||||
double precision :: alpha, beta, gama, delta
|
||||
double precision :: alpha, beta
|
||||
integer :: num_A,num_B
|
||||
double precision :: A_center(3),B_center(3),C_center(3)
|
||||
integer :: power_A(3),power_B(3)
|
||||
@ -125,7 +145,7 @@ BEGIN_PROVIDER [ double precision, ao_integrals_n_e_per_atom, (ao_num,ao_num,nuc
|
||||
! :math:`\langle \chi_i | -\frac{1}{|r-R_A|} | \chi_j \rangle`
|
||||
END_DOC
|
||||
implicit none
|
||||
double precision :: alpha, beta, gama, delta
|
||||
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)
|
||||
@ -268,6 +288,7 @@ double precision function NAI_pol_mult(A_center,B_center,power_A,power_B,alpha,b
|
||||
|
||||
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
|
||||
@ -579,61 +600,3 @@ double precision function V_r(n,alpha)
|
||||
end
|
||||
|
||||
|
||||
double precision function V_phi(n,m)
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Computes the angular $\phi$ part of the nuclear attraction integral:
|
||||
!
|
||||
! $\int_{0}^{2 \pi} \cos(\phi)^n \sin(\phi)^m d\phi$.
|
||||
END_DOC
|
||||
integer :: n,m, i
|
||||
double precision :: prod, Wallis
|
||||
prod = 1.d0
|
||||
do i = 0,shiftr(n,1)-1
|
||||
prod = prod/ (1.d0 + dfloat(m+1)/dfloat(n-i-i-1))
|
||||
enddo
|
||||
V_phi = 4.d0 * prod * Wallis(m)
|
||||
end
|
||||
|
||||
|
||||
double precision function V_theta(n,m)
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Computes the angular $\theta$ part of the nuclear attraction integral:
|
||||
!
|
||||
! $\int_{0}^{\pi} \cos(\theta)^n \sin(\theta)^m d\theta$
|
||||
END_DOC
|
||||
integer :: n,m,i
|
||||
double precision :: Wallis, prod
|
||||
include 'utils/constants.include.F'
|
||||
V_theta = 0.d0
|
||||
prod = 1.d0
|
||||
do i = 0,shiftr(n,1)-1
|
||||
prod = prod / (1.d0 + dfloat(m+1)/dfloat(n-i-i-1))
|
||||
enddo
|
||||
V_theta = (prod+prod) * Wallis(m)
|
||||
end
|
||||
|
||||
|
||||
double precision function Wallis(n)
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Wallis integral:
|
||||
!
|
||||
! $\int_{0}^{\pi} \cos(\theta)^n d\theta$.
|
||||
END_DOC
|
||||
double precision :: fact
|
||||
integer :: n,p
|
||||
include 'utils/constants.include.F'
|
||||
if(iand(n,1).eq.0)then
|
||||
Wallis = fact(shiftr(n,1))
|
||||
Wallis = pi * fact(n) / (dble(ibset(0_8,n)) * (Wallis+Wallis)*Wallis)
|
||||
else
|
||||
p = shiftr(n,1)
|
||||
Wallis = fact(p)
|
||||
Wallis = dble(ibset(0_8,p+p)) * Wallis*Wallis / fact(p+p+1)
|
||||
endif
|
||||
|
||||
end
|
||||
|
||||
|
||||
|
@ -4,13 +4,6 @@ doc: Read/Write |AO| integrals from/to disk [ Write | Read | None ]
|
||||
interface: ezfio,provider,ocaml
|
||||
default: None
|
||||
|
||||
[ao_integrals_threshold]
|
||||
type: Threshold
|
||||
doc: If | (pq|rs) | < `ao_integrals_threshold` then (pq|rs) is zero
|
||||
interface: ezfio,provider,ocaml
|
||||
default: 1.e-15
|
||||
ezfio_name: threshold_ao
|
||||
|
||||
[do_direct_integrals]
|
||||
type: logical
|
||||
doc: Compute integrals on the fly (very slow, only for debugging)
|
||||
|
@ -1,102 +1,123 @@
|
||||
double precision function ao_two_e_integral(i,j,k,l)
|
||||
implicit none
|
||||
|
||||
! ---
|
||||
|
||||
double precision function ao_two_e_integral(i, j, k, l)
|
||||
|
||||
BEGIN_DOC
|
||||
! integral of the AO basis <ik|jl> or (ij|kl)
|
||||
! i(r1) j(r1) 1/r12 k(r2) l(r2)
|
||||
END_DOC
|
||||
|
||||
integer,intent(in) :: i,j,k,l
|
||||
integer :: p,q,r,s
|
||||
double precision :: I_center(3),J_center(3),K_center(3),L_center(3)
|
||||
integer :: num_i,num_j,num_k,num_l,dim1,I_power(3),J_power(3),K_power(3),L_power(3)
|
||||
double precision :: integral
|
||||
implicit none
|
||||
include 'utils/constants.include.F'
|
||||
|
||||
integer, intent(in) :: i, j, k, l
|
||||
|
||||
integer :: p, q, r, s
|
||||
integer :: num_i,num_j,num_k,num_l,dim1,I_power(3),J_power(3),K_power(3),L_power(3)
|
||||
integer :: iorder_p(3), iorder_q(3)
|
||||
double precision :: I_center(3), J_center(3), K_center(3), L_center(3)
|
||||
double precision :: integral
|
||||
double precision :: P_new(0:max_dim,3),P_center(3),fact_p,pp
|
||||
double precision :: Q_new(0:max_dim,3),Q_center(3),fact_q,qq
|
||||
integer :: iorder_p(3), iorder_q(3)
|
||||
|
||||
double precision :: ao_two_e_integral_schwartz_accel
|
||||
|
||||
if (ao_prim_num(i) * ao_prim_num(j) * ao_prim_num(k) * ao_prim_num(l) > 1024 ) then
|
||||
ao_two_e_integral = ao_two_e_integral_schwartz_accel(i,j,k,l)
|
||||
else
|
||||
double precision :: ao_two_e_integral_cosgtos
|
||||
|
||||
dim1 = n_pt_max_integrals
|
||||
|
||||
num_i = ao_nucl(i)
|
||||
num_j = ao_nucl(j)
|
||||
num_k = ao_nucl(k)
|
||||
num_l = ao_nucl(l)
|
||||
ao_two_e_integral = 0.d0
|
||||
if(use_cosgtos) then
|
||||
!print *, ' use_cosgtos for ao_two_e_integral ?', use_cosgtos
|
||||
|
||||
if (num_i /= num_j .or. num_k /= num_l .or. num_j /= num_k)then
|
||||
do p = 1, 3
|
||||
I_power(p) = ao_power(i,p)
|
||||
J_power(p) = ao_power(j,p)
|
||||
K_power(p) = ao_power(k,p)
|
||||
L_power(p) = ao_power(l,p)
|
||||
I_center(p) = nucl_coord(num_i,p)
|
||||
J_center(p) = nucl_coord(num_j,p)
|
||||
K_center(p) = nucl_coord(num_k,p)
|
||||
L_center(p) = nucl_coord(num_l,p)
|
||||
enddo
|
||||
ao_two_e_integral = ao_two_e_integral_cosgtos(i, j, k, l)
|
||||
|
||||
double precision :: coef1, coef2, coef3, coef4
|
||||
double precision :: p_inv,q_inv
|
||||
double precision :: general_primitive_integral
|
||||
else
|
||||
|
||||
do p = 1, ao_prim_num(i)
|
||||
coef1 = ao_coef_normalized_ordered_transp(p,i)
|
||||
do q = 1, ao_prim_num(j)
|
||||
coef2 = coef1*ao_coef_normalized_ordered_transp(q,j)
|
||||
call give_explicit_poly_and_gaussian(P_new,P_center,pp,fact_p,iorder_p,&
|
||||
ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j), &
|
||||
I_power,J_power,I_center,J_center,dim1)
|
||||
p_inv = 1.d0/pp
|
||||
do r = 1, ao_prim_num(k)
|
||||
coef3 = coef2*ao_coef_normalized_ordered_transp(r,k)
|
||||
do s = 1, ao_prim_num(l)
|
||||
coef4 = coef3*ao_coef_normalized_ordered_transp(s,l)
|
||||
call give_explicit_poly_and_gaussian(Q_new,Q_center,qq,fact_q,iorder_q,&
|
||||
ao_expo_ordered_transp(r,k),ao_expo_ordered_transp(s,l), &
|
||||
K_power,L_power,K_center,L_center,dim1)
|
||||
q_inv = 1.d0/qq
|
||||
integral = general_primitive_integral(dim1, &
|
||||
P_new,P_center,fact_p,pp,p_inv,iorder_p, &
|
||||
Q_new,Q_center,fact_q,qq,q_inv,iorder_q)
|
||||
ao_two_e_integral = ao_two_e_integral + coef4 * integral
|
||||
enddo ! s
|
||||
enddo ! r
|
||||
enddo ! q
|
||||
enddo ! p
|
||||
if (ao_prim_num(i) * ao_prim_num(j) * ao_prim_num(k) * ao_prim_num(l) > 1024 ) then
|
||||
|
||||
ao_two_e_integral = ao_two_e_integral_schwartz_accel(i,j,k,l)
|
||||
|
||||
else
|
||||
|
||||
do p = 1, 3
|
||||
I_power(p) = ao_power(i,p)
|
||||
J_power(p) = ao_power(j,p)
|
||||
K_power(p) = ao_power(k,p)
|
||||
L_power(p) = ao_power(l,p)
|
||||
enddo
|
||||
double precision :: ERI
|
||||
dim1 = n_pt_max_integrals
|
||||
|
||||
do p = 1, ao_prim_num(i)
|
||||
coef1 = ao_coef_normalized_ordered_transp(p,i)
|
||||
do q = 1, ao_prim_num(j)
|
||||
coef2 = coef1*ao_coef_normalized_ordered_transp(q,j)
|
||||
do r = 1, ao_prim_num(k)
|
||||
coef3 = coef2*ao_coef_normalized_ordered_transp(r,k)
|
||||
do s = 1, ao_prim_num(l)
|
||||
coef4 = coef3*ao_coef_normalized_ordered_transp(s,l)
|
||||
integral = ERI( &
|
||||
ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j),ao_expo_ordered_transp(r,k),ao_expo_ordered_transp(s,l),&
|
||||
I_power(1),J_power(1),K_power(1),L_power(1), &
|
||||
I_power(2),J_power(2),K_power(2),L_power(2), &
|
||||
I_power(3),J_power(3),K_power(3),L_power(3))
|
||||
ao_two_e_integral = ao_two_e_integral + coef4 * integral
|
||||
enddo ! s
|
||||
enddo ! r
|
||||
enddo ! q
|
||||
enddo ! p
|
||||
num_i = ao_nucl(i)
|
||||
num_j = ao_nucl(j)
|
||||
num_k = ao_nucl(k)
|
||||
num_l = ao_nucl(l)
|
||||
ao_two_e_integral = 0.d0
|
||||
|
||||
if (num_i /= num_j .or. num_k /= num_l .or. num_j /= num_k)then
|
||||
do p = 1, 3
|
||||
I_power(p) = ao_power(i,p)
|
||||
J_power(p) = ao_power(j,p)
|
||||
K_power(p) = ao_power(k,p)
|
||||
L_power(p) = ao_power(l,p)
|
||||
I_center(p) = nucl_coord(num_i,p)
|
||||
J_center(p) = nucl_coord(num_j,p)
|
||||
K_center(p) = nucl_coord(num_k,p)
|
||||
L_center(p) = nucl_coord(num_l,p)
|
||||
enddo
|
||||
|
||||
double precision :: coef1, coef2, coef3, coef4
|
||||
double precision :: p_inv,q_inv
|
||||
double precision :: general_primitive_integral
|
||||
|
||||
do p = 1, ao_prim_num(i)
|
||||
coef1 = ao_coef_normalized_ordered_transp(p,i)
|
||||
do q = 1, ao_prim_num(j)
|
||||
coef2 = coef1*ao_coef_normalized_ordered_transp(q,j)
|
||||
call give_explicit_poly_and_gaussian(P_new,P_center,pp,fact_p,iorder_p,&
|
||||
ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j), &
|
||||
I_power,J_power,I_center,J_center,dim1)
|
||||
p_inv = 1.d0/pp
|
||||
do r = 1, ao_prim_num(k)
|
||||
coef3 = coef2*ao_coef_normalized_ordered_transp(r,k)
|
||||
do s = 1, ao_prim_num(l)
|
||||
coef4 = coef3*ao_coef_normalized_ordered_transp(s,l)
|
||||
call give_explicit_poly_and_gaussian(Q_new,Q_center,qq,fact_q,iorder_q,&
|
||||
ao_expo_ordered_transp(r,k),ao_expo_ordered_transp(s,l), &
|
||||
K_power,L_power,K_center,L_center,dim1)
|
||||
q_inv = 1.d0/qq
|
||||
integral = general_primitive_integral(dim1, &
|
||||
P_new,P_center,fact_p,pp,p_inv,iorder_p, &
|
||||
Q_new,Q_center,fact_q,qq,q_inv,iorder_q)
|
||||
ao_two_e_integral = ao_two_e_integral + coef4 * integral
|
||||
enddo ! s
|
||||
enddo ! r
|
||||
enddo ! q
|
||||
enddo ! p
|
||||
|
||||
else
|
||||
|
||||
do p = 1, 3
|
||||
I_power(p) = ao_power(i,p)
|
||||
J_power(p) = ao_power(j,p)
|
||||
K_power(p) = ao_power(k,p)
|
||||
L_power(p) = ao_power(l,p)
|
||||
enddo
|
||||
double precision :: ERI
|
||||
|
||||
do p = 1, ao_prim_num(i)
|
||||
coef1 = ao_coef_normalized_ordered_transp(p,i)
|
||||
do q = 1, ao_prim_num(j)
|
||||
coef2 = coef1*ao_coef_normalized_ordered_transp(q,j)
|
||||
do r = 1, ao_prim_num(k)
|
||||
coef3 = coef2*ao_coef_normalized_ordered_transp(r,k)
|
||||
do s = 1, ao_prim_num(l)
|
||||
coef4 = coef3*ao_coef_normalized_ordered_transp(s,l)
|
||||
integral = ERI( &
|
||||
ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j),ao_expo_ordered_transp(r,k),ao_expo_ordered_transp(s,l),&
|
||||
I_power(1),J_power(1),K_power(1),L_power(1), &
|
||||
I_power(2),J_power(2),K_power(2),L_power(2), &
|
||||
I_power(3),J_power(3),K_power(3),L_power(3))
|
||||
ao_two_e_integral = ao_two_e_integral + coef4 * integral
|
||||
enddo ! s
|
||||
enddo ! r
|
||||
enddo ! q
|
||||
enddo ! p
|
||||
|
||||
endif
|
||||
|
||||
endif
|
||||
|
||||
@ -104,6 +125,8 @@ double precision function ao_two_e_integral(i,j,k,l)
|
||||
|
||||
end
|
||||
|
||||
! ---
|
||||
|
||||
double precision function ao_two_e_integral_schwartz_accel(i,j,k,l)
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
@ -421,14 +444,17 @@ BEGIN_PROVIDER [ logical, ao_two_e_integrals_in_map ]
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ double precision, ao_two_e_integral_schwartz,(ao_num,ao_num) ]
|
||||
implicit none
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [ double precision, ao_two_e_integral_schwartz, (ao_num, ao_num) ]
|
||||
|
||||
BEGIN_DOC
|
||||
! Needed to compute Schwartz inequalities
|
||||
END_DOC
|
||||
|
||||
integer :: i,k
|
||||
double precision :: ao_two_e_integral,cpu_1,cpu_2, wall_1, wall_2
|
||||
implicit none
|
||||
integer :: i, k
|
||||
double precision :: ao_two_e_integral,cpu_1,cpu_2, wall_1, wall_2
|
||||
|
||||
ao_two_e_integral_schwartz(1,1) = ao_two_e_integral(1,1,1,1)
|
||||
!$OMP PARALLEL DO PRIVATE(i,k) &
|
||||
@ -445,6 +471,7 @@ BEGIN_PROVIDER [ double precision, ao_two_e_integral_schwartz,(ao_num,ao_num) ]
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
double precision function general_primitive_integral(dim, &
|
||||
P_new,P_center,fact_p,p,p_inv,iorder_p, &
|
||||
|
19
src/cosgtos_ao_int/EZFIO.cfg
Normal file
19
src/cosgtos_ao_int/EZFIO.cfg
Normal file
@ -0,0 +1,19 @@
|
||||
[ao_expoim_cosgtos]
|
||||
type: double precision
|
||||
doc: imag part for Exponents for each primitive of each cosGTOs |AO|
|
||||
size: (ao_basis.ao_num,ao_basis.ao_prim_num_max)
|
||||
interface: ezfio, provider
|
||||
|
||||
[use_cosgtos]
|
||||
type: logical
|
||||
doc: If true, use cosgtos for AO integrals
|
||||
interface: ezfio,provider,ocaml
|
||||
default: False
|
||||
|
||||
[ao_integrals_threshold]
|
||||
type: Threshold
|
||||
doc: If | (pq|rs) | < `ao_integrals_threshold` then (pq|rs) is zero
|
||||
interface: ezfio,provider,ocaml
|
||||
default: 1.e-15
|
||||
ezfio_name: threshold_ao
|
||||
|
4
src/cosgtos_ao_int/README.rst
Normal file
4
src/cosgtos_ao_int/README.rst
Normal file
@ -0,0 +1,4 @@
|
||||
==============
|
||||
cosgtos_ao_int
|
||||
==============
|
||||
|
210
src/cosgtos_ao_int/aos_cosgtos.irp.f
Normal file
210
src/cosgtos_ao_int/aos_cosgtos.irp.f
Normal file
@ -0,0 +1,210 @@
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [ double precision, ao_coef_norm_ord_transp_cosgtos, (ao_prim_num_max, ao_num) ]
|
||||
|
||||
implicit none
|
||||
integer :: i, j
|
||||
|
||||
do j = 1, ao_num
|
||||
do i = 1, ao_prim_num_max
|
||||
ao_coef_norm_ord_transp_cosgtos(i,j) = ao_coef_norm_ord_cosgtos(j,i)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [ complex*16, ao_expo_ord_transp_cosgtos, (ao_prim_num_max, ao_num) ]
|
||||
|
||||
implicit none
|
||||
integer :: i, j
|
||||
|
||||
do j = 1, ao_num
|
||||
do i = 1, ao_prim_num_max
|
||||
ao_expo_ord_transp_cosgtos(i,j) = ao_expo_ord_cosgtos(j,i)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [ double precision, ao_coef_norm_cosgtos, (ao_num, ao_prim_num_max) ]
|
||||
|
||||
implicit none
|
||||
|
||||
integer :: i, j, powA(3), nz
|
||||
double precision :: norm
|
||||
complex*16 :: overlap_x, overlap_y, overlap_z, C_A(3)
|
||||
complex*16 :: integ1, integ2, expo
|
||||
|
||||
nz = 100
|
||||
|
||||
C_A(1) = (0.d0, 0.d0)
|
||||
C_A(2) = (0.d0, 0.d0)
|
||||
C_A(3) = (0.d0, 0.d0)
|
||||
|
||||
ao_coef_norm_cosgtos = 0.d0
|
||||
|
||||
do i = 1, ao_num
|
||||
|
||||
powA(1) = ao_power(i,1)
|
||||
powA(2) = ao_power(i,2)
|
||||
powA(3) = ao_power(i,3)
|
||||
|
||||
! Normalization of the primitives
|
||||
if(primitives_normalized) then
|
||||
|
||||
do j = 1, ao_prim_num(i)
|
||||
|
||||
expo = ao_expo(i,j) + (0.d0, 1.d0) * ao_expoim_cosgtos(i,j)
|
||||
|
||||
call overlap_cgaussian_xyz(C_A, C_A, expo, expo, powA, powA, overlap_x, overlap_y, overlap_z, integ1, nz)
|
||||
call overlap_cgaussian_xyz(C_A, C_A, conjg(expo), expo, powA, powA, overlap_x, overlap_y, overlap_z, integ2, nz)
|
||||
|
||||
norm = 2.d0 * real( integ1 + integ2 )
|
||||
|
||||
ao_coef_norm_cosgtos(i,j) = ao_coef(i,j) / dsqrt(norm)
|
||||
enddo
|
||||
|
||||
else
|
||||
|
||||
do j = 1, ao_prim_num(i)
|
||||
ao_coef_norm_cosgtos(i,j) = ao_coef(i,j)
|
||||
enddo
|
||||
|
||||
endif
|
||||
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [ double precision, ao_coef_norm_ord_cosgtos, (ao_num, ao_prim_num_max) ]
|
||||
&BEGIN_PROVIDER [ complex*16 , ao_expo_ord_cosgtos, (ao_num, ao_prim_num_max) ]
|
||||
|
||||
implicit none
|
||||
integer :: i, j
|
||||
integer :: iorder(ao_prim_num_max)
|
||||
double precision :: d(ao_prim_num_max,3)
|
||||
|
||||
d = 0.d0
|
||||
|
||||
do i = 1, ao_num
|
||||
|
||||
do j = 1, ao_prim_num(i)
|
||||
iorder(j) = j
|
||||
d(j,1) = ao_expo(i,j)
|
||||
d(j,2) = ao_coef_norm_cosgtos(i,j)
|
||||
d(j,3) = ao_expoim_cosgtos(i,j)
|
||||
enddo
|
||||
|
||||
call dsort (d(1,1), iorder, ao_prim_num(i))
|
||||
call dset_order(d(1,2), iorder, ao_prim_num(i))
|
||||
call dset_order(d(1,3), iorder, ao_prim_num(i))
|
||||
|
||||
do j = 1, ao_prim_num(i)
|
||||
ao_expo_ord_cosgtos (i,j) = d(j,1) + (0.d0, 1.d0) * d(j,3)
|
||||
ao_coef_norm_ord_cosgtos(i,j) = d(j,2)
|
||||
enddo
|
||||
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [ double precision, ao_overlap_cosgtos, (ao_num, ao_num) ]
|
||||
&BEGIN_PROVIDER [ double precision, ao_overlap_cosgtos_x, (ao_num, ao_num) ]
|
||||
&BEGIN_PROVIDER [ double precision, ao_overlap_cosgtos_y, (ao_num, ao_num) ]
|
||||
&BEGIN_PROVIDER [ double precision, ao_overlap_cosgtos_z, (ao_num, ao_num) ]
|
||||
|
||||
implicit none
|
||||
integer :: i, j, n, l, dim1, power_A(3), power_B(3)
|
||||
double precision :: c, overlap, overlap_x, overlap_y, overlap_z
|
||||
complex*16 :: alpha, beta, A_center(3), B_center(3)
|
||||
complex*16 :: overlap1, overlap_x1, overlap_y1, overlap_z1
|
||||
complex*16 :: overlap2, overlap_x2, overlap_y2, overlap_z2
|
||||
|
||||
ao_overlap_cosgtos = 0.d0
|
||||
ao_overlap_cosgtos_x = 0.d0
|
||||
ao_overlap_cosgtos_y = 0.d0
|
||||
ao_overlap_cosgtos_z = 0.d0
|
||||
|
||||
dim1 = 100
|
||||
|
||||
!$OMP PARALLEL DO SCHEDULE(GUIDED) &
|
||||
!$OMP DEFAULT(NONE) &
|
||||
!$OMP PRIVATE( A_center, B_center, power_A, power_B, alpha, beta, i, j, n, l, c &
|
||||
!$OMP , overlap_x , overlap_y , overlap_z , overlap &
|
||||
!$OMP , overlap_x1, overlap_y1, overlap_z1, overlap1 &
|
||||
!$OMP , overlap_x2, overlap_y2, overlap_z2, overlap2 ) &
|
||||
!$OMP SHARED( nucl_coord, ao_power, ao_prim_num, ao_num, ao_nucl, dim1 &
|
||||
!$OMP , ao_overlap_cosgtos_x, ao_overlap_cosgtos_y, ao_overlap_cosgtos_z, ao_overlap_cosgtos &
|
||||
!$OMP , ao_coef_norm_ord_transp_cosgtos, ao_expo_ord_transp_cosgtos )
|
||||
|
||||
do j = 1, ao_num
|
||||
|
||||
A_center(1) = nucl_coord(ao_nucl(j),1) * (1.d0, 0.d0)
|
||||
A_center(2) = nucl_coord(ao_nucl(j),2) * (1.d0, 0.d0)
|
||||
A_center(3) = nucl_coord(ao_nucl(j),3) * (1.d0, 0.d0)
|
||||
power_A(1) = ao_power(j,1)
|
||||
power_A(2) = ao_power(j,2)
|
||||
power_A(3) = ao_power(j,3)
|
||||
|
||||
do i = 1, ao_num
|
||||
|
||||
B_center(1) = nucl_coord(ao_nucl(i),1) * (1.d0, 0.d0)
|
||||
B_center(2) = nucl_coord(ao_nucl(i),2) * (1.d0, 0.d0)
|
||||
B_center(3) = nucl_coord(ao_nucl(i),3) * (1.d0, 0.d0)
|
||||
power_B(1) = ao_power(i,1)
|
||||
power_B(2) = ao_power(i,2)
|
||||
power_B(3) = ao_power(i,3)
|
||||
|
||||
do n = 1, ao_prim_num(j)
|
||||
alpha = ao_expo_ord_transp_cosgtos(n,j)
|
||||
|
||||
do l = 1, ao_prim_num(i)
|
||||
c = ao_coef_norm_ord_transp_cosgtos(n,j) * ao_coef_norm_ord_transp_cosgtos(l,i)
|
||||
beta = ao_expo_ord_transp_cosgtos(l,i)
|
||||
|
||||
call overlap_cgaussian_xyz( A_center, B_center, alpha, beta, power_A, power_B &
|
||||
, overlap_x1, overlap_y1, overlap_z1, overlap1, dim1 )
|
||||
|
||||
call overlap_cgaussian_xyz( A_center, B_center, conjg(alpha), beta, power_A, power_B &
|
||||
, overlap_x2, overlap_y2, overlap_z2, overlap2, dim1 )
|
||||
|
||||
overlap_x = 2.d0 * real( overlap_x1 + overlap_x2 )
|
||||
overlap_y = 2.d0 * real( overlap_y1 + overlap_y2 )
|
||||
overlap_z = 2.d0 * real( overlap_z1 + overlap_z2 )
|
||||
overlap = 2.d0 * real( overlap1 + overlap2 )
|
||||
|
||||
ao_overlap_cosgtos(i,j) = ao_overlap_cosgtos(i,j) + c * overlap
|
||||
|
||||
if( isnan(ao_overlap_cosgtos(i,j)) ) then
|
||||
print*,'i, j', i, j
|
||||
print*,'l, n', l, n
|
||||
print*,'c, overlap', c, overlap
|
||||
print*, overlap_x, overlap_y, overlap_z
|
||||
stop
|
||||
endif
|
||||
|
||||
ao_overlap_cosgtos_x(i,j) = ao_overlap_cosgtos_x(i,j) + c * overlap_x
|
||||
ao_overlap_cosgtos_y(i,j) = ao_overlap_cosgtos_y(i,j) + c * overlap_y
|
||||
ao_overlap_cosgtos_z(i,j) = ao_overlap_cosgtos_z(i,j) + c * overlap_z
|
||||
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END PARALLEL DO
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
|
||||
|
7
src/cosgtos_ao_int/cosgtos_ao_int.irp.f
Normal file
7
src/cosgtos_ao_int/cosgtos_ao_int.irp.f
Normal file
@ -0,0 +1,7 @@
|
||||
program cosgtos_ao_int
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! TODO : Put the documentation of the program here
|
||||
END_DOC
|
||||
print *, 'Hello world'
|
||||
end
|
535
src/cosgtos_ao_int/one_e_Coul_integrals.irp.f
Normal file
535
src/cosgtos_ao_int/one_e_Coul_integrals.irp.f
Normal file
@ -0,0 +1,535 @@
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [ double precision, ao_integrals_n_e_cosgtos, (ao_num, ao_num)]
|
||||
|
||||
BEGIN_DOC
|
||||
!
|
||||
! Nucleus-electron interaction, in the cosgtos |AO| basis set.
|
||||
!
|
||||
! :math:`\langle \chi_i | -\sum_A \frac{1}{|r-R_A|} | \chi_j \rangle`
|
||||
!
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
integer :: num_A, num_B, power_A(3), power_B(3)
|
||||
integer :: i, j, k, l, n_pt_in, m
|
||||
double precision :: c, Z, A_center(3), B_center(3), C_center(3)
|
||||
complex*16 :: alpha, beta, c1, c2
|
||||
|
||||
complex*16 :: NAI_pol_mult_cosgtos
|
||||
|
||||
ao_integrals_n_e_cosgtos = 0.d0
|
||||
|
||||
!$OMP PARALLEL &
|
||||
!$OMP DEFAULT (NONE) &
|
||||
!$OMP PRIVATE ( i, j, k, l, m, alpha, beta, A_center, B_center, C_center &
|
||||
!$OMP , power_A, power_B, num_A, num_B, Z, c, c1, c2, n_pt_in ) &
|
||||
!$OMP SHARED ( ao_num, ao_prim_num, ao_nucl, nucl_coord, ao_power, nucl_num, nucl_charge &
|
||||
!$OMP , ao_expo_ord_transp_cosgtos, ao_coef_norm_ord_transp_cosgtos &
|
||||
!$OMP , n_pt_max_integrals, ao_integrals_n_e_cosgtos )
|
||||
|
||||
n_pt_in = n_pt_max_integrals
|
||||
|
||||
!$OMP DO SCHEDULE (dynamic)
|
||||
|
||||
do j = 1, ao_num
|
||||
num_A = ao_nucl(j)
|
||||
power_A(1:3) = ao_power(j,1:3)
|
||||
A_center(1:3) = nucl_coord(num_A,1:3)
|
||||
|
||||
do i = 1, ao_num
|
||||
num_B = ao_nucl(i)
|
||||
power_B(1:3) = ao_power(i,1:3)
|
||||
B_center(1:3) = nucl_coord(num_B,1:3)
|
||||
|
||||
do l = 1, ao_prim_num(j)
|
||||
alpha = ao_expo_ord_transp_cosgtos(l,j)
|
||||
|
||||
do m = 1, ao_prim_num(i)
|
||||
beta = ao_expo_ord_transp_cosgtos(m,i)
|
||||
|
||||
c = 0.d0
|
||||
do k = 1, nucl_num
|
||||
|
||||
Z = nucl_charge(k)
|
||||
|
||||
C_center(1:3) = nucl_coord(k,1:3)
|
||||
|
||||
!print *, ' '
|
||||
!print *, A_center, B_center, C_center, power_A, power_B
|
||||
!print *, real(alpha), real(beta)
|
||||
|
||||
c1 = NAI_pol_mult_cosgtos( A_center, B_center, power_A, power_B &
|
||||
, alpha, beta, C_center, n_pt_in )
|
||||
|
||||
!c2 = c1
|
||||
c2 = NAI_pol_mult_cosgtos( A_center, B_center, power_A, power_B &
|
||||
, conjg(alpha), beta, C_center, n_pt_in )
|
||||
|
||||
!print *, ' c1 = ', real(c1)
|
||||
!print *, ' c2 = ', real(c2)
|
||||
|
||||
c = c - Z * 2.d0 * real(c1 + c2)
|
||||
|
||||
enddo
|
||||
ao_integrals_n_e_cosgtos(i,j) = ao_integrals_n_e_cosgtos(i,j) &
|
||||
+ ao_coef_norm_ord_transp_cosgtos(l,j) &
|
||||
* ao_coef_norm_ord_transp_cosgtos(m,i) * c
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
complex*16 function NAI_pol_mult_cosgtos(A_center, B_center, power_A, power_B, alpha, beta, C_center, n_pt_in)
|
||||
|
||||
BEGIN_DOC
|
||||
!
|
||||
! Computes the electron-nucleus attraction with two primitves cosgtos.
|
||||
!
|
||||
! :math:`\langle g_i | \frac{1}{|r-R_c|} | g_j \rangle`
|
||||
!
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
include 'utils/constants.include.F'
|
||||
|
||||
integer, intent(in) :: n_pt_in, power_A(3), power_B(3)
|
||||
double precision, intent(in) :: C_center(3), A_center(3), B_center(3)
|
||||
complex*16, intent(in) :: alpha, beta
|
||||
|
||||
integer :: i, n_pt, n_pt_out
|
||||
double precision :: dist, const_mod
|
||||
complex*16 :: p, p_inv, rho, dist_integral, const, const_factor, coeff, factor
|
||||
complex*16 :: accu, P_center(3)
|
||||
complex*16 :: d(0:n_pt_in)
|
||||
|
||||
complex*16 :: V_n_e_cosgtos
|
||||
complex*16 :: crint
|
||||
|
||||
if ( (A_center(1)/=B_center(1)) .or. (A_center(2)/=B_center(2)) .or. (A_center(3)/=B_center(3)) .or. &
|
||||
(A_center(1)/=C_center(1)) .or. (A_center(2)/=C_center(2)) .or. (A_center(3)/=C_center(3)) ) then
|
||||
|
||||
continue
|
||||
|
||||
else
|
||||
|
||||
NAI_pol_mult_cosgtos = V_n_e_cosgtos( power_A(1), power_A(2), power_A(3) &
|
||||
, power_B(1), power_B(2), power_B(3) &
|
||||
, alpha, beta )
|
||||
return
|
||||
|
||||
endif
|
||||
|
||||
p = alpha + beta
|
||||
p_inv = (1.d0, 0.d0) / p
|
||||
rho = alpha * beta * p_inv
|
||||
|
||||
dist = 0.d0
|
||||
dist_integral = (0.d0, 0.d0)
|
||||
do i = 1, 3
|
||||
P_center(i) = (alpha * A_center(i) + beta * B_center(i)) * p_inv
|
||||
dist += (A_center(i) - B_center(i)) * (A_center(i) - B_center(i))
|
||||
dist_integral += (P_center(i) - C_center(i)) * (P_center(i) - C_center(i))
|
||||
enddo
|
||||
|
||||
const_factor = dist * rho
|
||||
const = p * dist_integral
|
||||
|
||||
const_mod = dsqrt(real(const_factor)*real(const_factor) + aimag(const_factor)*aimag(const_factor))
|
||||
if(const_mod > 80.d0) then
|
||||
NAI_pol_mult_cosgtos = (0.d0, 0.d0)
|
||||
return
|
||||
endif
|
||||
|
||||
factor = zexp(-const_factor)
|
||||
coeff = dtwo_pi * factor * p_inv
|
||||
|
||||
do i = 0, n_pt_in
|
||||
d(i) = (0.d0, 0.d0)
|
||||
enddo
|
||||
|
||||
n_pt = 2 * ( (power_A(1) + power_B(1)) + (power_A(2) + power_B(2)) + (power_A(3) + power_B(3)) )
|
||||
if(n_pt == 0) then
|
||||
NAI_pol_mult_cosgtos = coeff * crint(0, const)
|
||||
return
|
||||
endif
|
||||
|
||||
call give_cpolynomial_mult_center_one_e( A_center, B_center, alpha, beta &
|
||||
, power_A, power_B, C_center, n_pt_in, d, n_pt_out)
|
||||
|
||||
if(n_pt_out < 0) then
|
||||
NAI_pol_mult_cosgtos = (0.d0, 0.d0)
|
||||
return
|
||||
endif
|
||||
|
||||
accu = (0.d0, 0.d0)
|
||||
do i = 0, n_pt_out, 2
|
||||
accu += crint(shiftr(i, 1), const) * d(i)
|
||||
|
||||
! print *, shiftr(i, 1), real(const), real(d(i)), real(crint(shiftr(i, 1), const))
|
||||
enddo
|
||||
NAI_pol_mult_cosgtos = accu * coeff
|
||||
|
||||
end function NAI_pol_mult_cosgtos
|
||||
|
||||
! ---
|
||||
|
||||
subroutine give_cpolynomial_mult_center_one_e( A_center, B_center, alpha, beta &
|
||||
, power_A, power_B, C_center, n_pt_in, d, n_pt_out)
|
||||
|
||||
BEGIN_DOC
|
||||
! Returns the explicit polynomial in terms of the "t" variable of the following
|
||||
!
|
||||
! $I_{x1}(a_x, d_x,p,q) \times I_{x1}(a_y, d_y,p,q) \times I_{x1}(a_z, d_z,p,q)$.
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
|
||||
integer, intent(in) :: n_pt_in, power_A(3), power_B(3)
|
||||
double precision, intent(in) :: A_center(3), B_center(3), C_center(3)
|
||||
complex*16, intent(in) :: alpha, beta
|
||||
integer, intent(out) :: n_pt_out
|
||||
complex*16, intent(out) :: d(0:n_pt_in)
|
||||
|
||||
integer :: a_x, b_x, a_y, b_y, a_z, b_z
|
||||
integer :: n_pt1, n_pt2, n_pt3, dim, i, n_pt_tmp
|
||||
complex*16 :: p, P_center(3), rho, p_inv, p_inv_2
|
||||
complex*16 :: R1x(0:2), B01(0:2), R1xp(0:2),R2x(0:2)
|
||||
complex*16 :: d1(0:n_pt_in), d2(0:n_pt_in), d3(0:n_pt_in)
|
||||
|
||||
ASSERT (n_pt_in > 1)
|
||||
|
||||
p = alpha + beta
|
||||
p_inv = (1.d0, 0.d0) / p
|
||||
p_inv_2 = 0.5d0 * p_inv
|
||||
|
||||
do i = 1, 3
|
||||
P_center(i) = (alpha * A_center(i) + beta * B_center(i)) * p_inv
|
||||
enddo
|
||||
|
||||
do i = 0, n_pt_in
|
||||
d(i) = (0.d0, 0.d0)
|
||||
d1(i) = (0.d0, 0.d0)
|
||||
d2(i) = (0.d0, 0.d0)
|
||||
d3(i) = (0.d0, 0.d0)
|
||||
enddo
|
||||
|
||||
! ---
|
||||
|
||||
n_pt1 = n_pt_in
|
||||
|
||||
R1x(0) = (P_center(1) - A_center(1))
|
||||
R1x(1) = (0.d0, 0.d0)
|
||||
R1x(2) = -(P_center(1) - C_center(1))
|
||||
|
||||
R1xp(0) = (P_center(1) - B_center(1))
|
||||
R1xp(1) = (0.d0, 0.d0)
|
||||
R1xp(2) = -(P_center(1) - C_center(1))
|
||||
|
||||
R2x(0) = p_inv_2
|
||||
R2x(1) = (0.d0, 0.d0)
|
||||
R2x(2) = -p_inv_2
|
||||
|
||||
a_x = power_A(1)
|
||||
b_x = power_B(1)
|
||||
call I_x1_pol_mult_one_e_cosgtos(a_x, b_x, R1x, R1xp, R2x, d1, n_pt1, n_pt_in)
|
||||
|
||||
if(n_pt1 < 0) then
|
||||
n_pt_out = -1
|
||||
do i = 0, n_pt_in
|
||||
d(i) = (0.d0, 0.d0)
|
||||
enddo
|
||||
return
|
||||
endif
|
||||
|
||||
! ---
|
||||
|
||||
n_pt2 = n_pt_in
|
||||
|
||||
R1x(0) = (P_center(2) - A_center(2))
|
||||
R1x(1) = (0.d0, 0.d0)
|
||||
R1x(2) = -(P_center(2) - C_center(2))
|
||||
|
||||
R1xp(0) = (P_center(2) - B_center(2))
|
||||
R1xp(1) = (0.d0, 0.d0)
|
||||
R1xp(2) = -(P_center(2) - C_center(2))
|
||||
|
||||
a_y = power_A(2)
|
||||
b_y = power_B(2)
|
||||
call I_x1_pol_mult_one_e_cosgtos(a_y, b_y, R1x, R1xp, R2x, d2, n_pt2, n_pt_in)
|
||||
|
||||
if(n_pt2 < 0) then
|
||||
n_pt_out = -1
|
||||
do i = 0, n_pt_in
|
||||
d(i) = (0.d0, 0.d0)
|
||||
enddo
|
||||
return
|
||||
endif
|
||||
|
||||
! ---
|
||||
|
||||
n_pt3 = n_pt_in
|
||||
|
||||
R1x(0) = (P_center(3) - A_center(3))
|
||||
R1x(1) = (0.d0, 0.d0)
|
||||
R1x(2) = -(P_center(3) - C_center(3))
|
||||
|
||||
R1xp(0) = (P_center(3) - B_center(3))
|
||||
R1xp(1) = (0.d0, 0.d0)
|
||||
R1xp(2) = -(P_center(3) - C_center(3))
|
||||
|
||||
a_z = power_A(3)
|
||||
b_z = power_B(3)
|
||||
call I_x1_pol_mult_one_e_cosgtos(a_z, b_z, R1x, R1xp, R2x, d3, n_pt3, n_pt_in)
|
||||
|
||||
if(n_pt3 < 0) then
|
||||
n_pt_out = -1
|
||||
do i = 0, n_pt_in
|
||||
d(i) = (0.d0, 0.d0)
|
||||
enddo
|
||||
return
|
||||
endif
|
||||
|
||||
! ---
|
||||
|
||||
n_pt_tmp = 0
|
||||
call multiply_cpoly(d1, n_pt1, d2, n_pt2, d, n_pt_tmp)
|
||||
do i = 0, n_pt_tmp
|
||||
d1(i) = (0.d0, 0.d0)
|
||||
enddo
|
||||
|
||||
n_pt_out = 0
|
||||
call multiply_cpoly(d, n_pt_tmp, d3, n_pt3, d1, n_pt_out)
|
||||
do i = 0, n_pt_out
|
||||
d(i) = d1(i)
|
||||
enddo
|
||||
|
||||
end subroutine give_cpolynomial_mult_center_one_e
|
||||
|
||||
! ---
|
||||
|
||||
recursive subroutine I_x1_pol_mult_one_e_cosgtos(a, c, R1x, R1xp, R2x, d, nd, n_pt_in)
|
||||
|
||||
BEGIN_DOC
|
||||
! Recursive routine involved in the electron-nucleus potential
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
include 'utils/constants.include.F'
|
||||
|
||||
integer, intent(in) :: a, c, n_pt_in
|
||||
complex*16, intent(in) :: R1x(0:2), R1xp(0:2), R2x(0:2)
|
||||
integer, intent(inout) :: nd
|
||||
complex*16, intent(inout) :: d(0:n_pt_in)
|
||||
|
||||
integer :: nx, ix, dim, iy, ny
|
||||
complex*16 :: X(0:max_dim)
|
||||
complex*16 :: Y(0:max_dim)
|
||||
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: X, Y
|
||||
|
||||
dim = n_pt_in
|
||||
|
||||
if( (a==0) .and. (c==0)) then
|
||||
|
||||
nd = 0
|
||||
d(0) = (1.d0, 0.d0)
|
||||
return
|
||||
|
||||
elseif( (c < 0) .or. (nd < 0) ) then
|
||||
|
||||
nd = -1
|
||||
return
|
||||
|
||||
elseif((a == 0) .and. (c .ne. 0)) then
|
||||
|
||||
call I_x2_pol_mult_one_e_cosgtos(c, R1x, R1xp, R2x, d, nd, n_pt_in)
|
||||
|
||||
elseif(a == 1) then
|
||||
|
||||
nx = nd
|
||||
do ix = 0, n_pt_in
|
||||
X(ix) = (0.d0, 0.d0)
|
||||
Y(ix) = (0.d0, 0.d0)
|
||||
enddo
|
||||
|
||||
call I_x2_pol_mult_one_e_cosgtos(c-1, R1x, R1xp, R2x, X, nx, n_pt_in)
|
||||
|
||||
do ix = 0, nx
|
||||
X(ix) *= dble(c)
|
||||
enddo
|
||||
|
||||
call multiply_cpoly(X, nx, R2x, 2, d, nd)
|
||||
|
||||
ny = 0
|
||||
call I_x2_pol_mult_one_e_cosgtos(c, R1x, R1xp, R2x, Y, ny, n_pt_in)
|
||||
call multiply_cpoly(Y, ny, R1x, 2, d, nd)
|
||||
|
||||
else
|
||||
|
||||
nx = 0
|
||||
do ix = 0, n_pt_in
|
||||
X(ix) = (0.d0, 0.d0)
|
||||
Y(ix) = (0.d0, 0.d0)
|
||||
enddo
|
||||
|
||||
call I_x1_pol_mult_one_e_cosgtos(a-2, c, R1x, R1xp, R2x, X, nx, n_pt_in)
|
||||
|
||||
do ix = 0, nx
|
||||
X(ix) *= dble(a-1)
|
||||
enddo
|
||||
call multiply_cpoly(X, nx, R2x, 2, d, nd)
|
||||
|
||||
nx = nd
|
||||
do ix = 0, n_pt_in
|
||||
X(ix) = (0.d0, 0.d0)
|
||||
enddo
|
||||
|
||||
call I_x1_pol_mult_one_e_cosgtos(a-1, c-1, R1x, R1xp, R2x, X, nx, n_pt_in)
|
||||
do ix = 0, nx
|
||||
X(ix) *= dble(c)
|
||||
enddo
|
||||
|
||||
call multiply_cpoly(X, nx, R2x, 2, d, nd)
|
||||
|
||||
ny = 0
|
||||
call I_x1_pol_mult_one_e_cosgtos(a-1, c, R1x, R1xp, R2x, Y, ny, n_pt_in)
|
||||
call multiply_cpoly(Y, ny, R1x, 2, d, nd)
|
||||
|
||||
endif
|
||||
|
||||
end subroutine I_x1_pol_mult_one_e_cosgtos
|
||||
|
||||
! ---
|
||||
|
||||
recursive subroutine I_x2_pol_mult_one_e_cosgtos(c, R1x, R1xp, R2x, d, nd, dim)
|
||||
|
||||
BEGIN_DOC
|
||||
! Recursive routine involved in the electron-nucleus potential
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
include 'utils/constants.include.F'
|
||||
|
||||
integer, intent(in) :: dim, c
|
||||
complex*16, intent(in) :: R1x(0:2), R1xp(0:2), R2x(0:2)
|
||||
integer, intent(inout) :: nd
|
||||
complex*16, intent(out) :: d(0:max_dim)
|
||||
|
||||
integer :: i, nx, ix, ny
|
||||
complex*16 :: X(0:max_dim), Y(0:max_dim)
|
||||
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: X, Y
|
||||
|
||||
if(c == 0) then
|
||||
|
||||
nd = 0
|
||||
d(0) = (1.d0, 0.d0)
|
||||
return
|
||||
|
||||
elseif((nd < 0) .or. (c < 0)) then
|
||||
|
||||
nd = -1
|
||||
return
|
||||
|
||||
else
|
||||
|
||||
nx = 0
|
||||
do ix = 0, dim
|
||||
X(ix) = (0.d0, 0.d0)
|
||||
Y(ix) = (0.d0, 0.d0)
|
||||
enddo
|
||||
|
||||
call I_x1_pol_mult_one_e_cosgtos(0, c-2, R1x, R1xp, R2x, X, nx, dim)
|
||||
|
||||
do ix = 0, nx
|
||||
X(ix) *= dble(c-1)
|
||||
enddo
|
||||
|
||||
call multiply_cpoly(X, nx, R2x, 2, d, nd)
|
||||
|
||||
ny = 0
|
||||
do ix = 0, dim
|
||||
Y(ix) = (0.d0, 0.d0)
|
||||
enddo
|
||||
|
||||
call I_x1_pol_mult_one_e_cosgtos(0, c-1, R1x, R1xp, R2x, Y, ny, dim)
|
||||
|
||||
if(ny .ge. 0) then
|
||||
call multiply_cpoly(Y, ny, R1xp, 2, d, nd)
|
||||
endif
|
||||
|
||||
endif
|
||||
|
||||
end subroutine I_x2_pol_mult_one_e_cosgtos
|
||||
|
||||
! ---
|
||||
|
||||
complex*16 function V_n_e_cosgtos(a_x, a_y, a_z, b_x, b_y, b_z, alpha, beta)
|
||||
|
||||
BEGIN_DOC
|
||||
! Primitve nuclear attraction between the two primitves centered on the same atom.
|
||||
!
|
||||
! $p_1 = x^{a_x} y^{a_y} z^{a_z} \exp(-\alpha r^2)$
|
||||
!
|
||||
! $p_2 = x^{b_x} y^{b_y} z^{b_z} \exp(-\beta r^2)$
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
|
||||
integer, intent(in) :: a_x, a_y, a_z, b_x, b_y, b_z
|
||||
complex*16, intent(in) :: alpha, beta
|
||||
|
||||
double precision :: V_phi, V_theta
|
||||
complex*16 :: V_r_cosgtos
|
||||
|
||||
if( (iand(a_x + b_x, 1) == 1) .or. &
|
||||
(iand(a_y + b_y, 1) == 1) .or. &
|
||||
(iand(a_z + b_z, 1) == 1) ) then
|
||||
|
||||
V_n_e_cosgtos = (0.d0, 0.d0)
|
||||
|
||||
else
|
||||
|
||||
V_n_e_cosgtos = V_r_cosgtos(a_x + b_x + a_y + b_y + a_z + b_z + 1, alpha + beta) &
|
||||
* V_phi(a_x + b_x, a_y + b_y) &
|
||||
* V_theta(a_z + b_z, a_x + b_x + a_y + b_y + 1)
|
||||
endif
|
||||
|
||||
end function V_n_e_cosgtos
|
||||
|
||||
! ---
|
||||
|
||||
complex*16 function V_r_cosgtos(n, alpha)
|
||||
|
||||
BEGIN_DOC
|
||||
! Computes the radial part of the nuclear attraction integral:
|
||||
!
|
||||
! $\int_{0}^{\infty} r^n \exp(-\alpha r^2) dr$
|
||||
!
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
include 'utils/constants.include.F'
|
||||
|
||||
integer , intent(in) :: n
|
||||
complex*16, intent(in) :: alpha
|
||||
|
||||
double precision :: fact
|
||||
|
||||
if(iand(n, 1) .eq. 1) then
|
||||
V_r_cosgtos = 0.5d0 * fact(shiftr(n, 1)) / (alpha**(shiftr(n, 1) + 1))
|
||||
else
|
||||
V_r_cosgtos = sqpi * fact(n) / fact(shiftr(n, 1)) * (0.5d0/zsqrt(alpha))**(n+1)
|
||||
endif
|
||||
|
||||
end function V_r_cosgtos
|
||||
|
||||
! ---
|
||||
|
223
src/cosgtos_ao_int/one_e_kin_integrals.irp.f
Normal file
223
src/cosgtos_ao_int/one_e_kin_integrals.irp.f
Normal file
@ -0,0 +1,223 @@
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [ double precision, ao_deriv2_cosgtos_x, (ao_num, ao_num) ]
|
||||
&BEGIN_PROVIDER [ double precision, ao_deriv2_cosgtos_y, (ao_num, ao_num) ]
|
||||
&BEGIN_PROVIDER [ double precision, ao_deriv2_cosgtos_z, (ao_num, ao_num) ]
|
||||
|
||||
implicit none
|
||||
integer :: i, j, n, l, dim1, power_A(3), power_B(3)
|
||||
double precision :: c, deriv_tmp
|
||||
complex*16 :: alpha, beta, A_center(3), B_center(3)
|
||||
complex*16 :: overlap_x, overlap_y, overlap_z, overlap
|
||||
complex*16 :: overlap_x0_1, overlap_y0_1, overlap_z0_1
|
||||
complex*16 :: overlap_x0_2, overlap_y0_2, overlap_z0_2
|
||||
complex*16 :: overlap_m2_1, overlap_p2_1
|
||||
complex*16 :: overlap_m2_2, overlap_p2_2
|
||||
complex*16 :: deriv_tmp_1, deriv_tmp_2
|
||||
|
||||
|
||||
dim1 = 100
|
||||
|
||||
! -- Dummy call to provide everything
|
||||
|
||||
A_center(:) = (0.0d0, 0.d0)
|
||||
B_center(:) = (1.0d0, 0.d0)
|
||||
alpha = (1.0d0, 0.d0)
|
||||
beta = (0.1d0, 0.d0)
|
||||
power_A = 1
|
||||
power_B = 0
|
||||
call overlap_cgaussian_xyz( A_center, B_center, alpha, beta, power_A, power_B &
|
||||
, overlap_x0_1, overlap_y0_1, overlap_z0_1, overlap, dim1 )
|
||||
|
||||
! ---
|
||||
|
||||
!$OMP PARALLEL DO SCHEDULE(GUIDED) &
|
||||
!$OMP DEFAULT(NONE) &
|
||||
!$OMP PRIVATE( A_center, B_center, power_A, power_B, alpha, beta, i, j, l, n, c &
|
||||
!$OMP , deriv_tmp, deriv_tmp_1, deriv_tmp_2 &
|
||||
!$OMP , overlap_x, overlap_y, overlap_z, overlap &
|
||||
!$OMP , overlap_m2_1, overlap_p2_1, overlap_m2_2, overlap_p2_2 &
|
||||
!$OMP , overlap_x0_1, overlap_y0_1, overlap_z0_1, overlap_x0_2, overlap_y0_2, overlap_z0_2 ) &
|
||||
!$OMP SHARED( nucl_coord, ao_power, ao_prim_num, ao_num, ao_nucl, dim1 &
|
||||
!$OMP , ao_coef_norm_ord_transp_cosgtos, ao_expo_ord_transp_cosgtos &
|
||||
!$OMP , ao_deriv2_cosgtos_x, ao_deriv2_cosgtos_y, ao_deriv2_cosgtos_z )
|
||||
|
||||
do j = 1, ao_num
|
||||
A_center(1) = nucl_coord(ao_nucl(j),1) * (1.d0, 0.d0)
|
||||
A_center(2) = nucl_coord(ao_nucl(j),2) * (1.d0, 0.d0)
|
||||
A_center(3) = nucl_coord(ao_nucl(j),3) * (1.d0, 0.d0)
|
||||
power_A(1) = ao_power(j,1)
|
||||
power_A(2) = ao_power(j,2)
|
||||
power_A(3) = ao_power(j,3)
|
||||
|
||||
do i = 1, ao_num
|
||||
B_center(1) = nucl_coord(ao_nucl(i),1) * (1.d0, 0.d0)
|
||||
B_center(2) = nucl_coord(ao_nucl(i),2) * (1.d0, 0.d0)
|
||||
B_center(3) = nucl_coord(ao_nucl(i),3) * (1.d0, 0.d0)
|
||||
power_B(1) = ao_power(i,1)
|
||||
power_B(2) = ao_power(i,2)
|
||||
power_B(3) = ao_power(i,3)
|
||||
|
||||
ao_deriv2_cosgtos_x(i,j) = 0.d0
|
||||
ao_deriv2_cosgtos_y(i,j) = 0.d0
|
||||
ao_deriv2_cosgtos_z(i,j) = 0.d0
|
||||
|
||||
do n = 1, ao_prim_num(j)
|
||||
alpha = ao_expo_ord_transp_cosgtos(n,j)
|
||||
|
||||
do l = 1, ao_prim_num(i)
|
||||
c = ao_coef_norm_ord_transp_cosgtos(n,j) * ao_coef_norm_ord_transp_cosgtos(l,i)
|
||||
beta = ao_expo_ord_transp_cosgtos(l,i)
|
||||
|
||||
call overlap_cgaussian_xyz( A_center, B_center, alpha, beta, power_A, power_B &
|
||||
, overlap_x0_1, overlap_y0_1, overlap_z0_1, overlap, dim1 )
|
||||
|
||||
call overlap_cgaussian_xyz( A_center, B_center, alpha, conjg(beta), power_A, power_B &
|
||||
, overlap_x0_2, overlap_y0_2, overlap_z0_2, overlap, dim1 )
|
||||
|
||||
! ---
|
||||
|
||||
power_A(1) = power_A(1) - 2
|
||||
if(power_A(1) > -1) then
|
||||
call overlap_cgaussian_xyz( A_center, B_center, alpha, beta, power_A, power_B &
|
||||
, overlap_m2_1, overlap_y, overlap_z, overlap, dim1 )
|
||||
|
||||
call overlap_cgaussian_xyz( A_center, B_center, alpha, conjg(beta), power_A, power_B &
|
||||
, overlap_m2_2, overlap_y, overlap_z, overlap, dim1 )
|
||||
else
|
||||
overlap_m2_1 = (0.d0, 0.d0)
|
||||
overlap_m2_2 = (0.d0, 0.d0)
|
||||
endif
|
||||
|
||||
power_A(1) = power_A(1) + 4
|
||||
call overlap_cgaussian_xyz( A_center, B_center, alpha, beta, power_A, power_B &
|
||||
, overlap_p2_1, overlap_y, overlap_z, overlap, dim1 )
|
||||
|
||||
call overlap_cgaussian_xyz( A_center, B_center, alpha, conjg(beta), power_A, power_B &
|
||||
, overlap_p2_2, overlap_y, overlap_z, overlap, dim1 )
|
||||
|
||||
power_A(1) = power_A(1) - 2
|
||||
|
||||
deriv_tmp_1 = ( -2.d0 * alpha * (2.d0 * power_A(1) + 1.d0) * overlap_x0_1 &
|
||||
+ power_A(1) * (power_A(1) - 1.d0) * overlap_m2_1 &
|
||||
+ 4.d0 * alpha * alpha * overlap_p2_1 ) * overlap_y0_1 * overlap_z0_1
|
||||
|
||||
deriv_tmp_2 = ( -2.d0 * alpha * (2.d0 * power_A(1) + 1.d0) * overlap_x0_2 &
|
||||
+ power_A(1) * (power_A(1) - 1.d0) * overlap_m2_2 &
|
||||
+ 4.d0 * alpha * alpha * overlap_p2_2 ) * overlap_y0_2 * overlap_z0_2
|
||||
|
||||
deriv_tmp = 2.d0 * real(deriv_tmp_1 + deriv_tmp_2)
|
||||
|
||||
ao_deriv2_cosgtos_x(i,j) += c * deriv_tmp
|
||||
|
||||
! ---
|
||||
|
||||
power_A(2) = power_A(2) - 2
|
||||
if(power_A(2) > -1) then
|
||||
call overlap_cgaussian_xyz( A_center, B_center, alpha, beta, power_A, power_B &
|
||||
, overlap_x, overlap_m2_1, overlap_y, overlap, dim1 )
|
||||
|
||||
call overlap_cgaussian_xyz( A_center, B_center, alpha, conjg(beta), power_A, power_B &
|
||||
, overlap_x, overlap_m2_2, overlap_y, overlap, dim1 )
|
||||
else
|
||||
overlap_m2_1 = (0.d0, 0.d0)
|
||||
overlap_m2_2 = (0.d0, 0.d0)
|
||||
endif
|
||||
|
||||
power_A(2) = power_A(2) + 4
|
||||
call overlap_cgaussian_xyz( A_center, B_center, alpha, beta, power_A, power_B &
|
||||
, overlap_x, overlap_p2_1, overlap_y, overlap, dim1 )
|
||||
|
||||
call overlap_cgaussian_xyz( A_center, B_center, alpha, conjg(beta), power_A, power_B &
|
||||
, overlap_x, overlap_p2_2, overlap_y, overlap, dim1 )
|
||||
|
||||
power_A(2) = power_A(2) - 2
|
||||
|
||||
deriv_tmp_1 = ( -2.d0 * alpha * (2.d0 * power_A(2) + 1.d0) * overlap_y0_1 &
|
||||
+ power_A(2) * (power_A(2) - 1.d0) * overlap_m2_1 &
|
||||
+ 4.d0 * alpha * alpha * overlap_p2_1 ) * overlap_x0_1 * overlap_z0_1
|
||||
|
||||
deriv_tmp_2 = ( -2.d0 * alpha * (2.d0 * power_A(2) + 1.d0) * overlap_y0_2 &
|
||||
+ power_A(2) * (power_A(2) - 1.d0) * overlap_m2_2 &
|
||||
+ 4.d0 * alpha * alpha * overlap_p2_2 ) * overlap_x0_2 * overlap_z0_2
|
||||
|
||||
deriv_tmp = 2.d0 * real(deriv_tmp_1 + deriv_tmp_2)
|
||||
|
||||
ao_deriv2_cosgtos_y(i,j) += c * deriv_tmp
|
||||
|
||||
! ---
|
||||
|
||||
power_A(3) = power_A(3) - 2
|
||||
if(power_A(3) > -1) then
|
||||
call overlap_cgaussian_xyz( A_center, B_center, alpha, beta, power_A, power_B &
|
||||
, overlap_x, overlap_y, overlap_m2_1, overlap, dim1 )
|
||||
|
||||
call overlap_cgaussian_xyz( A_center, B_center, alpha, conjg(beta), power_A, power_B &
|
||||
, overlap_x, overlap_y, overlap_m2_2, overlap, dim1 )
|
||||
else
|
||||
overlap_m2_1 = (0.d0, 0.d0)
|
||||
overlap_m2_2 = (0.d0, 0.d0)
|
||||
endif
|
||||
|
||||
power_A(3) = power_A(3) + 4
|
||||
call overlap_cgaussian_xyz( A_center, B_center, alpha, beta, power_A, power_B &
|
||||
, overlap_x, overlap_y, overlap_p2_1, overlap, dim1 )
|
||||
|
||||
call overlap_cgaussian_xyz( A_center, B_center, alpha, conjg(beta), power_A, power_B &
|
||||
, overlap_x, overlap_y, overlap_p2_2, overlap, dim1 )
|
||||
|
||||
power_A(3) = power_A(3) - 2
|
||||
|
||||
deriv_tmp_1 = ( -2.d0 * alpha * (2.d0 * power_A(3) + 1.d0) * overlap_z0_1 &
|
||||
+ power_A(3) * (power_A(3) - 1.d0) * overlap_m2_1 &
|
||||
+ 4.d0 * alpha * alpha * overlap_p2_1 ) * overlap_x0_1 * overlap_y0_1
|
||||
|
||||
deriv_tmp_2 = ( -2.d0 * alpha * (2.d0 * power_A(3) + 1.d0) * overlap_z0_2 &
|
||||
+ power_A(3) * (power_A(3) - 1.d0) * overlap_m2_2 &
|
||||
+ 4.d0 * alpha * alpha * overlap_p2_2 ) * overlap_x0_2 * overlap_y0_2
|
||||
|
||||
deriv_tmp = 2.d0 * real(deriv_tmp_1 + deriv_tmp_2)
|
||||
|
||||
ao_deriv2_cosgtos_z(i,j) += c * deriv_tmp
|
||||
|
||||
! ---
|
||||
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END PARALLEL DO
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [double precision, ao_kinetic_integrals_cosgtos, (ao_num, ao_num)]
|
||||
|
||||
BEGIN_DOC
|
||||
!
|
||||
! Kinetic energy integrals in the cosgtos |AO| basis.
|
||||
!
|
||||
! $\langle \chi_i |\hat{T}| \chi_j \rangle$
|
||||
!
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
integer :: i, j
|
||||
|
||||
!$OMP PARALLEL DO DEFAULT(NONE) &
|
||||
!$OMP PRIVATE(i, j) &
|
||||
!$OMP SHARED(ao_num, ao_kinetic_integrals_cosgtos, ao_deriv2_cosgtos_x, ao_deriv2_cosgtos_y, ao_deriv2_cosgtos_z)
|
||||
do j = 1, ao_num
|
||||
do i = 1, ao_num
|
||||
ao_kinetic_integrals_cosgtos(i,j) = -0.5d0 * ( ao_deriv2_cosgtos_x(i,j) &
|
||||
+ ao_deriv2_cosgtos_y(i,j) &
|
||||
+ ao_deriv2_cosgtos_z(i,j) )
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END PARALLEL DO
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
1584
src/cosgtos_ao_int/two_e_Coul_integrals.irp.f
Normal file
1584
src/cosgtos_ao_int/two_e_Coul_integrals.irp.f
Normal file
File diff suppressed because it is too large
Load Diff
120
src/utils/cgtos_one_e.irp.f
Normal file
120
src/utils/cgtos_one_e.irp.f
Normal file
@ -0,0 +1,120 @@
|
||||
|
||||
! ---
|
||||
|
||||
complex*16 function overlap_cgaussian_x(A_center, B_center, alpha, beta, power_A, power_B, dim)
|
||||
|
||||
BEGIN_DOC
|
||||
!
|
||||
! \int_{-infty}^{+infty} (x-A_x)^ax (x-B_x)^bx exp(-alpha (x-A_x)^2) exp(- beta(x-B_X)^2) dx
|
||||
! with complex arguments
|
||||
!
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
include 'constants.include.F'
|
||||
|
||||
integer, intent(in) :: dim, power_A, power_B
|
||||
complex*16, intent(in) :: A_center, B_center, alpha, beta
|
||||
|
||||
integer :: i, iorder_p
|
||||
double precision :: fact_p_mod
|
||||
complex*16 :: P_new(0:max_dim), P_center, fact_p, p, inv_sq_p
|
||||
|
||||
complex*16 :: Fc_integral
|
||||
|
||||
|
||||
call give_explicit_cpoly_and_cgaussian_x( P_new, P_center, p, fact_p, iorder_p &
|
||||
, alpha, beta, power_A, power_B, A_center, B_center, dim)
|
||||
|
||||
fact_p_mod = dsqrt(real(fact_p)*real(fact_p) + aimag(fact_p)*aimag(fact_p))
|
||||
if(fact_p_mod .lt. 1.d-14) then
|
||||
overlap_cgaussian_x = (0.d0, 0.d0)
|
||||
return
|
||||
endif
|
||||
|
||||
|
||||
inv_sq_p = (1.d0, 0.d0) / zsqrt(p)
|
||||
|
||||
overlap_cgaussian_x = (0.d0, 0.d0)
|
||||
do i = 0, iorder_p
|
||||
overlap_cgaussian_x += P_new(i) * Fc_integral(i, inv_sq_p)
|
||||
enddo
|
||||
|
||||
overlap_cgaussian_x *= fact_p
|
||||
|
||||
end function overlap_cgaussian_x
|
||||
|
||||
! ---
|
||||
|
||||
subroutine overlap_cgaussian_xyz( A_center, B_center, alpha, beta, power_A, power_B &
|
||||
, overlap_x, overlap_y, overlap_z, overlap, dim )
|
||||
|
||||
BEGIN_DOC
|
||||
!
|
||||
! S_x = \int (x-A_x)^{a_x} exp(-\alpha(x-A_x)^2) (x-B_x)^{b_x} exp(-beta(x-B_x)^2) dx
|
||||
! S = S_x S_y S_z
|
||||
! for complex arguments
|
||||
!
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
include 'constants.include.F'
|
||||
|
||||
integer, intent(in) :: dim, power_A(3), power_B(3)
|
||||
complex*16, intent(in) :: A_center(3), B_center(3), alpha, beta
|
||||
complex*16, intent(out) :: overlap_x, overlap_y, overlap_z, overlap
|
||||
|
||||
integer :: i, nmax, iorder_p(3)
|
||||
double precision :: fact_p_mod
|
||||
complex*16 :: P_new(0:max_dim,3), P_center(3), fact_p, p, inv_sq_p
|
||||
complex*16 :: F_integral_tab(0:max_dim)
|
||||
|
||||
complex*16 :: Fc_integral
|
||||
|
||||
call give_explicit_cpoly_and_cgaussian(P_new, P_center, p, fact_p, iorder_p, alpha, beta, power_A, power_B, A_center, B_center, dim)
|
||||
|
||||
fact_p_mod = dsqrt(real(fact_p)*real(fact_p) + aimag(fact_p)*aimag(fact_p))
|
||||
if(fact_p_mod .lt. 1.d-14) then
|
||||
overlap_x = (1.d-10, 0.d0)
|
||||
overlap_y = (1.d-10, 0.d0)
|
||||
overlap_z = (1.d-10, 0.d0)
|
||||
overlap = (1.d-10, 0.d0)
|
||||
return
|
||||
endif
|
||||
|
||||
nmax = maxval(iorder_p)
|
||||
|
||||
inv_sq_p = (1.d0, 0.d0) / zsqrt(p)
|
||||
do i = 0, nmax
|
||||
F_integral_tab(i) = Fc_integral(i, inv_sq_p)
|
||||
enddo
|
||||
|
||||
overlap_x = P_new(0,1) * F_integral_tab(0)
|
||||
overlap_y = P_new(0,2) * F_integral_tab(0)
|
||||
overlap_z = P_new(0,3) * F_integral_tab(0)
|
||||
|
||||
do i = 1, iorder_p(1)
|
||||
overlap_x = overlap_x + P_new(i,1) * F_integral_tab(i)
|
||||
enddo
|
||||
call cgaussian_product_x(alpha, A_center(1), beta, B_center(1), fact_p, p, P_center(1))
|
||||
overlap_x *= fact_p
|
||||
|
||||
do i = 1, iorder_p(2)
|
||||
overlap_y = overlap_y + P_new(i,2) * F_integral_tab(i)
|
||||
enddo
|
||||
call cgaussian_product_x(alpha, A_center(2), beta, B_center(2), fact_p, p, P_center(2))
|
||||
overlap_y *= fact_p
|
||||
|
||||
do i = 1, iorder_p(3)
|
||||
overlap_z = overlap_z + P_new(i,3) * F_integral_tab(i)
|
||||
enddo
|
||||
call cgaussian_product_x(alpha, A_center(3), beta, B_center(3), fact_p, p, P_center(3))
|
||||
overlap_z *= fact_p
|
||||
|
||||
overlap = overlap_x * overlap_y * overlap_z
|
||||
|
||||
end subroutine overlap_cgaussian_xyz
|
||||
|
||||
! ---
|
||||
|
||||
|
780
src/utils/cgtos_utils.irp.f
Normal file
780
src/utils/cgtos_utils.irp.f
Normal file
@ -0,0 +1,780 @@
|
||||
|
||||
! ---
|
||||
|
||||
subroutine give_explicit_cpoly_and_cgaussian_x(P_new, P_center, p, fact_k, iorder, alpha, beta, a, b, A_center, B_center, dim)
|
||||
|
||||
BEGIN_DOC
|
||||
! Transform the product of
|
||||
! (x-x_A)^a (x-x_B)^b exp(-(r-A)^2 alpha) exp(-(r-B)^2 beta)
|
||||
! into
|
||||
! fact_k \sum_{i=0}^{iorder} (x-x_P)^i exp(-p(r-P)^2)
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
include 'constants.include.F'
|
||||
|
||||
integer, intent(in) :: dim
|
||||
integer, intent(in) :: a, b
|
||||
complex*16, intent(in) :: alpha, beta, A_center, B_center
|
||||
integer, intent(out) :: iorder
|
||||
complex*16, intent(out) :: p, P_center, fact_k
|
||||
complex*16, intent(out) :: P_new(0:max_dim)
|
||||
|
||||
integer :: n_new, i, j
|
||||
double precision :: tmp_mod
|
||||
complex*16 :: P_a(0:max_dim), P_b(0:max_dim)
|
||||
complex*16 :: p_inv, ab, d_AB, tmp
|
||||
|
||||
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: P_a, P_b
|
||||
|
||||
P_new = (0.d0, 0.d0)
|
||||
|
||||
! new exponent
|
||||
p = alpha + beta
|
||||
|
||||
! new center
|
||||
p_inv = (1.d0, 0.d0) / p
|
||||
ab = alpha * beta
|
||||
P_center = (alpha * A_center + beta * B_center) * p_inv
|
||||
|
||||
! get the factor
|
||||
d_AB = (A_center - B_center) * (A_center - B_center)
|
||||
tmp = ab * p_inv * d_AB
|
||||
tmp_mod = dsqrt(REAL(tmp)*REAL(tmp) + AIMAG(tmp)*AIMAG(tmp))
|
||||
if(tmp_mod .lt. 50.d0) then
|
||||
fact_k = zexp(-tmp)
|
||||
else
|
||||
fact_k = (0.d0, 0.d0)
|
||||
endif
|
||||
|
||||
! Recenter the polynomials P_a and P_b on P_center
|
||||
!DIR$ FORCEINLINE
|
||||
call recentered_cpoly2(P_a(0), A_center, P_center, a, P_b(0), B_center, P_center, b)
|
||||
n_new = 0
|
||||
|
||||
!DIR$ FORCEINLINE
|
||||
call multiply_cpoly(P_a(0), a, P_b(0), b, P_new(0), n_new)
|
||||
iorder = a + b
|
||||
|
||||
end subroutine give_explicit_cpoly_and_cgaussian_x
|
||||
|
||||
! ---
|
||||
|
||||
subroutine give_explicit_cpoly_and_cgaussian(P_new, P_center, p, fact_k, iorder, alpha, beta, a, b, A_center, B_center, dim)
|
||||
|
||||
BEGIN_DOC
|
||||
! Transforms the product of
|
||||
! (x-x_A)^a(1) (x-x_B)^b(1) (y-y_A)^a(2) (y-y_B)^b(2) (z-z_A)^a(3) (z-z_B)^b(3) exp(-(r-A)^2 alpha) exp(-(r-B)^2 beta)
|
||||
! into
|
||||
! fact_k * [ sum (l_x = 0,i_order(1)) P_new(l_x,1) * (x-P_center(1))^l_x ] exp (- p (x-P_center(1))^2 )
|
||||
! * [ sum (l_y = 0,i_order(2)) P_new(l_y,2) * (y-P_center(2))^l_y ] exp (- p (y-P_center(2))^2 )
|
||||
! * [ sum (l_z = 0,i_order(3)) P_new(l_z,3) * (z-P_center(3))^l_z ] exp (- p (z-P_center(3))^2 )
|
||||
!
|
||||
! WARNING ::: IF fact_k is too smal then:
|
||||
! returns a "s" function centered in zero
|
||||
! with an inifinite exponent and a zero polynom coef
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
include 'constants.include.F'
|
||||
|
||||
integer, intent(in) :: dim, a(3), b(3)
|
||||
complex*16, intent(in) :: alpha, beta, A_center(3), B_center(3)
|
||||
integer, intent(out) :: iorder(3)
|
||||
complex*16, intent(out) :: p, P_center(3), fact_k, P_new(0:max_dim,3)
|
||||
|
||||
integer :: n_new, i, j
|
||||
double precision :: tmp_mod
|
||||
complex*16 :: P_a(0:max_dim,3), P_b(0:max_dim,3)
|
||||
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: P_a, P_b
|
||||
|
||||
iorder(1) = 0
|
||||
iorder(2) = 0
|
||||
iorder(3) = 0
|
||||
|
||||
P_new(0,1) = (0.d0, 0.d0)
|
||||
P_new(0,2) = (0.d0, 0.d0)
|
||||
P_new(0,3) = (0.d0, 0.d0)
|
||||
|
||||
!DIR$ FORCEINLINE
|
||||
call cgaussian_product(alpha, A_center, beta, B_center, fact_k, p, P_center)
|
||||
|
||||
! IF fact_k is too smal then: returns a "s" function centered in zero
|
||||
! with an inifinite exponent and a zero polynom coef
|
||||
tmp_mod = dsqrt(REAL(fact_k)*REAL(fact_k) + AIMAG(fact_k)*AIMAG(fact_k))
|
||||
if(tmp_mod < 1d-14) then
|
||||
iorder = 0
|
||||
p = (1.d+14, 0.d0)
|
||||
fact_k = (0.d0 , 0.d0)
|
||||
P_new(0:max_dim,1:3) = (0.d0 , 0.d0)
|
||||
P_center(1:3) = (0.d0 , 0.d0)
|
||||
return
|
||||
endif
|
||||
|
||||
!DIR$ FORCEINLINE
|
||||
call recentered_cpoly2(P_a(0,1), A_center(1), P_center(1), a(1), P_b(0,1), B_center(1), P_center(1), b(1))
|
||||
iorder(1) = a(1) + b(1)
|
||||
do i = 0, iorder(1)
|
||||
P_new(i,1) = 0.d0
|
||||
enddo
|
||||
n_new = 0
|
||||
!DIR$ FORCEINLINE
|
||||
call multiply_cpoly(P_a(0,1), a(1), P_b(0,1), b(1), P_new(0,1), n_new)
|
||||
|
||||
!DIR$ FORCEINLINE
|
||||
call recentered_cpoly2(P_a(0,2), A_center(2), P_center(2), a(2), P_b(0,2), B_center(2), P_center(2), b(2))
|
||||
iorder(2) = a(2) + b(2)
|
||||
do i = 0, iorder(2)
|
||||
P_new(i,2) = 0.d0
|
||||
enddo
|
||||
n_new = 0
|
||||
!DIR$ FORCEINLINE
|
||||
call multiply_cpoly(P_a(0,2), a(2), P_b(0,2), b(2), P_new(0,2), n_new)
|
||||
|
||||
!DIR$ FORCEINLINE
|
||||
call recentered_cpoly2(P_a(0,3), A_center(3), P_center(3), a(3), P_b(0,3), B_center(3), P_center(3), b(3))
|
||||
iorder(3) = a(3) + b(3)
|
||||
do i = 0, iorder(3)
|
||||
P_new(i,3) = 0.d0
|
||||
enddo
|
||||
n_new = 0
|
||||
!DIR$ FORCEINLINE
|
||||
call multiply_cpoly(P_a(0,3), a(3), P_b(0,3), b(3), P_new(0,3), n_new)
|
||||
|
||||
end subroutine give_explicit_cpoly_and_cgaussian
|
||||
|
||||
! ---
|
||||
|
||||
!subroutine give_explicit_poly_and_gaussian_double(P_new,P_center,p,fact_k,iorder,alpha,beta,gama,a,b,A_center,B_center,Nucl_center,dim)
|
||||
! BEGIN_DOC
|
||||
! ! Transforms the product of
|
||||
! ! (x-x_A)^a(1) (x-x_B)^b(1) (x-x_A)^a(2) (y-y_B)^b(2) (z-z_A)^a(3) (z-z_B)^b(3)
|
||||
! ! exp(-(r-A)^2 alpha) exp(-(r-B)^2 beta) exp(-(r-Nucl_center)^2 gama
|
||||
! !
|
||||
! ! into
|
||||
! ! fact_k * [ sum (l_x = 0,i_order(1)) P_new(l_x,1) * (x-P_center(1))^l_x ] exp (- p (x-P_center(1))^2 )
|
||||
! ! * [ sum (l_y = 0,i_order(2)) P_new(l_y,2) * (y-P_center(2))^l_y ] exp (- p (y-P_center(2))^2 )
|
||||
! ! * [ sum (l_z = 0,i_order(3)) P_new(l_z,3) * (z-P_center(3))^l_z ] exp (- p (z-P_center(3))^2 )
|
||||
! END_DOC
|
||||
! implicit none
|
||||
! include 'constants.include.F'
|
||||
! integer, intent(in) :: dim
|
||||
! integer, intent(in) :: a(3),b(3) ! powers : (x-xa)**a_x = (x-A(1))**a(1)
|
||||
! double precision, intent(in) :: alpha, beta, gama ! exponents
|
||||
! double precision, intent(in) :: A_center(3) ! A center
|
||||
! double precision, intent(in) :: B_center (3) ! B center
|
||||
! double precision, intent(in) :: Nucl_center(3) ! B center
|
||||
! double precision, intent(out) :: P_center(3) ! new center
|
||||
! double precision, intent(out) :: p ! new exponent
|
||||
! double precision, intent(out) :: fact_k ! constant factor
|
||||
! double precision, intent(out) :: P_new(0:max_dim,3)! polynomial
|
||||
! integer , intent(out) :: iorder(3) ! i_order(i) = order of the polynomials
|
||||
!
|
||||
! double precision :: P_center_tmp(3) ! new center
|
||||
! double precision :: p_tmp ! new exponent
|
||||
! double precision :: fact_k_tmp,fact_k_bis ! constant factor
|
||||
! double precision :: P_new_tmp(0:max_dim,3)! polynomial
|
||||
! integer :: i,j
|
||||
! double precision :: binom_func
|
||||
!
|
||||
! ! First you transform the two primitives into a sum of primitive with the same center P_center_tmp and gaussian exponent p_tmp
|
||||
! call give_explicit_cpoly_and_cgaussian(P_new_tmp,P_center_tmp,p_tmp,fact_k_tmp,iorder,alpha,beta,a,b,A_center,B_center,dim)
|
||||
! ! Then you create the new gaussian from the product of the new one per the Nuclei one
|
||||
! call cgaussian_product(p_tmp,P_center_tmp,gama,Nucl_center,fact_k_bis,p,P_center)
|
||||
! fact_k = fact_k_bis * fact_k_tmp
|
||||
!
|
||||
! ! Then you build the coefficient of the new polynom
|
||||
! do i = 0, iorder(1)
|
||||
! P_new(i,1) = 0.d0
|
||||
! do j = i,iorder(1)
|
||||
! P_new(i,1) = P_new(i,1) + P_new_tmp(j,1) * binom_func(j,j-i) * (P_center(1) - P_center_tmp(1))**(j-i)
|
||||
! enddo
|
||||
! enddo
|
||||
! do i = 0, iorder(2)
|
||||
! P_new(i,2) = 0.d0
|
||||
! do j = i,iorder(2)
|
||||
! P_new(i,2) = P_new(i,2) + P_new_tmp(j,2) * binom_func(j,j-i) * (P_center(2) - P_center_tmp(2))**(j-i)
|
||||
! enddo
|
||||
! enddo
|
||||
! do i = 0, iorder(3)
|
||||
! P_new(i,3) = 0.d0
|
||||
! do j = i,iorder(3)
|
||||
! P_new(i,3) = P_new(i,3) + P_new_tmp(j,3) * binom_func(j,j-i) * (P_center(3) - P_center_tmp(3))**(j-i)
|
||||
! enddo
|
||||
! enddo
|
||||
!
|
||||
!end
|
||||
|
||||
! ---
|
||||
|
||||
subroutine cgaussian_product(a, xa, b, xb, k, p, xp)
|
||||
|
||||
BEGIN_DOC
|
||||
! complex Gaussian product
|
||||
! e^{-a (r-r_A)^2} e^{-b (r-r_B)^2} = k e^{-p (r-r_P)^2}
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
complex*16, intent(in) :: a, b, xa(3), xb(3)
|
||||
complex*16, intent(out) :: p, k, xp(3)
|
||||
|
||||
double precision :: tmp_mod
|
||||
complex*16 :: p_inv, xab(3), ab
|
||||
|
||||
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: xab
|
||||
|
||||
ASSERT (REAL(a) > 0.)
|
||||
ASSERT (REAL(b) > 0.)
|
||||
|
||||
! new exponent
|
||||
p = a + b
|
||||
|
||||
xab(1) = xa(1) - xb(1)
|
||||
xab(2) = xa(2) - xb(2)
|
||||
xab(3) = xa(3) - xb(3)
|
||||
|
||||
p_inv = (1.d0, 0.d0) / p
|
||||
ab = a * b * p_inv
|
||||
|
||||
k = ab * (xab(1)*xab(1) + xab(2)*xab(2) + xab(3)*xab(3))
|
||||
tmp_mod = dsqrt(REAL(k)*REAL(k) + AIMAG(k)*AIMAG(k))
|
||||
if(tmp_mod .gt. 40.d0) then
|
||||
k = (0.d0, 0.d0)
|
||||
xp(1:3) = (0.d0, 0.d0)
|
||||
return
|
||||
endif
|
||||
|
||||
k = zexp(-k)
|
||||
xp(1) = ( a * xa(1) + b * xb(1) ) * p_inv
|
||||
xp(2) = ( a * xa(2) + b * xb(2) ) * p_inv
|
||||
xp(3) = ( a * xa(3) + b * xb(3) ) * p_inv
|
||||
|
||||
end subroutine cgaussian_product
|
||||
|
||||
! ---
|
||||
|
||||
subroutine cgaussian_product_x(a, xa, b, xb, k, p, xp)
|
||||
|
||||
BEGIN_DOC
|
||||
! complex Gaussian product in 1D.
|
||||
! e^{-a (x-x_A)^2} e^{-b (x-x_B)^2} = K e^{-p (x-x_P)^2}
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
complex*16, intent(in) :: a, b, xa, xb
|
||||
complex*16, intent(out) :: p, k, xp
|
||||
|
||||
double precision :: tmp_mod
|
||||
complex*16 :: p_inv
|
||||
complex*16 :: xab, ab
|
||||
|
||||
ASSERT (REAL(a) > 0.)
|
||||
ASSERT (REAL(b) > 0.)
|
||||
|
||||
! new center
|
||||
p = a + b
|
||||
|
||||
xab = xa - xb
|
||||
|
||||
p_inv = (1.d0, 0.d0) / p
|
||||
ab = a * b * p_inv
|
||||
|
||||
k = ab * xab*xab
|
||||
tmp_mod = dsqrt(REAL(k)*REAL(k) + AIMAG(k)*AIMAG(k))
|
||||
if(tmp_mod > 40.d0) then
|
||||
k = (0.d0, 0.d0)
|
||||
xp = (0.d0, 0.d0)
|
||||
return
|
||||
endif
|
||||
|
||||
k = zexp(-k)
|
||||
xp = (a*xa + b*xb) * p_inv
|
||||
|
||||
end subroutine cgaussian_product_x
|
||||
|
||||
! ---
|
||||
|
||||
subroutine multiply_cpoly(b, nb, c, nc, d, nd)
|
||||
|
||||
BEGIN_DOC
|
||||
! Multiply two complex polynomials
|
||||
! D(t) += B(t) * C(t)
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
|
||||
integer, intent(in) :: nb, nc
|
||||
complex*16, intent(in) :: b(0:nb), c(0:nc)
|
||||
complex*16, intent(inout) :: d(0:nb+nc)
|
||||
integer, intent(out) :: nd
|
||||
|
||||
integer :: ndtmp, ib, ic
|
||||
double precision :: tmp_mod
|
||||
complex*16 :: tmp
|
||||
|
||||
if(ior(nc, nb) >= 0) then ! True if nc>=0 and nb>=0
|
||||
continue
|
||||
else
|
||||
return
|
||||
endif
|
||||
|
||||
ndtmp = nb + nc
|
||||
|
||||
do ic = 0, nc
|
||||
d(ic) = d(ic) + c(ic) * b(0)
|
||||
enddo
|
||||
|
||||
do ib = 1, nb
|
||||
d(ib) = d(ib) + c(0) * b(ib)
|
||||
do ic = 1, nc
|
||||
d(ib+ic) = d(ib+ic) + c(ic) * b(ib)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
do nd = ndtmp, 0, -1
|
||||
tmp = d(nd)
|
||||
tmp_mod = dsqrt(REAL(tmp)*REAL(tmp) + AIMAG(tmp)*AIMAG(tmp))
|
||||
if(tmp_mod .lt. 1.d-15) cycle
|
||||
exit
|
||||
enddo
|
||||
|
||||
end subroutine multiply_cpoly
|
||||
|
||||
! ---
|
||||
|
||||
subroutine add_cpoly(b, nb, c, nc, d, nd)
|
||||
|
||||
BEGIN_DOC
|
||||
! Add two complex polynomials
|
||||
! D(t) += B(t) + C(t)
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
complex*16, intent(in) :: b(0:nb), c(0:nc)
|
||||
integer, intent(inout) :: nb, nc
|
||||
integer, intent(out) :: nd
|
||||
complex*16, intent(out) :: d(0:nb+nc)
|
||||
|
||||
integer :: ib
|
||||
double precision :: tmp_mod
|
||||
complex*16 :: tmp
|
||||
|
||||
nd = nb + nc
|
||||
do ib = 0, max(nb, nc)
|
||||
d(ib) = d(ib) + c(ib) + b(ib)
|
||||
enddo
|
||||
|
||||
tmp = d(nd)
|
||||
tmp_mod = dsqrt(REAL(tmp)*REAL(tmp) + AIMAG(tmp)*AIMAG(tmp))
|
||||
do while( (tmp_mod .lt. 1.d-15) .and. (nd >= 0) )
|
||||
nd -= 1
|
||||
tmp = d(nd)
|
||||
tmp_mod = dsqrt(REAL(tmp)*REAL(tmp) + AIMAG(tmp)*AIMAG(tmp))
|
||||
if(nd < 0) exit
|
||||
enddo
|
||||
|
||||
end subroutine add_cpoly
|
||||
|
||||
! ---
|
||||
|
||||
subroutine add_cpoly_multiply(b, nb, cst, d, nd)
|
||||
|
||||
BEGIN_DOC
|
||||
! Add a complex polynomial multiplied by a complex constant
|
||||
! D(t) += cst * B(t)
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
|
||||
integer, intent(in) :: nb
|
||||
complex*16, intent(in) :: b(0:nb), cst
|
||||
integer, intent(inout) :: nd
|
||||
complex*16, intent(inout) :: d(0:max(nb, nd))
|
||||
|
||||
integer :: ib
|
||||
double precision :: tmp_mod
|
||||
complex*16 :: tmp
|
||||
|
||||
nd = max(nd, nb)
|
||||
if(nd /= -1) then
|
||||
|
||||
do ib = 0, nb
|
||||
d(ib) = d(ib) + cst * b(ib)
|
||||
enddo
|
||||
|
||||
tmp = d(nd)
|
||||
tmp_mod = dsqrt(REAL(tmp)*REAL(tmp) + AIMAG(tmp)*AIMAG(tmp))
|
||||
do while(tmp_mod .lt. 1.d-15)
|
||||
nd -= 1
|
||||
if(nd < 0) exit
|
||||
tmp = d(nd)
|
||||
tmp_mod = dsqrt(REAL(tmp)*REAL(tmp) + AIMAG(tmp)*AIMAG(tmp))
|
||||
enddo
|
||||
|
||||
endif
|
||||
|
||||
end subroutine add_cpoly_multiply
|
||||
|
||||
! ---
|
||||
|
||||
subroutine recentered_cpoly2(P_A, x_A, x_P, a, P_B, x_B, x_Q, b)
|
||||
|
||||
BEGIN_DOC
|
||||
!
|
||||
! write two complex polynomials (x-x_A)^a (x-x_B)^b
|
||||
! as P_A(x-x_P) and P_B(x-x_Q)
|
||||
!
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
|
||||
integer, intent(in) :: a, b
|
||||
complex*16, intent(in) :: x_A, x_P, x_B, x_Q
|
||||
complex*16, intent(out) :: P_A(0:a), P_B(0:b)
|
||||
|
||||
integer :: i, minab, maxab
|
||||
complex*16 :: pows_a(-2:a+b+4), pows_b(-2:a+b+4)
|
||||
|
||||
double precision :: binom_func
|
||||
|
||||
if((a<0) .or. (b<0)) return
|
||||
|
||||
maxab = max(a, b)
|
||||
minab = max(min(a, b), 0)
|
||||
|
||||
pows_a(0) = (1.d0, 0.d0)
|
||||
pows_a(1) = x_P - x_A
|
||||
|
||||
pows_b(0) = (1.d0, 0.d0)
|
||||
pows_b(1) = x_Q - x_B
|
||||
|
||||
do i = 2, maxab
|
||||
pows_a(i) = pows_a(i-1) * pows_a(1)
|
||||
pows_b(i) = pows_b(i-1) * pows_b(1)
|
||||
enddo
|
||||
|
||||
P_A(0) = pows_a(a)
|
||||
P_B(0) = pows_b(b)
|
||||
|
||||
do i = 1, min(minab, 20)
|
||||
P_A(i) = binom_transp(a-i,a) * pows_a(a-i)
|
||||
P_B(i) = binom_transp(b-i,b) * pows_b(b-i)
|
||||
enddo
|
||||
|
||||
do i = minab+1, min(a, 20)
|
||||
P_A(i) = binom_transp(a-i,a) * pows_a(a-i)
|
||||
enddo
|
||||
do i = minab+1, min(b, 20)
|
||||
P_B(i) = binom_transp(b-i,b) * pows_b(b-i)
|
||||
enddo
|
||||
|
||||
do i = 101, a
|
||||
P_A(i) = binom_func(a,a-i) * pows_a(a-i)
|
||||
enddo
|
||||
do i = 101, b
|
||||
P_B(i) = binom_func(b,b-i) * pows_b(b-i)
|
||||
enddo
|
||||
|
||||
end subroutine recentered_cpoly2
|
||||
|
||||
! ---
|
||||
|
||||
complex*16 function Fc_integral(n, inv_sq_p)
|
||||
|
||||
BEGIN_DOC
|
||||
! function that calculates the following integral
|
||||
! \int_{\-infty}^{+\infty} x^n \exp(-p x^2) dx
|
||||
! for complex valued p
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
include 'constants.include.F'
|
||||
|
||||
integer, intent(in) :: n
|
||||
complex*16, intent(in) :: inv_sq_p
|
||||
|
||||
! (n)!
|
||||
double precision :: fact
|
||||
|
||||
if(n < 0) then
|
||||
Fc_integral = (0.d0, 0.d0)
|
||||
return
|
||||
endif
|
||||
|
||||
! odd n
|
||||
if(iand(n, 1) .ne. 0) then
|
||||
Fc_integral = (0.d0, 0.d0)
|
||||
return
|
||||
endif
|
||||
|
||||
if(n == 0) then
|
||||
Fc_integral = sqpi * inv_sq_p
|
||||
return
|
||||
endif
|
||||
|
||||
Fc_integral = sqpi * 0.5d0**n * inv_sq_p**dble(n+1) * fact(n) / fact(shiftr(n, 1))
|
||||
|
||||
end function Fc_integral
|
||||
|
||||
! ---
|
||||
|
||||
complex*16 function crint(n, rho)
|
||||
|
||||
implicit none
|
||||
include 'constants.include.F'
|
||||
|
||||
integer, intent(in) :: n
|
||||
complex*16, intent(in) :: rho
|
||||
|
||||
integer :: i, mmax
|
||||
double precision :: rho_mod, rho_re, rho_im
|
||||
double precision :: sq_rho_re, sq_rho_im
|
||||
double precision :: n_tmp
|
||||
complex*16 :: sq_rho, rho_inv, rho_exp
|
||||
|
||||
complex*16 :: crint_smallz, cpx_erf
|
||||
|
||||
rho_re = REAL (rho)
|
||||
rho_im = AIMAG(rho)
|
||||
rho_mod = dsqrt(rho_re*rho_re + rho_im*rho_im)
|
||||
|
||||
if(rho_mod < 10.d0) then
|
||||
! small z
|
||||
|
||||
if(rho_mod .lt. 1.d-10) then
|
||||
crint = 1.d0 / dble(n + n + 1)
|
||||
else
|
||||
crint = crint_smallz(n, rho)
|
||||
endif
|
||||
|
||||
else
|
||||
! large z
|
||||
|
||||
if(rho_mod .gt. 40.d0) then
|
||||
|
||||
n_tmp = dble(n) + 0.5d0
|
||||
crint = 0.5d0 * gamma(n_tmp) / (rho**n_tmp)
|
||||
|
||||
else
|
||||
|
||||
! get \sqrt(rho)
|
||||
sq_rho_re = sq_op5 * dsqrt(rho_re + rho_mod)
|
||||
sq_rho_im = 0.5d0 * rho_im / sq_rho_re
|
||||
sq_rho = sq_rho_re + (0.d0, 1.d0) * sq_rho_im
|
||||
|
||||
rho_exp = 0.5d0 * zexp(-rho)
|
||||
rho_inv = (1.d0, 0.d0) / rho
|
||||
|
||||
crint = 0.5d0 * sqpi * cpx_erf(sq_rho_re, sq_rho_im) / sq_rho
|
||||
mmax = n
|
||||
if(mmax .gt. 0) then
|
||||
do i = 0, mmax-1
|
||||
crint = ((dble(i) + 0.5d0) * crint - rho_exp) * rho_inv
|
||||
enddo
|
||||
endif
|
||||
|
||||
! ***
|
||||
|
||||
endif
|
||||
|
||||
endif
|
||||
|
||||
! print *, n, real(rho), real(crint)
|
||||
|
||||
end function crint
|
||||
|
||||
! ---
|
||||
|
||||
complex*16 function crint_sum(n_pt_out, rho, d1)
|
||||
|
||||
implicit none
|
||||
include 'constants.include.F'
|
||||
|
||||
integer, intent(in) :: n_pt_out
|
||||
complex*16, intent(in) :: rho, d1(0:n_pt_out)
|
||||
|
||||
integer :: n, i, mmax
|
||||
double precision :: rho_mod, rho_re, rho_im
|
||||
double precision :: sq_rho_re, sq_rho_im
|
||||
complex*16 :: sq_rho, F0
|
||||
complex*16 :: rho_tmp, rho_inv, rho_exp
|
||||
complex*16, allocatable :: Fm(:)
|
||||
|
||||
complex*16 :: crint_smallz, cpx_erf
|
||||
|
||||
rho_re = REAL (rho)
|
||||
rho_im = AIMAG(rho)
|
||||
rho_mod = dsqrt(rho_re*rho_re + rho_im*rho_im)
|
||||
|
||||
if(rho_mod < 10.d0) then
|
||||
! small z
|
||||
|
||||
if(rho_mod .lt. 1.d-10) then
|
||||
|
||||
! print *, ' 111'
|
||||
! print *, ' rho = ', rho
|
||||
|
||||
crint_sum = d1(0)
|
||||
! print *, 0, 1
|
||||
|
||||
do i = 2, n_pt_out, 2
|
||||
|
||||
n = shiftr(i, 1)
|
||||
crint_sum = crint_sum + d1(i) / dble(n+n+1)
|
||||
|
||||
! print *, n, 1.d0 / dble(n+n+1)
|
||||
enddo
|
||||
|
||||
! ***
|
||||
|
||||
else
|
||||
|
||||
! print *, ' 222'
|
||||
! print *, ' rho = ', real(rho)
|
||||
! if(abs(aimag(rho)) .gt. 1d-15) then
|
||||
! print *, ' complex rho', rho
|
||||
! stop
|
||||
! endif
|
||||
|
||||
crint_sum = d1(0) * crint_smallz(0, rho)
|
||||
|
||||
! print *, 0, real(d1(0)), real(crint_smallz(0, rho))
|
||||
! if(abs(aimag(d1(0))) .gt. 1d-15) then
|
||||
! print *, ' complex d1(0)', d1(0)
|
||||
! stop
|
||||
! endif
|
||||
|
||||
do i = 2, n_pt_out, 2
|
||||
n = shiftr(i, 1)
|
||||
crint_sum = crint_sum + d1(i) * crint_smallz(n, rho)
|
||||
|
||||
! print *, n, real(d1(i)), real(crint_smallz(n, rho))
|
||||
! if(abs(aimag(d1(i))) .gt. 1d-15) then
|
||||
! print *, ' complex d1(i)', i, d1(i)
|
||||
! stop
|
||||
! endif
|
||||
|
||||
enddo
|
||||
|
||||
! print *, 'sum = ', real(crint_sum)
|
||||
! if(abs(aimag(crint_sum)) .gt. 1d-15) then
|
||||
! print *, ' complex crint_sum', crint_sum
|
||||
! stop
|
||||
! endif
|
||||
|
||||
! ***
|
||||
|
||||
endif
|
||||
|
||||
else
|
||||
! large z
|
||||
|
||||
if(rho_mod .gt. 40.d0) then
|
||||
|
||||
! print *, ' 333'
|
||||
! print *, ' rho = ', rho
|
||||
|
||||
rho_inv = (1.d0, 0.d0) / rho
|
||||
rho_tmp = 0.5d0 * sqpi * zsqrt(rho_inv)
|
||||
crint_sum = rho_tmp * d1(0)
|
||||
! print *, 0, rho_tmp
|
||||
|
||||
do i = 2, n_pt_out, 2
|
||||
n = shiftr(i, 1)
|
||||
rho_tmp = rho_tmp * (dble(n) + 0.5d0) * rho_inv
|
||||
crint_sum = crint_sum + rho_tmp * d1(i)
|
||||
! print *, n, rho_tmp
|
||||
enddo
|
||||
|
||||
! ***
|
||||
|
||||
else
|
||||
|
||||
! print *, ' 444'
|
||||
! print *, ' rho = ', rho
|
||||
|
||||
! get \sqrt(rho)
|
||||
sq_rho_re = sq_op5 * dsqrt(rho_re + rho_mod)
|
||||
sq_rho_im = 0.5d0 * rho_im / sq_rho_re
|
||||
sq_rho = sq_rho_re + (0.d0, 1.d0) * sq_rho_im
|
||||
!sq_rho = zsqrt(rho)
|
||||
|
||||
|
||||
F0 = 0.5d0 * sqpi * cpx_erf(sq_rho_re, sq_rho_im) / sq_rho
|
||||
crint_sum = F0 * d1(0)
|
||||
! print *, 0, F0
|
||||
|
||||
rho_exp = 0.5d0 * zexp(-rho)
|
||||
rho_inv = (1.d0, 0.d0) / rho
|
||||
|
||||
mmax = shiftr(n_pt_out, 1)
|
||||
if(mmax .gt. 0) then
|
||||
|
||||
allocate( Fm(mmax) )
|
||||
Fm(1:mmax) = (0.d0, 0.d0)
|
||||
|
||||
do n = 0, mmax-1
|
||||
F0 = ((dble(n) + 0.5d0) * F0 - rho_exp) * rho_inv
|
||||
Fm(n+1) = F0
|
||||
! print *, n, F0
|
||||
enddo
|
||||
|
||||
do i = 2, n_pt_out, 2
|
||||
n = shiftr(i, 1)
|
||||
crint_sum = crint_sum + Fm(n) * d1(i)
|
||||
enddo
|
||||
deallocate(Fm)
|
||||
endif
|
||||
|
||||
! ***
|
||||
|
||||
endif
|
||||
|
||||
endif
|
||||
|
||||
end function crint_sum
|
||||
|
||||
! ---
|
||||
|
||||
complex*16 function crint_smallz(n, rho)
|
||||
|
||||
BEGIN_DOC
|
||||
! Standard version of rint
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
integer, intent(in) :: n
|
||||
complex*16, intent(in) :: rho
|
||||
|
||||
integer, parameter :: kmax = 40
|
||||
double precision, parameter :: eps = 1.d-13
|
||||
|
||||
integer :: k
|
||||
double precision :: delta_mod
|
||||
complex*16 :: rho_k, ct, delta_k
|
||||
|
||||
ct = 0.5d0 * zexp(-rho) * gamma(dble(n) + 0.5d0)
|
||||
rho_k = (1.d0, 0.d0)
|
||||
crint_smallz = ct * rho_k / gamma(dble(n) + 1.5d0)
|
||||
|
||||
do k = 1, kmax
|
||||
|
||||
rho_k = rho_k * rho
|
||||
delta_k = ct * rho_k / gamma(dble(n+k) + 1.5d0)
|
||||
crint_smallz = crint_smallz + delta_k
|
||||
|
||||
delta_mod = dsqrt(REAL(delta_k)*REAL(delta_k) + AIMAG(delta_k)*AIMAG(delta_k))
|
||||
if(delta_mod .lt. eps) return
|
||||
enddo
|
||||
|
||||
if(delta_mod > eps) then
|
||||
write(*,*) ' pb in crint_smallz !'
|
||||
write(*,*) ' n, rho = ', n, rho
|
||||
write(*,*) ' delta_mod = ', delta_mod
|
||||
stop 1
|
||||
endif
|
||||
|
||||
end function crint_smallz
|
||||
|
||||
! ---
|
||||
|
204
src/utils/cpx_erf.irp.f
Normal file
204
src/utils/cpx_erf.irp.f
Normal file
@ -0,0 +1,204 @@
|
||||
|
||||
! ---
|
||||
|
||||
complex*16 function cpx_erf(x, y)
|
||||
|
||||
BEGIN_DOC
|
||||
!
|
||||
! compute erf(z) for z = x + i y
|
||||
!
|
||||
! REF: Abramowitz and Stegun
|
||||
!
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
|
||||
double precision, intent(in) :: x, y
|
||||
|
||||
double precision :: yabs
|
||||
complex*16 :: erf_tmp1, erf_tmp2, erf_tmp3, erf_tot
|
||||
|
||||
double precision :: erf_F
|
||||
complex*16 :: erf_E, erf_G, erf_H
|
||||
|
||||
yabs = dabs(y)
|
||||
|
||||
if(yabs .lt. 1.d-15) then
|
||||
|
||||
cpx_erf = (1.d0, 0.d0) * derf(x)
|
||||
return
|
||||
|
||||
else
|
||||
|
||||
erf_tmp1 = (1.d0, 0.d0) * derf(x)
|
||||
erf_tmp2 = erf_E(x, yabs) + erf_F(x, yabs)
|
||||
erf_tmp3 = zexp(-(0.d0, 2.d0) * x * yabs) * ( erf_G(x, yabs) + erf_H(x, yabs) )
|
||||
erf_tot = erf_tmp1 + erf_tmp2 - erf_tmp3
|
||||
|
||||
endif
|
||||
|
||||
if(y .gt. 0.d0) then
|
||||
cpx_erf = erf_tot
|
||||
else
|
||||
cpx_erf = CONJG(erf_tot)
|
||||
endif
|
||||
|
||||
end function cpx_erf
|
||||
|
||||
! ---
|
||||
|
||||
complex*16 function erf_E(x, yabs)
|
||||
|
||||
implicit none
|
||||
include 'constants.include.F'
|
||||
|
||||
double precision, intent(in) :: x, yabs
|
||||
|
||||
if( (dabs(x).gt.6.d0) .or. (x==0.d0) ) then
|
||||
erf_E = (0.d0, 0.d0)
|
||||
return
|
||||
endif
|
||||
|
||||
if(dabs(x) .lt. 1.d-7) then
|
||||
|
||||
erf_E = -inv_pi * (0.d0, 1.d0) * yabs
|
||||
|
||||
else
|
||||
|
||||
erf_E = 0.5d0 * inv_pi * dexp(-x*x) &
|
||||
* ((1.d0, 0.d0) - zexp(-(2.d0, 0.d0) * x * yabs)) / x
|
||||
|
||||
endif
|
||||
|
||||
end function erf_E
|
||||
|
||||
! ---
|
||||
|
||||
double precision function erf_F(x, yabs)
|
||||
|
||||
implicit none
|
||||
include 'constants.include.F'
|
||||
|
||||
double precision, intent(in) :: x, yabs
|
||||
|
||||
integer, parameter :: Nmax = 13
|
||||
|
||||
integer :: i
|
||||
double precision :: tmp1, tmp2, x2, ct
|
||||
|
||||
|
||||
if(dabs(x) .gt. 5.8d0) then
|
||||
|
||||
erf_F = 0.d0
|
||||
|
||||
else
|
||||
|
||||
x2 = x * x
|
||||
ct = x * inv_pi
|
||||
|
||||
erf_F = 0.d0
|
||||
do i = 1, Nmax
|
||||
|
||||
tmp1 = 0.25d0 * dble(i) * dble(i) + x2
|
||||
tmp2 = dexp(-tmp1) / tmp1
|
||||
erf_F = erf_F + tmp2
|
||||
|
||||
if(dabs(tmp2) .lt. 1d-15) exit
|
||||
enddo
|
||||
erf_F = ct * erf_F
|
||||
|
||||
endif
|
||||
|
||||
end function erf_F
|
||||
|
||||
! ---
|
||||
|
||||
complex*16 function erf_G(x, yabs)
|
||||
|
||||
implicit none
|
||||
include 'constants.include.F'
|
||||
|
||||
double precision, intent(in) :: x, yabs
|
||||
|
||||
integer, parameter :: Nmax = 13
|
||||
|
||||
integer :: i, tmpi, imin, imax
|
||||
double precision :: tmp0, tmp1, x2, idble
|
||||
complex*16 :: tmp2
|
||||
|
||||
if(x .eq. 0.d0) then
|
||||
erf_G = (0.d0, 0.d0)
|
||||
return
|
||||
endif
|
||||
|
||||
tmpi = int(2.d0 * yabs)
|
||||
imin = max(1, tmpi-Nmax)
|
||||
imax = tmpi + Nmax
|
||||
|
||||
x2 = x * x
|
||||
|
||||
erf_G = 0.d0
|
||||
do i = imin, imax
|
||||
|
||||
idble = dble(i)
|
||||
tmp0 = 0.5d0 * idble
|
||||
tmp1 = tmp0 * tmp0 + x2
|
||||
tmp2 = dexp( idble * yabs - tmp1 - dlog(tmp1) - dlog_2pi) * (x - (0.d0, 1.d0)*tmp0)
|
||||
|
||||
erf_G = erf_G + tmp2
|
||||
|
||||
enddo
|
||||
|
||||
end function erf_G
|
||||
|
||||
! ---
|
||||
|
||||
complex*16 function erf_H(x, yabs)
|
||||
|
||||
implicit none
|
||||
include 'constants.include.F'
|
||||
|
||||
double precision, intent(in) :: x, yabs
|
||||
|
||||
integer, parameter :: Nmax = 13
|
||||
|
||||
integer :: i
|
||||
double precision :: tmp0, tmp1, tmp_mod, x2, ct, idble
|
||||
complex*16 :: tmp2
|
||||
|
||||
if(x .eq. 0.d0) then
|
||||
erf_H = (0.d0, 0.d0)
|
||||
return
|
||||
endif
|
||||
|
||||
|
||||
if( (dabs(x) .lt. 10d0) .and. (yabs .lt. 6.1d0) ) then
|
||||
|
||||
x2 = x * x
|
||||
ct = 0.5d0 * inv_pi
|
||||
|
||||
erf_H = 0.d0
|
||||
do i = 1, Nmax
|
||||
|
||||
idble = dble(i)
|
||||
tmp0 = 0.5d0 * idble
|
||||
tmp1 = tmp0 * tmp0 + x2
|
||||
tmp2 = dexp(-tmp1-idble*yabs) * (x + (0.d0, 1.d0)*tmp0) / tmp1
|
||||
erf_H = erf_H + tmp2
|
||||
|
||||
tmp_mod = dsqrt(REAL(tmp2)*REAL(tmp2) + AIMAG(tmp2)*AIMAG(tmp2))
|
||||
if(tmp_mod .lt. 1d-15) exit
|
||||
enddo
|
||||
erf_H = ct * erf_H
|
||||
|
||||
else
|
||||
|
||||
erf_H = (0.d0, 0.d0)
|
||||
|
||||
endif
|
||||
|
||||
end function erf_H
|
||||
|
||||
! ---
|
||||
|
||||
|
@ -133,7 +133,7 @@ subroutine give_explicit_poly_and_gaussian_v(P_new, ldp, P_center, p, fact_k, io
|
||||
|
||||
BEGIN_DOC
|
||||
! Transforms the product of
|
||||
! (x-x_A)^a(1) (x-x_B)^b(1) (x-x_A)^a(2) (y-y_B)^b(2) (z-z_A)^a(3) (z-z_B)^b(3) exp(-(r-A)^2 alpha) exp(-(r-B)^2 beta)
|
||||
! (x-x_A)^a(1) (x-x_B)^b(1) (y-y_A)^a(2) (y-y_B)^b(2) (z-z_A)^a(3) (z-z_B)^b(3) exp(-(r-A)^2 alpha) exp(-(r-B)^2 beta)
|
||||
! into
|
||||
! fact_k * [ sum (l_x = 0,i_order(1)) P_new(l_x,1) * (x-P_center(1))^l_x ] exp (- p (x-P_center(1))^2 )
|
||||
! * [ sum (l_y = 0,i_order(2)) P_new(l_y,2) * (y-P_center(2))^l_y ] exp (- p (y-P_center(2))^2 )
|
||||
@ -427,6 +427,46 @@ subroutine gaussian_product_x(a,xa,b,xb,k,p,xp)
|
||||
end subroutine
|
||||
|
||||
|
||||
!-
|
||||
|
||||
subroutine gaussian_product_x_v(a,xa,b,xb,k,p,xp,n_points)
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Gaussian product in 1D with multiple xa
|
||||
! e^{-a (x-x_A)^2} e^{-b (x-x_B)^2} = K_{ab}^x e^{-p (x-x_P)^2}
|
||||
END_DOC
|
||||
|
||||
integer, intent(in) :: n_points
|
||||
double precision , intent(in) :: a,b ! Exponents
|
||||
double precision , intent(in) :: xa(n_points),xb ! Centers
|
||||
double precision , intent(out) :: p(n_points) ! New exponent
|
||||
double precision , intent(out) :: xp(n_points) ! New center
|
||||
double precision , intent(out) :: k(n_points) ! Constant
|
||||
|
||||
double precision :: p_inv
|
||||
integer :: ipoint
|
||||
|
||||
ASSERT (a>0.)
|
||||
ASSERT (b>0.)
|
||||
|
||||
double precision :: xab, ab
|
||||
|
||||
p = a+b
|
||||
p_inv = 1.d0/(a+b)
|
||||
ab = a*b*p_inv
|
||||
do ipoint = 1, n_points
|
||||
xab = xa(ipoint)-xb
|
||||
k(ipoint) = ab*xab*xab
|
||||
if (k(ipoint) > 40.d0) then
|
||||
k(ipoint)=0.d0
|
||||
cycle
|
||||
endif
|
||||
k(ipoint) = exp(-k(ipoint))
|
||||
xp(ipoint) = (a*xa(ipoint)+b*xb)*p_inv
|
||||
enddo
|
||||
end subroutine
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
@ -506,8 +546,10 @@ subroutine multiply_poly_v(b,nb,c,nc,d,nd,n_points)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
end
|
||||
|
||||
|
||||
subroutine add_poly(b,nb,c,nc,d,nd)
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
@ -1041,3 +1083,94 @@ double precision function rint1(n,rho)
|
||||
write(*,*)'pb in rint1 k too large!'
|
||||
stop 1
|
||||
end
|
||||
|
||||
! ---
|
||||
|
||||
double precision function V_phi(n, m)
|
||||
|
||||
BEGIN_DOC
|
||||
! Computes the angular $\phi$ part of the nuclear attraction integral:
|
||||
!
|
||||
! $\int_{0}^{2 \pi} \cos(\phi)^n \sin(\phi)^m d\phi$.
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
integer, intent(in) :: n, m
|
||||
|
||||
integer :: i
|
||||
double precision :: prod
|
||||
|
||||
double precision :: Wallis
|
||||
|
||||
prod = 1.d0
|
||||
do i = 0, shiftr(n, 1)-1
|
||||
prod = prod/ (1.d0 + dfloat(m+1)/dfloat(n-i-i-1))
|
||||
enddo
|
||||
V_phi = 4.d0 * prod * Wallis(m)
|
||||
|
||||
end function V_phi
|
||||
|
||||
! ---
|
||||
|
||||
double precision function V_theta(n, m)
|
||||
|
||||
BEGIN_DOC
|
||||
! Computes the angular $\theta$ part of the nuclear attraction integral:
|
||||
!
|
||||
! $\int_{0}^{\pi} \cos(\theta)^n \sin(\theta)^m d\theta$
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
include 'utils/constants.include.F'
|
||||
integer, intent(in) :: n, m
|
||||
|
||||
integer :: i
|
||||
double precision :: prod
|
||||
|
||||
double precision :: Wallis
|
||||
|
||||
V_theta = 0.d0
|
||||
prod = 1.d0
|
||||
do i = 0, shiftr(n, 1)-1
|
||||
prod = prod / (1.d0 + dfloat(m+1)/dfloat(n-i-i-1))
|
||||
enddo
|
||||
V_theta = (prod + prod) * Wallis(m)
|
||||
|
||||
end function V_theta
|
||||
|
||||
! ---
|
||||
|
||||
double precision function Wallis(n)
|
||||
|
||||
BEGIN_DOC
|
||||
! Wallis integral:
|
||||
!
|
||||
! $\int_{0}^{\pi} \cos(\theta)^n d\theta$.
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
include 'utils/constants.include.F'
|
||||
|
||||
integer, intent(in) :: n
|
||||
|
||||
integer :: p
|
||||
|
||||
double precision :: fact
|
||||
|
||||
if(iand(n, 1) .eq. 0) then
|
||||
|
||||
Wallis = fact(shiftr(n, 1))
|
||||
Wallis = pi * fact(n) / (dble(ibset(0_8, n)) * (Wallis + Wallis) * Wallis)
|
||||
|
||||
else
|
||||
|
||||
p = shiftr(n, 1)
|
||||
Wallis = fact(p)
|
||||
Wallis = dble(ibset(0_8, p+p)) * Wallis * Wallis / fact(p+p+1)
|
||||
|
||||
endif
|
||||
|
||||
end function Wallis
|
||||
|
||||
! ---
|
||||
|
||||
|
@ -32,9 +32,8 @@ double precision function overlap_gaussian_x(A_center,B_center,alpha,beta,power_
|
||||
end
|
||||
|
||||
|
||||
subroutine overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,&
|
||||
power_B,overlap_x,overlap_y,overlap_z,overlap,dim)
|
||||
implicit none
|
||||
subroutine overlap_gaussian_xyz(A_center, B_center, alpha, beta, power_A, power_B, overlap_x, overlap_y, overlap_z, overlap, dim)
|
||||
|
||||
BEGIN_DOC
|
||||
!.. math::
|
||||
!
|
||||
@ -42,7 +41,10 @@ subroutine overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,&
|
||||
! S = S_x S_y S_z
|
||||
!
|
||||
END_DOC
|
||||
|
||||
include 'constants.include.F'
|
||||
|
||||
implicit none
|
||||
integer,intent(in) :: dim ! dimension maximum for the arrays representing the polynomials
|
||||
double precision,intent(in) :: A_center(3),B_center(3) ! center of the x1 functions
|
||||
double precision, intent(in) :: alpha,beta
|
||||
@ -51,17 +53,18 @@ subroutine overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,&
|
||||
double precision :: P_new(0:max_dim,3),P_center(3),fact_p,p
|
||||
double precision :: F_integral_tab(0:max_dim)
|
||||
integer :: iorder_p(3)
|
||||
|
||||
call give_explicit_poly_and_gaussian(P_new,P_center,p,fact_p,iorder_p,alpha,beta,power_A,power_B,A_center,B_center,dim)
|
||||
if(fact_p.lt.1d-20)then
|
||||
overlap_x = 1.d-10
|
||||
overlap_y = 1.d-10
|
||||
overlap_z = 1.d-10
|
||||
overlap = 1.d-10
|
||||
return
|
||||
endif
|
||||
integer :: nmax
|
||||
double precision :: F_integral
|
||||
|
||||
call give_explicit_poly_and_gaussian(P_new, P_center, p, fact_p, iorder_p, alpha, beta, power_A, power_B, A_center, B_center, dim)
|
||||
if(fact_p.lt.1d-20)then
|
||||
overlap_x = 1.d-10
|
||||
overlap_y = 1.d-10
|
||||
overlap_z = 1.d-10
|
||||
overlap = 1.d-10
|
||||
return
|
||||
endif
|
||||
|
||||
nmax = maxval(iorder_p)
|
||||
do i = 0,nmax
|
||||
F_integral_tab(i) = F_integral(i,p)
|
||||
@ -93,40 +96,47 @@ subroutine overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,&
|
||||
|
||||
end
|
||||
|
||||
! ---
|
||||
|
||||
subroutine overlap_x_abs(A_center, B_center, alpha, beta, power_A, power_B, overlap_x, lower_exp_val, dx, nx)
|
||||
|
||||
subroutine overlap_x_abs(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,lower_exp_val,dx,nx)
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! .. math ::
|
||||
!
|
||||
! \int_{-infty}^{+infty} (x-A_center)^(power_A) * (x-B_center)^power_B * exp(-alpha(x-A_center)^2) * exp(-beta(x-B_center)^2) dx
|
||||
!
|
||||
END_DOC
|
||||
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
|
||||
|
||||
implicit none
|
||||
|
||||
integer, intent(in) :: power_A, power_B, nx
|
||||
double precision, intent(in) :: lower_exp_val, A_center, B_center, alpha, beta
|
||||
double precision, intent(out) :: overlap_x, dx
|
||||
|
||||
integer :: i, j, k, l
|
||||
double precision :: x_min, x_max, domain, x, factor, dist, p, p_inv, rho
|
||||
double precision :: P_center
|
||||
double precision :: tmp
|
||||
|
||||
if(power_A.lt.0 .or. power_B.lt.0) then
|
||||
overlap_x = 0.d0
|
||||
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 = 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
|
||||
if(rho*dist.gt.80.d0)then
|
||||
|
||||
if(rho*dist.gt.80.d0) then
|
||||
overlap_x= 0.d0
|
||||
return
|
||||
endif
|
||||
|
||||
factor = dexp(-rho * dist)
|
||||
|
||||
double precision :: tmp
|
||||
|
||||
tmp = dsqrt(lower_exp_val/p)
|
||||
x_min = P_center - tmp
|
||||
@ -143,7 +153,7 @@ subroutine overlap_x_abs(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,
|
||||
overlap_x = factor * dx * overlap_x
|
||||
end
|
||||
|
||||
|
||||
! ---
|
||||
|
||||
subroutine overlap_gaussian_xyz_v(A_center, B_center, alpha, beta, power_A, power_B, overlap, n_points)
|
||||
|
||||
@ -173,7 +183,7 @@ subroutine overlap_gaussian_xyz_v(A_center, B_center, alpha, beta, power_A, powe
|
||||
double precision :: F_integral
|
||||
double precision, allocatable :: P_new(:,:,:), P_center(:,:), fact_p(:)
|
||||
|
||||
ldp = maxval(power_A(1:3) + power_B(1:3))
|
||||
ldp = maxval( power_A(1:3) + power_B(1:3) )
|
||||
|
||||
allocate(P_new(n_points,0:ldp,3), P_center(n_points,3), fact_p(n_points))
|
||||
|
||||
|
@ -460,6 +460,33 @@ subroutine v2_over_x(v,x,res)
|
||||
|
||||
end
|
||||
|
||||
! ---
|
||||
|
||||
subroutine check_sym(A, n)
|
||||
|
||||
implicit none
|
||||
integer, intent(in) :: n
|
||||
double precision, intent(in) :: A(n,n)
|
||||
integer :: i, j
|
||||
double precision :: dev_sym, norm, tmp
|
||||
|
||||
dev_sym = 0.d0
|
||||
norm = 0.d0
|
||||
do i = 1, n
|
||||
do j = i+1, n
|
||||
tmp = A(j,i) - A(i,j)
|
||||
dev_sym += tmp * tmp
|
||||
norm += A(j,i) * A(j,i)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
print*, ' deviation from sym = ', dev_sym
|
||||
print*, ' norm = ', norm
|
||||
|
||||
end subroutine check_sym
|
||||
|
||||
! ---
|
||||
|
||||
subroutine sum_A_At(A, N)
|
||||
|
||||
!BEGIN_DOC
|
||||
|
Loading…
Reference in New Issue
Block a user