mirror of
https://github.com/QuantumPackage/qp2.git
synced 2024-08-24 21:21:47 +02:00
added grad2_jmu_manu.irp.f grad_lapl_jmu_manu.irp.f listj1b_sorted.irp.f
This commit is contained in:
parent
354ba6cb28
commit
5bd7b7ca6b
222
src/ao_many_one_e_ints/grad2_jmu_manu.irp.f
Normal file
222
src/ao_many_one_e_ints/grad2_jmu_manu.irp.f
Normal file
@ -0,0 +1,222 @@
|
||||
|
||||
BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2_test, (ao_num, ao_num, n_points_final_grid)]
|
||||
|
||||
BEGIN_DOC
|
||||
!
|
||||
! int dr2 phi_i(r2) phi_j(r2) 1s_j1b(r2)^2 u_12^mu [\grad_1 u_12^mu]
|
||||
!
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
integer :: i, j, ipoint, i_1s, i_fit
|
||||
double precision :: r(3), int_fit, expo_fit, coef_fit, coef_tmp
|
||||
double precision :: coef, beta, B_center(3), dist
|
||||
double precision :: alpha_1s, alpha_1s_inv, centr_1s(3), expo_coef_1s, tmp
|
||||
double precision :: wall0, wall1
|
||||
double precision, external :: NAI_pol_mult_erf_ao_with1s
|
||||
double precision :: j12_mu_r12
|
||||
double precision :: sigma_ij,dist_ij_ipoint,dsqpi_3_2
|
||||
dsqpi_3_2 = (dacos(-1.d0))**(3/2)
|
||||
|
||||
provide mu_erf final_grid_points j1b_pen ao_overlap_abs
|
||||
call wall_time(wall0)
|
||||
|
||||
|
||||
int2_u_grad1u_j1b2_test = 0.d0
|
||||
|
||||
!$OMP PARALLEL DEFAULT (NONE) &
|
||||
!$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, &
|
||||
!$OMP coef_fit, expo_fit, int_fit, tmp, alpha_1s, dist, &
|
||||
!$OMP sigma_ij, beta_ij, factor_ij_1s,center_ij_1s, dist_ij_ipoint, &
|
||||
!$OMP alpha_1s_inv, centr_1s, expo_coef_1s, coef_tmp) &
|
||||
!$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b3_size, &
|
||||
!$OMP final_grid_points, n_max_fit_slat, &
|
||||
!$OMP expo_gauss_j_mu_1_erf, coef_gauss_j_mu_1_erf, &
|
||||
!$OMP ao_prod_dist_grid, ao_prod_sigma, ao_overlap_abs_grid,ao_prod_center,dsqpi_3_2, &
|
||||
!$OMP List_all_comb_b3_coef, List_all_comb_b3_expo, &
|
||||
!$OMP List_all_comb_b3_cent, int2_u_grad1u_j1b2_test)
|
||||
!$OMP DO
|
||||
do ipoint = 1, n_points_final_grid
|
||||
do i = 1, ao_num
|
||||
do j = i, ao_num
|
||||
if(dabs(ao_overlap_abs_grid(j,i)).lt.1.d-10)cycle
|
||||
dist_ij_ipoint = ao_prod_dist_grid(j,i,ipoint) ! distance to the grid point for the distribution |chi_i(r)chi_j(r)|
|
||||
sigma_ij = ao_prod_sigma(j,i) ! typical spatial extension of the distribution |chi_i(r)chi_j(r)|
|
||||
r(1) = final_grid_points(1,ipoint)
|
||||
r(2) = final_grid_points(2,ipoint)
|
||||
r(3) = final_grid_points(3,ipoint)
|
||||
|
||||
tmp = 0.d0
|
||||
do i_1s = 1, List_all_comb_b3_size
|
||||
|
||||
coef = List_all_comb_b3_coef (i_1s)
|
||||
beta = List_all_comb_b3_expo (i_1s)
|
||||
! if(beta.gt.1.d3)cycle
|
||||
if(dabs(coef).lt.1.d-10)cycle
|
||||
B_center(1) = List_all_comb_b3_cent(1,i_1s)
|
||||
B_center(2) = List_all_comb_b3_cent(2,i_1s)
|
||||
B_center(3) = List_all_comb_b3_cent(3,i_1s)
|
||||
dist = (B_center(1) - r(1)) * (B_center(1) - r(1)) &
|
||||
+ (B_center(2) - r(2)) * (B_center(2) - r(2)) &
|
||||
+ (B_center(3) - r(3)) * (B_center(3) - r(3))
|
||||
sigma_ij = 1.d0/sigma_ij
|
||||
sigma_ij *= sigma_ij
|
||||
sigma_ij *= 0.5d0
|
||||
double precision :: beta_ij, factor_ij_1s, center_ij_1s(3)
|
||||
! call gaussian_product(sigma_ij,ao_prod_center(1:3,j,i),beta,B_center,factor_ij_1s,beta_ij,center_ij_1s)
|
||||
! if(factor_ij_1s*ao_overlap_abs_grid(j,i).lt.1.d-15)cycle
|
||||
! if(factor_ij_1s*dsqpi_3_2*(beta_ij)**(-3/2)*ao_overlap_abs_grid(j,i).lt.1.d-20)cycle
|
||||
|
||||
do i_fit = 1, n_max_fit_slat
|
||||
|
||||
expo_fit = expo_gauss_j_mu_1_erf(i_fit)
|
||||
call gaussian_product(expo_fit,r,beta,B_center,factor_ij_1s,beta_ij,center_ij_1s)
|
||||
if(factor_ij_1s*ao_overlap_abs_grid(j,i).lt.1.d-15)cycle
|
||||
coef_fit = coef_gauss_j_mu_1_erf(i_fit)
|
||||
|
||||
alpha_1s = beta + expo_fit
|
||||
alpha_1s_inv = 1.d0 / alpha_1s
|
||||
centr_1s(1) = alpha_1s_inv * (beta * B_center(1) + expo_fit * r(1))
|
||||
centr_1s(2) = alpha_1s_inv * (beta * B_center(2) + expo_fit * r(2))
|
||||
centr_1s(3) = alpha_1s_inv * (beta * B_center(3) + expo_fit * r(3))
|
||||
|
||||
expo_coef_1s = beta * expo_fit * alpha_1s_inv * dist
|
||||
if(expo_coef_1s .gt. 20.d0) cycle
|
||||
coef_tmp = coef * coef_fit * dexp(-expo_coef_1s)
|
||||
if(dabs(coef_tmp) .lt. 1d-08) cycle
|
||||
|
||||
int_fit = NAI_pol_mult_erf_ao_with1s(i, j, alpha_1s, centr_1s, 1.d+9, r)
|
||||
|
||||
tmp += coef_tmp * int_fit
|
||||
enddo
|
||||
enddo
|
||||
|
||||
int2_u_grad1u_j1b2_test(j,i,ipoint) = tmp
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
|
||||
do ipoint = 1, n_points_final_grid
|
||||
do i = 2, ao_num
|
||||
do j = 1, i-1
|
||||
int2_u_grad1u_j1b2_test(j,i,ipoint) = int2_u_grad1u_j1b2_test(i,j,ipoint)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
call wall_time(wall1)
|
||||
print*, ' wall time for int2_u_grad1u_j1b2_test', wall1 - wall0
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
|
||||
BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2_test_2, (ao_num, ao_num, n_points_final_grid)]
|
||||
|
||||
BEGIN_DOC
|
||||
!
|
||||
! int dr2 phi_i(r2) phi_j(r2) 1s_j1b(r2)^2 u_12^mu [\grad_1 u_12^mu]
|
||||
!
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
integer :: i, j, ipoint, i_1s, i_fit
|
||||
double precision :: r(3), int_fit, expo_fit, coef_fit, coef_tmp
|
||||
double precision :: coef, beta, B_center(3), dist
|
||||
double precision :: alpha_1s, alpha_1s_inv, centr_1s(3), expo_coef_1s, tmp
|
||||
double precision :: wall0, wall1
|
||||
double precision, external :: NAI_pol_mult_erf_ao_with1s
|
||||
double precision :: j12_mu_r12,int_j1b
|
||||
double precision :: sigma_ij,dist_ij_ipoint,dsqpi_3_2
|
||||
double precision :: beta_ij,center_ij_1s(3),factor_ij_1s
|
||||
dsqpi_3_2 = (dacos(-1.d0))**(3/2)
|
||||
|
||||
provide mu_erf final_grid_points j1b_pen ao_overlap_abs List_comb_thr_b3_cent
|
||||
call wall_time(wall0)
|
||||
|
||||
|
||||
int2_u_grad1u_j1b2_test_2 = 0.d0
|
||||
|
||||
!$OMP PARALLEL DEFAULT (NONE) &
|
||||
!$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, &
|
||||
!$OMP coef_fit, expo_fit, int_fit, tmp, alpha_1s, dist, &
|
||||
!$OMP beta_ij,center_ij_1s,factor_ij_1s, &
|
||||
!$OMP int_j1b,alpha_1s_inv, centr_1s, expo_coef_1s, coef_tmp) &
|
||||
!$OMP SHARED (n_points_final_grid, ao_num, List_comb_b3_size_thr, &
|
||||
!$OMP final_grid_points, n_max_fit_slat, &
|
||||
!$OMP expo_gauss_j_mu_1_erf, coef_gauss_j_mu_1_erf, &
|
||||
!$OMP ao_prod_dist_grid, ao_prod_sigma, ao_overlap_abs_grid,ao_prod_center,dsqpi_3_2, &
|
||||
!$OMP List_comb_thr_b3_coef, List_comb_thr_b3_expo, ao_abs_comb_b3_j1b, &
|
||||
!$OMP List_comb_thr_b3_cent, int2_u_grad1u_j1b2_test_2)
|
||||
!$OMP DO
|
||||
do ipoint = 1, n_points_final_grid
|
||||
do i = 1, ao_num
|
||||
do j = i, ao_num
|
||||
if(dabs(ao_overlap_abs_grid(j,i)).lt.1.d-10)cycle
|
||||
r(1) = final_grid_points(1,ipoint)
|
||||
r(2) = final_grid_points(2,ipoint)
|
||||
r(3) = final_grid_points(3,ipoint)
|
||||
|
||||
tmp = 0.d0
|
||||
do i_1s = 1, List_comb_b3_size_thr(j,i)
|
||||
|
||||
coef = List_comb_thr_b3_coef (i_1s,j,i)
|
||||
beta = List_comb_thr_b3_expo (i_1s,j,i)
|
||||
int_j1b = ao_abs_comb_b3_j1b(i_1s,j,i)
|
||||
if(dabs(coef)*dabs(int_j1b).lt.1.d-10)cycle
|
||||
B_center(1) = List_comb_thr_b3_cent(1,i_1s,j,i)
|
||||
B_center(2) = List_comb_thr_b3_cent(2,i_1s,j,i)
|
||||
B_center(3) = List_comb_thr_b3_cent(3,i_1s,j,i)
|
||||
dist = (B_center(1) - r(1)) * (B_center(1) - r(1)) &
|
||||
+ (B_center(2) - r(2)) * (B_center(2) - r(2)) &
|
||||
+ (B_center(3) - r(3)) * (B_center(3) - r(3))
|
||||
|
||||
do i_fit = 1, n_max_fit_slat
|
||||
|
||||
expo_fit = expo_gauss_j_mu_1_erf(i_fit)
|
||||
call gaussian_product(expo_fit,r,beta,B_center,factor_ij_1s,beta_ij,center_ij_1s)
|
||||
if(factor_ij_1s*dabs(coef*int_j1b)*dsqpi_3_2*beta_ij**(-3/2).lt.1.d-15)cycle
|
||||
coef_fit = coef_gauss_j_mu_1_erf(i_fit)
|
||||
|
||||
alpha_1s = beta + expo_fit
|
||||
alpha_1s_inv = 1.d0 / alpha_1s
|
||||
centr_1s(1) = alpha_1s_inv * (beta * B_center(1) + expo_fit * r(1))
|
||||
centr_1s(2) = alpha_1s_inv * (beta * B_center(2) + expo_fit * r(2))
|
||||
centr_1s(3) = alpha_1s_inv * (beta * B_center(3) + expo_fit * r(3))
|
||||
|
||||
expo_coef_1s = beta * expo_fit * alpha_1s_inv * dist
|
||||
if(expo_coef_1s .gt. 20.d0) cycle
|
||||
coef_tmp = coef * coef_fit * dexp(-expo_coef_1s)
|
||||
if(dabs(coef_tmp) .lt. 1d-08) cycle
|
||||
|
||||
int_fit = NAI_pol_mult_erf_ao_with1s(i, j, alpha_1s, centr_1s, 1.d+9, r)
|
||||
|
||||
tmp += coef_tmp * int_fit
|
||||
enddo
|
||||
enddo
|
||||
|
||||
int2_u_grad1u_j1b2_test_2(j,i,ipoint) = tmp
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
|
||||
do ipoint = 1, n_points_final_grid
|
||||
do i = 2, ao_num
|
||||
do j = 1, i-1
|
||||
int2_u_grad1u_j1b2_test_2(j,i,ipoint) = int2_u_grad1u_j1b2_test_2(i,j,ipoint)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
call wall_time(wall1)
|
||||
print*, ' wall time for int2_u_grad1u_j1b2_test_2', wall1 - wall0
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
287
src/ao_many_one_e_ints/grad_lapl_jmu_manu.irp.f
Normal file
287
src/ao_many_one_e_ints/grad_lapl_jmu_manu.irp.f
Normal file
@ -0,0 +1,287 @@
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_j1b_test, (ao_num, ao_num, n_points_final_grid)]
|
||||
|
||||
BEGIN_DOC
|
||||
!
|
||||
! int dr phi_i(r) phi_j(r) 1s_j1b(r) (erf(mu(R) |r - R| - 1) / |r - R|
|
||||
!
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
integer :: i, j, ipoint, i_1s
|
||||
double precision :: r(3), int_mu, int_coulomb
|
||||
double precision :: coef, beta, B_center(3)
|
||||
double precision :: tmp,int_j1b
|
||||
double precision :: wall0, wall1
|
||||
double precision, external :: NAI_pol_mult_erf_ao_with1s
|
||||
double precision :: sigma_ij,dist_ij_ipoint,dsqpi_3_2
|
||||
dsqpi_3_2 = (dacos(-1.d0))**(3/2)
|
||||
provide mu_erf final_grid_points j1b_pen
|
||||
call wall_time(wall0)
|
||||
|
||||
v_ij_erf_rk_cst_mu_j1b_test = 0.d0
|
||||
|
||||
!$OMP PARALLEL DEFAULT (NONE) &
|
||||
!$OMP PRIVATE (ipoint, i, j, i_1s, r, coef, beta, B_center, int_mu, int_coulomb, tmp, int_j1b)&
|
||||
!$OMP SHARED (n_points_final_grid, ao_num, List_comb_b2_size_thr, final_grid_points, &
|
||||
!$OMP List_comb_thr_b2_coef, List_comb_thr_b2_expo, List_comb_thr_b2_cent,ao_abs_comb_b2_j1b, &
|
||||
!$OMP v_ij_erf_rk_cst_mu_j1b_test, mu_erf, &
|
||||
!$OMP ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,dsqpi_3_2)
|
||||
!$OMP DO
|
||||
!do ipoint = 1, 10
|
||||
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)
|
||||
|
||||
do i = 1, ao_num
|
||||
do j = i, ao_num
|
||||
if(dabs(ao_overlap_abs_grid(j,i)).lt.1.d-20)cycle
|
||||
|
||||
tmp = 0.d0
|
||||
do i_1s = 1, List_comb_b2_size_thr(j,i)
|
||||
|
||||
coef = List_comb_thr_b2_coef (i_1s,j,i)
|
||||
beta = List_comb_thr_b2_expo (i_1s,j,i)
|
||||
int_j1b = ao_abs_comb_b2_j1b(i_1s,j,i)
|
||||
if(dabs(coef)*dabs(int_j1b).lt.1.d-10)cycle
|
||||
B_center(1) = List_comb_thr_b2_cent(1,i_1s,j,i)
|
||||
B_center(2) = List_comb_thr_b2_cent(2,i_1s,j,i)
|
||||
B_center(3) = List_comb_thr_b2_cent(3,i_1s,j,i)
|
||||
|
||||
int_mu = NAI_pol_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r)
|
||||
int_coulomb = NAI_pol_mult_erf_ao_with1s(i, j, beta, B_center, 1.d+9, r)
|
||||
|
||||
tmp += coef * (int_mu - int_coulomb)
|
||||
enddo
|
||||
|
||||
v_ij_erf_rk_cst_mu_j1b_test(j,i,ipoint) = tmp
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
|
||||
do ipoint = 1, n_points_final_grid
|
||||
do i = 2, ao_num
|
||||
do j = 1, i-1
|
||||
v_ij_erf_rk_cst_mu_j1b_test(j,i,ipoint) = v_ij_erf_rk_cst_mu_j1b_test(i,j,ipoint)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
call wall_time(wall1)
|
||||
print*, ' wall time for v_ij_erf_rk_cst_mu_j1b_test', wall1 - wall0
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_j1b_test, (ao_num, ao_num, n_points_final_grid, 3)]
|
||||
|
||||
BEGIN_DOC
|
||||
! int dr x phi_i(r) phi_j(r) 1s_j1b(r) (erf(mu(R) |r - R|) - 1)/|r - R|
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
integer :: i, j, ipoint
|
||||
double precision :: wall0, wall1
|
||||
|
||||
call wall_time(wall0)
|
||||
|
||||
do ipoint = 1, n_points_final_grid
|
||||
do i = 1, ao_num
|
||||
do j = 1, ao_num
|
||||
x_v_ij_erf_rk_cst_mu_j1b_test(j,i,ipoint,1) = x_v_ij_erf_rk_cst_mu_tmp_j1b(1,j,i,ipoint)
|
||||
x_v_ij_erf_rk_cst_mu_j1b_test(j,i,ipoint,2) = x_v_ij_erf_rk_cst_mu_tmp_j1b(2,j,i,ipoint)
|
||||
x_v_ij_erf_rk_cst_mu_j1b_test(j,i,ipoint,3) = x_v_ij_erf_rk_cst_mu_tmp_j1b(3,j,i,ipoint)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
call wall_time(wall1)
|
||||
print*, ' wall time for x_v_ij_erf_rk_cst_mu_j1b_test', wall1 - wall0
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_tmp_j1b_test, (3, ao_num, ao_num, n_points_final_grid)]
|
||||
|
||||
BEGIN_DOC
|
||||
! int dr x phi_i(r) phi_j(r) 1s_j1b(r) (erf(mu(R) |r - R|) - 1)/|r - R|
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
integer :: i, j, ipoint, i_1s
|
||||
double precision :: coef, beta, B_center(3), r(3), ints(3), ints_coulomb(3)
|
||||
double precision :: tmp_x, tmp_y, tmp_z
|
||||
double precision :: wall0, wall1
|
||||
double precision :: sigma_ij,dist_ij_ipoint,dsqpi_3_2,int_j1b
|
||||
dsqpi_3_2 = (dacos(-1.d0))**(3/2)
|
||||
|
||||
call wall_time(wall0)
|
||||
|
||||
x_v_ij_erf_rk_cst_mu_tmp_j1b_test = 0.d0
|
||||
|
||||
!$OMP PARALLEL DEFAULT (NONE) &
|
||||
!$OMP PRIVATE (ipoint, i, j, i_1s, r, coef, beta, B_center, ints, ints_coulomb, &
|
||||
!$OMP int_j1b, tmp_x, tmp_y, tmp_z) &
|
||||
!$OMP SHARED (n_points_final_grid, ao_num, List_comb_b2_size_thr, final_grid_points,&
|
||||
!$OMP List_comb_thr_b2_coef, List_comb_thr_b2_expo, List_comb_thr_b2_cent, &
|
||||
!$OMP x_v_ij_erf_rk_cst_mu_tmp_j1b_test, mu_erf,ao_abs_comb_b2_j1b, &
|
||||
!$OMP ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,dsqpi_3_2)
|
||||
!$OMP DO
|
||||
!do ipoint = 1, 10
|
||||
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)
|
||||
|
||||
do i = 1, ao_num
|
||||
do j = i, ao_num
|
||||
if(dabs(ao_overlap_abs_grid(j,i)).lt.1.d-20)cycle
|
||||
|
||||
tmp_x = 0.d0
|
||||
tmp_y = 0.d0
|
||||
tmp_z = 0.d0
|
||||
do i_1s = 1, List_comb_b2_size_thr(j,i)
|
||||
|
||||
coef = List_comb_thr_b2_coef (i_1s,j,i)
|
||||
beta = List_comb_thr_b2_expo (i_1s,j,i)
|
||||
int_j1b = ao_abs_comb_b2_j1b(i_1s,j,i)
|
||||
if(dabs(coef)*dabs(int_j1b).lt.1.d-10)cycle
|
||||
B_center(1) = List_comb_thr_b2_cent(1,i_1s,j,i)
|
||||
B_center(2) = List_comb_thr_b2_cent(2,i_1s,j,i)
|
||||
B_center(3) = List_comb_thr_b2_cent(3,i_1s,j,i)
|
||||
|
||||
call NAI_pol_x_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r, ints )
|
||||
call NAI_pol_x_mult_erf_ao_with1s(i, j, beta, B_center, 1.d+9, r, ints_coulomb)
|
||||
|
||||
tmp_x += coef * (ints(1) - ints_coulomb(1))
|
||||
tmp_y += coef * (ints(2) - ints_coulomb(2))
|
||||
tmp_z += coef * (ints(3) - ints_coulomb(3))
|
||||
enddo
|
||||
|
||||
x_v_ij_erf_rk_cst_mu_tmp_j1b_test(1,j,i,ipoint) = tmp_x
|
||||
x_v_ij_erf_rk_cst_mu_tmp_j1b_test(2,j,i,ipoint) = tmp_y
|
||||
x_v_ij_erf_rk_cst_mu_tmp_j1b_test(3,j,i,ipoint) = tmp_z
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
|
||||
do ipoint = 1, n_points_final_grid
|
||||
do i = 2, ao_num
|
||||
do j = 1, i-1
|
||||
x_v_ij_erf_rk_cst_mu_tmp_j1b_test(1,j,i,ipoint) = x_v_ij_erf_rk_cst_mu_tmp_j1b_test(1,i,j,ipoint)
|
||||
x_v_ij_erf_rk_cst_mu_tmp_j1b_test(2,j,i,ipoint) = x_v_ij_erf_rk_cst_mu_tmp_j1b_test(2,i,j,ipoint)
|
||||
x_v_ij_erf_rk_cst_mu_tmp_j1b_test(3,j,i,ipoint) = x_v_ij_erf_rk_cst_mu_tmp_j1b_test(3,i,j,ipoint)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
call wall_time(wall1)
|
||||
print*, ' wall time for x_v_ij_erf_rk_cst_mu_tmp_j1b_test', wall1 - wall0
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
! TODO analytically
|
||||
BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_test, (ao_num, ao_num, n_points_final_grid)]
|
||||
|
||||
BEGIN_DOC
|
||||
!
|
||||
! int dr2 phi_i(r2) phi_j(r2) 1s_j1b(r2) u(mu, r12)
|
||||
!
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
integer :: i, j, ipoint, i_1s, i_fit
|
||||
double precision :: r(3), int_fit, expo_fit, coef_fit
|
||||
double precision :: coef, beta, B_center(3)
|
||||
double precision :: tmp
|
||||
double precision :: wall0, wall1
|
||||
|
||||
double precision, external :: overlap_gauss_r12_ao_with1s
|
||||
double precision :: sigma_ij,dist_ij_ipoint,dsqpi_3_2,int_j1b
|
||||
dsqpi_3_2 = (dacos(-1.d0))**(3/2)
|
||||
|
||||
provide mu_erf final_grid_points j1b_pen
|
||||
call wall_time(wall0)
|
||||
|
||||
v_ij_u_cst_mu_j1b_test = 0.d0
|
||||
|
||||
!$OMP PARALLEL DEFAULT (NONE) &
|
||||
!$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, &
|
||||
!$OMP beta_ij_u, factor_ij_1s_u, center_ij_1s_u, &
|
||||
!$OMP coef_fit, expo_fit, int_fit, tmp,coeftot,int_j1b) &
|
||||
!$OMP SHARED (n_points_final_grid, ao_num, &
|
||||
!$OMP final_grid_points, n_max_fit_slat, &
|
||||
!$OMP expo_gauss_j_mu_x, coef_gauss_j_mu_x, &
|
||||
!$OMP List_comb_thr_b2_coef, List_comb_thr_b2_expo,List_comb_b2_size_thr, &
|
||||
!$OMP List_comb_thr_b2_cent, v_ij_u_cst_mu_j1b_test,ao_abs_comb_b2_j1b, &
|
||||
!$OMP ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,dsqpi_3_2)
|
||||
!$OMP DO
|
||||
!do ipoint = 1, 10
|
||||
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)
|
||||
|
||||
do i = 1, ao_num
|
||||
do j = i, ao_num
|
||||
if(dabs(ao_overlap_abs_grid(j,i)).lt.1.d-20)cycle
|
||||
|
||||
tmp = 0.d0
|
||||
do i_1s = 1, List_comb_b2_size_thr(j,i)
|
||||
|
||||
coef = List_comb_thr_b2_coef (i_1s,j,i)
|
||||
beta = List_comb_thr_b2_expo (i_1s,j,i)
|
||||
int_j1b = ao_abs_comb_b2_j1b(i_1s,j,i)
|
||||
if(dabs(coef)*dabs(int_j1b).lt.1.d-10)cycle
|
||||
B_center(1) = List_comb_thr_b2_cent(1,i_1s,j,i)
|
||||
B_center(2) = List_comb_thr_b2_cent(2,i_1s,j,i)
|
||||
B_center(3) = List_comb_thr_b2_cent(3,i_1s,j,i)
|
||||
|
||||
do i_fit = 1, n_max_fit_slat
|
||||
|
||||
expo_fit = expo_gauss_j_mu_x(i_fit)
|
||||
coef_fit = coef_gauss_j_mu_x(i_fit)
|
||||
coeftot = coef * coef_fit
|
||||
if(dabs(coeftot).lt.1.d-15)cycle
|
||||
double precision :: beta_ij_u, factor_ij_1s_u, center_ij_1s_u(3),coeftot
|
||||
call gaussian_product(beta,B_center,expo_fit,r,factor_ij_1s_u,beta_ij_u,center_ij_1s_u)
|
||||
if(factor_ij_1s_u*ao_overlap_abs_grid(j,i).lt.1.d-15)cycle
|
||||
int_fit = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j)
|
||||
|
||||
tmp += coef * coef_fit * int_fit
|
||||
enddo
|
||||
enddo
|
||||
|
||||
v_ij_u_cst_mu_j1b_test(j,i,ipoint) = tmp
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
|
||||
do ipoint = 1, n_points_final_grid
|
||||
do i = 2, ao_num
|
||||
do j = 1, i-1
|
||||
v_ij_u_cst_mu_j1b_test(j,i,ipoint) = v_ij_u_cst_mu_j1b_test(i,j,ipoint)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
call wall_time(wall1)
|
||||
print*, ' wall time for v_ij_u_cst_mu_j1b_test', wall1 - wall0
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
188
src/ao_many_one_e_ints/listj1b_sorted.irp.f
Normal file
188
src/ao_many_one_e_ints/listj1b_sorted.irp.f
Normal file
@ -0,0 +1,188 @@
|
||||
|
||||
BEGIN_PROVIDER [ integer, List_comb_b2_size_thr, (ao_num, ao_num)]
|
||||
&BEGIN_PROVIDER [ integer, max_List_comb_b2_size_thr]
|
||||
implicit none
|
||||
integer :: i_1s,i,j,ipoint
|
||||
double precision :: coef,beta,center(3),int_j1b,thr
|
||||
double precision :: r(3),weight,dist
|
||||
thr = 1.d-10
|
||||
List_comb_b2_size_thr = 0
|
||||
do i = 1, ao_num
|
||||
do j = i, ao_num
|
||||
do i_1s = 1, List_all_comb_b2_size
|
||||
coef = List_all_comb_b2_coef (i_1s)
|
||||
if(dabs(coef).lt.1.d-10)cycle
|
||||
beta = List_all_comb_b2_expo (i_1s)
|
||||
center(1:3) = List_all_comb_b2_cent(1:3,i_1s)
|
||||
int_j1b = 0.d0
|
||||
do ipoint = 1, n_points_final_grid
|
||||
r(1:3) = final_grid_points(1:3,ipoint)
|
||||
weight = final_weight_at_r_vector(ipoint)
|
||||
dist = ( center(1) - r(1) )*( center(1) - r(1) )
|
||||
dist += ( center(2) - r(2) )*( center(2) - r(2) )
|
||||
dist += ( center(3) - r(3) )*( center(3) - r(3) )
|
||||
int_j1b += dabs(aos_in_r_array_transp(ipoint,i) * aos_in_r_array_transp(ipoint,j))*dexp(-beta*dist) * weight
|
||||
enddo
|
||||
if(dabs(coef)*dabs(int_j1b).gt.thr)then
|
||||
List_comb_b2_size_thr(j,i) += 1
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
do i = 1, ao_num
|
||||
do j = 1, i-1
|
||||
List_comb_b2_size_thr(j,i) = List_comb_b2_size_thr(i,j)
|
||||
enddo
|
||||
enddo
|
||||
integer :: list(ao_num)
|
||||
do i = 1, ao_num
|
||||
list(i) = maxval(List_comb_b2_size_thr(:,i))
|
||||
enddo
|
||||
max_List_comb_b2_size_thr = maxval(list)
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ double precision, List_comb_thr_b2_coef, ( max_List_comb_b2_size_thr,ao_num, ao_num )]
|
||||
&BEGIN_PROVIDER [ double precision, List_comb_thr_b2_expo, ( max_List_comb_b2_size_thr,ao_num, ao_num )]
|
||||
&BEGIN_PROVIDER [ double precision, List_comb_thr_b2_cent, (3, max_List_comb_b2_size_thr,ao_num, ao_num )]
|
||||
&BEGIN_PROVIDER [ double precision, ao_abs_comb_b2_j1b, ( max_List_comb_b2_size_thr ,ao_num, ao_num)]
|
||||
implicit none
|
||||
integer :: i_1s,i,j,ipoint,icount
|
||||
double precision :: coef,beta,center(3),int_j1b,thr
|
||||
double precision :: r(3),weight,dist
|
||||
thr = 1.d-10
|
||||
ao_abs_comb_b2_j1b = 10000000.d0
|
||||
do i = 1, ao_num
|
||||
do j = i, ao_num
|
||||
icount = 0
|
||||
do i_1s = 1, List_all_comb_b2_size
|
||||
coef = List_all_comb_b2_coef (i_1s)
|
||||
if(dabs(coef).lt.1.d-10)cycle
|
||||
beta = List_all_comb_b2_expo (i_1s)
|
||||
center(1:3) = List_all_comb_b2_cent(1:3,i_1s)
|
||||
int_j1b = 0.d0
|
||||
do ipoint = 1, n_points_final_grid
|
||||
r(1:3) = final_grid_points(1:3,ipoint)
|
||||
weight = final_weight_at_r_vector(ipoint)
|
||||
dist = ( center(1) - r(1) )*( center(1) - r(1) )
|
||||
dist += ( center(2) - r(2) )*( center(2) - r(2) )
|
||||
dist += ( center(3) - r(3) )*( center(3) - r(3) )
|
||||
int_j1b += dabs(aos_in_r_array_transp(ipoint,i) * aos_in_r_array_transp(ipoint,j))*dexp(-beta*dist) * weight
|
||||
enddo
|
||||
if(dabs(coef)*dabs(int_j1b).gt.thr)then
|
||||
icount += 1
|
||||
List_comb_thr_b2_coef(icount,j,i) = coef
|
||||
List_comb_thr_b2_expo(icount,j,i) = beta
|
||||
List_comb_thr_b2_cent(1:3,icount,j,i) = center(1:3)
|
||||
ao_abs_comb_b2_j1b(icount,j,i) = int_j1b
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
do i = 1, ao_num
|
||||
do j = 1, i-1
|
||||
do icount = 1, List_comb_b2_size_thr(j,i)
|
||||
List_comb_thr_b2_coef(icount,j,i) = List_comb_thr_b2_coef(icount,i,j)
|
||||
List_comb_thr_b2_expo(icount,j,i) = List_comb_thr_b2_expo(icount,i,j)
|
||||
List_comb_thr_b2_cent(1:3,icount,j,i) = List_comb_thr_b2_cent(1:3,icount,i,j)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ integer, List_comb_b3_size_thr, (ao_num, ao_num)]
|
||||
&BEGIN_PROVIDER [ integer, max_List_comb_b3_size_thr]
|
||||
implicit none
|
||||
integer :: i_1s,i,j,ipoint
|
||||
double precision :: coef,beta,center(3),int_j1b,thr
|
||||
double precision :: r(3),weight,dist
|
||||
thr = 1.d-14
|
||||
List_comb_b3_size_thr = 0
|
||||
do i = 1, ao_num
|
||||
do j = i, ao_num
|
||||
do i_1s = 1, List_all_comb_b3_size
|
||||
coef = List_all_comb_b3_coef (i_1s)
|
||||
if(dabs(coef).lt.thr)cycle
|
||||
beta = List_all_comb_b3_expo (i_1s)
|
||||
center(1:3) = List_all_comb_b3_cent(1:3,i_1s)
|
||||
int_j1b = 0.d0
|
||||
do ipoint = 1, n_points_final_grid
|
||||
r(1:3) = final_grid_points(1:3,ipoint)
|
||||
weight = final_weight_at_r_vector(ipoint)
|
||||
dist = ( center(1) - r(1) )*( center(1) - r(1) )
|
||||
dist += ( center(2) - r(2) )*( center(2) - r(2) )
|
||||
dist += ( center(3) - r(3) )*( center(3) - r(3) )
|
||||
int_j1b += dabs(aos_in_r_array_transp(ipoint,i) * aos_in_r_array_transp(ipoint,j))*dexp(-beta*dist) * weight
|
||||
enddo
|
||||
if(dabs(coef)*dabs(int_j1b).gt.thr)then
|
||||
List_comb_b3_size_thr(j,i) += 1
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
do i = 1, ao_num
|
||||
do j = 1, i-1
|
||||
List_comb_b3_size_thr(j,i) = List_comb_b3_size_thr(i,j)
|
||||
enddo
|
||||
enddo
|
||||
integer :: list(ao_num)
|
||||
do i = 1, ao_num
|
||||
list(i) = maxval(List_comb_b3_size_thr(:,i))
|
||||
enddo
|
||||
max_List_comb_b3_size_thr = maxval(list)
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ double precision, List_comb_thr_b3_coef, ( max_List_comb_b3_size_thr,ao_num, ao_num )]
|
||||
&BEGIN_PROVIDER [ double precision, List_comb_thr_b3_expo, ( max_List_comb_b3_size_thr,ao_num, ao_num )]
|
||||
&BEGIN_PROVIDER [ double precision, List_comb_thr_b3_cent, (3, max_List_comb_b3_size_thr,ao_num, ao_num )]
|
||||
&BEGIN_PROVIDER [ double precision, ao_abs_comb_b3_j1b, ( max_List_comb_b3_size_thr ,ao_num, ao_num)]
|
||||
implicit none
|
||||
integer :: i_1s,i,j,ipoint,icount
|
||||
double precision :: coef,beta,center(3),int_j1b,thr
|
||||
double precision :: r(3),weight,dist
|
||||
thr = 1.d-14
|
||||
ao_abs_comb_b3_j1b = 10000000.d0
|
||||
do i = 1, ao_num
|
||||
do j = i, ao_num
|
||||
icount = 0
|
||||
do i_1s = 1, List_all_comb_b3_size
|
||||
coef = List_all_comb_b3_coef (i_1s)
|
||||
if(dabs(coef).lt.thr)cycle
|
||||
beta = List_all_comb_b3_expo (i_1s)
|
||||
center(1:3) = List_all_comb_b3_cent(1:3,i_1s)
|
||||
int_j1b = 0.d0
|
||||
do ipoint = 1, n_points_final_grid
|
||||
r(1:3) = final_grid_points(1:3,ipoint)
|
||||
weight = final_weight_at_r_vector(ipoint)
|
||||
dist = ( center(1) - r(1) )*( center(1) - r(1) )
|
||||
dist += ( center(2) - r(2) )*( center(2) - r(2) )
|
||||
dist += ( center(3) - r(3) )*( center(3) - r(3) )
|
||||
int_j1b += dabs(aos_in_r_array_transp(ipoint,i) * aos_in_r_array_transp(ipoint,j))*dexp(-beta*dist) * weight
|
||||
enddo
|
||||
if(dabs(coef)*dabs(int_j1b).gt.thr)then
|
||||
icount += 1
|
||||
List_comb_thr_b3_coef(icount,j,i) = coef
|
||||
List_comb_thr_b3_expo(icount,j,i) = beta
|
||||
List_comb_thr_b3_cent(1:3,icount,j,i) = center(1:3)
|
||||
ao_abs_comb_b3_j1b(icount,j,i) = int_j1b
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
do i = 1, ao_num
|
||||
do j = 1, i-1
|
||||
do icount = 1, List_comb_b3_size_thr(j,i)
|
||||
List_comb_thr_b3_coef(icount,j,i) = List_comb_thr_b3_coef(icount,i,j)
|
||||
List_comb_thr_b3_expo(icount,j,i) = List_comb_thr_b3_expo(icount,i,j)
|
||||
List_comb_thr_b3_cent(1,icount,j,i) = List_comb_thr_b3_cent(1,icount,i,j)
|
||||
List_comb_thr_b3_cent(2,icount,j,i) = List_comb_thr_b3_cent(2,icount,i,j)
|
||||
List_comb_thr_b3_cent(3,icount,j,i) = List_comb_thr_b3_cent(3,icount,i,j)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
|
||||
END_PROVIDER
|
@ -109,3 +109,8 @@ BEGIN_PROVIDER [ double precision, ao_prod_dist_grid, (ao_num, ao_num, n_points_
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
!BEGIN_PROVIDER [ double precision, ao_abs_prod_j1b, (ao_num, ao_num)]
|
||||
! implicit none
|
||||
!
|
||||
!END_PROVIDER
|
||||
|
@ -14,11 +14,40 @@ program test_ints
|
||||
my_n_pt_r_grid = 10 ! small grid for quick debug
|
||||
my_n_pt_a_grid = 26 ! small grid for quick debug
|
||||
touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
|
||||
call routine
|
||||
! call routine_int2_u_grad1u_j1b2
|
||||
! call routine_v_ij_erf_rk_cst_mu_j1b
|
||||
! call routine_x_v_ij_erf_rk_cst_mu_tmp_j1b
|
||||
call routine_v_ij_u_cst_mu_j1b
|
||||
! call routine_test_j1b
|
||||
|
||||
end
|
||||
|
||||
subroutine routine
|
||||
subroutine routine_test_j1b
|
||||
implicit none
|
||||
integer :: i,icount,j
|
||||
icount = 0
|
||||
! do i = 1, List_all_comb_b2_size
|
||||
! if(dabs(List_all_comb_b2_coef(i)).gt.1.d-10)then
|
||||
! icount += 1
|
||||
! endif
|
||||
! print*,i,List_all_comb_b2_expo(i),List_all_comb_b2_coef(i)
|
||||
! enddo
|
||||
! print*,'List_all_comb_b2_coef,icount = ',List_all_comb_b2_size
|
||||
do i = 1, ao_num
|
||||
do j = 1, ao_num
|
||||
do icount = 1, List_comb_b3_size_thr(j,i)
|
||||
print*,List_comb_thr_b3_cent(1:3,icount,j,i)
|
||||
! print*,'',j,i
|
||||
! print*,List_comb_b2_size_thr(j,i),List_comb_b3_size_thr(j,i),ao_overlap_abs_grid(j,i)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
print*,'max_List_comb_b2_size_thr = ',max_List_comb_b2_size_thr,List_all_comb_b2_size
|
||||
print*,'max_List_comb_b2_size_thr = ',max_List_comb_b3_size_thr,List_all_comb_b3_size
|
||||
|
||||
end
|
||||
|
||||
subroutine routine_int2_u_grad1u_j1b2
|
||||
implicit none
|
||||
integer :: i,j,ipoint,k,l
|
||||
double precision :: weight,accu_relat, accu_abs, contrib
|
||||
@ -47,7 +76,15 @@ subroutine routine
|
||||
do l = 1, ao_num
|
||||
do i = 1, ao_num
|
||||
do j = 1, ao_num
|
||||
array(j,i,l,k) += int2_u_grad1u_j1b2_test(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight
|
||||
array(j,i,l,k) += int2_u_grad1u_j1b2_test_2(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight
|
||||
! if(dabs(int2_u_grad1u_j1b2(j,i,ipoint)).gt.1.d-6)then
|
||||
! if(dabs(int2_u_grad1u_j1b2_test_2(j,i,ipoint)-int2_u_grad1u_j1b2(j,i,ipoint)).gt.1.d-6)then
|
||||
! print*,int2_u_grad1u_j1b2(j,i,ipoint), int2_u_grad1u_j1b2_test_2(j,i,ipoint),dabs(int2_u_grad1u_j1b2_test_2(j,i,ipoint)-int2_u_grad1u_j1b2(j,i,ipoint))
|
||||
! print*,i,j
|
||||
! print*,final_grid_points(:,i)
|
||||
! stop
|
||||
! endif
|
||||
! endif
|
||||
array_ref(j,i,l,k) += int2_u_grad1u_j1b2(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight
|
||||
enddo
|
||||
enddo
|
||||
@ -74,4 +111,183 @@ subroutine routine
|
||||
|
||||
|
||||
|
||||
end
|
||||
|
||||
subroutine routine_v_ij_erf_rk_cst_mu_j1b
|
||||
implicit none
|
||||
integer :: i,j,ipoint,k,l
|
||||
double precision :: weight,accu_relat, accu_abs, contrib
|
||||
double precision, allocatable :: array(:,:,:,:), array_ref(:,:,:,:)
|
||||
! print*,'ao_overlap_abs = '
|
||||
! do i = 1, ao_num
|
||||
! write(*,'(100(F10.5,X))')ao_overlap_abs(i,:)
|
||||
! enddo
|
||||
! print*,'center = '
|
||||
! do i = 1, ao_num
|
||||
! write(*,'(100(F10.5,X))')ao_prod_center(2,i,:)
|
||||
! enddo
|
||||
! print*,'sigma = '
|
||||
! do i = 1, ao_num
|
||||
! write(*,'(100(F10.5,X))')ao_prod_sigma(i,:)
|
||||
! enddo
|
||||
|
||||
|
||||
allocate(array(ao_num, ao_num, ao_num, ao_num))
|
||||
array = 0.d0
|
||||
allocate(array_ref(ao_num, ao_num, ao_num, ao_num))
|
||||
array_ref = 0.d0
|
||||
do ipoint = 1, n_points_final_grid
|
||||
weight = final_weight_at_r_vector(ipoint)
|
||||
do k = 1, ao_num
|
||||
do l = 1, ao_num
|
||||
do i = 1, ao_num
|
||||
do j = 1, ao_num
|
||||
array(j,i,l,k) += v_ij_erf_rk_cst_mu_j1b_test(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight
|
||||
array_ref(j,i,l,k) += v_ij_erf_rk_cst_mu_j1b(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
accu_relat = 0.d0
|
||||
accu_abs = 0.d0
|
||||
do k = 1, ao_num
|
||||
do l = 1, ao_num
|
||||
do i = 1, ao_num
|
||||
do j = 1, ao_num
|
||||
contrib = dabs(array(j,i,l,k) - array_ref(j,i,l,k))
|
||||
accu_abs += contrib
|
||||
if(dabs(array_ref(j,i,l,k)).gt.1.d-10)then
|
||||
accu_relat += contrib/dabs(array_ref(j,i,l,k))
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
print*,'accu_abs = ',accu_abs/dble(ao_num)**4
|
||||
print*,'accu_relat = ',accu_relat/dble(ao_num)**4
|
||||
|
||||
|
||||
|
||||
end
|
||||
|
||||
|
||||
subroutine routine_x_v_ij_erf_rk_cst_mu_tmp_j1b
|
||||
implicit none
|
||||
integer :: i,j,ipoint,k,l,m
|
||||
double precision :: weight,accu_relat, accu_abs, contrib
|
||||
double precision, allocatable :: array(:,:,:,:), array_ref(:,:,:,:)
|
||||
! print*,'ao_overlap_abs = '
|
||||
! do i = 1, ao_num
|
||||
! write(*,'(100(F10.5,X))')ao_overlap_abs(i,:)
|
||||
! enddo
|
||||
! print*,'center = '
|
||||
! do i = 1, ao_num
|
||||
! write(*,'(100(F10.5,X))')ao_prod_center(2,i,:)
|
||||
! enddo
|
||||
! print*,'sigma = '
|
||||
! do i = 1, ao_num
|
||||
! write(*,'(100(F10.5,X))')ao_prod_sigma(i,:)
|
||||
! enddo
|
||||
|
||||
|
||||
allocate(array(ao_num, ao_num, ao_num, ao_num))
|
||||
array = 0.d0
|
||||
allocate(array_ref(ao_num, ao_num, ao_num, ao_num))
|
||||
array_ref = 0.d0
|
||||
do ipoint = 1, n_points_final_grid
|
||||
weight = final_weight_at_r_vector(ipoint)
|
||||
do k = 1, ao_num
|
||||
do l = 1, ao_num
|
||||
do i = 1, ao_num
|
||||
do j = 1, ao_num
|
||||
do m = 1, 3
|
||||
array(j,i,l,k) += x_v_ij_erf_rk_cst_mu_tmp_j1b_test(m,j,i,ipoint) * aos_grad_in_r_array_transp(m,k,ipoint) * aos_in_r_array(l,ipoint) * weight
|
||||
array_ref(j,i,l,k) += x_v_ij_erf_rk_cst_mu_tmp_j1b(m,j,i,ipoint) * aos_grad_in_r_array_transp(m,k,ipoint) * aos_in_r_array(l,ipoint) * weight
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
accu_relat = 0.d0
|
||||
accu_abs = 0.d0
|
||||
do k = 1, ao_num
|
||||
do l = 1, ao_num
|
||||
do i = 1, ao_num
|
||||
do j = 1, ao_num
|
||||
contrib = dabs(array(j,i,l,k) - array_ref(j,i,l,k))
|
||||
accu_abs += contrib
|
||||
if(dabs(array_ref(j,i,l,k)).gt.1.d-10)then
|
||||
accu_relat += contrib/dabs(array_ref(j,i,l,k))
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
print*,'accu_abs = ',accu_abs/dble(ao_num)**4
|
||||
print*,'accu_relat = ',accu_relat/dble(ao_num)**4
|
||||
|
||||
|
||||
|
||||
end
|
||||
|
||||
|
||||
|
||||
subroutine routine_v_ij_u_cst_mu_j1b
|
||||
implicit none
|
||||
integer :: i,j,ipoint,k,l
|
||||
double precision :: weight,accu_relat, accu_abs, contrib
|
||||
double precision, allocatable :: array(:,:,:,:), array_ref(:,:,:,:)
|
||||
! print*,'ao_overlap_abs = '
|
||||
! do i = 1, ao_num
|
||||
! write(*,'(100(F10.5,X))')ao_overlap_abs(i,:)
|
||||
! enddo
|
||||
! print*,'center = '
|
||||
! do i = 1, ao_num
|
||||
! write(*,'(100(F10.5,X))')ao_prod_center(2,i,:)
|
||||
! enddo
|
||||
! print*,'sigma = '
|
||||
! do i = 1, ao_num
|
||||
! write(*,'(100(F10.5,X))')ao_prod_sigma(i,:)
|
||||
! enddo
|
||||
|
||||
|
||||
allocate(array(ao_num, ao_num, ao_num, ao_num))
|
||||
array = 0.d0
|
||||
allocate(array_ref(ao_num, ao_num, ao_num, ao_num))
|
||||
array_ref = 0.d0
|
||||
do ipoint = 1, n_points_final_grid
|
||||
weight = final_weight_at_r_vector(ipoint)
|
||||
do k = 1, ao_num
|
||||
do l = 1, ao_num
|
||||
do i = 1, ao_num
|
||||
do j = 1, ao_num
|
||||
array(j,i,l,k) += v_ij_u_cst_mu_j1b_test(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight
|
||||
array_ref(j,i,l,k) += v_ij_u_cst_mu_j1b(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
accu_relat = 0.d0
|
||||
accu_abs = 0.d0
|
||||
do k = 1, ao_num
|
||||
do l = 1, ao_num
|
||||
do i = 1, ao_num
|
||||
do j = 1, ao_num
|
||||
contrib = dabs(array(j,i,l,k) - array_ref(j,i,l,k))
|
||||
accu_abs += contrib
|
||||
if(dabs(array_ref(j,i,l,k)).gt.1.d-10)then
|
||||
accu_relat += contrib/dabs(array_ref(j,i,l,k))
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
print*,'accu_abs = ',accu_abs/dble(ao_num)**4
|
||||
print*,'accu_relat = ',accu_relat/dble(ao_num)**4
|
||||
|
||||
|
||||
|
||||
end
|
||||
|
Loading…
Reference in New Issue
Block a user