diff --git a/src/ao_one_e_ints/ao_overlap.irp.f b/src/ao_one_e_ints/ao_overlap.irp.f
index 11c95e42..d9061d67 100644
--- a/src/ao_one_e_ints/ao_overlap.irp.f
+++ b/src/ao_one_e_ints/ao_overlap.irp.f
@@ -54,6 +54,13 @@
         call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,overlap_y,overlap_z,overlap,dim1)
         c = ao_coef_normalized_ordered_transp(n,j) * ao_coef_normalized_ordered_transp(l,i)
         ao_overlap(i,j) += c * overlap
+        if(isnan(ao_overlap(i,j)))then
+         print*,'i,j',i,j
+         print*,'l,n',l,n
+         print*,'c,overlap',c,overlap
+         print*,overlap_x,overlap_y,overlap_z
+         stop
+        endif
         ao_overlap_x(i,j) += c * overlap_x
         ao_overlap_y(i,j) += c * overlap_y
         ao_overlap_z(i,j) += c * overlap_z
diff --git a/src/utils/one_e_integration.irp.f b/src/utils/one_e_integration.irp.f
index 2c597104..cacc3bf7 100644
--- a/src/utils/one_e_integration.irp.f
+++ b/src/utils/one_e_integration.irp.f
@@ -15,10 +15,10 @@ double precision function overlap_gaussian_x(A_center,B_center,alpha,beta,power_
   call give_explicit_poly_and_gaussian_x(P_new,P_center,p,fact_p,iorder_p,alpha,&
       beta,power_A,power_B,A_center,B_center,dim)
 
-! if(fact_p.lt.0.000001d0)then
-!   overlap_gaussian_x = 0.d0
-!   return
-! endif
+  if(fact_p.lt.1.d-20)then
+    overlap_gaussian_x = 0.d0
+    return
+  endif
 
   overlap_gaussian_x = 0.d0
   integer                        :: i