Model-based molecular phylogenetic analysis plays a key role in most comparisons of genomic data, and model selection is a critical step in all such analyses. We introduce ModelFinder, a novel model-selection method that: (a) incorporates a model of rate-heterogeneity across sites (RHAS) not previously considered in this context, and (b) allows the topology of the tree to be a variable during the search for an optimal model of evolution. Using alignments of amino acids obtained by simulation on a 100-tipped tree, ModelFinder is found to produce accurate estimates of RHAS, where other model-selection methods fail (e.g., when the model of evolution incorporates a bimodal distribution of RHAS). The benefits of ModelFinder are demonstrated using genomic data from a Tree-of-Life study involving Bacteria and Archaea. The results from our analysis of these data reveal that: (a) the optimal model of RHAS is trimodal, (b) the popular Γ-distribution-based model of RHAS provides a poor fit of the data, and (c) the optimal trees inferred under these two models of RHAS are very different. ModelFinder is implemented in the phylogenetic program IQ-TREE.