# How to get atomic orbital coefficients in PySCF

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).

• +1 Welcome to our new community and we hope to see much more of you in the future!!! Thank you for contributing your question here. How did you find us? I'm sure "matter modeling stack exchange" isn't the first thing that comes to mind for someone that wants to ask about PySCF! Jan 20 at 7:01
• @NikeDattani could just be for me, but the pyscf tag here is almost on the first page of Google search results for PySCF.
– Tyberius
Jan 20 at 13:40
• I googled "PySCF stack exchange" and this forum came up first :) Jan 20 at 16:24

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.

• It was getting a bit long, so I moved the conversation to chat. KobeGote, you will need to ping @SusiLehtola if you want them to see your comment.
– Tyberius
Jan 25 at 2:14

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).
• +1, just added the equation labels in case if someone else wants to refer to "Eq. 4 of this post" somewhere else. Is "i:th" a proper way to say $i^\textrm{th}$? You can revert that part back if you want! I didn't see it in this, but didn't look too far. As for the bullet points, that makes it look clearer to me, but feel free to revert if you want! Feb 2 at 0:14