anchor.diff:  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
  +  ! Returns the upper boundary of the degree of the polynomial 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) )
  +  n_pt_sup =  ishft( a_x+b_x+c_x+d_x + a_y+b_y+c_y+d_y + a_z+b_z+c_z+d_z,1 )
   end
   
   
  
positive.code:     p = alpha + beta
      q = delta + gama
      ASSERT (p+q >= 0.d0)
      n_pt = 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)
      coeff = pi_5_2 / (p * q * dsqrt(p+q))
      
      if (n_pt == 0) then
        ERI = coeff
  
negative.code:       L_center(p) = nucl_coord(num_l,p)
      enddo
      
      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)
          call give_explicit_poly_and_gaussian(Q_new,Q_center,qq,fact_q,iorder_q,&
              ao_expo_transp(r,k),ao_expo_transp(s,l),                 &
              K_power,L_power,K_center,L_center,dim1)
          q_inv = 1.d0/qq
          schwartz_kl(s,r) = general_primitive_integral(dim1,          &
              Q_new,Q_center,fact_q,qq,q_inv,iorder_q,                 &
              Q_new,Q_center,fact_q,qq,q_inv,iorder_q)                 &
              * 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)
  