Joint lensing and dynamical mass profile determinations of galaxy clusters are an excellent tool to constrain modification of gravity at cosmological scales. However, search for tiny departures from General Relativity calls for an accurate control of the systematics affecting the method. In this analysis we concentrate on the systematics in the reconstruction of mass profiles from the dynamics of cluster member galaxies, while assuming that lensing provides unbiased mass profile reconstructions. In particular, in the case study of linear $f(R)$ gravity, we aim at veryfying whether in realistic simulations of cluster formation a spurious detection of departure from GR can be detected due to violation of the main assumptions (e.g. dynamical equilibrium and spherical symmetry) on which the method is based. We aim at identifying and calibrating the impact of those systematics by analyzing a set of Dark Matter halos taken from $Lambda$CDM N-body cosmological simulations performed with the GADGET-3 code. [...] If no selection criteria are applied, $sim 60%$ of clusters in a $Lambda$CDM Universe (where GR is assumed) produce a spurious detection of modified gravity. We find that the probability of finding cluster in agreement with GR predictions $P_{GR}$ mainly depends on the properties of the halos projected phase-space and on shape orientation of the cluster along the line-of-sight projection. We define two observational criteria which correlate with the probability to find clusters in agreement with GR predictions and which can be used to select [...] those objects that are more suitable for the application of the proposed method. In particular, we find that according to these criteria the percentage of spurious detection can be lowered down to $sim 20%$ in the best case. Our results are relevant in view of data that will be available with the next generation surveys.