From c9938913789fcbd37290f020c19262ed05cf508e Mon Sep 17 00:00:00 2001 From: eginer Date: Wed, 2 Apr 2025 11:17:54 +0200 Subject: [PATCH] added test_e_pol_grid.irp.f --- plugins/local/extra_basis_int/NEED | 1 + .../extra_basis_int/test_e_pol_grid.irp.f | 33 +++++++ src/ao_extra_basis/aos.irp.f | 95 +++++++++++++++++++ src/ao_extra_basis/tuto/example_copy.sh | 2 + 4 files changed, 131 insertions(+) create mode 100644 plugins/local/extra_basis_int/test_e_pol_grid.irp.f diff --git a/plugins/local/extra_basis_int/NEED b/plugins/local/extra_basis_int/NEED index 76ae72cd..f32faa74 100644 --- a/plugins/local/extra_basis_int/NEED +++ b/plugins/local/extra_basis_int/NEED @@ -3,3 +3,4 @@ ao_one_e_ints ao_two_e_ints determinants becke_numerical_grid +dft_utils_in_r diff --git a/plugins/local/extra_basis_int/test_e_pol_grid.irp.f b/plugins/local/extra_basis_int/test_e_pol_grid.irp.f new file mode 100644 index 00000000..7e8882b1 --- /dev/null +++ b/plugins/local/extra_basis_int/test_e_pol_grid.irp.f @@ -0,0 +1,33 @@ +program test_e_pol + implicit none + io_extra_grid = "Read" + touch io_extra_grid + call routine +end + +subroutine routine + implicit none + integer :: i1,i2 + double precision :: integral_12, r1(3), r2(3), weight1, weight2 + double precision :: dm_r1_alpha,dm_r1_beta,dm_r1,dm_r2 + integral_12 = 0.d0 + do i1 = 1, n_points_final_grid + r1(1) = final_grid_points(1,i1) + r1(2) = final_grid_points(2,i1) + r1(3) = final_grid_points(3,i1) + weight1 = final_weight_at_r_vector(i1) + call dm_dft_alpha_beta_at_r(r1,dm_r1_alpha,dm_r1_beta) + dm_r1 = dm_r1_alpha+dm_r1_beta ! rhoA(r1) + do i2 = 1, n_points_extra_final_grid + r2(1) = final_grid_points_extra(1,i2) + r2(2) = final_grid_points_extra(2,i2) + r2(3) = final_grid_points_extra(3,i2) + weight2 = final_weight_at_r_vector_extra(i2) + dm_r2 = ao_extra_one_e_dm_at_extra_r(i2) ! rhoB(r2) + integral_12 += dm_r1 * dm_r2 * weight1 * weight2 + enddo + enddo + print*,'integral_12 = ',integral_12 + +end + diff --git a/src/ao_extra_basis/aos.irp.f b/src/ao_extra_basis/aos.irp.f index 56d6fb04..39a2e9b2 100644 --- a/src/ao_extra_basis/aos.irp.f +++ b/src/ao_extra_basis/aos.irp.f @@ -323,3 +323,98 @@ BEGIN_PROVIDER [ character*(4), ao_extra_l_char_space, (ao_extra_num) ] ao_extra_l_char_space(i) = give_ao_extra_character_space enddo END_PROVIDER + + + --- + +double precision function ao_extra_value(i, r) + + BEGIN_DOC + ! Returns the value of the i-th ao belonging the EXTRA BASIS at point $\textbf{r}$ + END_DOC + + implicit none + integer, intent(in) :: i + double precision, intent(in) :: r(3) + + integer :: m, num_ao + integer :: power_ao(3) + double precision :: center_ao(3) + double precision :: beta + double precision :: accu, dx, dy, dz, r2 + + num_ao = ao_extra_nucl(i) + power_ao(1:3) = ao_extra_power(i,1:3) + center_ao(1:3) = extra_nucl_coord(num_ao,1:3) + dx = r(1) - center_ao(1) + dy = r(2) - center_ao(2) + dz = r(3) - center_ao(3) + r2 = dx*dx + dy*dy + dz*dz + dx = dx**power_ao(1) + dy = dy**power_ao(2) + dz = dz**power_ao(3) + + accu = 0.d0 + do m = 1, ao_extra_prim_num(i) + beta = ao_extra_expo_ordered_transp(m,i) + accu += ao_extra_coef_normalized_ordered_transp(m,i) * dexp(-beta*r2) + enddo + ao_extra_value = accu * dx * dy * dz + +end + + +subroutine give_all_aos_extra_at_r(r, tmp_array) + implicit none + BEGIN_dOC + ! + ! input : r == r(1) = x and so on + ! + ! output : tmp_array(i) = EXTRA aos(i) evaluated in $\textbf{r}$ + ! + END_DOC + double precision, intent(in) :: r(3) + double precision, intent(out) :: tmp_array(ao_extra_num) + integer :: i + double precision :: ao_extra_value + do i = 1, ao_extra_num + tmp_array(i) = ao_extra_value(i, r) + enddo +end + +double precision function extra_density_at_r(r) + implicit none + BEGIN_dOC + ! + ! input : r == r(1) = x and so on + ! + ! output : density corresponding to the extra system at r + ! + END_DOC + double precision, intent(in) :: r(3) + integer :: mu,nu + double precision, allocatable :: tmp_array(:) + allocate(tmp_array(ao_extra_num)) + call give_all_aos_extra_at_r(r, tmp_array) + extra_density_at_r = 0.d0 + do nu = 1, ao_extra_num + do mu = 1, ao_extra_num + extra_density_at_r += ao_extra_one_e_dm(mu,nu,1) * tmp_array(mu) * tmp_array(nu) + enddo + enddo +end + +BEGIN_PROVIDER [ double precision, ao_extra_one_e_dm_at_extra_r, (n_points_extra_final_grid)] + implicit none + BEGIN_DOC +! ao_extra_one_e_dm_at_extra_r(i) = extra density on the extra grid points + END_DOC + integer :: i + double precision :: extra_density_at_r, r(3) + do i = 1, n_points_extra_final_grid + r(1) = final_grid_points_extra(1,i) + r(2) = final_grid_points_extra(2,i) + r(3) = final_grid_points_extra(3,i) + ao_extra_one_e_dm_at_extra_r(i) = extra_density_at_r(r) + enddo +END_PROVIDER diff --git a/src/ao_extra_basis/tuto/example_copy.sh b/src/ao_extra_basis/tuto/example_copy.sh index 0677b183..e08b2825 100755 --- a/src/ao_extra_basis/tuto/example_copy.sh +++ b/src/ao_extra_basis/tuto/example_copy.sh @@ -15,6 +15,8 @@ qp create_ezfio -b $basis_B $sys_B -o ${output_B} qp run scf # we save the density matrix in the EZFIO qp run save_one_e_dm +# we specify a small grid for the system B +qp set becke_numerical_grid extra_grid_type_sgn 0 # we create the system "A" qp create_ezfio -b $basis_A $sys_A -o ${output_A} # We perform an SCF calculation