Why does bond length decrease with increased basis set quality for SCF methods?

Question

From Reviews in Computational Chemistry Vol. 5, for $\ce{H2O}$ using HF, MP2, CCSD, and CCSD(T), all lengths in angstroms:

$$ \small \begin{array}{rcccc} \text{Basis} & \text{HF} & \text{MP2} & \text{CCSD} & \text{CCSD(T)} \\ \hline \text{STO-3G} & 0.989 & 1.013 & 1.028 & 1.028 \\ \text{3-21G} & 0.967 & 0.989 & 0.993 & 0.994 \\ \text{DZ} & 0.951 & 0.979 & 0.979 & 0.980 \\ \text{6-31G*} & 0.947 & 0.969 & 0.969 & 0.971 \\ \text{DZP} & 0.944 & 0.963 & 0.961 & 0.962 \\ \text{TZ2P} & 0.941 & 0.958 & 0.956 & 0.959 \\ \text{Experimental} & 0.957 \end{array} $$

In case one is skeptical due to effects from the lone pairs and the changing bond angle, here are results for $\ce{BH}$ using Hartree-Fock:

$$ \small \begin{array}{rc} \text{Basis} & r_{e} \\ \hline \text{cc-pVDZ} & 1.235950 \\ \text{cc-pVTZ} & 1.222062 \\ \text{cc-pVQZ} & 1.220306 \\ \text{cc-pV5Z} & 1.220010 \\ \text{cc-pV6Z} & 1.219877 \end{array} $$

and for $\ce{LiH}$ using Hartree-Fock:

$$ \small \begin{array}{rc} \text{Basis} & r_{e} \\ \hline \text{cc-pVDZ} & 1.619017 \\ \text{cc-pVTZ} & 1.607865 \\ \text{cc-pVQZ} & 1.605530 \\ \text{cc-pV5Z} & 1.605519 \end{array} $$

Going on the NIST CCCBDB and looking at calculated geometries, I cannot find a counterexample, at least when considering only HF, as adding dynamic correlation introduces an effect to counterbalance this bond shortening, since electrons can then properly repel each other to levels beyond the mean-field approximation, and DFT offers no insight beyond HF here.

One potential reason that doesn't receive much attention is intramolecular BSSE, where one part of a molecule borrows the basis functions from another. However, in the case of my calculations, I would argue that only cc-pVDZ could be considered a "deficient" basis, and the trend is still seen when leaving it out. I also believe that this is more plausible for molecules with distinct "domains" that can interact with each other, like one amino acid in a peptide forming a hydrogen bond with another AA, etc.

The current rationale I have is that as the basis set becomes more complete, and functions with higher angular momentum (polarization functions) are included, the overlap between the atomic wave functions is increased. This increased overlap results in higher electron density in the bonding region, which strengthens it, causing it to become shorter, like going from ethane to ethene to ethylene.

Specifically, I'm looking for a literature reference where an explanation is given for the phenomenon, not just examples with or without comparison to experiment.

Answer 1

I think that your rationale is right. I would think it in terms of energy instead to make a little more precise. Thinking about a diatomic molecule, the nuclei-nuclei repulsion function of distance does not change with the basis set used. So the difference is strictly in terms of electronic part. At large distance the atoms behave more similar to an isolated atom, so few spherical harmonics are enough (for a relative description). As atoms get closer it is needed more flexibility to represent the very different probability density that there is in the space region located between atoms. With a smaller basis set this wave function is approximated poorer. If so, it is expected that the energy rises in such region if one consider the variational principle.

Here I calculated for the oxygen molecule in the region close to equilibrium distance (about 1.2 A) the difference in electronic energy between HF/STO-3G and HF/6-31* basis set:

The behavior close to 1.4A should be just an artifact in convergence (I did not take care about that).

Sadly I never found a discussion about that in the literature. Up to my understanding there is not to much that can be said about that.