The energetically favorable spatial configuration of M(3+) ions and oxide-ion vacancies near a symmetrical grain boundary (GB) in cubic zirconia is determined for various trivalent species M(3+) (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 ∼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, Zr(4+) ions, is the dominant factor to lower the driving force for GBS at all concentrations. In contrast, the stabilization of M(3+) 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.