  Molybdenum clusters, characterised by their unique structure and intriguing
catalytic properties, have gained significant attention in recent years. In
several existing studies density functional theory (DFT) methods have been used
to find the lowest energy Mo clusters and explore their electronic and magnetic
structure. In all cases, with the exception of a single recent study, where a
genetic algorithm was employed, initial geometries of the clusters, prior to
geometry optimisation, were chosen using heuristic approaches based on symmetry
considerations and known structures. DFT calculations were performed using
different types of pseudopotentials, from hard to soft, and different types of
basis sets. However, no comprehensive study has yet been done in which a DFT
method with the best control on its precision would be complemented by a
reliable global minimum search method to find the lowest energy Mo clusters. In
this work, we employ a combination of a plane wave-based DFT method and
\emph{ab initio} random structure searching (AIRSS) technique to find the
lowest energy clusters of up to 10 Mo atoms. In each case, the search has been
performed for clusters with different spin multiplicities, which enabled us to
explore their magnetic structure. The results are compared for both hard and
soft pseudopotentials stressing the importance of treating more electrons
explicitly, in agreement with some of the previous studies. For most of the
low-energy magnetic structures found, we investigate the distribution of their
spin densities, and for all low energy clusters, we confirm their stability by
calculating their phonon structure. Finally, free energies of the Mo clusters,
within the quasi-harmonic approximation, are also calculated and discussed.
