  Large-scale applications of energy density functional (EDF) methods depend on
fast and reliable algorithms to solve the associated non-linear
self-consistency problem. When dealing with large single-particle variational
spaces, existing solvers can become very slow, and their performance dependent
on manual fine-tuning of numerical parameters. In addition, convergence can
sensitively depend on particularities of the EDF's parametrisation under
consideration. Using the widely-used Skyrme EDF as an example, we investigate
the impact of the parametrisation of the EDF, both in terms of the operator
structures present and the size of coupling constants, on the convergence of
numerical solvers. We focus on two aspects of the self-consistency cycle, which
are the diagonalisation of a fixed single-particle Hamiltonian on one hand and
the evolution of the mean-field densities and potentials on the other.
Throughout the article we use a coordinate-space representation, for which the
behaviour of algorithms can be straightforwardly analysed. We propose two
algorithmic improvements that are easily implementable in existing solvers,
heavy-ball dynamics and potential preconditioning. We demonstrate that these
methods can be made virtually parameter-free, requiring no manual fine-tuning
to achieve near-optimal performance except for isolated cases. The combination
of both methods decreases substantially the CPU time required to obtain
converged results. The improvements are illustrated for the MOCCa code that
solves the self-consistent HFB problem in a 3d coordinate space representation
for parametrisations of the standard Skyrme EDF at next-to-leading order in
gradients and its extension to next-to-next-to-leading order.
