How to get atomic orbital coefficients in PySCF

Question

I'm new to quantum chemistry software, so I'll try my best to explain my thought process. When running a DFT calculation, the software uses a linear combination of basis functions to approximate the electron density. The coefficients of the basis functions are optimized to minimize the energy. My question is: how do I get those coefficients out of PySCF?

I've provided some code below in case it is helpful to build off of:

from pyscf import gto, dft

# Create molecule and select basis set
mol = gto.Mole()
mol.atom = 'N 0 0 0; N 0 0 1.2'
mol.basis = 'STO-3G'
mol.build()

# Create DFT solver and xc functional
mf = dft.RKS(mol)
mf.xc = 'b3lyp'

# Run DFT calculation
mf.kernel()

In this example, I'm using the STO-3G basis set which requires 5 orbitals per nitrogen atom: 1s, 2s, and three 2p orbitals. Since there are 10 orbitals, I'm looking for the 10 values of the coefficients (several will be identical because of cylindrical symmetry).

Answer 1

Broadly, there are two flavors of Density Functional Theory (DFT): Orbital Free DFT and Kohn-Sham DFT (KS-DFT). In Orbital Free DFT, we do indeed directly optimize the density until we find the minimum energy. In KS-DFT, however, we instead optimize the coefficients of the KS orbitals, which are the solutions to the KS equations: $$\left(-\frac{1}{2}\nabla^2 + v(r)\right) \psi_i(r) = \epsilon_i \psi_i(r),\tag{1}$$ and which can be used to construct the density: $$\rho(r) = \sum_i^N |\psi_i(r)|^2.\tag{2}$$

The coefficients of the KS orbitals are readily available as a 2D array through mf.mo_coeff. If you read the quickstart guide, you can see them being manipulated in the section “Spatially Localized Molecular Orbitals”.

Now, as you mentioned we are using some finite basis for these calculations. In practice, these basis functions are used to represent the KS Orbitals, not the density: $$\psi_i(r) = \sum_j^M c_{ij} \phi_j(r).\tag{3}$$

You'll notice that if you try plugging this equation in to the equation for the density, you will end up squaring both the coefficients and the basis functions, so that it is not straightforward to compute the density in the original basis. Furthermore, just because the original basis can represent the KS orbitals well, there is no guarantee it can represent the density.

Nonetheless, we have to construct the density in order to compute the KS potential. Inside PySCF, in fact the density is constructed on a set of atom centered, radial grids. In the PySCF tutorial on DFT, you can see a snippet constructing the density on such a grid. You might also want to construct the density on a regularly spaced grid so that you can visualize it directly. This can be done using the cubegen module of PySCF.

Answer 2

Amending KobeGote's answer above, computing the density is in fact quite simple.

The molecular orbitals (MOs) are expanded in the basis set as, $$ \psi_i^\sigma({\bf r}) = \sum_\alpha C^\sigma_{\alpha i} \chi_\alpha({\bf r}), \tag{1}$$ where $\chi_i$ is the $i^{\textrm{th}}$ basis function and $\mathbf{C}^\sigma$ is the matrix of MO coefficients for spin $\sigma$. (Note that the MOs are stored by columns, and the number of MOs can be different from the number of basis functions: $\mathbf{C}^\sigma$ is typically not a square matrix.)

The spin-$\sigma$ electron density is, $$ n^\sigma({\bf r}) = \sum_{i \text{ occupied}} |\psi_i^\sigma({\bf r})|^2,\tag{2}$$ and the total electron density is simply, $$ n({\bf r}) = \sum_\sigma n^\sigma({\bf r}).\tag{3}$$

Evaluating the electron density can be done in two ways: either by expanding the square, $$ n^\sigma({\bf r}) = \sum_{i \alpha \beta} C^\sigma_{\alpha i} C^\sigma_{\beta i} \chi_\alpha({\bf r}) \chi_\beta({\bf r}) \equiv \sum_{\alpha \beta} P^\sigma_{\alpha \beta} \chi_\alpha({\bf r}) \chi_\beta({\bf r}) \tag{4}$$ where $$P^\sigma_{\alpha \beta} = \sum_i C^\sigma_{\alpha i} C^\sigma_{\beta i},\tag{5} $$ is the density matrix, or just by first computing the values of the occupied MOs and then squaring them and summing them together. If the basis set is large compared to the number of occupied orbitals (or one is interested in the density of only a single orbital), the latter approach is much more efficient.

Note that the density expansion is in the product basis $\chi_\alpha \chi_\beta$ which is:

• 1) different from the orbital basis $\chi_\alpha$ and,
• 2) pathologically linearly dependent. Due to the linear dependency, one can build a compact set of auxiliary basis functions $\chi_A$ in which the density can be expanded accurately as $n^\sigma({\bf r}) = \sum_A c^\sigma_A \chi_A({\bf r})$; this is the idea of the density-fitting (DF) approach which has been used to great success in quantum chemistry for several decades. I have recently proposed an approach for automatic formation of auxiliary basis sets for DF calculations in J. Chem. Theory Comput. 17, 6886 (2021).