From bd8fcdb03030af73d5f4bef74a6ccd0a8899702d Mon Sep 17 00:00:00 2001
From: Emmanuel Giner <giner.emmanuel@gmail.com>
Date: Mon, 2 Nov 2020 17:24:35 +0100
Subject: [PATCH] removed stupid bug in utils/one_e_integration.irp.f

---
 src/ao_one_e_ints/ao_overlap.irp.f |  7 +++++++
 src/utils/one_e_integration.irp.f  | 22 +++++++++++-----------
 2 files changed, 18 insertions(+), 11 deletions(-)

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 97eef89d..e7b3e051 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
@@ -53,13 +53,13 @@ subroutine overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,&
   integer                        :: iorder_p(3)
 
   call give_explicit_poly_and_gaussian(P_new,P_center,p,fact_p,iorder_p,alpha,beta,power_A,power_B,A_center,B_center,dim)
-!  if(fact_p.lt.1d-20)then
-!    overlap_x = 0.d0
-!    overlap_y = 0.d0
-!    overlap_z = 0.d0
-!    overlap = 0.d0
-!    return
-!  endif
+   if(fact_p.lt.1d-20)then
+     overlap_x = 0.d0
+     overlap_y = 0.d0
+     overlap_z = 0.d0
+     overlap = 0.d0
+     return
+   endif
   integer                        :: nmax
   double precision               :: F_integral
   nmax = maxval(iorder_p)