diff --git a/src/Pseudo_integrals/int.f90 b/src/Pseudo_integrals/int.f90
index 46e204eb..f7b35817 100644
--- a/src/Pseudo_integrals/int.f90
+++ b/src/Pseudo_integrals/int.f90
@@ -195,7 +195,7 @@ double precision, intent(in) :: v_kl(kmax,0:lmax),dz_kl(kmax,0:lmax)
 
 double precision :: fourpi,f,prod,prodp,binom,accu,bigR,bigI,ylm
 double precision :: theta_AC0,phi_AC0,theta_BC0,phi_BC0,ac,bc,big
-double precision :: areal,freal,breal,t1,t2,int_prod_bessel, int_prod_bessel_num_soph_p
+double precision :: areal,freal,breal,t1,t2,int_prod_bessel
 double precision :: arg
 
 integer :: ntot,ntotA,m,mu,mup,k1,k2,k3,ntotB,k1p,k2p,k3p,lambda,lambdap,ktot
@@ -276,7 +276,7 @@ if(ac.eq.0.d0.and.bc.eq.0.d0)then
     prod=bigI(0,0,l,m,n_a(1),n_a(2),n_a(3))
     prodp=bigI(0,0,l,m,n_b(1),n_b(2),n_b(3))
 
-      accu=accu+prod*prodp*v_kl(k,l)*int_prod_bessel_num_soph_p(ktot+2,g_a+g_b+dz_kl(k,l),0,0,areal,breal,arg)
+      accu=accu+prod*prodp*v_kl(k,l)*int_prod_bessel(ktot+2,g_a+g_b+dz_kl(k,l),0,0,areal,breal,arg)
 
    enddo
   enddo
@@ -309,7 +309,7 @@ else if(ac.ne.0.d0.and.bc.ne.0.d0)then
    do lambdap=0,lmax+ntotB
         do k=1,kmax
           do l=0,lmax
-          array_R(ktot,k,l,lambda,lambdap)= int_prod_bessel_num_soph_p(ktot+2,g_a+g_b+dz_kl(k,l),lambda,lambdap,areal,breal,arg)
+          array_R(ktot,k,l,lambda,lambdap)= int_prod_bessel(ktot+2,g_a+g_b+dz_kl(k,l),lambda,lambdap,areal,breal,arg)
           enddo
         enddo
    enddo
@@ -431,7 +431,7 @@ else if(ac.eq.0.d0.and.bc.ne.0.d0)then
     do k=1,kmax
       do l=0,lmax
 
-        array_R(ktot,k,l,0,lambdap)=  int_prod_bessel_num_soph_p(ktot+2,g_a+g_b+dz_kl(k,l),0,lambdap,areal,breal,arg)
+        array_R(ktot,k,l,0,lambdap)=  int_prod_bessel(ktot+2,g_a+g_b+dz_kl(k,l),0,lambdap,areal,breal,arg)
        enddo
     enddo
   enddo
@@ -517,7 +517,7 @@ else if(ac.ne.0.d0.and.bc.eq.0.d0)then
     do k=1,kmax
       do l=0,lmax
 
-      array_R(ktot,k,l,lambda,0)= int_prod_bessel_num_soph_p(ktot+2,g_a+g_b+dz_kl(k,l),lambda,0,areal,breal,arg)
+      array_R(ktot,k,l,lambda,0)= int_prod_bessel(ktot+2,g_a+g_b+dz_kl(k,l),lambda,0,areal,breal,arg)
       enddo
     enddo
   enddo
@@ -1898,14 +1898,32 @@ end
 !!
 !!    M_n(x) modified spherical bessel function 
 !!
-      double precision function int_prod_bessel(l,gam,n,m,a,b)
+      double precision function int_prod_bessel(l,gam,n,m,a,b,arg)
       implicit none
       integer n,k,m,q,l,kcp
       double precision gam,dblefact,fact,pi,a,b
-      double precision int,intold,sum,coef_nk,crochet
+      double precision int,intold,sum,coef_nk,crochet,u,int_prod_bessel_large,freal,arg
       logical done
 
+      u=(a+b)/(2.d0*dsqrt(gam))
+      freal=dexp(-arg)
+
+      if(a.eq.0.d0.and.b.eq.0.d0)then
+       if(n.ne.0.or.m.ne.0)then
+        int_prod_bessel=0.d0
+        return
+       endif
+       int_prod_bessel=crochet(l,gam)*freal
+       return
+      endif
+
+      if(u.gt.6.d0)then
+       int_prod_bessel=int_prod_bessel_large(l,gam,n,m,a,b,arg)
+      return
+      endif
+
       if(a.ne.0.d0.and.b.ne.0.d0)then
+
        q=0
        intold=-1.d0
        int=0.d0
@@ -1925,17 +1943,12 @@ end
          intold=int
         endif
        enddo
-       int_prod_bessel=int
-       if(kcp.gt.100) then
-          print*,"l",l
-          print*, "gam", gam
-          print*, "n", n
-          print*, "m", m
-          print*, "a", a
-          print*, "b", b
-          print*, "kcp", kcp
-          print*,'**WARNING** bad convergence in int_prod_bessel'
-        endif
+       int_prod_bessel=int*freal
+       if(kcp.gt.100)then
+        print*,'**WARNING** bad convergence in int_prod_bessel'
+!        else
+!        print*,'kcp=',kcp
+       endif
        return
       endif
 
@@ -1944,6 +1957,7 @@ end
         int_prod_bessel=0.d0
         return
        endif
+
        q=0
        intold=-1.d0
        int=0.d0
@@ -1959,7 +1973,7 @@ end
          intold=int
         endif
        enddo
-       int_prod_bessel=int
+       int_prod_bessel=int*freal
        if(kcp.gt.100)stop '**WARNING** bad convergence in int_prod_bessel'
        return
       endif
@@ -1969,6 +1983,7 @@ end
         int_prod_bessel=0.d0
         return
        endif
+
        q=0
        intold=-1.d0
        int=0.d0
@@ -1984,76 +1999,74 @@ end
          intold=int
         endif
        enddo
-       int_prod_bessel=int
+       int_prod_bessel=int*freal
        if(kcp.gt.100)stop '**WARNING** bad convergence in int_prod_bessel'
        return
-       endif
-
-      if(a.eq.0.d0.and.b.eq.0.d0)then
-       if(n.ne.0.or.m.ne.0)then
-        int_prod_bessel=0.d0
-        return
-       endif
-       int_prod_bessel=crochet(l,gam)
-       return
       endif
 
       stop 'pb in int_prod_bessel!!'
       end
 
-
-double precision function int_prod_bessel_num_soph_p(l,gam,n,m,a,b,arg)
+double precision function int_prod_bessel_large(l,gam,n,m,a,b,arg)
       implicit none
-      integer :: n,m,l
-      double precision :: gam,a,b,arg,arg_new
-      double precision :: bessel_mod,factor
-      logical :: not_done
-      double precision :: bigA,xold,x,dx,accu,intnew,intold,intold2,u,v,freal
-      integer :: iter, i, nI, n0
-      double precision :: eps
+      integer n,m,i,npts,l
+      double precision gam,a,b
+      double precision sum,x,bessel_mod,u,factor,arg
+      double precision xq(100),wq(100)
 
       u=(a+b)/(2.d0*dsqrt(gam))
-      arg_new=u**2-arg
-      freal=dexp(arg_new)
-      v=u/dsqrt(gam)
+      factor=dexp(u**2-arg)/dsqrt(gam)
 
-  bigA=v+dsqrt(-dlog(1.d-15)/gam)
-  n0=5
-  accu=0.d0
-  dx=bigA/(float(n0)-1.d0)
-  iter=0
-  do i=1,n0
-   x=(float(i)-1.d0)*dx
-   accu=accu+x**l*dexp(-gam*(x-v)**2)*bessel_mod(a*x,n)*bessel_mod(b*x,m)*dexp(-(a+b)*x)
-  enddo
+xq(1)= 5.38748089001123    
+xq(2)= 4.60368244955074   
+xq(3)= 3.94476404011563  
+xq(4)= 3.34785456738322 
+xq(5)= 2.78880605842813
+xq(6)= 2.25497400208928       
+xq(7)= 1.73853771211659      
+xq(8)= 1.23407621539532     
+xq(9)= 0.737473728545394   
+xq(10)= 0.245340708300901 
+xq(11)=-0.245340708300901
+xq(12)=-0.737473728545394    
+xq(13)=-1.23407621539532    
+xq(14)=-1.73853771211659   
+xq(15)=-2.25497400208928  
+xq(16)=-2.78880605842813 
+xq(17)=-3.34785456738322     
+xq(18)=-3.94476404011563    
+xq(19)=-4.60368244955074   
+xq(20)=-5.38748089001123  
+wq(1)=  2.229393645534151E-013
+wq(2)=  4.399340992273176E-010
+wq(3)=  1.086069370769280E-007
+wq(4)=  7.802556478532063E-006
+wq(5)=  2.283386360163528E-004
+wq(6)=  3.243773342237853E-003
+wq(7)=  2.481052088746362E-002
+wq(8)=  0.109017206020022
+wq(9)=  0.286675505362834
+wq(10)= 0.462243669600610
+wq(11)= 0.462243669600610
+wq(12)= 0.286675505362834
+wq(13)= 0.109017206020022
+wq(14)= 2.481052088746362E-002
+wq(15)= 3.243773342237853E-003
+wq(16)= 2.283386360163528E-004
+wq(17)= 7.802556478532063E-006
+wq(18)= 1.086069370769280E-007
+wq(19)= 4.399340992273176E-010
+wq(20)= 2.229393645534151E-013
 
-accu=accu*freal
-intold=accu*dx
+      npts=20
+!      call gauher(xq,wq,npts)
 
-eps=1.d-08
-nI=n0-1
-dx=dx/2.d0
-not_done=.true.
-
-do while(not_done)
- iter=iter+1
- accu=0.d0
- do i=1,nI
-  x=dx+(float(i)-1.d0)*2.d0*dx
-  accu=accu+dx*x**l*dexp(-gam*(x-v)**2)*bessel_mod(a*x,n)*bessel_mod(b*x,m)*dexp(-(a+b)*x)
- enddo
- accu=accu*freal
- intnew=intold/2.d0+accu
- if(iter.gt.1.and.dabs(intnew-intold).lt.eps.and.dabs(intnew-intold2).lt.eps)then
-  not_done=.false.
- else
-  intold2=intold
-  intold=intnew
-  dx=dx/2.d0
-  nI=2*nI
- endif
-enddo
-int_prod_bessel_num_soph_p=intold
+      sum=0.d0
+      do i=1,npts
+       x=(xq(i)+u)/dsqrt(gam)
+       sum=sum+wq(i)*x**l*bessel_mod(a*x,n)*bessel_mod(b*x,m)*dexp(-(a+b)*x)
+      enddo
+      int_prod_bessel_large=sum*factor
 end
 
 !!    Calculation of
diff --git a/tests/unit_test/unit_test.py b/tests/unit_test/unit_test.py
index 7ce1c0f6..b2881c92 100755
--- a/tests/unit_test/unit_test.py
+++ b/tests/unit_test/unit_test.py
@@ -169,7 +169,7 @@ def run_hf(geo, basis, mult=1, pseudo=False, remove_after_sucess=True):
 
     ref_energy["sto-3g"]["methane"] = Energy(-39.7267433402, None)
     ref_energy["vdz"]["SO2"] = Energy(None, -41.48912297776174)
-    ref_energy["vdz"]["HBO"] = Energy(None, -19.11982537379352)
+    ref_energy["vdz"]["HBO"] = Energy(None, -19.1198251160)
 
     # ~#~#~#~#~#~#~#~#~#~#~#~#~#~#~ #
     # G l o b a l _ v a r i a b l e #