I am trying to calculate the above mentioned RDMs from the MP2 calculation. For now, using psi4.energy I can just calculate the MP2 corrected energy and obtain the wavefunction object. But is it possible to obtain these RDMs as well after the MP2 calculation? Or a way to construct them from the MP2 wavefunction object? Any guidance would be great help. Thank you!
I need to ask a few technical questions before I know how to answer:
 What approximations do you want to use on top of MP2? Frozen core? Frozen virtual? Cholesky? Density fitting? DLPNO?
 What do you need the 3 and 4 RDM for? Psi doesn’t store these, as we have no need for them, and I’m not sure there’s even a standard way to define these. (The only thing I can think of is using MP2 as a source of RDMs for a multireference perturbation theory.)
 Do you want RDMs that include orbital relaxation terms or not?
 How big of a system do you want this for? You may want to create your own code for the 3RDM and 4RDM using Psi4Numpy, but I’d need to think on that more before recommending it conclusively.
Thank you for the quick reply. Actually I am not from chemistry background and hence my knowledge is not too deep but I will explain what exactly I am trying to calculate and I hope it be of help.
I am trying to calculate what are known as the one and two “orbital” reduced density matrices. I found out they can be calculated using 1,2,3 and 4 “electron” RDMs (Refer Orbital Entanglement in Quantum Chemistry). I want to make use of a basic post HartreeFock method as I will be using quantum computing later. Hence, I came to MP2. Below I have given the answers for your questions with my basic understanding of the subject.

I believe I am not looking for any particular approximation. If I can do just the MP2 calculation without any approximation, that would also work for me. However, if an approximation gives a better nbody wavefunction then I would prefer that. If it helps, the wavefunction is more important to me than the energy correction offered by these approximation methods.

Partly I have given the answer for this question in the above paragraph. Now for the definition of 3 and 4 RDMs, what I have in mind are the terms like <a_a^\dagger a_b^\dagger a_c^\dagger a_i a_j a_k > and <a_a^\dagger a_b^\dagger a_c^\dagger a_d^\dagger a_i a_j a_k a_l > respectively, where the expectation is wrt the nparticle wavefunction obtained by the postHF method and i,j,k,l and a,b,c,d are the occupied and virtual molecular orbitals respectively.

I am not exactly sure but if orbital relaxation gives a more accurate nbody ground state wavefunction then that would be preferable.

I am mostly working with small molecules like H2, LiH, H20 etc which can be solved on a quantum computer.
It would also work if I could just obtain the nbody wavefunction from MP2 calculation in the form of linear combination of slater determinants of HF state and its excitations (like in configuration interaction.) In this way I would be able to directly calculate the “orbital” RDMs without the requirement of 1,2,3 and 4 electron RDMs.
Hope this helps. If you want me to clarify anything further please let me know. Thank you.
You need to understand a bit more about electronic structure methods and reduced density matrices before you can proceed.
Many electronic structure methods (including MP2) do not determine the energy as the expectation value of some particular wavefunction. This can lead to them taking energies below FCI. The reduced density matrices of such methods are not defined as in the formula you gave. It sounds like the closest match to what you’re looking for would be to get the first order wavefunction correction computing during an MP2 computation. The energy of that wavefunction will be different from the MP2 energy, but that’s at least the simplest approximation to a wavefunction. Then you can take the RDMs of that.
Is that what you’re looking for?
I believe yes. If that is the closest I can get to the nbody wavefunction (expressed in a form similar to CI) from an MP2 calculation then I suppose it will work. And yes I will be able to calculate orbital RDMs directly from that.
Based on everything you’ve said, I’m going to recommend that you compute these via Psi4Numpy, using Psi4 for integrals but computing everything else with NumPy. This way, you understand what you’re doing and don’t need to go through the steep learning curve of working with Psi’s internal tensor library. The computational cost should be negligible for systems like you’re describing. See the equation between (1) and (2) here for the final energy expression. For the wavefunction expression, the reference determinant has coefficient 1. All other doubly excited determinants have weight given by removing one of the two twoelectronintegrals, i.e., the weight of the excitation from (i, j) to (a, b) is [iajb] / (e_i + e_j  e_a  e_b). That formula does not count (i, j) and (j, i) separately. You can count them separately if you wish if you include a factor of 1/2. The same applies for (a, b) and (b, a), which have their own factor of 1/2 if you want to count both.
Okay. I think I will be able to calculate what I want from the explanation you have given. Thank you.