Calculating the density matrix from a WFN file

Question

I want to calculate the density for the $\ce{HCN}$ molecule. This molecule has 14 electrons and thus 7 occupied molecular orbitals. I have a .wfn file which includes 8 parts, where each of them has 14 rows and 5 columns (due to 5 basis functions). The first part belongs to the basis functions (exponents) and the 7 other parts contain coefficients for each of seven MO's. How can I calculate the electron density using this file? Generally I have used the electron density formula based on the Roothaan expansion introduced in Szabo's textbook, namely:

$$ P_{\mu\nu} =2\sum_{a=1}^{N/2} C_{\mu a}C_{\nu a} $$

where $a$ runs over the number of electrons, though I have seen in several codes that $a$ is considered as the number of occupied MOs. $\mu$ and $\nu$ stand for the number of basis functions (5 here).

Should I use all 7 parts of the coefficients to obtain the electron density of the $\ce{HCN}$ ground state, or is the first part which has the minimum energy sufficient?

Answer 1

where a runs over the number of electrons, though I have seen in several codes that a is considered as the number of occupied MOs.

If it is closed shell (all spatial orbitals are either doubly occupied or empty), then running up to $N/2$ is the same as running over all occupied orbitals. Note that your formula is for closed shell systems described by a single determinant.

There are two representations of the first order reduced density matrix, or simply density matrix: the discrete one and the continuous one. What you have in the formula is the discrete one, represented in the space of primitive functions. The only reason why it is discrete is because the original Schrödinger equation has been projected on a finite dimensional space. You need to locate the Molecular Orbital Primitive Coefficients ($C_{\mu a}$) in the file. You could as well calculate the overlap matrix ($S$) of the primitives and then get the density matrix in the basis of molecular orbitals $$ D = C^t S P S C.$$ You will see that it is a diagonal matrix with 2's.

You could also represent in real space the general density matrix, which is continuous and spans $\mathbb{R}^3$. You could discretize the space and compute the values there. Then, the general formula $$ \rho(\pmb{r},\pmb{r}') = \sum_{a,b} D_{ab}\phi_a(\pmb r) \phi_b^*(\pmb r'),$$ where $\{\phi_a\}_{a=1}^M$ represent the spatial canonical molecular orbitals, and becomes $$ \rho(\pmb{r},\pmb{r}') = 2\sum_{a}^{N/2} \phi_a(\pmb r) \phi_a^*(\pmb r')$$ for a closed-shell single determinant wavefunction. Setting $\pmb r' = \pmb r$ you get the density.

If you plan to do it you will have to deal with some nuances, like the normalization used by the program that generated the .wfn file.