# How closely should the energy components of ab initio SCF results at a stationary nuclear configuration preserve the virial theorem?

When running a SCF geometry optimization for water molecule, I found that the energy components of the optimized geometrical structure appear to deviate from Virial theorem sightly:

Kinetic energy = 75.424942151953331
Total energy = -75.585962823645531

How should this discrepancy be interpreted? Is it an indication of inaccuracy in the geometrical optimization algorithm of the software being used?

• Please see my edits. I think the title is now a better reflection of the question, but please feel free to edit or roll back if you don't like it. – hBy2Py Aug 23 '17 at 2:53
• Thanks for the edits, the new title is better. If I read correctly this book, Basic Principles and Techniques of Molecular Quantum Mechanics By Ralph E. Christoffersen, pp 460, books.google.com/… virial theorem applies to non-exact wave function as long as the geometry is optimized, independent of the choice of basis set. – wang1908 Aug 23 '17 at 21:51

I hope I am not mistaken about this, but if I'm wrong then I'll learn so no big deal.

I believe this is a result of the fact that the electrons do not experience a true $r^{-1}$ potential due to electron-electron repulsion. Namely, we should expect that the electrons will repel each other which will cause the apparent attraction to the nuclei to be weaker than $r^{-1}$.

In other words, for a general potential, $$V(r)=r^{n}$$ we have, $$\frac{dV(r)}{dr}=nr^{n-1}$$

Because the Virial theorem can be stated as, $$\langle T\rangle=\frac12\langle\vec{r}\cdot\nabla V\rangle$$ But for a potential which is less than coulombic, we have, $$\langle T\rangle=-\frac{n}{2}\langle V\rangle=-nE_{tot}$$

That is to say, the electrons in water are feeling an effective potential of the form $V(r)=r^{-.99787}$ or so.

Or, in the language of the book you linked, the trial wavefunction of the form being tried for the hydrogen atom does not evaluate to a parameter $\zeta=1$.

Again, I am not completely confident that I nailed this answer, but for systems that have no obvious analytic form of the potential energy, there is no reason I see that we have should have $\langle T\rangle=-E_{tot}$ be exactly true. That is not to say the virial theorem is not satisfied, it certainly is, only that it is not as clean as we are used to.

I was a bit loose with some constants in my equations, please forgive me :)

• Thanks for the insightful answer. Although SCF is single-body mean field theory, the wave function thus derived is an approximation to the exact wave function of the all-electron Hamiltonian, in which case n = 1 in your formula, even after electron-electron repulsion is taken into account. The above deviation is only significant when a small basis set, say 3-21G, is used. For 6-31G deviation is rather negligible. You are probably right to the point that for smaller basis set, it does not evaluate to a parameter ζ=1, but I do not fully understand how. – wang1908 Aug 29 '17 at 22:19
• Yes I think you're correct. It makes sense. I'm completely guessing here but it could be that if smaller basis sets show a larger deviation, and they should follow the virial theorem exactly, that this is the result of some kind of intramolecular basis set superposition error... ? – jheindel Aug 30 '17 at 0:08
• @Wang-X-Y Ehhh... I ran some calculations on water with a number of different basis sets and there is very little change between them. In fact sometimes the smaller basis set had a ratio closer to one, but really the difference is quite small. I think I'll post these sometime tomorrow because pretty much my answer above isn't right. I found it very compelling, but the fact HF is a mean-field theory invalidates my reasoning. – jheindel Aug 30 '17 at 5:22
• Thanks for looking into this, it's very nice to have someone to compare notes:-) Here is the input of the calculation, for the purpose of the project I worked on UHF rather than RHF chosen, not sure where it caused the differences. molecule h2o { O H 1 0.96 H 1 0.96 2 104.5 } set basis 3-21g set scf_type pk set g_convergence GAU_VERYTIGHT set reference uhf E, wfn = optimize('scf',return_wfn=True) – wang1908 Aug 30 '17 at 11:47

I am posting a second answer rather than editing my first answer basically just because my first answer is wrong, and this answer is just going to present some data and then make the opposite claim of my first answer.

I ran some calculations using Molpro on $\ce{H2O}$ using a number of different basis sets. Molpro, conveniently, actually prints out the Virial ratio, and I believe it prints this as $\frac{\langle T\rangle}{\langle V\rangle}$ (as opposed to the other way). So, after optimizing the water using each basis set at the HF level, the ratio was found. Note that HF was used because finding this ratio for post-HF methods is much more complicated, and I believe that in the case of post-HF methods, my first answer actually would be correct. Namely, the ratio should no longer be -1.

\begin{array}{cccc} \hline \text{Basis}& \text{avdz} & \text{avtz} & \text{avqz} & \text{6-31g} & \text{6-31g**} & \text{6-31g++} & \text{STO-3g} \\ \hline \mathrm{Ratio} & -1.00080 & -1.00051 & -1.00008 & -0.99961 & -1.0008 & -0.99919 & -1.00600 \\ \mathrm{Deviation} & 0.00080 & 0.00051 & 0.00008 & -0.00039 & 0.0008 & -0.00081 & 0.00600 \\ \hline \end{array} Note that in the above avxz (x=d,t,q) is the aug-cc-pvxz basis set.

So, it does seem that with the correlation-consistent basis sets, the deviation from a ratio of one decreases as the basis set becomes more complete. The Pople basis sets don't follow any trend that I can see. STO-3g is the worst which makes sense cause it's the worst basis set.

I would suspect that these differences are not significant. I think this actually is a reasonable method of determining how converged a HF geometry is, but I am suspicious that something else is going on here. I don't know why there would be deviations. The only thing I can think of is intramolecular BSSE, which is rarely talked about because it usually doesn't matter, but maybe this measure is sensitive to it though. I kind of don't believe that though so I don't know what's going on. I guess one thing to try is to look at this ratio for a system which is known to suffer from intramolecular BSSE.

The only other thing I can imagine is the boring answer (but is probably right) that converging the total energy and the gradient, which is what are usually converged in geometry optimizations, does not converge the kinetic energy beyond the third-ish decimal place, so we see deviations.

As to your question of how much tolerance we should have for this, my guess is that the amount of tolerance we have now is sufficient because otherwise people would have figured out what was wrong and fixed it. That is, there is no obvious error in QM calculations that can be traced back to disobeying the Virial theorem. I would guess the reason for this is that usually post-HF methods are used anyways, so any effect of slightly disobeying the Virial theorem gets washed out because the geometry at higher levels is different from the HF geometry. Whether or not this same deviation is seen at higher levels of theory is something I don't know and would be trickier to answer.

• The numbers of deviations for different basis sets in your calculation shows trend similar to what I saw but the slightly different numerical values. These several explanation can all be valid, just difficult to evaluate to what extend they contribute, The question can be asked differently: seeing a deviation from Virial Theorem in one's optimized wave function, is it enough reason to infer that the basis set used is too small? – wang1908 Sep 1 '17 at 1:51

Originally, I was convinced that the molecular geometry has very limited influence on the virial ratio: only in so far as the effective basis set will change with the geometry as long as nuclei-centered Slater or Gauss (STO or GTO) basis functions are used.1

Hartree-Fock is accurate as full CI (within the given basis set and its postulates and earlier approximations) for one-electron systems. Let's look at $\ce{H}$: $$\begin{array}{lr} \mathrm{Basis\ set} & |VR-2| \\ \hline \mathrm{cc-pVDZ} & 0.00002250 \\ \mathrm{cc-pVTZ} & 0.00004588 \\ \mathrm{cc-pVQZ} & 0.00000153 \\ \hline \mathrm{def2-SVP} & 0.00000000 \\ \mathrm{def2-TZVP} & 0.00000049 \\ \mathrm{def2-QZVP} & 0.00000404 \\ \hline \end{array}$$ $VR$ is the virial ratio, whose value I expect to be 2. At this point, we can conclude a) that the basis set has an influence and b) it's more complicated than simple basis set size, even in reasonably constructed basis sets/basis set families.

Let's look at $\ce{H2+}$ at different $\ce{H}-\ce{H}$ distances, indicated by the factor, i.e. as multiples of the equilibrium distance. The equilibrium geometry was found at the computational level: HF/aug-cc-pVQZ. The energy is given relative to the energy of the equilibrium structure. $$\begin{array}{lll} \mathrm{factor} &\mathrm{rel.\ }E &|VR-2|\\ \hline 0.850 &0.120850254 &-0.11563802 \\ 0.900 &0.048912238 &-0.07368802 \\ 0.925 &0.026293496 &-0.05396290 \\ 0.950 &0.011177260 &-0.03503389 \\ 0.975 &0.002673918 &-0.01687637 \\ 1.000 &0.000000000 &+0.00053214 \\ 1.025 &0.002463880 &+0.01721247 \\ 1.050 &0.009453714 &+0.03318366 \\ 1.075 &0.020427030 &+0.04846357 \\ 1.100 &0.034901950 &+0.06306886 \\ 1.150 &0.072687778 &+0.09031731 \\ 1.200 &0.119906546 &+0.11504445 \\ 1.300 &0.233820556 &+0.15734817 \\ 1.400 &0.363035528 &+0.19075111 \\ \hline \end{array}$$ We can conclude that at least for this system, the influence from the geometry is a lot larger than from the basis set.

In the comments, the suspicion that integral neglect thresholds and similar settings influence $VR$ was voiced. I decided to rerun the calculations from the first table using the tightest preset defined in ORCA, EXTREMESCF, which is described as close to machine accuracy. Given that some of the thresholds are then below $10^{-15}$ au, this seems credible. The result is: $$\begin{array}{lr} \mathrm{Basis\ set} & |VR-2| \\ \hline \mathrm{cc-pVDZ} & 0.00002250 \\ \mathrm{cc-pVTZ} & 0.00004539 \\ \mathrm{cc-pVQZ} & 0.00000156 \\ \hline \mathrm{def2-SVP} & 0.00000000 \\ \mathrm{def2-TZVP} & 0.00000000 \\ \mathrm{def2-QZVP} & 0.00000000 \\ \hline \end{array}$$ The result surprised me, there actually is discernible influence on $VR$ from the thresholds. What's more, basis set families are different in this regard.

I also reran calculations from the second table when setting one of the atoms atomic charge to zero - i.e. with a ghost basis. The energy of the "dissociation" curve is not relevant for this crude investigation, and is thus omitted. $$\begin{array}{ll} \mathrm{factor} &|VR-2|\\ \hline 0.850 &0.00027735 \\ 0.900 &0.00027606 \\ 0.925 &0.00027525 \\ 0.950 &0.00027437 \\ 0.975 &0.00027345 \\ 1.000 &0.00027252 \\ 1.025 &0.00027160 \\ 1.050 &0.00027073 \\ 1.075 &0.00026994 \\ 1.100 &0.00026924 \\ 1.150 &0.00026819 \\ 1.200 &0.00026767 \\ 1.300 &0.00026785 \\ 1.400 &0.00026865 \\ \hline \end{array}$$ So the intramolecular BSSE can probably neglected for these investigations.

Further work would include extending post-HF methods to allow for calculation of the virial ratio. Application to methods that are regularly used for geometry optimization (e.g. DFT) is not trivial because exact calculation of the virial ratio is hard.2 A cursory evaluation of multi-electron systems on the HF level may also be useful.

1 I currently only have GTO programs at my disposal.
2 The kinetic part of the correlation energy is subsumed into the energy calculated by the XC functional [Jensen, Introduction to Computational Chemistry, 2nd Ed., Wiley 2007, p. 236]

• The ratio of 2 is kept only by stationary geometrical state. For non-stationary there is another term: 2<T>+<V> = R*dE/dR, in case of diatomic molecules. – wang1908 Sep 1 '17 at 1:55
• Thanks for confirming the effect of basis set on the deviation, particularly for the case of H atom. This is very convincing I think. As @jheindel mentioned in his answer, this 'complicated' issue might be related by BSSE. Or more likely this could be somehow related to basis set incompleteness or BSIE. It would be interesting the deviation can be somehow formulated into a measure of BSIE. – wang1908 Sep 1 '17 at 2:14
• @Wang-X-Y I was thinking the same thing. Truth be told I did those calculations because I thought it would be interesting if this could be formulated as a measure of intramolecular BSSE, but I think it has more to do with basis set incompleteness, and, honestly, how stringently you converge integrals. I would guess the likeliest thing is that this is from accumulation of error due to neglecting integrals. Although it seems weird this would manifest itself in the kinetic energy more than the potential energy. – jheindel Sep 1 '17 at 2:25
• @TAR86 you're first table is quite nice because it shows the trend does not necessarily depend on basis set size. The second table, however, is a bit artificial because the nuclear configuration needs to be at a minimum for VR-2 to theoretically equal zero. As Wang-X-Y says, another term must be included. – jheindel Sep 1 '17 at 2:28
• Also, I've become pretty convinced now that the only way we're gonna answer this question is if we get rid of all assumptions and do the calculation involving all integrals to be sure this isn't from screening out integrals. – jheindel Sep 1 '17 at 2:30