From dbbae1f990f2ec96beeecdb9513aab2d6f90e1f7 Mon Sep 17 00:00:00 2001 From: AbdAmmar Date: Mon, 17 Oct 2022 10:55:53 +0200 Subject: [PATCH] list2, list3, j1b, grad_j1, lapl_j1b ok --- src/ao_many_one_e_ints/grad2_jmu_modif.irp.f | 118 ---- .../grad_lapl_jmu_modif.irp.f | 105 --- .../listj1b_product_to_sum.irp.f | 227 +++++++ src/ao_one_e_ints/pot_ao_erf_ints.irp.f | 46 +- src/non_h_ints_mu/debug_integ_jmu_modif.irp.f | 339 +++++++++- src/non_h_ints_mu/grad_squared.irp.f | 3 +- src/non_h_ints_mu/j12_nucl_utils.irp.f | 598 ++++++++++++++++++ src/non_h_ints_mu/jmu_modif.irp.f | 266 -------- src/non_h_ints_mu/new_grad_tc.irp.f | 86 +-- src/non_h_ints_mu/numerical_integ.irp.f | 165 +++-- 10 files changed, 1253 insertions(+), 700 deletions(-) create mode 100644 src/ao_many_one_e_ints/listj1b_product_to_sum.irp.f create mode 100644 src/non_h_ints_mu/j12_nucl_utils.irp.f delete mode 100644 src/non_h_ints_mu/jmu_modif.irp.f 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 04d4efb9..163f6e2d 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,124 +1,6 @@ ! --- -BEGIN_PROVIDER [ integer, List_all_comb_b3_size] - - implicit none - - List_all_comb_b3_size = 3**nucl_num - -END_PROVIDER - -! --- - -BEGIN_PROVIDER [ integer, List_all_comb_b3, (nucl_num, List_all_comb_b3_size)] - - implicit none - integer :: i, j, ii, jj - integer, allocatable :: M(:,:), p(:) - - if(nucl_num .gt. 32) then - print *, ' nucl_num = ', nucl_num, '> 32' - stop - endif - - List_all_comb_b3(:,:) = 0 - List_all_comb_b3(:,List_all_comb_b3_size) = 2 - - allocate(p(nucl_num)) - p = 0 - - do i = 2, List_all_comb_b3_size-1 - do j = 1, nucl_num - - ii = 0 - do jj = 1, j-1, 1 - ii = ii + p(jj) * 3**(jj-1) - enddo - p(j) = modulo(i-1-ii, 3**j) / 3**(j-1) - - List_all_comb_b3(j,i) = p(j) - enddo - enddo - -END_PROVIDER - -! --- - - BEGIN_PROVIDER [ double precision, List_all_comb_b3_coef, ( List_all_comb_b3_size)] -&BEGIN_PROVIDER [ double precision, List_all_comb_b3_expo, ( List_all_comb_b3_size)] -&BEGIN_PROVIDER [ double precision, List_all_comb_b3_cent, (3, List_all_comb_b3_size)] - - implicit none - integer :: i, j, k, phase - double precision :: tmp_alphaj, tmp_alphak, facto - - provide j1b_pen - - List_all_comb_b3_coef = 0.d0 - List_all_comb_b3_expo = 0.d0 - List_all_comb_b3_cent = 0.d0 - - do i = 1, List_all_comb_b3_size - - do j = 1, nucl_num - tmp_alphaj = dble(List_all_comb_b3(j,i)) * j1b_pen(j) - - List_all_comb_b3_expo(i) += tmp_alphaj - List_all_comb_b3_cent(1,i) += tmp_alphaj * nucl_coord(j,1) - List_all_comb_b3_cent(2,i) += tmp_alphaj * nucl_coord(j,2) - List_all_comb_b3_cent(3,i) += tmp_alphaj * nucl_coord(j,3) - - enddo - - ASSERT(List_all_comb_b3_expo(i) .gt. 0d0) - if(List_all_comb_b3_expo(i) .lt. 1d-10) cycle - - List_all_comb_b3_cent(1,i) = List_all_comb_b3_cent(1,i) / List_all_comb_b3_expo(i) - List_all_comb_b3_cent(2,i) = List_all_comb_b3_cent(2,i) / List_all_comb_b3_expo(i) - List_all_comb_b3_cent(3,i) = List_all_comb_b3_cent(3,i) / List_all_comb_b3_expo(i) - enddo - - ! --- - - do i = 1, List_all_comb_b3_size - - do j = 2, nucl_num, 1 - tmp_alphaj = dble(List_all_comb_b3(j,i)) * j1b_pen(j) - do k = 1, j-1, 1 - tmp_alphak = dble(List_all_comb_b3(k,i)) * j1b_pen(k) - - List_all_comb_b3_coef(i) += tmp_alphaj * tmp_alphak * ( (nucl_coord(j,1) - nucl_coord(k,1)) * (nucl_coord(j,1) - nucl_coord(k,1)) & - + (nucl_coord(j,2) - nucl_coord(k,2)) * (nucl_coord(j,2) - nucl_coord(k,2)) & - + (nucl_coord(j,3) - nucl_coord(k,3)) * (nucl_coord(j,3) - nucl_coord(k,3)) ) - enddo - enddo - - if(List_all_comb_b3_expo(i) .lt. 1d-10) cycle - - List_all_comb_b3_coef(i) = List_all_comb_b3_coef(i) / List_all_comb_b3_expo(i) - enddo - - ! --- - - do i = 1, List_all_comb_b3_size - - facto = 1.d0 - phase = 0 - do j = 1, nucl_num - tmp_alphaj = dble(List_all_comb_b3(j,i)) - - facto *= 2.d0 / (gamma(tmp_alphaj+1.d0) * gamma(3.d0-tmp_alphaj)) - phase += List_all_comb_b3(j,i) - enddo - - List_all_comb_b3_coef(i) = (-1.d0)**dble(phase) * facto * dexp(-List_all_comb_b3_coef(i)) - enddo - -END_PROVIDER - -! --- - BEGIN_PROVIDER [ double precision, int2_grad1u_grad2u_j1b, (ao_num, ao_num, n_points_final_grid)] BEGIN_DOC 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 c789ad36..9dd715e2 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 @@ -1,110 +1,6 @@ ! --- -BEGIN_PROVIDER [ integer, List_all_comb_b2_size] - - implicit none - - List_all_comb_b2_size = 2**nucl_num - -END_PROVIDER - -! --- - -BEGIN_PROVIDER [ integer, List_all_comb_b2, (nucl_num, List_all_comb_b2_size)] - - implicit none - integer :: i, j - - if(nucl_num .gt. 32) then - print *, ' nucl_num = ', nucl_num, '> 32' - stop - endif - - List_all_comb_b2 = 0 - - do i = 0, List_all_comb_b2_size-1 - do j = 0, nucl_num-1 - if (btest(i,j)) then - List_all_comb_b2(j+1,i+1) = 1 - endif - enddo - enddo - -END_PROVIDER - -! --- - - BEGIN_PROVIDER [ double precision, List_all_comb_b2_coef, ( List_all_comb_b2_size)] -&BEGIN_PROVIDER [ double precision, List_all_comb_b2_expo, ( List_all_comb_b2_size)] -&BEGIN_PROVIDER [ double precision, List_all_comb_b2_cent, (3, List_all_comb_b2_size)] - - implicit none - integer :: i, j, k, phase - double precision :: tmp_alphaj, tmp_alphak - - provide j1b_pen - - List_all_comb_b2_coef = 0.d0 - List_all_comb_b2_expo = 0.d0 - List_all_comb_b2_cent = 0.d0 - - do i = 1, List_all_comb_b2_size - - do j = 1, nucl_num - tmp_alphaj = dble(List_all_comb_b2(j,i)) * j1b_pen(j) - - List_all_comb_b2_expo(i) += tmp_alphaj - List_all_comb_b2_cent(1,i) += tmp_alphaj * nucl_coord(j,1) - List_all_comb_b2_cent(2,i) += tmp_alphaj * nucl_coord(j,2) - List_all_comb_b2_cent(3,i) += tmp_alphaj * nucl_coord(j,3) - - enddo - - ASSERT(List_all_comb_b2_expo(i) .gt. 0d0) - if(List_all_comb_b2_expo(i) .lt. 1d-10) cycle - - List_all_comb_b2_cent(1,i) = List_all_comb_b2_cent(1,i) / List_all_comb_b2_expo(i) - List_all_comb_b2_cent(2,i) = List_all_comb_b2_cent(2,i) / List_all_comb_b2_expo(i) - List_all_comb_b2_cent(3,i) = List_all_comb_b2_cent(3,i) / List_all_comb_b2_expo(i) - enddo - - ! --- - - do i = 1, List_all_comb_b2_size - - do j = 2, nucl_num, 1 - tmp_alphaj = dble(List_all_comb_b2(j,i)) * j1b_pen(j) - do k = 1, j-1, 1 - tmp_alphak = dble(List_all_comb_b2(k,i)) * j1b_pen(k) - - List_all_comb_b2_coef(i) += tmp_alphaj * tmp_alphak * ( (nucl_coord(j,1) - nucl_coord(k,1)) * (nucl_coord(j,1) - nucl_coord(k,1)) & - + (nucl_coord(j,2) - nucl_coord(k,2)) * (nucl_coord(j,2) - nucl_coord(k,2)) & - + (nucl_coord(j,3) - nucl_coord(k,3)) * (nucl_coord(j,3) - nucl_coord(k,3)) ) - enddo - enddo - - if(List_all_comb_b2_expo(i) .lt. 1d-10) cycle - - List_all_comb_b2_coef(i) = List_all_comb_b2_coef(i) / List_all_comb_b2_expo(i) - enddo - - ! --- - - do i = 1, List_all_comb_b2_size - - phase = 0 - do j = 1, nucl_num - phase += List_all_comb_b2(j,i) - enddo - - List_all_comb_b2_coef(i) = (-1.d0)**dble(phase) * dexp(-List_all_comb_b2_coef(i)) - enddo - -END_PROVIDER - -! --- - BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_j1b, (ao_num, ao_num, n_points_final_grid)] BEGIN_DOC @@ -157,7 +53,6 @@ BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_j1b, (ao_num, ao_num, n_po tmp(j,i,ipoint) += coef * (int_mu - int_coulomb) enddo - enddo enddo enddo diff --git a/src/ao_many_one_e_ints/listj1b_product_to_sum.irp.f b/src/ao_many_one_e_ints/listj1b_product_to_sum.irp.f new file mode 100644 index 00000000..ff9f8ae5 --- /dev/null +++ b/src/ao_many_one_e_ints/listj1b_product_to_sum.irp.f @@ -0,0 +1,227 @@ + +! --- + +BEGIN_PROVIDER [ integer, List_all_comb_b2_size] + + implicit none + + List_all_comb_b2_size = 2**nucl_num + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [ integer, List_all_comb_b2, (nucl_num, List_all_comb_b2_size)] + + implicit none + integer :: i, j + + if(nucl_num .gt. 32) then + print *, ' nucl_num = ', nucl_num, '> 32' + stop + endif + + List_all_comb_b2 = 0 + + do i = 0, List_all_comb_b2_size-1 + do j = 0, nucl_num-1 + if (btest(i,j)) then + List_all_comb_b2(j+1,i+1) = 1 + endif + enddo + enddo + +END_PROVIDER + +! --- + + BEGIN_PROVIDER [ double precision, List_all_comb_b2_coef, ( List_all_comb_b2_size)] +&BEGIN_PROVIDER [ double precision, List_all_comb_b2_expo, ( List_all_comb_b2_size)] +&BEGIN_PROVIDER [ double precision, List_all_comb_b2_cent, (3, List_all_comb_b2_size)] + + implicit none + integer :: i, j, k, phase + double precision :: tmp_alphaj, tmp_alphak + double precision :: tmp_cent_x, tmp_cent_y, tmp_cent_z + + provide j1b_pen + + List_all_comb_b2_coef = 0.d0 + List_all_comb_b2_expo = 0.d0 + List_all_comb_b2_cent = 0.d0 + + do i = 1, List_all_comb_b2_size + + tmp_cent_x = 0.d0 + tmp_cent_y = 0.d0 + tmp_cent_z = 0.d0 + do j = 1, nucl_num + tmp_alphaj = dble(List_all_comb_b2(j,i)) * j1b_pen(j) + List_all_comb_b2_expo(i) += tmp_alphaj + tmp_cent_x += tmp_alphaj * nucl_coord(j,1) + tmp_cent_y += tmp_alphaj * nucl_coord(j,2) + tmp_cent_z += tmp_alphaj * nucl_coord(j,3) + enddo + + ASSERT(List_all_comb_b2_expo(i) .gt. 0d0) + if(List_all_comb_b2_expo(i) .lt. 1d-10) cycle + + List_all_comb_b2_cent(1,i) = tmp_cent_x / List_all_comb_b2_expo(i) + List_all_comb_b2_cent(2,i) = tmp_cent_y / List_all_comb_b2_expo(i) + List_all_comb_b2_cent(3,i) = tmp_cent_z / List_all_comb_b2_expo(i) + enddo + + ! --- + + do i = 1, List_all_comb_b2_size + + do j = 2, nucl_num, 1 + tmp_alphaj = dble(List_all_comb_b2(j,i)) * j1b_pen(j) + do k = 1, j-1, 1 + tmp_alphak = dble(List_all_comb_b2(k,i)) * j1b_pen(k) + + List_all_comb_b2_coef(i) += tmp_alphaj * tmp_alphak * ( (nucl_coord(j,1) - nucl_coord(k,1)) * (nucl_coord(j,1) - nucl_coord(k,1)) & + + (nucl_coord(j,2) - nucl_coord(k,2)) * (nucl_coord(j,2) - nucl_coord(k,2)) & + + (nucl_coord(j,3) - nucl_coord(k,3)) * (nucl_coord(j,3) - nucl_coord(k,3)) ) + enddo + enddo + + if(List_all_comb_b2_expo(i) .lt. 1d-10) cycle + + List_all_comb_b2_coef(i) = List_all_comb_b2_coef(i) / List_all_comb_b2_expo(i) + enddo + + ! --- + + do i = 1, List_all_comb_b2_size + + phase = 0 + do j = 1, nucl_num + phase += List_all_comb_b2(j,i) + enddo + + List_all_comb_b2_coef(i) = (-1.d0)**dble(phase) * dexp(-List_all_comb_b2_coef(i)) + enddo + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [ integer, List_all_comb_b3_size] + + implicit none + + List_all_comb_b3_size = 3**nucl_num + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [ integer, List_all_comb_b3, (nucl_num, List_all_comb_b3_size)] + + implicit none + integer :: i, j, ii, jj + integer, allocatable :: M(:,:), p(:) + + if(nucl_num .gt. 32) then + print *, ' nucl_num = ', nucl_num, '> 32' + stop + endif + + List_all_comb_b3(:,:) = 0 + List_all_comb_b3(:,List_all_comb_b3_size) = 2 + + allocate(p(nucl_num)) + p = 0 + + do i = 2, List_all_comb_b3_size-1 + do j = 1, nucl_num + + ii = 0 + do jj = 1, j-1, 1 + ii = ii + p(jj) * 3**(jj-1) + enddo + p(j) = modulo(i-1-ii, 3**j) / 3**(j-1) + + List_all_comb_b3(j,i) = p(j) + enddo + enddo + +END_PROVIDER + +! --- + + BEGIN_PROVIDER [ double precision, List_all_comb_b3_coef, ( List_all_comb_b3_size)] +&BEGIN_PROVIDER [ double precision, List_all_comb_b3_expo, ( List_all_comb_b3_size)] +&BEGIN_PROVIDER [ double precision, List_all_comb_b3_cent, (3, List_all_comb_b3_size)] + + implicit none + integer :: i, j, k, phase + double precision :: tmp_alphaj, tmp_alphak, facto + + provide j1b_pen + + List_all_comb_b3_coef = 0.d0 + List_all_comb_b3_expo = 0.d0 + List_all_comb_b3_cent = 0.d0 + + do i = 1, List_all_comb_b3_size + + do j = 1, nucl_num + tmp_alphaj = dble(List_all_comb_b3(j,i)) * j1b_pen(j) + + List_all_comb_b3_expo(i) += tmp_alphaj + List_all_comb_b3_cent(1,i) += tmp_alphaj * nucl_coord(j,1) + List_all_comb_b3_cent(2,i) += tmp_alphaj * nucl_coord(j,2) + List_all_comb_b3_cent(3,i) += tmp_alphaj * nucl_coord(j,3) + + enddo + + ASSERT(List_all_comb_b3_expo(i) .gt. 0d0) + if(List_all_comb_b3_expo(i) .lt. 1d-10) cycle + + List_all_comb_b3_cent(1,i) = List_all_comb_b3_cent(1,i) / List_all_comb_b3_expo(i) + List_all_comb_b3_cent(2,i) = List_all_comb_b3_cent(2,i) / List_all_comb_b3_expo(i) + List_all_comb_b3_cent(3,i) = List_all_comb_b3_cent(3,i) / List_all_comb_b3_expo(i) + enddo + + ! --- + + do i = 1, List_all_comb_b3_size + + do j = 2, nucl_num, 1 + tmp_alphaj = dble(List_all_comb_b3(j,i)) * j1b_pen(j) + do k = 1, j-1, 1 + tmp_alphak = dble(List_all_comb_b3(k,i)) * j1b_pen(k) + + List_all_comb_b3_coef(i) += tmp_alphaj * tmp_alphak * ( (nucl_coord(j,1) - nucl_coord(k,1)) * (nucl_coord(j,1) - nucl_coord(k,1)) & + + (nucl_coord(j,2) - nucl_coord(k,2)) * (nucl_coord(j,2) - nucl_coord(k,2)) & + + (nucl_coord(j,3) - nucl_coord(k,3)) * (nucl_coord(j,3) - nucl_coord(k,3)) ) + enddo + enddo + + if(List_all_comb_b3_expo(i) .lt. 1d-10) cycle + + List_all_comb_b3_coef(i) = List_all_comb_b3_coef(i) / List_all_comb_b3_expo(i) + enddo + + ! --- + + do i = 1, List_all_comb_b3_size + + facto = 1.d0 + phase = 0 + do j = 1, nucl_num + tmp_alphaj = dble(List_all_comb_b3(j,i)) + + facto *= 2.d0 / (gamma(tmp_alphaj+1.d0) * gamma(3.d0-tmp_alphaj)) + phase += List_all_comb_b3(j,i) + enddo + + List_all_comb_b3_coef(i) = (-1.d0)**dble(phase) * facto * dexp(-List_all_comb_b3_coef(i)) + enddo + +END_PROVIDER + +! --- + 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 d4ef4b28..e0d6254a 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 @@ -82,7 +82,7 @@ double precision function NAI_pol_mult_erf_ao_with1s(i_ao, j_ao, beta, B_center, double precision, external :: NAI_pol_mult_erf_with1s, NAI_pol_mult_erf_ao - ASSERT(beta .lt. 0.d0) + ASSERT(beta .ge. 0.d0) if(beta .lt. 1d-10) then NAI_pol_mult_erf_ao_with1s = NAI_pol_mult_erf_ao(i_ao, j_ao, mu_in, C_center) return @@ -234,17 +234,17 @@ double precision function NAI_pol_mult_erf_with1s( A1_center, A2_center, power_A double precision :: rint - ! e^{-alpha1 (r - A1)^2} e^{-alpha2 (r - A2)^2} = e^{K12} e^{-alpha12 (r - A12)^2} + ! e^{-alpha1 (r - A1)^2} e^{-alpha2 (r - A2)^2} = e^{-K12} e^{-alpha12 (r - A12)^2} alpha12 = alpha1 + alpha2 alpha12_inv = 1.d0 / alpha12 alpha12_inv_2 = 0.5d0 * alpha12_inv rho12 = alpha1 * alpha2 * alpha12_inv - - dist12 = 0.d0 - do i = 1, 3 - A12_center(i) = (alpha1 * A1_center(i) + alpha2 * A2_center(i)) * alpha12_inv - dist12 += (A1_center(i) - A2_center(i)) * (A1_center(i) - A2_center(i)) - enddo + A12_center(1) = (alpha1 * A1_center(1) + alpha2 * A2_center(1)) * alpha12_inv + A12_center(2) = (alpha1 * A1_center(2) + alpha2 * A2_center(2)) * alpha12_inv + A12_center(3) = (alpha1 * A1_center(3) + alpha2 * A2_center(3)) * alpha12_inv + dist12 = ( (A1_center(1) - A2_center(1)) * (A1_center(1) - A2_center(1)) & + + (A1_center(2) - A2_center(2)) * (A1_center(2) - A2_center(2)) & + + (A1_center(3) - A2_center(3)) * (A1_center(3) - A2_center(3)) ) const_factor12 = dist12 * rho12 if(const_factor12 > 80.d0) then @@ -254,19 +254,17 @@ double precision function NAI_pol_mult_erf_with1s( A1_center, A2_center, power_A ! --- - ! e^{K12} e^{-alpha12 (r - A12)^2} e^{-beta (r - B)^2} = e^{K} e^{-p (r - P)^2} - p = alpha12 + beta - p_inv = 1.d0 / p - p_inv_2 = 0.5d0 * p_inv - rho = alpha12 * beta * p_inv - - dist = 0.d0 - dist_integral = 0.d0 - do i = 1, 3 - P_center(i) = (alpha12 * A12_center(i) + beta * B_center(i)) * p_inv - dist += (A12_center(i) - B_center(i)) * (A12_center(i) - B_center(i)) - dist_integral += (P_center(i) - C_center(i)) * (P_center(i) - C_center(i)) - enddo + ! e^{-K12} e^{-alpha12 (r - A12)^2} e^{-beta (r - B)^2} = e^{-K} e^{-p (r - P)^2} + p = alpha12 + beta + p_inv = 1.d0 / p + p_inv_2 = 0.5d0 * p_inv + rho = alpha12 * beta * p_inv + P_center(1) = (alpha12 * A12_center(1) + beta * B_center(1)) * p_inv + P_center(2) = (alpha12 * A12_center(2) + beta * B_center(2)) * p_inv + P_center(3) = (alpha12 * A12_center(3) + beta * B_center(3)) * p_inv + dist = ( (A12_center(1) - B_center(1)) * (A12_center(1) - B_center(1)) & + + (A12_center(2) - B_center(2)) * (A12_center(2) - B_center(2)) & + + (A12_center(3) - B_center(3)) * (A12_center(3) - B_center(3)) ) const_factor = const_factor12 + dist * rho if(const_factor > 80.d0) then @@ -274,6 +272,12 @@ double precision function NAI_pol_mult_erf_with1s( A1_center, A2_center, power_A return endif + dist_integral = 0.d0 + do i = 1, 3 + dist_integral += (P_center(i) - C_center(i)) * (P_center(i) - C_center(i)) + enddo + + ! --- p_new = mu_in / dsqrt(p + mu_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 54026349..2bcd1358 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 @@ -6,15 +6,24 @@ program debug_integ_jmu_modif implicit none my_grid_becke = .True. - !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 = 30 + my_n_pt_a_grid = 50 + !my_n_pt_r_grid = 100 + !my_n_pt_a_grid = 170 touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid PROVIDE mu_erf j1b_pen - call test_grad_1_u_ij_mu() + call test_j1b_nucl() + call test_grad_j1b_nucl() + call test_lapl_j1b_nucl() + + call test_list_b2() + call test_list_b3() + + !call test_grad_1_u_ij_mu() + !call test_gradu_squared_u_ij_mu() end @@ -23,16 +32,13 @@ end subroutine test_grad_1_u_ij_mu() implicit none - integer :: i, j, ipoint - double precision :: acc_ij, acc_tot, eps_ij, i_exc, i_num - double precision, external :: num_grad_1_u_ij_mu_x - double precision, external :: num_grad_1_u_ij_mu_y - double precision, external :: num_grad_1_u_ij_mu_z + integer :: i, j, ipoint + double precision :: acc_ij, acc_tot, eps_ij, i_exc, i_num, normalz + double precision :: integ(3) print*, ' test_grad_1_u_ij_mu ...' - PROVIDE grad_1_u_ij_mu -! PROVIDE num_grad_1_u_ij_mu + PROVIDE grad_1_u_ij_mu eps_ij = 1d-6 acc_tot = 0.d0 @@ -41,9 +47,10 @@ subroutine test_grad_1_u_ij_mu() do j = 1, ao_num do i = 1, ao_num - i_exc = grad_1_u_ij_mu(i,j,ipoint,1) - !i_num = num_grad_1_u_ij_mu(i,j,ipoint,1) - i_num = num_grad_1_u_ij_mu_x(i, j, ipoint) + call num_grad_1_u_ij_mu(i, j, ipoint, integ) + + i_exc = grad_1_u_ij_mu(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 grad_1_u_ij_mu on', i, j, ipoint @@ -52,10 +59,10 @@ subroutine test_grad_1_u_ij_mu() print *, ' diff = ', acc_ij endif acc_tot += acc_ij + normalz += dabs(i_num) - i_exc = grad_1_u_ij_mu(i,j,ipoint,2) - !i_num = num_grad_1_u_ij_mu(i,j,ipoint,2) - i_num = num_grad_1_u_ij_mu_y(i, j, ipoint) + i_exc = grad_1_u_ij_mu(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 grad_1_u_ij_mu on', i, j, ipoint @@ -64,30 +71,316 @@ subroutine test_grad_1_u_ij_mu() print *, ' diff = ', acc_ij endif acc_tot += acc_ij + normalz += dabs(i_num) - i_exc = grad_1_u_ij_mu(i,j,ipoint,3) - !i_num = num_grad_1_u_ij_mu(i,j,ipoint,3) - i_num = num_grad_1_u_ij_mu_z(i, j, ipoint) + i_exc = grad_1_u_ij_mu(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 y part of grad_1_u_ij_mu on', i, j, ipoint + print *, ' problem in z part of grad_1_u_ij_mu 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 / dble(ao_num*ao_num*n_points_final_grid) + acc_tot = acc_tot / normalz print*, ' normalized acc = ', acc_tot + print*, ' normalz = ', normalz return end subroutine test_grad_1_u_ij_mu ! --- +subroutine test_gradu_squared_u_ij_mu() + implicit none + integer :: i, j, ipoint + double precision :: acc_ij, acc_tot, eps_ij, i_exc, i_num, normalz + double precision, external :: num_gradu_squared_u_ij_mu + + print*, ' test_gradu_squared_u_ij_mu ...' + + PROVIDE gradu_squared_u_ij_mu + + eps_ij = 1d-6 + acc_tot = 0.d0 + + do ipoint = 1, n_points_final_grid + do j = 1, ao_num + do i = 1, ao_num + + i_exc = gradu_squared_u_ij_mu(i,j,ipoint) + i_num = num_gradu_squared_u_ij_mu(i, j, ipoint) + acc_ij = dabs(i_exc - i_num) + if(acc_ij .gt. eps_ij) then + print *, ' problem in gradu_squared_u_ij_mu 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_gradu_squared_u_ij_mu + +! --- + +subroutine test_j1b_nucl() + + implicit none + integer :: ipoint + double precision :: acc_ij, acc_tot, eps_ij, i_exc, i_num, normalz + double precision :: r(3) + double precision, external :: j1b_nucl + + print*, ' test_j1b_nucl ...' + + PROVIDE v_1b + + eps_ij = 1d-7 + acc_tot = 0.d0 + + 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) + + i_exc = v_1b(ipoint) + i_num = j1b_nucl(r) + acc_ij = dabs(i_exc - i_num) + if(acc_ij .gt. eps_ij) then + print *, ' problem in v_1b on', ipoint + print *, ' analyt = ', i_exc + print *, ' numeri = ', i_num + print *, ' diff = ', acc_ij + endif + + acc_tot += acc_ij + normalz += dabs(i_num) + enddo + + acc_tot = acc_tot / normalz + print*, ' normalized acc = ', acc_tot + print*, ' normalz = ', normalz + + return +end subroutine test_j1b_nucl + +! --- + +subroutine test_grad_j1b_nucl() + + implicit none + integer :: ipoint + double precision :: acc_ij, acc_tot, eps_ij, i_exc, i_num, normalz + double precision :: r(3) + double precision, external :: grad_x_j1b_nucl + double precision, external :: grad_y_j1b_nucl + double precision, external :: grad_z_j1b_nucl + + print*, ' test_grad_j1b_nucl ...' + + PROVIDE v_1b_grad + + eps_ij = 1d-6 + acc_tot = 0.d0 + + 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) + + i_exc = v_1b_grad(1,ipoint) + i_num = grad_x_j1b_nucl(r) + acc_ij = dabs(i_exc - i_num) + if(acc_ij .gt. eps_ij) then + print *, ' problem in x of v_1b_grad on', ipoint + print *, ' analyt = ', i_exc + print *, ' numeri = ', i_num + print *, ' diff = ', acc_ij + endif + + i_exc = v_1b_grad(2,ipoint) + i_num = grad_y_j1b_nucl(r) + acc_ij = dabs(i_exc - i_num) + if(acc_ij .gt. eps_ij) then + print *, ' problem in y of v_1b_grad on', ipoint + print *, ' analyt = ', i_exc + print *, ' numeri = ', i_num + print *, ' diff = ', acc_ij + endif + + i_exc = v_1b_grad(3,ipoint) + i_num = grad_z_j1b_nucl(r) + acc_ij = dabs(i_exc - i_num) + if(acc_ij .gt. eps_ij) then + print *, ' problem in z of v_1b_grad on', ipoint + print *, ' analyt = ', i_exc + print *, ' numeri = ', i_num + print *, ' diff = ', acc_ij + endif + + acc_tot += acc_ij + normalz += dabs(i_num) + enddo + + acc_tot = acc_tot / normalz + print*, ' normalized acc = ', acc_tot + print*, ' normalz = ', normalz + + return +end subroutine test_grad_j1b_nucl + +! --- + +subroutine test_lapl_j1b_nucl() + + implicit none + integer :: ipoint + double precision :: acc_ij, acc_tot, eps_ij, i_exc, i_num, normalz + double precision :: r(3) + double precision, external :: lapl_j1b_nucl + + print*, ' test_lapl_j1b_nucl ...' + + PROVIDE v_1b_lapl + + eps_ij = 1d-5 + acc_tot = 0.d0 + + 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) + + i_exc = v_1b_lapl(ipoint) + i_num = lapl_j1b_nucl(r) + acc_ij = dabs(i_exc - i_num) + if(acc_ij .gt. eps_ij) then + print *, ' problem in v_1b_lapl on', ipoint + print *, ' analyt = ', i_exc + print *, ' numeri = ', i_num + print *, ' diff = ', acc_ij + endif + + acc_tot += acc_ij + normalz += dabs(i_num) + enddo + + acc_tot = acc_tot / normalz + print*, ' normalized acc = ', acc_tot + print*, ' normalz = ', normalz + + return +end subroutine test_lapl_j1b_nucl + +! --- + +subroutine test_list_b2() + + implicit none + integer :: ipoint + double precision :: acc_ij, acc_tot, eps_ij, i_exc, i_num, normalz + double precision :: r(3) + double precision, external :: j1b_nucl + + print*, ' test_list_b2 ...' + + PROVIDE v_1b_list_b2 + + eps_ij = 1d-7 + acc_tot = 0.d0 + + 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) + + i_exc = v_1b_list_b2(ipoint) + i_num = j1b_nucl(r) + acc_ij = dabs(i_exc - i_num) + if(acc_ij .gt. eps_ij) then + print *, ' problem in list_b2 on', ipoint + print *, ' analyt = ', i_exc + print *, ' numeri = ', i_num + print *, ' diff = ', acc_ij + endif + + acc_tot += acc_ij + normalz += dabs(i_num) + enddo + + acc_tot = acc_tot / normalz + print*, ' normalized acc = ', acc_tot + print*, ' normalz = ', normalz + + return +end subroutine test_list_b2 + +! --- + +subroutine test_list_b3() + + implicit none + integer :: ipoint + double precision :: acc_ij, acc_tot, eps_ij, i_exc, i_tmp, i_num, normalz + double precision :: r(3) + double precision, external :: j1b_nucl + + print*, ' test_list_b3 ...' + + PROVIDE v_1b_list_b3 + + eps_ij = 1d-7 + acc_tot = 0.d0 + + 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) + + i_exc = v_1b_list_b3(ipoint) + i_tmp = j1b_nucl(r) + i_num = i_tmp * i_tmp + acc_ij = dabs(i_exc - i_num) + if(acc_ij .gt. eps_ij) then + print *, ' problem in list_b3 on', ipoint + print *, ' analyt = ', i_exc + print *, ' numeri = ', i_num + print *, ' diff = ', acc_ij + endif + + acc_tot += acc_ij + normalz += dabs(i_num) + enddo + + acc_tot = acc_tot / normalz + print*, ' normalized acc = ', acc_tot + print*, ' normalz = ', normalz + + return +end subroutine test_list_b3 + +! --- diff --git a/src/non_h_ints_mu/grad_squared.irp.f b/src/non_h_ints_mu/grad_squared.irp.f index 08152ddf..6ef34118 100644 --- a/src/non_h_ints_mu/grad_squared.irp.f +++ b/src/non_h_ints_mu/grad_squared.irp.f @@ -27,9 +27,10 @@ BEGIN_PROVIDER [ double precision, gradu_squared_u_ij_mu, (ao_num, ao_num,n_poin PROVIDE j1b_type j1b_pen if(j1b_type .eq. 3) then + ! v1_1b^2 \int d2 \phi_i(2) \phi_j(2) \frac{-[1 - \erf(\mu r12)]^2}{4} v2_1b^2 do ipoint = 1, n_points_final_grid - tmp = fact3_j12(ipoint) + 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) diff --git a/src/non_h_ints_mu/j12_nucl_utils.irp.f b/src/non_h_ints_mu/j12_nucl_utils.irp.f new file mode 100644 index 00000000..344b1fa8 --- /dev/null +++ b/src/non_h_ints_mu/j12_nucl_utils.irp.f @@ -0,0 +1,598 @@ + +! --- + +BEGIN_PROVIDER [ double precision, v_1b, (n_points_final_grid)] + + implicit none + integer :: ipoint, i, j, phase + double precision :: x, y, z, dx, dy, dz + double precision :: a, d, e, fact_r + + do ipoint = 1, n_points_final_grid + + x = final_grid_points(1,ipoint) + y = final_grid_points(2,ipoint) + z = final_grid_points(3,ipoint) + + fact_r = 1.d0 + do j = 1, nucl_num + a = j1b_pen(j) + dx = x - nucl_coord(j,1) + dy = y - nucl_coord(j,2) + dz = z - nucl_coord(j,3) + d = dx*dx + dy*dy + dz*dz + e = 1.d0 - dexp(-a*d) + + fact_r = fact_r * e + enddo + + v_1b(ipoint) = fact_r + enddo + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [ double precision, v_1b_grad, (3, n_points_final_grid)] + + implicit none + integer :: ipoint, i, j, phase + double precision :: x, y, z, dx, dy, dz + double precision :: a, d, e + double precision :: fact_x, fact_y, fact_z + double precision :: ax_der, ay_der, az_der, a_expo + + do ipoint = 1, n_points_final_grid + + x = final_grid_points(1,ipoint) + y = final_grid_points(2,ipoint) + z = final_grid_points(3,ipoint) + + fact_x = 0.d0 + fact_y = 0.d0 + fact_z = 0.d0 + do i = 1, List_all_comb_b2_size + + phase = 0 + a_expo = 0.d0 + ax_der = 0.d0 + ay_der = 0.d0 + az_der = 0.d0 + do j = 1, nucl_num + a = dble(List_all_comb_b2(j,i)) * j1b_pen(j) + dx = x - nucl_coord(j,1) + dy = y - nucl_coord(j,2) + dz = z - nucl_coord(j,3) + + phase += List_all_comb_b2(j,i) + a_expo += a * (dx*dx + dy*dy + dz*dz) + ax_der += a * dx + ay_der += a * dy + az_der += a * dz + enddo + e = -2.d0 * (-1.d0)**dble(phase) * dexp(-a_expo) + + fact_x += e * ax_der + fact_y += e * ay_der + fact_z += e * az_der + enddo + + v_1b_grad(1,ipoint) = fact_x + v_1b_grad(2,ipoint) = fact_y + v_1b_grad(3,ipoint) = fact_z + enddo + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [ double precision, v_1b_lapl, (n_points_final_grid)] + + implicit none + integer :: ipoint, i, j, phase + double precision :: x, y, z, dx, dy, dz + double precision :: a, d, e, b + double precision :: fact_r + double precision :: ax_der, ay_der, az_der, a_expo + + do ipoint = 1, n_points_final_grid + + x = final_grid_points(1,ipoint) + y = final_grid_points(2,ipoint) + z = final_grid_points(3,ipoint) + + fact_r = 0.d0 + do i = 1, List_all_comb_b2_size + + phase = 0 + b = 0.d0 + a_expo = 0.d0 + ax_der = 0.d0 + ay_der = 0.d0 + az_der = 0.d0 + do j = 1, nucl_num + a = dble(List_all_comb_b2(j,i)) * j1b_pen(j) + dx = x - nucl_coord(j,1) + dy = y - nucl_coord(j,2) + dz = z - nucl_coord(j,3) + + phase += List_all_comb_b2(j,i) + b += a + a_expo += a * (dx*dx + dy*dy + dz*dz) + ax_der += a * dx + ay_der += a * dy + az_der += a * dz + enddo + + fact_r += (-1.d0)**dble(phase) * (-6.d0*b + 4.d0*(ax_der*ax_der + ay_der*ay_der + az_der*az_der) ) * dexp(-a_expo) + enddo + + v_1b_lapl(ipoint) = fact_r + enddo + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [ double precision, v_1b_list_b2, (n_points_final_grid)] + + implicit none + integer :: i, ipoint + double precision :: x, y, z, coef, expo, dx, dy, dz + double precision :: fact_r + + PROVIDE List_all_comb_b2_coef List_all_comb_b2_expo List_all_comb_b2_cent + + do ipoint = 1, n_points_final_grid + + x = final_grid_points(1,ipoint) + y = final_grid_points(2,ipoint) + z = final_grid_points(3,ipoint) + + fact_r = 0.d0 + do i = 1, List_all_comb_b2_size + + coef = List_all_comb_b2_coef(i) + expo = List_all_comb_b2_expo(i) + + dx = x - List_all_comb_b2_cent(1,i) + dy = y - List_all_comb_b2_cent(2,i) + dz = z - List_all_comb_b2_cent(3,i) + + fact_r += coef * dexp(-expo * (dx*dx + dy*dy + dz*dz)) + enddo + + v_1b_list_b2(ipoint) = fact_r + enddo + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [ double precision, v_1b_list_b3, (n_points_final_grid)] + + implicit none + integer :: i, ipoint + double precision :: x, y, z, coef, expo, dx, dy, dz + double precision :: fact_r + + PROVIDE List_all_comb_b3_coef List_all_comb_b3_expo List_all_comb_b3_cent + + do ipoint = 1, n_points_final_grid + + x = final_grid_points(1,ipoint) + y = final_grid_points(2,ipoint) + z = final_grid_points(3,ipoint) + + fact_r = 0.d0 + do i = 1, List_all_comb_b3_size + + coef = List_all_comb_b3_coef(i) + expo = List_all_comb_b3_expo(i) + + dx = x - List_all_comb_b3_cent(1,i) + dy = y - List_all_comb_b3_cent(2,i) + dz = z - List_all_comb_b3_cent(3,i) + + fact_r += coef * dexp(-expo * (dx*dx + dy*dy + dz*dz)) + enddo + + v_1b_list_b3(ipoint) = fact_r + enddo + +END_PROVIDER + +! --- + +double precision function jmu_modif(r1, r2) + + implicit none + double precision, intent(in) :: r1(3), r2(3) + double precision, external :: j12_mu, j12_nucl + + jmu_modif = j12_mu(r1, r2) * j12_nucl(r1, r2) + + return +end function jmu_modif + +! --- + +double precision function j12_mu(r1, r2) + + include 'constants.include.F' + + implicit none + double precision, intent(in) :: r1(3), r2(3) + double precision :: mu_r12, r12 + + 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)) ) + mu_r12 = mu_erf * r12 + + j12_mu = 0.5d0 * r12 * (1.d0 - derf(mu_r12)) - inv_sq_pi_2 * dexp(-mu_r12*mu_r12) / mu_erf + + return +end function j12_mu + +! --- + +double precision function j1b_nucl(r) + + implicit none + double precision, intent(in) :: r(3) + integer :: i + double precision :: a, d, e + + j1b_nucl = 1.d0 + + do i = 1, nucl_num + a = j1b_pen(i) + d = ( (r(1) - nucl_coord(i,1)) * (r(1) - nucl_coord(i,1)) & + + (r(2) - nucl_coord(i,2)) * (r(2) - nucl_coord(i,2)) & + + (r(3) - nucl_coord(i,3)) * (r(3) - nucl_coord(i,3)) ) + e = 1.d0 - exp(-a*d) + + j1b_nucl = j1b_nucl * e + enddo + + return +end function j1b_nucl + +! --- + +double precision function j12_nucl(r1, r2) + + implicit none + double precision, intent(in) :: r1(3), r2(3) + double precision, external :: j1b_nucl + + j12_nucl = j1b_nucl(r1) * j1b_nucl(r2) + + return +end function j12_nucl + +! --- + +! --------------------------------------------------------------------------------------- + +double precision function grad_x_j1b_nucl(r) + + implicit none + double precision, intent(in) :: r(3) + double precision :: r_eps(3), eps, fp, fm, delta + double precision, external :: j1b_nucl + + eps = 1d-6 + r_eps = r + delta = max(eps, dabs(eps*r(1))) + + r_eps(1) = r_eps(1) + delta + fp = j1b_nucl(r_eps) + r_eps(1) = r_eps(1) - 2.d0 * delta + fm = j1b_nucl(r_eps) + + grad_x_j1b_nucl = 0.5d0 * (fp - fm) / delta + + return +end function grad_x_j1b_nucl + +double precision function grad_y_j1b_nucl(r) + + implicit none + double precision, intent(in) :: r(3) + double precision :: r_eps(3), eps, fp, fm, delta + double precision, external :: j1b_nucl + + eps = 1d-6 + r_eps = r + delta = max(eps, dabs(eps*r(2))) + + r_eps(2) = r_eps(2) + delta + fp = j1b_nucl(r_eps) + r_eps(2) = r_eps(2) - 2.d0 * delta + fm = j1b_nucl(r_eps) + + grad_y_j1b_nucl = 0.5d0 * (fp - fm) / delta + + return +end function grad_y_j1b_nucl + +double precision function grad_z_j1b_nucl(r) + + implicit none + double precision, intent(in) :: r(3) + double precision :: r_eps(3), eps, fp, fm, delta + double precision, external :: j1b_nucl + + eps = 1d-6 + r_eps = r + delta = max(eps, dabs(eps*r(3))) + + r_eps(3) = r_eps(3) + delta + fp = j1b_nucl(r_eps) + r_eps(3) = r_eps(3) - 2.d0 * delta + fm = j1b_nucl(r_eps) + + grad_z_j1b_nucl = 0.5d0 * (fp - fm) / delta + + return +end function grad_z_j1b_nucl + +! --------------------------------------------------------------------------------------- + +! --- + +double precision function lapl_j1b_nucl(r) + + implicit none + double precision, intent(in) :: r(3) + double precision :: r_eps(3), eps, fp, fm, delta + double precision, external :: grad_x_j1b_nucl + double precision, external :: grad_y_j1b_nucl + double precision, external :: grad_z_j1b_nucl + + eps = 1d-5 + r_eps = r + + lapl_j1b_nucl = 0.d0 + + ! --- + + delta = max(eps, dabs(eps*r(1))) + r_eps(1) = r_eps(1) + delta + fp = grad_x_j1b_nucl(r_eps) + r_eps(1) = r_eps(1) - 2.d0 * delta + fm = grad_x_j1b_nucl(r_eps) + r_eps(1) = r_eps(1) + delta + + lapl_j1b_nucl += 0.5d0 * (fp - fm) / delta + + ! --- + + delta = max(eps, dabs(eps*r(2))) + r_eps(2) = r_eps(2) + delta + fp = grad_y_j1b_nucl(r_eps) + r_eps(2) = r_eps(2) - 2.d0 * delta + fm = grad_y_j1b_nucl(r_eps) + r_eps(2) = r_eps(2) + delta + + lapl_j1b_nucl += 0.5d0 * (fp - fm) / delta + + ! --- + + delta = max(eps, dabs(eps*r(3))) + r_eps(3) = r_eps(3) + delta + fp = grad_z_j1b_nucl(r_eps) + r_eps(3) = r_eps(3) - 2.d0 * delta + fm = grad_z_j1b_nucl(r_eps) + r_eps(3) = r_eps(3) + delta + + lapl_j1b_nucl += 0.5d0 * (fp - fm) / delta + + ! --- + + return +end function lapl_j1b_nucl + +! --- + +! --------------------------------------------------------------------------------------- + +double precision function grad1_x_jmu_modif(r1, r2) + + implicit none + double precision, intent(in) :: r1(3), r2(3) + double precision :: r1_eps(3), eps, fp, fm, delta + double precision, external :: jmu_modif + + eps = 1d-7 + r1_eps = r1 + delta = max(eps, dabs(eps*r1(1))) + + r1_eps(1) = r1_eps(1) + delta + fp = jmu_modif(r1_eps, r2) + r1_eps(1) = r1_eps(1) - 2.d0 * delta + fm = jmu_modif(r1_eps, r2) + + grad1_x_jmu_modif = 0.5d0 * (fp - fm) / delta + + return +end function grad1_x_jmu_modif + +double precision function grad1_y_jmu_modif(r1, r2) + + implicit none + double precision, intent(in) :: r1(3), r2(3) + double precision :: r1_eps(3), eps, fp, fm, delta + double precision, external :: jmu_modif + + eps = 1d-7 + r1_eps = r1 + delta = max(eps, dabs(eps*r1(2))) + + r1_eps(2) = r1_eps(2) + delta + fp = jmu_modif(r1_eps, r2) + r1_eps(2) = r1_eps(2) - 2.d0 * delta + fm = jmu_modif(r1_eps, r2) + + grad1_y_jmu_modif = 0.5d0 * (fp - fm) / delta + + return +end function grad1_y_jmu_modif + +double precision function grad1_z_jmu_modif(r1, r2) + + implicit none + double precision, intent(in) :: r1(3), r2(3) + double precision :: r1_eps(3), eps, fp, fm, delta + double precision, external :: jmu_modif + + eps = 1d-7 + r1_eps = r1 + delta = max(eps, dabs(eps*r1(3))) + + r1_eps(3) = r1_eps(3) + delta + fp = jmu_modif(r1_eps, r2) + r1_eps(3) = r1_eps(3) - 2.d0 * delta + fm = jmu_modif(r1_eps, r2) + + grad1_z_jmu_modif = 0.5d0 * (fp - fm) / delta + + return +end function grad1_z_jmu_modif + +! --------------------------------------------------------------------------------------- + +! --- + +! --------------------------------------------------------------------------------------- + +double precision function grad1_x_j12_mu_num(r1, r2) + + implicit none + double precision, intent(in) :: r1(3), r2(3) + double precision :: r1_eps(3), eps, fp, fm, delta + double precision, external :: j12_mu + + eps = 1d-7 + r1_eps = r1 + delta = max(eps, dabs(eps*r1(1))) + + r1_eps(1) = r1_eps(1) + delta + fp = j12_mu(r1_eps, r2) + r1_eps(1) = r1_eps(1) - 2.d0 * delta + fm = j12_mu(r1_eps, r2) + + grad1_x_j12_mu_num = 0.5d0 * (fp - fm) / delta + + return +end function grad1_x_j12_mu_num + +double precision function grad1_y_j12_mu_num(r1, r2) + + implicit none + double precision, intent(in) :: r1(3), r2(3) + double precision :: r1_eps(3), eps, fp, fm, delta + double precision, external :: j12_mu + + eps = 1d-7 + r1_eps = r1 + delta = max(eps, dabs(eps*r1(2))) + + r1_eps(2) = r1_eps(2) + delta + fp = j12_mu(r1_eps, r2) + r1_eps(2) = r1_eps(2) - 2.d0 * delta + fm = j12_mu(r1_eps, r2) + + grad1_y_j12_mu_num = 0.5d0 * (fp - fm) / delta + + return +end function grad1_y_j12_mu_num + +double precision function grad1_z_j12_mu_num(r1, r2) + + implicit none + double precision, intent(in) :: r1(3), r2(3) + double precision :: r1_eps(3), eps, fp, fm, delta + double precision, external :: j12_mu + + eps = 1d-7 + r1_eps = r1 + delta = max(eps, dabs(eps*r1(3))) + + r1_eps(3) = r1_eps(3) + delta + fp = j12_mu(r1_eps, r2) + r1_eps(3) = r1_eps(3) - 2.d0 * delta + fm = j12_mu(r1_eps, r2) + + grad1_z_j12_mu_num = 0.5d0 * (fp - fm) / delta + + return +end function grad1_z_j12_mu_num + +! --------------------------------------------------------------------------------------- + +! --- + +! --------------------------------------------------------------------------------------- + +double precision function grad1_x_j12_mu_exc(r1, r2) + + implicit none + double precision, intent(in) :: r1(3), r2(3) + double precision :: r12 + + grad1_x_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_x_j12_mu_exc = 0.5d0 * (1.d0 - derf(mu_erf * r12)) * (r1(1) - r2(1)) / r12 + + 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 + +! --------------------------------------------------------------------------------------- + +! --- + + diff --git a/src/non_h_ints_mu/jmu_modif.irp.f b/src/non_h_ints_mu/jmu_modif.irp.f deleted file mode 100644 index 59a4a104..00000000 --- a/src/non_h_ints_mu/jmu_modif.irp.f +++ /dev/null @@ -1,266 +0,0 @@ - -! --- - -double precision function jmu_modif(r1, r2) - - implicit none - double precision, intent(in) :: r1(3), r2(3) - double precision, external :: j12_mu, j12_nucl - - jmu_modif = j12_mu(r1, r2) * j12_nucl(r1, r2) - - return -end function jmu_modif - -! --- - -double precision function j12_mu(r1, r2) - - include 'constants.include.F' - - implicit none - double precision, intent(in) :: r1(3), r2(3) - double precision :: mu_r12, r12 - - 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)) ) - mu_r12 = mu_erf * r12 - - j12_mu = 0.5d0 * r12 * (1.d0 - derf(mu_r12)) - inv_sq_pi_2 * dexp(-mu_r12*mu_r12) / mu_erf - - return -end function j12_mu - -! --- - -double precision function j12_nucl(r1, r2) - - implicit none - double precision, intent(in) :: r1(3), r2(3) - integer :: i, j - double precision :: a1, d1, e1, a2, d2, e2 - - j12_nucl = 1.d0 - do i = 1, nucl_num - a1 = j1b_pen(i) - d1 = ( (r1(1) - nucl_coord(i,1)) * (r1(1) - nucl_coord(i,1)) & - + (r1(2) - nucl_coord(i,2)) * (r1(2) - nucl_coord(i,2)) & - + (r1(3) - nucl_coord(i,3)) * (r1(3) - nucl_coord(i,3)) ) - e1 = 1.d0 - exp(-a1*d1) - - do j = 1, nucl_num - a2 = j1b_pen(j) - d2 = ( (r2(1) - nucl_coord(j,1)) * (r2(1) - nucl_coord(j,1)) & - + (r2(2) - nucl_coord(j,2)) * (r2(2) - nucl_coord(j,2)) & - + (r2(3) - nucl_coord(j,3)) * (r2(3) - nucl_coord(j,3)) ) - e2 = 1.d0 - exp(-a2*d2) - - j12_nucl = j12_nucl * e1 * e2 - enddo - enddo - - return -end function j12_nucl - -! --- - -! --------------------------------------------------------------------------------------- - -double precision function grad1_x_jmu_modif(r1, r2) - - implicit none - double precision, intent(in) :: r1(3), r2(3) - double precision :: r1_eps(3), eps, fp, fm, delta - double precision, external :: jmu_modif - - eps = 1d-7 - r1_eps = r1 - delta = max(eps, dabs(eps*r1(1))) - - r1_eps(1) = r1_eps(1) + delta - fp = jmu_modif(r1_eps, r2) - r1_eps(1) = r1_eps(1) - 2.d0 * delta - fm = jmu_modif(r1_eps, r2) - - grad1_x_jmu_modif = 0.5d0 * (fp - fm) / delta - - return -end function grad1_x_jmu_modif - -double precision function grad1_y_jmu_modif(r1, r2) - - implicit none - double precision, intent(in) :: r1(3), r2(3) - double precision :: r1_eps(3), eps, fp, fm, delta - double precision, external :: jmu_modif - - eps = 1d-7 - r1_eps = r1 - delta = max(eps, dabs(eps*r1(2))) - - r1_eps(2) = r1_eps(2) + delta - fp = jmu_modif(r1_eps, r2) - r1_eps(2) = r1_eps(2) - 2.d0 * delta - fm = jmu_modif(r1_eps, r2) - - grad1_y_jmu_modif = 0.5d0 * (fp - fm) / delta - - return -end function grad1_y_jmu_modif - -double precision function grad1_z_jmu_modif(r1, r2) - - implicit none - double precision, intent(in) :: r1(3), r2(3) - double precision :: r1_eps(3), eps, fp, fm, delta - double precision, external :: jmu_modif - - eps = 1d-7 - r1_eps = r1 - delta = max(eps, dabs(eps*r1(3))) - - r1_eps(3) = r1_eps(3) + delta - fp = jmu_modif(r1_eps, r2) - r1_eps(3) = r1_eps(3) - 2.d0 * delta - fm = jmu_modif(r1_eps, r2) - - grad1_z_jmu_modif = 0.5d0 * (fp - fm) / delta - - return -end function grad1_z_jmu_modif - -! --------------------------------------------------------------------------------------- - -! --- - -! --------------------------------------------------------------------------------------- - -double precision function grad1_x_j12_mu_num(r1, r2) - - implicit none - double precision, intent(in) :: r1(3), r2(3) - double precision :: r1_eps(3), eps, fp, fm, delta - double precision, external :: j12_mu - - eps = 1d-7 - r1_eps = r1 - delta = max(eps, dabs(eps*r1(1))) - - r1_eps(1) = r1_eps(1) + delta - fp = j12_mu(r1_eps, r2) - r1_eps(1) = r1_eps(1) - 2.d0 * delta - fm = j12_mu(r1_eps, r2) - - grad1_x_j12_mu_num = 0.5d0 * (fp - fm) / delta - - return -end function grad1_x_j12_mu_num - -double precision function grad1_y_j12_mu_num(r1, r2) - - implicit none - double precision, intent(in) :: r1(3), r2(3) - double precision :: r1_eps(3), eps, fp, fm, delta - double precision, external :: j12_mu - - eps = 1d-7 - r1_eps = r1 - delta = max(eps, dabs(eps*r1(2))) - - r1_eps(2) = r1_eps(2) + delta - fp = j12_mu(r1_eps, r2) - r1_eps(2) = r1_eps(2) - 2.d0 * delta - fm = j12_mu(r1_eps, r2) - - grad1_y_j12_mu_num = 0.5d0 * (fp - fm) / delta - - return -end function grad1_y_j12_mu_num - -double precision function grad1_z_j12_mu_num(r1, r2) - - implicit none - double precision, intent(in) :: r1(3), r2(3) - double precision :: r1_eps(3), eps, fp, fm, delta - double precision, external :: j12_mu - - eps = 1d-7 - r1_eps = r1 - delta = max(eps, dabs(eps*r1(3))) - - r1_eps(3) = r1_eps(3) + delta - fp = j12_mu(r1_eps, r2) - r1_eps(3) = r1_eps(3) - 2.d0 * delta - fm = j12_mu(r1_eps, r2) - - grad1_z_j12_mu_num = 0.5d0 * (fp - fm) / delta - - return -end function grad1_z_j12_mu_num - -! --------------------------------------------------------------------------------------- - -! --- - -! --------------------------------------------------------------------------------------- - -double precision function grad1_x_j12_mu_exc(r1, r2) - - implicit none - double precision, intent(in) :: r1(3), r2(3) - double precision :: r12 - - grad1_x_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_x_j12_mu_exc = 0.5d0 * (1.d0 - derf(mu_erf * r12)) * (r1(1) - r2(1)) / r12 - - 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 - -! --------------------------------------------------------------------------------------- - -! --- - - diff --git a/src/non_h_ints_mu/new_grad_tc.irp.f b/src/non_h_ints_mu/new_grad_tc.irp.f index 26ed642c..9832e81d 100644 --- a/src/non_h_ints_mu/new_grad_tc.irp.f +++ b/src/non_h_ints_mu/new_grad_tc.irp.f @@ -1,81 +1,3 @@ - -! --- - - BEGIN_PROVIDER [ double precision, fact1_j12, ( n_points_final_grid)] -&BEGIN_PROVIDER [ double precision, fact2_j12, (3, n_points_final_grid)] -&BEGIN_PROVIDER [ double precision, fact3_j12, ( n_points_final_grid)] - - implicit none - integer :: ipoint, i, j, phase - double precision :: x, y, z, dx, dy, dz - double precision :: a, d, e, fact_r, fact_r_sq - double precision :: fact_x, fact_y, fact_z - double precision :: ax_der, ay_der, az_der, a_expo - - do ipoint = 1, n_points_final_grid - - x = final_grid_points(1,ipoint) - y = final_grid_points(2,ipoint) - z = final_grid_points(3,ipoint) - - ! --- - - fact_r = 1.d0 - fact_r_sq = 1.d0 - do j = 1, nucl_num - a = j1b_pen(j) - dx = x - nucl_coord(j,1) - dy = y - nucl_coord(j,2) - dz = z - nucl_coord(j,3) - d = x*x + y*y + z*z - e = 1.d0 - dexp(-a*d) - - fact_r = fact_r * e - fact_r_sq = fact_r_sq * e * e - enddo - fact1_j12(ipoint) = fact_r - fact3_j12(ipoint) = fact_r_sq - - ! --- - - fact_x = 0.d0 - fact_y = 0.d0 - fact_z = 0.d0 - do i = 1, List_all_comb_b2_size - - phase = 0 - a_expo = 0.d0 - ax_der = 0.d0 - ay_der = 0.d0 - az_der = 0.d0 - do j = 1, nucl_num - a = dble(List_all_comb_b2(j,i)) * j1b_pen(j) - dx = x - nucl_coord(j,1) - dy = y - nucl_coord(j,2) - dz = z - nucl_coord(j,3) - - phase += List_all_comb_b2(j,i) - a_expo += a * (dx*dx + dy*dy + dz*dz) - ax_der += a * dx - ay_der += a * dy - az_der += a * dz - enddo - e = -2.d0 * (-1.d0)**dble(phase) * dexp(-a_expo) - - fact_x += e * ax_der - fact_y += e * ay_der - fact_z += e * az_der - enddo - - fact2_j12(1,ipoint) = fact_x - fact2_j12(2,ipoint) = fact_y - fact2_j12(3,ipoint) = fact_z - - ! --- - - enddo - -END_PROVIDER ! --- @@ -103,10 +25,10 @@ BEGIN_PROVIDER [ double precision, grad_1_u_ij_mu, (ao_num, ao_num, n_points_fin y = final_grid_points(2,ipoint) z = final_grid_points(3,ipoint) - tmp0 = fact1_j12(ipoint) - tmp_x = fact2_j12(1,ipoint) - tmp_y = fact2_j12(2,ipoint) - tmp_z = fact2_j12(3,ipoint) + tmp0 = v_1b (ipoint) + tmp_x = v_1b_grad(1,ipoint) + tmp_y = v_1b_grad(2,ipoint) + tmp_z = v_1b_grad(3,ipoint) do j = 1, ao_num do i = 1, ao_num diff --git a/src/non_h_ints_mu/numerical_integ.irp.f b/src/non_h_ints_mu/numerical_integ.irp.f index abac1874..842851aa 100644 --- a/src/non_h_ints_mu/numerical_integ.irp.f +++ b/src/non_h_ints_mu/numerical_integ.irp.f @@ -4,132 +4,129 @@ ! ! \int dr2 [-1 * \grad_r1 u(r1,r2)] \phi_i(r2) \phi_j(r2) x 1s_j1b(r2) ! -BEGIN_PROVIDER [ double precision, num_grad_1_u_ij_mu, (ao_num, ao_num, n_points_final_grid, 3)] - implicit none - - integer :: i, j, ipoint, jpoint - double precision :: tmp, r1(3), r2(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 - - do j = 1, ao_num - do i = 1, ao_num - - 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) - - 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) - 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)) - enddo - - enddo - enddo - enddo - -END_PROVIDER +!BEGIN_PROVIDER [ double precision, num_grad_1_u_ij_mu, (ao_num, ao_num, n_points_final_grid, 3)] +! +! implicit none +! +! integer :: i, j, ipoint, jpoint +! double precision :: tmp, r1(3), r2(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 +! +! do j = 1, ao_num +! do i = 1, ao_num +! +! 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) +! +! 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) +! 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)) +! enddo +! +! enddo +! enddo +! enddo +! +!END_PROVIDER ! --- -double precision function num_grad_1_u_ij_mu_x(i, j, ipoint) +subroutine num_grad_1_u_ij_mu(i, j, ipoint, integ) implicit none - integer, intent(in) :: i, j, ipoint - integer :: jpoint - double precision :: tmp, r1(3), r2(3) - double precision, external :: ao_value - double precision, external :: j12_nucl - double precision, external :: grad1_x_j12_mu_num, grad1_x_j12_mu_exc - num_grad_1_u_ij_mu_x = 0.d0 + integer, intent(in) :: i, j, ipoint + double precision, intent(out) :: integ(3) + + integer :: jpoint + double precision :: tmp, r1(3), r2(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) 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) 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_x += tmp * (-1.d0 * grad1_x_j12_mu_exc(r1, r2)) + 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)) enddo -end function num_grad_1_u_ij_mu_x + integ(1) = tmp_x + integ(2) = tmp_y + integ(3) = tmp_z + + return +end subroutine num_grad_1_u_ij_mu ! --- -double precision function num_grad_1_u_ij_mu_y(i, j, ipoint) +double precision function num_gradu_squared_u_ij_mu(i, j, ipoint) implicit none + integer, intent(in) :: i, j, ipoint + integer :: jpoint - double precision :: tmp, r1(3), r2(3) + double precision :: tmp, r1(3), r2(3), r12 + double precision :: tmp_x, tmp_y, tmp_z, tmp1, tmp2 + double precision, external :: ao_value double precision, external :: j12_nucl - double precision, external :: grad1_y_j12_mu_num, grad1_y_j12_mu_exc - - num_grad_1_u_ij_mu_y = 0.d0 r1(1) = final_grid_points(1,ipoint) r1(2) = final_grid_points(2,ipoint) r1(3) = final_grid_points(3,ipoint) + num_gradu_squared_u_ij_mu = 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) - tmp = ao_value(i, r2) * ao_value(j, r2) * j12_nucl(r1, r2) * final_weight_at_r_vector(jpoint) + tmp_x = r1(1) - r2(1) + tmp_y = r1(2) - r2(2) + tmp_z = r1(3) - r2(3) + r12 = dsqrt( tmp_x*tmp_x + tmp_y*tmp_y + tmp_z*tmp_z ) + tmp1 = 1.d0 - derf(mu_erf * r12) + tmp2 = j12_nucl(r1, r2) + tmp = -0.25d0 * tmp1 * tmp1 * tmp2 * tmp2 * ao_value(i, r2) * ao_value(j, r2) * final_weight_at_r_vector(jpoint) - num_grad_1_u_ij_mu_y += tmp * (-1.d0 * grad1_y_j12_mu_exc(r1, r2)) + num_gradu_squared_u_ij_mu += tmp enddo -end function num_grad_1_u_ij_mu_y + return +end function num_gradu_squared_u_ij_mu ! --- -double precision function num_grad_1_u_ij_mu_z(i, j, ipoint) - - implicit none - integer, intent(in) :: i, j, ipoint - integer :: jpoint - double precision :: tmp, r1(3), r2(3) - double precision, external :: ao_value - double precision, external :: j12_nucl - double precision, external :: grad1_z_j12_mu_num, grad1_z_j12_mu_exc - - num_grad_1_u_ij_mu_z = 0.d0 - - r1(1) = final_grid_points(1,ipoint) - r1(2) = final_grid_points(2,ipoint) - r1(3) = final_grid_points(3,ipoint) - - 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) - 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_z += tmp * (-1.d0 * grad1_z_j12_mu_exc(r1, r2)) - enddo - -end function num_grad_1_u_ij_mu_z - -! ---