Calculating energy function with 1&2-RDM (reduced density matrix) in Qiskit

I am now working on designing quantum error mitigation with variational quantum eigensolvers (VQE). I need to get the 1&2-RDM from Qiskit, which has already been realized by Openchemist pakages for Qiskit to get RMDS. The question is how to get the energy function in the following formula (A14) from appendix A of The Google Hartree Fock paper. I need (A14) to be the objective function for gradient descent, then the VQE algorithm can be implemented.

Because the problem here is to calculate the ground state energy for a quantum chemistry problem, the first step is to extract the molecular information, which is contained in the tensors, h_ij and V_ijkl in (A14). Then, the objective function can be defined for gradient descent. Is there a way to realize the two step or there are other ways to solve this in Qiskit?