  We present different methods to increase the performance of Hybrid Monte
Carlo simulations of the Hubbard model in two-dimensions. Our simulations
concentrate on a hexagonal lattice, though can be easily generalized to other
lattices. It is found that best results can be achieved using a flexible GMRES
solver for matrix inversions and the second order Omelyan integrator with
Hasenbusch acceleration on different time scales for molecular dynamics. We
demonstrate how an arbitrary number of Hasenbusch mass terms can be included
into this geometry and find that the optimal speed depends weakly on the choice
of the number of Hasenbusch masses and their values. As such, the tuning of
these masses is amenable to automization and we present an algorithm for this
tuning that is based on the knowledge of the dependence of solver time and
forces on the Hasenbusch masses. We benchmark our algorithms to systems where
direct numerical diagonalization is feasible and find excellent agreement. We
also simulate systems with hexagonal lattice dimensions up to $102\times 102$
and $N_t=64$. We find that the Hasenbusch algorithm leads to a speed up of more
than an order of magnitude.
