We present a novel fitting formula for the halo concentration enhancement in chameleon $f(R)$ gravity relative to General Relativity (GR). The formula is derived by employing a large set of $N$-body simulations of the Hu-Sawicki $f(R)$ model which cover a wide range of model and cosmological parameters, resolutions and simulation box sizes. The complicated dependence of the concentration on halo mass $M$, redshift $z$, and the $f(R)$ and cosmological parameters can be combined into a simpler form that depends only on a rescaled mass $M/10^{p_2}$, with $p_2equiv1.5log_{10}left[|{bar{f}_R(z)}|/(1+z)right]+21.64$ and $bar{f}_R(z)$ the background scalar field at $z$, irrespective of the $f(R)$ model parameter. Our fitting formula can describe the concentration enhancement well for redshifts $zleq3$, nearly 7 orders of magnitude in $M/10^{p_2}$ and five decades in halo mass. This is part of a series of works which aims to provide a general framework for self-consistent and unbiased tests of gravity using present and upcoming galaxy cluster surveys. The fitting formula, which is the first quantitative model for the concentration enhancement due to chameleon type modified gravity, is an important part in this framework and will allow continuous exploration of the parameter space. It can also be used to model other statistics such as the matter power spectrum.