Efficient Method for Twist-Averaged Coupled Cluster Calculation of Gap Energy: Bulk Study of Stannic Oxide

We study gap energy of the semiconducting oxide SnO 2 through ab-initio calculations. DFT and coupled cluster calculations are presented and discussed. The efficiency of twist averaging in reducing the finite-size errors is tested through different functionals. We report an overestimation in the gap energy by finite-size scaling at the thermodynamic limit through equation-of-motion (EOM) CCSD calculations. To address one-body and many-body errors, we report a combination of the Kwee-Zhang-Krakauer (KZK) approach with twist averaging to explain twist averaged EOM-CCSD gap energy. In this approach, the correction to the gap energy originates from the difference between mean-field and many-body approximations and at the end the difference is added to the mean-filed gap of an infinite system to estimate the many-body gap in the thermodynamic limit. Unique single twist angles are introduced that yield cost-effective and accurate energies compared to the full twist averaging.


Introduction
SnO2 in the category of post-transition-metal oxides is used as transparent conducting oxides (TCOs) that are materials with specific features of high electrical conductivity and optical transparency in the visible range.These unique properties make SnO2 1 a material with a variety of applications including optoelectronic devices, 2 transparent transistors [3][4][5] and solid state gas sensors. 6,7 ue to such aspects, the electrical and optical properties of SnO2 have been subjects of many experimental and theoretical studies.Accordingly, owing to a correlation between transparency and band gap of SnO2, a detailed understanding of its electronic structure is essential.
In interpretation of adsorption measurements in SnO2, there have been controversies regarding the symmetry of top valence band and values of transitions from the valence band to the conduction band.To clear up the argument, Nagasawa et al, 8 Agekyan, 9 and Fröhlich et al. 10 established a direct-forbidden minimum gap of 3.6 eV at Γ, where the valence band maximum (VBM) had Γ !" symmetry.Also, experimental findings [11][12][13] and calculations 14 show conduction band minimum (CBM) has Γ # " symmetry.Therefore, there is a dipole forbidden transition between VBM and CBM.Although dipole transitions are symmetry forbidden in single photon optical absorption, there are dipole-allowed transitions in two-photon absorption, where the measured value of ~3.6 eV has been reported. 10,15 the study of SnO2, local density approximation (LDA), generalized gradient approximation (GGA), meta-GGAs, or van der Waals corrected functionals are known to underestimate the binding energy of the Sn d states.Hence, the calculated gaps based on the LDA [16][17][18][19][20][21][22][23] and the GGA 24- 26 are assumed to suffer from theoretical deficiencies.High level theories, e.g., hybrid functionals or GGA+U are known [27][28][29] to correct the position of the d states.Nonetheless, many factors such as the amount of exact exchange in hybrid functionals and the underlying GGA along with the value of U parameter in GGA+U might affect the electronic structure properties.Due to the full d states in SnO2, GGA+U 23,[30][31][32][33] does not provide correct description of the band gap.][36][37] Forms of Heyd-Scuseria-Ernzerhof (HSE) short-range screened hybrid 38,39 that incorporate a fraction of Hartree-Fock (HF) exchange are well-stablished and in widespread use.These functionals yield an agreement between the https://doi.org/10.26434/chemrxiv-2023-dks4g-v3ORCID: https://orcid.org/0000-0002-8259-8637Content not peer-reviewed by ChemRxiv.License: CC BY-NC-ND 4.0 accuracy and computational cost.Besides, full-range hybrids such as B3PW91 have shown adequate efficiency in band gap prediction of semiconductors. 40][44][45] As for SnO2, HSE03+G0W0 yields a fundamental gap of 3.65 eV close to the experimentally approved optical gap. 21A behavior that is assigned to the dipoleallowed direct transitions in the vicinity of the VBM at Γ. Despite of desirable accuracy and computational cost in the GW approximation, the errors due to the underlying DFT XC functional or self-consistency are difficult to assess. 467][48][49] While growing in popularity, the main focus for the widespread adoption of these methods has been controlling finite-size errors so that accuracy is preserved using small simulation cells (instead of supercells).This is of note in using coupled cluster theory with periodic boundary conditions to model extended systems.[52] Here, we focus on the lowest-order truncation in the coupled-cluster (CC) theory of the Green's function, namely, equation-of-motion CC theory with single and double excitations (EOM-CCSD) 53,54 to estimate the band gap energy of the SnO2.Several studies reported band gaps of various systems using periodic EOM-CCSD calculations. 46,47,55 Te computational cost of CCSD scales with the system size and the number of the k-points (N $ ) which necessitates extrapolations of the estimations in the thermodynamic (or infinite-size) limit (TDL). 47This limitation is enhanced by using larger basis sets to yield accurate calculations.To this end, twist averaging 56 is https://doi.org/10.26434/chemrxiv-2023-dks4g-v3ORCID: https://orcid.org/0000-0002-8259-8637Content not peer-reviewed by ChemRxiv.License: CC BY-NC-ND 4.0 used to reduce the finite-size errors and facilitate the extrapolation to the TDL. 50[58][59][60] The twist averaging is known to generate a smoother extrapolation to the TDL. 51Here, we introduce unique single twist angles yielding cost-effective and accurate results compared to the full twist averaging.For the many-body errors in the correlation energy arising from finite simulation cells subject to periodic boundary conditions we consider Kwee-Zhang-Krakauer (KZK) approach. 61These errors are distinct from the errors due to incomplete k-points in the Brillouin zone integration.The KZK approach considers finite-size corrections obtained from DFT calculations that is applied through post-processing steps to manybody calculations.At DFT level, we use the B3LYP, HSE06, PBE0, 62 B3PW91, and Perdew-Burke-Ernzerhof (PBE) 63 functionals, where the performance of short-range and full-range hybrids are evaluated.A combination of the KZK method with twist averaging 57 is introduced to address both one-body and many-body errors in the band gap calculations.
The rest of the paper is organized as follows.A description of the methodology is given in Sec.II.
In Sec III we present and analyze our calculations which focus on the gap energy of the SnO2.The reliability and validity of the applied methods are discussed in this section.The conclusion is given in Sec IV.

Computational details
Tetragonal rutile unit cell of SnO2 includes six atoms and is represented by two lattice parameters a and c and the internal parameter of u. [64][65][66][67][68] Each tin atom is surrounded by a distorted octahedron  Within the periodic CC theory, we calculate the band gap using the ionization potential (IP-EOM-CCSD) and the electron affinity (EA-EOM-CCSD).Here, we apply twist averaged boundary conditions (TABC) in order to reduce the single-particle errors.The twists can be chosen from a uniform Monkhorst-Pack grid, Baldereschi points, and high symmetry points, where offsetting from Γ may make other choices.Indeed, an insightful choice of the twists can reduce the cost while preserving the accuracy. 73e other finite-size errors, known as many-body contribution based on the DFT calculations to mimic the EOM-CCSD many-body errors.Indeed, the correction to the gap energy originates from the difference between mean-field and many-body approximations.This difference is then applied to the mean-filed gap of an infinite system to estimate the many-body gap in the thermodynamic limit. 77 address one-body and many-body contributions, we consider a combination of the KZK method with twist averaging. 57Here, we explain the twist averaged EOM-CCSD band gap energy augmented with KZK technique in an infinite simulation cell as follows.

A. Finite size errors
G '()*+,-,(/) where, the first term is the twist averaged CCSD energy for a simulation cell of size L and the second term is the finite-size correction based on DFT calculations.FS represents finite-size.
Hence, one has in which, N is the twist number, G 29' (∞) is the DFT gap computed using a fully converged kpoints.The finite-size DFT gap, G 29',91 (L,  & ), indicates twist averaged DFT gap using the same simulation cell and twists used in the CCSD twist averaged calculations.Our finite-size correction addresses all the finite-size errors.

III.
Results and discussions
In the Table I A comparison of the DFT gap values given in the Table I using 2 × 2 × 2 and 3 × 3 × 3 samplings with the different SBKJC and DZVP basis sets for Sn atoms through all the functionals, yields about 0.11 eV difference.As we noted, we could not afford DFT calculations using the DZVP basis set for Sn atoms with sampling that is beyond 3 × 3 × 3. Thus, we performed the extrapolations at this step.The extrapolation has the form of N $ +#/! to reach the thermodynamic (or infinite-size) limit (TDL). 46,47  $ is the number of the k-points in the mesh.To do the extrapolation we excluded the 1 × 1 × 1 values.As is seen in the Table I, the TDL values using the PBE0, B3LYP, and B3PW91 functionals and the DZVP basis set for Sn atom yields about 0.8 eV higher gap energy compared to the corresponding using 3 × 3 × 3 k-points calculations.Increasing k-points results in a convergence to the TDL limit.A switch to SBKJ for Sn atom allows us to increase the mesh to 5 × 5 × 5, so the interpolations (TDL) yield a lower average difference of 0.31 eV compared to the latest calculations with 5 × 5 × 5 k-points.The TDL value (Extrap,  > +#/! ) for each of the PBE0, https://doi.org/10.26434/chemrxiv-2023-dks4g-v3ORCID: https://orcid.org/0000-0002-8259-8637Content not peer-reviewed by ChemRxiv.License: CC BY-NC-ND 4.0 B3LYP, and B3PW91 functionals, using different basis sets for tin atoms, is given in the Table II.
For each functional, the difference between the corresponding extrapolated values is less than 0.2 eV.In the Table I Our results indicate that TABC with higher sampling is compatible with TDL and as the number of twists increase to infinity, the accuracy of the result is increased.
Besides, our results confirm that addressing finite-size errors in the mean-field approaches (as DFT) can be easily achieved through twist averaging as the number of twists tends to infinity.The performance of the extrapolation to reach TDL is shown in Figure 2. The solid and dashed lines (red, green, and purple) correspond to the gap energies within DFT using DZVP and SBKJC for Sn atoms, respectively, in PySCF package.The vertical red, green, and black lines show the calculated gap values using the PBE0, B3LYP, and HSE0 functionals, respectively, in the VASP package.Details are given in the Table S2 in the SI.
https://doi.org/10.26434/chemrxiv-2023-dks4g-v3ORCID: https://orcid.org/0000-0002-8259-8637Content not peer-reviewed by ChemRxiv.License: CC BY-NC-ND 4.0 From the figure, it is seen that irrespective of the type of applied basis set, all the estimations lie in the same range, as discussed above.Ab initio EOM-CCSD calculations in periodic systems are very computationally expensive.Indeed, it is not practical to go beyond 3 × 3 × 3 k-points in our system.Using 3 × 3 × 3 k-points, the EOM-CCSD gap is 3.32 eV and the corresponding TDL EOM-CCSD is 4.36 eV, see Table I.Given as the blue dashed line in the Figure 2, the least-squares extrapolation (N $ +#/! ) to the TDL limit is obtained by using 2 × 2 × 2 and 3 × 3 × 3 values (the precise uncertainty cannot be measured by two points).The TDL EOM-CCSD value is overestimated by about 0.8 eV that is attributed to lack of EOM-CCSD calculations using higher k-points, due to too demanding computational cost, and the many-body finite-size errors that could not be addressed here.As a result, the finite-size errors remain large.In the next section, the correction to this type of errors is discussed.

Twist averaged EOM-CCSD band gap energy augmented with KZK technique
As we discussed earlier, in order to reduce or eliminate the finite-size errors, it is important to perform the calculations using a supercell and a dense mesh.However, these solutions may not be achievable and generalized.In EOM-CCSD, one should notice that an increase in the k-points sampling results in a slower convergence and the cost exceeds the computational resources.For the same reason, it is not always practical to utilize the supercell for EOM-CCSD calculations to increase the accuracy and reliability in probing electronic structure properties.Employment of

IV. Conclusion
We presented ab initio study of the gap energy of SnO2.Our work provided various ways to apply periodic coupled cluster methods while there is a balance between cost and accuracy.We examined the efficiency of the twist averaging in reducing the finite-size errors through different functionals.
We showed that the gap obtained by the EOM-CCSD is overestimated in the thermodynamic extrapolation, indicating a large finite-size effect.We reported twist averaged EOM-CCSD band (∞) approach.Our found fundamental gap of 3.46 eV for the SnO2 is in close agreement with the two photon spectroscopy experiments.
of six oxygen atoms and each oxygen atom has three tins nearest neighbors at the corners of an almost equilateral triangle, see Figure 1.To obtain the bulk equilibrium structure, geometry is https://doi.org/10.26434/chemrxiv-2023-dks4g-v3ORCID: https://orcid.org/0000-0002-8259-8637Content not peer-reviewed by ChemRxiv.License: CC BY-NC-ND 4.0 relaxed with respect to the lattice a and internal parameter u using Vienna Ab-initio Simulation Package (VASP) 69 with PBE functional and kinetic energy cutoff of 400 eV and 16×16×16 Monkhorst-Pack 70 k-points mesh.Our relaxed lattice parameters of a = b = 4.827 and c = 3.246 A( , agree with the experimentally 67, 68 measured values of 4.737 and 3.186 A( , respectively, and generate experimental c/a ratio of 0.673.

Figure 1 .
Figure 1.Unit cell of SnO2.Red atoms are oxygen and grey atoms are tin.

Finite-size errors
due to finite simulation cells can be reduced by the Bloch theorem, where the Brillouin zone of the reciprocal lattice is spanned by enough boundary conditions or k-points to ensure the reliability of the probing quantities.Twist averaging is a well-established method to reduce the one-body (independent particle) finite-size errors (non-interacting kinetic, potential and Hartree energies, etc.) arising from incomplete k-points grid and facilitate the extrapolation to the TDL.Generally, integrating a periodic function over the Brillouin zone is performed at every kpoint.In the twist averaging, the k-points grid is offset from the origin by  & twist vectors, which is called twisted boundary conditions and then one takes average over all offsets, i.e., over the Bloch vectors  & in the first Brillouin zone of the simulation cell that is called twist averaging boundary condition (TABC).
, we present the direct gap values at Γ point within the DFT and EOM-CCSD formalisms using PySCF software package.For the DFT calculations we used the PBE0, B3LYP, and B3PW91 hybrid functionals.The band gap values are given with a n × n × n sampling of the https://doi.org/10.26434/chemrxiv-2023-dks4g-v3ORCID: https://orcid.org/0000-0002-8259-8637Content not peer-reviewed by ChemRxiv.License: CC BY-NC-ND 4.0 Brillouin zone, where n varies from 1 to 5. The utmost DFT calculations use a 3 × 3 × 3 sampling with the DZVP78 basis set for all atoms and 5 × 5 × 5 sampling with the DZVP basis set for the O atoms and the SBKJC 79 basis set for the Sn atoms.The reason lies in the fact that the cost of CCSD calculations is very high, so we opt for the SBKJC basis set for tin atoms to afford very demanding calculations.For the sake of comparison, we use the same basis set for the DFT calculations as well.Therefore, there is a chance of testing the reliability of the SBKJC basis set in comparison to the DZVP basis set in our estimations.We continue with this setting to obtain the EOM-CCSD gap.
, the TABC values are based on the 5 × 5 × 5 Monkhorst-Pack mesh and the Baldereschi point.It is seen that TABC results are in very close agreement with the TDL gap values.For instance, for a lower k-points, TABC calculations using 4 × 4 × 4 sampling and the Baldereschi point yields the same gap values of 4 × 4 × 4 k-points, without employing any twists.

Figure 2 .
Figure 2. Band gap extrapolation for the SnO2.The dashed lines and solid lines represent leastsquares extrapolation to the thermodynamic limit for each functional.The blue diamond show the CC gap values and the blue dashed line denotes least-squares extrapolation with the form  > +#/! .

(
∞) method based on the G J)K (∞), given in the Figure3, indicates that the PBE is not an appropriate functional in the DFT finite-size correction in the equation 2. Overall, our results indicate that twist averaged EOM-CCSD band gap energy augmented with KZK technique using PBE0 for the DFT finite-size correction is in close agreement with the gap energy based on two photon spectroscopies shown in TableI.As noted in the introduction, the lowest exciton transition of 3.56 eV measured in the two-photon absorption experiments 10 by adding the exciton binding energy of about 30 meV 9, 15, 80 yields a fundamental gap of 3.59 eV.Our twist averaged EOM-CCSD band gap energy augmented with KZK technique using PBE0 for the DFT finite-size correction, 3.46 eV, is in close agreement with this value.https://doi.org/10.26434/chemrxiv-2023-dks4g-v3ORCID: https://orcid.org/0000-0002-8259-8637Content not peer-reviewed by ChemRxiv.License: CC BY-NC-ND 4.0
gap energy augmented with KZK technique, G '()*+,-, **12 (∞) , using various functionals and samplings in the Brillouin zone to correct finite-size errors in the one-body and many-body context.The optimal sampling points are introduced to reduce the cost of the twist averaging.To https://doi.org/10.26434/chemrxiv-2023-dks4g-v3ORCID: https://orcid.org/0000-0002-8259-8637Content not peer-reviewed by ChemRxiv.License: CC BY-NC-ND 4.0 this end, Baldereschi point which is a special point in the Brillouin zone provided gap energies very close to the corresponding ones obtained by averaging over  & using Monkhorst IHIHI mesh in the Brillouin zone.Our results illustrated that we effectively addressed one-body and manybody errors in the gap energy estimation using G '()*+,-, **12 46,61,74,75, are not addressed by the twist averaging.Although, employment of large supercells can reduce such errors, it is impractical due to the cost and slow convergence.The KZK (Kwee-Zhang-Krakauer 61 ) approach is among the other techniques to reduce or cancel such errors; 76 this method includes finite-size corrections in the DFT calculations with finite-size functionals.We utilize this scheme to exploit corrections https://doi.org/10.26434/chemrxiv-2023-dks4g-v3ORCID: https://orcid.org/0000-0002-8259-8637Content not peer-reviewed by ChemRxiv.License: CC BY-NC-ND 4.0

TABLE IV
method suggests reliable gap energies but still the effect of the applied DFT functional in the G 29' (∞) term of the equation 2 is noticeable.