You might want to look at Folkner, et al theory paper,
"The Planetary and Lunar Ephemerides DE430 and DE431"
And also the JPL horizons ephermeides referenced by this paper, which is already implemented, does what you describe, and has a public interface so I believe you could personally access it. But learning to operate it appears to be quite technical - I know of the existence of the program but not how to operate it. Documentation is available on the JPL website, but appears quite technical.
Here is a list of what is considered taken from the above paper.
Thus we have are relativistic equations for point masses (based on linearized approximation of General relativity), a Newtonian analysis of what happens because of the "figure" (i.e. the non-spherical shape) of the various solar system bodies, and some tidal effects (the figure of the Earth's oceans doesn't stay constant, it changes due to the tides).
Goldstein's book, "Classical mechanics", has some discussion of the effect of figure on the Earth-moon system if you want more details of how the whole idea of "figure" is handled. But note it's a graduate level text (though it's introductory at the graduate level). Basically one expands the (Newtonian) gravitational potential in terms of series expansion via spherical harmonics.
Folkner points out that the moon needs a particularly complex model, due to lunar mass concentrations. So one needs a lot of the spherical and zonal harmonics, more so than for other bodies.
It's possible to do the analysis of the effects of "figure" in the context of full GR, but it's simpler to analyze it with Newtonian methods.
There are effects not covered in the JPL program that don't matter for the solar system models, but could matter for more exotic problems, effects such as the emission of gravitational waves due to the orbital motion of the bodies. Simulations that take into account more than linearized gravity are needed if one wants to include these sort of effects. Such simulations have been done for black hole mergers (binary inspirals), but I'm afraid I don't know the details, other than getting them working was a project for the entire scientific community.
JPL's program is rather old, so the way it is set up may not follow the most modern conventions. Probably with the right options, it can be configured to output the numbers according to modern conventions, but that's really just a guess.
An overview of the modern conventions would probably be the wiki article on the ICRS, the Inernational Celestial Reference System, see for instance
<<wiki link>>.
The technology behind the ICRS is very long baseline interferometry, which gives extremely precise measurements of the position of extragalactic radio sources, typically quasars.
A general observation - getting into the full level of detail of how we handle the solar system simulations is interesting, but to fully appreciate it one probably needs to understand General relativity first. It's not necessarily a good project to take on to learn about how GR works. The usual approach would be to learn the basics of GR, learn about approximations to GR such as PPN, then learn about the refinements to PPN needed to allow conversion between geocentric (earth-centered) and barycentric (solar system centered) coordinates.