Here is a quick answer. Apparent retrograde motion of a planet as seen from Earth means that its ecliptic longitude is decreasing. (Thanks to @barrycarter for correcting my previous statement.) You have to calculate the apparent angular position of the planet from the Earth in ecliptic coordinates (as explained here and here) and find places where the ecliptic longitude is decreasing.
You can subtract the positions of Earth from the position of the planet in either coordinates. Then you express that vector in spherical coordinates. It doesn't matter if it's heliocentric or geocentric. What matters is the direction that the axes of the coordinate system is pointing, not where the center is.
If you want RA and dec, use ICRF coordinates, and these can be transformed to ecliptic coordinates as well. An easy way to do this is to use the python package Skyfield.
Here's a sample script and a plot. This is Mars from 2010.0 to 2020.0 and the thick red line represents retrograde motion.
import numpy as np
import matplotlib.pyplot as plt
from skyfield.api import load
planets = load('de421.bsp')
earth = planets['earth']
mars = planets['mars']
year_zero = 2010
days = np.linspace(1, 3650, 10000)
years = year_zero + days / 365.2564
ts = load.timescale()
t = ts.utc(year_zero, 1, days)
# thanks to @barrycarter's comment, do it the right way!
eclat, eclon, ecd = earth.at(t).observe(mars).ecliptic_latlon()
eclondgs = (180./np.pi) * eclon.radians
eclondel = eclondgs[1:] - eclondgs[:-1]
eclondel[eclondel < -300] += 360. # this is a fudge for now
eclondel[eclondel > +300] -= 360. # this is a fudge for now
prograde = eclondel > 0.
eclon_prograde = eclondgs.copy()[:-1]
eclon_retrograde = eclondgs.copy()[:-1]
eclon_prograde[-prograde] = np.nan
eclon_retrograde[prograde] = np.nan
plt.plot(years[:-1], eclon_prograde, '-g', linewidth=1)
plt.plot(years[:-1], eclon_retrograde, '-r', linewidth=3)
plt.title('Mars Ecliptic longitude (degrees) AD 2010.0 to 2020.0',