How to get density matrix derivatives over nuclear positions in SCF calculations?

I need density matrix first derivatives with respect to nuclear positions in SCF (HF or DFT) calculations. I need that for developing of some constraint SCF related approach.

What is the easiest way to get it?
I realize that wfn.cphf_solve can be used to get density response if the hamiltonian perturbation is known. I figured out what wfn.cphf_solve has on input and output. It is handy and not difficult to use.

I am also aware that MintsHelper.ao_oei_deriv1(TYPE, atom) with TYPE = ‘KINETIC’, ‘POTENTIAL’ and ‘OVERLAP’ can give me the derivatives of one-electron hamiltonian and of overlap matrix.

But there is also derivative of two-electron J - K hamiltonian term. And it does not seem like it is easy to obtain. Could someone help with that? (or may be there is another route to density matrix derivatives)?

Thanks a lot.