Computer simulation has been increasingly present in the discovery and understanding of materials structure and properties. First-principles simulations, in particular, allow to meaningfully explore the material world without having to resource to experimental methods. Even though their extreme usefulness, first-principle methods tend to be very computer-intensive to the point of being impractical in many situations that are important for the understanding of materials’ behavior. As an alternative, computer simulations can be carried out using interatomic potentials (IPs), in which case an analytic mathematical model for the energy of the system is fitted to a set of crystal structure parameters and experimentally determined physical properties by minimizing a cost function. How close the properties calculated using an IP can get of those experimentally determined is limited both by its analytic form and the model’s parameters. For compounds with different atomic species, the multidimensional parameter space of the cost function most possibly exhibit a very complex landscape. In these cases, global minimization methods are often required. This work explores the use of a differential evolution (DE) algorithm for the parametrization of IPs for two compounds, namely berlinite (AlPO4) and zirconium tungstate (α-ZrW2O8). Several analytic interatomic potentials (including some previously proposed in the literature) were fitted to experimental lattice parameters, atomic positions, and elastic constants. Two-dimensional mappings of the interatomic potential parameters reveal the complex landscape of the cost function for berlinite and α-ZrW8O8, which exhibits several local minima, broad plateaus, discontinuities, and, in some cases, strong correlations between fitting parameters. The differential evolution algorithm was able to navigate the cost function landscape and, giving a sufficiently large population, was capable to found good candidates for the global minimum, despite all the mentioned difficulties. In all cases here explored, the interatomic potentials re-parametrized using DE yield calculated athermal properties in better agreement with the experimental observables as compared to previous results from the literature.