We present a specially designed evolutionary algorithm for the prediction of surface reconstructions. This technique allows one to automatically explore stable and low-energy metastable configurations with variable surface atoms and variable surface unit cells through the whole chemical potential range. The power of evolutionary search is demonstrated by the efficient identification of diamond 2 × 1 (100) and 2 × 1 (111)surface reconstructions with a fixed number of surface atoms and a fixed cell size. With further variation of surface unit cells, we study the reconstructions of the polar surface MgO (111). Experiment has detected an oxygen trimer (ozone) motif [Plass et al., Phys. Rev. Lett. 81, 4891 (1998)]. We predict another version of this motif which can be thermodynamically stable in extreme oxygen-rich conditions. Finally, we perform a variable stoichiometry search for a complex ternary system: semipolar GaN (1011) with and without adsorbed oxygen. The search yields a counterintuitive reconstruction based on N3 trimers. These examples demonstrate that an automated scheme to explore the energy landscape of surfaces will improve our understanding of surface reconstructions. The method presented in this paper can be generally applied to binary and multicomponent systems.