Reduced Scaling of Optimal Regional Orbital Localization via Sequential Exhaustion of the Single-Particle Space

Wannier functions have become a powerful tool in the electronic structure calculations of extended systems. The generalized Pipek-Mezey Wannier functions exhibit appealing characteristics (e.g., reaching an optimal localization and the separation of the σ–π orbitals) compared with other schemes. However, when applied to giant nanoscale systems, the orbital localization suffers from a large computational cost overhead if one is interested in localized states in a small fragment of the system. Herein, we present a swift, efficient, and robust approach for obtaining regionally localized orbitals of a subsystem within the generalized Pipek-Mezey scheme. The proposed algorithm introduces a reduced work space and sequentially exhausts the entire orbital space until the convergence of the localization functional. It tackles systems with ∼10000 electrons within 0.5 h with no loss in localization quality compared to the traditional approach. Regionally localized orbitals with a higher extent of localization are obtained via judiciously extending the subsystem’s size. Exemplifying on large bulk and a 4 nm wide slab of diamond with an NV– center, we demonstrate the methodology and discuss how the choice of the localization region affects the excitation energy of the defect. Furthermore, we show how the sequential algorithm is easily extended to stochastic methodologies that do not provide individual single-particle eigenstates. It is thus a promising tool to obtain regionally localized states for solving the electronic structure problems of a subsystem embedded in giant condensed systems.


Downfolded effective Hamiltonian
In large systems with a certain anisotropy (defects in semiconductors, molecules in solvent environments) all physical phenomena can be attributed to a small active space embedded in a host environment. Thus, it is common to map the problem onto the effective Hamiltonian, defined within an active space.
whereĉ † iσ andĉ i,σ are creation and annihilation operators in site i with spin σ andn † iσ is a particle number operator. The ε i , t ij are the on-site and hopping energies.
We extract the Hamiltonian parameters ε, t U i and V ij from the first-principles calcula-tions employing large supercells. To compute the onsite and hopping parameters we calculate the integral containing kinetic and ionic potential terms: The U i represents Coulomb on-site interactions of electrons with a different spin, while V ij is the Coulomb inter-site interaction, which we compute as: where, the V (r, r ) is the bare Coulomb interaction.
Excited states of the NV − center Table S1 shows the excited states of the NV − center computed in the basis of the Wannier functions that were obtained with different energy windows. The full space energy window is ∼ 24 eV below the Fermi energy. One can see that even 20 eV window results in an extremely underestimated result, while for 10 eV window the order of states is reversed. As a measure of the localization we report the value of the objective functional P (see main text). The P is set to 100% for case where the full space is used in the energy window. Preparation of stochastic basis using deterministic eigenstates The stochastic basis representing the complement (rest) space in our sF-PMWF calculations is prepared in a three-step manner. First, a random vector is constructed in the full space where m denotes the m th iteration in the outer-loop and |φ j is the eigenstate in the full space. The set of random coefficients {α m ij } are associated with the outer-loop step m, i.e., a different m corresponds to a different set of coefficients.
The second step is to perform Gram-Schmidt orthogonalization such that the stochastic basis is orthogonal to the core space where |ψ c k represents the state in the core space. The stochastic basis is then made mutually The last step is to normalize the stochastic basis such that and After these three steps, the construction of stochastic basis for the m th step is completed and it is ready to enter the work space.    Figure S3: The log of the time per outer-loop iteration (t outer ) as a function of the log of the number of states in the work space (N w ). The scaling of t outer with N w is derived from the slope of the linear fitting. Figure S4: Convergence of the functional P with respect to the outer-loop step m for the NV − center of the 215-atom system. Blue curve: localization performed with deterministic basis in the rest space. Orange curve: localization performed with stochastic basis in the rest space. The (16,32) combination is employed in both calculations.      Figure S7: Total wall time of orbital localization on each system with respect to the number of occupied states N s . Blue bar: F-PMWF using the full orbital space. Orange Bar: sF-PMWF using the work space.   Figure S8: Number of total SA iteration steps in sF-PMWF calculation relative to the number of total SA iteration steps in the F-PMWF calcuation for the 215-atom system using different N r . The N c is fixed at 16. Figure S9: Number of total SA iteration steps in sF-PMWF calculation relative to the number of total SA iteration steps in the F-PMWF calculation for the 511-atom system using different N r . The N c is fixed at 16. Figure S10: The log of the normalized time per iteration plotted as a function of the log of number of occupied states N s for the four investigated systems. The black line and square points represent the normalized t SA obtained from the F-PMWF method using the full orbital space. The red line and circle points represent the normalized t outer obtained from the sF-PMWF method using the constructed work space. The time per iteration is normalized to the largest grid (2303-atom system). The scaling is derived from the slope of each fitting.       Figure S15: The 4 regionally localized "p"-like states around the NV − center of the 511-atom system. The left 4 states are obtained from F-PMWF and the right 4 are obtained from sF-PMWF. The yellow and blue colors represent the phases of the single-particle wavefunction. The isosurface value is set at 0.05 for all the plots. Figure S16: The 4 regionally localized "p"-like states around the NV − center of the 999-atom system. The left 4 states are obtained from F-PMWF and the right 4 are obtained from sF-PMWF. The yellow and blue colors represent the phases of the single-particle wavefunction. The isosurface value is set at 0.05 for all the plots.

Supplementary Tables and Figures
Figure S17: Electron density constructed from the 4 regionally localized states on an arbitrary carbon of the four investigated systems: (a) 215-atom system; (b) 511-atom system; (c) 999-atom system; (d) 2303-atom system. The isosurface value is set at 0.01 for all the plots. Figure S18: The 4 regionally localized "p"-like states around the NV − center of the 215-atom system using different sizes of the fragment or using all the atoms. The last column shows the electron density constructed from these 4 states in each calculation. The isosurface value is set at 0.02 for all the plots.    Table S16: The spatial overlap between the two sets of "p-like" Wannier functions obtained from the sF-PMWF (ψ s ) method and the F-PMWF (ψ) method for the 2303-atom system.