I’ll start off by saying I don’t know how to get the MO array directly from Psi4. Hopefully this is easy to do with Psi4, but I haven’t seen it in the documentation.

However, I help develop a python module called AaronTools which can read FCHK files and give you MO values. I see no benefit to doing this with AaronTools if Psi4 can just give you the array. Here is a page describing how to get an FCHK file from Psi4. Here is the installation guide for AaronTools.

The version of AaronTools on PyPI is missing the convenience functions for constructing a reasonable grid of points, but these functions are present on the version currently on GitHub. The version on PyPI is still capable of giving you the MO values for an array of points, you will just have to construct the grid yourself.

The code with the GitHub version would look something like:

```
from AaronTools.fileIO import FileReader
from AaronTools.geometry import Geometry
# read the FCHK file and grab orbital stuff (not just the structure)
fr = FileReader("test.fchk", just_geom=False)
# Geometry is basically a molecule object
geom = Geometry(fr)
# get the Orbitals object from the parsed data
orbitals = fr.other["orbitals"]
# number of threads that will be used when evaluating an MO
# more threads = faster
# probably best to not exceed the number of cores on your CPU
n_threads = 8
# get grid primitives with the specified resolution
# the grid is centered on the molecule's centroid
# and oriented to minimize the volume
# n_ptsN - number of points along each axis
# vN - vector along each axis
# com - origin of the cube
# u - rotation matrix for vN relative to x, y, and z
# this is not on PyPI right now
n_pts1, n_pts2, n_pts3, v1, v2, v3, com, u = orbitals.get_cube_array(
geom,
padding=4, # extra space in angstroms around atoms
spacing=0.1, # target resolution in angstroms
standard_axes=False, # use True to get vN to be x, y, and z
)
# optional - check the max. memory that will be used when
# evaluating the MO if the result higher than your computer's
# free memory, maybe don't continue
# this is not on PyPI right now
n_val = n_pts1 * n_pts2 * n_pts3
mem_used = orbitals.memory_estimate(
"mo_value",
n_points=n_val,
n_atoms=len(fr.atoms),
n_jobs=n_threads,
)
# using 8 GB as an example
if mem_used > 8:
raise RuntimeError("insufficient memory: %.1fGB requested" % mem_used)
# get the actual points in the box
# points - array of coordinates (N_points x 3)
# n_list - number of points along each axis (some permutation
# of [n_pts1, n_pts2, n_pts3])
# this is mainly an accounting thing for writing .cube files
# sort=True is an accounting thing for writing cube file
# to get the points to increment
# along a particular axis in a predictable order
# sort=False means n_list will be [n_pts1, n_pts2, n_pts3]
# this is not on PyPI
points, n_list = orbitals.get_cube_points(
n_pts1, n_pts2, n_pts3, v1, v2, v3, com, sort=False
)
# get orbital values
# this is on PyPI
# first argument is index of MO (0-indexed) or an array
# of AO coefficients
# points is the (N x 3) array
# assuming closed shell...
homo_values = orbitals.mo_value(
orbitals.n_alpha - 1,
points,
n_jobs=n_threads
)
lumo_values = orbitals.mo_value(
orbitals.n_alpha,
points,
n_jobs=n_threads
)
```

AaronTools will print an error message if it encounters a Gaussian shell it can’t handle. It should be able to work with up to i functions for pure basis sets, and I think only up to f or g shells for Cartesian functions.