diff --git a/src/ao_many_one_e_ints/ao_erf_gauss.irp.f b/src/ao_many_one_e_ints/ao_erf_gauss.irp.f index 3d7fbe50..b1077161 100644 --- a/src/ao_many_one_e_ints/ao_erf_gauss.irp.f +++ b/src/ao_many_one_e_ints/ao_erf_gauss.irp.f @@ -212,9 +212,7 @@ subroutine NAI_pol_x_mult_erf_ao(i_ao, j_ao, mu_in, C_center, ints) ! Computes the following integral : ! ! $\int_{-\infty}^{infty} dr x * \chi_i(r) \chi_j(r) \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$. - ! ! $\int_{-\infty}^{infty} dr y * \chi_i(r) \chi_j(r) \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$. - ! ! $\int_{-\infty}^{infty} dr z * \chi_i(r) \chi_j(r) \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$. ! END_DOC @@ -279,9 +277,7 @@ subroutine NAI_pol_x_mult_erf_ao_v0(i_ao, j_ao, mu_in, C_center, LD_C, ints, LD_ ! Computes the following integral : ! ! $\int_{-\infty}^{infty} dr x * \chi_i(r) \chi_j(r) \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$. - ! ! $\int_{-\infty}^{infty} dr y * \chi_i(r) \chi_j(r) \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$. - ! ! $\int_{-\infty}^{infty} dr z * \chi_i(r) \chi_j(r) \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$. ! END_DOC @@ -1111,3 +1107,141 @@ end ! --- +subroutine NAI_pol_x2_mult_erf_ao_with1s(i_ao, j_ao, beta, B_center, mu_in, C_center, ints) + + BEGIN_DOC + ! + ! Computes the following integral : + ! + ! $\int_{-\infty}^{infty} dr x^2 * \chi_i(r) \chi_j(r) e^{-\beta (r - B_center)^2} \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$. + ! $\int_{-\infty}^{infty} dr y^2 * \chi_i(r) \chi_j(r) e^{-\beta (r - B_center)^2} \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$. + ! $\int_{-\infty}^{infty} dr z^2 * \chi_i(r) \chi_j(r) e^{-\beta (r - B_center)^2} \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$. + ! + END_DOC + + include 'utils/constants.include.F' + + implicit none + + integer, intent(in) :: i_ao, j_ao + double precision, intent(in) :: beta, B_center(3), mu_in, C_center(3) + double precision, intent(out) :: ints(3) + + integer :: i, j, power_Ai(3), power_Aj(3), n_pt_in, m + integer :: power_A1(3), power_A2(3) + double precision :: Ai_center(3), Aj_center(3), alphai, alphaj, coef, coefi + double precision :: integral0, integral1, integral2 + + double precision, external :: NAI_pol_mult_erf_with1s + + ASSERT(beta .ge. 0.d0) + if(beta .lt. 1d-10) then + call NAI_pol_x2_mult_erf_ao(i_ao, j_ao, mu_in, C_center, ints) + return + endif + + ints = 0.d0 + + power_Ai(1:3) = ao_power(i_ao,1:3) + power_Aj(1:3) = ao_power(j_ao,1:3) + + Ai_center(1:3) = nucl_coord(ao_nucl(i_ao),1:3) + Aj_center(1:3) = nucl_coord(ao_nucl(j_ao),1:3) + + n_pt_in = n_pt_max_integrals + + do i = 1, ao_prim_num(i_ao) + alphai = ao_expo_ordered_transp (i,i_ao) + coefi = ao_coef_normalized_ordered_transp(i,i_ao) + + do m = 1, 3 + + power_A1 = power_Ai + power_A1(m) += 1 + + power_A2 = power_Ai + power_A2(m) += 2 + + do j = 1, ao_prim_num(j_ao) + alphaj = ao_expo_ordered_transp (j,j_ao) + coef = coefi * ao_coef_normalized_ordered_transp(j,j_ao) + + integral0 = NAI_pol_mult_erf_with1s(Ai_center, Aj_center, power_Ai, power_Aj, alphai, alphaj, beta, B_center, C_center, n_pt_in, mu_in) + integral1 = NAI_pol_mult_erf_with1s(Ai_center, Aj_center, power_A1, power_Aj, alphai, alphaj, beta, B_center, C_center, n_pt_in, mu_in) + integral2 = NAI_pol_mult_erf_with1s(Ai_center, Aj_center, power_A2, power_Aj, alphai, alphaj, beta, B_center, C_center, n_pt_in, mu_in) + + ints(m) += coef * (integral2 + Ai_center(m) * (2.d0*integral1 + Ai_center(m)*integral0)) + enddo + enddo + enddo + +end subroutine NAI_pol_x2_mult_erf_ao_with1s + +! --- + +subroutine NAI_pol_x2_mult_erf_ao(i_ao, j_ao, mu_in, C_center, ints) + + BEGIN_DOC + ! + ! Computes the following integral : + ! + ! $\int_{-\infty}^{infty} dr x^2 * \chi_i(r) \chi_j(r) \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$. + ! $\int_{-\infty}^{infty} dr y^2 * \chi_i(r) \chi_j(r) \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$. + ! $\int_{-\infty}^{infty} dr z^2 * \chi_i(r) \chi_j(r) \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$. + ! + END_DOC + + include 'utils/constants.include.F' + + implicit none + + integer, intent(in) :: i_ao, j_ao + double precision, intent(in) :: mu_in, C_center(3) + double precision, intent(out) :: ints(3) + + integer :: i, j, num_A, num_B, power_A(3), power_B(3), n_pt_in, m + integer :: power_A1(3), power_A2(3) + double precision :: A_center(3), B_center(3), alpha, beta, coef + double precision :: integral0, integral1, integral2 + + double precision :: NAI_pol_mult_erf + + ints = 0.d0 + + num_A = ao_nucl(i_ao) + power_A(1:3) = ao_power(i_ao,1:3) + A_center(1:3) = nucl_coord(num_A,1:3) + num_B = ao_nucl(j_ao) + power_B(1:3) = ao_power(j_ao,1:3) + B_center(1:3) = nucl_coord(num_B,1:3) + + n_pt_in = n_pt_max_integrals + + do i = 1, ao_prim_num(i_ao) + alpha = ao_expo_ordered_transp(i,i_ao) + + do m = 1, 3 + + power_A1 = power_A + power_A1(m) += 1 + + power_A2 = power_A + power_A2(m) += 2 + + do j = 1, ao_prim_num(j_ao) + beta = ao_expo_ordered_transp(j,j_ao) + coef = ao_coef_normalized_ordered_transp(j,j_ao) * ao_coef_normalized_ordered_transp(i,i_ao) + + integral0 = NAI_pol_mult_erf(A_center, B_center, power_A , power_B, alpha, beta, C_center, n_pt_in, mu_in) + integral1 = NAI_pol_mult_erf(A_center, B_center, power_A1, power_B, alpha, beta, C_center, n_pt_in, mu_in) + integral2 = NAI_pol_mult_erf(A_center, B_center, power_A2, power_B, alpha, beta, C_center, n_pt_in, mu_in) + + ints(m) += coef * (integral2 + A_center(m) * (2.d0*integral1 + A_center(m)*integral0)) + enddo + enddo + enddo + +end subroutine NAI_pol_x2_mult_erf_ao + +! --- + diff --git a/src/ao_many_one_e_ints/grad_lapl_jmu_modif.irp.f b/src/ao_many_one_e_ints/grad_lapl_jmu_modif.irp.f index 25bb2f8b..9d34e1d7 100644 --- a/src/ao_many_one_e_ints/grad_lapl_jmu_modif.irp.f +++ b/src/ao_many_one_e_ints/grad_lapl_jmu_modif.irp.f @@ -195,7 +195,6 @@ END_PROVIDER ! --- -! TODO analytically BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b, (ao_num, ao_num, n_points_final_grid)] BEGIN_DOC @@ -217,6 +216,8 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b, (ao_num, ao_num, n_points_ call wall_time(wall0) provide mu_erf final_grid_points j1b_pen + PROVIDE ng_fit_jast expo_gauss_j_mu_x coef_gauss_j_mu_x + PROVIDE List_all_comb_b2_size List_all_comb_b2_coef List_all_comb_b2_expo List_all_comb_b2_cent v_ij_u_cst_mu_j1b = 0.d0 @@ -229,7 +230,6 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b, (ao_num, ao_num, n_points_ !$OMP List_all_comb_b2_coef, List_all_comb_b2_expo, & !$OMP List_all_comb_b2_cent, v_ij_u_cst_mu_j1b) !$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) @@ -240,10 +240,13 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b, (ao_num, ao_num, n_points_ tmp = 0.d0 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) +! do i_fit = ng_fit_jast, ng_fit_jast +! expo_fit = 5.0d0 +! coef_fit = 1.0d0 + ! --- coef = List_all_comb_b2_coef (1) @@ -253,7 +256,6 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b, (ao_num, ao_num, n_points_ B_center(3) = List_all_comb_b2_cent(3,1) int_fit = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j) -! if(dabs(int_fit*coef) .lt. 1d-12) cycle tmp += coef * coef_fit * int_fit @@ -298,3 +300,137 @@ END_PROVIDER ! --- +BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_an, (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 + + include 'constants.include.F' + + implicit none + integer :: i, j, ipoint, i_1s + double precision :: r(3), r1_2 + double precision :: int_c1, int_e1, int_o + double precision :: int_c2(3), int_e2(3) + double precision :: int_c3(3), int_e3(3) + double precision :: coef, beta, B_center(3) + double precision :: tmp, ct + double precision :: wall0, wall1 + + double precision, external :: overlap_gauss_r12_ao_with1s + double precision, external :: NAI_pol_mult_erf_ao_with1s + + print*, ' providing v_ij_u_cst_mu_j1b_an ...' + call wall_time(wall0) + + provide mu_erf final_grid_points j1b_pen + PROVIDE ng_fit_jast expo_gauss_j_mu_x coef_gauss_j_mu_x + PROVIDE List_all_comb_b2_size List_all_comb_b2_coef List_all_comb_b2_expo List_all_comb_b2_cent + + ct = inv_sq_pi_2 / mu_erf + + v_ij_u_cst_mu_j1b_an = 0.d0 + + !$OMP PARALLEL DEFAULT (NONE) & + !$OMP PRIVATE (ipoint, i, j, i_1s, r, coef, beta, B_center, & + !$OMP r1_2, tmp, int_c1, int_e1, int_o, int_c2, & + !$OMP int_e2, int_c3, int_e3) & + !$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b2_size, & + !$OMP final_grid_points, mu_erf, ct, & + !$OMP expo_gauss_j_mu_x, coef_gauss_j_mu_x, & + !$OMP List_all_comb_b2_coef, List_all_comb_b2_expo, & + !$OMP List_all_comb_b2_cent, v_ij_u_cst_mu_j1b_an) + !$OMP DO + 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) + r1_2 = 0.5d0 * (r(1)*r(1) + r(2)*r(2) + r(3)*r(3)) + + do i = 1, ao_num + do j = i, ao_num + + ! --- + + coef = List_all_comb_b2_coef (1) + beta = List_all_comb_b2_expo (1) + B_center(1) = List_all_comb_b2_cent(1,1) + B_center(2) = List_all_comb_b2_cent(2,1) + B_center(3) = List_all_comb_b2_cent(3,1) + + int_c1 = NAI_pol_mult_erf_ao_with1s(i, j, beta, B_center, 1.d+9, r) + int_e1 = NAI_pol_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r) + + call NAI_pol_x_mult_erf_ao_with1s(i, j, beta, B_center, 1.d+9, r, int_c2) + call NAI_pol_x_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r, int_e2) + + call NAI_pol_x2_mult_erf_ao_with1s(i, j, beta, B_center, 1.d+9, r, int_c3) + call NAI_pol_x2_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r, int_e3) + + int_o = overlap_gauss_r12_ao_with1s(B_center, beta, r, mu_erf*mu_erf, i, j) + + tmp = coef & + * ( r1_2 * (int_c1 - int_e1) & + - r(1) * (int_c2(1) - int_e2(1)) - r(2) * (int_c2(2) - int_e2(2)) - r(3) * (int_c2(3) - int_e2(3)) & + + 0.5d0 * (int_c3(1) + int_c3(2) + int_c3(3) - int_e3(1) - int_e3(2) - int_e3(3)) & + - ct * int_o & + ) + + ! --- + + do i_1s = 2, List_all_comb_b2_size + + coef = List_all_comb_b2_coef (i_1s) + beta = List_all_comb_b2_expo (i_1s) + B_center(1) = List_all_comb_b2_cent(1,i_1s) + B_center(2) = List_all_comb_b2_cent(2,i_1s) + B_center(3) = List_all_comb_b2_cent(3,i_1s) + + int_c1 = NAI_pol_mult_erf_ao_with1s(i, j, beta, B_center, 1.d+9, r) + int_e1 = NAI_pol_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r) + + call NAI_pol_x_mult_erf_ao_with1s(i, j, beta, B_center, 1.d+9, r, int_c2) + call NAI_pol_x_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r, int_e2) + + call NAI_pol_x2_mult_erf_ao_with1s(i, j, beta, B_center, 1.d+9, r, int_c3) + call NAI_pol_x2_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r, int_e3) + + int_o = overlap_gauss_r12_ao_with1s(B_center, beta, r, mu_erf*mu_erf, i, j) + + tmp = tmp + coef & + * ( r1_2 * (int_c1 - int_e1) & + - r(1) * (int_c2(1) - int_e2(1)) - r(2) * (int_c2(2) - int_e2(2)) - r(3) * (int_c2(3) - int_e2(3)) & + + 0.5d0 * (int_c3(1) + int_c3(2) + int_c3(3) - int_e3(1) - int_e3(2) - int_e3(3)) & + - ct * int_o & + ) + + enddo + + ! --- + + v_ij_u_cst_mu_j1b_an(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_an(j,i,ipoint) = v_ij_u_cst_mu_j1b_an(i,j,ipoint) + enddo + enddo + enddo + + call wall_time(wall1) + print*, ' wall time for v_ij_u_cst_mu_j1b_an', wall1 - wall0 + +END_PROVIDER + +! --- + diff --git a/src/ao_tc_eff_map/fit_j.irp.f b/src/ao_tc_eff_map/fit_j.irp.f index 4730d003..0fc3da2f 100644 --- a/src/ao_tc_eff_map/fit_j.irp.f +++ b/src/ao_tc_eff_map/fit_j.irp.f @@ -36,16 +36,25 @@ END_PROVIDER END_PROVIDER BEGIN_PROVIDER [ double precision, expo_j_xmu, (n_fit_1_erf_x) ] - implicit none - BEGIN_DOC - ! F(x) = x * (1 - erf(x)) - 1/sqrt(pi) * exp(-x**2) is fitted with a gaussian and a Slater - ! - ! \approx - 1/sqrt(pi) * exp(-alpha * x ) exp(-beta * x**2) - ! - ! where alpha = expo_j_xmu(1) and beta = expo_j_xmu(2) - END_DOC - expo_j_xmu(1) = 1.7477d0 - expo_j_xmu(2) = 0.668662d0 + + BEGIN_DOC + ! F(x) = x * (1 - erf(x)) - 1/sqrt(pi) * exp(-x**2) is fitted with a gaussian and a Slater + ! + ! \approx - 1/sqrt(pi) * exp(-alpha * x ) exp(-beta * x**2) + ! + ! where alpha = expo_j_xmu(1) and beta = expo_j_xmu(2) + END_DOC + + implicit none + + !expo_j_xmu(1) = 1.7477d0 + !expo_j_xmu(2) = 0.668662d0 + + !expo_j_xmu(1) = 1.74766377595541d0 + !expo_j_xmu(2) = 0.668719925486403d0 + + expo_j_xmu(1) = 1.74770446934522d0 + expo_j_xmu(2) = 0.668659706559979d0 END_PROVIDER diff --git a/src/non_h_ints_mu/tc_integ.irp.f b/src/non_h_ints_mu/tc_integ.irp.f index b2c0df31..ce65b203 100644 --- a/src/non_h_ints_mu/tc_integ.irp.f +++ b/src/non_h_ints_mu/tc_integ.irp.f @@ -70,14 +70,14 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_ao, (ao_num, ao_num, n_points_f elseif((j1b_type .eq. 3) .or. (j1b_type .eq. 4)) then - PROVIDE v_1b_grad v_ij_erf_rk_cst_mu_j1b v_ij_u_cst_mu_j1b x_v_ij_erf_rk_cst_mu_j1b + PROVIDE v_1b_grad v_ij_erf_rk_cst_mu_j1b v_ij_u_cst_mu_j1b_an x_v_ij_erf_rk_cst_mu_j1b int2_grad1_u12_ao = 0.d0 !$OMP PARALLEL & !$OMP DEFAULT (NONE) & !$OMP PRIVATE (ipoint, i, j, x, y, z, tmp0, tmp1, tmp2, tmp_x, tmp_y, tmp_z) & !$OMP SHARED ( ao_num, n_points_final_grid, final_grid_points, v_1b, v_1b_grad & - !$OMP , v_ij_erf_rk_cst_mu_j1b, v_ij_u_cst_mu_j1b, x_v_ij_erf_rk_cst_mu_j1b, int2_grad1_u12_ao) + !$OMP , v_ij_erf_rk_cst_mu_j1b, v_ij_u_cst_mu_j1b_an, x_v_ij_erf_rk_cst_mu_j1b, int2_grad1_u12_ao) !$OMP DO SCHEDULE (static) do ipoint = 1, n_points_final_grid x = final_grid_points(1,ipoint) @@ -90,7 +90,7 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_ao, (ao_num, ao_num, n_points_f do j = 1, ao_num do i = 1, ao_num tmp1 = tmp0 * v_ij_erf_rk_cst_mu_j1b(i,j,ipoint) - tmp2 = v_ij_u_cst_mu_j1b(i,j,ipoint) + tmp2 = v_ij_u_cst_mu_j1b_an(i,j,ipoint) int2_grad1_u12_ao(i,j,ipoint,1) = tmp1 * x - tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,1) - tmp2 * tmp_x int2_grad1_u12_ao(i,j,ipoint,2) = tmp1 * y - tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,2) - tmp2 * tmp_y int2_grad1_u12_ao(i,j,ipoint,3) = tmp1 * z - tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,3) - tmp2 * tmp_z @@ -100,7 +100,7 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_ao, (ao_num, ao_num, n_points_f !$OMP END DO !$OMP END PARALLEL - FREE v_ij_erf_rk_cst_mu_j1b v_ij_u_cst_mu_j1b x_v_ij_erf_rk_cst_mu_j1b + FREE v_ij_erf_rk_cst_mu_j1b v_ij_u_cst_mu_j1b_an x_v_ij_erf_rk_cst_mu_j1b elseif(j1b_type .ge. 100) then diff --git a/src/tc_bi_ortho/print_tc_energy.irp.f b/src/tc_bi_ortho/print_tc_energy.irp.f index 2f667a48..522f4cd7 100644 --- a/src/tc_bi_ortho/print_tc_energy.irp.f +++ b/src/tc_bi_ortho/print_tc_energy.irp.f @@ -9,11 +9,11 @@ program print_tc_energy print *, 'Hello world' 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 = 100 + my_n_pt_a_grid = 170 !my_n_pt_r_grid = 100 !my_n_pt_a_grid = 266 diff --git a/src/utils/integration.irp.f b/src/utils/integration.irp.f index b548b18a..72029c73 100644 --- a/src/utils/integration.irp.f +++ b/src/utils/integration.irp.f @@ -418,7 +418,7 @@ subroutine gaussian_product_x(a,xa,b,xb,k,p,xp) xab = xa-xb ab = ab*p_inv k = ab*xab*xab - if (k > 40.d0) then + if (k > 400.d0) then k=0.d0 return endif