When doing a frequency calculation (psi4.frequency(return_wfn="on")), the force constants are printed on the output. Are they stored anywhere on the wfn object? I can extract the frequencies and the hessian matrices from the Wavefunction object, but I don’t see a method for the force constants. Is it possible to extract them into an numpy array?
What about the normal modes?
Are you’re referring to the “force constants” given for each vibrational frequency in the output? To my knowledge, those values are not provided in the wfn object, but they can be calculated from the Hessian and geometry. The process of that calculation is long, depending on your background and interest level, but the details are provided here for the Gaussian code, though the approach is basically the same for all quantum chemistry programs.
It might help to simplify things if you could explain what you want to do with the force constants once you have them. No guarantees, though! 
Yes, for example the force constants printed in the following output when you do a frequency calculation:
Vibration 7 8 9
Freq [cm^-1] 1842.5361 3876.1906 3950.9309
Irrep A1 A1 B2
Reduced mass [u] 1.0804 1.0473 1.0798
Force const [mDyne/A] 2.1612 9.2707 9.9312
Turning point v=0 [a0] 0.2459 0.1722 0.1680
RMS dev v=0 [a0 u^1/2] 0.1808 0.1246 0.1234
IR activ [km/mol] 75.1412 19.5437 46.6616
Char temp [K] 2650.9992 5576.9753 5684.5099
So that means the force constant are calculated internally, but I guess they’re not stored anywhere that can be conveniently accessed through the python API.
I did find a function that sort of does this: psi4.driver.qcdb.vib.harmonic_analysis. It’s not as convenient to use, but one of the dictionary entries it outputs are the force constants.
I am using QC-derived force constants for force-field development.