Input should be the AO density matrix and the AOs on the grid. The output is a rank 3 tensor that is "number of grid points" by "number of AOs" by "number of AOs".
I would just code up "X alpha" or the Slater functional (see if either is in GauXC so you can compute the correct answer).