mirror of
https://github.com/QuantumPackage/qp2.git
synced 2024-08-30 04:33:41 +02:00
integrals over r2 tested
This commit is contained in:
parent
dbbae1f990
commit
b13a315cc1
@ -1,7 +1,7 @@
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [ double precision, int2_grad1u_grad2u_j1b, (ao_num, ao_num, n_points_final_grid)]
|
||||
BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2, (ao_num, ao_num, n_points_final_grid)]
|
||||
|
||||
BEGIN_DOC
|
||||
!
|
||||
@ -21,7 +21,7 @@ BEGIN_PROVIDER [ double precision, int2_grad1u_grad2u_j1b, (ao_num, ao_num, n_po
|
||||
provide mu_erf final_grid_points j1b_pen
|
||||
call wall_time(wall0)
|
||||
|
||||
int2_grad1u_grad2u_j1b = 0.d0
|
||||
int2_grad1u2_grad2u2_j1b2 = 0.d0
|
||||
|
||||
!$OMP PARALLEL DEFAULT (NONE) &
|
||||
!$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, &
|
||||
@ -30,12 +30,13 @@ BEGIN_PROVIDER [ double precision, int2_grad1u_grad2u_j1b, (ao_num, ao_num, n_po
|
||||
!$OMP final_grid_points, n_max_fit_slat, &
|
||||
!$OMP expo_gauss_1_erf_x_2, coef_gauss_1_erf_x_2, &
|
||||
!$OMP List_all_comb_b3_coef, List_all_comb_b3_expo, &
|
||||
!$OMP List_all_comb_b3_cent, int2_grad1u_grad2u_j1b)
|
||||
!$OMP List_all_comb_b3_cent, int2_grad1u2_grad2u2_j1b2)
|
||||
|
||||
allocate( tmp(ao_num,ao_num,n_points_final_grid) )
|
||||
tmp = 0.d0
|
||||
|
||||
!$OMP DO
|
||||
!do ipoint = 1, 10
|
||||
do ipoint = 1, n_points_final_grid
|
||||
do i = 1, ao_num
|
||||
do j = i, ao_num
|
||||
@ -69,7 +70,7 @@ BEGIN_PROVIDER [ double precision, int2_grad1u_grad2u_j1b, (ao_num, ao_num, n_po
|
||||
do ipoint = 1, n_points_final_grid
|
||||
do i = 1, ao_num
|
||||
do j = i, ao_num
|
||||
int2_grad1u_grad2u_j1b(j,i,ipoint) += tmp(j,i,ipoint)
|
||||
int2_grad1u2_grad2u2_j1b2(j,i,ipoint) += tmp(j,i,ipoint)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
@ -81,19 +82,19 @@ BEGIN_PROVIDER [ double precision, int2_grad1u_grad2u_j1b, (ao_num, ao_num, n_po
|
||||
do ipoint = 1, n_points_final_grid
|
||||
do i = 1, ao_num
|
||||
do j = 1, i-1
|
||||
int2_grad1u_grad2u_j1b(j,i,ipoint) = int2_grad1u_grad2u_j1b(i,j,ipoint)
|
||||
int2_grad1u2_grad2u2_j1b2(j,i,ipoint) = int2_grad1u2_grad2u2_j1b2(i,j,ipoint)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
call wall_time(wall1)
|
||||
print*, ' wall time for int2_grad1u_grad2u_j1b', wall1 - wall0
|
||||
print*, ' wall time for int2_grad1u2_grad2u2_j1b2', wall1 - wall0
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [ double precision, int2_u2_j1b, (ao_num, ao_num, n_points_final_grid)]
|
||||
BEGIN_PROVIDER [ double precision, int2_u2_j1b2, (ao_num, ao_num, n_points_final_grid)]
|
||||
|
||||
BEGIN_DOC
|
||||
!
|
||||
@ -113,7 +114,7 @@ BEGIN_PROVIDER [ double precision, int2_u2_j1b, (ao_num, ao_num, n_points_final_
|
||||
provide mu_erf final_grid_points j1b_pen
|
||||
call wall_time(wall0)
|
||||
|
||||
int2_u2_j1b = 0.d0
|
||||
int2_u2_j1b2 = 0.d0
|
||||
|
||||
!$OMP PARALLEL DEFAULT (NONE) &
|
||||
!$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, &
|
||||
@ -122,12 +123,13 @@ BEGIN_PROVIDER [ double precision, int2_u2_j1b, (ao_num, ao_num, n_points_final_
|
||||
!$OMP final_grid_points, n_max_fit_slat, &
|
||||
!$OMP expo_gauss_j_mu_x_2, coef_gauss_j_mu_x_2, &
|
||||
!$OMP List_all_comb_b3_coef, List_all_comb_b3_expo, &
|
||||
!$OMP List_all_comb_b3_cent, int2_u2_j1b)
|
||||
!$OMP List_all_comb_b3_cent, int2_u2_j1b2)
|
||||
|
||||
allocate( tmp(ao_num,ao_num,n_points_final_grid) )
|
||||
tmp = 0.d0
|
||||
|
||||
!$OMP DO
|
||||
!do ipoint = 1, 10
|
||||
do ipoint = 1, n_points_final_grid
|
||||
do i = 1, ao_num
|
||||
do j = i, ao_num
|
||||
@ -161,7 +163,7 @@ BEGIN_PROVIDER [ double precision, int2_u2_j1b, (ao_num, ao_num, n_points_final_
|
||||
do ipoint = 1, n_points_final_grid
|
||||
do i = 1, ao_num
|
||||
do j = i, ao_num
|
||||
int2_u2_j1b(j,i,ipoint) += tmp(j,i,ipoint)
|
||||
int2_u2_j1b2(j,i,ipoint) += tmp(j,i,ipoint)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
@ -173,13 +175,13 @@ BEGIN_PROVIDER [ double precision, int2_u2_j1b, (ao_num, ao_num, n_points_final_
|
||||
do ipoint = 1, n_points_final_grid
|
||||
do i = 1, ao_num
|
||||
do j = 1, i-1
|
||||
int2_u2_j1b(j,i,ipoint) = int2_u2_j1b(i,j,ipoint)
|
||||
int2_u2_j1b2(j,i,ipoint) = int2_u2_j1b2(i,j,ipoint)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
call wall_time(wall1)
|
||||
print*, ' wall time for int2_u2_j1b', wall1 - wall0
|
||||
print*, ' wall time for int2_u2_j1b2', wall1 - wall0
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
@ -297,7 +299,7 @@ END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b, (ao_num, ao_num, n_points_final_grid)]
|
||||
BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2, (ao_num, ao_num, n_points_final_grid)]
|
||||
|
||||
BEGIN_DOC
|
||||
!
|
||||
@ -317,7 +319,7 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b, (ao_num, ao_num, n_points_
|
||||
provide mu_erf final_grid_points j1b_pen
|
||||
call wall_time(wall0)
|
||||
|
||||
int2_u_grad1u_j1b = 0.d0
|
||||
int2_u_grad1u_j1b2 = 0.d0
|
||||
|
||||
!$OMP PARALLEL DEFAULT (NONE) &
|
||||
!$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, &
|
||||
@ -327,7 +329,7 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b, (ao_num, ao_num, n_points_
|
||||
!$OMP final_grid_points, n_max_fit_slat, &
|
||||
!$OMP expo_gauss_j_mu_1_erf, coef_gauss_j_mu_1_erf, &
|
||||
!$OMP List_all_comb_b3_coef, List_all_comb_b3_expo, &
|
||||
!$OMP List_all_comb_b3_cent, int2_u_grad1u_j1b)
|
||||
!$OMP List_all_comb_b3_cent, int2_u_grad1u_j1b2)
|
||||
|
||||
allocate( tmp(ao_num,ao_num,n_points_final_grid) )
|
||||
tmp = 0.d0
|
||||
@ -380,7 +382,7 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b, (ao_num, ao_num, n_points_
|
||||
do ipoint = 1, n_points_final_grid
|
||||
do i = 1, ao_num
|
||||
do j = i, ao_num
|
||||
int2_u_grad1u_j1b(j,i,ipoint) += tmp(j,i,ipoint)
|
||||
int2_u_grad1u_j1b2(j,i,ipoint) += tmp(j,i,ipoint)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
@ -392,13 +394,13 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b, (ao_num, ao_num, n_points_
|
||||
do ipoint = 1, n_points_final_grid
|
||||
do i = 1, ao_num
|
||||
do j = 1, i-1
|
||||
int2_u_grad1u_j1b(j,i,ipoint) = int2_u_grad1u_j1b(i,j,ipoint)
|
||||
int2_u_grad1u_j1b2(j,i,ipoint) = int2_u_grad1u_j1b2(i,j,ipoint)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
call wall_time(wall1)
|
||||
print*, ' wall time for int2_u_grad1u_j1b', wall1 - wall0
|
||||
print*, ' wall time for int2_u_grad1u_j1b2', wall1 - wall0
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
@ -33,6 +33,7 @@ BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_j1b, (ao_num, ao_num, n_po
|
||||
tmp = 0.d0
|
||||
|
||||
!$OMP DO
|
||||
!do ipoint = 1, 10
|
||||
do ipoint = 1, n_points_final_grid
|
||||
do i = 1, ao_num
|
||||
do j = i, ao_num
|
||||
@ -141,6 +142,7 @@ BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_tmp_j1b, (3, ao_num, ao_
|
||||
tmp = 0.d0
|
||||
|
||||
!$OMP DO
|
||||
!do ipoint = 1, 10
|
||||
do ipoint = 1, n_points_final_grid
|
||||
do i = 1, ao_num
|
||||
do j = i, ao_num
|
||||
@ -235,6 +237,7 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b, (ao_num, ao_num, n_points_
|
||||
tmp = 0.d0
|
||||
|
||||
!$OMP DO
|
||||
!do ipoint = 1, 10
|
||||
do ipoint = 1, n_points_final_grid
|
||||
do i = 1, ao_num
|
||||
do j = i, ao_num
|
||||
|
@ -103,7 +103,7 @@ double precision function NAI_pol_mult_erf_ao_with1s(i_ao, j_ao, beta, B_center,
|
||||
alpha2 = ao_expo_ordered_transp(j,j_ao)
|
||||
|
||||
coef12 = ao_coef_normalized_ordered_transp(j,j_ao) * ao_coef_normalized_ordered_transp(i,i_ao)
|
||||
if(coef12 .lt. 1d-14) cycle
|
||||
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 )
|
||||
|
@ -7,20 +7,32 @@ program debug_integ_jmu_modif
|
||||
|
||||
my_grid_becke = .True.
|
||||
|
||||
my_n_pt_r_grid = 30
|
||||
my_n_pt_a_grid = 50
|
||||
!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 = 150
|
||||
my_n_pt_a_grid = 194
|
||||
touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
|
||||
|
||||
PROVIDE mu_erf j1b_pen
|
||||
|
||||
call test_j1b_nucl()
|
||||
call test_grad_j1b_nucl()
|
||||
call test_lapl_j1b_nucl()
|
||||
!call test_j1b_nucl()
|
||||
!call test_grad_j1b_nucl()
|
||||
!call test_lapl_j1b_nucl()
|
||||
|
||||
call test_list_b2()
|
||||
call test_list_b3()
|
||||
!call test_list_b2()
|
||||
!call test_list_b3()
|
||||
|
||||
!call test_fit_u()
|
||||
call test_fit_u2()
|
||||
!call test_fit_ugradu()
|
||||
|
||||
!call test_v_ij_u_cst_mu_j1b()
|
||||
!call test_v_ij_erf_rk_cst_mu_j1b()
|
||||
!call test_x_v_ij_erf_rk_cst_mu_j1b()
|
||||
!call test_int2_u2_j1b2()
|
||||
!call test_int2_grad1u2_grad2u2_j1b2()
|
||||
|
||||
!call test_grad_1_u_ij_mu()
|
||||
!call test_gradu_squared_u_ij_mu()
|
||||
@ -29,6 +41,252 @@ end
|
||||
|
||||
! ---
|
||||
|
||||
subroutine test_v_ij_u_cst_mu_j1b()
|
||||
|
||||
implicit none
|
||||
integer :: i, j, ipoint
|
||||
double precision :: acc_ij, acc_tot, eps_ij, i_exc, i_num, normalz
|
||||
double precision, external :: num_v_ij_u_cst_mu_j1b
|
||||
|
||||
print*, ' test_v_ij_u_cst_mu_j1b ...'
|
||||
|
||||
PROVIDE v_ij_u_cst_mu_j1b
|
||||
|
||||
eps_ij = 1d-8
|
||||
acc_tot = 0.d0
|
||||
|
||||
!do ipoint = 1, 10
|
||||
do ipoint = 1, n_points_final_grid
|
||||
do j = 1, ao_num
|
||||
do i = 1, ao_num
|
||||
|
||||
i_exc = v_ij_u_cst_mu_j1b(i,j,ipoint)
|
||||
i_num = num_v_ij_u_cst_mu_j1b(i,j,ipoint)
|
||||
acc_ij = dabs(i_exc - i_num)
|
||||
if(acc_ij .gt. eps_ij) then
|
||||
print *, ' problem in v_ij_u_cst_mu_j1b 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_v_ij_u_cst_mu_j1b
|
||||
|
||||
! ---
|
||||
|
||||
subroutine test_v_ij_erf_rk_cst_mu_j1b()
|
||||
|
||||
implicit none
|
||||
integer :: i, j, ipoint
|
||||
double precision :: acc_ij, acc_tot, eps_ij, i_exc, i_num, normalz
|
||||
double precision, external :: num_v_ij_erf_rk_cst_mu_j1b
|
||||
|
||||
print*, ' test_v_ij_erf_rk_cst_mu_j1b ...'
|
||||
|
||||
PROVIDE v_ij_erf_rk_cst_mu_j1b
|
||||
|
||||
eps_ij = 1d-8
|
||||
acc_tot = 0.d0
|
||||
|
||||
!do ipoint = 1, 10
|
||||
do ipoint = 1, n_points_final_grid
|
||||
do j = 1, ao_num
|
||||
do i = 1, ao_num
|
||||
|
||||
i_exc = v_ij_erf_rk_cst_mu_j1b(i,j,ipoint)
|
||||
i_num = num_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint)
|
||||
acc_ij = dabs(i_exc - i_num)
|
||||
if(acc_ij .gt. eps_ij) then
|
||||
print *, ' problem in v_ij_erf_rk_cst_mu_j1b 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_v_ij_erf_rk_cst_mu_j1b
|
||||
|
||||
! ---
|
||||
|
||||
subroutine test_x_v_ij_erf_rk_cst_mu_j1b()
|
||||
|
||||
implicit none
|
||||
integer :: i, j, ipoint
|
||||
double precision :: acc_ij, acc_tot, eps_ij, i_exc, i_num, normalz
|
||||
double precision :: integ(3)
|
||||
|
||||
print*, ' test_x_v_ij_erf_rk_cst_mu_j1b ...'
|
||||
|
||||
PROVIDE x_v_ij_erf_rk_cst_mu_j1b
|
||||
|
||||
eps_ij = 1d-8
|
||||
acc_tot = 0.d0
|
||||
|
||||
!do ipoint = 1, 10
|
||||
do ipoint = 1, n_points_final_grid
|
||||
do j = 1, ao_num
|
||||
do i = 1, ao_num
|
||||
|
||||
call num_x_v_ij_erf_rk_cst_mu_j1b(i, j, ipoint, integ)
|
||||
|
||||
i_exc = x_v_ij_erf_rk_cst_mu_j1b(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 x_v_ij_erf_rk_cst_mu_j1b 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)
|
||||
|
||||
i_exc = x_v_ij_erf_rk_cst_mu_j1b(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 x_v_ij_erf_rk_cst_mu_j1b 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)
|
||||
|
||||
i_exc = x_v_ij_erf_rk_cst_mu_j1b(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 z part of x_v_ij_erf_rk_cst_mu_j1b 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_x_v_ij_erf_rk_cst_mu_j1b
|
||||
|
||||
! ---
|
||||
|
||||
subroutine test_int2_u2_j1b2()
|
||||
|
||||
implicit none
|
||||
integer :: i, j, ipoint
|
||||
double precision :: acc_ij, acc_tot, eps_ij, i_exc, i_num, normalz
|
||||
double precision, external :: num_int2_u2_j1b2
|
||||
|
||||
print*, ' test_int2_u2_j1b2 ...'
|
||||
|
||||
PROVIDE int2_u2_j1b2
|
||||
|
||||
eps_ij = 1d-8
|
||||
acc_tot = 0.d0
|
||||
|
||||
!do ipoint = 1, 10
|
||||
do ipoint = 1, n_points_final_grid
|
||||
do j = 1, ao_num
|
||||
do i = 1, ao_num
|
||||
|
||||
i_exc = int2_u2_j1b2(i,j,ipoint)
|
||||
i_num = num_int2_u2_j1b2(i,j,ipoint)
|
||||
acc_ij = dabs(i_exc - i_num)
|
||||
if(acc_ij .gt. eps_ij) then
|
||||
print *, ' problem in int2_u2_j1b2 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_int2_u2_j1b2
|
||||
|
||||
! ---
|
||||
|
||||
subroutine test_int2_grad1u2_grad2u2_j1b2()
|
||||
|
||||
implicit none
|
||||
integer :: i, j, ipoint
|
||||
double precision :: acc_ij, acc_tot, eps_ij, i_exc, i_num, normalz
|
||||
double precision, external :: num_int2_grad1u2_grad2u2_j1b2
|
||||
|
||||
print*, ' test_int2_grad1u2_grad2u2_j1b2 ...'
|
||||
|
||||
PROVIDE int2_grad1u2_grad2u2_j1b2
|
||||
|
||||
eps_ij = 1d-8
|
||||
acc_tot = 0.d0
|
||||
|
||||
!do ipoint = 1, 10
|
||||
do ipoint = 1, n_points_final_grid
|
||||
do j = 1, ao_num
|
||||
do i = 1, ao_num
|
||||
|
||||
i_exc = int2_grad1u2_grad2u2_j1b2(i,j,ipoint)
|
||||
i_num = num_int2_grad1u2_grad2u2_j1b2(i,j,ipoint)
|
||||
acc_ij = dabs(i_exc - i_num)
|
||||
if(acc_ij .gt. eps_ij) then
|
||||
print *, ' problem in int2_grad1u2_grad2u2_j1b2 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_int2_grad1u2_grad2u2_j1b2
|
||||
|
||||
! ---
|
||||
|
||||
subroutine test_grad_1_u_ij_mu()
|
||||
|
||||
implicit none
|
||||
@ -384,3 +642,206 @@ end subroutine test_list_b3
|
||||
|
||||
! ---
|
||||
|
||||
subroutine test_fit_ugradu()
|
||||
|
||||
implicit none
|
||||
|
||||
integer :: ipoint, i
|
||||
double precision :: i_exc, i_fit, i_num, x2
|
||||
double precision :: r1(3), r2(3), grad(3)
|
||||
double precision :: eps_ij, acc_tot, acc_ij, normalz, coef, expo
|
||||
|
||||
double precision, external :: j12_mu
|
||||
|
||||
print*, ' test_fit_ugradu ...'
|
||||
|
||||
eps_ij = 1d-7
|
||||
acc_tot = 0.d0
|
||||
|
||||
r2 = 0.d0
|
||||
|
||||
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)
|
||||
x2 = r1(1) * r1(1) + r1(2) * r1(2) + r1(3) * r1(3)
|
||||
if(x2 .lt. 1d-10) cycle
|
||||
|
||||
i_fit = 0.d0
|
||||
do i = 1, n_max_fit_slat
|
||||
expo = expo_gauss_j_mu_1_erf(i)
|
||||
coef = coef_gauss_j_mu_1_erf(i)
|
||||
i_fit += coef * dexp(-expo*x2)
|
||||
enddo
|
||||
i_fit = i_fit / dsqrt(x2)
|
||||
|
||||
call grad1_j12_mu_exc(r1, r2, grad)
|
||||
|
||||
! ---
|
||||
|
||||
i_exc = j12_mu(r1, r2) * grad(1)
|
||||
i_num = i_fit * r1(1)
|
||||
acc_ij = dabs(i_exc - i_num)
|
||||
if(acc_ij .gt. eps_ij) then
|
||||
print *, ' problem on x in test_fit_ugradu on', ipoint
|
||||
print *, ' analyt = ', i_exc
|
||||
print *, ' numeri = ', i_num
|
||||
print *, ' diff = ', acc_ij
|
||||
endif
|
||||
acc_tot += acc_ij
|
||||
normalz += dabs(i_exc)
|
||||
|
||||
! ---
|
||||
|
||||
i_exc = j12_mu(r1, r2) * grad(2)
|
||||
i_num = i_fit * r1(2)
|
||||
acc_ij = dabs(i_exc - i_num)
|
||||
if(acc_ij .gt. eps_ij) then
|
||||
print *, ' problem on y in test_fit_ugradu on', ipoint
|
||||
print *, ' analyt = ', i_exc
|
||||
print *, ' numeri = ', i_num
|
||||
print *, ' diff = ', acc_ij
|
||||
endif
|
||||
acc_tot += acc_ij
|
||||
normalz += dabs(i_exc)
|
||||
|
||||
! ---
|
||||
|
||||
i_exc = j12_mu(r1, r2) * grad(3)
|
||||
i_num = i_fit * r1(3)
|
||||
acc_ij = dabs(i_exc - i_num)
|
||||
if(acc_ij .gt. eps_ij) then
|
||||
print *, ' problem on z in test_fit_ugradu on', ipoint
|
||||
print *, ' analyt = ', i_exc
|
||||
print *, ' numeri = ', i_num
|
||||
print *, ' diff = ', acc_ij
|
||||
endif
|
||||
acc_tot += acc_ij
|
||||
normalz += dabs(i_exc)
|
||||
|
||||
! ---
|
||||
|
||||
enddo
|
||||
|
||||
acc_tot = acc_tot / normalz
|
||||
print*, ' normalized acc = ', acc_tot
|
||||
print*, ' normalz = ', normalz
|
||||
|
||||
return
|
||||
end subroutine test_fit_ugradu
|
||||
|
||||
! ---
|
||||
|
||||
subroutine test_fit_u()
|
||||
|
||||
implicit none
|
||||
|
||||
integer :: ipoint, i
|
||||
double precision :: i_exc, i_fit, i_num, x2
|
||||
double precision :: r1(3), r2(3)
|
||||
double precision :: eps_ij, acc_tot, acc_ij, normalz, coef, expo
|
||||
|
||||
double precision, external :: j12_mu
|
||||
|
||||
print*, ' test_fit_u ...'
|
||||
|
||||
eps_ij = 1d-7
|
||||
acc_tot = 0.d0
|
||||
|
||||
r2 = 0.d0
|
||||
|
||||
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)
|
||||
x2 = r1(1) * r1(1) + r1(2) * r1(2) + r1(3) * r1(3)
|
||||
if(x2 .lt. 1d-10) cycle
|
||||
|
||||
i_fit = 0.d0
|
||||
do i = 1, n_max_fit_slat
|
||||
expo = expo_gauss_j_mu_x(i)
|
||||
coef = coef_gauss_j_mu_x(i)
|
||||
i_fit += coef * dexp(-expo*x2)
|
||||
enddo
|
||||
|
||||
i_exc = j12_mu(r1, r2)
|
||||
i_num = i_fit
|
||||
acc_ij = dabs(i_exc - i_num)
|
||||
if(acc_ij .gt. eps_ij) then
|
||||
print *, ' problem in test_fit_u on', ipoint
|
||||
print *, ' analyt = ', i_exc
|
||||
print *, ' numeri = ', i_num
|
||||
print *, ' diff = ', acc_ij
|
||||
endif
|
||||
|
||||
acc_tot += acc_ij
|
||||
normalz += dabs(i_exc)
|
||||
enddo
|
||||
|
||||
acc_tot = acc_tot / normalz
|
||||
print*, ' normalized acc = ', acc_tot
|
||||
print*, ' normalz = ', normalz
|
||||
|
||||
return
|
||||
end subroutine test_fit_u
|
||||
|
||||
! ---
|
||||
|
||||
subroutine test_fit_u2()
|
||||
|
||||
implicit none
|
||||
|
||||
integer :: ipoint, i
|
||||
double precision :: i_exc, i_fit, i_num, x2
|
||||
double precision :: r1(3), r2(3)
|
||||
double precision :: eps_ij, acc_tot, acc_ij, normalz, coef, expo
|
||||
|
||||
double precision, external :: j12_mu
|
||||
|
||||
print*, ' test_fit_u2 ...'
|
||||
|
||||
eps_ij = 1d-7
|
||||
acc_tot = 0.d0
|
||||
|
||||
r2 = 0.d0
|
||||
|
||||
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)
|
||||
x2 = r1(1) * r1(1) + r1(2) * r1(2) + r1(3) * r1(3)
|
||||
if(x2 .lt. 1d-10) cycle
|
||||
|
||||
i_fit = 0.d0
|
||||
do i = 1, n_max_fit_slat
|
||||
expo = expo_gauss_j_mu_x_2(i)
|
||||
coef = coef_gauss_j_mu_x_2(i)
|
||||
i_fit += coef * dexp(-expo*x2)
|
||||
enddo
|
||||
|
||||
i_exc = j12_mu(r1, r2) * j12_mu(r1, r2)
|
||||
i_num = i_fit
|
||||
acc_ij = dabs(i_exc - i_num)
|
||||
if(acc_ij .gt. eps_ij) then
|
||||
print *, ' problem in test_fit_u2 on', ipoint
|
||||
print *, ' analyt = ', i_exc
|
||||
print *, ' numeri = ', i_num
|
||||
print *, ' diff = ', acc_ij
|
||||
endif
|
||||
|
||||
acc_tot += acc_ij
|
||||
normalz += dabs(i_exc)
|
||||
enddo
|
||||
|
||||
acc_tot = acc_tot / normalz
|
||||
print*, ' normalized acc = ', acc_tot
|
||||
print*, ' normalz = ', normalz
|
||||
|
||||
return
|
||||
end subroutine test_fit_u2
|
||||
|
||||
! ---
|
||||
|
||||
|
@ -70,12 +70,20 @@ END_PROVIDER
|
||||
integer :: i
|
||||
double precision :: tmp
|
||||
double precision :: expos(n_max_fit_slat), alpha, beta
|
||||
double precision :: alpha_opt, beta_opt
|
||||
|
||||
!alpha_opt = 2.d0 * expo_j_xmu(1)
|
||||
!beta_opt = 2.d0 * expo_j_xmu(2)
|
||||
|
||||
! direct opt
|
||||
alpha_opt = 3.52751759d0
|
||||
beta_opt = 1.26214809d0
|
||||
|
||||
tmp = 0.25d0 / (mu_erf * mu_erf * dacos(-1.d0))
|
||||
|
||||
alpha = 2.d0 * expo_j_xmu(1) * mu_erf
|
||||
alpha = alpha_opt * mu_erf
|
||||
call expo_fit_slater_gam(alpha, expos)
|
||||
beta = 2.d0 * expo_j_xmu(2) * mu_erf * mu_erf
|
||||
beta = beta_opt * mu_erf * mu_erf
|
||||
|
||||
do i = 1, n_max_fit_slat
|
||||
expo_gauss_j_mu_x_2(i) = expos(i) + beta
|
||||
@ -101,12 +109,20 @@ END_PROVIDER
|
||||
integer :: i
|
||||
double precision :: tmp
|
||||
double precision :: expos(n_max_fit_slat), alpha, beta
|
||||
double precision :: alpha_opt, beta_opt
|
||||
|
||||
!alpha_opt = expo_j_xmu(1) + expo_gauss_1_erf_x(1)
|
||||
!beta_opt = expo_j_xmu(2) + expo_gauss_1_erf_x(2)
|
||||
|
||||
! direct opt
|
||||
alpha_opt = 2.87875632d0
|
||||
beta_opt = 1.34801003d0
|
||||
|
||||
tmp = -0.25d0 / (mu_erf * dsqrt(dacos(-1.d0)))
|
||||
|
||||
alpha = (expo_j_xmu(1) + expo_gauss_1_erf_x(1)) * mu_erf
|
||||
alpha = alpha_opt * mu_erf
|
||||
call expo_fit_slater_gam(alpha, expos)
|
||||
beta = (expo_j_xmu(2) + expo_gauss_1_erf_x(2)) * mu_erf * mu_erf
|
||||
beta = beta_opt * mu_erf * mu_erf
|
||||
|
||||
do i = 1, n_max_fit_slat
|
||||
expo_gauss_j_mu_1_erf(i) = expos(i) + beta
|
||||
@ -162,8 +178,8 @@ double precision function j_mu_fit_gauss(x)
|
||||
j_mu_fit_gauss = 0.d0
|
||||
do i = 1, n_max_fit_slat
|
||||
alpha = expo_gauss_j_mu_x(i)
|
||||
coef = coef_gauss_j_mu_x(i)
|
||||
j_mu_fit_gauss += coef_gauss_j_mu_x(i) * dexp(-expo_gauss_j_mu_x(i)*x*x)
|
||||
coef = coef_gauss_j_mu_x(i)
|
||||
j_mu_fit_gauss += coef * dexp(-alpha*x*x)
|
||||
enddo
|
||||
|
||||
end
|
||||
|
@ -33,7 +33,7 @@ BEGIN_PROVIDER [ double precision, gradu_squared_u_ij_mu, (ao_num, ao_num,n_poin
|
||||
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)
|
||||
gradu_squared_u_ij_mu(j,i,ipoint) += tmp * int2_grad1u2_grad2u2_j1b2(i,j,ipoint)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
@ -237,6 +237,30 @@ end function j12_mu
|
||||
|
||||
! ---
|
||||
|
||||
double precision function j12_mu_gauss(r1, r2)
|
||||
|
||||
implicit none
|
||||
double precision, intent(in) :: r1(3), r2(3)
|
||||
integer :: i
|
||||
double precision :: r12, coef, expo
|
||||
|
||||
r12 = (r1(1) - r2(1)) * (r1(1) - r2(1)) &
|
||||
+ (r1(2) - r2(2)) * (r1(2) - r2(2)) &
|
||||
+ (r1(3) - r2(3)) * (r1(3) - r2(3))
|
||||
|
||||
j12_mu_gauss = 0.d0
|
||||
do i = 1, n_max_fit_slat
|
||||
expo = expo_gauss_j_mu_x(i)
|
||||
coef = coef_gauss_j_mu_x(i)
|
||||
|
||||
j12_mu_gauss += coef * dexp(-expo*r12)
|
||||
enddo
|
||||
|
||||
return
|
||||
end function j12_mu_gauss
|
||||
|
||||
! ---
|
||||
|
||||
double precision function j1b_nucl(r)
|
||||
|
||||
implicit none
|
||||
@ -535,63 +559,30 @@ end function grad1_z_j12_mu_num
|
||||
|
||||
! ---
|
||||
|
||||
! ---------------------------------------------------------------------------------------
|
||||
|
||||
double precision function grad1_x_j12_mu_exc(r1, r2)
|
||||
subroutine grad1_j12_mu_exc(r1, r2, grad)
|
||||
|
||||
implicit none
|
||||
double precision, intent(in) :: r1(3), r2(3)
|
||||
double precision :: r12
|
||||
double precision, intent(in) :: r1(3), r2(3)
|
||||
double precision, intent(out) :: grad(3)
|
||||
double precision :: dx, dy, dz, r12, tmp
|
||||
|
||||
grad1_x_j12_mu_exc = 0.d0
|
||||
grad = 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)) )
|
||||
dx = r1(1) - r2(1)
|
||||
dy = r1(2) - r2(2)
|
||||
dz = r1(3) - r2(3)
|
||||
|
||||
r12 = dsqrt( dx * dx + dy * dy + dz * dz )
|
||||
if(r12 .lt. 1d-10) return
|
||||
|
||||
grad1_x_j12_mu_exc = 0.5d0 * (1.d0 - derf(mu_erf * r12)) * (r1(1) - r2(1)) / r12
|
||||
tmp = 0.5d0 * (1.d0 - derf(mu_erf * r12)) / r12
|
||||
|
||||
grad(1) = tmp * dx
|
||||
grad(2) = tmp * dy
|
||||
grad(3) = tmp * dz
|
||||
|
||||
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
|
||||
|
||||
! ---------------------------------------------------------------------------------------
|
||||
end subroutine grad1_j12_mu_exc
|
||||
|
||||
! ---
|
||||
|
||||
|
@ -10,13 +10,10 @@
|
||||
! implicit none
|
||||
!
|
||||
! integer :: i, j, ipoint, jpoint
|
||||
! double precision :: tmp, r1(3), r2(3)
|
||||
! double precision :: tmp, r1(3), r2(3), grad(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
|
||||
!
|
||||
@ -34,9 +31,11 @@
|
||||
! 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))
|
||||
! call grad1_j12_mu_exc(r1, r2, grad)
|
||||
!
|
||||
! num_grad_1_u_ij_mu(i,j,ipoint,1) += tmp * (-1.d0 * grad(1))
|
||||
! num_grad_1_u_ij_mu(i,j,ipoint,2) += tmp * (-1.d0 * grad(2))
|
||||
! num_grad_1_u_ij_mu(i,j,ipoint,3) += tmp * (-1.d0 * grad(3))
|
||||
! enddo
|
||||
!
|
||||
! enddo
|
||||
@ -47,6 +46,249 @@
|
||||
|
||||
! ---
|
||||
|
||||
double precision function num_v_ij_u_cst_mu_j1b(i, j, ipoint)
|
||||
|
||||
BEGIN_DOC
|
||||
!
|
||||
! \int dr2 u12 \phi_i(r2) \phi_j(r2) x v_1b(r2)
|
||||
!
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
|
||||
integer, intent(in) :: i, j, ipoint
|
||||
|
||||
integer :: jpoint
|
||||
double precision :: r1(3), r2(3)
|
||||
|
||||
double precision, external :: ao_value
|
||||
double precision, external :: j12_mu, j1b_nucl, j12_mu_gauss
|
||||
|
||||
r1(1) = final_grid_points(1,ipoint)
|
||||
r1(2) = final_grid_points(2,ipoint)
|
||||
r1(3) = final_grid_points(3,ipoint)
|
||||
|
||||
num_v_ij_u_cst_mu_j1b = 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)
|
||||
|
||||
num_v_ij_u_cst_mu_j1b += ao_value(i, r2) * ao_value(j, r2) * j12_mu_gauss(r1, r2) * j1b_nucl(r2) * final_weight_at_r_vector(jpoint)
|
||||
enddo
|
||||
|
||||
return
|
||||
end function num_v_ij_u_cst_mu_j1b
|
||||
|
||||
! ---
|
||||
|
||||
double precision function num_int2_u2_j1b2(i, j, ipoint)
|
||||
|
||||
BEGIN_DOC
|
||||
!
|
||||
! \int dr2 u12^2 \phi_i(r2) \phi_j(r2) x v_1b(r2)^2
|
||||
!
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
|
||||
integer, intent(in) :: i, j, ipoint
|
||||
|
||||
integer :: jpoint, i_fit
|
||||
double precision :: r1(3), r2(3)
|
||||
double precision :: dx, dy, dz, r12, x2, tmp1, tmp2, tmp3, coef, expo
|
||||
|
||||
double precision, external :: ao_value
|
||||
double precision, external :: j1b_nucl
|
||||
|
||||
r1(1) = final_grid_points(1,ipoint)
|
||||
r1(2) = final_grid_points(2,ipoint)
|
||||
r1(3) = final_grid_points(3,ipoint)
|
||||
|
||||
num_int2_u2_j1b2 = 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)
|
||||
dx = r1(1) - r2(1)
|
||||
dy = r1(2) - r2(2)
|
||||
dz = r1(3) - r2(3)
|
||||
x2 = dx * dx + dy * dy + dz * dz
|
||||
r12 = dsqrt(x2)
|
||||
|
||||
tmp1 = j1b_nucl(r2)
|
||||
tmp2 = tmp1 * tmp1 * ao_value(i, r2) * ao_value(j, r2) * final_weight_at_r_vector(jpoint)
|
||||
|
||||
tmp3 = 0.d0
|
||||
do i_fit = 1, n_max_fit_slat
|
||||
expo = expo_gauss_j_mu_x_2(i_fit)
|
||||
coef = coef_gauss_j_mu_x_2(i_fit)
|
||||
|
||||
tmp3 += coef * dexp(-expo*x2)
|
||||
enddo
|
||||
|
||||
num_int2_u2_j1b2 += tmp2 * tmp3
|
||||
enddo
|
||||
|
||||
return
|
||||
end function num_int2_u2_j1b2
|
||||
|
||||
! ---
|
||||
|
||||
double precision function num_int2_grad1u2_grad2u2_j1b2(i, j, ipoint)
|
||||
|
||||
BEGIN_DOC
|
||||
!
|
||||
! \int dr2 \frac{-[erf(mu r12) -1]^2}{4} \phi_i(r2) \phi_j(r2) x v_1b(r2)^2
|
||||
!
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
|
||||
integer, intent(in) :: i, j, ipoint
|
||||
|
||||
integer :: jpoint, i_fit
|
||||
double precision :: r1(3), r2(3)
|
||||
double precision :: dx, dy, dz, r12, x2, tmp1, tmp2, tmp3, coef, expo
|
||||
|
||||
double precision, external :: ao_value
|
||||
double precision, external :: j1b_nucl
|
||||
|
||||
r1(1) = final_grid_points(1,ipoint)
|
||||
r1(2) = final_grid_points(2,ipoint)
|
||||
r1(3) = final_grid_points(3,ipoint)
|
||||
|
||||
num_int2_grad1u2_grad2u2_j1b2 = 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)
|
||||
dx = r1(1) - r2(1)
|
||||
dy = r1(2) - r2(2)
|
||||
dz = r1(3) - r2(3)
|
||||
x2 = dx * dx + dy * dy + dz * dz
|
||||
r12 = dsqrt(x2)
|
||||
|
||||
tmp1 = j1b_nucl(r2)
|
||||
tmp2 = tmp1 * tmp1 * ao_value(i, r2) * ao_value(j, r2) * final_weight_at_r_vector(jpoint)
|
||||
|
||||
tmp3 = 0.d0
|
||||
do i_fit = 1, n_max_fit_slat
|
||||
expo = expo_gauss_1_erf_x_2(i_fit)
|
||||
coef = coef_gauss_1_erf_x_2(i_fit)
|
||||
|
||||
tmp3 += coef * dexp(-expo*x2)
|
||||
enddo
|
||||
tmp3 = -0.25d0 * tmp3
|
||||
|
||||
num_int2_grad1u2_grad2u2_j1b2 += tmp2 * tmp3
|
||||
enddo
|
||||
|
||||
return
|
||||
end function num_int2_grad1u2_grad2u2_j1b2
|
||||
|
||||
! ---
|
||||
|
||||
double precision function num_v_ij_erf_rk_cst_mu_j1b(i, j, ipoint)
|
||||
|
||||
BEGIN_DOC
|
||||
!
|
||||
! \int dr2 [erf(mu r12) -1]/r12 \phi_i(r2) \phi_j(r2) x v_1b(r2)
|
||||
!
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
|
||||
integer, intent(in) :: i, j, ipoint
|
||||
|
||||
integer :: jpoint
|
||||
double precision :: r1(3), r2(3)
|
||||
double precision :: dx, dy, dz, r12, tmp1, tmp2
|
||||
|
||||
double precision, external :: ao_value
|
||||
double precision, external :: j1b_nucl
|
||||
|
||||
r1(1) = final_grid_points(1,ipoint)
|
||||
r1(2) = final_grid_points(2,ipoint)
|
||||
r1(3) = final_grid_points(3,ipoint)
|
||||
|
||||
num_v_ij_erf_rk_cst_mu_j1b = 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)
|
||||
dx = r1(1) - r2(1)
|
||||
dy = r1(2) - r2(2)
|
||||
dz = r1(3) - r2(3)
|
||||
r12 = dsqrt( dx * dx + dy * dy + dz * dz )
|
||||
if(r12 .lt. 1d-10) cycle
|
||||
|
||||
tmp1 = (derf(mu_erf * r12) - 1.d0) / r12
|
||||
tmp2 = tmp1 * ao_value(i, r2) * ao_value(j, r2) * j1b_nucl(r2) * final_weight_at_r_vector(jpoint)
|
||||
|
||||
num_v_ij_erf_rk_cst_mu_j1b += tmp2
|
||||
enddo
|
||||
|
||||
return
|
||||
end function num_v_ij_erf_rk_cst_mu_j1b
|
||||
|
||||
! ---
|
||||
|
||||
subroutine num_x_v_ij_erf_rk_cst_mu_j1b(i, j, ipoint, integ)
|
||||
|
||||
BEGIN_DOC
|
||||
!
|
||||
! \int dr2 [erf(mu r12) -1]/r12 \phi_i(r2) \phi_j(r2) x v_1b(r2) x r2
|
||||
!
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
|
||||
integer, intent(in) :: i, j, ipoint
|
||||
double precision, intent(out) :: integ(3)
|
||||
|
||||
integer :: jpoint
|
||||
double precision :: r1(3), r2(3), grad(3)
|
||||
double precision :: dx, dy, dz, r12, tmp1, tmp2
|
||||
double precision :: tmp_x, tmp_y, tmp_z
|
||||
|
||||
double precision, external :: ao_value
|
||||
double precision, external :: j1b_nucl
|
||||
|
||||
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)
|
||||
dx = r1(1) - r2(1)
|
||||
dy = r1(2) - r2(2)
|
||||
dz = r1(3) - r2(3)
|
||||
r12 = dsqrt( dx * dx + dy * dy + dz * dz )
|
||||
if(r12 .lt. 1d-10) cycle
|
||||
|
||||
tmp1 = (derf(mu_erf * r12) - 1.d0) / r12
|
||||
tmp2 = tmp1 * ao_value(i, r2) * ao_value(j, r2) * j1b_nucl(r2) * final_weight_at_r_vector(jpoint)
|
||||
|
||||
tmp_x += tmp2 * r2(1)
|
||||
tmp_y += tmp2 * r2(2)
|
||||
tmp_z += tmp2 * r2(3)
|
||||
enddo
|
||||
|
||||
integ(1) = tmp_x
|
||||
integ(2) = tmp_y
|
||||
integ(3) = tmp_z
|
||||
|
||||
return
|
||||
end subroutine num_x_v_ij_erf_rk_cst_mu_j1b
|
||||
|
||||
! ---
|
||||
|
||||
subroutine num_grad_1_u_ij_mu(i, j, ipoint, integ)
|
||||
|
||||
implicit none
|
||||
@ -55,14 +297,11 @@ subroutine num_grad_1_u_ij_mu(i, j, ipoint, integ)
|
||||
double precision, intent(out) :: integ(3)
|
||||
|
||||
integer :: jpoint
|
||||
double precision :: tmp, r1(3), r2(3)
|
||||
double precision :: tmp, r1(3), r2(3), grad(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)
|
||||
@ -77,9 +316,11 @@ subroutine num_grad_1_u_ij_mu(i, j, ipoint, integ)
|
||||
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 += 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))
|
||||
call grad1_j12_mu_exc(r1, r2, grad)
|
||||
|
||||
tmp_x += tmp * (-1.d0 * grad(1))
|
||||
tmp_y += tmp * (-1.d0 * grad(2))
|
||||
tmp_z += tmp * (-1.d0 * grad(3))
|
||||
enddo
|
||||
|
||||
integ(1) = tmp_x
|
||||
|
Loading…
Reference in New Issue
Block a user