I have generated MOs of a molecular system and I need to recreate its electronic density via $\rho(r)=\sum_{A_i}|\psi_i(\mathbf{r})|^2$ where $A$ is a set of occupied MOs. I am trying to do this via extracting each MO individually in Gaussian16 Cubegen utility:
cubegen 16 MO=${i} input.fchk MO${i}.cube
I am looking for post-processing codes which can manipulate gaussian cube data, i.e., square and summation, preferably in a batch-processing fashion. For the codes I know, Gaussian's Cubman can add only two cubes files at a time; another code has issues with Gaussian's MO cube file format. I tried to understand PySCF implementation, but I am not able to fully grasp it. If someone has a simpler, standalone implementation for such a code, please let me know in the answers below.