With the aid of ab initio calculations, a realistic interatomic potential was constructed for the Mg-Cu-Y ternary system under the proposed formalism of smoothed and long-range second-moment approximation of tight-binding. Taking the potential as the starting base, an atomistic computation/simulation route was developed for designing favored and optimized compositions for Mg-Cu-Y metallic glass formation. Simulations revealed that the physical origin of metallic glass formation is the collapse of crystalline lattice when solute concentration exceeds a critical value, thus leading to predict a hexagonal region in the Mg-Cu-Y composition triangle, within which metallic glass formation is energetically favored. It is proposed that the hexagonal region can be defined as the intrinsic glass formation region, or quantitative glass formation ability of the system. Inside the hexagonal region, the driving force for formation of each specific glassy alloy was further calculated and correlated with its forming ability in practice. Calculations pinpointed the optimized stoichiometry in the Mg-Cu-Y system to be Mg64Cu16Y20, at which the formation driving force reaches its maximum, suggesting that metallic glasses designed to have compositions around Mg64Cu16Y20 are most stable or easiest to obtain. The predictions derived directly from the atomistic simulations are supported by experimental observations reported so far in the literature. Furthermore, Honeycutt-Anderson analysis indicated that pentagonal bipyramids (although not aggregating to form icosahedra) dominate in the local structure of the Mg-Cu-Y metallic glasses. A microscopic picture of the medium-range packing can then be described as an extended network of the pentagonal bipyramids, entangled with the fourfold and sixfold disclination lines, jointly fulfilling the space of the metallic glasses.