{'fraction_non_alphanumeric': 0.03253078478220915, 'fraction_numerical': 0.01180849108619739, 'mean_word_length': 4.538167938931298, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
First-principles LCPAO Approach for Insulators under Finite Electric Field

We propose a linear-combination-of-pseudo-atomic-orbitals scheme for a finite electric field method based on the modern theory of polarization. We derive the matrix elements of the effective potential for the field and the corresponding terms of the forces on atoms. In addition, we successfully evaluated the dielectric constants and Born effective charges of typical semiconducting and insulating materials. Our formalism will aid in the study of materials under electric fields.


Introduction
The effects of electric fields on insulators are of significant research interest, and have been investigated experimentally and theoretically for several applications. In theoretical works based on quantum mechanics, dielectric properties including responses for electric fields have been investigated intensively since the past fifty years [1][2][3]. In the modern theory of polarization (MTP) [4], first-principles calculations can be used on insulating solids to determine the electric polarization as a well-defined bulk quantity with Berry phases. In addition, polarization and Berry phases can be leveraged to probe the electronic structure of solids under electric fields [5].
MTP addresses the difficulty in describing electric polarization in bulk insulators with Berry phases obtained from overlap matrices between the periodic parts of Bloch orbitals at a k-point and the next k-point, instead of explicit evaluation of the position operator. The difficulty arises from the periodic boundary condition applied in electronic structure calculations based on the density functional theory (DFT), and it cannot be determined how electrons cross the cell boundary. This is because the existence of electrons is described by continuous electronic densities. Considering electric fields in the DFT framework, the boundary condition leads to another difficulty in that the periodicity can yield unphysical fields. Additional potentials for finite electric fields include terms arising from the product of the electric field and position operator E · r, but it is not acceptable when there are electrons at the boundary (i.e. saw-tooth potential in bulk systems). The position can be substituted with the polarization P to consider the potential for electric fields in bulk insulators with electric-field-polarized Wannier functions [6], because MTP shows that polarization is independent of the boundary, except for arbitrary modules of polarization. The first DFT calculations using E · P were reported in 1998, but large localization regions were required to get converged field-polarized Wannier functions [7]. In 2002, practical DFT schemes for E · P were proposed [5], and they were implemented within the plane wave method [5] or Car-Parrinello method [8]. The electronic structures of III-V or II-VI semiconductors were successfully evaluated under finite electric fields, and dielectric properties such as dielectric constants and the Born effective charges were also calculated. Such a scheme was also implemented in a projector augmented wave framework [9].
For realistic conditions, large scale calculations are generally required, and localized orbitals such as the linear combination of atomic orbitals (LCAO) method are appropriate in terms of scalability. Interesting models to consider the electric field effect on vacancies, impurities, and interfacial systems may require hundreds or thousands of atoms. Moreover, MTP for the LCAO scheme demands a conversion from Bloch orbitals to atomic orbitals [10]. Generally, atomic orbitals do not comprise the complete basis set, and therefore, Pulay's correction is required. These requirements make implementation complicated, which is a difficulty in the LCAO scheme. In the timedependent schrödinger equation, a time-evolving tight-binding scheme for finite electric fields was proposed based on MTP [11]. It was implemented in the time-dependent DFT framework of the SIESTA code, a linear combination of pseudo atomic orbitals (LCPAO) first-principle code [12,13] (the implementation details have not been reported). However, with regard to time-independent schemes, only a density matrix approach in the LCAO scheme for electric field effects has been reported thus far [14]. The LCAO implementation of effective potential and the forces on atoms for finite electric fields practically applicable to ordinal DFT is expected to expand the application range to consider even large scale or realistic models (e.g. vacancies in the diamond nitrogen-vacancy center) or complicated magnetic systems (e.g. spin spirals with the spin-orbit interaction (SOI)).
In this study, we developed a first LCAO implementation of the finite electric field scheme including the forces on atoms compatible with self-consistent field (SCF) loops in the DFT framework based on MTP. In addition, we evaluated electronic structures under uniform electric fields. Moreover, we calculated the dielectric properties of insulators and semiconductors: dielectric constants and Born effective charges of III-V and II-VI semiconductors and group IV insulating materials. The calculated dielectric properties agree with previously reported computational values. The LCAO implementation should serve as a powerful tool to investigate electric field effects on large scale systems.

LCPAO method and modern theory of electric polarization
The LCPAO method uses basis sets localized on atoms, and the basis functions are based on atomic orbitals [10]. The required number of basis functions is significantly smaller than that in the plane wave method, which enables solving eigenvalue problems by direct diagonalization. The electronic contribution to electric polarization is described via the Berry phase using Bloch orbitals. While Bloch orbitals can be expressed easily by plane wave basis sets, using localized basis sets of PAOs Bloch orbitals in the LCPAO method are given as where r is the position, N is the number of cells considered in calculations, "i" is the imaginary unit, k is a wave number, R n is the lattice vector for cell n, c (σk) iα,µ is an LCPAO coefficient connecting a PAO for orbital α belonging to atom i with a Bloch orbital of spin σ and state µ, τ i is the position of an atom i, and φ R iα is a PAO for orbital α belonging to atom i in a cell moved by R from the original cell, and it is expanded and can be solved by direct diagonalization to obtain the eigenvalues of energies ǫ (σk) . According to MTP [4], the electronic contribution to electric polarization P e in an insulating system with periodic cells is obtained from G a · P e , G b · P e , and G c · P e , which can be expressed by the Berry phase as where f is a weighting coefficient of the spin degeneracy per state, G a , G b and G c are the reciprocal lattice vectors for the cell vectors v a , v b , and v c , respectively, e is the elementary charge, Ω is the cell volume, N a , N b and N c are the numbers of k-points to discretize the first Brillouin zone along G a , is an overlap matrix of the periodic parts u σk µ of occupied Bloch orbitals between adjacent k-points, and k(I a , I b , The manner of expressing electric polarization by the Berry phase has been successful, since it provides evaluations based on well-defined quantities instead of direct evaluation of the position r. M (σ,Ia,I b ,Ic) a can be calculated by where we introduced T Na ·r φ Rn jβ , and µ and ν run over occupied states. Although M (σ,k(Ia,I b ,Ic),k(Ia+1,I b ,Ic)) itself is gaugedependent, the dependence vanishes when considering the sum over the first Brillouin zone, and the electric polarization P e can be evaluated. By introducing T [10] becomes a matrix form of Eq. 4, which is most suitable for computation in terms of efficiency and simplicity, especially for iterative calculations such as SCF calculations. This is because it suffices to prepare T (k) aiα,jβ once at the beginning of each SCF loop.

Effective potential for finite electric fields
We start with an expression of the electric enthalpy functional F , which is the total energy for an insulating system under a finite electric field: Here, E is the electric field, E KS is the ordinary counterpart of Kohn-Sham total energy described by field-polarized orbitals, and the polarization P comprises the ionic and electronic contributions: P = P I +P e . The solution that minimizes F [E] is desirable, such that diagonalization can be used to get solutions for insulators under electric fields in the same manner as in variational approaches to the ordinary Kohn-Sham equation. Here, considering an additional term of the effective potential V = −ΩE · P is required. In the plane wave method, Souza et al. suggested the conjugategradient method using a gradient δF /δ u (σk) µ | [5]. For the LCPAO method, we propose direct diagonalization by adding matrix elements of the effective potential into those of the ordinary Kohn-Sham Hamiltonian H 0 . The matrix representation of the effective potential A can be extracted from the relation Moreover, we expand Eq. 5 with the cell vectors v a , v b , and v c , and obtain where Then, we get Therefore, the matrix elements of A are given as Finally, in the case of a finite electric field, the Hamiltonian H (σk) = H (σk) 0 + A (σk) can be evaluated to obtain a solution through Eq. 2. In addition, this formalism is compatible with non-collinear spin density functionals, and is applicable to cases considering the SOI.

Forces under finite electric fields
For the forces on atoms induced by the effective potential, the ionic contribution is ΩE · ∂P I ∂τ i = Z i eE, where Z i is the core charge of atom i, and the electronic contribution in the case of the complete basis set is taken into account without additional terms, owing to the Hellmann-Feynman theorem. In the LCAO representation, the additional terms, i.e. Pulay forces, are necessary, so that the electronic contribution is iα,jβ = φ 0 iα |φ Rn jβ [10]. After SCF convergence, only the c † ∂T ∂τ i c term needs to be calculated explicitly, while evaluation of other terms in the electronic contribution is replaced with that of the ∂c † ∂τ i Hc + c † H ∂c ∂τ i This is realized using ε (σRn) determined after SCF convergence with the derivative of s (Rn) .
The norm-conserving pseudopotentials we used included the 2s and 2p electrons for C and O as valence electrons, while the 2p and 3s for Mg; 3s and 3p for Al, Si, P and S; 3d and 4s for Zn; 3d, 4s and 4p for Ga, Ge and As; 4s and 4p for Se; 4d and 5s for Cd; and 4d, 5s and 5p for In and Sb were also considered. We used a 32 × 32 × 32 regular k-point grid and real space grids corresponding to energy cutoffs larger than 300 Ry to obtain the converged results of the dielectric constants and Born effective charges. We also used the following standard and precise PAO basis sets ofXr-sn s pn p dn d fn f , where X, r and n s , n p , n d and n f denote the element and cutoff radius in the unit of Bohr and numbers of s-, p-, d-and f-orbital sets, respectively: for the We used computational models where the unit cell is the primitive cell of face centered cubic or zinc-blende structures with two atoms for computation. We applied an electric field of 0.1 GV/m along the a-axis of the conventional cubic cell including eight atoms. We optimized the lattice constants with the precise PAOs under the LDA, and w neglected the changes due to an external electric field. The optimized lattice constants were a = 2.4989 Å for C, a = 3.8219 Å for Si, a = 3.8471 Å for AlP, a = 3.9845 Å for AlAs, a = 3.8367 Å for GaP, a = 3.9998 Å for GaAs, a = 4.3184 Å for AlSb, a = 4.1413 Å for InP, a = 3.7529 Å for ZnS, a = 3.9356 Å for ZnSe, a = 4.2526 Å for ZnTe, a = 4.2446 Å for CdSe, a = 4.5433 Å for CdTe, and a = 2.9410 Å for MgO. The dielectric constant tensor ε was evaluated from where ǫ 0 and δ αβ are the electric constant and Kronecker delta, respectively. When the electric polarization induced by the electric field ∆P is given with and without atomic relaxation, ε αβ is equal to the static and electronic dielectric constants ε sαβ , ε ∞αβ , respectively, and the following can be derived where ∆P α = P α (E = E βêβ ) − P α (E = 0). In addition, through the forces F on atoms without the atomic relaxation, the Born effective charge tensor Z * of an atom was also evaluated from We then evaluated ε ∞xx , ε sxx and Z * xx . Here, we drop the subscripts and denote them as ε ∞ , ε s and Z * . The finite electric field was considered after SCF convergence under the zero field of each SCF loop, i.e., its effective potential was added after the first SCF convergence. The matrix elements of the effective potential was calculated from the Hamiltonian matrix in the previous SCF iteration and the overlap matrices. The overlap matrices were prepared at the beginning of each SCF loop. We noted that the forces converged when the magnitudes were less than 1 × 10 −5 Hartree/Bohr. This is because a Z * value of 1 yields a force of approximately 2 × 10 −4 Hartree/Bohr under an electric field of 0.1 GV/m, and highly precise evaluation of the forces is essential to describe the electronic structures under electric fields with ε s and Z * . In the LCPAO method, the so-called egg box effect leads to numerical errors arising from the real space grid, and it is critical in evaluating the forces, dipole moment, and stress tensor. The egg box effect can be resolved by grid cell sampling and iterative evaluation of such quantities with grid shifting and frozen density matrix for the LCPAO [10]. The density gradient in the GGA is sensitive to the egg box effect, especially since it is estimated from the electronic density on the real space grid. Therefore, we applied grid cell sampling to evaluate the forces precisely by introducing a fine real space grid of at least 100 × 100 × 100 after SCF convergence in the GGA cases. Another workaround is to consider a stronger electric field to enhance contributions from the electric field to the forces, but it sometimes violates the limitation of the electric field strength due to the Zener breakdown [5], and we did not treat the workaround in this work. To evaluate ε s , the atomic relaxation needs to be considered. However, the computational cost of the grid cell sampling with such fine grid was significantly high, because at least eight evaluations are required (e.g. 2 × 2 × 2 grid shifts). Therefore, we used the "one-shot" grid cell sampling, i.e., the forces based on the fine grid were evaluated from the density matrix obtained from the rough grid without the grid shifts.  [20] were obtained from density functional perturbation theory with GGA, while values of Ref. [5] and Ref. [9] were calculated from DFT with LDA and GGA, respectively, and the finite electric field methods based on MTP. Our

Results and discussion
To validate our formalism, we calculated the electronic dielectric constants ε ∞ of III-V, II-VI semiconductors and group IV insulators. First, we calculated ε ∞ with LDA to compare our values to those from the plane wave formalism based on the MTP by Souza et al. [5]. Our calculated values are shown in Table 1. Our values of the AlP and AlAs based on the precise PAO basis set coincide with their values, while the ones of the GaP and GaAs have more than 10% and 20% errors. Although our LCPAO formalism is equivalent to Souza et al.'s formalism, this difference might arise from different conditions such as pseudopotentials. Indeed, our measured values are in good agreement with recently reported values under GGA based on the density functional perturbation theory. Next, we calculated ε ∞ under the GGA, and the values were improved (i.e. closer to experimental values) for all the systems compared to LDA. We confirmed that a moderately large number of k-points are necessary to obtain the converged values as reported by Souza et al. [5]. We also investigated the PAO quality dependence, and the values for the standard basis sets are smaller than those for the precise ones, except the InP. The standard basis sets yield PAOs around atom centers in a more narrow domain than the precise basis sets. A more localized description of electronic states through the standard basis sets can constrain the electronic distribution more strongly, and suppress ε ∞ . These results demonstrated that our LCPAO implementation of the effective potential for the finite electric field was accurate.
To investigate the influence of the SOI, we focused on ZnTe and CdTe including tellurium, which is the heaviest element in our systems. When considering the SOI with the precise basis set under GGA, the ε ∞ of ZnTe and CdTe was 8.37 and 8.10, and increased by 0.13 and 0.21, respectively. The strong SOI can drastically change the electronic structure so that these increases occur. In fact, the band gap of ZnTe and CdTe changed from 1.73 eV and 1.04 eV to 1.60 eV and 0.77 eV after including the SOI, respectively (see also Appendix A for the relationship between the band gap and ε ∞ ). The ε ∞ of AlSb also increased by 0.11, and it became 11.01. Our results illustrated that for systems including heavy-elements, the SOI effect on ε ∞ should be considered.
As shown in Table 2, we also calculated the Born effective charges Z * of cations except for the group IV systems to validate our formalism of the forces on atoms under electric fields. Our calculated results were robust with respect to exchange correlation functionals, and each Z * under the LDA is approximately equal to that under GGA. Our values were in good agreement with not only the theoretical values, but also experimental values in previous studies. The mean absolute error of Z * for our values and the experimental values with respect to the systems was~6%. Therefore, we concluded that our formalism of forces under finite electric fields was also correct. In including the SOI, however, only negligible changes appeared. Although a fine real space grid allowed us to get converged values, we confirmed that even a rough real space grid could provide good values, which suggests that the computational cost can be reduced in calculating Z * .
Then, we evaluated the atomic relaxation under finite electric fields through the static dielectric constants ε s . As shown in Table 3, the values based on LDA and GGA were successfully reproduced, except for GaAs. The LDA value of GaAs was overestimated compared to the results of previous studies, because our calculated band gap of 0.25 eV was underestimated compared to the experimental value of 1.579 eV [31], which made more  [22][23][24] were obtained from DFT with LDA and the finite difference method based on the MTP, while those of Refs. [5,8] and Ref. [9] were calculated from DFT with LDA and GGA, respectively, and the finite electric field methods based on the MTP. Our calculated values under two types of conditions for the exchange correlation functional and the quality of PAOs are shown. The value of Ref. [24] is that of the anion.  delocalized pictures of electrons to give the overestimation. The GGA values were closer to the experimental ones than those of LDA, and the AlP, AlAs, and GaP values had errors within 7%, while the errors for GaAs and MgO were within~30%. Although we reproduced the LDA value of MgO reported in the previous theoretical study [8], there was a mismatch of the MgO values between the experimental and theoretical results. With an experimental lattice constant of 4.212Å [32] for its conventional cell, the ε s of MgO became 9.56, i.e., it approached the experimental ε ∞ of 9.8 [21], although our lattice constant was 1% smaller compared to the experimental value. In the atomic relaxation, the precision of the forces on atoms is of importance, and even minor errors in the magnitude of 10 −4 Hartree/Bohr can reflect large errors in ε s . As a result, it is suggested that not only a fine k-space grid, but also a fine real space grid is required to describe the movements of atoms in solids under a finite electric field. Although the computational cost of the ordinary grid cell sampling with a fine real space grid was significantly high to continue relaxation, we confirmed that the "one shot" grid cell sampling was efficient and effective in avoiding the numerical errors induced by the real space grid or density gradient. Finally, it was proved that our formalism for the forces was correct, and that it can be applied to both atomic relaxation and molecular dynamics under finite electric fields.

Conclusion
We developed a first-principle LCPAO scheme of a finite electric field method based on MTP to yield self-consistent solutions, including the forces on atoms in the DFT framework. Moreover, we investigated the PAO dependence and exchange correlation dependence by performing systematic calculations. Our implementation successfully reproduced the electronic and static dielectric constants and Born effective charges of III-V and II-VI semiconductors and group IV insulating materials. We confirmed that the implementation for the SOI is necessary for systems including heavy elements. Our LCPAO implementation is expected to be advantageous over conventional formalisms such as the plane wave method in terms of efficiency and extensibility, especially for large scale calculations modeling realistic systems through the O(N) method, and the atomic or orbital decomposition of the effective potential. Our formalism will expand the problem size such that complicated systems such as devices can be considered, and the molecular dynamics of complicated systems such as batteries under a finite electric field can be realized with a lower computational cost.