From 02968f569e6d3321c0ce54026dbb5fe431f69fa1 Mon Sep 17 00:00:00 2001
From: Emmanuel Giner <giner.emmanuel@gmail.com>
Date: Wed, 27 Mar 2019 12:56:32 +0100
Subject: [PATCH] fixed bug for dummy atoms X

---
 src/becke_numerical_grid/grid_becke.irp.f | 20 ++++++++++++++++++++
 src/dft_utils_in_r/dm_in_r.irp.f          |  2 ++
 src/nuclei/atomic_radii.irp.f             |  6 +++---
 3 files changed, 25 insertions(+), 3 deletions(-)

diff --git a/src/becke_numerical_grid/grid_becke.irp.f b/src/becke_numerical_grid/grid_becke.irp.f
index 4da9f4c9..38d4053f 100644
--- a/src/becke_numerical_grid/grid_becke.irp.f
+++ b/src/becke_numerical_grid/grid_becke.irp.f
@@ -195,6 +195,20 @@ 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
+         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
     enddo
   enddo
@@ -221,6 +235,12 @@ BEGIN_PROVIDER [double precision, final_weight_at_r, (n_points_integration_angul
         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
         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
       enddo
     enddo
   enddo
diff --git a/src/dft_utils_in_r/dm_in_r.irp.f b/src/dft_utils_in_r/dm_in_r.irp.f
index cc8dc4b4..6a19fed0 100644
--- a/src/dft_utils_in_r/dm_in_r.irp.f
+++ b/src/dft_utils_in_r/dm_in_r.irp.f
@@ -227,6 +227,8 @@ END_PROVIDER
 
  BEGIN_PROVIDER [double precision, one_e_dm_alpha_at_r, (n_points_final_grid,N_states) ]
 &BEGIN_PROVIDER [double precision, one_e_dm_beta_at_r, (n_points_final_grid,N_states) ]
+&BEGIN_PROVIDER [double precision, elec_beta_num_grid_becke , (N_states) ]
+&BEGIN_PROVIDER [double precision, elec_alpha_num_grid_becke , (N_states) ]
  implicit none
  BEGIN_DOC
 ! one_e_dm_alpha_at_r(i,istate) = n_alpha(r_i,istate)
diff --git a/src/nuclei/atomic_radii.irp.f b/src/nuclei/atomic_radii.irp.f
index 7210980d..82487b9d 100644
--- a/src/nuclei/atomic_radii.irp.f
+++ b/src/nuclei/atomic_radii.irp.f
@@ -1,4 +1,4 @@
-BEGIN_PROVIDER [ double precision, slater_bragg_radii, (100)]
+BEGIN_PROVIDER [ double precision, slater_bragg_radii, (0:100)]
  implicit none
  BEGIN_DOC
  ! atomic radii in Angstrom defined in table I of JCP 41, 3199 (1964) Slater
@@ -54,10 +54,10 @@ BEGIN_PROVIDER [ double precision, slater_bragg_radii, (100)]
 
 END_PROVIDER
 
-BEGIN_PROVIDER [double precision, slater_bragg_radii_ua, (100)]
+BEGIN_PROVIDER [double precision, slater_bragg_radii_ua, (0:100)]
  implicit none
  integer :: i
- do i = 1, 100
+ do i = 0, 100
   slater_bragg_radii_ua(i) = slater_bragg_radii(i) * 1.889725989d0
  enddo
 END_PROVIDER