diff --git a/src/ao_two_e_ints/map_integrals.irp.f b/src/ao_two_e_ints/map_integrals.irp.f index 01b923e5..c6b01b07 100644 --- a/src/ao_two_e_ints/map_integrals.irp.f +++ b/src/ao_two_e_ints/map_integrals.irp.f @@ -324,6 +324,52 @@ subroutine get_ao_two_e_integrals_non_zero_jl(j,l,thresh,sze_max,sze,out_val,out end +subroutine get_ao_two_e_integrals_non_zero_jl_from_list(j,l,thresh,list,n_list,sze_max,out_val,out_val_index,non_zero_int) + use map_module + implicit none + BEGIN_DOC + ! Gets multiple AO bi-electronic integral from the AO map . + ! All non-zero i are retrieved for j,k,l fixed. + END_DOC + double precision, intent(in) :: thresh + integer, intent(in) :: j,l, n_list,list(2,sze_max),sze_max + real(integral_kind), intent(out) :: out_val(sze_max) + integer, intent(out) :: out_val_index(2,sze_max),non_zero_int + + integer :: i,k + integer(key_kind) :: hash + double precision :: tmp + + PROVIDE ao_two_e_integrals_in_map + non_zero_int = 0 + if (ao_overlap_abs(j,l) < thresh) then + out_val = 0.d0 + return + endif + + non_zero_int = 0 + integer :: kk + do kk = 1, n_list + k = list(1,kk) + i = list(2,kk) + integer, external :: ao_l4 + double precision, external :: ao_two_e_integral + !DIR$ FORCEINLINE + if (ao_two_e_integral_schwartz(i,k)*ao_two_e_integral_schwartz(j,l) < thresh) then + cycle + endif + call two_e_integrals_index(i,j,k,l,hash) + call map_get(ao_integrals_map, hash,tmp) + if (dabs(tmp) < thresh ) cycle + non_zero_int = non_zero_int+1 + out_val_index(1,non_zero_int) = i + out_val_index(2,non_zero_int) = k + out_val(non_zero_int) = tmp + enddo + +end + + function get_ao_map_size() diff --git a/src/becke_numerical_grid/EZFIO.cfg b/src/becke_numerical_grid/EZFIO.cfg index ed89428c..ca2100a1 100644 --- a/src/becke_numerical_grid/EZFIO.cfg +++ b/src/becke_numerical_grid/EZFIO.cfg @@ -8,3 +8,9 @@ default: 2 type: integer doc: Total number of grid points interface: ezfio + +[thresh_grid] +type: double precision +doc: threshold on the weight of a given grid point +interface: ezfio,provider,ocaml +default: 1.e-20 diff --git a/src/becke_numerical_grid/grid_becke_per_atom.irp.f b/src/becke_numerical_grid/grid_becke_per_atom.irp.f new file mode 100644 index 00000000..6026350b --- /dev/null +++ b/src/becke_numerical_grid/grid_becke_per_atom.irp.f @@ -0,0 +1,53 @@ + + + BEGIN_PROVIDER [integer, n_pts_per_atom, (nucl_num)] +&BEGIN_PROVIDER [integer, n_pts_max_per_atom] + BEGIN_DOC + ! Number of points which are non zero + END_DOC + integer :: i,j,k,l + n_pts_per_atom = 0 + do j = 1, nucl_num + do i = 1, n_points_radial_grid -1 + do k = 1, n_points_integration_angular + if(dabs(final_weight_at_r(k,i,j)) < thresh_grid)then + cycle + endif + n_pts_per_atom(j) += 1 + enddo + enddo + enddo + n_pts_max_per_atom = maxval(n_pts_per_atom) +END_PROVIDER + + BEGIN_PROVIDER [double precision, final_grid_points_per_atom, (3,n_pts_max_per_atom,nucl_num)] +&BEGIN_PROVIDER [double precision, final_weight_at_r_vector_per_atom, (n_pts_max_per_atom,nucl_num) ] +&BEGIN_PROVIDER [integer, index_final_points_per_atom, (3,n_pts_max_per_atom,nucl_num) ] +&BEGIN_PROVIDER [integer, index_final_points_per_atom_reverse, (n_points_integration_angular,n_points_radial_grid,nucl_num) ] + implicit none + integer :: i,j,k,l,i_count(nucl_num) + double precision :: r(3) + i_count = 0 + do j = 1, nucl_num + do i = 1, n_points_radial_grid -1 + do k = 1, n_points_integration_angular + if(dabs(final_weight_at_r(k,i,j)) < thresh_grid)then + cycle + endif + i_count(j) += 1 + final_grid_points_per_atom(1,i_count(j),j) = grid_points_per_atom(1,k,i,j) + final_grid_points_per_atom(2,i_count(j),j) = grid_points_per_atom(2,k,i,j) + final_grid_points_per_atom(3,i_count(j),j) = grid_points_per_atom(3,k,i,j) + final_weight_at_r_vector_per_atom(i_count(j),j) = final_weight_at_r(k,i,j) + index_final_points_per_atom(1,i_count(j),j) = k + index_final_points_per_atom(2,i_count(j),j) = i + index_final_points_per_atom(3,i_count(j),j) = j + index_final_points_per_atom_reverse(k,i,j) = i_count(j) + enddo + enddo + enddo + + + + +END_PROVIDER diff --git a/src/becke_numerical_grid/grid_becke_vector.irp.f b/src/becke_numerical_grid/grid_becke_vector.irp.f index a595cd0b..9856bcbd 100644 --- a/src/becke_numerical_grid/grid_becke_vector.irp.f +++ b/src/becke_numerical_grid/grid_becke_vector.irp.f @@ -8,9 +8,9 @@ BEGIN_PROVIDER [integer, n_points_final_grid] do j = 1, nucl_num do i = 1, n_points_radial_grid -1 do k = 1, n_points_integration_angular -! if(dabs(final_weight_at_r(k,i,j)) < 1.d-30)then -! cycle -! endif + if(dabs(final_weight_at_r(k,i,j)) < tresh_grid)then + cycle + endif n_points_final_grid += 1 enddo enddo @@ -39,9 +39,9 @@ END_PROVIDER do j = 1, nucl_num do i = 1, n_points_radial_grid -1 do k = 1, n_points_integration_angular - !if(dabs(final_weight_at_r(k,i,j)) < 1.d-30)then - ! cycle - !endif + if(dabs(final_weight_at_r(k,i,j)) < thresh_grid)then + cycle + endif i_count += 1 final_grid_points(1,i_count) = grid_points_per_atom(1,k,i,j) final_grid_points(2,i_count) = grid_points_per_atom(2,k,i,j) diff --git a/src/dft_utils_in_r/ao_in_r.irp.f b/src/dft_utils_in_r/ao_in_r.irp.f index 17892832..767f329c 100644 --- a/src/dft_utils_in_r/ao_in_r.irp.f +++ b/src/dft_utils_in_r/ao_in_r.irp.f @@ -121,3 +121,26 @@ enddo END_PROVIDER + + BEGIN_PROVIDER[double precision, aos_in_r_array_per_atom, (ao_num,n_pts_max_per_atom,nucl_num)] +&BEGIN_PROVIDER[double precision, aos_in_r_array_per_atom_transp, (n_pts_max_per_atom,ao_num,nucl_num)] + implicit none + BEGIN_DOC + ! aos_in_r_array_per_atom(i,j,k) = value of the ith ao on the jth grid point attached on the kth atom + END_DOC + integer :: i,j,k + double precision :: aos_array(ao_num), r(3) + do k = 1, nucl_num + do i = 1, n_pts_per_atom(k) + r(1) = final_grid_points_per_atom(1,i,k) + r(2) = final_grid_points_per_atom(2,i,k) + r(3) = final_grid_points_per_atom(3,i,k) + call give_all_aos_at_r(r,aos_array) + do j = 1, ao_num + aos_in_r_array_per_atom(j,i,k) = aos_array(j) + aos_in_r_array_per_atom_transp(i,j,k) = aos_array(j) + enddo + enddo + enddo + END_PROVIDER +