From c9f579483af138575cb7a37b9cf1fdb39f22c2fa Mon Sep 17 00:00:00 2001 From: Abdallah Ammar Date: Mon, 15 May 2023 00:31:28 +0200 Subject: [PATCH] added Gill grid --- src/ao_many_one_e_ints/listj1b.irp.f | 24 +- src/becke_numerical_grid/EZFIO.cfg | 12 + src/becke_numerical_grid/extra_grid.irp.f | 278 +++++++++++----- .../extra_grid_vector.irp.f | 15 +- src/becke_numerical_grid/grid_becke.irp.f | 314 ++++++++++++------ .../grid_becke_vector.irp.f | 35 +- .../integration_radial.irp.f | 142 ++++---- .../step_function_becke.irp.f | 43 ++- src/non_h_ints_mu/j12_nucl_utils.irp.f | 4 +- 9 files changed, 570 insertions(+), 297 deletions(-) diff --git a/src/ao_many_one_e_ints/listj1b.irp.f b/src/ao_many_one_e_ints/listj1b.irp.f index 93ac459e..02963605 100644 --- a/src/ao_many_one_e_ints/listj1b.irp.f +++ b/src/ao_many_one_e_ints/listj1b.irp.f @@ -7,17 +7,17 @@ BEGIN_PROVIDER [integer, List_all_comb_b2_size] PROVIDE j1b_type - if(j1b_type .eq. 3) then + if((j1b_type .eq. 3) .or. (j1b_type .eq. 103)) then List_all_comb_b2_size = 2**nucl_num - elseif(j1b_type .eq. 4) then + elseif((j1b_type .eq. 4) .or. (j1b_type .eq. 104)) then List_all_comb_b2_size = nucl_num + 1 else - print *, 'j1b_type = ', j1b_pen, 'is not implemented' + print *, 'j1b_type = ', j1b_type, 'is not implemented' stop endif @@ -67,7 +67,7 @@ END_PROVIDER List_all_comb_b2_expo = 0.d0 List_all_comb_b2_cent = 0.d0 - if(j1b_type .eq. 3) then + if((j1b_type .eq. 3) .or. (j1b_type .eq. 103)) then do i = 1, List_all_comb_b2_size @@ -121,7 +121,7 @@ END_PROVIDER List_all_comb_b2_coef(i) = (-1.d0)**dble(phase) * dexp(-List_all_comb_b2_coef(i)) enddo - elseif(j1b_type .eq. 4) then + elseif((j1b_type .eq. 4) .or. (j1b_type .eq. 104)) then List_all_comb_b2_coef( 1) = 1.d0 List_all_comb_b2_expo( 1) = 0.d0 @@ -136,7 +136,7 @@ END_PROVIDER else - print *, 'j1b_type = ', j1b_pen, 'is not implemented' + print *, 'j1b_type = ', j1b_type, 'is not implemented' stop endif @@ -156,18 +156,18 @@ BEGIN_PROVIDER [ integer, List_all_comb_b3_size] implicit none double precision :: tmp - if(j1b_type .eq. 3) then + if((j1b_type .eq. 3) .or. (j1b_type .eq. 103)) then List_all_comb_b3_size = 3**nucl_num - elseif(j1b_type .eq. 4) then + elseif((j1b_type .eq. 4) .or. (j1b_type .eq. 104)) 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' + print *, 'j1b_type = ', j1b_type, 'is not implemented' stop endif @@ -230,7 +230,7 @@ END_PROVIDER List_all_comb_b3_expo = 0.d0 List_all_comb_b3_cent = 0.d0 - if(j1b_type .eq. 3) then + if((j1b_type .eq. 3) .or. (j1b_type .eq. 103)) then do i = 1, List_all_comb_b3_size @@ -287,7 +287,7 @@ END_PROVIDER 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 + elseif((j1b_type .eq. 4) .or. (j1b_type .eq. 104)) then ii = 1 List_all_comb_b3_coef( ii) = 1.d0 @@ -347,7 +347,7 @@ END_PROVIDER else - print *, 'j1b_type = ', j1b_pen, 'is not implemented' + print *, 'j1b_type = ', j1b_type, 'is not implemented' stop endif diff --git a/src/becke_numerical_grid/EZFIO.cfg b/src/becke_numerical_grid/EZFIO.cfg index 4083e0e7..7861f074 100644 --- a/src/becke_numerical_grid/EZFIO.cfg +++ b/src/becke_numerical_grid/EZFIO.cfg @@ -64,3 +64,15 @@ doc: Number of angular extra_grid points given from input. Warning, this number interface: ezfio,provider,ocaml default: 1202 +[rad_grid_type] +type: character*(32) +doc: method used to sample the radial space. Possible choices are [KNOWLES | GILL] +interface: ezfio,provider,ocaml +default: KNOWLES + +[extra_rad_grid_type] +type: character*(32) +doc: method used to sample the radial space. Possible choices are [KNOWLES | GILL] +interface: ezfio,provider,ocaml +default: KNOWLES + diff --git a/src/becke_numerical_grid/extra_grid.irp.f b/src/becke_numerical_grid/extra_grid.irp.f index db691e55..9bd24f22 100644 --- a/src/becke_numerical_grid/extra_grid.irp.f +++ b/src/becke_numerical_grid/extra_grid.irp.f @@ -1,96 +1,149 @@ +! --- + BEGIN_PROVIDER [integer, n_points_extra_radial_grid] &BEGIN_PROVIDER [integer, n_points_extra_integration_angular] - implicit none - BEGIN_DOC - ! n_points_extra_radial_grid = number of radial grid points_extra per atom - ! - ! n_points_extra_integration_angular = number of angular grid points_extra per atom - ! - ! These numbers are automatically set by setting the grid_type_sgn parameter - END_DOC -if(.not.my_extra_grid_becke)then - select case (extra_grid_type_sgn) - case(0) - n_points_extra_radial_grid = 23 - n_points_extra_integration_angular = 170 - case(1) - n_points_extra_radial_grid = 50 - n_points_extra_integration_angular = 194 - case(2) - n_points_extra_radial_grid = 75 - n_points_extra_integration_angular = 302 - case(3) - n_points_extra_radial_grid = 99 - n_points_extra_integration_angular = 590 - case default - write(*,*) '!!! Quadrature grid not available !!!' - stop - end select -else - n_points_extra_radial_grid = my_n_pt_r_extra_grid - n_points_extra_integration_angular = my_n_pt_a_extra_grid -endif + + BEGIN_DOC + ! n_points_extra_radial_grid = number of radial grid points_extra per atom + ! + ! n_points_extra_integration_angular = number of angular grid points_extra per atom + ! + ! These numbers are automatically set by setting the grid_type_sgn parameter + END_DOC + + implicit none + + if(.not.my_extra_grid_becke)then + select case (extra_grid_type_sgn) + case(0) + n_points_extra_radial_grid = 23 + n_points_extra_integration_angular = 170 + case(1) + n_points_extra_radial_grid = 50 + n_points_extra_integration_angular = 194 + case(2) + n_points_extra_radial_grid = 75 + n_points_extra_integration_angular = 302 + case(3) + n_points_extra_radial_grid = 99 + n_points_extra_integration_angular = 590 + case default + write(*,*) '!!! Quadrature grid not available !!!' + stop + end select + else + n_points_extra_radial_grid = my_n_pt_r_extra_grid + n_points_extra_integration_angular = my_n_pt_a_extra_grid + endif + END_PROVIDER +! --- + BEGIN_PROVIDER [integer, n_points_extra_grid_per_atom] - implicit none + BEGIN_DOC ! Number of grid points_extra per atom END_DOC + + implicit none n_points_extra_grid_per_atom = n_points_extra_integration_angular * n_points_extra_radial_grid END_PROVIDER +! --- + BEGIN_PROVIDER [double precision, grid_points_extra_radial, (n_points_extra_radial_grid)] &BEGIN_PROVIDER [double precision, dr_radial_extra_integral] - implicit none BEGIN_DOC ! points_extra in [0,1] to map the radial integral [0,\infty] END_DOC + + implicit none + integer :: i + dr_radial_extra_integral = 1.d0/dble(n_points_extra_radial_grid-1) - integer :: i do i = 1, n_points_extra_radial_grid grid_points_extra_radial(i) = dble(i-1) * dr_radial_extra_integral enddo END_PROVIDER +! --- + BEGIN_PROVIDER [double precision, grid_points_extra_per_atom, (3,n_points_extra_integration_angular,n_points_extra_radial_grid,nucl_num)] + BEGIN_DOC ! x,y,z coordinates of grid points_extra used for integration in 3d space END_DOC + implicit none - integer :: i,j,k - double precision :: dr,x_ref,y_ref,z_ref - double precision :: knowles_function - do i = 1, nucl_num - x_ref = nucl_coord(i,1) - y_ref = nucl_coord(i,2) - z_ref = nucl_coord(i,3) - do j = 1, n_points_extra_radial_grid-1 - double precision :: x,r - ! x value for the mapping of the [0, +\infty] to [0,1] - x = grid_points_extra_radial(j) + integer :: i, j, k + double precision :: dr, x_ref, y_ref, z_ref + double precision :: x, r, tmp + double precision, external :: knowles_function - ! value of the radial coordinate for the integration - r = knowles_function(alpha_knowles(grid_atomic_number(i)),m_knowles,x) + grid_points_extra_per_atom = 0.d0 - ! explicit values of the grid points_extra centered around each atom - do k = 1, n_points_extra_integration_angular - grid_points_extra_per_atom(1,k,j,i) = & - x_ref + angular_quadrature_points_extra(k,1) * r - grid_points_extra_per_atom(2,k,j,i) = & - y_ref + angular_quadrature_points_extra(k,2) * r - grid_points_extra_per_atom(3,k,j,i) = & - z_ref + angular_quadrature_points_extra(k,3) * r + PROVIDE extra_rad_grid_type + if(extra_rad_grid_type .eq. "KNOWLES") then + + do i = 1, nucl_num + x_ref = nucl_coord(i,1) + y_ref = nucl_coord(i,2) + z_ref = nucl_coord(i,3) + do j = 1, n_points_extra_radial_grid-1 + + ! x value for the mapping of the [0, +\infty] to [0,1] + x = grid_points_extra_radial(j) + ! value of the radial coordinate for the integration + r = knowles_function(alpha_knowles(grid_atomic_number(i)), m_knowles, x) + + ! explicit values of the grid points_extra centered around each atom + do k = 1, n_points_extra_integration_angular + grid_points_extra_per_atom(1,k,j,i) = x_ref + angular_quadrature_points_extra(k,1) * r + grid_points_extra_per_atom(2,k,j,i) = y_ref + angular_quadrature_points_extra(k,2) * r + grid_points_extra_per_atom(3,k,j,i) = z_ref + angular_quadrature_points_extra(k,3) * r + enddo enddo enddo - enddo + + elseif(extra_rad_grid_type .eq. "GILL") then + ! GILL & CHIEN, 2002 + + do i = 1, nucl_num + x_ref = nucl_coord(i,1) + y_ref = nucl_coord(i,2) + z_ref = nucl_coord(i,3) + do j = 1, n_points_extra_radial_grid-1 + + r = R_gill * dble(j-1)**2 / dble(n_points_extra_radial_grid-j+1)**2 + + ! explicit values of the grid points_extra centered around each atom + do k = 1, n_points_extra_integration_angular + grid_points_extra_per_atom(1,k,j,i) = x_ref + angular_quadrature_points_extra(k,1) * r + grid_points_extra_per_atom(2,k,j,i) = y_ref + angular_quadrature_points_extra(k,2) * r + grid_points_extra_per_atom(3,k,j,i) = z_ref + angular_quadrature_points_extra(k,3) * r + enddo + enddo + enddo + + else + + print*, " extra_rad_grid_type = ", extra_rad_grid_type, ' is not implemented' + stop + + endif + + END_PROVIDER +! --- + BEGIN_PROVIDER [double precision, weight_at_r_extra, (n_points_extra_integration_angular,n_points_extra_radial_grid,nucl_num) ] + BEGIN_DOC ! Weight function at grid points_extra : w_n(r) according to the equation (22) ! of Becke original paper (JCP, 88, 1988) @@ -99,11 +152,14 @@ BEGIN_PROVIDER [double precision, weight_at_r_extra, (n_points_extra_integration ! represented by the last dimension and the points_extra are labelled by the ! other dimensions. END_DOC + implicit none - integer :: i,j,k,l,m - double precision :: r(3) - double precision :: accu,cell_function_becke - double precision :: tmp_array(nucl_num) + integer :: i, j, k, l, m + double precision :: r(3) + double precision :: accu + double precision :: tmp_array(nucl_num) + double precision, external :: cell_function_becke + ! run over all points_extra in space ! that are referred to each atom do j = 1, nucl_num @@ -114,6 +170,7 @@ BEGIN_PROVIDER [double precision, weight_at_r_extra, (n_points_extra_integration r(1) = grid_points_extra_per_atom(1,l,k,j) r(2) = grid_points_extra_per_atom(2,l,k,j) r(3) = grid_points_extra_per_atom(3,l,k,j) + accu = 0.d0 ! For each of these points_extra in space, ou need to evaluate the P_n(r) do i = 1, nucl_num @@ -124,18 +181,19 @@ BEGIN_PROVIDER [double precision, weight_at_r_extra, (n_points_extra_integration enddo accu = 1.d0/accu weight_at_r_extra(l,k,j) = tmp_array(j) * accu + if(isnan(weight_at_r_extra(l,k,j)))then - print*,'isnan(weight_at_r_extra(l,k,j))' - print*,l,k,j - accu = 0.d0 - do i = 1, nucl_num - ! function defined for each atom "i" by equation (13) and (21) with k == 3 - tmp_array(i) = cell_function_becke(r,i) ! P_n(r) - print*,i,tmp_array(i) - ! Then you compute the summ the P_n(r) function for each of the "r" points_extra - accu += tmp_array(i) - enddo - write(*,'(100(F16.10,X))')tmp_array(j) , accu + print*,'isnan(weight_at_r_extra(l,k,j))' + print*,l,k,j + accu = 0.d0 + do i = 1, nucl_num + ! function defined for each atom "i" by equation (13) and (21) with k == 3 + tmp_array(i) = cell_function_becke(r,i) ! P_n(r) + print*,i,tmp_array(i) + ! Then you compute the summ the P_n(r) function for each of the "r" points_extra + accu += tmp_array(i) + enddo + write(*,'(100(F16.10,X))')tmp_array(j) , accu stop endif enddo @@ -144,35 +202,73 @@ BEGIN_PROVIDER [double precision, weight_at_r_extra, (n_points_extra_integration END_PROVIDER +! --- + BEGIN_PROVIDER [double precision, final_weight_at_r_extra, (n_points_extra_integration_angular,n_points_extra_radial_grid,nucl_num) ] + BEGIN_DOC ! Total weight on each grid point which takes into account all Lebedev, Voronoi and radial weights. END_DOC + implicit none - integer :: i,j,k,l,m - double precision :: r(3) - double precision :: accu,cell_function_becke - double precision :: tmp_array(nucl_num) - double precision :: contrib_integration,x - double precision :: derivative_knowles_function,knowles_function - ! run over all points_extra in space - do j = 1, nucl_num ! that are referred to each atom - do i = 1, n_points_extra_radial_grid -1 !for each radial grid attached to the "jth" atom - x = grid_points_extra_radial(i) ! x value for the mapping of the [0, +\infty] to [0,1] - do k = 1, n_points_extra_integration_angular ! for each angular point attached to the "jth" atom - contrib_integration = derivative_knowles_function(alpha_knowles(grid_atomic_number(j)),m_knowles,x)& - *knowles_function(alpha_knowles(grid_atomic_number(j)),m_knowles,x)**2 - final_weight_at_r_extra(k,i,j) = weights_angular_points_extra(k) * weight_at_r_extra(k,i,j) * contrib_integration * dr_radial_extra_integral - if(isnan(final_weight_at_r_extra(k,i,j)))then - print*,'isnan(final_weight_at_r_extra(k,i,j))' - print*,k,i,j - write(*,'(100(F16.10,X))')weights_angular_points_extra(k) , weight_at_r_extra(k,i,j) , contrib_integration , dr_radial_extra_integral - stop - endif + integer :: i, j, k, l, m + double precision :: r(3) + double precision :: tmp_array(nucl_num) + double precision :: contrib_integration, x, tmp + double precision, external :: derivative_knowles_function, knowles_function + + PROVIDE extra_rad_grid_type + if(extra_rad_grid_type .eq. "KNOWLES") then + + ! run over all points_extra in space + do j = 1, nucl_num ! that are referred to each atom + do i = 1, n_points_extra_radial_grid -1 !for each radial grid attached to the "jth" atom + x = grid_points_extra_radial(i) ! x value for the mapping of the [0, +\infty] to [0,1] + do k = 1, n_points_extra_integration_angular ! for each angular point attached to the "jth" atom + contrib_integration = derivative_knowles_function(alpha_knowles(grid_atomic_number(j)),m_knowles,x)& + * knowles_function(alpha_knowles(grid_atomic_number(j)),m_knowles,x)**2 + final_weight_at_r_extra(k,i,j) = weights_angular_points_extra(k) * weight_at_r_extra(k,i,j) * contrib_integration * dr_radial_extra_integral + if(isnan(final_weight_at_r_extra(k,i,j)))then + print*,'isnan(final_weight_at_r_extra(k,i,j))' + print*,k,i,j + write(*,'(100(F16.10,X))')weights_angular_points_extra(k) , weight_at_r_extra(k,i,j) , contrib_integration , dr_radial_extra_integral + stop + endif + enddo enddo enddo - enddo + + elseif(extra_rad_grid_type .eq. "GILL") then + ! GILL & CHIEN, 2002 + + PROVIDE R_gill + tmp = 2.d0 * R_gill * R_gill * R_gill * dble(n_points_extra_radial_grid) + + ! run over all points_extra in space + do j = 1, nucl_num ! that are referred to each atom + do i = 1, n_points_extra_radial_grid -1 !for each radial grid attached to the "jth" atom + contrib_integration = tmp * dble(i-1)**5 / dble(n_points_extra_radial_grid-i+1)**7 + + do k = 1, n_points_extra_integration_angular ! for each angular point attached to the "jth" atom + final_weight_at_r_extra(k,i,j) = weights_angular_points_extra(k) * weight_at_r_extra(k,i,j) * contrib_integration + if(isnan(final_weight_at_r_extra(k,i,j)))then + print*,'isnan(final_weight_at_r_extra(k,i,j))' + print*,k,i,j + write(*,'(100(F16.10,X))') weights_angular_points_extra(k), weight_at_r_extra(k,i,j), contrib_integration + stop + endif + enddo + enddo + enddo + + else + + print*, " extra_rad_grid_type = ", extra_rad_grid_type, ' is not implemented' + stop + + endif + END_PROVIDER diff --git a/src/becke_numerical_grid/extra_grid_vector.irp.f b/src/becke_numerical_grid/extra_grid_vector.irp.f index 3a5e6d3c..e4fc03b5 100644 --- a/src/becke_numerical_grid/extra_grid_vector.irp.f +++ b/src/becke_numerical_grid/extra_grid_vector.irp.f @@ -1,26 +1,35 @@ +! --- + BEGIN_PROVIDER [integer, n_points_extra_final_grid] - implicit none + BEGIN_DOC ! Number of points_extra which are non zero END_DOC - integer :: i,j,k,l + + implicit none + integer :: i, j, k, l + n_points_extra_final_grid = 0 + do j = 1, nucl_num do i = 1, n_points_extra_radial_grid -1 do k = 1, n_points_extra_integration_angular - if(dabs(final_weight_at_r_extra(k,i,j)) < thresh_extra_grid)then + if(dabs(final_weight_at_r_extra(k,i,j)) < thresh_extra_grid) then cycle endif n_points_extra_final_grid += 1 enddo enddo enddo + print*,'n_points_extra_final_grid = ',n_points_extra_final_grid print*,'n max point = ',n_points_extra_integration_angular*(n_points_extra_radial_grid*nucl_num - 1) ! call ezfio_set_becke_numerical_grid_n_points_extra_final_grid(n_points_extra_final_grid) END_PROVIDER +! --- + BEGIN_PROVIDER [double precision, final_grid_points_extra, (3,n_points_extra_final_grid)] &BEGIN_PROVIDER [double precision, final_weight_at_r_vector_extra, (n_points_extra_final_grid) ] &BEGIN_PROVIDER [integer, index_final_points_extra, (3,n_points_extra_final_grid) ] diff --git a/src/becke_numerical_grid/grid_becke.irp.f b/src/becke_numerical_grid/grid_becke.irp.f index 79f15c9a..21b9f98d 100644 --- a/src/becke_numerical_grid/grid_becke.irp.f +++ b/src/becke_numerical_grid/grid_becke.irp.f @@ -1,103 +1,174 @@ + +! --- + BEGIN_PROVIDER [integer, n_points_radial_grid] &BEGIN_PROVIDER [integer, n_points_integration_angular] - implicit none - BEGIN_DOC - ! n_points_radial_grid = number of radial grid points per atom - ! - ! n_points_integration_angular = number of angular grid points per atom - ! - ! These numbers are automatically set by setting the grid_type_sgn parameter - END_DOC -if(.not.my_grid_becke)then - select case (grid_type_sgn) - case(0) - n_points_radial_grid = 23 - n_points_integration_angular = 170 - case(1) - n_points_radial_grid = 50 - n_points_integration_angular = 194 - case(2) - n_points_radial_grid = 75 - n_points_integration_angular = 302 - case(3) - n_points_radial_grid = 99 - n_points_integration_angular = 590 - case default - write(*,*) '!!! Quadrature grid not available !!!' - stop - end select -else - n_points_radial_grid = my_n_pt_r_grid - n_points_integration_angular = my_n_pt_a_grid -endif + + BEGIN_DOC + ! n_points_radial_grid = number of radial grid points per atom + ! + ! n_points_integration_angular = number of angular grid points per atom + ! + ! These numbers are automatically set by setting the grid_type_sgn parameter + END_DOC + + implicit none + + if(.not.my_grid_becke)then + select case (grid_type_sgn) + case(0) + n_points_radial_grid = 23 + n_points_integration_angular = 170 + case(1) + n_points_radial_grid = 50 + n_points_integration_angular = 194 + case(2) + n_points_radial_grid = 75 + n_points_integration_angular = 302 + case(3) + n_points_radial_grid = 99 + n_points_integration_angular = 590 + case default + write(*,*) '!!! Quadrature grid not available !!!' + stop + end select + else + n_points_radial_grid = my_n_pt_r_grid + n_points_integration_angular = my_n_pt_a_grid + endif + END_PROVIDER +! --- + BEGIN_PROVIDER [integer, n_points_grid_per_atom] - implicit none + BEGIN_DOC ! Number of grid points per atom END_DOC + + implicit none + n_points_grid_per_atom = n_points_integration_angular * n_points_radial_grid END_PROVIDER -BEGIN_PROVIDER [integer , m_knowles] - implicit none +! --- + +BEGIN_PROVIDER [integer, m_knowles] + BEGIN_DOC ! value of the "m" parameter in the equation (7) of the paper of Knowles (JCP, 104, 1996) END_DOC + + implicit none + m_knowles = 3 + END_PROVIDER +! --- + +BEGIN_PROVIDER [double precision, R_gill] + + implicit none + + R_gill = 3.d0 + +END_PROVIDER + +! --- + BEGIN_PROVIDER [double precision, grid_points_radial, (n_points_radial_grid)] &BEGIN_PROVIDER [double precision, dr_radial_integral] - implicit none BEGIN_DOC ! points in [0,1] to map the radial integral [0,\infty] END_DOC - dr_radial_integral = 1.d0/dble(n_points_radial_grid-1) - integer :: i + + implicit none + integer :: i + + dr_radial_integral = 1.d0 / dble(n_points_radial_grid-1) + do i = 1, n_points_radial_grid grid_points_radial(i) = dble(i-1) * dr_radial_integral enddo END_PROVIDER +! --- + BEGIN_PROVIDER [double precision, grid_points_per_atom, (3,n_points_integration_angular,n_points_radial_grid,nucl_num)] + BEGIN_DOC ! x,y,z coordinates of grid points used for integration in 3d space END_DOC + implicit none - integer :: i,j,k - double precision :: dr,x_ref,y_ref,z_ref - double precision :: knowles_function - do i = 1, nucl_num - x_ref = nucl_coord(i,1) - y_ref = nucl_coord(i,2) - z_ref = nucl_coord(i,3) - do j = 1, n_points_radial_grid-1 - double precision :: x,r - ! x value for the mapping of the [0, +\infty] to [0,1] - x = grid_points_radial(j) + integer :: i, j, k + double precision :: dr, x_ref, y_ref, z_ref + double precision :: x, r, tmp + double precision, external :: knowles_function - ! value of the radial coordinate for the integration - r = knowles_function(alpha_knowles(grid_atomic_number(i)),m_knowles,x) + grid_points_per_atom = 0.d0 - ! explicit values of the grid points centered around each atom - do k = 1, n_points_integration_angular - grid_points_per_atom(1,k,j,i) = & - x_ref + angular_quadrature_points(k,1) * r - grid_points_per_atom(2,k,j,i) = & - y_ref + angular_quadrature_points(k,2) * r - grid_points_per_atom(3,k,j,i) = & - z_ref + angular_quadrature_points(k,3) * r + PROVIDE rad_grid_type + if(rad_grid_type .eq. "KNOWLES") then + + do i = 1, nucl_num + x_ref = nucl_coord(i,1) + y_ref = nucl_coord(i,2) + z_ref = nucl_coord(i,3) + do j = 1, n_points_radial_grid-1 + + ! x value for the mapping of the [0, +\infty] to [0,1] + x = grid_points_radial(j) + ! value of the radial coordinate for the integration + r = knowles_function(alpha_knowles(grid_atomic_number(i)), m_knowles, x) + + ! explicit values of the grid points centered around each atom + do k = 1, n_points_integration_angular + grid_points_per_atom(1,k,j,i) = x_ref + angular_quadrature_points(k,1) * r + grid_points_per_atom(2,k,j,i) = y_ref + angular_quadrature_points(k,2) * r + grid_points_per_atom(3,k,j,i) = z_ref + angular_quadrature_points(k,3) * r + enddo enddo enddo - enddo + + elseif(rad_grid_type .eq. "GILL") then + ! GILL & CHIEN, 2002 + + do i = 1, nucl_num + x_ref = nucl_coord(i,1) + y_ref = nucl_coord(i,2) + z_ref = nucl_coord(i,3) + do j = 1, n_points_radial_grid-1 + + r = R_gill * dble(j-1)**2 / dble(n_points_radial_grid-j+1)**2 + + ! explicit values of the grid points centered around each atom + do k = 1, n_points_integration_angular + grid_points_per_atom(1,k,j,i) = x_ref + angular_quadrature_points(k,1) * r + grid_points_per_atom(2,k,j,i) = y_ref + angular_quadrature_points(k,2) * r + grid_points_per_atom(3,k,j,i) = z_ref + angular_quadrature_points(k,3) * r + enddo + enddo + enddo + + else + + print*, " rad_grid_type = ", rad_grid_type, ' is not implemented' + stop + + endif + END_PROVIDER -BEGIN_PROVIDER [double precision, weight_at_r, (n_points_integration_angular,n_points_radial_grid,nucl_num) ] +! --- + +BEGIN_PROVIDER [double precision, weight_at_r, (n_points_integration_angular,n_points_radial_grid,nucl_num)] + BEGIN_DOC ! Weight function at grid points : w_n(r) according to the equation (22) ! of Becke original paper (JCP, 88, 1988) @@ -106,11 +177,13 @@ BEGIN_PROVIDER [double precision, weight_at_r, (n_points_integration_angular,n_p ! represented by the last dimension and the points are labelled by the ! other dimensions. END_DOC + implicit none - integer :: i,j,k,l,m - double precision :: r(3) - double precision :: accu,cell_function_becke - double precision :: tmp_array(nucl_num) + integer :: i, j, k, l, m + double precision :: r(3), accu + double precision :: tmp_array(nucl_num) + double precision, external :: cell_function_becke + ! run over all points in space ! that are referred to each atom do j = 1, nucl_num @@ -121,28 +194,30 @@ BEGIN_PROVIDER [double precision, weight_at_r, (n_points_integration_angular,n_p r(1) = grid_points_per_atom(1,l,k,j) r(2) = grid_points_per_atom(2,l,k,j) r(3) = grid_points_per_atom(3,l,k,j) + accu = 0.d0 ! For each of these points in space, ou need to evaluate the P_n(r) do i = 1, nucl_num ! function defined for each atom "i" by equation (13) and (21) with k == 3 - tmp_array(i) = cell_function_becke(r,i) ! P_n(r) + tmp_array(i) = cell_function_becke(r, i) ! P_n(r) ! Then you compute the summ the P_n(r) function for each of the "r" points accu += tmp_array(i) enddo accu = 1.d0/accu weight_at_r(l,k,j) = tmp_array(j) * accu - if(isnan(weight_at_r(l,k,j)))then - print*,'isnan(weight_at_r(l,k,j))' - print*,l,k,j - accu = 0.d0 - do i = 1, nucl_num - ! function defined for each atom "i" by equation (13) and (21) with k == 3 - tmp_array(i) = cell_function_becke(r,i) ! P_n(r) - print*,i,tmp_array(i) - ! Then you compute the summ the P_n(r) function for each of the "r" points - accu += tmp_array(i) - enddo - write(*,'(100(F16.10,X))')tmp_array(j) , accu + + if(isnan(weight_at_r(l,k,j))) then + print*,'isnan(weight_at_r(l,k,j))' + print*,l,k,j + accu = 0.d0 + do i = 1, nucl_num + ! function defined for each atom "i" by equation (13) and (21) with k == 3 + tmp_array(i) = cell_function_becke(r,i) ! P_n(r) + print*,i,tmp_array(i) + ! Then you compute the summ the P_n(r) function for each of the "r" points + accu += tmp_array(i) + enddo + write(*,'(100(F16.10,X))')tmp_array(j) , accu stop endif enddo @@ -151,35 +226,76 @@ BEGIN_PROVIDER [double precision, weight_at_r, (n_points_integration_angular,n_p END_PROVIDER +! --- + +BEGIN_PROVIDER [double precision, final_weight_at_r, (n_points_integration_angular,n_points_radial_grid,nucl_num)] -BEGIN_PROVIDER [double precision, final_weight_at_r, (n_points_integration_angular,n_points_radial_grid,nucl_num) ] BEGIN_DOC - ! Total weight on each grid point which takes into account all Lebedev, Voronoi and radial weights. + ! Total weight on each grid point which takes into account all Lebedev, Voronoi and radial weights. END_DOC + implicit none - integer :: i,j,k,l,m - double precision :: r(3) - double precision :: accu,cell_function_becke - double precision :: tmp_array(nucl_num) - double precision :: contrib_integration,x - double precision :: derivative_knowles_function,knowles_function - ! run over all points in space - do j = 1, nucl_num ! that are referred to each atom - do i = 1, n_points_radial_grid -1 !for each radial grid attached to the "jth" atom - x = grid_points_radial(i) ! x value for the mapping of the [0, +\infty] to [0,1] - do k = 1, n_points_integration_angular ! for each angular point attached to the "jth" atom - contrib_integration = derivative_knowles_function(alpha_knowles(grid_atomic_number(j)),m_knowles,x)& - *knowles_function(alpha_knowles(grid_atomic_number(j)),m_knowles,x)**2 - final_weight_at_r(k,i,j) = weights_angular_points(k) * weight_at_r(k,i,j) * contrib_integration * dr_radial_integral - if(isnan(final_weight_at_r(k,i,j)))then - print*,'isnan(final_weight_at_r(k,i,j))' - print*,k,i,j - write(*,'(100(F16.10,X))')weights_angular_points(k) , weight_at_r(k,i,j) , contrib_integration , dr_radial_integral - stop - endif + integer :: i, j, k, l, m + double precision :: r(3) + double precision :: tmp_array(nucl_num) + double precision :: contrib_integration, x, tmp + double precision, external :: derivative_knowles_function, knowles_function + + final_weight_at_r = 0.d0 + + PROVIDE rad_grid_type + if(rad_grid_type .eq. "KNOWLES") then + + ! run over all points in space + do j = 1, nucl_num ! that are referred to each atom + do i = 1, n_points_radial_grid -1 !for each radial grid attached to the "jth" atom + x = grid_points_radial(i) ! x value for the mapping of the [0, +\infty] to [0,1] + + do k = 1, n_points_integration_angular ! for each angular point attached to the "jth" atom + contrib_integration = derivative_knowles_function(alpha_knowles(grid_atomic_number(j)), m_knowles, x) & + * knowles_function(alpha_knowles(grid_atomic_number(j)), m_knowles, x)**2 + + final_weight_at_r(k,i,j) = weights_angular_points(k) * weight_at_r(k,i,j) * contrib_integration * dr_radial_integral + + if(isnan(final_weight_at_r(k,i,j))) then + print*,'isnan(final_weight_at_r(k,i,j))' + print*,k,i,j + write(*,'(100(F16.10,X))') weights_angular_points(k), weight_at_r(k,i,j), contrib_integration + stop + endif + enddo enddo enddo - enddo + + elseif(rad_grid_type .eq. "GILL") then + ! GILL & CHIEN, 2002 + + tmp = 2.d0 * R_gill * R_gill * R_gill * dble(n_points_radial_grid) + + ! run over all points in space + do j = 1, nucl_num ! that are referred to each atom + do i = 1, n_points_radial_grid - 1 !for each radial grid attached to the "jth" atom + contrib_integration = tmp * dble(i-1)**5 / dble(n_points_radial_grid-i+1)**7 + do k = 1, n_points_integration_angular ! for each angular point attached to the "jth" atom + final_weight_at_r(k,i,j) = weights_angular_points(k) * weight_at_r(k,i,j) * contrib_integration + + if(isnan(final_weight_at_r(k,i,j))) then + print*,'isnan(final_weight_at_r(k,i,j))' + print*,k,i,j + write(*,'(100(F16.10,X))') weights_angular_points(k), weight_at_r(k,i,j), contrib_integration, dr_radial_integral + stop + endif + enddo + enddo + enddo + + else + + print*, " rad_grid_type = ", rad_grid_type, ' is not implemented' + stop + + endif 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 343bd054..fd185641 100644 --- a/src/becke_numerical_grid/grid_becke_vector.irp.f +++ b/src/becke_numerical_grid/grid_becke_vector.irp.f @@ -21,22 +21,27 @@ BEGIN_PROVIDER [integer, n_points_final_grid] call ezfio_set_becke_numerical_grid_n_points_final_grid(n_points_final_grid) END_PROVIDER - BEGIN_PROVIDER [double precision, final_grid_points, (3,n_points_final_grid)] -&BEGIN_PROVIDER [double precision, final_weight_at_r_vector, (n_points_final_grid) ] -&BEGIN_PROVIDER [integer, index_final_points, (3,n_points_final_grid) ] -&BEGIN_PROVIDER [integer, index_final_points_reverse, (n_points_integration_angular,n_points_radial_grid,nucl_num) ] - implicit none +! --- + + BEGIN_PROVIDER [double precision, final_grid_points, (3,n_points_final_grid)] +&BEGIN_PROVIDER [double precision, final_weight_at_r_vector, (n_points_final_grid)] +&BEGIN_PROVIDER [integer, index_final_points, (3,n_points_final_grid)] +&BEGIN_PROVIDER [integer, index_final_points_reverse, (n_points_integration_angular,n_points_radial_grid,nucl_num)] + BEGIN_DOC -! final_grid_points(1:3,j) = (/ x, y, z /) of the jth grid point -! -! final_weight_at_r_vector(i) = Total weight function of the ith grid point which contains the Lebedev, Voronoi and radial weights contributions -! -! index_final_points(1:3,i) = gives the angular, radial and atomic indices associated to the ith grid point -! -! index_final_points_reverse(i,j,k) = index of the grid point having i as angular, j as radial and l as atomic indices + ! final_grid_points(1:3,j) = (/ x, y, z /) of the jth grid point + ! + ! final_weight_at_r_vector(i) = Total weight function of the ith grid point which contains the Lebedev, Voronoi and radial weights contributions + ! + ! index_final_points(1:3,i) = gives the angular, radial and atomic indices associated to the ith grid point + ! + ! index_final_points_reverse(i,j,k) = index of the grid point having i as angular, j as radial and l as atomic indices END_DOC - integer :: i,j,k,l,i_count - double precision :: r(3) + + implicit none + integer :: i, j, k, l, i_count + double precision :: r(3) + i_count = 0 do j = 1, nucl_num do i = 1, n_points_radial_grid -1 @@ -59,6 +64,8 @@ END_PROVIDER END_PROVIDER +! --- + BEGIN_PROVIDER [double precision, final_grid_points_transp, (n_points_final_grid,3)] implicit none BEGIN_DOC diff --git a/src/becke_numerical_grid/integration_radial.irp.f b/src/becke_numerical_grid/integration_radial.irp.f index 44c83070..3de151ab 100644 --- a/src/becke_numerical_grid/integration_radial.irp.f +++ b/src/becke_numerical_grid/integration_radial.irp.f @@ -1,71 +1,93 @@ - double precision function knowles_function(alpha,m,x) - implicit none - BEGIN_DOC -! Function proposed by Knowles (JCP, 104, 1996) for distributing the radial points : -! the Log "m" function ( equation (7) in the paper ) - END_DOC - double precision, intent(in) :: alpha,x - integer, intent(in) :: m -!print*, x - knowles_function = -alpha * dlog(1.d0-x**m) - end - double precision function derivative_knowles_function(alpha,m,x) - implicit none - BEGIN_DOC -! Derivative of the function proposed by Knowles (JCP, 104, 1996) for distributing the radial points - END_DOC - double precision, intent(in) :: alpha,x - integer, intent(in) :: m - double precision :: f - f = x**(m-1) - derivative_knowles_function = alpha * dble(m) * f / (1.d0 - x*f) - end +! --- - BEGIN_PROVIDER [double precision, alpha_knowles, (100)] - implicit none - integer :: i - BEGIN_DOC -! Recommended values for the alpha parameters according to the paper of Knowles (JCP, 104, 1996) -! as a function of the nuclear charge - END_DOC +double precision function knowles_function(alpha, m, x) - ! H-He - alpha_knowles(1) = 5.d0 - alpha_knowles(2) = 5.d0 + BEGIN_DOC + ! Function proposed by Knowles (JCP, 104, 1996) for distributing the radial points : + ! the Log "m" function ( equation (7) in the paper ) + END_DOC + + implicit none + double precision, intent(in) :: alpha, x + integer, intent(in) :: m - ! Li-Be - alpha_knowles(3) = 7.d0 - alpha_knowles(4) = 7.d0 + !print*, x + knowles_function = -alpha * dlog(1.d0-x**m) - ! B-Ne - do i = 5, 10 - alpha_knowles(i) = 5.d0 - enddo + return +end - ! Na-Mg - do i = 11, 12 - alpha_knowles(i) = 7.d0 - enddo +! --- - ! Al-Ar - do i = 13, 18 - alpha_knowles(i) = 5.d0 - enddo +double precision function derivative_knowles_function(alpha, m, x) - ! K-Ca - do i = 19, 20 - alpha_knowles(i) = 7.d0 - enddo + BEGIN_DOC + ! Derivative of the function proposed by Knowles (JCP, 104, 1996) for distributing the radial points + END_DOC - ! Sc-Zn - do i = 21, 30 - alpha_knowles(i) = 5.d0 - enddo + implicit none + double precision, intent(in) :: alpha, x + integer, intent(in) :: m + double precision :: f - ! Ga-Kr - do i = 31, 100 - alpha_knowles(i) = 7.d0 - enddo + f = x**(m-1) + derivative_knowles_function = alpha * dble(m) * f / (1.d0 - x*f) + + return +end + +! --- + +BEGIN_PROVIDER [double precision, alpha_knowles, (100)] + + BEGIN_DOC + ! Recommended values for the alpha parameters according to the paper of Knowles (JCP, 104, 1996) + ! as a function of the nuclear charge + END_DOC + + implicit none + integer :: i + + ! H-He + alpha_knowles(1) = 5.d0 + alpha_knowles(2) = 5.d0 + + ! Li-Be + alpha_knowles(3) = 7.d0 + alpha_knowles(4) = 7.d0 + + ! B-Ne + do i = 5, 10 + alpha_knowles(i) = 5.d0 + enddo + + ! Na-Mg + do i = 11, 12 + alpha_knowles(i) = 7.d0 + enddo + + ! Al-Ar + do i = 13, 18 + alpha_knowles(i) = 5.d0 + enddo + + ! K-Ca + do i = 19, 20 + alpha_knowles(i) = 7.d0 + enddo + + ! Sc-Zn + do i = 21, 30 + alpha_knowles(i) = 5.d0 + enddo + + ! Ga-Kr + do i = 31, 100 + alpha_knowles(i) = 7.d0 + enddo + +END_PROVIDER + +! --- - END_PROVIDER diff --git a/src/becke_numerical_grid/step_function_becke.irp.f b/src/becke_numerical_grid/step_function_becke.irp.f index 2905c6c0..6048c35f 100644 --- a/src/becke_numerical_grid/step_function_becke.irp.f +++ b/src/becke_numerical_grid/step_function_becke.irp.f @@ -20,31 +20,42 @@ double precision function f_function_becke(x) f_function_becke = 1.5d0 * x - 0.5d0 * x*x*x end -double precision function cell_function_becke(r,atom_number) - implicit none - double precision, intent(in) :: r(3) - integer, intent(in) :: atom_number +! --- + +double precision function cell_function_becke(r, atom_number) + BEGIN_DOC -! atom_number :: atom on which the cell function of Becke (1988, JCP,88(4)) + ! atom_number :: atom on which the cell function of Becke (1988, JCP,88(4)) ! r(1:3) :: x,y,z coordinantes of the current point END_DOC - double precision :: mu_ij,nu_ij - double precision :: distance_i,distance_j,step_function_becke - integer :: j - distance_i = (r(1) - nucl_coord_transp(1,atom_number) ) * (r(1) - nucl_coord_transp(1,atom_number)) + + implicit none + double precision, intent(in) :: r(3) + integer, intent(in) :: atom_number + integer :: j + double precision :: mu_ij, nu_ij + double precision :: distance_i, distance_j, step_function_becke + + distance_i = (r(1) - nucl_coord_transp(1,atom_number) ) * (r(1) - nucl_coord_transp(1,atom_number)) distance_i += (r(2) - nucl_coord_transp(2,atom_number) ) * (r(2) - nucl_coord_transp(2,atom_number)) distance_i += (r(3) - nucl_coord_transp(3,atom_number) ) * (r(3) - nucl_coord_transp(3,atom_number)) - distance_i = dsqrt(distance_i) + distance_i = dsqrt(distance_i) + cell_function_becke = 1.d0 do j = 1, nucl_num - if(j==atom_number)cycle - distance_j = (r(1) - nucl_coord_transp(1,j) ) * (r(1) - nucl_coord_transp(1,j)) - distance_j+= (r(2) - nucl_coord_transp(2,j) ) * (r(2) - nucl_coord_transp(2,j)) - distance_j+= (r(3) - nucl_coord_transp(3,j) ) * (r(3) - nucl_coord_transp(3,j)) - distance_j = dsqrt(distance_j) - mu_ij = (distance_i - distance_j)*nucl_dist_inv(atom_number,j) + if(j==atom_number) cycle + + distance_j = (r(1) - nucl_coord_transp(1,j) ) * (r(1) - nucl_coord_transp(1,j)) + distance_j += (r(2) - nucl_coord_transp(2,j) ) * (r(2) - nucl_coord_transp(2,j)) + distance_j += (r(3) - nucl_coord_transp(3,j) ) * (r(3) - nucl_coord_transp(3,j)) + distance_j = dsqrt(distance_j) + + mu_ij = (distance_i - distance_j) * nucl_dist_inv(atom_number,j) nu_ij = mu_ij + slater_bragg_type_inter_distance_ua(atom_number,j) * (1.d0 - mu_ij*mu_ij) + cell_function_becke *= step_function_becke(nu_ij) enddo + + return end 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 079cb388..9b91a8ed 100644 --- a/src/non_h_ints_mu/j12_nucl_utils.irp.f +++ b/src/non_h_ints_mu/j12_nucl_utils.irp.f @@ -59,7 +59,7 @@ BEGIN_PROVIDER [ double precision, v_1b, (n_points_final_grid)] else - print*, 'j1b_type = ', j1b_pen, 'is not implemented for v_1b' + print*, 'j1b_type = ', j1b_type, 'is not implemented for v_1b' stop endif @@ -158,7 +158,7 @@ BEGIN_PROVIDER [double precision, v_1b_grad, (3, n_points_final_grid)] else - print*, 'j1b_type = ', j1b_pen, 'is not implemented' + print*, 'j1b_type = ', j1b_type, 'is not implemented' stop endif