This thesis is a theoretical study of configurational disorder in icosahedral boron-rich solids, in particular boron carbide, including also the development of a methodological framework for treating configurational disorder in such materials, namely superatom-special quasirandom structure (SA-SQS). In terms of its practical implementations, the SA-SQS method is demonstrated to be capable of efficiently modeling configurational disorder in icosahedral boron-rich solids, whiles the thermodynamic stability as well as the properties of the configurationally disordered icosahedral boron-rich solids, modeled from the SA-SQS method, can be directly investigated, using the density functional theory (DFT). In case of boron carbide, especially B4C and B13C2 compositions, the SA-SQS method is used for modeling configurational disorder, arising from a high concentration of low-energy B/C substitutional defects. The results, obtained from the DFT-based calculations, demonstrate that configurational disorder of B and C atoms in boron carbide is not only thermodynamically favored at high temperature, but it also plays an important role in altering the properties of boron carbide − for example, restoration of higher rhombohedral symmetry of B4C, a metal-to-nonmetal transition and a drastic increase in the elastic moduli of B13C2. The configurational disorder can also explain large discrepancies, regarding the proper- ties of boron carbide, between experiments and previous theoretical calculations, having been a long standing controversial issue in the field of icosahedral boron- rich solids, as the calculated properties of the disordered boron carbides are found to be in qualitatively good agreement with those, observed in experiments. In order to investigate the configurational evolution of B4C as a function of temperature, beyond the SA-SQS level, a brute-force cluster-expansion method in combination with Monte Carlo simulations is implemented. The results demonstrate that configurational disorder in B4C indeed essentially takes place within the icosahedra in a way that justifies the focus on lowenergy defect patterns of the superatom picture. The investigation of the thermodynamic stability of icosahedral carbon-rich boron carbides beyond the believed solubility limit of carbon (20 at.% C) demonstrates that, apart from B4C generally addressed in the literature, B2.5C represented by B10Cp2(CC) is predicted to be thermodynamically stable with respect to B4C as well as pure boron and carbon under high pressure, ranging between 40 and 67 GPa, and also at elevated temperature. B2.5C is expected to be metastable at ambient pressure, as indicated by its dynamical and mechanical stabilities at 0 GPa. A possible synthesis route of B2.5C and a fingerprint for its characterization from the simulations of x-ray powder diffraction pattern are suggested. Besides modeling configurational disorder in boron carbide, the SA-SQS method also opens up for theoretical studies of new alloys between different icosahedral boron-rich solids − for example, (B6O)1−x(B13C2)x and B12(As1−xPx)2. As for the pseudo-binary (B6O)1−x(B13C2)x alloy, it is predicted to display a miscibility gap resulting in B6O-rich and either ordered or disordered B13C2-rich domains for intermediate global compositions at all temperatures up to melting points of the materials. However, some intermixing of B6O and B13C2 to form solid solutions is also predicted at high temperature. A noticeable mutual solubility of icosahedral B12As2 and B12P2 in each other to form B12(As1−xPx)2 disordered alloy is predicted even at room temperature, and a complete closure of a pseudo-binary miscibility gap is achieved at around 900 K. Apart from B12(As1−xPx)2, the thermodynamic stability of other compounds and alloys in the ternary B-As-P system is also investigated. For the binary B-As system, zincblende BAs is found to be thermodynamically unstable with respect to icosahedral B12As2 and gray arsenic at 0 K and increasingly so at higher temperature, indicating that BAs may merely exist as a metastable phase. This is in contrast to the binary B-P system, in which zinc-blende BP and icosahedral B12P2 are both predicted to be stable. Owing to the instability of BAs with respect to B12As2 and gray arsenic, only a tiny amount of BAs is predicted to be able to dissolve in BP to form BAs1−xPx disordered alloy at elevated temperature. For example, less than 5% BAs can dissolve in BP at 1000 K. As for the binary As-P system, As1−xPx disordered alloys are predicted at elevated temperature − for example, a disordered solid solution of up to ∼75% As in black phosphorus as well as a small solubility of ∼1% P in gray arsenic at 750 K, together with the presence of miscibility gaps. The thermodynamic stability of three different compositions of α-rhombohedral boron-like boron subnitride, having been proposed so far in the literature, is investigated. Those are, B6N, B13N2, and B38N6, represented respectively by B12(N-N), B12(NBN), and [B12(N-N)]0.33[B12(NBN)]0.67. It is found that, out of these sub- nitrides, only B38N6 is thermodynamically stable from 0 GPa up to ∼7.5 GPa, depending on the temperature, and is thus concluded as a stable composition of α-rhombohedral boron-like boron subnitride.