Magnetic fields are expected to play a key role in the dynamics and the ejection mechanisms that accompany the merger of two neutron stars. General relativistic magnetohydrodynamic (MHD) simulations offer a unique opportunity to unravel the details of the ongoing physical processes. Nevertheless, current numerical studies are severely limited by the fact that any affordable resolution remains insufficient to fully capture the small-scale dynamo, initially triggered by the Kelvin-Helmholtz instability, and later sourced by several MHD processes involving differential rotation. Here, we alleviate this limitation by using explicit large-eddy simulations, a technique where the unresolved dynamics occurring at the sub-grid scales (SGS) is modeled by extra terms, which are functions of the resolved fields and their derivatives. The combination of high-order numerical schemes, high resolutions, and the gradient SGS model allow us to capture the small-scale dynamos produced during the binary neutron star mergers. Here we follow the first 50 milliseconds after the merger and, for the first time, we find numerical convergence on the magnetic field amplification, in terms of integrated energy and spectral distribution over spatial scales. We also find that the average intensity of the magnetic field in the remnant saturates at \(\sim 10^{16}\)~G around \(5\)~ms after the merger. After \(20-30\)~ms, both toroidal and poloidal magnetic field components grow continuously, fed by the winding mechanism that provides a slow inverse cascade. We find no clear hints for magneto-rotational instabilities, and no significant impact of the magnetic field on the redistribution of angular momentum in the remnant in our simulations, probably due to the very turbulent and dynamical topology of the magnetic field at all stages, with small-scale components largely dominating over the large-scale ones.