PTEROSOR/QuantumPackage
10
0
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-06-21 20:52:28 +02:00
Code Releases Activity

list2, list3, j1b, grad_j1, lapl_j1b ok

Browse Source
This commit is contained in:
AbdAmmar 2022-10-17 10:55:53 +02:00
parent 4f0a0f68fc
commit dbbae1f990
10 changed files with 1253 additions and 700 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

118
src/ao_many_one_e_ints/grad2_jmu_modif.irp.f
View File

@ -1,124 +1,6 @@
! ---
BEGIN_PROVIDER [ integer, List_all_comb_b3_size]
implicit none
List_all_comb_b3_size = 3**nucl_num
END_PROVIDER
! ---
BEGIN_PROVIDER [ integer, List_all_comb_b3, (nucl_num, List_all_comb_b3_size)]
implicit none
integer :: i, j, ii, jj
integer, allocatable :: M(:,:), p(:)
if(nucl_num .gt. 32) then
print *, ' nucl_num = ', nucl_num, '> 32'
stop
endif
List_all_comb_b3(:,:) = 0
List_all_comb_b3(:,List_all_comb_b3_size) = 2
allocate(p(nucl_num))
p = 0
do i = 2, List_all_comb_b3_size-1
do j = 1, nucl_num
ii = 0
do jj = 1, j-1, 1
ii = ii + p(jj) * 3**(jj-1)
enddo
p(j) = modulo(i-1-ii, 3**j) / 3**(j-1)
List_all_comb_b3(j,i) = p(j)
enddo
enddo
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, List_all_comb_b3_coef, ( List_all_comb_b3_size)]
&BEGIN_PROVIDER [ double precision, List_all_comb_b3_expo, ( List_all_comb_b3_size)]
&BEGIN_PROVIDER [ double precision, List_all_comb_b3_cent, (3, List_all_comb_b3_size)]
implicit none
integer :: i, j, k, phase
double precision :: tmp_alphaj, tmp_alphak, facto
provide j1b_pen
List_all_comb_b3_coef = 0.d0
List_all_comb_b3_expo = 0.d0
List_all_comb_b3_cent = 0.d0
do i = 1, List_all_comb_b3_size
do j = 1, nucl_num
tmp_alphaj = dble(List_all_comb_b3(j,i)) * j1b_pen(j)
List_all_comb_b3_expo(i) += tmp_alphaj
List_all_comb_b3_cent(1,i) += tmp_alphaj * nucl_coord(j,1)
List_all_comb_b3_cent(2,i) += tmp_alphaj * nucl_coord(j,2)
List_all_comb_b3_cent(3,i) += tmp_alphaj * nucl_coord(j,3)
enddo
ASSERT(List_all_comb_b3_expo(i) .gt. 0d0)
if(List_all_comb_b3_expo(i) .lt. 1d-10) cycle
List_all_comb_b3_cent(1,i) = List_all_comb_b3_cent(1,i) / List_all_comb_b3_expo(i)
List_all_comb_b3_cent(2,i) = List_all_comb_b3_cent(2,i) / List_all_comb_b3_expo(i)
List_all_comb_b3_cent(3,i) = List_all_comb_b3_cent(3,i) / List_all_comb_b3_expo(i)
enddo
! ---
do i = 1, List_all_comb_b3_size
do j = 2, nucl_num, 1
tmp_alphaj = dble(List_all_comb_b3(j,i)) * j1b_pen(j)
do k = 1, j-1, 1
tmp_alphak = dble(List_all_comb_b3(k,i)) * j1b_pen(k)
List_all_comb_b3_coef(i) += tmp_alphaj * tmp_alphak * ( (nucl_coord(j,1) - nucl_coord(k,1)) * (nucl_coord(j,1) - nucl_coord(k,1)) &
+ (nucl_coord(j,2) - nucl_coord(k,2)) * (nucl_coord(j,2) - nucl_coord(k,2)) &
+ (nucl_coord(j,3) - nucl_coord(k,3)) * (nucl_coord(j,3) - nucl_coord(k,3)) )
enddo
enddo
if(List_all_comb_b3_expo(i) .lt. 1d-10) cycle
List_all_comb_b3_coef(i) = List_all_comb_b3_coef(i) / List_all_comb_b3_expo(i)
enddo
! ---
do i = 1, List_all_comb_b3_size
facto = 1.d0
phase = 0
do j = 1, nucl_num
tmp_alphaj = dble(List_all_comb_b3(j,i))
facto *= 2.d0 / (gamma(tmp_alphaj+1.d0) * gamma(3.d0-tmp_alphaj))
phase += List_all_comb_b3(j,i)
enddo
List_all_comb_b3_coef(i) = (-1.d0)**dble(phase) * facto * dexp(-List_all_comb_b3_coef(i))
enddo
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, int2_grad1u_grad2u_j1b, (ao_num, ao_num, n_points_final_grid)]
BEGIN_DOC

105
src/ao_many_one_e_ints/grad_lapl_jmu_modif.irp.f
View File

@ -1,110 +1,6 @@
! ---
BEGIN_PROVIDER [ integer, List_all_comb_b2_size]
implicit none
List_all_comb_b2_size = 2**nucl_num
END_PROVIDER
! ---
BEGIN_PROVIDER [ integer, List_all_comb_b2, (nucl_num, List_all_comb_b2_size)]
implicit none
integer :: i, j
if(nucl_num .gt. 32) then
print *, ' nucl_num = ', nucl_num, '> 32'
stop
endif
List_all_comb_b2 = 0
do i = 0, List_all_comb_b2_size-1
do j = 0, nucl_num-1
if (btest(i,j)) then
List_all_comb_b2(j+1,i+1) = 1
endif
enddo
enddo
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, List_all_comb_b2_coef, ( List_all_comb_b2_size)]
&BEGIN_PROVIDER [ double precision, List_all_comb_b2_expo, ( List_all_comb_b2_size)]
&BEGIN_PROVIDER [ double precision, List_all_comb_b2_cent, (3, List_all_comb_b2_size)]
implicit none
integer :: i, j, k, phase
double precision :: tmp_alphaj, tmp_alphak
provide j1b_pen
List_all_comb_b2_coef = 0.d0
List_all_comb_b2_expo = 0.d0
List_all_comb_b2_cent = 0.d0
do i = 1, List_all_comb_b2_size
do j = 1, nucl_num
tmp_alphaj = dble(List_all_comb_b2(j,i)) * j1b_pen(j)
List_all_comb_b2_expo(i) += tmp_alphaj
List_all_comb_b2_cent(1,i) += tmp_alphaj * nucl_coord(j,1)
List_all_comb_b2_cent(2,i) += tmp_alphaj * nucl_coord(j,2)
List_all_comb_b2_cent(3,i) += tmp_alphaj * nucl_coord(j,3)
enddo
ASSERT(List_all_comb_b2_expo(i) .gt. 0d0)
if(List_all_comb_b2_expo(i) .lt. 1d-10) cycle
List_all_comb_b2_cent(1,i) = List_all_comb_b2_cent(1,i) / List_all_comb_b2_expo(i)
List_all_comb_b2_cent(2,i) = List_all_comb_b2_cent(2,i) / List_all_comb_b2_expo(i)
List_all_comb_b2_cent(3,i) = List_all_comb_b2_cent(3,i) / List_all_comb_b2_expo(i)
enddo
! ---
do i = 1, List_all_comb_b2_size
do j = 2, nucl_num, 1
tmp_alphaj = dble(List_all_comb_b2(j,i)) * j1b_pen(j)
do k = 1, j-1, 1
tmp_alphak = dble(List_all_comb_b2(k,i)) * j1b_pen(k)
List_all_comb_b2_coef(i) += tmp_alphaj * tmp_alphak * ( (nucl_coord(j,1) - nucl_coord(k,1)) * (nucl_coord(j,1) - nucl_coord(k,1)) &
+ (nucl_coord(j,2) - nucl_coord(k,2)) * (nucl_coord(j,2) - nucl_coord(k,2)) &
+ (nucl_coord(j,3) - nucl_coord(k,3)) * (nucl_coord(j,3) - nucl_coord(k,3)) )
enddo
enddo
if(List_all_comb_b2_expo(i) .lt. 1d-10) cycle
List_all_comb_b2_coef(i) = List_all_comb_b2_coef(i) / List_all_comb_b2_expo(i)
enddo
! ---
do i = 1, List_all_comb_b2_size
phase = 0
do j = 1, nucl_num
phase += List_all_comb_b2(j,i)
enddo
List_all_comb_b2_coef(i) = (-1.d0)**dble(phase) * dexp(-List_all_comb_b2_coef(i))
enddo
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_j1b, (ao_num, ao_num, n_points_final_grid)]
BEGIN_DOC
@ -157,7 +53,6 @@ BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_j1b, (ao_num, ao_num, n_po
tmp(j,i,ipoint) += coef * (int_mu - int_coulomb)
enddo
enddo
enddo
enddo

227
src/ao_many_one_e_ints/listj1b_product_to_sum.irp.f Normal file
View File

@ -0,0 +1,227 @@
! ---
BEGIN_PROVIDER [ integer, List_all_comb_b2_size]
implicit none
List_all_comb_b2_size = 2**nucl_num
END_PROVIDER
! ---
BEGIN_PROVIDER [ integer, List_all_comb_b2, (nucl_num, List_all_comb_b2_size)]
implicit none
integer :: i, j
if(nucl_num .gt. 32) then
print *, ' nucl_num = ', nucl_num, '> 32'
stop
endif
List_all_comb_b2 = 0
do i = 0, List_all_comb_b2_size-1
do j = 0, nucl_num-1
if (btest(i,j)) then
List_all_comb_b2(j+1,i+1) = 1
endif
enddo
enddo
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, List_all_comb_b2_coef, ( List_all_comb_b2_size)]
&BEGIN_PROVIDER [ double precision, List_all_comb_b2_expo, ( List_all_comb_b2_size)]
&BEGIN_PROVIDER [ double precision, List_all_comb_b2_cent, (3, List_all_comb_b2_size)]
implicit none
integer :: i, j, k, phase
double precision :: tmp_alphaj, tmp_alphak
double precision :: tmp_cent_x, tmp_cent_y, tmp_cent_z
provide j1b_pen
List_all_comb_b2_coef = 0.d0
List_all_comb_b2_expo = 0.d0
List_all_comb_b2_cent = 0.d0
do i = 1, List_all_comb_b2_size
tmp_cent_x = 0.d0
tmp_cent_y = 0.d0
tmp_cent_z = 0.d0
do j = 1, nucl_num
tmp_alphaj = dble(List_all_comb_b2(j,i)) * j1b_pen(j)
List_all_comb_b2_expo(i) += tmp_alphaj
tmp_cent_x += tmp_alphaj * nucl_coord(j,1)
tmp_cent_y += tmp_alphaj * nucl_coord(j,2)
tmp_cent_z += tmp_alphaj * nucl_coord(j,3)
enddo
ASSERT(List_all_comb_b2_expo(i) .gt. 0d0)
if(List_all_comb_b2_expo(i) .lt. 1d-10) cycle
List_all_comb_b2_cent(1,i) = tmp_cent_x / List_all_comb_b2_expo(i)
List_all_comb_b2_cent(2,i) = tmp_cent_y / List_all_comb_b2_expo(i)
List_all_comb_b2_cent(3,i) = tmp_cent_z / List_all_comb_b2_expo(i)
enddo
! ---
do i = 1, List_all_comb_b2_size
do j = 2, nucl_num, 1
tmp_alphaj = dble(List_all_comb_b2(j,i)) * j1b_pen(j)
do k = 1, j-1, 1
tmp_alphak = dble(List_all_comb_b2(k,i)) * j1b_pen(k)
List_all_comb_b2_coef(i) += tmp_alphaj * tmp_alphak * ( (nucl_coord(j,1) - nucl_coord(k,1)) * (nucl_coord(j,1) - nucl_coord(k,1)) &
+ (nucl_coord(j,2) - nucl_coord(k,2)) * (nucl_coord(j,2) - nucl_coord(k,2)) &
+ (nucl_coord(j,3) - nucl_coord(k,3)) * (nucl_coord(j,3) - nucl_coord(k,3)) )
enddo
enddo
if(List_all_comb_b2_expo(i) .lt. 1d-10) cycle
List_all_comb_b2_coef(i) = List_all_comb_b2_coef(i) / List_all_comb_b2_expo(i)
enddo
! ---
do i = 1, List_all_comb_b2_size
phase = 0
do j = 1, nucl_num
phase += List_all_comb_b2(j,i)
enddo
List_all_comb_b2_coef(i) = (-1.d0)**dble(phase) * dexp(-List_all_comb_b2_coef(i))
enddo
END_PROVIDER
! ---
BEGIN_PROVIDER [ integer, List_all_comb_b3_size]
implicit none
List_all_comb_b3_size = 3**nucl_num
END_PROVIDER
! ---
BEGIN_PROVIDER [ integer, List_all_comb_b3, (nucl_num, List_all_comb_b3_size)]
implicit none
integer :: i, j, ii, jj
integer, allocatable :: M(:,:), p(:)
if(nucl_num .gt. 32) then
print *, ' nucl_num = ', nucl_num, '> 32'
stop
endif
List_all_comb_b3(:,:) = 0
List_all_comb_b3(:,List_all_comb_b3_size) = 2
allocate(p(nucl_num))
p = 0
do i = 2, List_all_comb_b3_size-1
do j = 1, nucl_num
ii = 0
do jj = 1, j-1, 1
ii = ii + p(jj) * 3**(jj-1)
enddo
p(j) = modulo(i-1-ii, 3**j) / 3**(j-1)
List_all_comb_b3(j,i) = p(j)
enddo
enddo
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, List_all_comb_b3_coef, ( List_all_comb_b3_size)]
&BEGIN_PROVIDER [ double precision, List_all_comb_b3_expo, ( List_all_comb_b3_size)]
&BEGIN_PROVIDER [ double precision, List_all_comb_b3_cent, (3, List_all_comb_b3_size)]
implicit none
integer :: i, j, k, phase
double precision :: tmp_alphaj, tmp_alphak, facto
provide j1b_pen
List_all_comb_b3_coef = 0.d0
List_all_comb_b3_expo = 0.d0
List_all_comb_b3_cent = 0.d0
do i = 1, List_all_comb_b3_size
do j = 1, nucl_num
tmp_alphaj = dble(List_all_comb_b3(j,i)) * j1b_pen(j)
List_all_comb_b3_expo(i) += tmp_alphaj
List_all_comb_b3_cent(1,i) += tmp_alphaj * nucl_coord(j,1)
List_all_comb_b3_cent(2,i) += tmp_alphaj * nucl_coord(j,2)
List_all_comb_b3_cent(3,i) += tmp_alphaj * nucl_coord(j,3)
enddo
ASSERT(List_all_comb_b3_expo(i) .gt. 0d0)
if(List_all_comb_b3_expo(i) .lt. 1d-10) cycle
List_all_comb_b3_cent(1,i) = List_all_comb_b3_cent(1,i) / List_all_comb_b3_expo(i)
List_all_comb_b3_cent(2,i) = List_all_comb_b3_cent(2,i) / List_all_comb_b3_expo(i)
List_all_comb_b3_cent(3,i) = List_all_comb_b3_cent(3,i) / List_all_comb_b3_expo(i)
enddo
! ---
do i = 1, List_all_comb_b3_size
do j = 2, nucl_num, 1
tmp_alphaj = dble(List_all_comb_b3(j,i)) * j1b_pen(j)
do k = 1, j-1, 1
tmp_alphak = dble(List_all_comb_b3(k,i)) * j1b_pen(k)
List_all_comb_b3_coef(i) += tmp_alphaj * tmp_alphak * ( (nucl_coord(j,1) - nucl_coord(k,1)) * (nucl_coord(j,1) - nucl_coord(k,1)) &
+ (nucl_coord(j,2) - nucl_coord(k,2)) * (nucl_coord(j,2) - nucl_coord(k,2)) &
+ (nucl_coord(j,3) - nucl_coord(k,3)) * (nucl_coord(j,3) - nucl_coord(k,3)) )
enddo
enddo
if(List_all_comb_b3_expo(i) .lt. 1d-10) cycle
List_all_comb_b3_coef(i) = List_all_comb_b3_coef(i) / List_all_comb_b3_expo(i)
enddo
! ---
do i = 1, List_all_comb_b3_size
facto = 1.d0
phase = 0
do j = 1, nucl_num
tmp_alphaj = dble(List_all_comb_b3(j,i))
facto *= 2.d0 / (gamma(tmp_alphaj+1.d0) * gamma(3.d0-tmp_alphaj))
phase += List_all_comb_b3(j,i)
enddo
List_all_comb_b3_coef(i) = (-1.d0)**dble(phase) * facto * dexp(-List_all_comb_b3_coef(i))
enddo
END_PROVIDER
! ---

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

@ -82,7 +82,7 @@ double precision function NAI_pol_mult_erf_ao_with1s(i_ao, j_ao, beta, B_center,
double precision, external :: NAI_pol_mult_erf_with1s, NAI_pol_mult_erf_ao
ASSERT(beta .lt. 0.d0)
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
@ -234,17 +234,17 @@ double precision function NAI_pol_mult_erf_with1s( A1_center, A2_center, power_A
double precision :: rint
! e^{-alpha1 (r - A1)^2} e^{-alpha2 (r - A2)^2} = e^{K12} e^{-alpha12 (r - A12)^2}
! 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
dist12 = 0.d0
do i = 1, 3
A12_center(i) = (alpha1 * A1_center(i) + alpha2 * A2_center(i)) * alpha12_inv
dist12 += (A1_center(i) - A2_center(i)) * (A1_center(i) - A2_center(i))
enddo
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
@ -254,19 +254,17 @@ double precision function NAI_pol_mult_erf_with1s( A1_center, A2_center, power_A
! ---
! 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
dist = 0.d0
dist_integral = 0.d0
do i = 1, 3
P_center(i) = (alpha12 * A12_center(i) + beta * B_center(i)) * p_inv
dist += (A12_center(i) - B_center(i)) * (A12_center(i) - B_center(i))
dist_integral += (P_center(i) - C_center(i)) * (P_center(i) - C_center(i))
enddo
! 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
@ -274,6 +272,12 @@ double precision function NAI_pol_mult_erf_with1s( A1_center, A2_center, power_A
return
endif
dist_integral = 0.d0
do i = 1, 3
dist_integral += (P_center(i) - C_center(i)) * (P_center(i) - C_center(i))
enddo
! ---
p_new = mu_in / dsqrt(p + mu_in * mu_in)

339
src/non_h_ints_mu/debug_integ_jmu_modif.irp.f
View File

@ -6,15 +6,24 @@ program debug_integ_jmu_modif
implicit none
my_grid_becke = .True.
!my_n_pt_r_grid = 30
!my_n_pt_a_grid = 50
my_n_pt_r_grid = 100
my_n_pt_a_grid = 170
my_n_pt_r_grid = 30
my_n_pt_a_grid = 50
!my_n_pt_r_grid = 100
!my_n_pt_a_grid = 170
touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
PROVIDE mu_erf j1b_pen
call test_grad_1_u_ij_mu()
call test_j1b_nucl()
call test_grad_j1b_nucl()
call test_lapl_j1b_nucl()
call test_list_b2()
call test_list_b3()
!call test_grad_1_u_ij_mu()
!call test_gradu_squared_u_ij_mu()
end
@ -23,16 +32,13 @@ end
subroutine test_grad_1_u_ij_mu()
implicit none
integer :: i, j, ipoint
double precision :: acc_ij, acc_tot, eps_ij, i_exc, i_num
double precision, external :: num_grad_1_u_ij_mu_x
double precision, external :: num_grad_1_u_ij_mu_y
double precision, external :: num_grad_1_u_ij_mu_z
integer :: i, j, ipoint
double precision :: acc_ij, acc_tot, eps_ij, i_exc, i_num, normalz
double precision :: integ(3)
print*, ' test_grad_1_u_ij_mu ...'
PROVIDE grad_1_u_ij_mu
! PROVIDE num_grad_1_u_ij_mu
PROVIDE grad_1_u_ij_mu
eps_ij = 1d-6
acc_tot = 0.d0
@ -41,9 +47,10 @@ subroutine test_grad_1_u_ij_mu()
do j = 1, ao_num
do i = 1, ao_num
i_exc = grad_1_u_ij_mu(i,j,ipoint,1)
!i_num = num_grad_1_u_ij_mu(i,j,ipoint,1)
i_num = num_grad_1_u_ij_mu_x(i, j, ipoint)
call num_grad_1_u_ij_mu(i, j, ipoint, integ)
i_exc = grad_1_u_ij_mu(i,j,ipoint,1)
i_num = integ(1)
acc_ij = dabs(i_exc - i_num)
if(acc_ij .gt. eps_ij) then
print *, ' problem in x part of grad_1_u_ij_mu on', i, j, ipoint
@ -52,10 +59,10 @@ subroutine test_grad_1_u_ij_mu()
print *, ' diff = ', acc_ij
endif
acc_tot += acc_ij
normalz += dabs(i_num)
i_exc = grad_1_u_ij_mu(i,j,ipoint,2)
!i_num = num_grad_1_u_ij_mu(i,j,ipoint,2)
i_num = num_grad_1_u_ij_mu_y(i, j, ipoint)
i_exc = grad_1_u_ij_mu(i,j,ipoint,2)
i_num = integ(2)
acc_ij = dabs(i_exc - i_num)
if(acc_ij .gt. eps_ij) then
print *, ' problem in y part of grad_1_u_ij_mu on', i, j, ipoint
@ -64,30 +71,316 @@ subroutine test_grad_1_u_ij_mu()
print *, ' diff = ', acc_ij
endif
acc_tot += acc_ij
normalz += dabs(i_num)
i_exc = grad_1_u_ij_mu(i,j,ipoint,3)
!i_num = num_grad_1_u_ij_mu(i,j,ipoint,3)
i_num = num_grad_1_u_ij_mu_z(i, j, ipoint)
i_exc = grad_1_u_ij_mu(i,j,ipoint,3)
i_num = integ(3)
acc_ij = dabs(i_exc - i_num)
if(acc_ij .gt. eps_ij) then
print *, ' problem in y part of grad_1_u_ij_mu on', i, j, ipoint
print *, ' problem in z part of grad_1_u_ij_mu on', i, j, ipoint
print *, ' analyt integ = ', i_exc
print *, ' numeri integ = ', i_num
print *, ' diff = ', acc_ij
endif
acc_tot += acc_ij
normalz += dabs(i_num)
enddo
enddo
enddo
acc_tot = acc_tot / dble(ao_num*ao_num*n_points_final_grid)
acc_tot = acc_tot / normalz
print*, ' normalized acc = ', acc_tot
print*, ' normalz = ', normalz
return
end subroutine test_grad_1_u_ij_mu
! ---
subroutine test_gradu_squared_u_ij_mu()
implicit none
integer :: i, j, ipoint
double precision :: acc_ij, acc_tot, eps_ij, i_exc, i_num, normalz
double precision, external :: num_gradu_squared_u_ij_mu
print*, ' test_gradu_squared_u_ij_mu ...'
PROVIDE gradu_squared_u_ij_mu
eps_ij = 1d-6
acc_tot = 0.d0
do ipoint = 1, n_points_final_grid
do j = 1, ao_num
do i = 1, ao_num
i_exc = gradu_squared_u_ij_mu(i,j,ipoint)
i_num = num_gradu_squared_u_ij_mu(i, j, ipoint)
acc_ij = dabs(i_exc - i_num)
if(acc_ij .gt. eps_ij) then
print *, ' problem in gradu_squared_u_ij_mu on', i, j, ipoint
print *, ' analyt integ = ', i_exc
print *, ' numeri integ = ', i_num
print *, ' diff = ', acc_ij
endif
acc_tot += acc_ij
normalz += dabs(i_num)
enddo
enddo
enddo
acc_tot = acc_tot / normalz
print*, ' normalized acc = ', acc_tot
print*, ' normalz = ', normalz
return
end subroutine test_gradu_squared_u_ij_mu
! ---
subroutine test_j1b_nucl()
implicit none
integer :: ipoint
double precision :: acc_ij, acc_tot, eps_ij, i_exc, i_num, normalz
double precision :: r(3)
double precision, external :: j1b_nucl
print*, ' test_j1b_nucl ...'
PROVIDE v_1b
eps_ij = 1d-7
acc_tot = 0.d0
do ipoint = 1, n_points_final_grid
r(1) = final_grid_points(1,ipoint)
r(2) = final_grid_points(2,ipoint)
r(3) = final_grid_points(3,ipoint)
i_exc = v_1b(ipoint)
i_num = j1b_nucl(r)
acc_ij = dabs(i_exc - i_num)
if(acc_ij .gt. eps_ij) then
print *, ' problem in v_1b on', ipoint
print *, ' analyt = ', i_exc
print *, ' numeri = ', i_num
print *, ' diff = ', acc_ij
endif
acc_tot += acc_ij
normalz += dabs(i_num)
enddo
acc_tot = acc_tot / normalz
print*, ' normalized acc = ', acc_tot
print*, ' normalz = ', normalz
return
end subroutine test_j1b_nucl
! ---
subroutine test_grad_j1b_nucl()
implicit none
integer :: ipoint
double precision :: acc_ij, acc_tot, eps_ij, i_exc, i_num, normalz
double precision :: r(3)
double precision, external :: grad_x_j1b_nucl
double precision, external :: grad_y_j1b_nucl
double precision, external :: grad_z_j1b_nucl
print*, ' test_grad_j1b_nucl ...'
PROVIDE v_1b_grad
eps_ij = 1d-6
acc_tot = 0.d0
do ipoint = 1, n_points_final_grid
r(1) = final_grid_points(1,ipoint)
r(2) = final_grid_points(2,ipoint)
r(3) = final_grid_points(3,ipoint)
i_exc = v_1b_grad(1,ipoint)
i_num = grad_x_j1b_nucl(r)
acc_ij = dabs(i_exc - i_num)
if(acc_ij .gt. eps_ij) then
print *, ' problem in x of v_1b_grad on', ipoint
print *, ' analyt = ', i_exc
print *, ' numeri = ', i_num
print *, ' diff = ', acc_ij
endif
i_exc = v_1b_grad(2,ipoint)
i_num = grad_y_j1b_nucl(r)
acc_ij = dabs(i_exc - i_num)
if(acc_ij .gt. eps_ij) then
print *, ' problem in y of v_1b_grad on', ipoint
print *, ' analyt = ', i_exc
print *, ' numeri = ', i_num
print *, ' diff = ', acc_ij
endif
i_exc = v_1b_grad(3,ipoint)
i_num = grad_z_j1b_nucl(r)
acc_ij = dabs(i_exc - i_num)
if(acc_ij .gt. eps_ij) then
print *, ' problem in z of v_1b_grad on', ipoint
print *, ' analyt = ', i_exc
print *, ' numeri = ', i_num
print *, ' diff = ', acc_ij
endif
acc_tot += acc_ij
normalz += dabs(i_num)
enddo
acc_tot = acc_tot / normalz
print*, ' normalized acc = ', acc_tot
print*, ' normalz = ', normalz
return
end subroutine test_grad_j1b_nucl
! ---
subroutine test_lapl_j1b_nucl()
implicit none
integer :: ipoint
double precision :: acc_ij, acc_tot, eps_ij, i_exc, i_num, normalz
double precision :: r(3)
double precision, external :: lapl_j1b_nucl
print*, ' test_lapl_j1b_nucl ...'
PROVIDE v_1b_lapl
eps_ij = 1d-5
acc_tot = 0.d0
do ipoint = 1, n_points_final_grid
r(1) = final_grid_points(1,ipoint)
r(2) = final_grid_points(2,ipoint)
r(3) = final_grid_points(3,ipoint)
i_exc = v_1b_lapl(ipoint)
i_num = lapl_j1b_nucl(r)
acc_ij = dabs(i_exc - i_num)
if(acc_ij .gt. eps_ij) then
print *, ' problem in v_1b_lapl on', ipoint
print *, ' analyt = ', i_exc
print *, ' numeri = ', i_num
print *, ' diff = ', acc_ij
endif
acc_tot += acc_ij
normalz += dabs(i_num)
enddo
acc_tot = acc_tot / normalz
print*, ' normalized acc = ', acc_tot
print*, ' normalz = ', normalz
return
end subroutine test_lapl_j1b_nucl
! ---
subroutine test_list_b2()
implicit none
integer :: ipoint
double precision :: acc_ij, acc_tot, eps_ij, i_exc, i_num, normalz
double precision :: r(3)
double precision, external :: j1b_nucl
print*, ' test_list_b2 ...'
PROVIDE v_1b_list_b2
eps_ij = 1d-7
acc_tot = 0.d0
do ipoint = 1, n_points_final_grid
r(1) = final_grid_points(1,ipoint)
r(2) = final_grid_points(2,ipoint)
r(3) = final_grid_points(3,ipoint)
i_exc = v_1b_list_b2(ipoint)
i_num = j1b_nucl(r)
acc_ij = dabs(i_exc - i_num)
if(acc_ij .gt. eps_ij) then
print *, ' problem in list_b2 on', ipoint
print *, ' analyt = ', i_exc
print *, ' numeri = ', i_num
print *, ' diff = ', acc_ij
endif
acc_tot += acc_ij
normalz += dabs(i_num)
enddo
acc_tot = acc_tot / normalz
print*, ' normalized acc = ', acc_tot
print*, ' normalz = ', normalz
return
end subroutine test_list_b2
! ---
subroutine test_list_b3()
implicit none
integer :: ipoint
double precision :: acc_ij, acc_tot, eps_ij, i_exc, i_tmp, i_num, normalz
double precision :: r(3)
double precision, external :: j1b_nucl
print*, ' test_list_b3 ...'
PROVIDE v_1b_list_b3
eps_ij = 1d-7
acc_tot = 0.d0
do ipoint = 1, n_points_final_grid
r(1) = final_grid_points(1,ipoint)
r(2) = final_grid_points(2,ipoint)
r(3) = final_grid_points(3,ipoint)
i_exc = v_1b_list_b3(ipoint)
i_tmp = j1b_nucl(r)
i_num = i_tmp * i_tmp
acc_ij = dabs(i_exc - i_num)
if(acc_ij .gt. eps_ij) then
print *, ' problem in list_b3 on', ipoint
print *, ' analyt = ', i_exc
print *, ' numeri = ', i_num
print *, ' diff = ', acc_ij
endif
acc_tot += acc_ij
normalz += dabs(i_num)
enddo
acc_tot = acc_tot / normalz
print*, ' normalized acc = ', acc_tot
print*, ' normalz = ', normalz
return
end subroutine test_list_b3
! ---

3
src/non_h_ints_mu/grad_squared.irp.f
View File

@ -27,9 +27,10 @@ BEGIN_PROVIDER [ double precision, gradu_squared_u_ij_mu, (ao_num, ao_num,n_poin
PROVIDE j1b_type j1b_pen
if(j1b_type .eq. 3) then
! v1_1b^2 \int d2 \phi_i(2) \phi_j(2) \frac{-[1 - \erf(\mu r12)]^2}{4} v2_1b^2
do ipoint = 1, n_points_final_grid
tmp = fact3_j12(ipoint)
tmp = v_1b(ipoint) * v_1b(ipoint)
do j = 1, ao_num
do i = 1, ao_num
gradu_squared_u_ij_mu(j,i,ipoint) += tmp * int2_grad1u_grad2u_j1b(i,j,ipoint)

598
src/non_h_ints_mu/j12_nucl_utils.irp.f Normal file
View File

@ -0,0 +1,598 @@
! ---
BEGIN_PROVIDER [ double precision, v_1b, (n_points_final_grid)]
implicit none
integer :: ipoint, i, j, phase
double precision :: x, y, z, dx, dy, dz
double precision :: a, d, e, fact_r
do ipoint = 1, n_points_final_grid
x = final_grid_points(1,ipoint)
y = final_grid_points(2,ipoint)
z = final_grid_points(3,ipoint)
fact_r = 1.d0
do j = 1, nucl_num
a = j1b_pen(j)
dx = x - nucl_coord(j,1)
dy = y - nucl_coord(j,2)
dz = z - nucl_coord(j,3)
d = dx*dx + dy*dy + dz*dz
e = 1.d0 - dexp(-a*d)
fact_r = fact_r * e
enddo
v_1b(ipoint) = fact_r
enddo
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, v_1b_grad, (3, n_points_final_grid)]
implicit none
integer :: ipoint, i, j, phase
double precision :: x, y, z, dx, dy, dz
double precision :: a, d, e
double precision :: fact_x, fact_y, fact_z
double precision :: ax_der, ay_der, az_der, a_expo
do ipoint = 1, n_points_final_grid
x = final_grid_points(1,ipoint)
y = final_grid_points(2,ipoint)
z = final_grid_points(3,ipoint)
fact_x = 0.d0
fact_y = 0.d0
fact_z = 0.d0
do i = 1, List_all_comb_b2_size
phase = 0
a_expo = 0.d0
ax_der = 0.d0
ay_der = 0.d0
az_der = 0.d0
do j = 1, nucl_num
a = dble(List_all_comb_b2(j,i)) * j1b_pen(j)
dx = x - nucl_coord(j,1)
dy = y - nucl_coord(j,2)
dz = z - nucl_coord(j,3)
phase += List_all_comb_b2(j,i)
a_expo += a * (dx*dx + dy*dy + dz*dz)
ax_der += a * dx
ay_der += a * dy
az_der += a * dz
enddo
e = -2.d0 * (-1.d0)**dble(phase) * dexp(-a_expo)
fact_x += e * ax_der
fact_y += e * ay_der
fact_z += e * az_der
enddo
v_1b_grad(1,ipoint) = fact_x
v_1b_grad(2,ipoint) = fact_y
v_1b_grad(3,ipoint) = fact_z
enddo
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, v_1b_lapl, (n_points_final_grid)]
implicit none
integer :: ipoint, i, j, phase
double precision :: x, y, z, dx, dy, dz
double precision :: a, d, e, b
double precision :: fact_r
double precision :: ax_der, ay_der, az_der, a_expo
do ipoint = 1, n_points_final_grid
x = final_grid_points(1,ipoint)
y = final_grid_points(2,ipoint)
z = final_grid_points(3,ipoint)
fact_r = 0.d0
do i = 1, List_all_comb_b2_size
phase = 0
b = 0.d0
a_expo = 0.d0
ax_der = 0.d0
ay_der = 0.d0
az_der = 0.d0
do j = 1, nucl_num
a = dble(List_all_comb_b2(j,i)) * j1b_pen(j)
dx = x - nucl_coord(j,1)
dy = y - nucl_coord(j,2)
dz = z - nucl_coord(j,3)
phase += List_all_comb_b2(j,i)
b += a
a_expo += a * (dx*dx + dy*dy + dz*dz)
ax_der += a * dx
ay_der += a * dy
az_der += a * dz
enddo
fact_r += (-1.d0)**dble(phase) * (-6.d0*b + 4.d0*(ax_der*ax_der + ay_der*ay_der + az_der*az_der) ) * dexp(-a_expo)
enddo
v_1b_lapl(ipoint) = fact_r
enddo
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, v_1b_list_b2, (n_points_final_grid)]
implicit none
integer :: i, ipoint
double precision :: x, y, z, coef, expo, dx, dy, dz
double precision :: fact_r
PROVIDE List_all_comb_b2_coef List_all_comb_b2_expo List_all_comb_b2_cent
do ipoint = 1, n_points_final_grid
x = final_grid_points(1,ipoint)
y = final_grid_points(2,ipoint)
z = final_grid_points(3,ipoint)
fact_r = 0.d0
do i = 1, List_all_comb_b2_size
coef = List_all_comb_b2_coef(i)
expo = List_all_comb_b2_expo(i)
dx = x - List_all_comb_b2_cent(1,i)
dy = y - List_all_comb_b2_cent(2,i)
dz = z - List_all_comb_b2_cent(3,i)
fact_r += coef * dexp(-expo * (dx*dx + dy*dy + dz*dz))
enddo
v_1b_list_b2(ipoint) = fact_r
enddo
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, v_1b_list_b3, (n_points_final_grid)]
implicit none
integer :: i, ipoint
double precision :: x, y, z, coef, expo, dx, dy, dz
double precision :: fact_r
PROVIDE List_all_comb_b3_coef List_all_comb_b3_expo List_all_comb_b3_cent
do ipoint = 1, n_points_final_grid
x = final_grid_points(1,ipoint)
y = final_grid_points(2,ipoint)
z = final_grid_points(3,ipoint)
fact_r = 0.d0
do i = 1, List_all_comb_b3_size
coef = List_all_comb_b3_coef(i)
expo = List_all_comb_b3_expo(i)
dx = x - List_all_comb_b3_cent(1,i)
dy = y - List_all_comb_b3_cent(2,i)
dz = z - List_all_comb_b3_cent(3,i)
fact_r += coef * dexp(-expo * (dx*dx + dy*dy + dz*dz))
enddo
v_1b_list_b3(ipoint) = fact_r
enddo
END_PROVIDER
! ---
double precision function jmu_modif(r1, r2)
implicit none
double precision, intent(in) :: r1(3), r2(3)
double precision, external :: j12_mu, j12_nucl
jmu_modif = j12_mu(r1, r2) * j12_nucl(r1, r2)
return
end function jmu_modif
! ---
double precision function j12_mu(r1, r2)
include 'constants.include.F'
implicit none
double precision, intent(in) :: r1(3), r2(3)
double precision :: mu_r12, r12
r12 = dsqrt( (r1(1) - r2(1)) * (r1(1) - r2(1)) &
+ (r1(2) - r2(2)) * (r1(2) - r2(2)) &
+ (r1(3) - r2(3)) * (r1(3) - r2(3)) )
mu_r12 = mu_erf * r12
j12_mu = 0.5d0 * r12 * (1.d0 - derf(mu_r12)) - inv_sq_pi_2 * dexp(-mu_r12*mu_r12) / mu_erf
return
end function j12_mu
! ---
double precision function j1b_nucl(r)
implicit none
double precision, intent(in) :: r(3)
integer :: i
double precision :: a, d, e
j1b_nucl = 1.d0
do i = 1, nucl_num
a = j1b_pen(i)
d = ( (r(1) - nucl_coord(i,1)) * (r(1) - nucl_coord(i,1)) &
+ (r(2) - nucl_coord(i,2)) * (r(2) - nucl_coord(i,2)) &
+ (r(3) - nucl_coord(i,3)) * (r(3) - nucl_coord(i,3)) )
e = 1.d0 - exp(-a*d)
j1b_nucl = j1b_nucl * e
enddo
return
end function j1b_nucl
! ---
double precision function j12_nucl(r1, r2)
implicit none
double precision, intent(in) :: r1(3), r2(3)
double precision, external :: j1b_nucl
j12_nucl = j1b_nucl(r1) * j1b_nucl(r2)
return
end function j12_nucl
! ---
! ---------------------------------------------------------------------------------------
double precision function grad_x_j1b_nucl(r)
implicit none
double precision, intent(in) :: r(3)
double precision :: r_eps(3), eps, fp, fm, delta
double precision, external :: j1b_nucl
eps = 1d-6
r_eps = r
delta = max(eps, dabs(eps*r(1)))
r_eps(1) = r_eps(1) + delta
fp = j1b_nucl(r_eps)
r_eps(1) = r_eps(1) - 2.d0 * delta
fm = j1b_nucl(r_eps)
grad_x_j1b_nucl = 0.5d0 * (fp - fm) / delta
return
end function grad_x_j1b_nucl
double precision function grad_y_j1b_nucl(r)
implicit none
double precision, intent(in) :: r(3)
double precision :: r_eps(3), eps, fp, fm, delta
double precision, external :: j1b_nucl
eps = 1d-6
r_eps = r
delta = max(eps, dabs(eps*r(2)))
r_eps(2) = r_eps(2) + delta
fp = j1b_nucl(r_eps)
r_eps(2) = r_eps(2) - 2.d0 * delta
fm = j1b_nucl(r_eps)
grad_y_j1b_nucl = 0.5d0 * (fp - fm) / delta
return
end function grad_y_j1b_nucl
double precision function grad_z_j1b_nucl(r)
implicit none
double precision, intent(in) :: r(3)
double precision :: r_eps(3), eps, fp, fm, delta
double precision, external :: j1b_nucl
eps = 1d-6
r_eps = r
delta = max(eps, dabs(eps*r(3)))
r_eps(3) = r_eps(3) + delta
fp = j1b_nucl(r_eps)
r_eps(3) = r_eps(3) - 2.d0 * delta
fm = j1b_nucl(r_eps)
grad_z_j1b_nucl = 0.5d0 * (fp - fm) / delta
return
end function grad_z_j1b_nucl
! ---------------------------------------------------------------------------------------
! ---
double precision function lapl_j1b_nucl(r)
implicit none
double precision, intent(in) :: r(3)
double precision :: r_eps(3), eps, fp, fm, delta
double precision, external :: grad_x_j1b_nucl
double precision, external :: grad_y_j1b_nucl
double precision, external :: grad_z_j1b_nucl
eps = 1d-5
r_eps = r
lapl_j1b_nucl = 0.d0
! ---
delta = max(eps, dabs(eps*r(1)))
r_eps(1) = r_eps(1) + delta
fp = grad_x_j1b_nucl(r_eps)
r_eps(1) = r_eps(1) - 2.d0 * delta
fm = grad_x_j1b_nucl(r_eps)
r_eps(1) = r_eps(1) + delta
lapl_j1b_nucl += 0.5d0 * (fp - fm) / delta
! ---
delta = max(eps, dabs(eps*r(2)))
r_eps(2) = r_eps(2) + delta
fp = grad_y_j1b_nucl(r_eps)
r_eps(2) = r_eps(2) - 2.d0 * delta
fm = grad_y_j1b_nucl(r_eps)
r_eps(2) = r_eps(2) + delta
lapl_j1b_nucl += 0.5d0 * (fp - fm) / delta
! ---
delta = max(eps, dabs(eps*r(3)))
r_eps(3) = r_eps(3) + delta
fp = grad_z_j1b_nucl(r_eps)
r_eps(3) = r_eps(3) - 2.d0 * delta
fm = grad_z_j1b_nucl(r_eps)
r_eps(3) = r_eps(3) + delta
lapl_j1b_nucl += 0.5d0 * (fp - fm) / delta
! ---
return
end function lapl_j1b_nucl
! ---
! ---------------------------------------------------------------------------------------
double precision function grad1_x_jmu_modif(r1, r2)
implicit none
double precision, intent(in) :: r1(3), r2(3)
double precision :: r1_eps(3), eps, fp, fm, delta
double precision, external :: jmu_modif
eps = 1d-7
r1_eps = r1
delta = max(eps, dabs(eps*r1(1)))
r1_eps(1) = r1_eps(1) + delta
fp = jmu_modif(r1_eps, r2)
r1_eps(1) = r1_eps(1) - 2.d0 * delta
fm = jmu_modif(r1_eps, r2)
grad1_x_jmu_modif = 0.5d0 * (fp - fm) / delta
return
end function grad1_x_jmu_modif
double precision function grad1_y_jmu_modif(r1, r2)
implicit none
double precision, intent(in) :: r1(3), r2(3)
double precision :: r1_eps(3), eps, fp, fm, delta
double precision, external :: jmu_modif
eps = 1d-7
r1_eps = r1
delta = max(eps, dabs(eps*r1(2)))
r1_eps(2) = r1_eps(2) + delta
fp = jmu_modif(r1_eps, r2)
r1_eps(2) = r1_eps(2) - 2.d0 * delta
fm = jmu_modif(r1_eps, r2)
grad1_y_jmu_modif = 0.5d0 * (fp - fm) / delta
return
end function grad1_y_jmu_modif
double precision function grad1_z_jmu_modif(r1, r2)
implicit none
double precision, intent(in) :: r1(3), r2(3)
double precision :: r1_eps(3), eps, fp, fm, delta
double precision, external :: jmu_modif
eps = 1d-7
r1_eps = r1
delta = max(eps, dabs(eps*r1(3)))
r1_eps(3) = r1_eps(3) + delta
fp = jmu_modif(r1_eps, r2)
r1_eps(3) = r1_eps(3) - 2.d0 * delta
fm = jmu_modif(r1_eps, r2)
grad1_z_jmu_modif = 0.5d0 * (fp - fm) / delta
return
end function grad1_z_jmu_modif
! ---------------------------------------------------------------------------------------
! ---
! ---------------------------------------------------------------------------------------
double precision function grad1_x_j12_mu_num(r1, r2)
implicit none
double precision, intent(in) :: r1(3), r2(3)
double precision :: r1_eps(3), eps, fp, fm, delta
double precision, external :: j12_mu
eps = 1d-7
r1_eps = r1
delta = max(eps, dabs(eps*r1(1)))
r1_eps(1) = r1_eps(1) + delta
fp = j12_mu(r1_eps, r2)
r1_eps(1) = r1_eps(1) - 2.d0 * delta
fm = j12_mu(r1_eps, r2)
grad1_x_j12_mu_num = 0.5d0 * (fp - fm) / delta
return
end function grad1_x_j12_mu_num
double precision function grad1_y_j12_mu_num(r1, r2)
implicit none
double precision, intent(in) :: r1(3), r2(3)
double precision :: r1_eps(3), eps, fp, fm, delta
double precision, external :: j12_mu
eps = 1d-7
r1_eps = r1
delta = max(eps, dabs(eps*r1(2)))
r1_eps(2) = r1_eps(2) + delta
fp = j12_mu(r1_eps, r2)
r1_eps(2) = r1_eps(2) - 2.d0 * delta
fm = j12_mu(r1_eps, r2)
grad1_y_j12_mu_num = 0.5d0 * (fp - fm) / delta
return
end function grad1_y_j12_mu_num
double precision function grad1_z_j12_mu_num(r1, r2)
implicit none
double precision, intent(in) :: r1(3), r2(3)
double precision :: r1_eps(3), eps, fp, fm, delta
double precision, external :: j12_mu
eps = 1d-7
r1_eps = r1
delta = max(eps, dabs(eps*r1(3)))
r1_eps(3) = r1_eps(3) + delta
fp = j12_mu(r1_eps, r2)
r1_eps(3) = r1_eps(3) - 2.d0 * delta
fm = j12_mu(r1_eps, r2)
grad1_z_j12_mu_num = 0.5d0 * (fp - fm) / delta
return
end function grad1_z_j12_mu_num
! ---------------------------------------------------------------------------------------
! ---
! ---------------------------------------------------------------------------------------
double precision function grad1_x_j12_mu_exc(r1, r2)
implicit none
double precision, intent(in) :: r1(3), r2(3)
double precision :: r12
grad1_x_j12_mu_exc = 0.d0
r12 = dsqrt( (r1(1) - r2(1)) * (r1(1) - r2(1)) &
+ (r1(2) - r2(2)) * (r1(2) - r2(2)) &
+ (r1(3) - r2(3)) * (r1(3) - r2(3)) )
if(r12 .lt. 1d-10) return
grad1_x_j12_mu_exc = 0.5d0 * (1.d0 - derf(mu_erf * r12)) * (r1(1) - r2(1)) / r12
return
end function grad1_x_j12_mu_exc
double precision function grad1_y_j12_mu_exc(r1, r2)
implicit none
double precision, intent(in) :: r1(3), r2(3)
double precision :: r12
grad1_y_j12_mu_exc = 0.d0
r12 = dsqrt( (r1(1) - r2(1)) * (r1(1) - r2(1)) &
+ (r1(2) - r2(2)) * (r1(2) - r2(2)) &
+ (r1(3) - r2(3)) * (r1(3) - r2(3)) )
if(r12 .lt. 1d-10) return
grad1_y_j12_mu_exc = 0.5d0 * (1.d0 - derf(mu_erf * r12)) * (r1(2) - r2(2)) / r12
return
end function grad1_y_j12_mu_exc
double precision function grad1_z_j12_mu_exc(r1, r2)
implicit none
double precision, intent(in) :: r1(3), r2(3)
double precision :: r12
grad1_z_j12_mu_exc = 0.d0
r12 = dsqrt( (r1(1) - r2(1)) * (r1(1) - r2(1)) &
+ (r1(2) - r2(2)) * (r1(2) - r2(2)) &
+ (r1(3) - r2(3)) * (r1(3) - r2(3)) )
if(r12 .lt. 1d-10) return
grad1_z_j12_mu_exc = 0.5d0 * (1.d0 - derf(mu_erf * r12)) * (r1(3) - r2(3)) / r12
return
end function grad1_z_j12_mu_exc
! ---------------------------------------------------------------------------------------
! ---

266
src/non_h_ints_mu/jmu_modif.irp.f
View File

@ -1,266 +0,0 @@
! ---
double precision function jmu_modif(r1, r2)
implicit none
double precision, intent(in) :: r1(3), r2(3)
double precision, external :: j12_mu, j12_nucl
jmu_modif = j12_mu(r1, r2) * j12_nucl(r1, r2)
return
end function jmu_modif
! ---
double precision function j12_mu(r1, r2)
include 'constants.include.F'
implicit none
double precision, intent(in) :: r1(3), r2(3)
double precision :: mu_r12, r12
r12 = dsqrt( (r1(1) - r2(1)) * (r1(1) - r2(1)) &
+ (r1(2) - r2(2)) * (r1(2) - r2(2)) &
+ (r1(3) - r2(3)) * (r1(3) - r2(3)) )
mu_r12 = mu_erf * r12
j12_mu = 0.5d0 * r12 * (1.d0 - derf(mu_r12)) - inv_sq_pi_2 * dexp(-mu_r12*mu_r12) / mu_erf
return
end function j12_mu
! ---
double precision function j12_nucl(r1, r2)
implicit none
double precision, intent(in) :: r1(3), r2(3)
integer :: i, j
double precision :: a1, d1, e1, a2, d2, e2
j12_nucl = 1.d0
do i = 1, nucl_num
a1 = j1b_pen(i)
d1 = ( (r1(1) - nucl_coord(i,1)) * (r1(1) - nucl_coord(i,1)) &
+ (r1(2) - nucl_coord(i,2)) * (r1(2) - nucl_coord(i,2)) &
+ (r1(3) - nucl_coord(i,3)) * (r1(3) - nucl_coord(i,3)) )
e1 = 1.d0 - exp(-a1*d1)
do j = 1, nucl_num
a2 = j1b_pen(j)
d2 = ( (r2(1) - nucl_coord(j,1)) * (r2(1) - nucl_coord(j,1)) &
+ (r2(2) - nucl_coord(j,2)) * (r2(2) - nucl_coord(j,2)) &
+ (r2(3) - nucl_coord(j,3)) * (r2(3) - nucl_coord(j,3)) )
e2 = 1.d0 - exp(-a2*d2)
j12_nucl = j12_nucl * e1 * e2
enddo
enddo
return
end function j12_nucl
! ---
! ---------------------------------------------------------------------------------------
double precision function grad1_x_jmu_modif(r1, r2)
implicit none
double precision, intent(in) :: r1(3), r2(3)
double precision :: r1_eps(3), eps, fp, fm, delta
double precision, external :: jmu_modif
eps = 1d-7
r1_eps = r1
delta = max(eps, dabs(eps*r1(1)))
r1_eps(1) = r1_eps(1) + delta
fp = jmu_modif(r1_eps, r2)
r1_eps(1) = r1_eps(1) - 2.d0 * delta
fm = jmu_modif(r1_eps, r2)
grad1_x_jmu_modif = 0.5d0 * (fp - fm) / delta
return
end function grad1_x_jmu_modif
double precision function grad1_y_jmu_modif(r1, r2)
implicit none
double precision, intent(in) :: r1(3), r2(3)
double precision :: r1_eps(3), eps, fp, fm, delta
double precision, external :: jmu_modif
eps = 1d-7
r1_eps = r1
delta = max(eps, dabs(eps*r1(2)))
r1_eps(2) = r1_eps(2) + delta
fp = jmu_modif(r1_eps, r2)
r1_eps(2) = r1_eps(2) - 2.d0 * delta
fm = jmu_modif(r1_eps, r2)
grad1_y_jmu_modif = 0.5d0 * (fp - fm) / delta
return
end function grad1_y_jmu_modif
double precision function grad1_z_jmu_modif(r1, r2)
implicit none
double precision, intent(in) :: r1(3), r2(3)
double precision :: r1_eps(3), eps, fp, fm, delta
double precision, external :: jmu_modif
eps = 1d-7
r1_eps = r1
delta = max(eps, dabs(eps*r1(3)))
r1_eps(3) = r1_eps(3) + delta
fp = jmu_modif(r1_eps, r2)
r1_eps(3) = r1_eps(3) - 2.d0 * delta
fm = jmu_modif(r1_eps, r2)
grad1_z_jmu_modif = 0.5d0 * (fp - fm) / delta
return
end function grad1_z_jmu_modif
! ---------------------------------------------------------------------------------------
! ---
! ---------------------------------------------------------------------------------------
double precision function grad1_x_j12_mu_num(r1, r2)
implicit none
double precision, intent(in) :: r1(3), r2(3)
double precision :: r1_eps(3), eps, fp, fm, delta
double precision, external :: j12_mu
eps = 1d-7
r1_eps = r1
delta = max(eps, dabs(eps*r1(1)))
r1_eps(1) = r1_eps(1) + delta
fp = j12_mu(r1_eps, r2)
r1_eps(1) = r1_eps(1) - 2.d0 * delta
fm = j12_mu(r1_eps, r2)
grad1_x_j12_mu_num = 0.5d0 * (fp - fm) / delta
return
end function grad1_x_j12_mu_num
double precision function grad1_y_j12_mu_num(r1, r2)
implicit none
double precision, intent(in) :: r1(3), r2(3)
double precision :: r1_eps(3), eps, fp, fm, delta
double precision, external :: j12_mu
eps = 1d-7
r1_eps = r1
delta = max(eps, dabs(eps*r1(2)))
r1_eps(2) = r1_eps(2) + delta
fp = j12_mu(r1_eps, r2)
r1_eps(2) = r1_eps(2) - 2.d0 * delta
fm = j12_mu(r1_eps, r2)
grad1_y_j12_mu_num = 0.5d0 * (fp - fm) / delta
return
end function grad1_y_j12_mu_num
double precision function grad1_z_j12_mu_num(r1, r2)
implicit none
double precision, intent(in) :: r1(3), r2(3)
double precision :: r1_eps(3), eps, fp, fm, delta
double precision, external :: j12_mu
eps = 1d-7
r1_eps = r1
delta = max(eps, dabs(eps*r1(3)))
r1_eps(3) = r1_eps(3) + delta
fp = j12_mu(r1_eps, r2)
r1_eps(3) = r1_eps(3) - 2.d0 * delta
fm = j12_mu(r1_eps, r2)
grad1_z_j12_mu_num = 0.5d0 * (fp - fm) / delta
return
end function grad1_z_j12_mu_num
! ---------------------------------------------------------------------------------------
! ---
! ---------------------------------------------------------------------------------------
double precision function grad1_x_j12_mu_exc(r1, r2)
implicit none
double precision, intent(in) :: r1(3), r2(3)
double precision :: r12
grad1_x_j12_mu_exc = 0.d0
r12 = dsqrt( (r1(1) - r2(1)) * (r1(1) - r2(1)) &
+ (r1(2) - r2(2)) * (r1(2) - r2(2)) &
+ (r1(3) - r2(3)) * (r1(3) - r2(3)) )
if(r12 .lt. 1d-10) return
grad1_x_j12_mu_exc = 0.5d0 * (1.d0 - derf(mu_erf * r12)) * (r1(1) - r2(1)) / r12
return
end function grad1_x_j12_mu_exc
double precision function grad1_y_j12_mu_exc(r1, r2)
implicit none
double precision, intent(in) :: r1(3), r2(3)
double precision :: r12
grad1_y_j12_mu_exc = 0.d0
r12 = dsqrt( (r1(1) - r2(1)) * (r1(1) - r2(1)) &
+ (r1(2) - r2(2)) * (r1(2) - r2(2)) &
+ (r1(3) - r2(3)) * (r1(3) - r2(3)) )
if(r12 .lt. 1d-10) return
grad1_y_j12_mu_exc = 0.5d0 * (1.d0 - derf(mu_erf * r12)) * (r1(2) - r2(2)) / r12
return
end function grad1_y_j12_mu_exc
double precision function grad1_z_j12_mu_exc(r1, r2)
implicit none
double precision, intent(in) :: r1(3), r2(3)
double precision :: r12
grad1_z_j12_mu_exc = 0.d0
r12 = dsqrt( (r1(1) - r2(1)) * (r1(1) - r2(1)) &
+ (r1(2) - r2(2)) * (r1(2) - r2(2)) &
+ (r1(3) - r2(3)) * (r1(3) - r2(3)) )
if(r12 .lt. 1d-10) return
grad1_z_j12_mu_exc = 0.5d0 * (1.d0 - derf(mu_erf * r12)) * (r1(3) - r2(3)) / r12
return
end function grad1_z_j12_mu_exc
! ---------------------------------------------------------------------------------------
! ---

86
src/non_h_ints_mu/new_grad_tc.irp.f
View File

@ -1,81 +1,3 @@
! ---
BEGIN_PROVIDER [ double precision, fact1_j12, ( n_points_final_grid)]
&BEGIN_PROVIDER [ double precision, fact2_j12, (3, n_points_final_grid)]
&BEGIN_PROVIDER [ double precision, fact3_j12, ( n_points_final_grid)]
implicit none
integer :: ipoint, i, j, phase
double precision :: x, y, z, dx, dy, dz
double precision :: a, d, e, fact_r, fact_r_sq
double precision :: fact_x, fact_y, fact_z
double precision :: ax_der, ay_der, az_der, a_expo
do ipoint = 1, n_points_final_grid
x = final_grid_points(1,ipoint)
y = final_grid_points(2,ipoint)
z = final_grid_points(3,ipoint)
! ---
fact_r = 1.d0
fact_r_sq = 1.d0
do j = 1, nucl_num
a = j1b_pen(j)
dx = x - nucl_coord(j,1)
dy = y - nucl_coord(j,2)
dz = z - nucl_coord(j,3)
d = x*x + y*y + z*z
e = 1.d0 - dexp(-a*d)
fact_r = fact_r * e
fact_r_sq = fact_r_sq * e * e
enddo
fact1_j12(ipoint) = fact_r
fact3_j12(ipoint) = fact_r_sq
! ---
fact_x = 0.d0
fact_y = 0.d0
fact_z = 0.d0
do i = 1, List_all_comb_b2_size
phase = 0
a_expo = 0.d0
ax_der = 0.d0
ay_der = 0.d0
az_der = 0.d0
do j = 1, nucl_num
a = dble(List_all_comb_b2(j,i)) * j1b_pen(j)
dx = x - nucl_coord(j,1)
dy = y - nucl_coord(j,2)
dz = z - nucl_coord(j,3)
phase += List_all_comb_b2(j,i)
a_expo += a * (dx*dx + dy*dy + dz*dz)
ax_der += a * dx
ay_der += a * dy
az_der += a * dz
enddo
e = -2.d0 * (-1.d0)**dble(phase) * dexp(-a_expo)
fact_x += e * ax_der
fact_y += e * ay_der
fact_z += e * az_der
enddo
fact2_j12(1,ipoint) = fact_x
fact2_j12(2,ipoint) = fact_y
fact2_j12(3,ipoint) = fact_z
! ---
enddo
END_PROVIDER
! ---
@ -103,10 +25,10 @@ BEGIN_PROVIDER [ double precision, grad_1_u_ij_mu, (ao_num, ao_num, n_points_fin
y = final_grid_points(2,ipoint)
z = final_grid_points(3,ipoint)
tmp0 = fact1_j12(ipoint)
tmp_x = fact2_j12(1,ipoint)
tmp_y = fact2_j12(2,ipoint)
tmp_z = fact2_j12(3,ipoint)
tmp0 = v_1b (ipoint)
tmp_x = v_1b_grad(1,ipoint)
tmp_y = v_1b_grad(2,ipoint)
tmp_z = v_1b_grad(3,ipoint)
do j = 1, ao_num
do i = 1, ao_num

165
src/non_h_ints_mu/numerical_integ.irp.f
View File

@ -4,132 +4,129 @@
!
! \int dr2 [-1 * \grad_r1 u(r1,r2)] \phi_i(r2) \phi_j(r2) x 1s_j1b(r2)
!
BEGIN_PROVIDER [ double precision, num_grad_1_u_ij_mu, (ao_num, ao_num, n_points_final_grid, 3)]
implicit none
integer :: i, j, ipoint, jpoint
double precision :: tmp, r1(3), r2(3)
double precision, external :: ao_value
double precision, external :: j12_nucl
double precision, external :: grad1_x_j12_mu_num, grad1_x_j12_mu_exc
double precision, external :: grad1_y_j12_mu_num, grad1_y_j12_mu_exc
double precision, external :: grad1_z_j12_mu_num, grad1_z_j12_mu_exc
num_grad_1_u_ij_mu = 0.d0
do j = 1, ao_num
do i = 1, ao_num
do ipoint = 1, n_points_final_grid
r1(1) = final_grid_points(1,ipoint)
r1(2) = final_grid_points(2,ipoint)
r1(3) = final_grid_points(3,ipoint)
do jpoint = 1, n_points_final_grid
r2(1) = final_grid_points(1,jpoint)
r2(2) = final_grid_points(2,jpoint)
r2(3) = final_grid_points(3,jpoint)
tmp = ao_value(i, r2) * ao_value(j, r2) * j12_nucl(r1, r2) * final_weight_at_r_vector(jpoint)
num_grad_1_u_ij_mu(i,j,ipoint,1) += tmp * (-1.d0 * grad1_x_j12_mu_exc(r1, r2))
num_grad_1_u_ij_mu(i,j,ipoint,2) += tmp * (-1.d0 * grad1_y_j12_mu_exc(r1, r2))
num_grad_1_u_ij_mu(i,j,ipoint,3) += tmp * (-1.d0 * grad1_z_j12_mu_exc(r1, r2))
enddo
enddo
enddo
enddo
END_PROVIDER
!BEGIN_PROVIDER [ double precision, num_grad_1_u_ij_mu, (ao_num, ao_num, n_points_final_grid, 3)]
!
! implicit none
!
! integer :: i, j, ipoint, jpoint
! double precision :: tmp, r1(3), r2(3)
!
! double precision, external :: ao_value
! double precision, external :: j12_nucl
! double precision, external :: grad1_x_j12_mu_num, grad1_x_j12_mu_exc
! double precision, external :: grad1_y_j12_mu_num, grad1_y_j12_mu_exc
! double precision, external :: grad1_z_j12_mu_num, grad1_z_j12_mu_exc
!
! num_grad_1_u_ij_mu = 0.d0
!
! do j = 1, ao_num
! do i = 1, ao_num
!
! do ipoint = 1, n_points_final_grid
! r1(1) = final_grid_points(1,ipoint)
! r1(2) = final_grid_points(2,ipoint)
! r1(3) = final_grid_points(3,ipoint)
!
! do jpoint = 1, n_points_final_grid
! r2(1) = final_grid_points(1,jpoint)
! r2(2) = final_grid_points(2,jpoint)
! r2(3) = final_grid_points(3,jpoint)
! tmp = ao_value(i, r2) * ao_value(j, r2) * j12_nucl(r1, r2) * final_weight_at_r_vector(jpoint)
!
! num_grad_1_u_ij_mu(i,j,ipoint,1) += tmp * (-1.d0 * grad1_x_j12_mu_exc(r1, r2))
! num_grad_1_u_ij_mu(i,j,ipoint,2) += tmp * (-1.d0 * grad1_y_j12_mu_exc(r1, r2))
! num_grad_1_u_ij_mu(i,j,ipoint,3) += tmp * (-1.d0 * grad1_z_j12_mu_exc(r1, r2))
! enddo
!
! enddo
! enddo
! enddo
!
!END_PROVIDER
! ---
double precision function num_grad_1_u_ij_mu_x(i, j, ipoint)
subroutine num_grad_1_u_ij_mu(i, j, ipoint, integ)
implicit none
integer, intent(in) :: i, j, ipoint
integer :: jpoint
double precision :: tmp, r1(3), r2(3)
double precision, external :: ao_value
double precision, external :: j12_nucl
double precision, external :: grad1_x_j12_mu_num, grad1_x_j12_mu_exc
num_grad_1_u_ij_mu_x = 0.d0
integer, intent(in) :: i, j, ipoint
double precision, intent(out) :: integ(3)
integer :: jpoint
double precision :: tmp, r1(3), r2(3)
double precision :: tmp_x, tmp_y, tmp_z
double precision, external :: ao_value
double precision, external :: j12_nucl
double precision, external :: grad1_x_j12_mu_num, grad1_x_j12_mu_exc
double precision, external :: grad1_y_j12_mu_num, grad1_y_j12_mu_exc
double precision, external :: grad1_z_j12_mu_num, grad1_z_j12_mu_exc
r1(1) = final_grid_points(1,ipoint)
r1(2) = final_grid_points(2,ipoint)
r1(3) = final_grid_points(3,ipoint)
tmp_x = 0.d0
tmp_y = 0.d0
tmp_z = 0.d0
do jpoint = 1, n_points_final_grid
r2(1) = final_grid_points(1,jpoint)
r2(2) = final_grid_points(2,jpoint)
r2(3) = final_grid_points(3,jpoint)
tmp = ao_value(i, r2) * ao_value(j, r2) * j12_nucl(r1, r2) * final_weight_at_r_vector(jpoint)
num_grad_1_u_ij_mu_x += tmp * (-1.d0 * grad1_x_j12_mu_exc(r1, r2))
tmp_x += tmp * (-1.d0 * grad1_x_j12_mu_exc(r1, r2))
tmp_y += tmp * (-1.d0 * grad1_y_j12_mu_exc(r1, r2))
tmp_z += tmp * (-1.d0 * grad1_z_j12_mu_exc(r1, r2))
enddo
end function num_grad_1_u_ij_mu_x
integ(1) = tmp_x
integ(2) = tmp_y
integ(3) = tmp_z
return
end subroutine num_grad_1_u_ij_mu
! ---
double precision function num_grad_1_u_ij_mu_y(i, j, ipoint)
double precision function num_gradu_squared_u_ij_mu(i, j, ipoint)
implicit none
integer, intent(in) :: i, j, ipoint
integer :: jpoint
double precision :: tmp, r1(3), r2(3)
double precision :: tmp, r1(3), r2(3), r12
double precision :: tmp_x, tmp_y, tmp_z, tmp1, tmp2
double precision, external :: ao_value
double precision, external :: j12_nucl
double precision, external :: grad1_y_j12_mu_num, grad1_y_j12_mu_exc
num_grad_1_u_ij_mu_y = 0.d0
r1(1) = final_grid_points(1,ipoint)
r1(2) = final_grid_points(2,ipoint)
r1(3) = final_grid_points(3,ipoint)
num_gradu_squared_u_ij_mu = 0.d0
do jpoint = 1, n_points_final_grid
r2(1) = final_grid_points(1,jpoint)
r2(2) = final_grid_points(2,jpoint)
r2(3) = final_grid_points(3,jpoint)
tmp = ao_value(i, r2) * ao_value(j, r2) * j12_nucl(r1, r2) * final_weight_at_r_vector(jpoint)
tmp_x = r1(1) - r2(1)
tmp_y = r1(2) - r2(2)
tmp_z = r1(3) - r2(3)
r12 = dsqrt( tmp_x*tmp_x + tmp_y*tmp_y + tmp_z*tmp_z )
tmp1 = 1.d0 - derf(mu_erf * r12)
tmp2 = j12_nucl(r1, r2)
tmp = -0.25d0 * tmp1 * tmp1 * tmp2 * tmp2 * ao_value(i, r2) * ao_value(j, r2) * final_weight_at_r_vector(jpoint)
num_grad_1_u_ij_mu_y += tmp * (-1.d0 * grad1_y_j12_mu_exc(r1, r2))
num_gradu_squared_u_ij_mu += tmp
enddo
end function num_grad_1_u_ij_mu_y
return
end function num_gradu_squared_u_ij_mu
! ---
double precision function num_grad_1_u_ij_mu_z(i, j, ipoint)
implicit none
integer, intent(in) :: i, j, ipoint
integer :: jpoint
double precision :: tmp, r1(3), r2(3)
double precision, external :: ao_value
double precision, external :: j12_nucl
double precision, external :: grad1_z_j12_mu_num, grad1_z_j12_mu_exc
num_grad_1_u_ij_mu_z = 0.d0
r1(1) = final_grid_points(1,ipoint)
r1(2) = final_grid_points(2,ipoint)
r1(3) = final_grid_points(3,ipoint)
do jpoint = 1, n_points_final_grid
r2(1) = final_grid_points(1,jpoint)
r2(2) = final_grid_points(2,jpoint)
r2(3) = final_grid_points(3,jpoint)
tmp = ao_value(i, r2) * ao_value(j, r2) * j12_nucl(r1, r2) * final_weight_at_r_vector(jpoint)
num_grad_1_u_ij_mu_z += tmp * (-1.d0 * grad1_z_j12_mu_exc(r1, r2))
enddo
end function num_grad_1_u_ij_mu_z
! ---