added Gill grid

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

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
+

src/becke_numerical_grid/extra_grid.irp.f

@@ -1,7 +1,9 @@
 
+! ---
+
  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
   !
@@ -9,6 +11,9 @@
   !
   ! 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)
@@ -31,66 +36,114 @@ 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
-  dr_radial_extra_integral = 1.d0/dble(n_points_extra_radial_grid-1)
+
+  implicit none
   integer :: i
+
+  dr_radial_extra_integral = 1.d0/dble(n_points_extra_radial_grid-1)
   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
+  double precision           :: x, r, tmp
+  double precision, external :: knowles_function
+
+  grid_points_extra_per_atom = 0.d0
+
+  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
-      double precision               :: x,r
+
         ! 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) =                              &
+          grid_points_extra_per_atom(1,k,j,i) = x_ref + angular_quadrature_points_extra(k,1) * r
-            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(2,k,j,i) =                              &
+          grid_points_extra_per_atom(3,k,j,i) = z_ref + angular_quadrature_points_extra(k,3) * r
-            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
+
+  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           :: 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,6 +181,7 @@ 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
@@ -144,18 +202,25 @@ 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           :: contrib_integration, x, tmp
-  double precision               :: derivative_knowles_function,knowles_function
+  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
@@ -174,5 +239,36 @@ BEGIN_PROVIDER [double precision, final_weight_at_r_extra, (n_points_extra_integ
       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
 

src/becke_numerical_grid/extra_grid_vector.irp.f

@@ -1,11 +1,17 @@
 
+! ---
+
 BEGIN_PROVIDER [integer, n_points_extra_final_grid]
-  implicit none
+
   BEGIN_DOC
   ! Number of points_extra which are non zero
   END_DOC
+
+  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
@@ -16,11 +22,14 @@ BEGIN_PROVIDER [integer, n_points_extra_final_grid]
       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) ]

src/becke_numerical_grid/grid_becke.irp.f

@@ -1,6 +1,9 @@
+
+! ---
+
  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
   !
@@ -8,6 +11,9 @@
   !
   ! 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)
@@ -30,74 +36,139 @@ 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_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)
+
+  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
+  double precision           :: x, r, tmp
+  double precision, external :: knowles_function
+
+  grid_points_per_atom = 0.d0
+
+  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
-      double precision               :: x,r
+
         ! 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) =                              &
+          grid_points_per_atom(1,k,j,i) = x_ref + angular_quadrature_points(k,1) * r
-            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(2,k,j,i) =                              &
+          grid_points_per_atom(3,k,j,i) = z_ref + angular_quadrature_points(k,3) * r
-            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
+
+  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_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           :: r(3), accu
-  double precision               :: accu,cell_function_becke
   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,6 +194,7 @@ 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
@@ -131,6 +205,7 @@ BEGIN_PROVIDER [double precision, weight_at_r, (n_points_integration_angular,n_p
         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
@@ -151,26 +226,59 @@ 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_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           :: contrib_integration, x, tmp
-  double precision               :: derivative_knowles_function,knowles_function
+  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
+
+  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
@@ -181,5 +289,13 @@ BEGIN_PROVIDER [double precision, final_weight_at_r, (n_points_integration_angul
       enddo
     enddo
 
+  else
+
+    print*, " rad_grid_type = ", rad_grid_type, ' is not implemented'
+    stop
+
+  endif
+
 END_PROVIDER
 
+

src/becke_numerical_grid/grid_becke_vector.irp.f

@@ -21,11 +21,13 @@ 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_DOC
   !  final_grid_points(1:3,j) = (/ x, y, z /) of the jth grid point
   !
@@ -35,8 +37,11 @@ END_PROVIDER
   !
   ! 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
+
+  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

src/becke_numerical_grid/integration_radial.irp.f

@@ -1,35 +1,54 @@
+
+! ---
+
 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
+ 
+  implicit none
   double precision, intent(in) :: alpha, x
   integer,          intent(in) :: m
+
   !print*, x
   knowles_function = -alpha * dlog(1.d0-x**m)
+
+  return
 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
+
+  implicit none
   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)
+
+  return
 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
 
+  implicit none
+  integer :: i
+
   ! H-He
   alpha_knowles(1) = 5.d0
   alpha_knowles(2) = 5.d0
@@ -69,3 +88,6 @@
   enddo
 
 END_PROVIDER
+
+! ---
+

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
   BEGIN_DOC
   ! 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
+
+  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
-  integer                        :: j
+
   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)
+
   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)
     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
 

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