With aoforce you get the harmonic frequencies. They are usually kind of ok for intramolecular stretch and bend modes that appear in the
infrared region.
In THz spectroscopy your are typically interested in intermolecular modes, typically large amplitude motions, e.g. in the solvent shells.
For these modes the harmonic approximation doesn't work well.
A more appropriate approach for this regime is to compute the spectrum with a MD appraoch via the dipole autocorrelation function.
