We model neutron star cooling with several microscopic nuclear equations of state based on different nucleon-nucleon interactions and three-body forces, and compatible with the recent GW170817 neutron star merger event. They all feature strong direct Urca processes. We find that all models are able to describe well the current set of cooling data for isolated neutron stars, provided that large and extended proton 1S0 gaps and no neutron 3PF2 gaps are active in the stellar matter. We then analyze the neutron star mass distributions predicted by the different models and single out the preferred ones.