Calculating the density matrix from a WFN file
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.