The classification of mineralized areas into different groups based on mineral grade and prospectivity is a practical problem in the area of optimal risk, time, and cost management of exploration projects. The purpose of this paper was to present a new approach for optimizing the grade classification model of an orebody. That is to say, through hybridizing machine learning with a metaheuristic algorithm called Harmony Search (HS), a proper model for the spatial distribution of the grade classes was obtained, while improving the computational cost of the traditional classification methods. The HS is an algorithm inspired by the simulation of the process where a composer tries to harmonize a piece of music. By interpolating the dataset of Cu and Mo concentrations in the surface rock samples taken from Kooh-Panj deposit district, the grid data of the two elements were extracted. To estimate the true number of groups in the dataset, five popular indices were used in this regard, which in turn determined two classes as the optimal number of groups. Harmony Search Learning (HSL) was used to classify the grid dataset of Cu and Mo. The comparison of the results of the proposed approach with the conventional k-means clustering suggested that the use of HSL method significantly reduced the cost function of the problem (up to 13%). The adaptation of the mineralization class derived from HSL and k-means clustering to borehole locations proved that the results of the HSL were more successful in the accurate estimation of the economic mineralization class identified by the exploratory excavations. In this respect, the HSL technique could significantly improve the k-means performance by 25%. Furthermore, the results of the HSL were more consistent with the lithological units and alteration zones involved in the ore-forming processes. The use of the HS-based learning rectified the disadvantages arising from the typical clustering methods regarding the entrapment in local optimums. It also led to the extraction of weak mineralization signals, numerically laid in boundary conditions. This approach can be extended to more than two geochemical variables and can be a valuable tool for the classification of the mineralized areas to design and optimally manage the mineral exploration projects.