anchor.diff:                  ao_expo_transp(p,i),ao_expo_transp(q,j),ao_expo_transp(p,i),ao_expo_transp(q,j),&
                   I_power(1),J_power(1),I_power(1),J_power(1),         &
                   I_power(2),J_power(2),I_power(2),J_power(2),         &
  -                I_power(3),J_power(3),I_power(3),J_power(3))
  +                I_power(3),J_power(3),I_power(3),J_power(3))*coef2*coef2
  +        if (schwartz_kl(0,0)*schwartz_ij < thresh) then
  +           cycle
  +        endif
           do r = 1, ao_prim_num(k)
  +          if (schwartz_kl(0,r)*schwartz_ij < thresh) then
  +             cycle
  +          endif
             coef3 = coef2*ao_coef_transp(r,k)
             do s = 1, ao_prim_num(l)
  -            coef4 = coef3*ao_coef_transp(s,l)
  -            if (schwartz_kl(r,s)*schwartz_ij*coef4*coef4 < thresh) then
  +            if (schwartz_kl(s,r)*schwartz_ij < thresh) then
                  cycle
               endif
  +            coef4 = coef3*ao_coef_transp(s,l)
               integral = ERI(                                          &
                   ao_expo_transp(p,i),ao_expo_transp(q,j),ao_expo_transp(r,k),ao_expo_transp(s,l),&
                   I_power(1),J_power(1),K_power(1),L_power(1),         &
  
positive.code: integer function n_pt_sup(a_x,b_x,c_x,d_x,a_y,b_y,c_y,d_y,a_z,b_z,c_z,d_z)
    implicit none
    BEGIN_DOC
    ! Returns the upper boundary of the degree of the polynom involved in the
    ! bielctronic integral :
    !       Ix(a_x,b_x,c_x,d_x) * Iy(a_y,b_y,c_y,d_y) * Iz(a_z,b_z,c_z,d_z)
    END_DOC
    integer                        :: a_x,b_x,c_x,d_x,a_y,b_y,c_y,d_y,a_z,b_z,c_z,d_z
    n_pt_sup =  2 * ( (a_x+b_x+c_x+d_x) + (a_y+b_y+c_y+d_y) + (a_z+b_z+c_z+d_z) )
  end
  
  
  
negative.code:     enddo
      double  precision              :: ERI
  
      schwartz_kl(0,0) = 0.d0
      do r = 1, ao_prim_num(k)
        coef1 = ao_coef_transp(r,k)*ao_coef_transp(r,k)
        schwartz_kl(0,r) = 0.d0
        do s = 1, ao_prim_num(l)
          coef2 = coef1*ao_coef_transp(s,l)*ao_coef_transp(s,l)
          schwartz_kl(s,r) = ERI(                                      &
              ao_expo_transp(r,k),ao_expo_transp(s,l),ao_expo_transp(r,k),ao_expo_transp(s,l),&
              K_power(1),L_power(1),K_power(1),L_power(1),             &
              K_power(2),L_power(2),K_power(2),L_power(2),             &
              K_power(3),L_power(3),K_power(3),L_power(3)) * &
              coef2
          schwartz_kl(0,r) = max(schwartz_kl(0,r),schwartz_kl(s,r))
        enddo
        schwartz_kl(0,0) = max(schwartz_kl(0,r),schwartz_kl(0,0))
      enddo
  
      do p = 1, ao_prim_num(i)
  