We simulate 4dSU(N) pure-gauge theories at large N using a parallel tempering scheme that combines simulations with open and periodic boundary conditions, implementing the algorithm originally proposed by Martin Hasenbusch for 2dCPN−1 models. That allows to dramatically suppress the topological freezing suffered from standard local algorithms, reducing the autocorrelation time of Q2 up to two orders of magnitude. Using this algorithm in combination with simulations at non-zero imaginary θ we are able to refine state-of-the-art results for the large-N behavior of the quartic coefficient of the θ-dependence of the vacuum energy b2, reaching an accuracy comparable with that of the large-N limit of the topological susceptibility.