From f985af03953b2d25e6706cac2412eb38cf544075 Mon Sep 17 00:00:00 2001 From: Abdallah Ammar Date: Sun, 7 May 2023 15:07:51 +0200 Subject: [PATCH] jast 4 added --- src/ao_many_one_e_ints/NEED | 1 + src/ao_many_one_e_ints/grad2_jmu_modif.irp.f | 62 ++-- .../grad_lapl_jmu_modif.irp.f | 62 ++-- src/ao_many_one_e_ints/listj1b.irp.f | 296 +++++++++++++----- src/non_h_ints_mu/grad_squared.irp.f | 2 +- src/non_h_ints_mu/j12_nucl_utils.irp.f | 173 +++++++--- src/non_h_ints_mu/tc_integ.irp.f | 4 +- 7 files changed, 404 insertions(+), 196 deletions(-) diff --git a/src/ao_many_one_e_ints/NEED b/src/ao_many_one_e_ints/NEED index 0d08442c..c57219cd 100644 --- a/src/ao_many_one_e_ints/NEED +++ b/src/ao_many_one_e_ints/NEED @@ -3,3 +3,4 @@ ao_two_e_ints becke_numerical_grid mo_one_e_ints dft_utils_in_r +tc_keywords 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 722ff2ff..fb4d71f3 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 @@ -25,11 +25,11 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2, (ao_num, ao_num, n_poin int2_grad1u2_grad2u2 = 0.d0 - !$OMP PARALLEL DEFAULT (NONE) & - !$OMP PRIVATE (ipoint, i, j, i_fit, r, coef_fit, expo_fit, tmp) & - !$OMP SHARED (n_points_final_grid, ao_num, final_grid_points, ng_fit_jast, & - !$OMP expo_gauss_1_erf_x_2, coef_gauss_1_erf_x_2,int2_grad1u2_grad2u2) - !$OMP DO + !$OMP PARALLEL DEFAULT (NONE) & + !$OMP PRIVATE (ipoint, i, j, i_fit, r, coef_fit, expo_fit, tmp) & + !$OMP SHARED (n_points_final_grid, ao_num, final_grid_points, ng_fit_jast, & + !$OMP expo_gauss_1_erf_x_2, coef_gauss_1_erf_x_2,int2_grad1u2_grad2u2) + !$OMP DO do ipoint = 1, n_points_final_grid r(1) = final_grid_points(1,ipoint) r(2) = final_grid_points(2,ipoint) @@ -51,8 +51,8 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2, (ao_num, ao_num, n_poin enddo enddo enddo - !$OMP END DO - !$OMP END PARALLEL + !$OMP END DO + !$OMP END PARALLEL do ipoint = 1, n_points_final_grid do i = 2, ao_num @@ -94,15 +94,15 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2, (ao_num, ao_num, n int2_grad1u2_grad2u2_j1b2 = 0.d0 - !$OMP PARALLEL DEFAULT (NONE) & - !$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, & - !$OMP coef_fit, expo_fit, int_fit, tmp) & - !$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b3_size, & - !$OMP final_grid_points, ng_fit_jast, & - !$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_grad1u2_grad2u2_j1b2) - !$OMP DO + !$OMP PARALLEL DEFAULT (NONE) & + !$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, & + !$OMP coef_fit, expo_fit, int_fit, tmp) & + !$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b3_size, & + !$OMP final_grid_points, ng_fit_jast, & + !$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_grad1u2_grad2u2_j1b2) + !$OMP DO do ipoint = 1, n_points_final_grid r(1) = final_grid_points(1,ipoint) r(2) = final_grid_points(2,ipoint) @@ -121,7 +121,7 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2, (ao_num, ao_num, n int_fit = overlap_gauss_r12_ao(r, expo_fit, i, j) tmp += -0.25d0 * coef_fit * int_fit -! if(dabs(coef_fit*int_fit) .lt. 1d-12) cycle +! if(dabs(coef_fit*int_fit) .lt. 1d-12) cycle ! --- @@ -146,8 +146,8 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2, (ao_num, ao_num, n enddo enddo enddo - !$OMP END DO - !$OMP END PARALLEL + !$OMP END DO + !$OMP END PARALLEL do ipoint = 1, n_points_final_grid do i = 2, ao_num @@ -188,15 +188,15 @@ BEGIN_PROVIDER [ double precision, int2_u2_j1b2, (ao_num, ao_num, n_points_final int2_u2_j1b2 = 0.d0 - !$OMP PARALLEL DEFAULT (NONE) & - !$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, & - !$OMP coef_fit, expo_fit, int_fit, tmp) & - !$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b3_size, & - !$OMP final_grid_points, ng_fit_jast, & - !$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_j1b2) - !$OMP DO + !$OMP PARALLEL DEFAULT (NONE) & + !$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, & + !$OMP coef_fit, expo_fit, int_fit, tmp) & + !$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b3_size, & + !$OMP final_grid_points, ng_fit_jast, & + !$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_j1b2) + !$OMP DO do ipoint = 1, n_points_final_grid r(1) = final_grid_points(1,ipoint) r(2) = final_grid_points(2,ipoint) @@ -215,7 +215,7 @@ BEGIN_PROVIDER [ double precision, int2_u2_j1b2, (ao_num, ao_num, n_points_final int_fit = overlap_gauss_r12_ao(r, expo_fit, i, j) tmp += coef_fit * int_fit -! if(dabs(coef_fit*int_fit) .lt. 1d-12) cycle +! if(dabs(coef_fit*int_fit) .lt. 1d-12) cycle ! --- @@ -240,8 +240,8 @@ BEGIN_PROVIDER [ double precision, int2_u2_j1b2, (ao_num, ao_num, n_points_final enddo enddo enddo - !$OMP END DO - !$OMP END PARALLEL + !$OMP END DO + !$OMP END PARALLEL do ipoint = 1, n_points_final_grid do i = 2, ao_num 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 fc30cd83..25bb2f8b 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 @@ -24,12 +24,12 @@ BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_j1b, (ao_num, ao_num, n_po v_ij_erf_rk_cst_mu_j1b = 0.d0 - !$OMP PARALLEL DEFAULT (NONE) & - !$OMP PRIVATE (ipoint, i, j, i_1s, r, coef, beta, B_center, int_mu, int_coulomb, tmp) & - !$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b2_size, final_grid_points, & - !$OMP List_all_comb_b2_coef, List_all_comb_b2_expo, List_all_comb_b2_cent, & - !$OMP v_ij_erf_rk_cst_mu_j1b, mu_erf) - !$OMP DO + !$OMP PARALLEL DEFAULT (NONE) & + !$OMP PRIVATE (ipoint, i, j, i_1s, r, coef, beta, B_center, int_mu, int_coulomb, tmp) & + !$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b2_size, final_grid_points, & + !$OMP List_all_comb_b2_coef, List_all_comb_b2_expo, List_all_comb_b2_cent, & + !$OMP v_ij_erf_rk_cst_mu_j1b, mu_erf) + !$OMP DO !do ipoint = 1, 10 do ipoint = 1, n_points_final_grid r(1) = final_grid_points(1,ipoint) @@ -51,7 +51,7 @@ BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_j1b, (ao_num, ao_num, n_po int_mu = NAI_pol_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r) int_coulomb = NAI_pol_mult_erf_ao_with1s(i, j, beta, B_center, 1.d+9, r) -! if(dabs(coef)*dabs(int_mu - int_coulomb) .lt. 1d-12) cycle +! if(dabs(coef)*dabs(int_mu - int_coulomb) .lt. 1d-12) cycle tmp += coef * (int_mu - int_coulomb) @@ -77,8 +77,8 @@ BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_j1b, (ao_num, ao_num, n_po enddo enddo enddo - !$OMP END DO - !$OMP END PARALLEL + !$OMP END DO + !$OMP END PARALLEL do ipoint = 1, n_points_final_grid do i = 2, ao_num @@ -112,13 +112,13 @@ BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_j1b, (ao_num, ao_num, n_ x_v_ij_erf_rk_cst_mu_j1b = 0.d0 - !$OMP PARALLEL DEFAULT (NONE) & - !$OMP PRIVATE (ipoint, i, j, i_1s, r, coef, beta, B_center, ints, ints_coulomb, & - !$OMP tmp_x, tmp_y, tmp_z) & - !$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b2_size, final_grid_points,& - !$OMP List_all_comb_b2_coef, List_all_comb_b2_expo, List_all_comb_b2_cent, & - !$OMP x_v_ij_erf_rk_cst_mu_j1b, mu_erf) - !$OMP DO + !$OMP PARALLEL DEFAULT (NONE) & + !$OMP PRIVATE (ipoint, i, j, i_1s, r, coef, beta, B_center, ints, ints_coulomb, & + !$OMP tmp_x, tmp_y, tmp_z) & + !$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b2_size, final_grid_points,& + !$OMP List_all_comb_b2_coef, List_all_comb_b2_expo, List_all_comb_b2_cent, & + !$OMP x_v_ij_erf_rk_cst_mu_j1b, mu_erf) + !$OMP DO !do ipoint = 1, 10 do ipoint = 1, n_points_final_grid r(1) = final_grid_points(1,ipoint) @@ -143,7 +143,7 @@ BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_j1b, (ao_num, ao_num, n_ call NAI_pol_x_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r, ints ) call NAI_pol_x_mult_erf_ao_with1s(i, j, beta, B_center, 1.d+9, r, ints_coulomb) -! if( dabs(coef)*(dabs(ints(1)-ints_coulomb(1)) + dabs(ints(2)-ints_coulomb(2)) + dabs(ints(3)-ints_coulomb(3))) .lt. 3d-10) cycle +! if( dabs(coef)*(dabs(ints(1)-ints_coulomb(1)) + dabs(ints(2)-ints_coulomb(2)) + dabs(ints(3)-ints_coulomb(3))) .lt. 3d-10) cycle tmp_x += coef * (ints(1) - ints_coulomb(1)) tmp_y += coef * (ints(2) - ints_coulomb(2)) @@ -175,8 +175,8 @@ BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_j1b, (ao_num, ao_num, n_ enddo enddo enddo - !$OMP END DO - !$OMP END PARALLEL + !$OMP END DO + !$OMP END PARALLEL do ipoint = 1, n_points_final_grid do i = 2, ao_num @@ -220,15 +220,15 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b, (ao_num, ao_num, n_points_ v_ij_u_cst_mu_j1b = 0.d0 - !$OMP PARALLEL DEFAULT (NONE) & - !$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, & - !$OMP coef_fit, expo_fit, int_fit, tmp) & - !$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b2_size, & - !$OMP final_grid_points, ng_fit_jast, & - !$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) - !$OMP DO + !$OMP PARALLEL DEFAULT (NONE) & + !$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, & + !$OMP coef_fit, expo_fit, int_fit, tmp) & + !$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b2_size, & + !$OMP final_grid_points, ng_fit_jast, & + !$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) + !$OMP DO !do ipoint = 1, 10 do ipoint = 1, n_points_final_grid r(1) = final_grid_points(1,ipoint) @@ -253,7 +253,7 @@ 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 +! if(dabs(int_fit*coef) .lt. 1d-12) cycle tmp += coef * coef_fit * int_fit @@ -280,8 +280,8 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b, (ao_num, ao_num, n_points_ enddo enddo enddo - !$OMP END DO - !$OMP END PARALLEL + !$OMP END DO + !$OMP END PARALLEL do ipoint = 1, n_points_final_grid do i = 2, ao_num diff --git a/src/ao_many_one_e_ints/listj1b.irp.f b/src/ao_many_one_e_ints/listj1b.irp.f index 4698cb27..c06d64bb 100644 --- a/src/ao_many_one_e_ints/listj1b.irp.f +++ b/src/ao_many_one_e_ints/listj1b.irp.f @@ -1,17 +1,34 @@ ! --- -BEGIN_PROVIDER [ integer, List_all_comb_b2_size] +BEGIN_PROVIDER [integer, List_all_comb_b2_size] implicit none - List_all_comb_b2_size = 2**nucl_num + PROVIDE j1b_type + + if(j1b_type .eq. 3) then + + List_all_comb_b2_size = 2**nucl_num + + elseif(j1b_type .eq. 4) then + + List_all_comb_b2_size = nucl_num + 1 + + else + + print *, 'j1b_type = ', j1b_pen, 'is not implemented' + stop + + endif + + print *, ' nb of linear terms in the envelope is ', List_all_comb_b2_size END_PROVIDER ! --- -BEGIN_PROVIDER [ integer, List_all_comb_b2, (nucl_num, List_all_comb_b2_size)] +BEGIN_PROVIDER [integer, List_all_comb_b2, (nucl_num, List_all_comb_b2_size)] implicit none integer :: i, j @@ -50,57 +67,79 @@ END_PROVIDER List_all_comb_b2_expo = 0.d0 List_all_comb_b2_cent = 0.d0 - do i = 1, List_all_comb_b2_size + if(j1b_type .eq. 3) then - 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 + do i = 1, List_all_comb_b2_size - 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)) ) + 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 + + 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 - 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 - ! --- + 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) - do i = 1, List_all_comb_b2_size + 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 - phase = 0 - do j = 1, nucl_num - phase += List_all_comb_b2(j,i) + 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 - List_all_comb_b2_coef(i) = (-1.d0)**dble(phase) * dexp(-List_all_comb_b2_coef(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 + + elseif(j1b_type .eq. 4) then + + List_all_comb_b2_coef( 1) = 1.d0 + List_all_comb_b2_expo( 1) = 0.d0 + List_all_comb_b2_cent(1:3,1) = 0.d0 + do i = 1, nucl_num + List_all_comb_b2_coef( i+1) = -1.d0 + List_all_comb_b2_expo( i+1) = j1b_pen( i) + List_all_comb_b2_cent(1,i+1) = nucl_coord(i,1) + List_all_comb_b2_cent(2,i+1) = nucl_coord(i,2) + List_all_comb_b2_cent(3,i+1) = nucl_coord(i,3) + enddo + + else + + print *, 'j1b_type = ', j1b_pen, 'is not implemented' + stop + + endif !print *, ' coeff, expo & cent of list b2' !do i = 1, List_all_comb_b2_size @@ -115,14 +154,31 @@ END_PROVIDER BEGIN_PROVIDER [ integer, List_all_comb_b3_size] implicit none + double precision :: tmp - List_all_comb_b3_size = 3**nucl_num + if(j1b_type .eq. 3) then + + List_all_comb_b3_size = 3**nucl_num + + elseif(j1b_type .eq. 4) then + + tmp = 0.5d0 * dble(nucl_num) * (dble(nucl_num) + 3.d0) + List_all_comb_b3_size = int(tmp) + 1 + + else + + print *, 'j1b_type = ', j1b_pen, 'is not implemented' + stop + + endif + + print *, ' nb of linear terms in the square of the envelope is ', List_all_comb_b3_size END_PROVIDER ! --- -BEGIN_PROVIDER [ integer, List_all_comb_b3, (nucl_num, List_all_comb_b3_size)] +BEGIN_PROVIDER [integer, List_all_comb_b3, (nucl_num, List_all_comb_b3_size)] implicit none integer :: i, j, ii, jj @@ -162,7 +218,11 @@ END_PROVIDER implicit none integer :: i, j, k, phase + integer :: ii double precision :: tmp_alphaj, tmp_alphak, facto + double precision :: tmp1, tmp2, tmp3, tmp4 + double precision :: xi, yi, zi, xj, yj, zj + double precision :: dx, dy, dz, r2 provide j1b_pen @@ -170,60 +230,126 @@ END_PROVIDER List_all_comb_b3_expo = 0.d0 List_all_comb_b3_cent = 0.d0 - do i = 1, List_all_comb_b3_size + if(j1b_type .eq. 3) then - 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) + 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 + + if(List_all_comb_b3_expo(i) .lt. 1d-10) cycle + ASSERT(List_all_comb_b3_expo(i) .gt. 0d0) + + 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 - if(List_all_comb_b3_expo(i) .lt. 1d-10) cycle - ASSERT(List_all_comb_b3_expo(i) .gt. 0d0) + ! --- - 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) - do i = 1, List_all_comb_b3_size + 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 - 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) + if(List_all_comb_b3_expo(i) .lt. 1d-10) cycle - 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)) ) + 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 + + elseif(j1b_type .eq. 4) then + + ii = 1 + List_all_comb_b3_coef( ii) = 1.d0 + List_all_comb_b3_expo( ii) = 0.d0 + List_all_comb_b3_cent(1:3,ii) = 0.d0 + + do i = 1, nucl_num + ii = ii + 1 + List_all_comb_b3_coef( ii) = -2.d0 + List_all_comb_b3_expo( ii) = j1b_pen( i) + List_all_comb_b3_cent(1,ii) = nucl_coord(i,1) + List_all_comb_b3_cent(2,ii) = nucl_coord(i,2) + List_all_comb_b3_cent(3,ii) = nucl_coord(i,3) + enddo + + do i = 1, nucl_num + ii = ii + 1 + List_all_comb_b3_coef( ii) = 1.d0 + List_all_comb_b3_expo( ii) = 2.d0 * j1b_pen(i) + List_all_comb_b3_cent(1,ii) = nucl_coord(i,1) + List_all_comb_b3_cent(2,ii) = nucl_coord(i,2) + List_all_comb_b3_cent(3,ii) = nucl_coord(i,3) + enddo + + do i = 1, nucl_num-1 + + tmp1 = j1b_pen(i) + + xi = nucl_coord(i,1) + yi = nucl_coord(i,2) + zi = nucl_coord(i,3) + + do j = i+1, nucl_num + + tmp2 = j1b_pen(j) + tmp3 = tmp1 + tmp2 + tmp4 = 1.d0 / tmp3 + + xj = nucl_coord(j,1) + yj = nucl_coord(j,2) + zj = nucl_coord(j,3) + + dx = xi - xj + dy = yi - yj + dz = zi - zj + r2 = dx*dx + dy*dy + dz*dz + + ii = ii + 1 + List_all_comb_b3_coef( ii) = dexp(-tmp1*tmp2*tmp4*r2) + List_all_comb_b3_expo( ii) = tmp3 + List_all_comb_b3_cent(1,ii) = tmp4 * (tmp1 * xi + tmp2 * xj) + List_all_comb_b3_cent(2,ii) = tmp4 * (tmp1 * yi + tmp2 * yj) + List_all_comb_b3_cent(3,ii) = tmp4 * (tmp1 * zi + tmp2 * zj) enddo enddo - if(List_all_comb_b3_expo(i) .lt. 1d-10) cycle + else - List_all_comb_b3_coef(i) = List_all_comb_b3_coef(i) / List_all_comb_b3_expo(i) - enddo + print *, 'j1b_type = ', j1b_pen, 'is not implemented' + stop - ! --- - - 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 + endif !print *, ' coeff, expo & cent of list b3' !do i = 1, List_all_comb_b3_size diff --git a/src/non_h_ints_mu/grad_squared.irp.f b/src/non_h_ints_mu/grad_squared.irp.f index ea4bd36c..c93f62de 100644 --- a/src/non_h_ints_mu/grad_squared.irp.f +++ b/src/non_h_ints_mu/grad_squared.irp.f @@ -267,7 +267,7 @@ END_PROVIDER ! --- -BEGIN_PROVIDER [ double precision, u12sq_j1bsq, (ao_num, ao_num, n_points_final_grid) ] +BEGIN_PROVIDER [double precision, u12sq_j1bsq, (ao_num, ao_num, n_points_final_grid)] implicit none integer :: ipoint, i, j 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 c9f62b18..63a9db4b 100644 --- a/src/non_h_ints_mu/j12_nucl_utils.irp.f +++ b/src/non_h_ints_mu/j12_nucl_utils.irp.f @@ -8,79 +8,160 @@ BEGIN_PROVIDER [ double precision, v_1b, (n_points_final_grid)] double precision :: x, y, z, dx, dy, dz double precision :: a, d, e, fact_r - do ipoint = 1, n_points_final_grid + if(j1b_type .eq. 3) then - x = final_grid_points(1,ipoint) - y = final_grid_points(2,ipoint) - z = final_grid_points(3,ipoint) + ! v(r) = \Pi_{a} [1 - \exp(-\alpha_a (r - r_a)^2)] - 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) + do ipoint = 1, n_points_final_grid - fact_r = fact_r * e + 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 - v_1b(ipoint) = fact_r - enddo + elseif(j1b_type .eq. 4) then + + ! v(r) = 1 - \sum_{a} \exp(-\alpha_a (r - r_a)^2) + + 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 + + fact_r = fact_r - dexp(-a*d) + enddo + + v_1b(ipoint) = fact_r + enddo + + else + + print*, 'j1b_type = ', j1b_pen, 'is not implemented' + stop + + endif END_PROVIDER ! --- -BEGIN_PROVIDER [ double precision, v_1b_grad, (3, n_points_final_grid)] +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 :: x, y, z, dx, dy, dz, r2 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 + PROVIDE j1b_type - x = final_grid_points(1,ipoint) - y = final_grid_points(2,ipoint) - z = final_grid_points(3,ipoint) + if(j1b_type .eq. 3) then - fact_x = 0.d0 - fact_y = 0.d0 - fact_z = 0.d0 - do i = 1, List_all_comb_b2_size + ! v(r) = \Pi_{a} [1 - \exp(-\alpha_a (r - r_a)^2)] - phase = 0 - a_expo = 0.d0 - ax_der = 0.d0 - ay_der = 0.d0 - az_der = 0.d0 + 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 + + elseif(j1b_type .eq. 4) then + + ! v(r) = 1 - \sum_{a} \exp(-\alpha_a (r - r_a)^2) + + 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) + + 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) + r2 = dx*dx + dy*dy + dz*dz - fact_x += e * ax_der - fact_y += e * ay_der - fact_z += e * az_der + a = j1b_pen(j) + e = a * dexp(-a * r2) + + ax_der += e * dx + ay_der += e * dy + az_der += e * dz + enddo + + v_1b_grad(1,ipoint) = 2.d0 * ax_der + v_1b_grad(2,ipoint) = 2.d0 * ay_der + v_1b_grad(3,ipoint) = 2.d0 * 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 + else + + print*, 'j1b_type = ', j1b_pen, 'is not implemented' + stop + + endif 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 1a86a2e7..17c24d4b 100644 --- a/src/non_h_ints_mu/tc_integ.irp.f +++ b/src/non_h_ints_mu/tc_integ.irp.f @@ -68,7 +68,7 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_ao, (ao_num, ao_num, n_points_f !$OMP END DO !$OMP END PARALLEL - elseif(j1b_type .eq. 3) then + 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 @@ -219,7 +219,7 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_square_ao, (ao_num, ao_num, n_p !$OMP END DO !$OMP END PARALLEL - elseif(j1b_type .eq. 3) then + elseif((j1b_type .eq. 3) .or. (j1b_type .eq. 4)) then PROVIDE u12sq_j1bsq u12_grad1_u12_j1b_grad1_j1b grad12_j12