diff --git a/src/Integrals_Monoelec/pot_ao_pseudo_ints.irp.f b/src/Integrals_Monoelec/pot_ao_pseudo_ints.irp.f
index f3efaa4c..615ed127 100644
--- a/src/Integrals_Monoelec/pot_ao_pseudo_ints.irp.f
+++ b/src/Integrals_Monoelec/pot_ao_pseudo_ints.irp.f
@@ -15,6 +15,7 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integral_local, (ao_num_align,ao_nu
   BEGIN_DOC
   ! Local pseudo-potential
   END_DOC
+  include 'Utils/constants.include.F'
   double precision               :: alpha, beta, gama, delta
   integer                        :: num_A,num_B
   double precision               :: A_center(3),B_center(3),C_center(3)
@@ -68,7 +69,7 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integral_local, (ao_num_align,ao_nu
           c = 0.d0
           
           if (dabs(ao_coef_normalized_ordered_transp(l,j)*ao_coef_normalized_ordered_transp(m,i))&
-                < ao_integrals_threshold) then
+                < thresh) then
             cycle
           endif
           do  k = 1, nucl_num
@@ -112,6 +113,7 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integral_local, (ao_num_align,ao_nu
   BEGIN_DOC
   ! Local pseudo-potential
   END_DOC
+  include 'Utils/constants.include.F'
   double precision               :: alpha, beta, gama, delta
   integer                        :: num_A,num_B
   double precision               :: A_center(3),B_center(3),C_center(3)
@@ -165,7 +167,7 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integral_local, (ao_num_align,ao_nu
           c = 0.d0
           
           if (dabs(ao_coef_normalized_ordered_transp(l,j)*ao_coef_normalized_ordered_transp(m,i))&
-                < ao_integrals_threshold) then
+                < thresh) then
             cycle
           endif
 
diff --git a/src/Utils/constants.include.F b/src/Utils/constants.include.F
index b55880c3..991ef80a 100644
--- a/src/Utils/constants.include.F
+++ b/src/Utils/constants.include.F
@@ -8,3 +8,4 @@ double precision, parameter :: dfour_pi =  4.d0*dacos(-1.d0)
 double precision, parameter :: dtwo_pi =  2.d0*dacos(-1.d0)
 double precision, parameter :: inv_sq_pi =  1.d0/dsqrt(dacos(-1.d0))
 double precision, parameter :: inv_sq_pi_2 = 0.5d0/dsqrt(dacos(-1.d0))
+double precision, parameter :: thresh = 1.d-15
diff --git a/src/Utils/integration.irp.f b/src/Utils/integration.irp.f
index ade130f7..ad57c52d 100644
--- a/src/Utils/integration.irp.f
+++ b/src/Utils/integration.irp.f
@@ -78,7 +78,7 @@ subroutine give_explicit_poly_and_gaussian(P_new,P_center,p,fact_k,iorder,alpha,
   
   !DEC$ FORCEINLINE
   call gaussian_product(alpha,A_center,beta,B_center,fact_k,p,P_center)
-  if (fact_k < ao_integrals_threshold) then
+  if (fact_k < thresh) then
     fact_k = 0.d0
     return
   endif