Computing the kinetic energy of a molecular orbital

Question

I want to compute the kinetic energy of the molecular orbitals computed by Gaussian 09, the variable U in this paper. Can it be computed from the orbital energy using a formula or does it need to be computed as part of the computation? How to do that in Gaussian 09?

Answer 1

Looking up the above paper,

The value of the kinetic energy U for each orbital in the initial state (usually the ground state) of the target is a theoretical quantity evaluated in any atomic or molecular wave function code that calculates the total energy.

Given the LCAO expansion $\psi_{i} ({\bf r}) = C_{\alpha i} \chi_\alpha({\bf r})$, what you therefore need to calculate is $U_i = \text{Tr } {\bf C}_i^\text{T} {\bf T} \bf{C}_i$, where ${\bf T}$ is the kinetic energy matrix with elements $T_{\alpha \beta} = \langle \chi_\alpha | -\frac 1 2 \nabla^2 | \chi_\beta \rangle$.

As far as I know, Gaussian does not compute these elements individually, only their sum, which is the total kinetic energy. To compute these variables, you need access to the kinetic energy and orbital coefficients matrices, and I am not sure how easy it is to access these, especially if you don't have access to the program's sources.

However, as I answered elsewhere just the day before, these quantities are easily accessible in the open-source Python program PySCF.

For example, with the code given below, I get the following orbital kinetic energies for N2:

Orbital 1 kinetic energy 601.73 eV
Orbital 2 kinetic energy 602.64 eV
Orbital 3 kinetic energy 70.41 eV
Orbital 4 kinetic energy 63.42 eV
Orbital 5 kinetic energy 54.89 eV
Orbital 6 kinetic energy 44.08 eV
Orbital 7 kinetic energy 44.08 eV

The last five values can be compared to 71.13 eV, 63.18 eV, 54.91 eV, and 44.30 eV given in the above paper; I do not know what bond length they used for the molecule.

The code to run the calculation is here

import pyscf

mol = pyscf.M(
    atom = 'N 0 0 0; N 0 0 1.098',
    basis = '6-311G',
    symmetry = True,
)

mf = mol.HF()
mf.kernel()

occs = mf.mo_occ
C = mf.mo_coeff
T = mol.intor_symmetric('int1e_kin')
for icol, occ in enumerate(occs):
    if occ>0.0:
        Ekin = C[:,icol].T @ T @ C[:,icol]
        print('Orbital {} kinetic energy {:.2f} eV'.format(icol+1,Ekin*27.2114))