Here are couple of references that describe professional uses of a post-Newtonian formalism to model the planets and the Earth's Moon:
Standish, et al. "Orbital ephemerides of the Sun, Moon, and planets," Explanatory Supplement to the Astronomical Almanac (1992): 279-323.
The relevant equation is 8-1 on page 3.
Petit and Luzum (eds.), "IERS Technical Note No. 36, IERS Conventions (2010)," International Earth Rotation and Reference Systems Service, Frankfurt, Germany (2010).
You'll want to look at chapter 10. This describes the Moon, but not the planets. The IERS's focus is the Earth.
Note that the IERS document does describes relativistic time scales. Mixing and matching the non-relativistic time scales we use to keep time on the surface of the Earth with a post-Newtonian formalism doesn't buy much. JPL (the first reference) uses its own relativistic time scale, Teph. This is now very close to one of the officially defined time scales, Barycentric Dynamical Time (TDB). The IERS document uses a different relativistic time scale that doesn't tick with Earth-based clocks.
(Note: The acronym doesn't match the name. That's intentional. The official acronym is short for the French name "Temps Dynamique Barycentrique", regardless of what language you use for the full name.)
Edit: Regarding the simplified version specified in the question
My old, old copy of Marion, Classical Dynamics has something very close to the formula in the question. It extends the Newtonian relation
$$ \frac{d^2}{d\theta^2}\left(\frac 1 r\right) + \frac 1 r
= -\frac m {l^2}r^2 F(r)
= GM \frac {m^2}{l^2}$$
to
$$ \frac{d^2}{d\theta^2}\left(\frac 1 r\right) + \frac 1 r
= GM \frac {m^2}{l^2} + \frac {3GM} {r^2 c^2} $$
from which one the force modified $F(r)$ can be expressed as
$$F(r) = -\frac {GMm}{r^2}\left( 1 + \frac{3 \, l^2}{m^2c^2 r^2} \right)$$
Using $\vec l \equiv \vec r \times (m \vec v)$, this can be rewritten as
$$F(r) = -\frac {GMm}{r^2}
\left( 1 + \frac{3 \, ||\vec r \times \vec v||^2}{c^2 r^2} \right) $$