From 436b8815807c3ad04dc5035244c6ed9e7c4f2580 Mon Sep 17 00:00:00 2001 From: eginer Date: Tue, 28 Mar 2023 11:21:19 +0200 Subject: [PATCH] added thresh_cycle --- src/ao_many_one_e_ints/grad2_jmu_manu.irp.f | 48 ++-- .../grad_lapl_jmu_manu.irp.f | 46 ++-- src/ao_many_one_e_ints/listj1b_sorted.irp.f | 42 ++-- src/non_h_ints_mu/total_tc_int.irp.f | 19 ++ src/tc_keywords/EZFIO.cfg | 6 + src/tc_scf/tc_scf.irp.f | 4 +- src/tc_scf/test_int.irp.f | 214 +++++++++++------- src/tools/sort_wf.irp.f | 2 +- 8 files changed, 212 insertions(+), 169 deletions(-) diff --git a/src/ao_many_one_e_ints/grad2_jmu_manu.irp.f b/src/ao_many_one_e_ints/grad2_jmu_manu.irp.f index 8e253d75..14170ede 100644 --- a/src/ao_many_one_e_ints/grad2_jmu_manu.irp.f +++ b/src/ao_many_one_e_ints/grad2_jmu_manu.irp.f @@ -38,7 +38,7 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test, (ao_num, ao_n !$OMP expo_gauss_1_erf_x_2, coef_gauss_1_erf_x_2, & !$OMP List_comb_thr_b3_coef, List_comb_thr_b3_expo, & !$OMP List_comb_thr_b3_cent, int2_grad1u2_grad2u2_j1b2_test, ao_abs_comb_b3_j1b, & - !$OMP ao_overlap_abs,sq_pi_3_2) + !$OMP ao_overlap_abs,sq_pi_3_2,thrsh_cycle_tc) !$OMP DO SCHEDULE(dynamic) do ipoint = 1, n_points_final_grid r(1) = final_grid_points(1,ipoint) @@ -46,7 +46,7 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test, (ao_num, ao_n r(3) = final_grid_points(3,ipoint) do i = 1, ao_num do j = i, ao_num - if(ao_overlap_abs(j,i) .lt. 1.d-12) then + if(ao_overlap_abs(j,i) .lt. thrsh_cycle_tc) then cycle endif @@ -58,7 +58,7 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test, (ao_num, ao_n do i_fit = 1, ng_fit_jast expo_fit = expo_gauss_1_erf_x_2(i_fit) coef_fit = -0.25d0 * coef_gauss_1_erf_x_2(i_fit) - if(dabs(coef_fit*int_j1b*sq_pi_3_2*(expo_fit)**(-1.5d0)).lt.1.d-10)cycle +! if(dabs(coef_fit*int_j1b*sq_pi_3_2*(expo_fit)**(-1.5d0)).lt.thrsh_cycle_tc)cycle int_gauss = overlap_gauss_r12_ao(r, expo_fit, i, j) int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint) += coef_fit * int_gauss enddo @@ -81,8 +81,7 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test, (ao_num, ao_n !DIR$ FORCEINLINE call gaussian_product(expo_fit,r,beta,B_center,factor_ij_1s,beta_ij,center_ij_1s) coef_fit = -0.25d0 * coef_gauss_1_erf_x_2(i_fit) * coef -! if(dabs(coef_fit*factor_ij_1s*int_j1b).lt.1.d-10)cycle ! old version - if(dabs(coef_fit*factor_ij_1s*int_j1b*sq_pi_3_2*(beta_ij)**(-1.5d0)).lt.1.d-10)cycle +! if(dabs(coef_fit*factor_ij_1s*int_j1b*sq_pi_3_2*(beta_ij)**(-1.5d0)).lt.thrsh_cycle_tc)cycle ! call overlap_gauss_r12_ao_with1s_v(B_center, beta, final_grid_points_transp, & ! expo_fit, i, j, int_fit_v, n_points_final_grid) int_gauss = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j) @@ -145,14 +144,14 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test_v, (ao_num, ao !$OMP expo_gauss_1_erf_x_2, coef_gauss_1_erf_x_2, & !$OMP List_comb_thr_b3_coef, List_comb_thr_b3_expo, & !$OMP List_comb_thr_b3_cent, big_array,& - !$OMP ao_abs_comb_b3_j1b,ao_overlap_abs) + !$OMP ao_abs_comb_b3_j1b,ao_overlap_abs,thrsh_cycle_tc) ! allocate(int_fit_v(n_points_final_grid)) !$OMP DO SCHEDULE(dynamic) do i = 1, ao_num do j = i, ao_num - if(ao_overlap_abs(j,i) .lt. 1.d-12) then + if(ao_overlap_abs(j,i) .lt. thrsh_cycle_tc) then cycle endif @@ -161,7 +160,6 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test_v, (ao_num, ao 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-15)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) @@ -243,7 +241,7 @@ BEGIN_PROVIDER [ double precision, int2_u2_j1b2_test, (ao_num, ao_num, n_points_ !$OMP final_grid_points, ng_fit_jast, & !$OMP expo_gauss_j_mu_x_2, coef_gauss_j_mu_x_2, & !$OMP List_comb_thr_b3_coef, List_comb_thr_b3_expo,sq_pi_3_2, & - !$OMP List_comb_thr_b3_cent, int2_u2_j1b2_test,ao_abs_comb_b3_j1b) + !$OMP List_comb_thr_b3_cent, int2_u2_j1b2_test,ao_abs_comb_b3_j1b,thrsh_cycle_tc) !$OMP DO do ipoint = 1, n_points_final_grid r(1) = final_grid_points(1,ipoint) @@ -260,11 +258,11 @@ BEGIN_PROVIDER [ double precision, int2_u2_j1b2_test, (ao_num, ao_num, n_points_ ! --- --- --- int_j1b = ao_abs_comb_b3_j1b(1,j,i) - if(dabs(int_j1b).lt.1.d-10) cycle + if(dabs(int_j1b).lt.thrsh_cycle_tc) cycle do i_fit = 1, ng_fit_jast expo_fit = expo_gauss_j_mu_x_2(i_fit) coef_fit = coef_gauss_j_mu_x_2(i_fit) - if(dabs(coef_fit*int_j1b*sq_pi_3_2*(expo_fit)**(-1.5d0)).lt.1.d-10)cycle +! if(dabs(coef_fit*int_j1b*sq_pi_3_2*(expo_fit)**(-1.5d0)).lt.thrsh_cycle_tc)cycle int_fit = overlap_gauss_r12_ao(r, expo_fit, i, j) tmp += coef_fit * int_fit enddo @@ -278,7 +276,7 @@ BEGIN_PROVIDER [ double precision, int2_u2_j1b2_test, (ao_num, ao_num, n_points_ 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 +! if(dabs(coef)*dabs(int_j1b).lt.thrsh_cycle_tc)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) @@ -288,8 +286,7 @@ BEGIN_PROVIDER [ double precision, int2_u2_j1b2_test, (ao_num, ao_num, n_points_ coef_fit = coef_gauss_j_mu_x_2(i_fit) !DIR$ FORCEINLINE call gaussian_product(expo_fit,r,beta,B_center,factor_ij_1s,beta_ij,center_ij_1s) -! if(dabs(coef_fit*coef*factor_ij_1s*int_j1b).lt.1.d-10)cycle ! old version - if(dabs(coef_fit*coef*factor_ij_1s*int_j1b*sq_pi_3_2*(beta_ij)**(-1.5d0)).lt.1.d-10)cycle +! if(dabs(coef_fit*coef*factor_ij_1s*int_j1b*sq_pi_3_2*(beta_ij)**(-1.5d0)).lt.thrsh_cycle_tc)cycle int_fit = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j) tmp += coef * coef_fit * int_fit enddo @@ -350,7 +347,7 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_x_j1b2_test, (ao_num, ao_num, n !$OMP final_grid_points, ng_fit_jast, & !$OMP expo_gauss_j_mu_1_erf, coef_gauss_j_mu_1_erf, & !$OMP List_comb_thr_b3_coef, List_comb_thr_b3_expo, & - !$OMP List_comb_thr_b3_cent, int2_u_grad1u_x_j1b2_test,ao_abs_comb_b3_j1b,sq_pi_3_2) + !$OMP List_comb_thr_b3_cent, int2_u_grad1u_x_j1b2_test,ao_abs_comb_b3_j1b,sq_pi_3_2,thrsh_cycle_tc) !$OMP DO do ipoint = 1, n_points_final_grid @@ -369,7 +366,7 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_x_j1b2_test, (ao_num, ao_num, n 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 + if(dabs(coef)*dabs(int_j1b).lt.thrsh_cycle_tc)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) @@ -392,8 +389,7 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_x_j1b2_test, (ao_num, ao_num, n expo_coef_1s = beta * expo_fit * alpha_1s_inv * dist coef_tmp = coef * coef_fit * dexp(-expo_coef_1s) sq_alpha = alpha_1s_inv * dsqrt(alpha_1s_inv) -! if(dabs(coef_tmp*int_j1b) .lt. 1d-10) cycle ! old version - if(dabs(coef_tmp*int_j1b*sq_pi_3_2*sq_alpha) .lt. 1d-10) cycle +! if(dabs(coef_tmp*int_j1b*sq_pi_3_2*sq_alpha) .lt. thrsh_cycle_tc) cycle call NAI_pol_x_mult_erf_ao_with1s(i, j, alpha_1s, centr_1s, 1.d+9, r, int_fit) @@ -470,13 +466,13 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2_test, (ao_num, ao_num, n_p !$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) + !$OMP List_comb_thr_b3_cent, int2_u_grad1u_j1b2_test,thrsh_cycle_tc) !$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 + if(dabs(ao_overlap_abs_grid(j,i)).lt.thrsh_cycle_tc) cycle r(1) = final_grid_points(1,ipoint) r(2) = final_grid_points(2,ipoint) @@ -489,10 +485,10 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2_test, (ao_num, ao_num, n_p ! --- --- --- int_j1b = ao_abs_comb_b3_j1b(1,j,i) - if(dabs(int_j1b).lt.1.d-10) cycle +! if(dabs(int_j1b).lt.thrsh_cycle_tc) cycle do i_fit = 1, ng_fit_jast expo_fit = expo_gauss_j_mu_1_erf(i_fit) - if(dabs(int_j1b)*dsqpi_3_2*expo_fit**(-1.5d0).lt.1.d-15) cycle +! if(dabs(int_j1b)*dsqpi_3_2*expo_fit**(-1.5d0).lt.thrsh_cycle_tc) cycle coef_fit = coef_gauss_j_mu_1_erf(i_fit) int_fit = NAI_pol_mult_erf_ao_with1s(i, j, expo_fit, r, 1.d+9, r) tmp += coef_fit * int_fit @@ -507,7 +503,7 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2_test, (ao_num, ao_num, n_p 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 +! if(dabs(coef)*dabs(int_j1b).lt.thrsh_cycle_tc)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) @@ -517,7 +513,7 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2_test, (ao_num, ao_num, n_p do i_fit = 1, ng_fit_jast 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**(-1.5d0).lt.1.d-15)cycle +! if(factor_ij_1s*dabs(coef*int_j1b)*dsqpi_3_2*beta_ij**(-1.5d0).lt.thrsh_cycle_tc)cycle coef_fit = coef_gauss_j_mu_1_erf(i_fit) alpha_1s = beta + expo_fit @@ -527,9 +523,9 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2_test, (ao_num, ao_num, n_p 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 +! 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 +! 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) diff --git a/src/ao_many_one_e_ints/grad_lapl_jmu_manu.irp.f b/src/ao_many_one_e_ints/grad_lapl_jmu_manu.irp.f index 5c9f81e9..66a2b961 100644 --- a/src/ao_many_one_e_ints/grad_lapl_jmu_manu.irp.f +++ b/src/ao_many_one_e_ints/grad_lapl_jmu_manu.irp.f @@ -31,7 +31,7 @@ BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_j1b_test, (ao_num, ao_num, !$OMP SHARED (n_points_final_grid, ao_num, List_comb_thr_b2_size, 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 ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,dsqpi_3_2,thrsh_cycle_tc) !$OMP DO !do ipoint = 1, 10 do ipoint = 1, n_points_final_grid @@ -41,7 +41,7 @@ BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_j1b_test, (ao_num, ao_num, do i = 1, ao_num do j = i, ao_num - if(dabs(ao_overlap_abs_grid(j,i)).lt.1.d-20)cycle + if(dabs(ao_overlap_abs_grid(j,i)).lt.thrsh_cycle_tc)cycle tmp = 0.d0 do i_1s = 1, List_comb_thr_b2_size(j,i) @@ -49,7 +49,7 @@ BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_j1b_test, (ao_num, ao_num, 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 +! if(dabs(coef)*dabs(int_j1b).lt.thrsh_cycle_tc)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) @@ -110,7 +110,7 @@ BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_j1b_test, (ao_num, ao_nu !$OMP SHARED (n_points_final_grid, ao_num, List_comb_thr_b2_size, 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_j1b_test, mu_erf,ao_abs_comb_b2_j1b, & - !$OMP ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma) + !$OMP ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,thrsh_cycle_tc) ! !$OMP ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,dsqpi_3_2,expo_erfc_mu_gauss) !$OMP DO do ipoint = 1, n_points_final_grid @@ -120,7 +120,7 @@ BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_j1b_test, (ao_num, ao_nu do i = 1, ao_num do j = i, ao_num - if(dabs(ao_overlap_abs_grid(j,i)).lt.1.d-10)cycle + if(dabs(ao_overlap_abs_grid(j,i)).lt.thrsh_cycle_tc)cycle tmp_x = 0.d0 tmp_y = 0.d0 @@ -130,19 +130,11 @@ BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_j1b_test, (ao_num, ao_nu 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 + ! if(dabs(coef)*dabs(int_j1b).lt.thrsh_cycle_tc)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) -! if(ao_prod_center(1,j,i).ne.10000.d0)then -! ! approximate 1 - erf(mu r12) by a gaussian * 10 -! !DIR$ FORCEINLINE -! call gaussian_product(expo_erfc_mu_gauss,r, & -! ao_prod_sigma(j,i),ao_prod_center(1,j,i), & -! factor_ij_1s,beta_ij,center_ij_1s) -! if(dabs(coef * factor_ij_1s*int_j1b*10.d0 * dsqpi_3_2 * beta_ij**(-1.5d0)).lt.1.d-10)cycle -! endif 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) @@ -216,7 +208,7 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_test, (ao_num, ao_num, n_po !$OMP expo_gauss_j_mu_x, coef_gauss_j_mu_x, & !$OMP List_comb_thr_b2_coef, List_comb_thr_b2_expo,List_comb_thr_b2_size, & !$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 ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,dsqpi_3_2,thrsh_cycle_tc) !$OMP DO do ipoint = 1, n_points_final_grid r(1) = final_grid_points(1,ipoint) @@ -225,7 +217,7 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_test, (ao_num, ao_num, n_po do i = 1, ao_num do j = i, ao_num - if(dabs(ao_overlap_abs_grid(j,i)).lt.1.d-20)cycle + if(dabs(ao_overlap_abs_grid(j,i)).lt.thrsh_cycle_tc)cycle tmp = 0.d0 @@ -234,11 +226,11 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_test, (ao_num, ao_num, n_po ! --- --- --- int_j1b = ao_abs_comb_b2_j1b(1,j,i) - if(dabs(int_j1b).lt.1.d-10) cycle + ! if(dabs(int_j1b).lt.thrsh_cycle_tc) cycle do i_fit = 1, ng_fit_jast expo_fit = expo_gauss_j_mu_x(i_fit) coef_fit = coef_gauss_j_mu_x(i_fit) - if(ao_overlap_abs_grid(j,i).lt.1.d-15) cycle + ! if(ao_overlap_abs_grid(j,i).lt.thrsh_cycle_tc) cycle int_fit = overlap_gauss_r12_ao(r, expo_fit, i, j) tmp += coef_fit * int_fit enddo @@ -251,7 +243,7 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_test, (ao_num, ao_num, n_po 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 +! if(dabs(coef)*dabs(int_j1b).lt.thrsh_cycle_tc)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) @@ -259,9 +251,9 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_test, (ao_num, ao_num, n_po 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 +! if(dabs(coeftot).lt.thrsh_cycle_tc)cycle 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 +! if(factor_ij_1s_u*ao_overlap_abs_grid(j,i).lt.thrsh_cycle_tc)cycle int_fit = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j) tmp += coef * coef_fit * int_fit enddo @@ -325,7 +317,7 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_ng_1_test, (ao_num, ao_num, !$OMP expo_gauss_j_mu_x, coef_gauss_j_mu_x, & !$OMP List_comb_thr_b2_coef, List_comb_thr_b2_expo,List_comb_thr_b2_size, & !$OMP List_comb_thr_b2_cent, v_ij_u_cst_mu_j1b_ng_1_test,ao_abs_comb_b2_j1b, & - !$OMP ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,dsqpi_3_2) + !$OMP ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,dsqpi_3_2,thrsh_cycle_tc) !$OMP DO do ipoint = 1, n_points_final_grid r(1) = final_grid_points(1,ipoint) @@ -334,7 +326,7 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_ng_1_test, (ao_num, ao_num, do i = 1, ao_num do j = i, ao_num - if(dabs(ao_overlap_abs_grid(j,i)).lt.1.d-20)cycle + if(dabs(ao_overlap_abs_grid(j,i)).lt.thrsh_cycle_tc)cycle tmp = 0.d0 @@ -343,7 +335,7 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_ng_1_test, (ao_num, ao_num, ! --- --- --- int_j1b = ao_abs_comb_b2_j1b(1,j,i) - if(dabs(int_j1b).lt.1.d-10) cycle +! if(dabs(int_j1b).lt.thrsh_cycle_tc) cycle expo_fit = expo_good_j_mu_1gauss int_fit = overlap_gauss_r12_ao(r, expo_fit, i, j) tmp += int_fit @@ -356,7 +348,7 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_ng_1_test, (ao_num, ao_num, 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 +! if(dabs(coef)*dabs(int_j1b).lt.thrsh_cycle_tc)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) @@ -364,9 +356,9 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_ng_1_test, (ao_num, ao_num, expo_fit = expo_good_j_mu_1gauss coef_fit = 1.d0 coeftot = coef * coef_fit - if(dabs(coeftot).lt.1.d-15)cycle + if(dabs(coeftot).lt.thrsh_cycle_tc)cycle 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 + if(factor_ij_1s_u*ao_overlap_abs_grid(j,i).lt.thrsh_cycle_tc)cycle int_fit = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j) tmp += coef * coef_fit * int_fit ! enddo diff --git a/src/ao_many_one_e_ints/listj1b_sorted.irp.f b/src/ao_many_one_e_ints/listj1b_sorted.irp.f index bf493fbb..9bcce449 100644 --- a/src/ao_many_one_e_ints/listj1b_sorted.irp.f +++ b/src/ao_many_one_e_ints/listj1b_sorted.irp.f @@ -3,15 +3,16 @@ &BEGIN_PROVIDER [ integer, max_List_comb_thr_b2_size] implicit none integer :: i_1s,i,j,ipoint - double precision :: coef,beta,center(3),int_j1b,thr + double precision :: coef,beta,center(3),int_j1b double precision :: r(3),weight,dist - thr = 1.d-15 List_comb_thr_b2_size = 0 + print*,'List_all_comb_b2_size = ',List_all_comb_b2_size +! pause 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-15)cycle + if(dabs(coef).lt.thrsh_cycle_tc)cycle beta = List_all_comb_b2_expo (i_1s) beta = max(beta,1.d-12) center(1:3) = List_all_comb_b2_cent(1:3,i_1s) @@ -24,7 +25,7 @@ dist += ( center(3) - r(3) )*( center(3) - r(3) ) int_j1b += dabs(aos_in_r_array_extra_transp(ipoint,i) * aos_in_r_array_extra_transp(ipoint,j))*dexp(-beta*dist) * weight enddo - if(dabs(coef)*dabs(int_j1b).gt.thr)then + if(dabs(coef)*dabs(int_j1b).gt.thrsh_cycle_tc)then List_comb_thr_b2_size(j,i) += 1 endif enddo @@ -40,6 +41,7 @@ list(i) = maxval(List_comb_thr_b2_size(:,i)) enddo max_List_comb_thr_b2_size = maxval(list) + print*,'max_List_comb_thr_b2_size = ',max_List_comb_thr_b2_size END_PROVIDER @@ -49,16 +51,15 @@ END_PROVIDER &BEGIN_PROVIDER [ double precision, ao_abs_comb_b2_j1b, ( max_List_comb_thr_b2_size ,ao_num, ao_num)] implicit none integer :: i_1s,i,j,ipoint,icount - double precision :: coef,beta,center(3),int_j1b,thr + double precision :: coef,beta,center(3),int_j1b double precision :: r(3),weight,dist - thr = 1.d-15 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-12)cycle + if(dabs(coef).lt.thrsh_cycle_tc)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 @@ -70,7 +71,7 @@ END_PROVIDER dist += ( center(3) - r(3) )*( center(3) - r(3) ) int_j1b += dabs(aos_in_r_array_extra_transp(ipoint,i) * aos_in_r_array_extra_transp(ipoint,j))*dexp(-beta*dist) * weight enddo - if(dabs(coef)*dabs(int_j1b).gt.thr)then + if(dabs(coef)*dabs(int_j1b).gt.thrsh_cycle_tc)then icount += 1 List_comb_thr_b2_coef(icount,j,i) = coef List_comb_thr_b2_expo(icount,j,i) = beta @@ -98,17 +99,17 @@ END_PROVIDER &BEGIN_PROVIDER [ integer, max_List_comb_thr_b3_size] implicit none integer :: i_1s,i,j,ipoint - double precision :: coef,beta,center(3),int_j1b,thr + double precision :: coef,beta,center(3),int_j1b double precision :: r(3),weight,dist - thr = 1.d-15 List_comb_thr_b3_size = 0 + print*,'List_all_comb_b3_size = ',List_all_comb_b3_size do i = 1, ao_num do j = 1, ao_num 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) center(1:3) = List_all_comb_b3_cent(1:3,i_1s) - if(dabs(coef).lt.thr)cycle + if(dabs(coef).lt.thrsh_cycle_tc)cycle int_j1b = 0.d0 do ipoint = 1, n_points_extra_final_grid r(1:3) = final_grid_points_extra(1:3,ipoint) @@ -118,7 +119,7 @@ END_PROVIDER dist += ( center(3) - r(3) )*( center(3) - r(3) ) int_j1b += dabs(aos_in_r_array_extra_transp(ipoint,i) * aos_in_r_array_extra_transp(ipoint,j))*dexp(-beta*dist) * weight enddo - if(dabs(coef)*dabs(int_j1b).gt.thr)then + if(dabs(coef)*dabs(int_j1b).gt.thrsh_cycle_tc)then List_comb_thr_b3_size(j,i) += 1 endif enddo @@ -144,9 +145,8 @@ END_PROVIDER &BEGIN_PROVIDER [ double precision, ao_abs_comb_b3_j1b, ( max_List_comb_thr_b3_size ,ao_num, ao_num)] implicit none integer :: i_1s,i,j,ipoint,icount - double precision :: coef,beta,center(3),int_j1b,thr + double precision :: coef,beta,center(3),int_j1b double precision :: r(3),weight,dist - thr = 1.d-15 ao_abs_comb_b3_j1b = 10000000.d0 do i = 1, ao_num do j = 1, ao_num @@ -156,7 +156,7 @@ END_PROVIDER beta = List_all_comb_b3_expo (i_1s) beta = max(beta,1.d-12) center(1:3) = List_all_comb_b3_cent(1:3,i_1s) - if(dabs(coef).lt.thr)cycle + if(dabs(coef).lt.thrsh_cycle_tc)cycle int_j1b = 0.d0 do ipoint = 1, n_points_extra_final_grid r(1:3) = final_grid_points_extra(1:3,ipoint) @@ -166,7 +166,7 @@ END_PROVIDER dist += ( center(3) - r(3) )*( center(3) - r(3) ) int_j1b += dabs(aos_in_r_array_extra_transp(ipoint,i) * aos_in_r_array_extra_transp(ipoint,j))*dexp(-beta*dist) * weight enddo - if(dabs(coef)*dabs(int_j1b).gt.thr)then + if(dabs(coef)*dabs(int_j1b).gt.thrsh_cycle_tc)then icount += 1 List_comb_thr_b3_coef(icount,j,i) = coef List_comb_thr_b3_expo(icount,j,i) = beta @@ -177,15 +177,5 @@ END_PROVIDER enddo enddo -! do i = 1, ao_num -! do j = 1, i-1 -! do icount = 1, List_comb_thr_b3_size(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:3,icount,j,i) = List_comb_thr_b3_cent(1:3,icount,i,j) -! enddo -! enddo -! enddo - END_PROVIDER diff --git a/src/non_h_ints_mu/total_tc_int.irp.f b/src/non_h_ints_mu/total_tc_int.irp.f index 2fd2719c..4f8dc74d 100644 --- a/src/non_h_ints_mu/total_tc_int.irp.f +++ b/src/non_h_ints_mu/total_tc_int.irp.f @@ -68,7 +68,26 @@ BEGIN_PROVIDER [double precision, ao_tc_int_chemist, (ao_num, ao_num, ao_num, ao END_PROVIDER +BEGIN_PROVIDER [double precision, ao_tc_int_chemist_no_cycle, (ao_num, ao_num, ao_num, ao_num)] ! --- + implicit none + integer :: i, j, k, l + double precision :: wall1, wall0 + print *, ' providing ao_tc_int_chemist_no_cycle ...' + call wall_time(wall0) + do j = 1, ao_num + do l = 1, ao_num + do i = 1, ao_num + do k = 1, ao_num + ao_tc_int_chemist_no_cycle(k,i,l,j) = tc_grad_square_ao(k,i,l,j) + tc_grad_and_lapl_ao(k,i,l,j) + ao_two_e_coul(k,i,l,j) +! ao_tc_int_chemist(k,i,l,j) = ao_two_e_coul(k,i,l,j) + enddo + enddo + enddo + enddo + call wall_time(wall1) + print *, ' wall time for ao_tc_int_chemist_no_cycle ', wall1 - wall0 +END_PROVIDER BEGIN_PROVIDER [double precision, ao_tc_int_chemist_test, (ao_num, ao_num, ao_num, ao_num)] diff --git a/src/tc_keywords/EZFIO.cfg b/src/tc_keywords/EZFIO.cfg index e397e700..62b6d2bf 100644 --- a/src/tc_keywords/EZFIO.cfg +++ b/src/tc_keywords/EZFIO.cfg @@ -46,6 +46,12 @@ doc: Thresholds on the energy for iterative Davidson used in TC interface: ezfio,provider,ocaml default: 1.e-5 +[thrsh_cycle_tc] +type: Threshold +doc: Thresholds to cycle the integrals with the envelop +interface: ezfio,provider,ocaml +default: 1.e-15 + [max_it_dav] type: integer doc: nb max of iteration in Davidson used in TC diff --git a/src/tc_scf/tc_scf.irp.f b/src/tc_scf/tc_scf.irp.f index 85389f30..ae3b609b 100644 --- a/src/tc_scf/tc_scf.irp.f +++ b/src/tc_scf/tc_scf.irp.f @@ -11,8 +11,8 @@ program tc_scf print *, ' starting ...' my_grid_becke = .True. - my_n_pt_r_grid = 30 - my_n_pt_a_grid = 50 + my_n_pt_r_grid = 60 + my_n_pt_a_grid = 110 ! 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 diff --git a/src/tc_scf/test_int.irp.f b/src/tc_scf/test_int.irp.f index 0866cdaf..b9287d58 100644 --- a/src/tc_scf/test_int.irp.f +++ b/src/tc_scf/test_int.irp.f @@ -21,25 +21,22 @@ program test_ints touch my_extra_grid_becke my_n_pt_r_extra_grid my_n_pt_a_extra_grid !! OK -!call routine_int2_u_grad1u_j1b2 -!! OK -!call routine_v_ij_erf_rk_cst_mu_j1b -!! OK +! call routine_int2_u_grad1u_j1b2 +! OK +! call routine_v_ij_erf_rk_cst_mu_j1b +! OK ! call routine_x_v_ij_erf_rk_cst_mu_j1b -!! OK -! call routine_v_ij_u_cst_mu_j1b - -!! OK -!call routine_int2_u2_j1b2 - -!! OK -!call routine_int2_u_grad1u_x_j1b2 - -!! OK +! OK +! call routine_int2_u2_j1b2 +! OK +! call routine_int2_u_grad1u_x_j1b2 +! OK ! call routine_int2_grad1u2_grad2u2_j1b2 ! call routine_int2_u_grad1u_j1b2 ! call test_total_grad_lapl ! call test_total_grad_square +! call test_int2_grad1_u12_ao_test +! call routine_v_ij_u_cst_mu_j1b_test ! call test_ao_tc_int_chemist ! call test_grid_points_ao ! call test_tc_scf @@ -53,12 +50,12 @@ program test_ints !call test_two_e_tc_non_hermit_integral() - call test_tc_grad_square_ao_test() - - PROVIDE TC_HF_energy VARTC_HF_energy - print *, ' TC_HF_energy = ', TC_HF_energy - print *, ' VARTC_HF_energy = ', VARTC_HF_energy +! call test_tc_grad_square_ao_test() +!!PROVIDE TC_HF_energy VARTC_HF_energy +!!print *, ' TC_HF_energy = ', TC_HF_energy +!!print *, ' VARTC_HF_energy = ', VARTC_HF_energy + call test_old_ints end ! --- @@ -157,6 +154,9 @@ subroutine routine_int2_u_grad1u_j1b2 enddo enddo enddo + print*,'******' + print*,'******' + print*,'routine_int2_u_grad1u_j1b2' print*,'accu_abs = ',accu_abs/dble(ao_num)**4 print*,'accu_relat = ',accu_relat/dble(ao_num)**4 @@ -169,20 +169,6 @@ subroutine routine_v_ij_erf_rk_cst_mu_j1b 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)) @@ -215,6 +201,9 @@ subroutine routine_v_ij_erf_rk_cst_mu_j1b enddo enddo enddo + print*,'******' + print*,'******' + print*,'routine_v_ij_erf_rk_cst_mu_j1b' print*,'accu_abs = ',accu_abs/dble(ao_num)**4 print*,'accu_relat = ',accu_relat/dble(ao_num)**4 @@ -228,20 +217,6 @@ subroutine routine_x_v_ij_erf_rk_cst_mu_j1b 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)) @@ -276,6 +251,10 @@ subroutine routine_x_v_ij_erf_rk_cst_mu_j1b enddo enddo enddo + + print*,'******' + print*,'******' + print*,'routine_x_v_ij_erf_rk_cst_mu_j1b' print*,'accu_abs = ',accu_abs/dble(ao_num)**4 print*,'accu_relat = ',accu_relat/dble(ao_num)**4 @@ -290,20 +269,6 @@ subroutine routine_v_ij_u_cst_mu_j1b_test 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)) @@ -336,6 +301,9 @@ subroutine routine_v_ij_u_cst_mu_j1b_test enddo enddo enddo + print*,'******' + print*,'******' + print*,'routine_v_ij_u_cst_mu_j1b_test' print*,'accu_abs = ',accu_abs/dble(ao_num)**4 print*,'accu_relat = ',accu_relat/dble(ao_num)**4 @@ -427,19 +395,6 @@ subroutine routine_int2_u2_j1b2 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 @@ -473,6 +428,9 @@ subroutine routine_int2_u2_j1b2 enddo enddo enddo + print*,'******' + print*,'******' + print*,'routine_int2_u2_j1b2' print*,'accu_abs = ',accu_abs/dble(ao_num)**4 print*,'accu_relat = ',accu_relat/dble(ao_num)**4 @@ -486,19 +444,6 @@ subroutine routine_int2_u_grad1u_x_j1b2 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 @@ -534,6 +479,9 @@ subroutine routine_int2_u_grad1u_x_j1b2 enddo enddo enddo + print*,'******' + print*,'******' + print*,'routine_int2_u_grad1u_x_j1b2' print*,'accu_abs = ',accu_abs/dble(ao_num)**4 print*,'accu_relat = ',accu_relat/dble(ao_num)**4 @@ -579,6 +527,9 @@ subroutine routine_v_ij_u_cst_mu_j1b enddo enddo enddo + print*,'******' + print*,'******' + print*,'routine_v_ij_u_cst_mu_j1b' print*,'accu_abs = ',accu_abs/dble(ao_num)**4 print*,'accu_relat = ',accu_relat/dble(ao_num)**4 @@ -742,6 +693,9 @@ subroutine test_total_grad_lapl enddo enddo enddo + print*,'******' + print*,'******' + print*,' test_total_grad_lapl' print*,'accu_abs = ',accu_abs/dble(ao_num)**4 print*,'accu_relat = ',accu_relat/dble(ao_num)**4 @@ -767,6 +721,9 @@ subroutine test_total_grad_square enddo enddo enddo + print*,'******' + print*,'******' + print*,'test_total_grad_square' print*,'accu_abs = ',accu_abs/dble(ao_num)**4 print*,'accu_relat = ',accu_relat/dble(ao_num)**4 @@ -1057,3 +1014,86 @@ end ! --- + +subroutine test_old_ints + implicit none + integer :: i,j,k,l + double precision :: old, new, contrib, get_ao_tc_sym_two_e_pot + double precision :: integral_sym , integral_nsym,accu + PROVIDE ao_tc_sym_two_e_pot_in_map + accu = 0.d0 + do j = 1, ao_num + do l= 1, ao_num + do i = 1, ao_num + do k = 1, ao_num +! integral_sym = get_ao_tc_sym_two_e_pot(i, j, k, l, ao_tc_sym_two_e_pot_map) + ! ao_non_hermit_term_chemist(k,i,l,j) = < k l | [erf( mu r12) - 1] d/d_r12 | i j > on the AO basis +! integral_nsym = ao_non_hermit_term_chemist(k,i,l,j) +! old = integral_sym + integral_nsym +! old = tc_grad_square_ao(k,i,l,j) + tc_grad_and_lapl_ao(k,i,l,j) + ao_two_e_coul(k,i,l,j) + new = ao_tc_int_chemist_test(k,i,l,j) + old = ao_tc_int_chemist_no_cycle(k,i,l,j) + contrib = dabs(old - new) + if(contrib.gt.1.d-6)then + print*,'problem !!' + print*,i,j,k,l + print*,old, new, contrib + endif + accu += contrib + enddo + enddo + enddo + enddo + print*,'******' + print*,'******' + print*,'in test_old_ints' + print*,'accu = ',accu/dble(ao_num**4) + +end + +subroutine test_int2_grad1_u12_ao_test + implicit none + integer :: i,j,ipoint,m,k,l + double precision :: weight,accu_relat, accu_abs, contrib + double precision, allocatable :: array(:,:,:,:), array_ref(:,:,:,:) + 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 m = 1, 3 + 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) += int2_grad1_u12_ao_test(j,i,ipoint,m) * aos_grad_in_r_array_transp(m,k,ipoint) * aos_in_r_array(l,ipoint) * weight + array_ref(j,i,l,k) += int2_grad1_u12_ao(j,i,ipoint,m) * 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*,'******' + print*,'******' + print*,'test_int2_grad1_u12_ao_test' + print*,'accu_abs = ',accu_abs/dble(ao_num)**4 + print*,'accu_relat = ',accu_relat/dble(ao_num)**4 +end diff --git a/src/tools/sort_wf.irp.f b/src/tools/sort_wf.irp.f index 038b24c8..95b1a964 100644 --- a/src/tools/sort_wf.irp.f +++ b/src/tools/sort_wf.irp.f @@ -14,7 +14,7 @@ subroutine routine output=trim(ezfio_filename)//'.wf_sorted' i_unit_output = getUnitAndOpen(output,'w') do i = 1, N_det - write(i_unit_output, *)i,dabs(psi_coef_sorted(i,1)) + write(i_unit_output, *)i,dabs(psi_coef_sorted(i,1))/dabs(psi_coef_sorted(1,1)) enddo end