Cauchy-Schwarz inequality for integral derivatives
This is indeed a way that pre-screening of two-electron integrals is done, and if you thought of this yourself, then I'm very impressed! Though I doubt that is too much of a compliment.

Unfortunately, Tyberius is correct that you cannot be sure this inequality will hold after taking the derivative. Fortunately, the Cauchy-Schwartz inequality holds for any arbitrary vectors which have an inner-product. This means that you can simply switch the order you are doing things and it will work.

For what comes below, I will basically replicate some work from a paper by Reinhart Alrichs et al., and will use their notation as well [1].

In what follows, a superscipt $i$ represents a derivative in any of the cartesian $x$, $y$, or $z$ directions.

The Cauchy-Schwarz inequality states that,
$$\tag{1}\label{eq:1}
|(ab|cd)|\le Q_{ab}Q_{cd}
$$
where,
$$\tag{2}\label{eq:2}
Q_{ab}=(ab|ab)^{\frac{1}{2}}
$$
and the same for $Q_{cd}$.

If we take the derivative in the $i$ direction of $(\ref{eq:1})$, we have,
$$\tag{3}\label{eq:3}
(ab|cd)^i=(a^ib|cd)+(ab^i|cd)+(ab|c^id)+(ab|cd^i)
$$
which is contained in the work you have done, but with some extra stuff from the original inequality.

Each of these terms, however, are still a part of the original space and hence still have a well-defined inner-product. This means we can apply the Cauchy-Schwarz inequality directly to these derivatives as a means of screening the integrals.

That is,
$$\tag{4}\label{eq:4}
|(a^ib|cd)|\le Q_{a^ib}Q_{cd}
$$
and so one for the other derivatives.

The advantage of this is that one can give an upper bound for the size of these integrals based on calculations which have to be done for the other integrals anyways, and some simple to compute auxiliary quantities. See [1] for details on these auxiliary quantities (eqs. 10-19). I do not fully understand this part of the procedure, so I'm not restating it here.

Next, you asked where one should do the screening. The authors of [1] advocate what they call a "global screening" methodology. This means that all of the quantities required for calculating the upper-bound are pre-computed (before any integral loops but after any checks for symmetry), and compared to some pre-specified threshold, say $10^{-10}$. Although we do have to compute some new things, these are much cheaper than the integrals themselves and thus save quite a bit of time.

Ref. [1] also goes into detail about some test cases as well as doing the same screening for second-derivatives and hence the application of this pre-screening to Hessians.



A quite different approach which is very interesting to read is given in ref. [2], which came after and improves on [1]. This method, quite naturally I think, takes advantage of a generalization of Cauchy-Schwarz known as Holder's Inequality. For our purposes, the relevant form of Holder's inequality is (see the section titled Lebesgue Measure or ref. [2]),
$$\tag{5}
\int|f(x)g(x)|dx\le\left(\int|f(x)|^mdx\right)^{\frac{1}{m}} \left(\int|g(x)|^n dx\right)^{\frac{1}{n}}
$$
for $\frac{1}{m}+\frac{1}{n}=1$.

Notice, however, that while this strengthens the inequality by giving a class of possibilities for different $m$ and $n$, it weakens the inequality in that we may only use one variable $x$. Thus, to apply this to 2-electron integrals, the authors of ref. [2] first integrate out one of the variables, which gives us the same inequality for our physical situation as,
$$\tag{6}
\int|(ab|cd)|\le\left(\int|V_{ab}|^md\textbf{r}\right)^{\frac{1}{m}} \left(\int|(cd|cd)|^n d\textbf{r}\right)^{\frac{1}{n}}
$$
where $V_{ab}$ is the potential due to electron 1. Thus, in this case the screening must always take place after the first coordinate is integrated-out. This method, however, is very flexible in that it apparently gives tighter bounds (a smaller maximum), and also has many possible solutions. Ref. [2] discusses four of these possible solutions, and compares to the Schwarz pre-screening method, which is a special case of what they derive. Again, some special quantities must be pre-computed, but this is not such a big deal for the time it saves.

Ref. [2] does not discuss derivatives, though I assume the same method ought to possible as in the first part, only things might be complicated by the fact we need we need to integrate out a coordinate. Maybe not in the first derivative case, but in the second derivative case, this would present a problem as far as I can see.



[1] Horn, H., Weiß, H., Haser, M., Ehrig, M., &amp; Ahlrichs, R. (1991). Prescreening of two‐electron integral derivatives in SCF gradient and Hessian calculations. Journal of computational chemistry, 12(9), 1058-1064.

[2] Gill, P. M., Johnson, B. G., &amp; Pople, J. A. (1994). A simple yet powerful upper bound for Coulomb integrals. Chemical physics letters, 217(1-2), 65-68.