Does the Hartree Fock energy of a virtual orbital fulfill the virial theorem?

The virial theorem does not hold for orbital energies separately. Remember that

$$ 2\langle T\rangle=k\langle V\rangle $$

holds only for homogeneous degree $k$ potentials of form $V\sim x^k$ or $V\sim r^k$ ($k=-1$ for the Coulomb potential). However, the potential felt by the separate electrons in the HF scheme is absolutely not of this form (it is not even a local potential in general).

As a simple example, think of a He-like atom, whose occupied spatial orbital satisfies

$$ \left[ -\frac{\hbar^2}{2m}\nabla^2 -\frac{Z_{\text{eff}}(r)e^2}{4\pi\epsilon_0r} \right]\chi_{\text{HF}}(r)=\varepsilon_{\text{HF}}\chi_{\text{HF}}(r) \ , $$ where $$ Z_{\text{eff}}(r) = Z-r\int_0^{\infty}\mathrm{d}r'r'^2 \frac{|\chi_{\text{HF}}(r')|^2}{r_>} \ , $$ and $r_{>}=\text{max}(r,r')$. In this simple case of a two-electron singlet state, the HF potential is actually local (as there is no exchange interaction between different spins), but its form is not of $r^k$, so the virial theorem does not say anything useful here (you could show that $Z_{\text{eff}}(r)$ is roughly an exponential-like function interpolating between $Z_{\text{eff}}(0)=Z$ and $Z_{\text{eff}}(\infty)=Z-1$).

Some numerical results for He ($Z=2$): $$ E_{\text{HF}}\approx -2.862 \, E_h \ , $$ $$ \langle T\rangle_{\text{HF}}\approx 2.862 \, E_h \ , $$ $$ \varepsilon_{\text{HF}}\approx-0.918 \, E_h \ . $$ As you can see, the virial theorem holds for the total energy and kinetic energy, but $\varepsilon_{\text{HF}}\neq-\frac{1}{2}\langle T\rangle_{\text{HF}}$ (the kinetic energy of a single electron is half of the total kinetic energy, hence the factor of $1/2$). Based on the above, there is no reason to expect anything different for virtual orbitals either.

Also, please note that the HF energy is $not$ the sum of occupied orbital energies, as half of the electron-electron interaction must be subtracted in order to avoid double counting.

One more remark: your question seems to imply that (total) energies obtained from a variational procedure always satisfy the virial theorem, which, however, is not true. Such results typically hold for coordinate-scaling variational calculations, where the energy calculated with trial function $\Phi(\vec{r}_1,...,\vec{r}_N)$ is improved by the following scaling: $$ \Phi(\vec{r}_1,...,\vec{r}_N) \rightarrow \Phi_{\lambda}(\vec{r}_1,...,\vec{r}_N)= \Phi(\lambda\vec{r}_1,...,\lambda\vec{r}_N) \ \ \ \ (\lambda>0) \ , $$ $$ E(\lambda)=\frac{\langle\Phi_{\lambda}|H|\Phi_{\lambda}\rangle} {\langle\Phi_{\lambda}|\Phi_{\lambda}\rangle} \ , $$ $$ \frac{\partial}{\partial\lambda}E(\lambda)\Big|_{\lambda=\lambda_{\text{opt}}}=0 \ \ \ , \ \ \ E_{\text{opt}}=E(\lambda_{\text{opt}}) \ . $$ Other kinds of trial functions (like e.g. a linear variational Ansatz) in general do not satisfy the virial theorem. This is the reason why the "exact" HF energy (that is, HF in the complete basis set limit) satisfies the virial theorem, but the LCAO HF energy (calculated with a finite number of fixed basis functions) does not.

(For molecules in the Born-Oppenheimer approximation, the virial theorem only holds in the equilibrium geometry.)