The energetically favorable spatial configuration of M3+ ions and oxide-ion vacancies near a symmetrical grain boundary (GB) in cubic zirconia is determined for various trivalent species M3+ (M = Al, Sc, Y, Gd, La), and the driving force for grain boundary segregation (GBS) quantitatively examined using atomistic Monte Carlo simulations in conjunction with static lattice calculations. For a high concentration of similar to 10 mol %, it is found that point defects near a GB plane preferentially occupy specific sites to minimize total lattice energy, rather than being randomly distributed. Systematic analysis shows that energetically stable configurations of segregants vary depending on their ionic radii. Analysis of the driving force for GBS as a function of dopant concentration reveals that three important factors govern GBS. First, occupation of specific sites by point defects is necessary to minimize the total lattice energy; enrichment of point defects near the GB plane with random configuration does not decrease the total lattice energy significantly because of strong Coulombic interactions. Second, the factors governing GBS change with increasing dopant concentration. At dilute concentrations, relief of bond strain is the dominant factor, while at high concentrations Coulombic interactions, which depend strongly on the specific arrangement of defects, become another dominant factor. Third, the stabilization of matrix cations, Zr4+ ions, is the dominant factor to lower the driving force for GBS at all concentrations. In contrast, the stabilization of M3+ ions does not necessarily contribute to GBS of point defects at high concentrations. These findings suggest practical ways to control GBS to enhance materials' properties or minimize detrimental effects.