Density Functional Theory for Buckyballs within Symmetrized Icosahedral Basis

We have developed a highly efficient computation method based on density functional theory (DFT) within a set of fully symmetrized basis functions for the C60 buckyball, which possesses the icosahedral (Ih) point-group symmetry with 120 symmetry operations. We demonstrate that our approach is much more efficient than the conventional approach based on three-dimensional plane waves. When applied to the calculation of optical transitions, our method is more than one order of magnitude faster than the existing DFT package with a conventional plane-wave basis. This makes it very convenient for modeling optical and transport properties of quantum devices related to buckyball crystals. The method introduced here can be easily extended to other fullerene-like materials.

Buckyballs can also form single crystals and they have decent mobility for device applications [13,14]. Many fullerenes can display superconductivity at relatively high temperatures. It was observed that Cs-doped C 60 single crystals have a superconducting transition temperature at T C = 40 K [15] and several alkali metal-doped C 60 compounds exhibit T C in the range between 19 K and 47 K [16]. The existence of high-Tc superconductivity in fullerenes is likely caused by the strong electron-phonon interaction, but detailed microscopic theory for understanding this is still not available. The theoretical development for such an important problem is mainly hindered by the complexity of the system, which requires heavy computation and highly sophisticated theoretical analyses. Thus, an effort to significantly reduce the computation effort for calculating the electronic states in such systems is warranted.
For the C 60 buckyball, which is a truncated icosahedron [3,17,18], it is natural to choose a basis set consisting of products of spherical harmonic functions [Y lm (Ω)] and localized basis functions along the radial direction. The point group theory can be used to find proper linear combinations of spherical harmonics to form symmetry-adapted basis functions (SABFs) which transform according to irreducible representations (IRs) of the Since we aim to demonstrate the usefulness of the symmetrized icosahedral basis to facilitate the DFT calculation of high-symmetry systems, such as C60 and related crystals, we choose to start with a simpler scheme.
In the GTH approach, ( ) can be well described by a simple analytic form [26] ( ) = − erf ( ) + ∑ (4) where erf denotes the error function and is the ionic charge. and are fitting parameters for the C atom. The Fourier transform of ( ) is also given by a simple form [26] ( ) = − + ∑ / / (5) where is the sample volume.
( ) in Equation (2) denotes the self-consistent Hartree potential and the last term in Equation (1) denotes the exchange-correlation potential, which is deduced from the Monte Carlo results calculated by Ceperley and Alder [28] and parametrized by Perdew and Zunger [29]. The DFT effective potential is determined selfconsistently until its root-mean-square change is less than 10 −6 Ry. The nonlocal pseudopotential ( ) in the GTH approach is given by where Figure 1. Schematic diagram of buckyball [25].
Nanomaterials 2023, 13,1912 3 of 17 In density functional theory (DFT), the Kohn-Sham Hamiltonian for an electron in the C 60 molecule is written as: where we have adopted the atomic units throughout the paper with energy measured in Rydberg (Ry) and distance in bohr. In the right-hand side of Equation (1), the first, second, and third term describes the kinetic energy, the local pseudopotential, and the nonlocal pseudopotential, respectively. The local pseudopotential consists of three terms: V loc (r) = V ion (r) + V H (r) + V xc (r) (2) where the first term describes the ionic local potential with in which τ σ denotes the position of different C atoms in the buckyball. For simplicity, we adopt the norm-conserving pseudopotential (NCPP) developed by Goedecker, Teter, and Hutter (GTH) [26]. The ultrasoft pseudopotential (USPP) developed by Vanderbilt [27] can also be adopted in our current approach. However, the implementation of the projector augmented part used to reduce the number of plane waves (or spherical harmonics here) in the basis to achieve faster convergence in the calculation will require more effort. Since we aim to demonstrate the usefulness of the symmetrized icosahedral basis to facilitate the DFT calculation of high-symmetry systems, such as C 60 and related crystals, we choose to start with a simpler scheme.
In the GTH approach, V σ I (r) can be well described by a simple analytic form [26] V I (r) = − Z ion r erf( √ α 0 r) + ∑ 3 j=0 D j r 2j e −α 0 r 2 (4) where erf denotes the error function and Z ion is the ionic charge. α 0 and D j are fitting parameters for the C atom. The Fourier transform of V I (r) is also given by a simple form [26] where V is the sample volume. V H (r) in Equation (2) denotes the self-consistent Hartree potential and the last term in Equation (1) denotes the exchange-correlation potential, which is deduced from the Monte Carlo results calculated by Ceperley and Alder [28] and parametrized by Perdew and Zunger [29]. The DFT effective potential is determined self-consistently until its root-mean-square change is less than 10 −6 Ry. The nonlocal pseudopotential (V nl in the GTH approach is given bŷ for l = 0.1 and |m| ≤ l. s = 1.2 denotes different β functions used for each angular momentum, l. The normalized p l s (r) functions have simple analytic forms as given in [26]. h l s are energy parameters. The projection of β s lm function in the wave-vector space reads which can also be expressed in a simple analytical form [26].

Generation of Symmetrized Angular Basis Functions
To take advantage of the full point-group symmetry, we first construct the C 60 -adapted symmetrized angular basis functions via a suitable linear combination of the spherical harmonics Y lm (r). We define the symmetrized angular functions (called icosahedral harmonics) with symmetry type ( Γ, v) compatible with angular momentum l as where and D (l) Here, Λ denotes the 120 symmetry operations in the I h group and n(Γ) is the dimension of the irreducible representation Γ and h is the order of the point group. The index v labels the degenerate partners of an IR. One straightforward way to find the symmetrization coefficients C Γv lm is to use the projection method, commonly described in textbooks [30]. However, it requires the knowledge of transformation matrices Γ i (Λ) of the spherical harmonics for each symmetry operator Λ associated with the i-th IR of the point group. For I h group, the above procedure can be quite tedious. Although these coefficients for the I h group can be obtained from the projection method [31], it still requires a lot of effort to implement these coefficients in the current DFT code.
Here, we adopt a more practical approach. We obtain these symmetrization coefficients via diagonalization of an effective Hamiltonian matrix within a minimum basis set of the form { e −αr Y lm (Ω)} for each fixed index of angular momentum, l. The effective potential is taken to be of the form: where τ σ denotes the positions of 60 carbon atoms in C 60 , V I (r ) can be any model potential.
Here, we simply choose V I (r ) to be the GTH atomic local pseudopotential adopted in the current approach. The matrix element of V e f f (r) in the subspace { Y lm (Ω)} with a fixed l is given by Obviously, V e f f obeys the same point group symmetry as C 60 . Thus, the eigenvectors of V e f f defined in the subspace { Y lm (Ω)} with a fixed l will transform according to the IRs of the I h group. From the degeneracy (g) of the corresponding eigenvalues, we can identify the possible irreducible representations. For example, if g = 1, 4, 5, the states must belong to the Γ 1 (A), Γ 4 (G), and Γ 5 (H) representations, respectively. If g = 3, the states must belong to either the Γ 2 (F 1 ) or Γ 3 (F 2 ) representation. To pin down the precise symmetry type of each eigenstate, we simply evaluate the coupling matrix element between the state and a known basis state belonging to the possible representations as given in [32,33]. When the matrix element is nonzero, the symmetry type is identified. The above method can be easily applied to systems with any other point group. The symmetrized basis functions with lowest l for the five irreducible representations of the I h group with either even or odd parity are listed in Table 1. Table 1. The symmetrized basis functions with the lowest l for the five irreducible representations of the I h group with either even (g) or odd (u) parity.

The Choice of Basis Functions
To simplify the computation effort, we choose the B-spline [34] augmented icosahedral harmonics (BAIHs) as basis functions for the DFT calculation of C 60 buckyball. The BAIHs are defined as with where B n (r) denotes the B-spline functions of a suitable order [29].
where I nl (k) = drj l (kr)r 2 B n (r) (19) For each fixed B-spline basis B n (r), we can take linear combinations of Y lm (Ω) to construct symmetrized basis functions that transform according to the irreducible representations of an icosahedral group, I (or I h = i ⊗ I, where i denotes inversion).
To make sure that the degenerate partners in each irreducible representation transform in the same way as in the corresponding basis functions for lower l, we do the following. Let K Γv l 1 (Ω) denote the symmetrized basis function for the v-th degenerate partner for the Γ representation, and l 1 denotes the lowest possible l for this representation. We can take linear combinations of states |l; Γv such that Thus, Z matrix is proportional to the inverse of the n(Γ)-dimensional matrix with ( v, v ) elements given by l, Γv H e f f l 1 , Γv . Using this simple procedure, we can obtain coefficients for symmetrized basis states for all higher l needed. With this approach, we obtained the symmetry coefficients C Γv lm for the 5 irreducible representations Γ 1 (A), Γ 2 (F 1 ), Γ 3 (F 2 ), Γ 4 (G), and Γ 5 (H) of I h group for l up to 55.

The Matrix Elements of Local Pseudopotential
Since V ion (r) should be invariant under all symmetry operations of the I h group, we where For the ideal buckyball, the magnitude of τ σ is the same for all 60 carbon atoms. Note that there is an overflow problem in j l 2iα j τ σ r when α j is large, since where we have used j 1 (z) = sin(z)/z/z − cos(z)/z.
Using the recursion relation for spherical Bessel functions, we have independent of n. V H (r) = d 3 r ρ(r )/|r − r | denotes the Hartree potential, where ρ(r) denotes the charge density. Since the wavefunctions are expanded in terms of BAIHs, we can write Here, the charge density transforms according to the Γ 1g representation. Thus, we L (Ω)ρ(r). Due to the high symmetry, the integral only has to cover 1/120 of the whole solid angle as shown by half the surface area enclosed by green lines in Figure 1. Using the expansion where r < = min(r, r )and r > = max(r, r ), and Here, V xc (r) denotes the exchange-correlation potential. In each iteration, we shall expand V xc (r) in terms of spherical harmonics in the form where the integral can be done efficiently by using the symmetry property.
Then, the matrix elements read and LM can be evaluated for all L (up to 55) before the DFT calculation.

Nonlocal Pseudopotential
The matrix elements ofV nl within the BAIH basis can be evaluated efficiently due to its separable form. We have where I nl (k) is given in Equation (19) and (9). A total of 18 B-splines, defined over a range of 12 a.u., are used to expand the radial wavefunction. The maximum angular moment L of the symmetry-adapted basis function (SABF) used is 45. The numbers of SABFs for these 10 IR's are listed in Table 2. Due to the small dimension of the Hamiltonian matrix for each symmetry type considered, the diagonalization of the Hamiltonian can be done efficiently with a direct solver.  spline  18  18  18  18  18  18  18  18  18  18   SABF  23  13  46  60  46  60  69  72  92  84   Total  414  234  828  1080  828  1080  1242  1296 1656 1512

Energy Levels
We found that there are 32 distinct energy levels (not including the degeneracy factor) for the occupied levels (labeled by E v,j for the C 60 molecule. Table 3 displays the number of occurrences (n) for each symmetry type or the occupied energy levels. We also list the corresponding values of n × g in the last row for each symmetry type to account for the level degeneracy ( g). Note that the sum of n × g for all occupied levels is equal to 120 since each C atom contributes two electrons to the occupied levels (or valence states). Table 3. The number of occurrences (n) and the number of basis states (n × g ) for each symmetry type for occupied energy levels. 3  0  3  12  3  12  16 16 35 20 For comparison purposes, we also performed calculations of the C 60 buckyball by using the Quantum ESPRESSO (QE) plane-wave-based package [35] with the NCPP option; the exchange-correlation functional parametrized by Perdew and Zunger [29] (same as the one used in the current approach) was used. In the QE calculation, a cubic supercell with 20Å along each side is chosen for the calculation. Thus, we are modeling an artificial buckyball crystal with the QE approach, instead of a single molecule as considered in our current code. We have checked the suitability of the vacuum length used in the QE calculation and found that the results concerned here are not significantly altered when the cell size is varied between 15Å and 20Å. Since C 60 is not electrically polarized, based on previous calculations on graphene nanoribbons [36], a vacuum space of~10Å is enough. So, a cell length of 16-20Å for C 60 is typical (10-14Å of vacuum + 6Å for the buckyball). The energy cutoff of 70 Ry was used (typical for the NCPP adopted).
We show the comparison of results obtained by the current method and those by using both the QE package [35] and Gaussian 16 package [20] for the inter-level energy spacings of the highest 20 occupied (valence) levels and lowest 3 unoccupied (conduction) levels in Table 4 The HOMO (highest occupied molecular orbital) level obtained by the current calculation is at −2.71 eV. We define the energy-level differences ∆c j = E c,j+1 − E c,j and ∆v j = E v,j+1 − E v,j , for the inter-level energy spacings between two consecutive unoccupied (conduction) levels and occupied (valence) levels, respectively, while the band gap energy is given by E g = E c,1 − E v,1 . Here, j = 1 denotes the topmost valence (lowest conduction) energy level, and higher j indicates energetically decreasing (increasing) levels for valence (conduction) levels. The corresponding symmetry types (IRs) are also indicated in parentheses.
As can be seen from Table 4, the results obtained for the bound states in the molecule (which include all occupied levels and some low-lying conduction levels) by the present method agree quite well with those results obtained by QE and reasonably well with Gaussian 16. In the calculation with Gaussian orbitals, we selected the VWN option (also within the local density approximation) [37,38]. We note that the Gaussian package uses the all-electron approach instead of the pseudopotential method. Thus, there is more deviation between our results and Gaussian 16 results. For the high-lying conduction states (which correspond to unbound states of the molecule), we see some deviation of our results from those obtained by QE and Gaussian 16. Table 4. Energy level differences (in eV) of C 60 molecules. See text for details. E g means the energy difference between the lowest unoccupied molecular orbital (LUMO) (with symmetry Γ 2u ) and the highest occupied molecular orbital (HOMO) (with symmetry Γ 5u ). 33 This discrepancy is due to the different boundary conditions imposed in different approaches. In our approach and Gaussian 16, we consider only an isolated C 60 molecule, while an artificial C 60 solid was considered in the supercell approach used by QE. For the unbound conduction states, there will be a strong overlap between states derived from neighboring C 60 molecules in the artificial solid. Thus, these states will have significant dispersion. Namely, these energy levels are k-dependent, where k is the wavevector of the C 60 solid. For the current approach, we choose a finite range for the B-spline basis functions. This effectively introduces a quantum barrier for these unbound states. The discrete levels we obtained for these unbound states represent a discretized sampling of the continuum states. Thus, the energy spacing of these unbound states will depend on the range of the knot sequence chosen for the B-spline basis functions. However, this finite-size sampling of continuum states can still give a reasonable description of the optical excitation spectrum even into the continuum region, as long as the energy spacings are smaller than the line broadening used to mimic the absorption spectra. That means the sampling is dense enough to capture the main features of the optical excitation spectrum. For Gaussian 16, there is no rigid boundary. However, the energy spacings between unbound states depend sensitively on the number of Gaussian orbitals chosen and the values of exponents used.

Transitions BAIH QE Gaussian
To do a bench-mark comparison in computation speed we run both the current code and the QE package with a single processor. The CPU time needed to calculate all eigenstates for the C 60 molecule is~300 s with the current code, while it would take~1000 s to get only the 120 occupied levels by QE. We note that the diagonalization procedure in QE was done via the conjugate gradient (CG) method. To study the optical properties of C 60 , we need to include many conduction states. If we also calculate 300 unoccupied levels by QE, the CPU time needed will increase to~20 h (on a single processor). Therefore, our method is more than two orders of magnitude faster than the well-optimized QE package for such an application. With further optimization, the current code can be made even more efficient.

Optical Absorption Spectrum
In this section, we calculate the optical absorption spectrum of the C 60 buckyball, while neglecting the excitonic effect and compare it with the corresponding results obtained by QE. The optical absorption spectrum is proportional to the imaginary part of the dielectric response function as given by [39]

Optical Absorption Spectrum
In this section, we calculate the optical absorption spectrum of the C60 buckyball, while neglecting the excitonic effect and compare it with the corresponding results obtained by QE. The optical absorption spectrum is proportional to the imaginary part of the dielectric response function as given by [39] ( ħ ) ∝ ∑ |⟨ , | • | , ⟩ , | ħ + , − , .
Here denotes the polarization vector of the photon, | , ⟩ and | , ⟩ denote the ith valence (occupied) state with energy , and j-th conduction (unoccupied) state with energy , , respectively. denotes the momentum operator, and ħ is the photon energy. Since the symmetry types of the eigenstates obtained by the current code are already known, we can apply the selection rules and significantly reduce the computation effort for calculating the dielectric response function, (ħ ). In this work, all eigenstates of the C60 buckyball are localized functions, and it is convenient to use the commutator relation = i ℏ [ , ] and obtain Let  ( ) and  ( ) denote the molecular orbitals (MOs) in a C60 buckyball with symmetry type  and  for occupied and unoccupied levels, respectively.
The photon transforms like = 1 spherical harmonics which has symmetry under the group. The selection rule imposed by applying the group theory indicates that the transition  →  is forbidden when the vector coupling coefficient

Optical Absorption Spectrum
In this section, we calculate the optical absorption spectrum o while neglecting the excitonic effect and compare it with the corre tained by QE. The optical absorption spectrum is proportional to t the dielectric response function as given by [39] ( ħ ) ∝ ∑ |⟨ , | • | , ⟩ , | ħ + , − Here denotes the polarization vector of the photon, | , ⟩ an th valence (occupied) state with energy , and j-th conduction (un energy , , respectively. denotes the momentum operator, and ergy. Since the symmetry types of the eigenstates obtained by the ready known, we can apply the selection rules and significantly red effort for calculating the dielectric response function, (ħ ).
In this of the C60 buckyball are localized functions, and it is convenient to relation = i ℏ [ , ] and obtain The photon transforms like = 1 spherical harmonics which has the group. The selection rule imposed by applying the group the transition  →  is forbidden when the vector c Hereê denotes the polarization vector of the photon, |v, i and |c, j denote the i-th valence (occupied) state with energy E v,i and j-th conduction (unoccupied) state with energy E c,j , respectively. p denotes the momentum operator, and

Optical Absorption Spectrum
In this section, we calculate the optical absorption spectr while neglecting the excitonic effect and compare it with the c tained by QE. The optical absorption spectrum is proportiona the dielectric response function as given by [39] ( ħ ) ∝ ∑ |⟨ , | • | , ⟩ , | ħ + Here denotes the polarization vector of the photon, | , th valence (occupied) state with energy , and j-th conductio energy , , respectively. denotes the momentum operator, ergy. Since the symmetry types of the eigenstates obtained b ready known, we can apply the selection rules and significantl effort for calculating the dielectric response function, (ħ ).
In of the C60 buckyball are localized functions, and it is convenie relation = i ℏ [ , ] and obtain The photon transforms like = 1 spherical harmonics which the group. The selection rule imposed by applying the group transition  →  is forbidden when the vecto Ω is the photon energy. Since the symmetry types of the eigenstates obtained by the current code are already known, we can apply the selection rules and significantly reduce the computation effort for calculating the dielectric response function, ε 2 (

Optical Absorption Spectrum
In this section, we calculate the optical absorption spectrum of the C60 buckyb while neglecting the excitonic effect and compare it with the corresponding results tained by QE. The optical absorption spectrum is proportional to the imaginary part the dielectric response function as given by [39] ( ħ ) ∝ ∑ |⟨ , | • | , ⟩ , | ħ + , − , . ( Here denotes the polarization vector of the photon, | , ⟩ and | , ⟩ denote th th valence (occupied) state with energy , and j-th conduction (unoccupied) state w energy , , respectively. denotes the momentum operator, and ħ is the photon ergy. Since the symmetry types of the eigenstates obtained by the current code are ready known, we can apply the selection rules and significantly reduce the computat effort for calculating the dielectric response function, (ħ ).
In this work, all eigensta of the C60 buckyball are localized functions, and it is convenient to use the commuta relation = i ℏ [ , ] and obtain

Optical Absorption Spectrum
In this section, we calculate the optical absorption spectrum of the C60 buckyball, while neglecting the excitonic effect and compare it with the corresponding results obtained by QE. The optical absorption spectrum is proportional to the imaginary part of the dielectric response function as given by [39] ( ħ ) ∝ ∑ |⟨ , | • | , ⟩ , | ħ + , − , .
Here denotes the polarization vector of the photon, | , ⟩ and | , ⟩ denote the ith valence (occupied) state with energy , and j-th conduction (unoccupied) state with energy , , respectively. denotes the momentum operator, and ħ is the photon energy. Since the symmetry types of the eigenstates obtained by the current code are already known, we can apply the selection rules and significantly reduce the computation effort for calculating the dielectric response function, (ħ ).
In this work, all eigenstates of the C60 buckyball are localized functions, and it is convenient to use the commutator relation = i ℏ [ , ] and obtain Let  ( ) and  ( ) denote the molecular orbitals (MOs) in a C60 buckyball with symmetry type  and  for occupied and unoccupied levels, respectively.
The photon transforms like = 1 spherical harmonics which has symmetry under the group. The selection rule imposed by applying the group theory indicates that the transition  →  is forbidden when the vector coupling coefficient

Optical Absorption Spectrum
In this section, we calculate the optical absorption spectrum of while neglecting the excitonic effect and compare it with the corresp tained by QE. The optical absorption spectrum is proportional to the the dielectric response function as given by [39] Here denotes the polarization vector of the photon, | , ⟩ and th valence (occupied) state with energy , and j-th conduction (unoc energy , , respectively. denotes the momentum operator, and ħ ergy. Since the symmetry types of the eigenstates obtained by the c ready known, we can apply the selection rules and significantly reduc effort for calculating the dielectric response function, (ħ ).
In this w of the C60 buckyball are localized functions, and it is convenient to us relation = i Let ϕ Γ h v h i (r h ) and ϕ Γ e v e j (r e ) denote the molecular orbitals (MOs) in a C 60 buckyball with symmetry type Γ h v h and Γ e v e for occupied and unoccupied levels, respectively. The photon transforms like l = 1 spherical harmonics which has Γ 2u symmetry under the I h group. The selection rule imposed by applying the group theory indicates that the transition Γ h v h → Γ e v e is forbidden when the vector coupling coefficient C(Γ h v h , Γ 2u v, ; Γ e v e ) = 0. Here, the vector coupling coefficient plays the same role as the Clebesh-Gordan coefficient for the products of two spherical harmonics. Namely, for a given valence state with symmetry type Γ h , the final conduction state must belong to an IR, compatible with the symmetry of the product Γ h ⊗ Γ 2u . The selection rules can be worked out by using the group theory, and they are listed in Table 5. Table 5. Selection rules for dipole allowed transitions for a C 60 buckyball.

Occupied
Unoccupied Based on the Wigner-Ekart theorem [40], the dipole matrix elements v, i|r|c, j can be written as v, i|r|c, j = v, i||r||c, j C(Γ h v h , Γ 2u v, ; Γ e v e ) (31) where v, i||r||c, j is called the reduced matrix element, which is independent of the indices of the degenerate partners in the initial state ( v h ) and final state ( v e ), as well as the polarization of the photon (indexed by v). Thus, for each allowed transition between two manifolds, we only have to evaluate the reduced matrix element once and immediately obtain all related matrix elements with saving in computation time of n v × n c × 3 fold, where n v and n c are the dimensions of the Γ h and Γ e IRs, respectively. We have worked out the vector-coupling coefficients for the I h group based on group theory. The results are listed in the Appendix A.
In the calculation of the imaginary part of dielectric response function, ε 2 (Ω), we have introduced a broadening parameter Γ. Namely, the delta function in Equation (29)

Optical Absorption Spectrum
In this section, we calculate the optical absorption spectrum of the C60 buckyball while neglecting the excitonic effect and compare it with the corresponding results obtained by QE. The optical absorption spectrum is proportional to the imaginary part of the dielectric response function as given by [39] ( ħ ) ∝ ∑ |⟨ , | • | , ⟩ , | ħ + , − , .
Here denotes the polarization vector of the photon, | , ⟩ and | , ⟩ denote the ith valence (occupied) state with energy , and j-th conduction (unoccupied) state with energy , , respectively. denotes the momentum operator, and ħ is the photon energy. Since the symmetry types of the eigenstates obtained by the current code are already known, we can apply the selection rules and significantly reduce the computation effort for calculating the dielectric response function, (ħ ). In this work, all eigenstates of the C60 buckyball are localized functions, and it is convenient to use the commutator relation = i ℏ [ , ] and obtain Let  ( ) and  ( ) denote the molecular orbitals (MOs) in a C60 buckyball with symmetry type  and  for occupied and unoccupied levels, respectively The photon transforms like = 1 spherical harmonics which has symmetry under the group. The selection rule imposed by applying the group theory indicates that the transition  →  is forbidden when the vector coupling coefficient . ε 2 (Ω) of a C 60 buckyball calculated by the current method is shown in Figure 2

Optical Absorption Spectrum
In this section, we calculate the optical absorption spectrum of the C60 buckyball, while neglecting the excitonic effect and compare it with the corresponding results obtained by QE. The optical absorption spectrum is proportional to the imaginary part of the dielectric response function as given by [39] Here denotes the polarization vector of the photon, | , ⟩ and | , ⟩ denote the ith valence (occupied) state with energy , and j-th conduction (unoccupied) state with energy , , respectively. denotes the momentum operator, and ħ is the photon energy. Since the symmetry types of the eigenstates obtained by the current code are already known, we can apply the selection rules and significantly reduce the computation effort for calculating the dielectric response function, (ħ ). is forbidden when the vector coupling coefficient Ω < 6 eV. For photon energies higher than 6 eV, we still get similar spectral features with roughly the same average oscillator strengths, but the details are somewhat different. The main deviation is due to the difference in boundary conditions used between our approach and the QE package. Here, we consider an isolated C 60 buckyball confined by an infinite potential at a radius of 12 bohrs (imposed by the cutoff of B-spline functions), while the QE package adopts a supercell with periodic boundary condition. Different boundary conditions will lead to different dielectric response functions at high photon energies [41,42]. Since all the eigenstates of well-defined symmetry have been obtained, the computation of the dipole strength based on Equation (31) can be calculated very efficiently (with less than 10 seconds) when the selection rule and Wigner-Ekart theorem are adopted. If we perform the calculation by a brute-force method without considering the symmetry, it will take much longer. The same concept can be applied to the calculation of the excitonic effect for the C 60 buckyball by solving the Bethe-Salpeter equation [21]. It is expected that the use of symmetrized basis can also speed up the computation significantly in comparison to the brute-force method.

Conclusions
Using the symmetrized-basis approach, we have implemented a highly efficient DFT code for the C60 buckyball. The energy levels calculated by this method are in close agreement with those obtained by using the Quantum Espresso (QE) package and in fair agreement with results obtained by the all-electron calculation with Gaussian 16 package. The computation time needed to obtain the self-consistent charge density of C60 buckyball is about 1/3 of that by using the QE package. Note that our code is not yet fully optimized, and therefore, it has the potential to speed up further. Once the self-consistent charge den-

Optical Absorption Spectrum
In this section, we calculate the optical absorption spectrum while neglecting the excitonic effect and compare it with the corre tained by QE. The optical absorption spectrum is proportional to t the dielectric response function as given by [39]

Conclusions
Using the symmetrized-basis approach, we have implemented a highly efficient DFT code for the C 60 buckyball. The energy levels calculated by this method are in close agreement with those obtained by using the Quantum Espresso (QE) package and in fair agreement with results obtained by the all-electron calculation with Gaussian 16 package.
The computation time needed to obtain the self-consistent charge density of C 60 buckyball is about 1/3 of that by using the QE package. Note that our code is not yet fully optimized, and therefore, it has the potential to speed up further. Once the self-consistent charge density is obtained, the computation of 120 occupied levels and 300 unoccupied levels with the current code takes only about 300 s, while it would take more than 100 times longer to do the same by using the QE package. For the calculation of the optical excitation spectrum from the 120 occupied levels to 300 unoccupied levels (not including the excitonic effect), the CPU time needed is less than 10 s after obtaining the eigenstates.
The method can be readily extended to other fullerenes such as C 70 and C 80 within the geometry-adapted symmetrized basis set (with reduced numbers of symmetry operations). Here, we have introduced a simple scheme to utilize the computation method to extract the coefficients in the symmetrized basis that transform according to the IRs of the underlying point group as described in Section 2.3. Thus, applying the same idea to other fullerenes can be conveniently achieved. This method can also be extended to study fullerene-like crystals and fullerene-related quantum devices. For such application, we will calculate the coupling matrix describing the interaction of the fullerene with neighboring objects based on the first principles in the framework of the linear combination of molecular orbitals (LCMO) approach. Since we use symmetry-adapted harmonics augmented by B-splines as basis functions, all eigenfunctions of the systems considered are localized near the fullerene surface, this can be done conveniently, and the resulting Hamiltonian matrix will be sparse.
For optoelectronic properties, the exciton plays a significant role. A model calculation of excitonic states in C 60 crystals based on the LCMO method has been reported, in which the MOs are deduced from a DFT-GW calculation [43]. The overlap integrals of MOs between two adjacent C 60 molecules have been neglected. Here, with the use of icosahedral harmonics augmented by B-splines as basis functions, the intra-molecular optical transition matrix and electron-hole Coulomb scattering matrix can be computed with very little effort. Thus, we can calculate the effect of inter-molecular overlap integrals on the excitonic states efficiently with the current approach.
Author Contributions: C.-Y.R. did the group-theory analysis for the selection rule and vectorcoupling coefficients, and wrote the code for the current DFT calculation. R.K.P. assisted in finding the initial set of icosahedral harmonics and carried out DFT calculations based on the QE and Gaussian 16 packages for comparison purposes. Y.-C.C. initiated the idea, developed the underlying formulation, wrote the code for the part to generate the symmetrized basis by a model Hamiltonian, and supervised the whole project. C.-Y.R. and Y.-C.C. wrote the major part of the manuscript and R.K.P. proofread it and provided corrections. All authors have read and agreed to the published version of the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.