Apparent inversion of oxidation state of sulfate ion using minimal basis set

Question

It is well known that the sulfate ion has essentially no 3d-2p bonding; hence I used the minimal basis set(STO-3G) for the NBO calculation I did via Gaussian, just to check it.

However, I got a surprising result- all four S-O bonds had higher coefficients on the sulfur than oxygen, which contradicts the basic fact that the oxidation state of sulfur in the sulfate ion is +VI; the four S-O bonds must all be "assigned to the oxygen", i.e. have higher coefficients on the oxygen.

At first I thought this was because I used the Hartree-Fock method, which is often too rough for precise calculations. So I tried using B3LYP and the density=current command. Same results still.

Now comes the main question- what would be the reason for this contradiction, and what methods can be used to solve this?

(All S-O bonds were set to be 150pm with perfect tetrahedral symmetry)

The basis set is too small, it is probably missing flexibility. — Martin - マーチン, Sep 09 '21 at 22:09
I thought that way until I realised that larger basis sets give sulfur an NBO charge of about +2.5, which is clearly less similar to the known empirical value(+1.7) than is the value from my calculation(+1.3). — Kanghun Kim, Sep 10 '21 at 02:36

score 4 · Accepted Answer · answered Sep 13 '21 at 19:27

The chosen structure and level of theory is not remotely close to a reasonable approximation.
You are calculating an anion with charge -2, but there are no polarisation or diffuse functions in STO-3G. Hartree Fock behaves quite poorly, especially for those systems.

This is quite obvious when you look at the bond lengths. For lack of time I'll refer to values from Wikipedia (sulfate): $d(\ce{S-O}) = \pu{149 pm}$. At the HF/STO-3G level of theory that value is $d_\mathrm{HF}(\ce{S-O}) = \pu{176 pm}$.

Running a bonding analysis on an approximate wave-function for a structure, which is not even close to experimental values will not produce any meaningful results.
Gaussian even warns you about this:

Warning!  S  atom    1 may be hypervalent but has no d functions.

So while there is little involvement of the d-orbitals, they still make a huge difference. Remember the correlation energy is only less than 1% of the total energy, but it is exactly that tiny bit you have to get right to make reasonable approximations.

With that in mind it is easy to spot that the NBO analysis is horrible, too. The charge on sulfur is way too low, as is the same for oxygen.

                                     Natural Population
             Natural    ---------------------------------------------
  Atom No    Charge        Core      Valence    Rydberg      Total
 --------------------------------------------------------------------
    S  1    0.50743     10.00000     5.49257    0.00000    15.49257
    O  2   -0.62686      2.00000     6.62686    0.00000     8.62686
    O  3   -0.62686      2.00000     6.62686    0.00000     8.62686
    O  4   -0.62686      2.00000     6.62686    0.00000     8.62686
    O  5   -0.62686      2.00000     6.62686    0.00000     8.62686
 ====================================================================
 * Total * -2.00000     18.00000    32.00000    0.00000    50.00000

If you force the structure, you'll get a wave-function, which won't tell you anything remotely approximate. So you might get all sorts of artifacts, even such that apparently lead to a contradiction.

For the above the bond orbital then is quite within reason, meaning that it is obviously also wrong:

  22. (1.97736) BD ( 1) S  1- O  2
               ( 68.31%)   0.8265* S  1 s( 25.00%)p 3.00( 75.00%)
                                         0.0000  0.0000  0.5000  0.0000  0.8660
                                         0.0000  0.0000  0.0000  0.0000
               ( 31.69%)   0.5629* O  2 s(  3.15%)p30.79( 96.85%)
                                         0.0000  0.1774 -0.9841  0.0000  0.0000

You can see how terrible this approximation is as oxygen uses an almost pure p-orbital for bonding.

If you increase the level of theory to DF-BP86/def2-SVP, you'll find $d_\mathrm{BP86}(\ce{S-O}) = \pu{154 pm}$, which is still not good, but at least closer.

The charges make a lot more sense:

                                     Natural Population
             Natural    ---------------------------------------------
  Atom No    Charge        Core      Valence    Rydberg      Total
 --------------------------------------------------------------------
    S  1    2.34734      9.99995     3.43215    0.22057    13.65266
    O  2   -1.08684      1.99997     7.07959    0.00727     9.08684
    O  3   -1.08683      1.99997     7.07959    0.00727     9.08683
    O  4   -1.08683      1.99997     7.07959    0.00727     9.08683
    O  5   -1.08683      1.99997     7.07959    0.00727     9.08683
 ====================================================================
 * Total * -2.00000     17.99984    31.75051    0.24965    50.00000

So does the bonding orbital:

  22. (1.97871) BD ( 1) S  1- O  2
               ( 34.07%)   0.5837* S  1 s( 25.00%)p 2.95( 73.81%)d 0.05(  1.19%)
                                         0.0000  0.0000  0.5000 -0.0018  0.0000
                                         0.8557  0.0772  0.0000  0.0000  0.0000
                                         0.0000  0.0000  0.0000  0.0000  0.0000
                                         0.0000  0.0944 -0.0545
               ( 65.93%)   0.8119* O  2 s( 23.20%)p 3.30( 76.64%)d 0.01(  0.16%)
                                         0.0000  0.4812  0.0210 -0.8754 -0.0102
                                         0.0000  0.0000  0.0000  0.0000  0.0000
                                         0.0000  0.0000  0.0348 -0.0201

With all this keep in mind: This is a calculation for an anion at virtually $\pu{0 K}$ in vacuum. There are no counter charges. That in itself is a poor approximation. You can easily see this in the eigenvalues for the molecular orbitals $\varepsilon_\mathrm{gas}(\mathrm{HOMO}) = \color{red}{+}\pu{7 eV}$.

This might be circumvented by using a solvent model (PCM, water), at least the energy is now more reasonable $\varepsilon(\mathrm{HOMO}) = \pu{-3 eV}$. The structure is virtually unchanged and only slightly shortened bonds.

As a general takeaway: Anions need at least polarisation functions, better are diffuse functions. Split valence basis sets are for heavy duty, work horse calculations only. For analysis and energies always use at least triple zeta basis sets. Always calibrate your DFA calculations (DFT Functional Selection Criteria). The smaller the molecule the better the level of theory needs to be.

"As a general takeaway: Anions need at least polarisation functions, better are diffuse functions. Split valence basis sets are for heavy duty, work horse calculations only." Does this hold for cations as well, may I ask? — Kanghun Kim, Sep 14 '21 at 00:09
@KanghunKim For cations you will not have the problem with the orbital energies. The rest applies to every calculation. Furthermore, split valence basis sets usually have polarisation functions already. Benchmarking is always the only way to get a reasonable approximation. — Martin - マーチン, Sep 14 '21 at 17:44

Apparent inversion of oxidation state of sulfate ion using minimal basis set

1 Answers1