From b13a315cc11c3223d21f7ee7093ce682b2a8a27a Mon Sep 17 00:00:00 2001 From: AbdAmmar Date: Mon, 17 Oct 2022 17:51:46 +0200 Subject: [PATCH] integrals over r2 tested --- src/ao_many_one_e_ints/grad2_jmu_modif.irp.f | 38 +- .../grad_lapl_jmu_modif.irp.f | 3 + ...j1b_product_to_sum.irp.f => listj1b.irp.f} | 0 src/ao_one_e_ints/pot_ao_erf_ints.irp.f | 2 +- src/non_h_ints_mu/debug_integ_jmu_modif.irp.f | 475 +++++++++++++++++- src/non_h_ints_mu/fit_j.irp.f | 28 +- src/non_h_ints_mu/grad_squared.irp.f | 2 +- src/non_h_ints_mu/j12_nucl_utils.irp.f | 89 ++-- src/non_h_ints_mu/numerical_integ.irp.f | 269 +++++++++- 9 files changed, 810 insertions(+), 96 deletions(-) rename src/ao_many_one_e_ints/{listj1b_product_to_sum.irp.f => listj1b.irp.f} (100%) diff --git a/src/ao_many_one_e_ints/grad2_jmu_modif.irp.f b/src/ao_many_one_e_ints/grad2_jmu_modif.irp.f index 163f6e2d..7e08bd97 100644 --- a/src/ao_many_one_e_ints/grad2_jmu_modif.irp.f +++ b/src/ao_many_one_e_ints/grad2_jmu_modif.irp.f @@ -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 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 9dd715e2..b847a630 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 @@ -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 diff --git a/src/ao_many_one_e_ints/listj1b_product_to_sum.irp.f b/src/ao_many_one_e_ints/listj1b.irp.f similarity index 100% rename from src/ao_many_one_e_ints/listj1b_product_to_sum.irp.f rename to src/ao_many_one_e_ints/listj1b.irp.f diff --git a/src/ao_one_e_ints/pot_ao_erf_ints.irp.f b/src/ao_one_e_ints/pot_ao_erf_ints.irp.f index e0d6254a..1d2d8faf 100644 --- a/src/ao_one_e_ints/pot_ao_erf_ints.irp.f +++ b/src/ao_one_e_ints/pot_ao_erf_ints.irp.f @@ -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 ) diff --git a/src/non_h_ints_mu/debug_integ_jmu_modif.irp.f b/src/non_h_ints_mu/debug_integ_jmu_modif.irp.f index 2bcd1358..e59b5f7a 100644 --- a/src/non_h_ints_mu/debug_integ_jmu_modif.irp.f +++ b/src/non_h_ints_mu/debug_integ_jmu_modif.irp.f @@ -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 + +! --- + diff --git a/src/non_h_ints_mu/fit_j.irp.f b/src/non_h_ints_mu/fit_j.irp.f index 34f3a31a..d53209d0 100644 --- a/src/non_h_ints_mu/fit_j.irp.f +++ b/src/non_h_ints_mu/fit_j.irp.f @@ -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 diff --git a/src/non_h_ints_mu/grad_squared.irp.f b/src/non_h_ints_mu/grad_squared.irp.f index 6ef34118..b9c98ea9 100644 --- a/src/non_h_ints_mu/grad_squared.irp.f +++ b/src/non_h_ints_mu/grad_squared.irp.f @@ -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 diff --git a/src/non_h_ints_mu/j12_nucl_utils.irp.f b/src/non_h_ints_mu/j12_nucl_utils.irp.f index 344b1fa8..a6dd0939 100644 --- a/src/non_h_ints_mu/j12_nucl_utils.irp.f +++ b/src/non_h_ints_mu/j12_nucl_utils.irp.f @@ -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 ! --- diff --git a/src/non_h_ints_mu/numerical_integ.irp.f b/src/non_h_ints_mu/numerical_integ.irp.f index 842851aa..17b666aa 100644 --- a/src/non_h_ints_mu/numerical_integ.irp.f +++ b/src/non_h_ints_mu/numerical_integ.irp.f @@ -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