As the earliest relics of star formation episodes in the history of the Universe, the most massive galaxies are the key to our understanding of the evolution of the stellar population, cosmic structure, and supermassive black holes (SMBHs). However, the details of their formation histories remain uncertain. We address these problems by planning a large survey sample of 101 ultramassive galaxies (𝑧 ≤ 0.3, |𝛿 + 24 • | < 45 • , |𝑏| > 8 • ), including 76% ellipticals (E), 17% lenticulars (S0), and 7% spirals (S) brighter than 𝑀 𝐾 ≤ −27 mag (stellar mass 2 × 10 12 < 𝑀 ★ < 5 × 10 12 M ) with the High Angular Resolution Monolithic Optical and Near-infrared Integral field spectrograph (HARMONI) instrument on the Extremely Large Telescope (ELT). Our sample comprises diverse galaxy environments ranging from isolated to dense-clusters galaxies that have not been studied systematically to date. The primary goals of the project are to (1) explore the stellar dynamics deep inside galaxy nuclei and weigh the central supermassive black holes (SMBHs), (2) constrain the black hole scaling relations at the highest-galactic mass ladder, and (3) probe the late-time assembly of these most massive galaxies through the stellar population and kinematical gradients. In this paper, we describe our most massive galaxies survey, discuss the distinct demographics and environmental properties of the selected galaxies, and simulate their HARMONI 𝐼 𝑧 , 𝐼 𝑧 + 𝐽, and 𝐻 + 𝐾-band Integral Field Spectroscopy (IFS) observations via combining the inferred stellar-mass models from the images of Pan-STARRS observations, a synthetic spectrum of the assumed stellar population, and an SMBH with a specific mass estimated based on different black hole scaling relations. We find HARMONI can produce excellent state-of-art integral field unit (IFU) datacubes even for the targets beyond the limits of the current suite of 8m class ground-based telescopes with adaptive optic assisted like Gemini and Very Large Telescope in a relatively short exposure time. We extract the stellar kinematics from the simulated IFU data using the stellar-absorption CaT and some stellar features in (𝐼 𝑧 , 𝐼 𝑧 + 𝐽) and 𝐻 + 𝐾 band, respectively; then used them to reconstruct the SMBH mass and its error using the Markov Chain Monte Carlo (MCMC) simulation. The reconstructed stellar kinematics were achieved with a median random noise of Δ𝑉 rms 1.5% for all nine simulated galaxies, which is representative for our MMBH survey, across the HARMONI field-of-view. Thus, our simulations and modelings can be used as benchmarks to evaluate the instrument models and pipelines dedicated to HARMONI, appearing to be promising tools to exploit the unprecedented capabilities of ELT in future works.