# Calculating the density matrix from a WFN file

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?

• What do you mean by "the first part which has the minimum energy"? And welcome on Chem.SE! Apr 2 '18 at 11:31
• Thanks. My mean was in my original .wfn file that the energy of each part has been written on the first line of it. The most negative energy belongs to first MO and it increases as the MO's goes on. Apr 2 '18 at 15:39
• Calculating a density matrix by hand from a wfn/wfx file is possible, but even for small molecules a time consuming work. The wfn-file (AIM traditional wavefunction data) is better read in by a program of choice (e.g. multiwfn, see sobereva.com/multiwfn) which then can calculate the desired properties. Apr 3 '18 at 20:11
• Do you just want to have a density matrix of your wavefunction or do you want to have a real electron density? Could you clarify this in your question, please. Apr 3 '18 at 20:22
• @awvwgk Thanks for your suggestion but notice I don't want to calculate electron density value for a specific coordinate rather I wanna to obtain the electron density as a function of distance (R). Otherwise I could simply obtain electron density using Softwares such as AIMALL or AIM2000. Apr 4 '18 at 3:41

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.