Predicting linear dependencies in a basis set in advance of calculating any integrals

Question

I've been doing some calculations for a water molecule, with a big uncontracted basis set: For oxygen I am using aug-cc-pV9Z supplemented with the "tight" functions from cc-pCV7Z. The s-type exponents are listed below, with all numbers up to 0.04456 constituting the s-type exponents of aug-cc-pV9Z, and the rest of the numbers constituting the "tight" s-type exponents of the CV7Z correction:

14977011.0 2218105.60 497972.050 136123.290 42655.7170 15004.6890 5776.15000 2375.75410 1017.44250 448.248580 203.162260 94.8087090 45.4553660 22.3219240 11.1990240 5.70236100 2.88556500 1.45704800 0.72489200 0.36172700 0.18245000 0.90164000E-01 0.04456 496.301272333 283.450755017 161.886207025 92.4574853342 52.8049100131 30.1582777362

After several hours of calculating integrals, I was told by the program that there were two near-linear-dependencies (determined by noticing that two eigenvalues of the overlap matrix were smaller than the default tolerance). The corresponding eigenvectors were removed automatically by the program, and my Hartree-Fock energy was higher than what it was, for the much smaller: aug-cc-pV9Z basis set (without the extra "tight" functions).

I decided to try to manually remove two basis functions to avoid any issue of near-linear-dependencies. But which ones should be removed?

I could plot all of the functions and see which ones look the most similar, but it would be tedious, and I do not know of any measure to determine which pairs of functions are most similar for this purpose, so I would merely be attempting to eye-ball it.

I suppose I could remove all of the exponents except for the ones that I think might lead to linear dependencies (based on the exponents being similar in magnitude), and then calculate the overlap matrix (which would be very fast since I would only have 4 functions in my basis set), and see if the small eigenvalues (near-linear-dependencies) still appear?

The two exponents most similar to each other are unsurprisingly, the two most diffuse functions:

0.90164000 0.04456 

but they surely won't lead to any problems, since they exist in the uncontracted aug-cc-pV9Z basis set, which does not lead to any issues.

So next I looked at the exponents that were most similar to each other percentage-wise:

94.8087090 92.4574853342

and it sure turned out that this guess was correct! The overlap matrix after removing one of them, now only had one eigenvalue that was too small (instead of two!). I then guessed to remove a function from the following second pair, which was the next closest to each other percentage-wise:

45.4553660 52.8049100131

and miraculously, the overlap matrix now had no eigenvalues below the tolerance (i.e. no near-linear-dependencies!), and the Hartree-Fock energy was lower than what I obtained with a vanilla uncontracted aug-cc-pV9Z (as expected).

Is it always safe to just look for the N pairs of exponents that are most similar to each other percentage-wise, and remove one from each pair, to cure N overly low eigenvalues? If so, why is it that none of the mainstream electronic structure packages have been able to implement an a priori test to "predict" near-linear-dependencies before spending several hours doing the integrals? I suppose when the geometry gets much more complicated, so will the procedure for predicting near-linear-dependencies in advance, but for diatomics and triatomics like water, is there ever a case in which the guessing procedure I used here fails?

Answer 1

You might be interested to know that I've recently presented a general solution to this problem in Curing basis set overcompleteness with pivoted Cholesky decompositions, J. Chem. Phys. 151, 241102 (2019).

The method is amazingly versatile, it also works if you have nuclei that are "unphysically" close to each other as I showed in Accurate reproduction of strongly repulsive interatomic potentials, Phys. Rev. A 101, 032504 (2020)

As pointed out in 1, the method can be implemented in two ways: you can either modify the usual orthonormalization routine, in which case your program might still waste a lot of time computing integrals that don't even appear in the calculation, or you can use the method to generate a customized basis set for the system you want to study by completely removing all shells that don't appear in the calculation.

As in the note in 2, the only thing you need to implement this method is the overlap matrix, which is very cheap to compute. Implementations of the method are presently available in my ERKALE code, as well as Psi4 and PySCF.