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?
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 :)
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.
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]
