PSI4 is making use of point group symmetry of water, i.e C2V which has 4 irreps, A1, A2, B1, B2. The following is the occupancy of different irreps: A1 - 4, A2 - 0, B1 - 1, B2 - 2
Since the wavefunction transforms as fully symmetric irrep, i.e A1, you will get a 4x4, 0x0, 1x1 and 2x2 matrix, corresponding to A1 * A1, A2 * A2, B1 * B1, B2 * B2 combinations.
You can get the 5x7 matrix if you set the symmetry as c1 in your input file:
h2o = psi4.geometry(""“
H 1 0.96
H 1 0.96 2 104.5
and then use
Cocc = wfn.Ca_subset("AO", "OCC")