fixed bugs with dummy atom and becke grid

src/becke_numerical_grid/atomic_number.irp.f

@@ -0,0 +1,9 @@
+BEGIN_PROVIDER [ integer, grid_atomic_number, (nucl_num) ]
+ implicit none
+ BEGIN_DOC
+ ! Atomic number used to adjust the grid
+ END_DOC
+ grid_atomic_number(:) = max(1,int(nucl_charge(:)))
+
+END_PROVIDER
+

src/becke_numerical_grid/grid_becke.irp.f

@@ -146,7 +146,7 @@ BEGIN_PROVIDER [double precision, grid_points_per_atom, (3,n_points_integration_
       x = grid_points_radial(j)
 
       ! value of the radial coordinate for the integration
-      r = knowles_function(alpha_knowles(int(nucl_charge(i))),m_knowles,x)
+      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
@@ -232,8 +232,8 @@ BEGIN_PROVIDER [double precision, final_weight_at_r, (n_points_integration_angul
     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(int(nucl_charge(j))),m_knowles,x)&
-            *knowles_function(alpha_knowles(int(nucl_charge(j))),m_knowles,x)**2
+        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))' 

src/becke_numerical_grid/grid_becke_vector.irp.f

@@ -1,5 +1,6 @@
 
 BEGIN_PROVIDER [integer, n_points_final_grid]
+  implicit none
   BEGIN_DOC
   ! Number of points which are non zero
   END_DOC
@@ -8,7 +9,7 @@ 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)) < tresh_grid)then
+        if(dabs(final_weight_at_r(k,i,j)) < thresh_grid)then
           cycle
         endif
         n_points_final_grid += 1

src/becke_numerical_grid/step_function_becke.irp.f

@@ -31,10 +31,6 @@ double precision function cell_function_becke(r,atom_number)
   double precision               :: mu_ij,nu_ij
   double precision               :: distance_i,distance_j,step_function_becke
   integer                        :: j
-  if(int(nucl_charge(atom_number))==0)then
-   cell_function_becke = 0.d0
-   return
-  endif
   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))
@@ -42,7 +38,6 @@ double precision function cell_function_becke(r,atom_number)
   cell_function_becke = 1.d0
   do j = 1, nucl_num
     if(j==atom_number)cycle
-    if(int(nucl_charge(j))==0)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))

src/nuclei/atomic_radii.irp.f

@@ -66,7 +66,7 @@ BEGIN_PROVIDER [double precision, slater_bragg_radii_per_atom, (nucl_num)]
  implicit none
  integer :: i
  do i = 1, nucl_num
-  slater_bragg_radii_per_atom(i) = slater_bragg_radii(int(nucl_charge(i)))
+  slater_bragg_radii_per_atom(i) = slater_bragg_radii(max(1,int(nucl_charge(i))))
  enddo
 END_PROVIDER
 
@@ -74,7 +74,7 @@ BEGIN_PROVIDER [double precision, slater_bragg_radii_per_atom_ua, (nucl_num)]
  implicit none
  integer :: i
  do i = 1, nucl_num
-  slater_bragg_radii_per_atom_ua(i) = slater_bragg_radii_ua(int(nucl_charge(i)))
+  slater_bragg_radii_per_atom_ua(i) = slater_bragg_radii_ua(max(1,int(nucl_charge(i))))
  enddo
 END_PROVIDER