anchor.diff:              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)
  -            if (abs(fact_p*fact_q) < 1.d-8) then
  -              exit
  -            endif
               q_inv = 1.d0/qq
               integral = general_primitive_integral(dim1,              &
                   P_new,P_center,fact_p,pp,p_inv,iorder_p,             &
                   Q_new,Q_center,fact_q,qq,q_inv,iorder_q)
               ao_bielec_integral +=  coef4 * integral
  -!            if (abs(integral) < thresh) then
  -!              exit
  -!            endif
  -            do_r = .True.
  -            do_q = .True.
  -            do_p = .True.
             enddo ! s
           enddo  ! r
         enddo   ! q
  
positive.code:       L_power(p) = ao_power(l,p)
      enddo
      double  precision              :: ERI
      do_p = .True.
      do p = 1, ao_prim_num(i)
        if (.not.do_p) exit
        do_p = .False.
        do_q = .True.
        coef1 = ao_coef_transp(p,i)
        do q = 1, ao_prim_num(j)
          if (.not.do_q) exit
          do_q = .False.
          do_r = .True.
          coef2 = coef1*ao_coef_transp(q,j)
          do r = 1, ao_prim_num(k)
            if (.not.do_r) exit
            do_r = .False.
            coef3 = coef2*ao_coef_transp(r,k)
            do s = 1, ao_prim_num(l)
              coef4 = coef3*ao_coef_transp(s,l)
  
negative.code:   num_k = ao_nucl(k)
    num_l = ao_nucl(l)
    ao_bielec_integral = 0.d0
    
    if (num_i /= num_j .or. num_k /= num_l .or. num_j /= num_k)then
      do p = 1, 3
  