Second-principles method including electron and lattice degrees of freedom

We present a first-principles-based (second-principles) scheme that permits large-scale materials simulations including both atomic and electronic degrees of freedom on the same footing. The method is based on a predictive quantum-mechanical theory, e.g., Density Functional Theory, and its accuracy can be systematically improved at a very modest computational cost. Our approach is based on dividing the electron density of the system into a reference part - typically corresponding to the system's neutral, geometry-dependent ground state - and a deformation part - defined as the difference between the actual and reference densities. We then take advantage of the fact that the bulk part of the system's energy depends on the reference density alone; this part can be efficiently and accurately described by a force field, thus avoiding explicit consideration of the electrons. Then, the effects associated to the difference density can be treated perturbatively with good precision by working in a suitably chosen Wannier function basis. Further, the electronic model can be restricted to the bands of interest. All these features combined yield a very flexible and computationally very efficient scheme. Here we present the basic formulation of this approach, as well as a practical strategy to compute model parameters for realistic materials. We illustrate the accuracy and scope of the proposed method with two case studies, namely, the relative stability of various spin arrangements in NiO and the formation of a two-dimensional electron gas at the interface between band insulators LaAlO$_3$ and SrTiO$_3$. We conclude by discussing ways to overcome the limitations of the present approach (most notably, the assumption of a fixed bonding topology), as well as its many envisioned possibilities and future extensions.


I. INTRODUCTION
Over the past two decades first-principles methods, in particular those based on efficient schemes like density functional theory (DFT) [1][2][3][4][5], have become an indispensable tool in applied and fundamental studies of molecules, nanostructures, and solids. Modern DFT implementations make it possible to compute the energy and properties (vibrational, electronic, magnetic) of a compound from elementary information about its structure and composition. Hence, in DFT investigations the experimental input can usually be reduced to a minimum (the number of atoms of the different chemical species, and a first guess for the atomic positions and unit cell lattice vectors). Further, the behavior of hypothetical materials can be readily investigated, which turns the methods into the ultimate predictive tool for application, e.g., in materials design problems.
However, interpreting or predicting the results of experiments requires, in many cases, to go beyond the time and length scales that the most efficient DFT methods can reach today. This becomes a very stringent limitation when, as it frequently happens, the experiments are performed in conditions that are out of the comfort zone of DFT calculations, i.e., at ambient temperature, under applied time-dependent external fields, out of equilibrium, under the presence of (charged) defects, etc.
The development of efficient schemes to tackle such challenging situations, which are of critical importance in areas ranging from biophysics to condensed matter physics and materials science, constitutes a very active research field. Especially promising are QM/MM multiscale approaches in which different parts of the system are treated at different levels of description: The most computationally intensive methods [based on quantum mechanics (QM), as for example DFT itself] are applied to a region involving a relatively small number of atoms and electrons, while a large embedding region is treated in a less accurate molecular mechanics (MM) way (e.g., by using one of many available semiempirical schemes).
Today's multiscale implementations tend to rely on semiempirical methods-like tight-binding [6,7] and forcefield [8,9] schemes-that were first introduced decades ago. In some cases, such schemes are designed to retain DFTlike accuracy and flexibility as much as possible. One relevant example are the self-consistent-charge density-functional tight-binding (DFTB) techniques [10][11][12], and related approaches [13][14][15], which retain the electronic description and permit an essentially complete treatment of the compounds. Another relevant example are the effective Hamiltonians developed to describe ferroelectric phase transitions and other functional effects [16][17][18]; these are purely lattice models (i.e., without an explicit treatment of the electrons) based on a physically-motivated coarse-grained representation of the material, and have been shown to be very useful even in nontrivial situations involving chemical disorder [19] and magnetoelectric effects [20], among others. Such methods have demonstrated their ability to tackle many important problems (see, e.g., Refs. [15,[21][22][23] for the DFTB approach), and constitute very powerful tools. Nevertheless, they are limited when it comes to treating situations in which the key interactions involve minute energy differences (of the order of meV's per atom) and a great accuracy is needed, or where a complete atomistic description of the material is required.
Another aspect in which many approximate approaches fail is in the simultaneous treatment, at a similar level of accuracy and completeness, of electronic and lattice degrees of freedom. Most methods in the literature are strongly biased towards either the electronic [24][25][26] or the lattice [8,9,[16][17][18] properties. Further, the few schemes that attempt a realistic, simultaneous treatment of both types of variables usually involve very coarse-grained representations [27][28][29].
Here we introduce a new scheme to tackle the problem of simulating both atomic and electronic degrees of freedom on the same footing, with arbitrarily high accuracy, and at a modest computational cost. In essence, we will combine (i) an accurate model potential to describe the lattice-dynamical properties of the system with (ii) a tight-binding-like approach to describe the relevant electronic degrees of freedom and electron-phonon interactions. Our scheme will be limited to problems in which it is possible to identify an underlying lattice or bonding topology that is not broken during the course of the simulation. As we will show below, such a fixed-topology hypothesis permits drastic simplifications in the description of the system, yielding a computationally efficient scheme whose accuracy can be systematically improved to match that of a DFT calculation, if needed. Note that, while our assumption of an underlying lattice may seem very restrictive at first, in fact it is not. There are myriads of problems of great current interest-ranging from electronic and thermal transport phenomena to functional (dielectric, ferroelectric, piezoelectric, magnetoelectric) effects and most optical properties-that are perfectly compatible with it. Further, this restriction can be greatly alleviated by combining our potentials with DFT calculations in a multiscale scheme, a task for which our models are ideally suited.
At its core, our new scheme relies on the usage of a force field to treat interatomic interactions, capable of providing a very accurate description of the lattice-dynamical properties of the material of interest. In particular, the scheme recently introduced by some of us in Ref. [30] constitutes an excellent choice for our purposes, as it takes advantage of the aforementioned fixed-topology condition to yield physically transparent models whose ability to match DFT results can be systematically improved.
Then, a critical feature of our approach is to identify such a lattice-dynamical model with the description of the material in the Born-Oppenheimer surface, i.e., with the DFT solution of the neutral system in its electronic ground state. Since the force fields of Ref. [30] do not treat electrons explicitly, this identification implies that our models will not tackle the description of electronic bonding, as DFTB schemes do. In other words, we will not be concerned with modeling the interactions responsible for the cohesive energy of the material, or for the occurrence of a certain basic lattice topology and structural features. Within our scheme, all such properties are simply taken for granted, and constitute the starting point of our models.
Instead, our models focus on the description of electronic states that differ from the ground state. These are the truly relevant configurations for the analysis of excitations, transport, competing magnetic orders, etc. By focusing on them, and by adopting a description based on material-(and topology-)specific electronic wave functions, we can afford a very accurate treatment of the electronic part while keeping the models relatively simple and computationally light.
As we will see, while it bears similarities with DFTB schemes, the present approach is ultimately more closely related to Hubbard-like methods. Yet, at variance with the usual semiempirical Hubbard Hamiltonians, our models are firmly based on a higher-level first-principles theory, treating all lattice degrees of freedom, and the relevant electronic ones, with similarly high (full DFT at the asymptotic limit) accuracy. The term "second-principles," used in the title of this paper, is meant to emphasize such a solid first-principles foundation.

II. OVERVIEW OF THE METHOD
In this paper we introduce the formal framework of our new computational scheme to perform large-scale simulations. The approach is methodologically based on DFT and obtains all the necessary information to simulate a material from this technique, so we have named it second-principles densityfunctional-theory (SP-DFT). Hence, the paper contains two main sections that describe (i) the development of the theory our models are based on as a systematic approximation to DFT (Sec. III) and (ii) the approach to calculate the model parameters from DFT (Sec. IV).
Given the generality of the approach and the many magnitudes that need to be calculated, we devote this section to (i) enumerate the most important foundations on which the method is based, (ii) point to the corresponding sections of this paper where the reader can find a full description of the various approximations and the physical phenomena behind them, (iii) present a comparison with previous methodologies based on similar ideas, and (iv) describe how to improve on the approximations made in a systematic way to achieve full DFT accuracy as the asymptotic limit.
(1) As in DFTB approaches, the starting point of our method is the expansion of the DFT total energy (Sec. III B) with respect to charge density fluctuations around a reference electron density (Sec. III A). Our first approximation involves limiting the expansion to second order. We choose our expansion so that the zeroth order term is the total DFT energy for the reference electron density, which we represent using an accurate and efficient force field. The first and second order terms involve small corrections to this energy and require explicit electronic structure calculations. The key to the success of this approach is that, contrary to other electronic structure calculation methods, the correction (and the self-consistent equations, Sec. III H) only depends on the difference density, and not on the total density. This magnitude involves few electrons and does not require knowledge of all the electronic states of the system, and can thus be handled with relatively modest computational resources.
(2) The electronic wave functions are expanded in a minimal basis set made of localized Wannier functions (Sec. III C). The Hamiltonian matrix elements are represented in this basis set, and the resultant one-and two-electron integrals become the main parameters of the calculation (Sec. III C). While these parameters can, in principle, be computed from DFT simulations, our approach here is to fit them so that the model reproduces a training set of first-principles data. An extension of the expressions to treat magnetic systems is developed in Sec. III D, where we see that Hubbard-and Stoner-like parameters appear in a natural way.
(3) Although our basis set is localized in space, and therefore the relevant real space matrix elements are restricted to only relatively close neighbors, their values can be affected by the electrostatic interactions due to charges and electric dipoles located at distant regions of the material. Thus, we split our one-and two-electron integrals into near-field and far-field contributions. The latter are computed using the multipole expansion described in Sec. III E, while the former are obtained from DFT using the recipe provided in Sec. IV.
(4) We introduce an explicit dependence of the electronic Hamiltonian on atomic positions by expanding the oneelectron parameters as a function of distortions of a reference structure. We thus account for electron-lattice couplings (Sec. III F). At this level a new approximation is introduced to neglect the dependence of the two-electron integrals with the geometry.
After finishing with the foundations of the method and gathering together all the previous ingredients, we find an expression for the total energy of the system (Sec. III G). Its minimization with respect to the coefficients of the wave functions in a Wannier basis gives rise to a set of Kohn-Sham-like equations that need to be solved self-consistently (Sec. III H). After convergence, the forces on all the atoms and the stresses on the unit cell lattice vectors are computed (Sec. III I), allowing for structural-relaxation or moleculardynamics simulations over very large length and time scales. To end the section on the method itself, we give in Sec. III J a brief account of some practicalities concerning an efficient implementation of our approach.
In Sec. IV we propose a tentative approach for a systematic calculation of the model variables from first principles. Note that the development of a systematic-and automaticstrategy for the construction of models with predefined accuracy is a technically challenging task that remains for future work.
In Sec. V we provide the details about the SCALE-UP code package developed in the course of this work, stressing the computational efficiency of the simulations based on the new methodology. This is followed by Sec. VI where we compare and highlight the novelty of our approach to two closely related techniques, the effective-Hamiltonian and the DFTB methods.
Then, we describe in Sec. VII a couple of nontrivial applications that were selected to highlight the flexibility of the models, the great physical insight that they provide, and their ability to account for complex properties with DFT-like accuracy. The two chosen systems present interactions of very different origin: (i) the magnetic Mott-Hubbard insulator NiO and (ii) the two-dimensional electron gas (2DEG) that appears at the interface between band insulators LaAlO 3 and SrTiO 3 . Note that, while accuracy will be highlighted, in these initial applications we have focused on testing the ability of our scheme to tackle challenging situations from the physics standpoint, and to produce reasonable predictions beyond the information including in the training set. Finally we present our conclusions and a brief panoramic overview of future extensions of the method and possible fields of application in Sec. VIII.

A. Basic definitions
As customarily done in most first-principles schemes, we assume the Born-Oppenheimer approximation to separate the dynamics of nuclei and electrons. Hence, we consider the positions of the nuclei as fixed parameters of the electronic Hamiltonian. Our approach will give us access to the potential energy surface (PES), i.e., for each configuration of the nuclei, the total energy of the system will be computed.
Our goal is to describe the electrons in the system, and the relevant electronic interactions, in the simplest possible way. Hence, we will typically focus on valence and conduction states, and will thus work with a lattice of ionic cores comprised by the nuclei and the corresponding core electrons, which will not be modeled explicitly. Here we use indistinctively the terms atoms, ions, and nuclei to refer to such ionic cores.
Our method relies on the following key concepts: the reference atomic geometry (RAG henceforth) and the reference electronic density (RED in the following). As in the recent development of model potentials for lattice-dynamical studies [30], the first step towards the construction of our model is the choice of a RAG, that is, one particular configuration of the nuclei that we will use as a reference to describe any other configuration. In principle, no restrictions are imposed on the choice of RAG. However, it is usually convenient to employ the ground state structure or, alternatively, a suitably chosen high-symmetry configuration. Note that these choices correspond to extrema of the PES, so that the corresponding forces on the atoms and stresses on the cell are zero. Further, the higher the symmetry of the RAG, the fewer the coupling terms needed to describe the system and, in turn, the number of parameters to be determined from first principles.
To describe the atomic configuration of the system we shall adopt a notation similar to that of Ref. [30]. In what follows, all the magnitudes related with the atomic structure will be labeled using Greek subindices. For the sake of simplicity, we shall assume a periodic three dimensional infinite crystal, with the lattice cells denoted by uppercase letters ( , , . . .) and the atoms in the cell by lowercase letters (λ, δ, . . .). In this manner, the lattice vector of cell is ⃗ R , and the reference position of atom λ is ⃗ τ λ . In order to allow for a more compact notation, a cell/atom pair will sometimes be represented by a lowercase bold subindex, so that ⃗ R λ ↔ λ. Any possible crystal configuration can be specified by expressing the atomic positions, ⃗ r λ , as a distortion of the RAG, as where 1 is the identity matrix, ← → η is the homogeneous strain tensor, and ⃗ u λ is the absolute displacement of atom λ in cell with respect to the strained reference structure.
The second step is the definition of a RED, n 0 (⃗ r), for each possible atomic configuration. Our method relies on the fact that, in most cases, the self-consistent electron density, n(⃗ r), will be very close to the RED, so that changes in physical properties can be described by the small [with respect to n 0 (⃗ r)] deformation density, δn(⃗ r), defined as n(⃗ r) = n 0 (⃗ r) + δn(⃗ r). ( Several remarks are in order about Eq. (2). First, with n(⃗ r) we denote the electron density that integrates to the number of electrons (i.e., it is positive). It is trivially related with the charge density (in atomic units it just requires making it negative due to the sign of the electronic charge).
Second, this separation of the charge density into reference and deformation contributions is similar to what is commonly found in DFTB schemes, and even in first-principles methods [31]. However, this parallelism may be misleading. Indeed, it is important to note that we make no assumption on the form of the RED. In most cases, e.g., nonmagnetic insulators, it will be most sensible to identify the RED with the ground state density of the neutral system. Nevertheless, as will be illustrated in Sec. VI B for Mott insulator NiO, other choices are also possible and very convenient in some situations.
Third, our RED will typically be an actual solution of the electronic problem, as opposed to some approximate density, e.g., a sum of spherical atomic-like densities, possibly taken from the isolated-atom solution, as used in some DFTB schemes [15,32]. Fourth, the concepts of RAG and RED are completely independent: In our formalism, we define a RED for every atomic structure accessible by the system, and not only for the reference atomic geometry. Finally, let us remark, in advance to Sec. III J, that our method does not require an explicit calculation of n 0 (⃗ r) (or any other function in real space, for that matter), a feature that allows us to reduce the computational cost significantly.
In order to further clarify the concept of RED, let us discuss the application of our method to the relevant case of a doped semiconductor. As sketched in Fig. 1, our hypothetical semiconductor is made of two different types of atoms (represented by large green and small red balls, respectively) in a square planar geometry with a three-atom repeated cell. The RAG corresponds to the high-symmetry configuration in which the large atom is located at the center of the square, while the small atoms lie at the centers of the sides. In the neutral (undoped) case, a self-consistent DFT calculation of the RAG would yield an electronic configuration with all the valence bands occupied and all the conduction bands empty, as illustrated in Fig. 1(e). The associated electron density would be our RED, n 0 (⃗ r), represented by the green clouds in Fig. 1(b); the associated energy would be E (0) , using the notation that will be introduced in Sec. III B. Now, if we dope the neutral system by adding or removing electrons, there will be a response of the electronic cloud, which will tend to screen the field caused by the extra charge. The doping electron (respectively, hole) will occupy the states at the bottom of the conduction band (respectively, top of the valence band). The doping-induced charge redistribution can be viewed as resulting from an admixture of occupied and unoccupied states of the reference neutral configuration.
FIG. 1. Schematic cartoon that represents the key physical concepts for the development of the second-principles models: the reference atomic structure and the reference and deformation electron densities. Panels (a)-(c): the meaning of the balls (which represent the position of the atoms in a hypothetical semiconductor), and the green clouds (which represent charge densities) are explained in the main text. Panels (d)-(f): the horizontal lines represent the one-electron energy levels obtained at the corresponding atomic structures and for the reference electronic configuration (neutral ground state). Full green circles represent full occupation of a given state by electrons. Half filled orange/green circles indicate partial occupation of a particular level. The notations E (0) , E (1) , and E (2) for the energies are introduced in Sec. III B. The parameters γ , U , and, I are defined in Secs. III C and III D. Only the case of doping with electrons is sketched. Doping with holes would lead to an equivalent picture. The resulting state, described by the total charge density n(⃗ r), is sketched in Figs. 1(a) and 1(d). The difference between the total electronic density and the RED is the deformation density δn(⃗ r). Such a deformation density, which is the key quantity in our scheme, captures both the doping and the system's response to it, as sketched in Figs. 1(c) and 1(f).
Finally, let us further stress the independence between RAG and RED. Note that all three quantities n(⃗ r), n 0 (⃗ r), and δn(⃗ r) are in fact parametric functions of the atomic positions. This is illustrated in Fig. 2, which sketches a case in which one atom is displaced from the RAG.

B. Approximate expression for the energy
Let us consider an atomic geometry characterized by the homogeneous strain tensor ← → η and the individual atomic displacements {⃗ u λ } as described in Eq. (1). Our main objective is to find a functional form that permits an accurate approximation of the DFT total energy at a low computational cost. The DFT energy can be written as FIG. 2. Schematic cartoon emphasizing that the division of the electron density into reference and deformation parts is carried out for any geometrical configuration of the system, as defined by the strain ← → η and atomic displacements {⃗ u λ }. The distortion of the reference atomic geometry is illustrated by the off-centering of one atom (indicated with a black arrow). The atomic distortion results in a modified n 0 (⃗ r) [panel (b)], as well as in additional changes depicted in panel (c), where the green and orange clouds denote positive and negative variations in the electronic density. Symbols have the same meaning as in Fig. 1.
In this expression, the first term on the right-hand side includes the kinetic energy of a collection of noninteracting electrons as computed through the one-particle kinetic energy operator,t; this first term also includes the action of an external potential, v ext , which gathers contributions from the nuclei (or ionic cores) and, possibly, other external fields. The second term is the Coulomb electrostatic energy, which in the context of quantum mechanics of condensed matter systems is also referred to as the Hartree term. The third term, E xc [n], is the so-called exchange and correlation functional, which contains the correlation contribution to the kinetic energy in the interacting electron system, as well as any electron-electron interaction effect beyond the classic Coulomb repulsion. The last term, E nn , is the nucleus-nucleus electrostatic energy. Note that Eq. (3) is written in atomic units, which are used throughout the paper. (|e| = m e = = a B = 1, where |e| is the magnitude of the electron charge, m e is the electronic mass, and a B is the Bohr radius).
As already mentioned, we assume that the Born-Oppenheimer approximation applies, so that the positions of the nuclei can be considered as parameters of the Hamiltonian. We also assume periodic boundary conditions. (Finite systems can be trivially considered by, e.g., adopting a supercell approach [33].) Within periodic boundary conditions, the eigenfunctions of the one-particle Kohn-Sham equations, |ψ j ⃗ k ⟩, can be written as Bloch states characterized by the wave vector ⃗ k and the band index j , with the occupation of a state given by o j ⃗ k . Note that Eq. (3) is valid for any geometric structure of the system, and we implicitly assume that the total energy (E DFT ), the one-particle eigenstates (|ψ j ⃗ k ⟩), and all derived magnitudes (such as the electron densities n, n 0 , and δn) depend on the structural parameters ← → η and {⃗ u λ }.
The total energy of Eq. (3) is a functional of the density which, as described in Eq. (2), can be written as the sum of a reference part n 0 (⃗ r) and a deformation part δn(⃗ r).
where we have introduced functional derivatives of E xc . In principle, Eq. (4) is exact if we consider all the orders in the expansion. (Expansions like this one are frequently found in the formulations of the adiabatic density functional perturbation theory [34,35].) In practice, under the assumption of a small deformation density, the expansion can be cut at second order. As we shall show in Secs. III D and III H, this approximation includes as a particular case the full Hartree-Fock-theory; hence, we expect it to be accurate enough for our current purposes. Within the previous approximation, we can write the total energy as a sum of terms coming from the contributions of the deformation density at zeroth (reference density), first, and second orders. Formally we write where the individual terms have the following form [36]. For the zeroth-order term E (0) , we get The above equation corresponds, without approximation, to the exact DFT energy for the reference density n 0 . We can choose the RED so that E (0) is the dominant contribution to the total energy of the system, and here comes a key advantage of our approach: We can compute E (0) by employing a model potential that depends only on the atomic positions, where the electrons (assumed to remain on the Born-Oppenheimer surface) are integrated out. This represents a huge gain with respect to other schemes that, like the typical DFTB methods, require an explicit and accurate treatment of the electronic interactions yielding the RED as well as solving numerically for E (0) and n 0 for each atomic configuration considered in the simulation. The first-order term involves the one-electron excitations as captured by the deformation density, Here,ĥ 0 is the Kohn-Sham [1] one-electron Hamiltonian defined for the RED, where v H (n 0 ; ⃗ r) and v xc [n 0 ; ⃗ r] are, respectively, the reference and exchange-correlation, potentials. It is important to note that Eq. (7) is different from the one usually employed in DFTB methods (see, for example, Refs. [10] and [12]): while typical DFTB schemes include a plain sum of one-electron energies, here we deal with the difference between the value of this quantity for the actual system and for the reference one [see sketch in Fig. 1(f)]. Such a difference is a much smaller quantity, more amenable to accurate calculations. Finally, the two-electron contribution from the deformation density, E (2) , is given by where the screened electron-electron interaction operator, g(⃗ r,⃗ r ′ ), is Here, δ 2 E xc /δn(⃗ r)δn(⃗ r ′ ) captures the effective screening of the two-electron interactions due to exchange and correlation. The latter magnitude is particularly important in chemistry, as it is related to the hardness of a material [37]. In summary, in this section we have shown how a particular splitting of the total density, into reference and a deformation parts, allows us to expand the DFT energy around n 0 and as a function of δn. This expansion can be truncated at second order while keeping a high accuracy; nevertheless, this approach can be systematically improved by including higher-order terms in δn, in analogy to what is done, e.g., in the so-called DFTB3 method [21]. While the general idea is reminiscent of DFTB methods in the literature, [10][11][12] our scheme has two distinct advantages. On one hand, the zeroth-order term can be conveniently parametrized by means of a lattice model potential, so that it can be evaluated in a fast and accurate way without explicit consideration of the electrons. On the other hand, the first-order term is much smaller, and can be calculated more accurately, than in usual DFTB schemes, as it takes the form of a perturbative correction.

Choice of Wannier functions
Our formulation requires the computation of the matrix elements of the Kohn-Sham one-electron Hamiltonian, as defined for the RED, in terms of Bloch wave functions [Eq. (7)], as well as various integrals involving the deformation charge density [Eq. (11)]. To compute these terms, we will expand the Bloch waves on a basis of Wannier-like functions (WFs), |χ a ⟩, in the spirit of Ref. [38]. There are several reasons for our choice of a Wannier basis set over the atomic orbitals most commonly used in DFTB formulations [10][11][12]15].
First, the Wannier orbitals are naturally adapted to the specific material under investigation. In fact, they will be typically obtained from a full first-principles simulation of the band structure of the target material, which permits a more accurate parametrization of the system while retaining a minimal basis set.
Second, the Wannier functions are spatially localized, and several localization schemes are available in the literature [38][39][40][41][42][43][44]. The localization will be exploited in our secondprinciples method to restrict the real-space matrix elements to those involving relatively close neighbors, as will be explained in Sec. IV.
Third, the localized Wannier functions can be chosen to be orthogonal. Note that methods with nonorthogonal basis functions require the calculation of the overlap integrals that have a nontrivial behavior as a function of the geometry of the system. Moreover, the one-particle Kohn-Sham equations in matrix form become a generalized eigenvalue problem, whose solution requires a computationally demanding inversion of the overlap matrix. The use of orthogonal Wannier functions allows us to bypass these shortcomings.
Fourth, the Wannier functions enable a very flexible description of the electronic band structure, as they can be constructed to span the space corresponding to a specific set of bands [38,45]. Therefore, the electronic states can be efficiently split into: (i) an active set playing an important role in the properties under study; and (ii) a background set that will be integrated out from the explicit treatment. For instance, if the problem of interest involves the formation of low-energy electron-hole excitons, our active set would be comprised by the top-valence and bottom-conduction bands, and we would use the corresponding Wannier functions as a basis set.
Typically, we will start from a set of Bloch-like Hamiltonian eigenstates |ψ (0) n ⃗ k ⟩ that define a manifold of J bands associated to the RED. Then, following, e.g., the recipe of Ref. [42], we have where the Wannier function | ⃗ R A a⟩ is labeled by the cell A at which it is centered (associated to the lattice vector ⃗ R A ) and by a discrete index a. Note that we use Latin subindices to label all physical quantities related with the electrons; to alleviate the notation, we group in the bold symbol a both the cell and the discrete index, so that a ↔ ⃗ R A a. In Eq. (13), V is the volume of the primitive unit cell, the integral is carried out over the whole Brillouin zone (BZ), the index m runs over all the J bands of the manifold, and the T ( ⃗ k) matrices represent unitary transformations among the J Bloch orbitals at a given wave vector. Figure 3 shows three paradigmatic examples: a nonmagnetic insulator [bulk SrTiO 3 , Fig. 3(a)], a nonmagnetic metal [bulk Cu, Fig. 3(b)], and an antiferromagnetic insulator [bulk NiO, Fig. 3(c)]. In the first case, the valence bands are well separated in energy from other bands; further, they have a well-defined character strongly reminiscent of the corresponding atomic orbitals. More precisely, three isolated manifolds corresponding to the occupied valence bands-with dominant O-2s, Sr-4p, and O-2p character, respectively-are clearly visible; these bands are centered around 17 eV, 15 eV, and 3 eV below the valence-band top, respectively. The Bloch eigenstates for these bands can be directly used to compute the corresponding localized Wannier functions following the scheme of Ref. [42] or similar ones. In contrast, the bottom conduction bands of SrTiO 3 have a dominant Ti-t 2g character, but overlap in energy with higher-lying (Ti-e g ) conduction bands. The situation is even more complicated in the cases of Figs. 3(b) and 3(c), where the critical bands-i.e., those around the Fermi energy in the case of Cu, and those comprising the Ni-3d manifold in the case of NiO-are strongly entangled with other states. In such cases, we may need to use a disentanglement method-like, e.g., the one proposed in Ref. [38]-to identify a minimal active manifold.
Note that the inverse transformation from Wannier to Bloch functions reads where the connection between the c (0) ja ⃗ k coefficients and the transformation matrices in Eq. (13) is given in Ref. [36]. The Wannier functions corresponding to the RED [ Fig. 1(b)] form a complete basis of the Hilbert space. Hence, we can use them to represent any perturbed electronic configuration of the system [e.g., the one sketched in Fig. 1

(a)] as
where the sum can be extended to as many bands as needed to accurately describe the phenomenon of interest. (As in the example of Fig. 1, this might be the addition of an electron and the associated screening.) Finally, note that the Wannier function basis is implicitly dependent on the structural parameters ← → η and {⃗ u λ }, and it should be recomputed for every new RED corresponding to varying atomic positions. Ultimately, our models will capture all such effects implicitly in the electron-lattice coupling terms, whose calculation is described in Sec. III F. Also, henceforth we will assume that each and every one of the WFs in our basis can be unambiguously associated with a particular atom at (around) which it is centered. Further, we will use the notation a ∈ λ to refer to all the WFs associated to atom λ in cell , an identification that will be necessary when discussing our treatment of electrostatic couplings in Sec. III E.

Equations in a Wannier basis
Using Eq. (15), we can write the electron density n(⃗ r) in terms of the Wannier functions, We can assume we will work with real Wannier functions [45] and therefore drop the complex conjugates in our equations. In Eq. (16) we have introduced a reduced density matrix, which, following the nomenclature of Ref. [46], will be referred to as the occupation matrix for the WFs. This occupation matrix has the usual properties, including periodicity when the Wannier functions are displaced by the same lattice vector in real space. Equation (16) can similarly be applied to the RED, where the calculation of the occupation matrix is performed with the coefficients of the Bloch functions that define the reference electronic density, c (0) j α ⃗ k , as in Eq. (14). In order to quantify the difference between the two densities defined in Eqs. (16) and (18), we introduce a deformation occupation matrix, which will be the central magnitude in our calculations. Now the deformation density can be written as Using these definitions, we can rewrite the E (1) and E (2) energy terms. Introducing Eq. (19) into Eqs. (7) and (11) we get, after some algebra [36] and respectively, where γ ab and U aba ′ b ′ are the primary parameters that define our electronic model. These parameters can be obtained, respectively, from the integrals of the one-and twoelectron operators computed in DFT simulations, as and Alternatively, they can be fitted so that the model reproduces a training set of first-principles data.

D. Magnetic systems
The above expressions are valid for systems without spin polarization. The procedure to construct the energy for magnetic cases is very similar, but there are subtleties pertaining to the choice of RED.
In principle, one could use a RED corresponding to a particular realization of the spin order, e.g., the antiferromagnetic ground state for a typical magnetic insulator, or the ferromagnetic ground state for a typical magnetic metal. However, such a choice is likely to result in a less accurate description of other spin arrangements, which would hamper the application of the model to investigate certain phenomena (e.g., a spin-ordering transition).
Alternatively, one might adopt a nonmagnetic RED around which to construct the model. Such a RED might correspond to an actual computable state: For example, it could be obtained from a nonmagnetic DFT simulation in which a perfect pairing of spin-up and spin-down electrons is imposed. Further, as we will see below for the case of NiO, in some cases it is possible and convenient to consider a virtual RED whose character can be inspected a posteriori. This latter option follows the spirit of the usual approach to the construction of spin-phonon effective Hamiltonians [47], where the parameters defining the reference state cannot be computed directly from DFT, but are effectively fitted by requesting the model to reproduce the properties of specific spin arrangements.
In the following we assume a nonmagnetic RED, and present an otherwise general formulation. The E (0) and E (1) terms thus describe the lattice and one-electron energetics corresponding to the nonmagnetic RED, and do not capture any effect related with the spin polarization. In contrast, the screened electron-electron interaction operator [Eq. (12)] is spin dependent and equal to where s and s ′ are spin indices that can take "up" or "down" values which we denote, respectively, by ↑ and ↓ symbols. This distinction in the screened electron-electron operator leads us to introduce two kinds of U parameters, and which describe, respectively, the interactions between electrons with parallel (U par ) and antiparallel (U anti ) spins. As a consequence, E (2) in spin-polarized systems is where and D s ab is the deformation occupation matrix for the s spin channel, defined for the up and down spins as (30) and respectively. Replacing Eqs. (29)-(31) into Eq. (28), For physical clarity, and to establish the link of Eqs. (26) and (27) with Eq. (24), it is convenient to write U par and U anti in terms of Hubbard-(U ) and Stoner-(I )-like parameters: It is also convenient to introduce and so that Eq. (32) can be rewritten as: Note that the value of U in Eqs. (33) and (34) is consistent with the one in Eq. (24) if we consider a non-spin-polarized density In addition, note that the newly introduced constant I aba ′ b ′ only plays a role in spin-polarized systems and is necessarily responsible for magnetism.

Connection with other schemes
The two-electron interaction constants-U and I defined in Eqs. (33) and (34), respectively-are formally similar to the four-index integrals typically found in Hartree-Fock theory [4] and can be chosen to completely match this approach. However, one should note that the electron-electron interaction in our Hubbard-like and Stoner-like constants is not the bare one, but is screened by the exchange-correlation potential associated to the reference density, n 0 [see Eq. (25)]. This fact brings our formulation closer to the so-called DFT+U [48,49] and GW [50,51] methods.
Looking in more detail at our expressions for U and I , we find that they are very similar to those of U par [Eq. (26)] and U anti [Eq. (27)], except that the operator involved in the double integral is, respectively, and Thus, we see that U contains the classical Hartree interactions, screened by exchange and correlation. Moreover, from Eq. (39) we see that U , as used here, is related with the deformation occupation matrix D U , that captures the total change of the electron density (i.e., the sum of the deformation occupation matrix for both components of spins). Therefore, it is consistent with the usual definition U = d 2 E/dn 2 , i.e., it quantifies the energy needed to add or remove electrons. On the other hand, Stoner's [52][53][54] I only includes terms with quantum origin. In particular g I provides the difference in interaction between electrons with parallel and antiparallel spins.

One-electron parameters
The matrix element γ ab [Eq. (23)] gathers Coulomb interactions associated to the electrostatic potential created by both electrons and nuclei, which acts on the WFs χ a and χ b . Note that these are the only long-ranged interactions in the system, since all other contributions (kinetic, exchange-correlation, external applied fields) can be considered local or semilocal. In the following we discuss the detailed form of this electrostatic part of γ ab , which we denote γ elec ab . Let us first consider the part of γ elec ab associated to the electrostatic potential created by the electrons, γ elec,e ab . We have The expression of the reference electron density in terms of the occupation of Wanniers, o (0) c , and squares of Wannier functions in the reference state will be described in more detail in Sec. III J. Following the criteria of Ref. [55], the one-electron matrix elements related with the Coulomb electron-electron interaction can be split into two categories (see Fig. 4): (i) the near-field regime, where the two WFs (a and b) significantly FIG. 4. Schematic representation of the near-and far-field interactions. The shape of the orbitals (represented here by two t 2g -like WFs labeled a and b) is important in the determination of the short-range part of the γ and U interactions. In addition, the diagonal terms like γ aa and U aabb also include far-field effects due to charges and dipoles at distant regions of the material (see WF c in the figure). As regards the far-field interactions, the precise shape of the charge distributions generating the potential is not critical (illustrated by the diffuse orbital at point c), and can be approximated by a multipole expansion.
overlap with the third WF (c) that creates the electrostatic potential, and (ii) the far-field regime, where this overlap is negligible.
In the far-field regime, the electrostatic potential outside the region where a source charge χ c is located can be expressed as a multipole expansion (see Chapter 4 of Ref. [56]). More precisely, we can write the far-field (FF) potential created by the charge distribution given by which applies to ⃗ r points for which χ c (⃗ r) ≈ 0. Now, let λ label the atom-located at the RAG reference position ⃗ Then, assuming that |⃗ r ′′ | ≪ |⃗ r − ⃗ τ λ | and using the superscript T to indicate the transpose operation, as necessary to compute inner dot products, we get Now, substituting Eq.
The coefficient of the first term is the total charge (i.e., the monopole), and it is given by The coefficient of the second term is the electric dipole moment associated to χ c , which amounts to where ⃗ r c represents the centroid of χ c . Quadrupole and higherorder moments follow in the expansion, but here we assume they can be neglected. Finally, the full FF potential created by the electrons at point ⃗ r is simply given by where the prime indicates that we sum only over WF's such that χ c (⃗ r) ≈ 0.
Let us now consider the part of γ elec ab associated to the potential created by the nuclei, which we call γ elec,n ab . In analogy with the electronic case, we write the FF electrostatic potential created by the nuclei at point ⃗ r as where the primed sums run only over atoms λ whose associated WFs a ∈ λ satisfy χ a (⃗ r) ≈ 0.
Then, adding all far-field contributions to γ elec ab , and assigning each WF to its associated nucleus, we get where is the charge of ion λ, while ⃗ p λ is the local dipole associated to that very ion. Note that we add together the contributions from electrons and nuclei, which allows us to talk about ions in a strict sense [57]. We can further approximate this local dipole using the Born charge tensor ← → Z * λ , to obtain In order to get the final expression for the FF potential, we note that the electrostatic interactions described above do not take place in vacuum, but in the material at its reference electronic density. Thus, we need to take into account that the RED will react to screen such interactions, and that such a screening can be modelled by the high-frequency dielectric tensor of the material at its RED. Thus, the far-field potential at the center of WF where ⃗ e λa is a unitary vector parallel to ⃗ τ λ − ⃗ r a , ← → ϵ ∞ is the highfrequency dielectic tensor, and the primed sums are restricted in the usual way.
We can now divide γ ab in long-range (lr) and short-range (sr) contributions (see Fig. 4). Considering that χ a and χ b are strongly localized and orthogonal to each other, we define γ lr ab as Then, we effectively define the short-range part of γ ab as Note that the short-range interactions defined in this way include electrostatic effects as well as others associated to chemical bonding, orbital hybridization, etc. These interactions do not have a simple analytic form; hence, in order to construct our models, they will generally be fitted to reproduce DFT results.
It is important to note that the above derivation, and decomposition in long-and short-range parts, is exact and does not involve any approximation, except for: (i) the truncation of the multipole expansion and (ii) the analytic form introduced for the long-range electrostatic interactions, which strictly speaking only applies to homogeneous materials with a band gap [58].
Finally, note that the γ ab couplings can be expected to be short in range, as they involve WFs χ a and χ b that are strongly localized in space and decay exponentially as we move away from their centers. Hence, the γ matrix can be expected to be sparse, which will result in more efficient calculations. It is important to note that this short-ranged character of the γ ab couplings is expected despite the fact that the interactions contributing to γ lr ab are electrostatic and thus long ranged.

Two-electron integrals
In a similar vein, we can split U in short-and long-range contributions, so that where the long-range part will contain the classical FF interaction between electrons that can be approximated analytically, while the short-range part will contain all other interactions including many-body effects.
Similarly to the case above, we expect that (i) long-range two-electron integrals should be very small unless a overlaps with b and a ′ overlaps with b ′ , (ii) we can truncate the multipolar expansion at the monopole level, and (iii) the electrostatic interactions take place in a medium characterized by the high-frequency dielectric tensor of the material at the RED. Under these conditions we choose U lr In order to avoid the divergence of this term, we assume that all one-body integrals (a = a ′ ) are fully included in the shortrange part, U sr aba ′ b ′ . Assigning the Wannier functions to their closest nucleus, and summing over all the atoms in the lattice, we find that the total two-electron long-range energy adds to where q λ is the change in charge of the atom λ when compared to the RED state [Eq. (54)]. Thus, the long-range part of U simply updates the one-electron FF potentials due to the charge transfers between atoms. We would like to stress again that the separation in longand short-range parts does not involve any approximation; indeed, effects usually considered important in many physical phenomena, like, e.g., the anisotropy of the Wannier orbitals at short distances [59], are included in U sr .

F. Electron-lattice coupling
The system's geometry determines the reference density n 0 (⃗ r) as well as the corresponding Hamiltonian. In our scheme, such a dependence of the model parameters on the atomic configuration is captured by the electron-lattice coupling terms.
Let us consider the lattice dependence of the one-electron integrals γ [Eq. (23)]. In Sec. III E 1, these parameters were split in short-and long-range contributions. The explicit dependence of the long-range part with the distortion of the lattice is clearly seen in Eq. (55), where the electric dipole that enters in the multipole expansion of the far field potential [Eq. (56)] depends linearly with the atomic displacements, as computed with respect to the RAG. As regards the dependence of γ sr ab on the atomic configuration (see Fig. 5), we include it by expanding where quantifies the relative displacement of atoms λ and υ. In addition, ⃗ f and ← → g are the first-and second-rank tensors that characterize the electron-lattice coupling, closely related to the concept of vibronic constants [60].
We have checked that including quadratic constants is enough to describe typical changes in the value of γ with the geometry. For example, we have inspected the γ parameters associated with the oxygen 2p-like WFs of SrTiO 3 , i.e., with the valence band of the material, and plotted in Fig. 6 the three that are most sensitive to structural deformations: They correspond to the diagonal elements of the σ and π functions centered on the oxygen ions [see Fig. 3(a)] and a π − π off-diagonal term. We find that, if we use a quadratic expansion to describe such a structural dependence, the errors are smaller than 1% over a wide range of distortions, up to 0.3Å. Hence, given the strong changes occurring in the hybridization of ferroelectriclike materials like SrTiO 3 , we consider that the approximation employed in Eq. (62) should be reasonable for most systems.
Moreover, in the cases studied so far, we have found that the quadratic constants are typically much smaller than the linear ones; further, among the quadratic constants, the diagonal ones are clearly dominant. Hence, in Eq. (62) we restrict the expansion to two-atom terms, so that The physical meaning of ⃗ f ab,λυ is particularly obvious when a = b: it represents the force created by an electron occupying the WF χ a over the surrounding atoms [see Figs. 5(a) and 5(b)]. Such a parameter is key to quantify phenomena like the Jahn-Teller effect in solids [60] or polaron formation [61].
Off-diagonal terms in ⃗ f describe the mixing of two WFs upon an atomic distortion and thus quantify changes in covalency [see Figs. 5(c) and 5(d)]. They can be identified with the pseudo Jahn-Teller vibronic constants and are thus relevant to a wide variety of phenomena including ferroelectricity [60], spin-crossover [62], and spin-phonon coupling [63].
Finally, the geometrical dependence of the two-electron parameters, U and I , can be included in our model in a similar way. Nevertheless, since these terms are not explicitly dependent on the potential created by the ions, their value can be expected to be less sensitive to changes in the atomic configuration. Hence, in this paper, and in analogy to what is customarily done in model Hamiltonian and DFT+U approaches [48], we will neglect such effects keeping in mind that these may be introduced, if necessary, in the future.

G. Total energy
Replacing the expressions for the one-electron [Eq. (21)], and two-electron [Eq. (39)] integrals into Eq. (5) for the total energy, we get or, equivalently, in terms of the spin-up and spin-down densities, Now we introduce the decomposition of the γ [Eqs. (57) and (58)] and U [Eqs. (59) and (60)] parameters into long and short-range parts and gather together all the long-range terms to obtain Note that for the case of a non-spin-polarized system (D I ab = 0) the expression for the total energy reduces to

H. Self-consistent equations
As it is clearly seen in Eq. (67), the total energy in our formalism depends on the deformation occupation matrix defined in Eq. (19) and later generalized for the case of spin-polarized systems in Eqs. (30) and (31). This quantity is directly related with the deformation charge density, i.e., with the difference between the total charge density and the reference electronic density. It can be computed from the coefficients of the Bloch wave functions in the basis of Wannier functions, which are thus the only variational parameters of the method.
Solving for the ground state amounts to finding a point at which the energy is stationary upon variations in the electronic density, n(⃗ r). Following a textbook procedure [2][3][4][5], we obtain a set of self-consistent conditions analogous to the Kohn-Sham equations [ where, as defined above, ε s j ⃗ k , is the j th band energy at wave vector ⃗ k for the spin channel s. The corresponding Hamiltonian matrix h s where h s ab is the real-space Hamiltonian Note that this is a mean-field problem fully equivalent to that of the Hartree-Fock approach, and it must be solved selfconsistently. The practical procedure for finding the solution is straightforward: Given an initial guess for the deformation occupation matrix (D s ab ), we compute the corresponding mean-field Hamiltonian (h s ab ); from the diagonalization of this matrix we obtain a new deformation occupation matrix, and the procedure is iterated until reaching self-consistency. Note that electrostatic effects are accounted for by computing the long-range part of the γ and U parameters; this is our scheme's equivalent to solving Poisson's equation, as customarily done in DFT and other approaches. Finally, note that in cases in which the system does not present any electron excitation, i.e., whenever the full density is equal to the reference density and we have D s ab = 0, no self-consistent procedure is needed to obtain the solution.

I. Forces and stresses
Forces and stresses can be computed by direct derivation of the total energy [Eq. (65)] with respect to the atomic positions and cell strains. After some algebra, the result for the forces is where λ denotes a specific atom in a certain cell; here we assume that electron-lattice couplings are restricted to the oneelectron terms. The derivative of E (0) can be computed directly and exactly from the force-field on which our model is based. The deformation occupation matrix D U ab depends on the eigenvector coefficients and occupations. However, its derivative with respect to the atomic displacement is not required, since the energy is stationary with respect to these coefficients and occupations on the Born-Oppenheimer surface, and the Hellman-Feynman theorem guarantees that their first-order variation will not modify the total energy, and therefore will not affect the forces. Moreover, due to the orthogonality of the basis set used, no orthogonality forces need to be included, as it is the case when using a basis of nonorthogonal atomic orbitals (see Appendix A of Ref. [31]).
It is interesting to further inspect the similarity between the second term in our forces and the Hellmann-Feynman result [3], as (via a Fourier transform) ⃗ ∇ λ γ ab is analogous to ⟨ψ j ⃗ k | ⃗ ∇ λĥ0 |ψ j ⃗ k ⟩, and D U ab plays the role of the occupations o j ⃗ k . This connection should be considered with caution, though, as our forces have a dominant contribution from the RED state, which is also included in the Hellmann-Feynman expression. It is also interesting to note that, if we included the dependence of U (and I ) on the nuclear positions, we would have a Pulay term in Eq. (72) [64], reflecting the change of the WF basis set with the atomic displacements. Now we calculate the stress tensor in an analogous way. We adopt the standard definition [3] where V is the volume of the real-space cell and ∂ ′ denotes derivative keeping the fractional coordinates of the atoms in the system constant. We notice that there are only three terms in the energy that depend explicitly on the strain tensor ← → η , namely, the RED energy E (0) , the short-range one-electron term, γ sr ab , and the electrostatic energy. Thus, we have where E elec corresponds with the electrostatic energy as written in the fourth contribution to the total energy in Eq. (67). As in the case of the forces, the E (0) derivative is computed from the force field that describes the RED state. Similarly, the calculation of the last, electrostatic term can be achieved via Ewald summation techniques (see, e.g., Ref. [65]). The only term that requires further manipulation is the derivative of γ sr ab , Eq. (62), with respect to the strain, which yields As in the case of the forces, the similarity between this result and the Hellmann-Feynman expression is apparent.
To end this section, let us stress that only excited electrons/holes, which render D ab ̸ = 0, create forces and stresses not included in the underlying force-field described by E (0) . In fact, in the typical case, the dominant contribution to both forces and stresses will come from the derivative of E (0) , with corrections that are linear in the difference occupation matrix.

J. Practical considerations
So far we have introduced a method for the simulation of materials at a large scale. We have presented the basic physical ingredients (reference atomic geometry, reference electronic density, deformation density, etc.), that allow us to approximate the DFT total energy, forces and stresses.
In this section we discuss some practicalities involved in the implementation of this method in a computer code to perform actual calculations. Of course, different implementations are in principle possible; here we briefly describe some details pertaining to our specific choices, which should be illustrative of the technical issues that need to be tackled.

Definition of the RED
The formulation above is written in terms of differences between the actual and reference states of the system in a completely general way. However, from a practical point of view, an appropriate choice of the RED, n 0 (⃗ r), is a necessary first step towards an efficient implementation of our method.
The most important ingredient to define n 0 is the reference occupation matrix that, following Eq. (17), amounts to where o (0) j ⃗ k and c (0) ja ⃗ k characterize the RED. While it would be possible to use the d (0) ab result computed from first principles to perform second-principles simulations, in the following we shall simplify this expression in order to obtain a more convenient form. Note that the reference occupation matrix satisfies d (0) ab = 0 for a and b belonging to different band manifolds. (By definition, if a and b belong to different bands, they cannot appear simultaneously in the expansion of a particular Bloch state [Eq. (14)], and the corresponding d (0) ab [Eq. (77)] will vanish.) It is thus possible to rewrite Eq. (77) and split the sum over states in two, one over manifolds J and a second one over bands within a manifold.
After having established this property, we impose that all the bands j that belong to the same manifold J have the same occupation in the RED where ω ⃗ k is the weight of each ⃗ k-point in the BZ. As we assume an homogeneous sampling in reciprocal space, ω ⃗ k = N −1 k , where N k is the total number of points in our BZ mesh. Thus, for example, in a diamagnetic insulator [see Fig. 3

(a)]
(where the valence and conduction band always belong to different manifolds) we would choose the occupation for the reference states so that all the valence bands are fully occupied while all conduction bands are completely empty. In this way the reference electronic density for a diamagnetic insulator is simply the ground state density.
For metals [ Fig. 3(b)], and magnetic insulators [ Fig. 3(c)], where a disentanglement procedure [38] has to be carried out to separate the desired bands from others with which they are hybridized in a given energy window, the choice is not so simple. In such cases we distribute all the electrons of the entangled bands equally among the bands in the manifold. For example, in the case of metallic copper [ Fig. 3(b)], which has an electronic configuration 3d 10 4s 1 where the five 3d functions cross with the 4s-like band, we would distribute the eleven electrons over the six bands taking o J = 11/6. On the other hand, for NiO where we have used the fact that o (0) a = o J for all WFs in the manifold.
We would like to stress that this approach simply allows for a more efficient computational method since we can still retrieve the full electron distribution in space by substituting the first-principles WFs into Eqs. (16) and (18).

Deformation electron density
From the definition of the charge density in terms of the density matrix, Eq. (16), and the orthogonality of the Wannier basis functions, it trivially follows that the total number of electrons in the system is the trace of the density matrix, Therefore, the trace of the deformation matrix gives the number of extra electrons or holes that dope the system. From the very definition of the deformation matrix, we can deduce that if its diagonal element, is negative (positive), that means that we are creating holes (inserting electrons) on that particular state, as illustrated in Fig. 1. The fact that most of these electron/hole excitations take place around the Fermi energy has a very important consequence with regards to the efficiency of the method. In order to calculate D aa , and the total energy [see Eq. (67)], we do not need to obtain all the eigenvalues of the oneelectron Hamiltonian, but just those around the Fermi energy. This opens up the possibility to use efficient diagonalization techniques that allow a fast calculation of a few relevant eigenvalues, e.g., Lanczos. This approach allows us to speed up the calculations in a very significant manner. (The diagonalization of the full Hamiltonian matrix is one of the main computational bottlenecks in electronic structure methods.) Along these lines, the possibility of obtaining linear scaling within our method will be discussed in a future publication.

IV. PARAMETER CALCULATION
The method presented above allows for the simulation of very large systems under operation conditions assuming that a few parameters describing one-electron and two-electron interactions, as well as the electron-lattice couplings, are known beforehand. For the sake of preserving predicting power, it is important to compute those parameters from first principles.
All the electronic parameters of our models have welldefined expressions [see Eq. (23) for γ ab , Eq. (40) for U aba ′ b ′ , and Eq. (41) for I aba ′ b ′ ], whose computation requires only the knowledge of the Wannier functions, the one-electron Hamiltonian, and the operators involved in the double integrals, all of them defined for the RED. Since the chosen basis functions are localized in space, the required calculations could be performed on small supercells. Such a direct approach to obtain the model parameters is thus, in principle, feasible.
Note that there has been significant work to calculate related integrals from first principles, as can be found, e.g., in Refs. [42,48,49,51,[66][67][68][69][70][71][72][73]]. Yet, we feel that most of these approaches are too restrictive for the more general task that we pursue in this paper. For instance, the focus in the previous references is placed on strongly correlated electrons in a single center, while we are also interested in multicenter integrals.
A significant effort would thus be required to implement the calculation of the more complex interactions, including all the potentially relevant ones, and developing tools to derive minimal models that retain only the dominant parameters and capture the main physical effects. Note that, in a typical system, the number of potentially relevant integrals will be very large. In fact, the presence of four-index integrals like U aba ′ b ′ and I aba ′ b ′ is the reason why Hartree-Fock schemes scale much worse than DFT methods with respect to the number of basis functions in the calculation [∼O(N 5 ) vs ∼O(N 3 ), respectively]. Hence, at the present stage we have not attempted a direct first-principles calculation of the parameters, which is a challenge that remains for the future. Instead, we have devised a practical scheme to fit our models to relevant first-principles data.

A. Parameter fitting
Our procedure comprises several steps.

Training set
First, we identify a training set (TS) of representative atomic and electronic configurations from which the relevant model parameters will be identified and computed. For example, the training set for a magnetic system should contain simulations for several spin arrangements, so that the mechanisms responsible for the magnetic couplings can be captured. Additionally, if we want to study a system whose bands are very sensitive to the atomic structure, the training set should contain calculations for different geometries so that this effect is captured. Alternatively, if we want to describe how doping affects the physical properties of a material, then different DFT simulations on charged systems should be carried out [74], etc.
Let us note that it is typically possible to restrict the TS to atomic and electronic configurations compatible with small simulation boxes. This translates into (and is consistent with) the fact that, when expressed in a basis of localized WFs, the nonelectrostatic interactions in most materials are short ranged.
We will use N TS to denote the total number of TS elements, noting that we will run a single-point first-principles calculation for each of them. Further, N RAG is the number of TS configurations that correspond to the reference atomic geometry.

Filtering the training set
where the T matrices are unitary transformations that convert the first-principles eigenstates into Bloch-like waves associated to specific (localized) WFs. These transformations can be obtained by employing codes like WANNIER90 [66], which implements a particular localization scheme, i.e., a particular way to compute optimum T matrices [42,45].
Once the h s ab (i) Hamiltonians are known, we can identify the pairs of WFs with a large enough interaction and which need to be retained in the fitting procedure. In practice, we introduce a cutoff energy δε h such that h s ab (i) > δε h , for at least one i in the TS, defines the Hamiltonian matrix elements to be retained.
(Diagonal elements, h s aa , are always considered independently of their value.) This condition allows us to identify the WF pairs (a, b) to be included in the fitting procedure, regardless of the geometry or spin arrangement.
In Fig. 7 we compare the full first-principles bands for SrTiO 3 and NiO with those obtained from models corresponding to different energy cutoffs. For all the δε h values considered, we also indicate the number of parameters in the resulting models. This allows us to estimate the size of the model (and associated computational cost) needed to achieve an acceptable description of the band structure.

Identifying most relevant model interactions
Our models, even though we truncate them at second order of the expansion in Eq. (5), contain a daunting number of electron-electron interaction parameters. Constructing an actual model usually involves further approximations regarding the spatial range of the interactions, the maximum number of different bodies (WFs) involved, etc. Hence, we need a procedure to identify the simplest models that can reproduce the first-principles TS data with an accuracy that is sufficient for our purposes.
The scheme we have implemented is based on a very simple logic: We start from a certain complete model that may contain, in principle, all possible one-electron, two-electron, and electron-lattice parameters. We can then fit such a model to reproduce the one-electron Hamiltonians {h s ab (i)} of our TS within a certain accuracy. Typically, by doing so, and by systematically exploring different combinations of parameters in the model, we can identify the simplest interactions (i.e., those that are shortest in range, involve fewest WFs, etc.) sufficient to achieve the desired level of accuracy; in other words, in this way we can identify noncritical couplings that just render the fitting problem underdetermined, do not improve the accuracy of the model, and can thus be disregarded. Naturally, this split between relevant and irrelevant couplings is strongly dependent on the choice of the training set, which should be complete enough to capture the physical effects of interest.
To better understand how the scheme works, consider the one-electron integrals γ ab in the case of a nonmagnetic material like SrTiO 3 . These parameters will be the only ones entering the description of the band structure of the RAG in the RED state. Hence, we can fit them directly by requiring our model to reproduce the Hamiltonian h ab of this particular, reference state with a certain accuracy.
More generally, the one-electron Hamiltonian corresponding to a TS configuration will reflect electronic excitations departing from the RED state. More precisely, we can recall Eq. (71) to write where we restrict ourselves to TS configurations at the RAG, so that no electron-lattice term appears in this equation. It is convenient to isolate the U contribution by defining and its average over all the RAG configurations in the training set,h Analogously, the antisymmetrization of the Hamiltonian matrix elements with respect to the spin yields We expect that the most important U aba ′ b ′ and I aba ′ b ′ parameters will be, respectively, those involving WF pairs whose corresponding h U ab (i) and h I ab (i) are most strongly dependent on the TS state. Hence, we introduce the twoelectron cutoff energy, δε ee , and retain only the (a,b) pairs that satisfy, for at least one TS configuration, at least one of the following conditions: or Note that we gauge the h U ab matrix elements with respect to thē h U ab average values so that the corresponding cutoff condition does not depend on the one-electron couplings γ .
Once we have selected all the {(a,b)} pairs that fulfill such criteria, we can build the list of potentially relevant U and I constants to be considered in the fit. Note that the number of free parameters is usually reduced by the fact that the U aba ′ b ′ and I aba ′ b ′ integrals are invariant upon permutations of the (a,b,a ′ ,b ′ ) indexes. In some cases, and in spite of the reduction of parameters due to symmetry, the list of relevant interactions is excessively long and needs to be further trimmed to successfully carry out the fitting. In such situations we introduce a third cutoff, δD, that operates over the difference occupation matrix to select the interactions associated to important changes of the electron density. When doing so we only accept U constants for which at least one pair of the associated indexes fulfills and the corresponding expression for I Let us also note that the most relevant γ ab parameters are trivially identified when we filter the TS one-electron Hamiltonians as described above.

Fitting the RAG model
Once our list of relevant γ , U , and I parameters is complete, we fit them to reproduce the {h s ab } matrix elements above the δε h energy cutoff introduced previously. We have found it convenient to perform the fit in several steps, so that different types of parameters are computed separately. More precisely, we first fit the U parameters by requesting that our model reproduces the h U ab (i) −h U ab matrices [Eqs. (86) and (87)]. Analogously, we obtain the I constants by fitting to the h I ab (i) matrices [Eq. (88)]. Importantly, both of these fits are independent of the one-electron integrals, and have typically yielded well-posed, overdetermined systems of equations in the cases we have so far considered. Finally, we obtain the γ parameters from the fitted U 's directly from Eq. (86). Direct comparison of the modeled bands with those obtained from the full first-principles {h s ab (i)} set provides an estimate of the goodness of the model (see the example in Sec. VI B and, particularly, Fig. 9).
Note that, alternatively, one could try a direct fit of all the γ , U , and I parameters to the real-space Hamiltonians, using Eq. (71). However, we typically find that this strategy leads to nearly-singular problems in which very different solutions lead to comparably good results. In the general case, such a difficulty may be mitigated by extending the TS. However, here we adopted the simple and practical procedure described above, which permits a numerically stable method that yields accurate and physically sound models.
To end with this section we would like to stress that the γ ab constants obtained with this procedure contain both the shortand long-range contributions described above [Eq. (58)]. In order to isolate γ sr ab , we simply subtract the corresponding electrostatic contribution [Eq. (57)] from the determined, full γ ab value. In order to calculate the electrostatic contribution [Eqs. (55)-(56)] we need first-principles results for the Born charge tensor, ← → Z * λ , and the high-frequency dielectric tensor, ← → ϵ ∞ , that can routinely be obtained for systems where the RED is insulating [58].

Relevant electron-lattice interactions
As above, we assume that the deviations from the RAG only affect the one-electron integrals γ , and not the U and I parameters. We further assume that such a dependence on the atomic structure is given by the linear and quadratic electronlattice constants ⃗ f and ← → g introduced in Eq. (62).
The selection of the most important electron-lattice couplings is performed by observing how much a particular matrix element h s ab changes with a particular distortion of the lattice with respect to the RAG. To quantify this change, we need to compare pairs of configurations i and i ′ that correspond to the same electronic state (e.g., to the same spin arrangement, to the same amount of electron/hole doping, etc.) but differ in their atomic structure. More precisely, one of the configurations must correspond to the RAG state (i), while the other one (i ′ ) is characterized by a distortion given by {⃗ u λ (i ′ )}. For simplicity, here we will restrict to distortions involving only one displacement component of one atom, so that we only have one specific u λα (i ′ ) ̸ = 0, where α labels the spatial direction. We then consider that a particular atom λ participates in the electron lattice affecting the h s ab element if where δf e−l is a new cutoff. Note that, for a large enough distortion |u λα |, this condition pertains to both the linear ( ⃗ f ) and quadratic ( ← → g ) electron-lattice interactions. Yet, since we activate a single atomic displacement at a time, in the case of ← → g we are only probing the diagonal elements. Restricting ourselves to the diagonal elements of ← → g is justified by the observation that, in the systems we have so far studied, those are the only significant ones. At any rate, the scheme can be trivially extended to check a possible contribution of offdiagonal terms.

V. IMPLEMENTATION OF THE ALGORITHM: THE SCALE-UP CODE
We have implemented this new method in the SCALE-UP (Second-principles Computational Approach for Lattice and Electrons) package, written in Fortran 90 and parallelized using Message Passing Interface (MPI). Presently, this code can perform single-point calculations, geometry optimizations, and Born-Oppenheimer molecular dynamics using either full diagonalization or the Lanczos scheme mentioned above.
The energy of the reference state, E (0) , is obtained from model potentials like those introduced by Wojdeł et al. [30], which are interfaced with SCALE-UP. We have also developed an auxiliary toolbox (MODELMAKER) for the calculation of all the parameters defining E (1) and E (2) , using as input DFT results for one-electron Hamiltonians in the format of WANNIER90 [66]. As shown in Sec. VI, these implementations can be used to create models that match the accuracy of the DFT calculations at an enormously reduced computational cost, opening the door to large-scale simulations (up to tens of thousands of atoms) of systems with a complex electronic structure, using modest computational resources.
The input to the code is based on the flexible data format (fdf) library used in SIESTA [31] and contains several pythonbased tools to plot bands, density of states, geometries and other properties.

VI. EXAMPLES OF APPLICATION
In order to illustrate the method, we will discuss its application to two nontrivial systems with interactions of very different origin. The first example consists in the calculation of the energy of a Mott-Hubbard insulator, NiO, for different magnetic phases. Our goal here is to show that the method can be used to deal accurately with complicated electronic structures including phenomena like magnetism in transition metal oxides. In this example we will also show how the method can tackle rather large systems (2000 atoms) that are at the limit of what can be done with first-principles methods today, reducing the computational burden by orders of magnitude.
The second application involves the two-dimensional electron gas (2DEG) that appears at the interface between band insulators LaAlO 3 and SrTiO 3 [75]. We will not discuss here the origin of the 2DEG, which has been treated in great detail in the bibliography [28,[75][76][77]. Rather, we will check whether our approach can predict the redistribution of the conduction electrons at the LaAlO 3 /SrTiO 3 interface, and the accompanying lattice distortion, as obtained from first principles. Thus, this example will showcase the treatment of electron-lattice couplings and electrostatics within our approach.

A. Details of the first-principles simulations
We construct our models following the recipes described in Sec. IV, and the first-principles data are obtained from smallscale calculations with the VASP package [78][79][80]. The local density approximation (LDA) to density-functional theory is used to create the TS data for SrTiO 3 . The calculations for NiO are also based on the LDA, but in this case an extra Hubbard-U term is included to account for the strong electron correlations [81], as will be discussed below. We employ the projectoraugmented wave (PAW) scheme [80] to treat the atomic cores, solving explicitly for the following electrons: Ni's 3s, 3p, 3d, and 4s; O's 2s and 2p; Sr's 3s, 3p, and 4s; and Ti's 3s, 3p, 3d, and 4s. The electronic wave functions are described with a plane-wave basis truncated at 300 eV for NiO and at 400 eV for SrTiO 3 . The integrals in reciprocal space are carried out using -centered 4 × 4 × 4 k-point meshes in both cases.

B. NiO magnetic couplings
Transition metal oxides are very interesting as they present optical, magnetic, and structural properties that are, very often, tightly coupled with each other. This fact, together with the large variety of functional properties that they can display, makes them a big focus of attention in both basic and applied materials science. From a theoretical point of view their study is complicated, mostly because of the strong correlations associated to the electrons in the compact d shell (especially, those of first-row 3d transition metal ions) and the frequent presence of many competing magnetic phases. Naively, one may expect most of these oxides to be metallic due to their open-shell nature while, in fact, many are insulators. This problem strongly affects computational methods; in fact, most common approaches, like DFT with local or semilocal exchange-correlation functionals, often fail to correctly reproduce the magnetic or insulating properties of these compounds. Simulations at this level of theory yield too diffuse states with underestimated interactions. The key to simulate successfully these materials lies in the way electron-electron interactions are handled. A panoply of methods have been developed to deal with this issue from first principles, ranging from the inclusion of a Hubbard-like term in the Hamiltonian in the so-called DFT+U methods, to more sophisticated schemes with dynamically screened interactions, such as the GW approximations [50,51].
For this first application of our scheme, we have chosen a simple transition metal oxide, NiO, a staple example of many new electronic structure simulation methods. Our goal is to show how our second-principles scheme can handily be used to compute the properties of strongly correlated materials, based on parameters obtained from first-principles LDA+U simulations. In particular we will study the band structure and magnetism of this archetypical binary oxide.

Model parameters
The geometric structure of NiO shows some subtle and not fully understood distortions associated to its magnetism [82,83], but this issue is beyond the scope of the present paper. In order to keep the model as simple and illustrative as possible, we neglect the lattice degrees of freedom in this case. As regards the spin order, neutron diffraction experiments [84] show that the ground state corresponds to the so-called AF2 phase, where planes of spin-up and spin-down polarized nickels alternate along the ⟨111⟩ direction of the conventional cell [see Fig. 8(c)]. Further experiments [85,86] evidence that the AF2 to paramagnetic transition is of second order and occurs at a Nèel temperature T N = 524 K.
The simulations of NiO are carried out in the experimental rocksalt cell, with a lattice constant of 4.17Å [87]. This cell is compatible with several spin arrangements ranging from fully ferromagnetic [FM; Fig. 8(a)] to various antiferromagnetic (AF) ones [ Figs. 8(b) and 8(c)]. Our TS includes the ground state spin arrangement (AF2) as well as the ferromagnetic one, which we choose because it represents a relevant limit case for spin-spin interactions. We use the LDA+U approach introduced by Dudarev et al., [81] with an effective U value of 7 eV, applied only on the 3d orbitals of Ni. These calculations indicate that the AF2 solution is more stable than the FM one by 89 meV per formula unit (f.u.).
Ligand-field theory predicts the magnetism in this lattice to be the result of the half-filled e g shell of the octahedrally coordinated Ni 2+ ion [see Fig. 8(d)]. Thus, we expect the levels around the Fermi energy to have this character. After calculating the electronic structure from first principles within the LDA+U, we find that the top valence and bottom conduction bands are composed of several strongly entangled states, as shown in Fig. 3(c). Thus, we project our WFs seeking to disentangle orbitals participating in the valence band [Ni(t 2g ), Ni(e g ), and O(p)] from others in the conduction band; to do this we use the tools provided within the WANNIER90 package. A graphical representation of the resulting orbitals [ Fig. 3(c)] clearly shows that we are able to isolate bands with the expected chemical character: The isosurfaces of the maximally-localized WFs (MLWFs) at the right hand side of Fig. 3(c) clearly resemble the shape of the O(p), Ni(d xy ), and Ni(d 3z 2 −r 2 ) orbitals for the valence band, and of the Ni(d 3z 2 −r 2 ) orbital for the bottom of the conduction band. Given the strong entanglement of these bands, we use a reference occupation for the construction of our model that is obtained by populating equally all of them. As discussed in Sec. III J 1, this amounts to assuming o J = 7/8 = 0.875 electrons per band and spin channel.
At this point we start with the analysis of the Hamiltonian as described in Sec. IV, using the {h s ab (i)} set obtained after the disentanglement procedure. First, we seek to choose a reasonable value of δε h that allows us to describe accurately the system's bands without including an excessive number of γ ab terms. As can be seen in Fig. 7, δε h = 0.05 eV is a reasonable choice; this involves the use of 71 γ terms per f.u.
In order to decide the values of δε ee and δD that will determine the U and I parameters considered in the fit, we first examine the occupation difference of each of the WFs, D aa . Let us recall that the diagonal elements of the deformation occupation matrices are defined as the difference between the reduced density matrix computed at the LDA+U level for the corresponding configuration in the TS and the reference one. In the FM phase, the bands with e g character for the majority spin channel are expected to be fully occupied (d aa = 1), while they should be empty for the minority spin (d aa = 0). Therefore, for WFs with e g character, we expect to have D aa = 0.125 for the majority spins and D aa = −0.875 for the minority spins, for both Ni atoms. For the antiferromagnetic configuration, as the spin in one of the Ni atoms is reversed, we expect the same  behavior for D aa as in the FM solution for one of the nickels, while the character of the majority and minority spins must be exchanged for the second Ni. In both cases (FM and AF2), the t 2g bands are fully occupied, so that the corresponding D aa should be 0.125. These tendencies are well reproduced in our TS configurations as can be seen in Table I, where the occupations are obtained using the {h s ab (i)} set obtained after applying δε h = 0.05 eV filter. The differences with respect to the ideal ionic values are due to the chemical bonding between Ni and O.
If the results for the FM and AF2 configurations are compared, we can see that the only significant change pertains to the occupation of the e g -like orbitals of the Ni ion whose spin is flipped: The majority (0.125) and minority (−0.756) values of the difference occupation are basically exchanged as we move from the FM to the AF2 calculation, as expected from the localization of the magnetic moment over these orbitals in Ni 2+ ions (see Fig. 8). This is an indication that, in order to capture the magnetic interactions in this system, only electron-electron interactions involving e g -type WFs are necessary. Moreover, we can observe how, on one hand, the average orbital occupations vary very little (essentially by 0.005, 0.001, and 0.003 electrons for the Ni(e g ), Ni(t 2g ), and O(p)-like WF, respectively) between spin configurations. This fact translates into a similar value of D U [Eq. (37)]. As can be seen in Eq. (67), if D U is the same for the different configurations of the TS, its contribution to the total energy is constant, i.e., it does not play any role in the calculations of the relevant energy differences. On the other hand, the spinup/spin-down differences of occupation are strongly changing between the FM and AF2 configurations. This indicates that only Stoner-type (I ) interactions are relevant to describe the relative stability of the magnetic phases in the training set. As can be seen in Table II, by playing with δε ee and δD it is straightforward to confirm that the e g −e g interactions drive magnetism in this system, while Ni(t 2g )-like and O(p)-like energy levels have a secondary role. Nevertheless, including the latter interactions is necessary to accurately describe the bands.
In Fig. 9 we show the second-principles computed bands for two different set of parameters: (i) the first one, obtained after a filtering the electron-electron interactions with a threshold of δε ee = 1.10 eV, was selected to include only couplings between Ni(e g )-like WFs. The results obtained with this set of parameters are labeled SP-LDAU-Ni(e g ); (ii) the second one, obtained with a threshold of δε ee = 0.20 eV, corresponds to a case in which interactions between Ni(e g )-like and Ni(t 2g )-like  of the table) and second-principles calculations (lower part of the table). The latter have been modeled after the LDA+U calculations (highlighted) and, as can be seen, converge towards the results obtained with this method when reducing δε ee . HSE stands for the hybrid exchange and correlation functional proposed by Heyd-Scuseria-Ernzerhof [90], PSIC for pseudo self-interaction-correction, ASIC for atomic self-interactioncorrection, and GGA for generalized gradient approximation. WFs at the same atom, as well as between Ni(e g )-like and nearest-neighboring O(p)-like WFs, are also included. The results are labeled SP-LDAU-Ni+O. As can be seen, the bands for the FM state are better reproduced in the second case, because of the correction of the diagonal spin-up/spindown Ni(t 2g ) and O(p) energies. These diagonal Hamiltonian matrix elements, which determine the center of mass of the corresponding bands, vary with the Ni(e g ) occupation as expressed in Eq. (71). Indeed, we can estimate the maximum error for the h s ab terms as

Method
where h s ab (i) is a matrix element directly obtained from the first-principles TS and h s ab is computed from Eq. (71) for a given set of parameters. This maximum error reduces from 0.651 eV in the SP-LDAU-Ni(e g ) case to 0.132 eV for SP-LDAU-Ni+O. We also considered a third, intermediate case with δε ee = 0.5 eV, where Ni(e g )-Ni(t 2g ) interactions are included but those with oxygen orbitals are neglected [SP-LDA-Ni(3d) in Table II]. The associated maximum error is 0.29 eV for such a choice.
It is interesting to note that our TS does not contain enough information to fit reliably all the Ni(e g )-Ni(e g ) interactions compatible with our filters. More precisely, we find that there are only two relevant I ab,a ′ b ′ constants: one related to the self-energy of the e g states (I aa,aa , where a is a e g -like basis function), and a second one quantifying the interaction between the two e g states in the same atom (essentially, the exchange interaction known as Hund's coupling). It is clear that our TS is not suitable to distinguish between such interactions. In both the FM and AF2 phases, the Ni 2+ ions display a S = 1 spin configuration; yet, the interplay between self-interaction and Hund coupling only appears when trying to differentiate between the high-(S = 1) and low-(S = 0) spin intra-atomic states. We checked whether such an indeterminacy affects the energy difference between the FM and AF2 phases as obtained from the model. To do so, we first add a Hubbard-U constant associated to the self-energy of the e g WFs, and make it equal to the corresponding I constant (i.e., we impose I aaaa = U aaaa ). In this way the self-interaction of an electron placed in one of these orbitals is U − I = 0, while the interaction between two electrons that only differ in their spin is U + I = 2I = 2U . Then, we vary this self-interaction parameter between 0 eV and 6 eV, optimize the interaction between different e g orbitals to reproduce the bands, and calculate the FM-AF2 energy gap. We observe that the energy difference is quite insensitive to the value of I aaaa , varying by less than 5% in the explored range. Hence, we simply take I aaaa = 2 eV to fix the indeterminacy in the model.

Results
Magnetism in rocksalt structures is usually described using a Heisenberg Hamiltonian with coupling constants between first-(J 1 ) and second-(J 2 ) nearest neighbors [86,[88][89][90][91][92]. These constants can be obtained from the energy differences between different spin arrangements by solving the equation system: which involves the spin arrangements of Fig. 8. E ref stands for the energy of a reference phase. We employ our models to compute the energy of these phases, using a 2 × 2 × 2 supercell containing 16 atoms, as sketched in Fig. 10. After converging the calculations, we plot the spatial distribution of D I (using the spatial representation of the WFs in our basis) to check that the obtained solutions correctly correspond to the FM, AF1, and AF2 spin arrangements. These plots are shown in Fig. 10, where it can be seen that the electron distribution in the simulations perfectly matches the spin orderings sketched in Fig. 8. We can now check the numerical results for the phase energies obtained with the three second-principles parametrizations [denoted, respectively, SP-LDAU-Ni(e g ), SP-LDAU-Ni(3d), and SP-LDAU-Ni+O in Table II] and compare them to our full DFT+U result and data in the literature (see, e.g., Ref. [89]).
We find that the coupling constants computed from our models compare quite well with the first-principles results. Indeed, we find the J 2 , running along the 180 • Ni-O-Ni bridge, to be much stronger than the J 1 coupling along the 90 • Ni-O-Ni path. It is worth noting that a parametrization as simple as that of the SP-LDAU-Ni(e g ) model captures this essential feature already. Then, when we include I couplings between Ni(e g ) and Ni(t 2g ) WFs, we obtain a very similar result, with a very small J 1 = −0.04 meV. Finally, when we include electronelectron couplings with the oxygen orbitals, we get a value for J 1 that is very close to the first-principles result.

C. Electron gas at the LaAlO 3 /SrTiO 3 interface
Now we tackle the well-studied electron gas appearing at the interface between LaAlO 3 (LAO) and SrTiO 3 (STO). The origin of the 2DEG has been intensively debated in the literature [77,93]. Here we are going to consider an idealized defect-free interface in which the driving force for the 2DEG is the so-called polar catastrophe that was proposed originally [75,76], which arises from the charge discontinuity between LaAlO 3 and SrTiO 3 when the bilayer is grown along the (001) pseudocubic direction of the perovskite lattice [94]. In such a case, the occurrence of the metallic state strongly depends on the electrostatic boundary conditions on each side of the interface. [28] Let us look at them in some detail to establish the basic elements of the calculation.
From simple electrostatic arguments we know that where D LAO and D STO are the normal components of the displacement field in LaAlO 3 and SrTiO 3 , respectively, and σ free is the free charge density at the interface between the materials. Hence, depending on the particular values of D LAO and D STO (which can be controlled in a simulation by varying the charges at the open surfaces of the layers [28]), a 2DEG appears at the interface according to Eq. (96). Figure 11 This setup is ideal to test our method, since the main physical effects are related to the negative doping of SrTiO 3 , such a doping being controlled by the electrostatic boundary conditions. Further, the properties of the 2DEG (e.g., spatial extension) depend essentially on the ability of SrTiO 3 to screen these additional charges, which in turn involves the electron-lattice couplings in our models.
We simulate the LaAlO 3 /SrTiO 3 interface [ Fig. 11(c)] by considering a slab of N = N LAO + N STO 5-atom perovskite unit cells, where the first N LAO cells are occupied by LaAlO 3 and the following N STO cells by SrTiO 3 . Following Stengel [28], we do not consider the electronic details of the interface, as these were found to be of little relevance to describe the main physical features of the 2DEG. In fact, as regards the construction of our model, we treat the entire slab as if it was made of SrTiO 3 , but introducing the following modifications: (i) on the LaAlO 3 side, the levels of the conduction band are shifted up (by appropriately modifying the γ aa self-energies) so that they do not interact with those on the SrTiO 3 side, and (ii) to account for the large disparity between the LaAlO 3 and SrTiO 3 dielectric constants (the latter is around 25 times larger than the former at room temperature) we simply freeze the coordinates of the atoms on the LaAlO 3 layer at the RAG, to prevent atomic displacements from screening electric fields.
As regards the electrostatic boundary conditions, we consider that D STO = 0 and D LAO = −N e /S, i.e., we have N e electrons per unit area S doping the slab. To impose such conditions, we first freeze into the centrosymmetric structure the atomic positions of N freeze unit cells at the end of the SrTiO 3 side of the slab [see Fig. 11(c)]. Secondly, we use the image-charge method [56] [see Fig. 11(c)] to introduce an electric field from the LaAlO 3 side of the interface consistent with D LAO = −N e /S.

Model parameters
We now describe how we obtain the parameters for the SrTiO 3 layer. As already mentioned, SrTiO 3 is a nonmagnetic insulator and the RED corresponds to the ground state of the undoped system. This allows us to take the lattice potential for pure SrTiO 3 described in Ref. [30] as the E (0) term of our model [see Eqs. (5) and (6)], using the LDA-relaxed cubic phase as RAG. [We slightly modified the force field of Ref. [30], by tuning the interaction between first-nearestneighboring Ti and O pairs, to exactly reproduce a dielectric constant of 500 for the cubic phase (see Fig. 12), as obtained from LDA calculations in Ref. [28]].
We then extended the model to include the electronic states associated with the bottom of the conduction band of SrTiO 3 , which present a dominant Ti(t 2g ) character [see Fig. 3(a)]. We followed the recipe in Sec. IV to extract the γ parameters describing these bands.
Note that our focus in this application was to capture the electron-lattice effects that determine the properties of the 2DEG, and we were not concerned with electron-electron couplings beyond the LDA. Thus, we did not include U or I terms in our model, and used a TS that contains the RAG and distorted structures (with individual atoms displaced by 0.05Å, 0.10Å, and 0.15Å from their RAG positions), all of which where assumed to be in the RED state.
We then found all γ , ⃗ f , and ← → g parameters [Eq. (62)] compatible with the choices δε h = 0.05 eV and δf e−l = 1.0 eV/Å [95]. We observe that the electron-lattice constants associated to diagonal one-electron terms, γ aa , are much more sensitive to the displacement of the ions than the off-diagonal ones. The distortions that induce larger changes involve the Ti-O bond, as expected according to the long literature on covalency in ferroelectric oxides and related materials [60,[96][97][98]. Finally, we took the Born charges and high-frequency dielectric tensor used in our lattice model, E (0) , to compute the electrostatic energy associated to the electronic degrees of freedom.

Results
We now compare the results for the LaAlO 3 /SrTiO 3 interface obtained with the above-described model and the LDA results of Ref. [28], where the calculations were performed for N LAO = 5 and N STO = 12. Using these values, we carry out geometry optimizations with the constraints illustrated in Fig. 11(c). The results of these calculations are shown in Figs. 13(a) and 13(b), where we compare the electron densities from first-principles LDA and second-principles simulations for N e = 0.3 and N e = 0.5, respectively. Moreover, in Figs. 13(c) and 13(d) we show the obtained lattice distortions, in terms of the layer-by-layer rumpling, for the same cases.
We can observe that the second-principles and LDA results match well for both atomic and electronic structure. Moreover, following the discussion in Ref. [28], we checked that our model captures correctly the influence that various physical parameters (e.g., linear and nonlinear dielectric response of the lattice, etc.) have in the final result. As regards relatively small errors in the shape of the electronic density profiles, we attribute them to technical differences (pseudopotentials, etc.) in the LDA calculations of Ref. [28] and those performed to construct our models.

D. Additional considerations
To finish this section we would like to give some estimations of the computer time required to carry out the secondprinciples calculations for the systems discussed in the present work. We will focus first on NiO, as we carried out both the DFT and second-principles computations for the same unit cells and on the same computational platform, so a reliable comparison of timings should be feasible. In Table III we show the time necessary to perform a single-point calculation using a single CPU, using the same reciprocal space sampling in the LDA+U and second-principles simulations. The values for 4-and 16-ion cells clearly show the very large speed-up of our second-principles simulations when compared to standard DFT even in small cells. Looking at Table III we see that there are very small differences between the timing results between the different parameterizations associated to different levels of description of the electron-electron interactions in NiO. Finally, in order to give an estimation of the scaling of the method as it stands now (i.e., at an early stage of implementation), we carried out a calculation of a 10 × 10 × 10 periodic supercell that contains 2000 atoms. This is approximately the size limit for single-point DFT calculations, and it would require significant computational resources and a highly parallelized code. However, this simulation at the SP-LDAU-Ni(e g ) level took 6.6 hours of a single CPU, suggesting that calculations including tens of thousands of atoms are within reach using our models. Turning now to the simulation of the LaAlO 3 /SrTiO 3 interface we note that each of the geometry optimizations involving a 85-atom supercell took about 13 minutes using a single desktop CPU.

VII. RELATION WITH OTHER METHODS
As already highlighted in the Introduction, various methods have been proposed during the last few decades to bridge over time and length scales while keeping DFT accuracy. They range from empirical potentials [99], to efficient linearscaling DFT implementations permitting million atom simulations [100], or the use of the non-self-consistent (e.g., Harris) functionals [13], among others. A critical discussion of all of them is out of the scope of the present paper. Instead, we will focus on the two well-established methods most closely related to the present scheme, namely, the self-consistent charge density functional tight-binding (SCC-DFTB), and the effective Hamiltonians for lattice dynamical studies of ferroelectrics and related materials. The former puts the emphasis on the description of the electronic structure, while the latter is a purely lattice model without an explicit treatment of the electrons. Detailed connections of different aspects of the methodology have already been made at the corresponding sections where the main features of our model have been presented. Yet, once a complete view of our new scheme has been given, it seems appropriate to further compare the present work with these two previous proposals, recalling the basic features of those methods, pointing to the specific sections where the equivalent approaches have been discussed (enclosed in parenthesis below), and emphasizing the new aspects and how we can go beyond the scope of the previous schemes.
Our model is based on a simple and computationallyefficient, electron-free description of the lattice dynamical or vibrational properties of the material of interest. The first step on this direction was taken in the 90's works of Vanderbilt, Rabe, and others, who introduced first-principles model potentials (which they called effective Hamiltonians) to describe the ferroelectric phase transitions of perovskite oxides like BaTiO 3 and PbTiO 3 [16][17][18]. The effective-Hamiltonian approach involves a drastic simplification of the material, which is coarse-grained to retain only the lattice degrees of freedom associated with the ferroelectric properties. The material is then described in term of a reduced number of variables (i.e., local polar modes and cell strains) whose iterations are parametrized by means of a Taylor series of the energy around a suitably chosen reference structure. Such a scheme is physically-motivated, computationally very efficient, and its precision can be improved, to some extent, in a well-defined way. The application of the original scheme to increasingly complex oxides has shown its generality, the good transferability of the interatomic couplings among dissimilar chemical environments, and the reliability and predictive power of the models. More recently, this scheme has been generalized by some of us [30] to retain all the atomic degrees of freedom, thus removing the coarse-grain approximation. These models can be trivially formulated for any material, their accuracy is systematically improvable, their interpretation is physically transparent, and they lend themselves to very efficient schemes to compute the model parameters from first principles.
In the approach of Ref. [30], it is implicitly assumed that electrons follow the atoms adiabatically, so that, for any atomic configuration, the electrons are in their ground state configuration. Hence, such models allow us to perform Born-Oppenheimer molecular-dynamics simulations, exactly in the same way as most DFT codes do. Nevertheless, if we wanted to monitor the electronic properties during such a simulation, or consider situations in which the electrons are not in their ground state (e.g., excitations, additional carriers, etc.), we obviously need to extend the models. This is precisely the step we have undertaken in this work, explicitly including the effect of the electronic band structure using a tight-binding-like description.
Among the most efficient implementations for large scale atomistic simulations, the SCC-DFTB has got a remarkable success. Several approximations are shared with our new scheme. Among them, the second-order expansion of the DFT total energy with respect a charge density fluctuation (Sec. III B), the use of a localized minimal basis set to expand the electronic wave functions (Sec. III C), or the approximation (Sec. III C 2) and parametrization (Sec. IV) of interaction integrals. Within this framework, SCC-DFTB is comparable in speed with semiempirical methods, roughly 2-3 orders of magnitude faster than standard DFT [32]. Therefore, it has become a successful technique in the study of large organic and biological molecules.
However, there are also some remarkable differences between SCC-DFTB and the present approach. In the SCC-DFTB method, the charge density fluctuations are quantified with respect a reference electron density that is defined as the superposition of neutral atom densities; hence, the reference density is not related to the actual electronic structure of the specific material being simulated. Therefore, the fluctuations include the difference in the charge density between the ground state charge density and the reference one due to the chemical bonding. As a consequence, all the valence band orbitals must be included in the calculation for an accurate enough description of δρ. In contrast, we use a difference reference electron density in our case: We usually take the self-consistent solution for the ground state for a representative atomic configuration of the material of interest. Our deformation charge density thus represents a local deficiency/excess of electrons that occur in excited or perturbed electronic states. In other words, our model focuses on the description of electronic excitations and is not hampered by the need to account for chemical bonding. This opens the door to the treatment of excitons, polarons, transport properties, or the analysis of competing magnetic orders, in a very accurate way. In contrast, our method does not allow us to compute cohesive energies and other basic quantities that are accessible within a SCC-DFTB scheme.
The second difference is in the choice of the minimal basis set. The SCC-DFTB method employs atom-centered localized orbitals [101], typically the product of a Slater orbital for the radial part times a spherical harmonic for the angular part. At least one radial function for each valence shell occupied in the isolated atom must be included. Moreover, the basis functions are nonorthogonal, giving rise to an overlap matrix and to a generalized eigenvalue problem.
In contrast, in our scheme we can select the electronic bands that really play a role in the description of the properties under analysis, and use a basis of Wannier functions coming from unitary transformations of the corresponding manifold. This results in a basis set that is perfectly adapted to the specific material and property under investigation, and which may reduce by orders of magnitude the number of basis functions and dimension of the the associated Hamiltonian matrices. For example, in the case of SrTiO 3 treated above in Sec. VI C, instead of requiring 30 atomic orbitals per formula unit [9 atomic orbitals for Sr and Ti (1 s, 3 p, and 5 d) and 4 for the O (1 s and 3 p)], we can consider only the 3-t 2g orbitals of the conduction band of SrTiO 3 . Since we use orthogonal Wannier functions, the mathematical problem of diagonalizing the Hamiltonian is more amenable and does not require the costly inversion of overlap matrices. Naturally, the reduction of the computational burden translates into an increase of the size of the systems that can be simulated.
The third difference comes from the evaluation of the matrix elements. In the SCC-DFTB method, only two-center Hamiltonian and overlap matrix elements are treated and explicitly evaluated, neglecting several three-center contributions to the corresponding matrix elements. Moreover, (i) the diagonal elements of the Hamiltonian are taken from the eigenvalues of the free atom [12], so the crystal field terms are not considered; and (ii) intra-atomic electron-electron interactions are averaged, without considering differences shell-byshell [32]. In our case we explicitly retain many three-and four-center integrals in the parametrization process, and the diagonal matrix elements are sensitive to electrostatic effects, as explained in Sec. III E. This allows us to treat transition metal systems that traditionally have been challenging for SCC-DFTB approaches [102,103].
The last important difference is the treatment of the lattice and the interatomic interactions. In the SCC-DFTB method, the term that involves the ion-ion repulsion, together with the DFT double counting terms, are included in a short-range repulsive potential. It is approximated as a sum of two-body potentials, fitted from the difference of self-consistent DFT calculations and the corresponding tight-binding band-structure energy for suitable reference systems at various interatomic distances. In the new proposed scheme, we deal with the lattice dynamical part by using a model potential (force field) that is directly fitted to a training set of relevant DFT data, achieving very good agreement with the first-principles energies (typically, accuracies below 1 meV per atom can be achieved for the relevant part of the energy surface). Hence, the generality of the employed model potential, and the flexibility in its definition and truncation, makes it possible to have a very accurate description of the Born-Oppenheimer surface.

VIII. CONCLUDING REMARKS
In this paper we have presented a first-principles-based multiscale method, which we denominate second-principles, that makes it possible to compute the properties of materials at an atomic level, with an accuracy essentially equal to DFT, and at a very reduced computational cost. Our approach is based on dividing the electron density of the system into a reference part, usually corresponding to its neutral ground state at any geometry, and a deformation part, defined as the difference between the actual and reference densities. We take advantage of the fact that the largest part of the system's energy depends on the reference density and can be efficiently and accurately described by a force field with no explicit consideration of the electrons. Then, the effects associated to the difference density can be treated perturbatively with good precision by working in the Wannier function basis corresponding to the reference state. Further, the electronic description can be restricted to the bands of interest, which renders a computationally very efficient scheme.
Conceptually, the present approach constitutes a fresh look at the problem of how to describe lattice and electronic degrees of freedom simultaneously and effectively, introducing a convenient partition of the energy that permits an accurate treatment of both types of variables and their mutual interactions. In our view, our method constitutes a significant step beyond the usual techniques-ranging from molecular-mechanics to tight-binding and quantum-mechanics/molecular-mechanics schemes-towards a more unified model.
As illustrated by the examples described here, the present approach allows us to obtain DFT-like accuracy in the analysis of subtle physical effects, like those determining the relative stability of the magnetic phases of NiO, or those involved in the structural relaxations and screening processes associated to the two-dimensional electron gas formed at the interface of LaAlO 3 and SrTiO 3 . Note that these problems-which involve electron correlation effects, transition-metal ions, etc.-are usually hard to treat within DFTB schemes [103].
As currently formulated, our approach has only one essential limitation: It is restricted to systems in which it is possible to (loosely) define an underlying bonding topology that is to be preserved. Hence, while the method allows the system undergo significant structural modifications, e.g., like those involved in typical ferroelectric or ferroelastic phase transitions, it is not possible to study full-blown bond breaking directly with it. Nevertheless, this limitation can be overcome by using our method in multiscale simulations that permit a more detailed treatment (e.g., with DFT) of the regions of the material in which the constant-topology condition is not satisfied.
It is also important to note that the constant-topology condition is perfectly compatible with the study of many structurally nontrivial cases, such as nanostructured materials, surfaces, chemically-disordered solid solutions, coexistence of different structural and electronic phases, etc. Hence, the application scope of our scheme is enormous.
Let us also note that the present method can be extended to cover physical effects not mentioned here. For example, it is possible to expand it to treat relativistic phenomena (as spinorbit effects) or time-dependent nonequilibrium situations (as resulting from the interaction with light) in essentially the same way as the initial DFT implementations were extended to do so (by implementing a noncolinear treatment of magnetism, time-dependent DFT, etc.). Further, since our method provides us with a Hamiltonian for the interacting system, one could imagine solving the electronic problem in ways that go beyond the mean-field approach adopted here, and thus better account for many body effects.
The ability to simulate systems with thousands of atoms, treating both lattice and electrons accurately, may permit for the first time predictive investigations of a variety of intriguing phenomena-e.g., Mott transitions, coupled spin-lattice dynamics, charged and conducting domain walls, polaron transport, etc.-in realistic conditions of temperature, applied fields, etc. We thus believe that the present method has the potential to significantly advance our understanding of some of today's most interesting problems in condensed-matter physics and material science.