anchor.diff:              ao_expo_transp(p,i),ao_expo_transp(q,j),                 &
               I_power,J_power,I_center,J_center,dim1)
           p_inv = 1.d0/pp
  -        schwartz_ij = general_primitive_integral(dim1,              &
  -                P_new,P_center,fact_p,pp,p_inv,iorder_p,             &
  -                P_new,P_center,fact_p,pp,p_inv,iorder_p)
  -        schwartz_ij = schwartz_ij
  +        schwartz_ij = general_primitive_integral(dim1,               &
  +            P_new,P_center,fact_p,pp,p_inv,iorder_p,                 &
  +            P_new,P_center,fact_p,pp,p_inv,iorder_p) *               &
  +            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
             double precision               :: coef3
             coef3 = coef2*ao_coef_transp(r,k)
             do s = 1, ao_prim_num(l)
               double precision               :: coef4
  -            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)
               double precision               :: general_primitive_integral
               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),             &
  
positive.code:                 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))
          do r = 1, ao_prim_num(k)
            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
                 cycle
              endif
              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),         &
  
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)
  