We simulate Nf = 2 + 1 QCD at the physical point combining open and periodic boundary conditions in a parallel tempering framework, following the original proposal by M. Hasenbusch for 2d CPN−1 models, which has been recently implemented and widely employed in 4d SU(N) pure Yang-Mills theories too. We show that using this algorithm it is possible to achieve a sizable reduction of the auto-correlation time of the topological charge in dynamical fermions simulations both at zero and finite temperature, allowing to avoid topology freezing down to lattice spacings as fine as a ∼ 0.02 fm. Therefore, this implementation of the Parallel Tempering on Boundary Conditions algorithm has the potential to substantially push forward the investigation of the QCD vacuum properties by means of lattice simulations.