Higher-order derivatives not working -- Libint2 limitation?

[chertianser]

I'm attempting to reproduce the example after following the installation instructions provided in the README. Quax runs fine for the first and second order derivative calls, but orders above two will fail:

import quax 
import psi4

molecule = psi4.geometry('''
                         0 1
                         H 0.0 0.0 -0.80000000000
                         H 0.0 0.0  0.80000000000
                         units bohr
                         ''')

dz1_dz1_dz2_dz2 = quax.core.geom_deriv(molecule, '6-31g', 'ccsd(t)', deriv_order=4, partial=(2,2,5,5))

python3: /h/292/ctser/functor/libint/BUILD/libint-2.9.0/PREFIX/include/libint2/./engine.impl.h:483: const target_ptr_vec& libint2::Engine::compute1(const libint2::Shell&, const libint2::Shell&): Assertion `deriv_order_ <= 2 && "feature not implemented"' failed.
/var/spool/slurmd/job216925/slurm_script: line 15: 1964948 Aborted                 (core dumped) python3 run_quax.py

The check for derivative order stated at the end of issue #1 have been since removed from quax/methods/ints.py. Could this be reinstated somehow? Thanks!