  (Abridged) Understanding the details of the formation process of massive
(i.e. M<8-10M$_\odot$) stars is a long-standing problem in astrophysics. [...]
We present a method to derive accurate timescales of the different evolutionary
phases of the high-mass star formation process. We model a representative
number of massive clumps of the ATLASGAL-TOP100 sample which cover all the
evolutionary stages. The models describe an isothermal collapse and the
subsequent warm-up phase, for which we follow their chemical evolution. The
timescale of each phase is derived by comparing the results of the models with
the properties of the sources of the ATLASGAL-TOP100 sample, taking into
account the mass and luminosity of the clumps, and the column densities of
methyl acetylene (CH$_3$CCH), acetonitrile (CH$_3$CN), formaldehyde (H$_2$CO)
and methanol (CH$_3$OH). We find that the chosen molecular tracers are affected
by the thermal evolution of the clumps, showing steep ice evaporation gradients
from 10$^3$ to 10$^5$ AU during the warm-up phase. We succeed in reproducing
the observed column densities of CH$_3$CCH and CH$_3$CN, while H$_2$CO and
CH$_3$OH show a poorer agreement with the observed values. The total (massive)
star formation time is found to be $\sim5.2\times10^5$ yr, which is defined by
the timescales of the individual evolutionary phases of the ATLASGAL-TOP100
sample: $\sim5\times10^4$ yr for 70-$\mu$m weak, $\sim1.2\times10^5$ yr for
mid-IR weak, $\sim2.4\times10^5$ yr for mid-IR bright and $\sim1.1\times10^5$
yr for HII-regions phases. Our models, with an appropriate selection of
molecular tracers that can act as chemical clocks, allow to get robust
estimates of the duration of the individual phases of the high-mass star
formation process, with the advantage of being capable to include additional
tracers aimed at increasing the accuracy of the estimated timescales.
