Ab Initio Absorption Spectrum of Nio Combining Molecular Dynamics with the Embedded Cluster Approach in a Discrete Reaction Field

We developed a procedure that combines three complementary computational methodologies to improve the theoretical description of the electronic structure of nickel oxide. The starting point is a Car-Parrinello molecular dynamics simulation to incorporate vibrorotational degrees of freedom into the material model. By means of complete active space self-consistent field second-order perturbation theory (CASPT2) calculations on embedded clusters extracted from the resulting trajectory, we describe localized spectroscopic phenomena on NiO with an efficient treatment of electron correlation. The inclusion of thermal motion into the theoretical description allows us to study electronic transitions that, otherwise, would be dipole forbidden in the ideal structure and results in a natural reproduction of the band broadening. Moreover, we improved the embedded cluster model by incorporating self-consistently at the complete active space self-consistent field (CASSCF) level a discrete (or direct) reaction field (DRF) in the cluster surroundings. The DRF approach offers an efficient treatment of electric response effects of the crystalline embedding to the electronic transitions localized in the cluster. We offer accurate theoretical estimates of the absorption spectrum and the density of states around the Fermi level of NiO, and a comprehensive explanation of the source of the broadening and the relaxation of the charge transfer states due to the adaptation of the environment.


I. INTRODUCTION
Nickel oxide is a prototype transition-metal oxide with a relatively simple crystalline structure but complex electric and magnetic properties.Understanding the electronic structure of this basic material is an important step to achieve a better knowledge about more sophisticated systems such as superconducting oxides and manganites with magnetoresistance.][6] The clear nonconductivity observed in NiO is in contradiction to what can be expected from simple arguments regarding the partially filled 3d band.Various spectroscopical techniques have been employed since the 1950s to study the electronic structure of NiO and understand the nature of its nonconductivity.The most widely used techniques have been based on photoabsorption spectroscopy or the combination of photoemission with inverse-photoemission spectroscopy.These studies invariably showed that there is an energy gap between the insulating electronic ground state of the system and the lowest conducting excited states.This gap is the source of nonconductivity in NiO and is measured to be ∼4 eV.However, the interpretation offered by those studies about the nature of the conductivity gap in NiO is diverse.Each technique measures different electronic processes, which cause the analysis to be dependent on the type of spectroscopy employed, even though they coincide in the measure of the conductivity gap.The difficulty in encountering a satisfactory and univocal explanation of the insulating behavior of NiO is directly related to the complexity of its electronic structure that can be observed from the spectroscopic data available.
Ultraviolet-visible absorption spectroscopy (UAS) offers direct measurements of electronic transition energies.The UAS spectrum of NiO shows a high absorption plateau above ∼4 eV and some features with low optical activity below it. 7he high absorption plateau extends constantly to 8 eV with some signals appearing along it. 8Thus, the optical gap in NiO can be defined as the optically almost inactive region below the absorption edge at ∼4 eV.The signals observed below the optical gap are very similar to the spectrum of Ni:MgO. 9Therefore, these weak absorptions are assigned to dipole-forbidden intra-atomic Ni d-d transitions 7 revealing the high localization of Ni(3d) orbitals on the Ni site. 8The Ni 2+ cation can be considered in an octahedral coordination environment with six O 2− anions where the crystal-field strength (10Dq) splits the Ni(3d) orbitals in e g and t 2g molecular orbitals (MO).The ground state of the Ni 2+ ion in this octahedral environment is a 3 A 2g (t 6 2g e 2 g ) and the lowest d-d excited states have configurations of the form (t 5 2g e 3 g ).The electron-hole pair formed by these d-d transitions is a so-called Frenkel exciton localized on the Ni atom and does not contribute to conductivity.][12][13][14][15] They all confirm the original assignment of Newman and Chrenko, although not all techniques are equally well suited to resolve all the different electronic states of the Ni-3d 8 manifold.
Early interpretations of the nature of the absorption edge at ∼4 eV were more differing than the assignments of the weak features in the optical gap. 7,8,16The resonant Raman scattering experiments by Merlin and collaborators for NiO (Refs.17 and 18) gave an empirical proof pointing toward a ligand-to-metal charge transfer (LMCT), from the ligand O(2p) to metal Ni(3d) orbitals, as the source of the absorption band.However, it is not clear if this charge-transfer process is a localized d 8 → d 9 L transition (here L means a hole in the ligand) or an itinerant intersite d 8 + d 8 → d 8 L + d 9 transition, with no interaction between electron and hole.Evidence of photoconductivity slightly above the optical gap indicates the presence of free carriers in the system as a result of the absorption process. 16,19Since the electrons in the localized Ni(3d) orbitals will be trapped, the only available carriers formed by the LMCT are the holes in the ligand O(2p) orbitals, which by hopping can produce polaronic conductivity.Hence, this observation points toward an intersite LMCT as the source of the high absorption plateau.Nevertheless, the fact that photoconductivity appears slightly above the optical gap indicates that there are some excitonic effects involved in the absorption process that point to a localized intrasite LMCT instead, where the electron-hole pair will be moderately trapped at the site of absorption.Therefore, the conductivity gap of NiO should be slightly larger than the observed optical gap.
Photoemission spectroscopy (PES) measures the binding energy of valence electrons revealing the valence density of states (DOS) in relation to the Fermi level (E f ).The x-ray photoemission spectroscopy (XPS), 20 resonant XPS, 21,22 and angle-resolved XPS (Refs.23 and 24) have been used to study the valence electronic structure of NiO.Furthermore, inverse PES (IPES) measures the binding energy of electrons in the lowest unoccupied levels, allowing the study of conducting bands and complementing PES data.Thus, the PES/IPES combination provides the complete DOS that covers both valence and conducting levels of the system.The XPS combined with bremsstrahlung-isochromat spectroscopy (BIS) shows a gap between the top of the valence band and the bottom of the conducting band of 4.3 eV for NiO. 19It must be noted that the exact value of this conductivity gap is subject to the band position chosen to define the highest valence and lowest conducting levels, obtaining values that range from 3.20 to 5.67 eV. 25 The experimental results show that the lowest conduction band has mainly d 9 character, and the upper part of the valence band has contributions from both d 7 and d 8 L levels.Nowadays, it is widely accepted that NiO is a charge-transfer insulator with a final d 8 L + d 9 state because the d 8 L levels are found at slightly lower binding energies than the d 7 levels.
However, the descriptions of the nature of the conductivity gap present many subtle differences regarding the degree of localization of the transitions involved 25 and the contribution and hybridization of the Ni(3d) and O(2p) levels in the valence band.Moreover, the indirect measure of the electronic levels offered by PES/IPES techniques are inherently delocalized, giving no information on local transitions which could be important to get the complete picture of the electronic structure of NiO.A study with site-specific x-ray photoelectron spectroscopy (SSXPS) localized on Ni and O sites showed that the top of the valence band has a higher O(2p)-Ni(3d) mixing than previously thought 26 and the comparison of these techniques with x-ray absorption spectroscopy (XAS) proved that this kind of data obtained for the valence band is dependent on the methodology applied. 27][30] An ongoing challenge during the last 60 years has been the development of a theoretical model that could offer an accurate description of both the electric and magnetic phenomena observed on NiO.2][33] The source of the conductivity gap in this model is the large electron correlation in NiO, which is caused by the high localization of 3d states on Ni sites and generates a large onsite Coulomb interaction (U ) for the electrons in the 3d orbitals.5][36] U is defined as the energy of the intersite d 8 + d 8 → d 7 + d 9 process, which makes PES/IPES techniques very convenient to measure it.However, as shown by the spectroscopic data, LMCT processes also play an important role in NiO.The incorporation of the charge-transfer energy ( ) into the theoretical model was performed by Zaanen et al., 37 who also defined the ratio U/ as a parameter to classify the source of nonconductivity in different materials.According to this model, NiO is an intermediate insulator with U and on the same order of magnitude.This behavior makes U and difficult to measure experimentally due to the hybridization between Ni(3d) and O(2p) states.Therefore, much theoretical research has been performed in the last decades to obtain accurate estimates of these parameters.
Density functional theory (DFT) calculations with the local-density approximation (LDA) functional applied to a periodic model of the system is known to be a solid approach to the description of Bloch itinerant states in metallic systems.However, the band-gap values obtained with LDA are clearly underestimated for NiO (∼0.5 eV) (Ref.38) and other ionic insulators.The latter is a consequence of the shortcomings of LDA to correctly treat the electron correlation and properly describe the localized Ni(3d) states in the vicinity of E f .Subsequent improvements over the LDA functional offered more accurate methods to describe highly correlated materials.Some examples are the self-interaction correction (SIC-LDA), 39 the dynamical mean field theory taking into account many-body effects, 40 treating local states with Wannier functions, 41 or by means of hybrid DFT functionals. 42Other approaches add explicitly the onsite Coulomb interaction through various implementations of LDA + U (Refs.43 and  44) and, additionally, the electron correlation can also be accounted for through Green's functions (GW). 45,46The DOS of NiO calculated with a recent application of GW@LDA + U (Green's functions applied perturbatively to LDA + U) is in very good agreement with experimental PES/IPES data, obtaining gaps for = 3.75 eV and U = 5.2 eV, 47 which enforces the description of NiO having ≈ U but being a charge-transfer insulator.An important drawback of the LDA + U method is that it requires U as a starting parameter, which can be either estimated from experimental data or calculated specifically with some external method such as the constrained DFT formalism. 48n alternative to the periodic calculations is the embedded cluster approach.The model of the system defined by this approach is composed by few atoms that form the active site of the material (cluster) which is immersed into a finite approximation of its environment (embedding).Hence, the cluster can be described "all-electron" by a wave-functionbased method capable to treat the strong electron correlation found in systems such as NiO, without any a priori assumption about the effect of correlation on the electronic structure.On the other hand, several implementations have been reported in the literature to treat the remainder of the crystal, i.e., the most appropriate way to embed the cluster.Many of these approaches are based on the theory of separability of McWeeny, 49 the subsystem formulation of DFT by Cortona, 50 or the incremental scheme due to Fulde and Stoll. 51Either atomic group of atoms or periodic calculations can be used to construct a potential that represents the part of the crystal outside the cluster.][54][55][56][57][58][59][60][61][62][63][64] Especially relevant for this work are the recently reported studies of the d-d transitions in cuprates combining advanced embedding techniques with high-level ab initio calculations by Hozoi and collaborators. 65,66ere, we follow a well-established approach in which the embedding is formed by a combination of ab initio modelpotentials 67 (AIMP) with a point charge distribution in order to reproduce the long-range electrostatic interactions.Albeit this approach is not able to deal with delocalized states such as the Bloch states of periodic systems, in highly correlated materials it has proven very accurate in the description of the ground and excited states where the localized Ni(3d) states play an important role. 68,69Fujimori et al. 70 gave a reinterpretation of the PES/IPES spectrum based on parametrized configuration interaction (CI) calculations on a NiO 10− 6 cluster pointing to the charge transfer as the source of the conductivity gap in NiO.They also obtained d-d transitions in good agreement with the UAS data but slightly overestimated LMCT energies (7-8 eV).Furthermore, a multireference CI (MR-CI) based method applied to a NiO 10− 6 cluster 71 resulted in better values for the LMCT process (6.2 eV).Ionization and electron affinity calculations were also performed using a NiO 9− 6 and a NiO 11− 6 cluster, respectively, which resulted in estimations of U = 11.6 eV and = 10.3 eV.These results were further improved in the same study by including semiempirical corrections to the correlation energy and bulk polarization effects through the discrete (or direct) reaction field (DRF) method obtaining values in better agreement with the experimental band gap, = 4.4 eV and U = 4.6 eV.The aim of this work is to continue developing the ab initio approach offered by the embedded cluster model and provide a more complete description of the electronic structure of NiO.By combining ab initio molecular dynamics (MD) simulations with highly correlated electronic-structure wavefunction calculations, we obtain not only the transition energies of local excitations, but also their intensities and bandwidths.
We use Car-Parrinello MD to explore the vibrational degrees of freedom of the ground-state potential energy surface of NiO.3][74] The combination of Car-Parrinello MD and CASPT2 calculations has proven very powerful in the description of the spectroscopic properties of the cytosine molecule in water. 75Hence, this work presents the application of an analog approach to crystalline materials.Moreover, the embedded cluster model has been improved in comparison to previous studies by a self-consistent incorporation of the aforementioned DRF method to take into account the extra cluster polarization produced by the electronic excitations.We give estimations of the low-energy region of the UAS spectrum of NiO, including d-d transitions and LMCT processes, and of its DOS nearby E f , from which we extract the U and parameters, and compare the results obtained with a static embedding and the DRF approach.

II. COMPUTATIONAL SCHEME
The comparison of the optical absorption spectra at 77 and 300 K (Ref.7) shows that the position of the features in the optical gap do not change with the temperature.However, the principal goal of our molecular dynamics simulations is not to investigate the temperature dependence of the transition energies, but we aim to investigate to what extent the thermal movement of the nuclei and the resulting deviations from the ideal cubic structure influence the electronic structure and the interpretation of the measured spectra.
The proposed ab initio procedure to study the electronic structure combines Car-Parrinello MD (Ref.76) simulations at a given temperature with the highly efficient treatment of electron correlation and excited states offered by the CASPT2 method. 77,78Car-Parrinello MD explores the classically allowed vibrational levels accessible in the electronic ground state, and CASPT2 calculates the vertical electronic excitations in a limited number of snapshots of the simulation.Hence, the effect of the movement of the ions on the electronic structure can now be incorporated.The extraction of these conformational structures from the MD trajectory is not trivial because it requires the transformation of the employed material model, from the periodic description of the molecular dynamics to the embedded cluster model of the electronicstructure calculations.The extraction step was automated by a tool developed with FORTRAN95 specifically for this purpose.This algorithm selects multiple distorted structures from the simulation using as reference the ideal structure of the cluster and performs the required transformations on the model.Furthermore, it is important to consider extra-cluster polarization effects induced by the electronic excitations considered, which can produce large electronic distortions on the cluster 79 and were shown to have a significant role in NiO (Refs.71 and 80) and other related compounds. 81,82ence, once the clusters are generated, the polarization of the embedding is incorporated in a self-consistent scheme that has been developed integrating the DRF method [83][84][85] with the electronic-structure calculations.

A. Combination of Car-Parrinello MD with CASPT2
The simulation box in the Car-Parrinello MD consisted of a NiO supercell of 15.7372 Bohr side length, which is equivalent to a 2 × 2 × 2 crystallographic unit cell and contains 32 Ni atoms and 32 O atoms.The core electrons are described by Troullier-Martins norm-conserving pseudopotentials 86 and, in the case of Ni, nonlinear core corrections were added.The electronic potential was calculated by means of the Becke-Lee-Yang-Parr (BLYP) density functional 87,88 using a plane-waves basis set with a cutoff of 80 Ry.The system was equilibrated at 300 K and then maintained at a steady temperature with a Nosé-Hoover chain thermostat. 89,90The simulation was performed in time steps of 0.096 fs with a fictitious electron mass μ = 700 a.u.Finally, the conformational structures were extracted from a trajectory of 2 ps.This MD simulation was carried out with the CPMD software version 3.11.1.
The number of extracted conformations varied depending on the type of calculation performed on them.The main factor in choosing the size of the set of conformations is the computational cost of the subsequent electronic-structure calculations to be performed on each structure.The simpler calculations of small clusters in a static embedding were performed over a set of 100 conformations, but the more complex calculations of bigger clusters obliged us to use a smaller set of 32 conformations.To check the sample's representativeness of the simulation, the thermal distribution of Ni-O bond lengths of the two sets is compared in Fig. 1.The close resemblance of the two density curves shows that both sets used in this work can be considered effectively equivalent in terms of simulation sampling.
The embedded cluster model of each conformation was built from a single supercell equivalent to a 10 × 10 × 10 NiO crystallographic unit cell with a total of 8000 atoms, which   92 where the inner point charges have the corresponding formal charge for Ni 2+ and O 2− ions, and the charges on the outside of the cube have fractional charges to ensure the overall charge neutrality.][95][96][97] The basis sets used throughout the electronic-structure calculations of all clusters are of the atomic natural orbitals type including scalar relativistic effects (ANO-RCC).The contracted basis functions for Ni were (5s, 4p, 2d, 1f ) (Ref.98) and for the O atoms were (3s, 2p, 1d). 99These basis sets were chosen as a balance between computational cost and flexibility to describe electron correlation and polarization phenomena in the cluster.The reference wave function of the CASPT2 calculation is constructed following the complete active space self-consistent field (CASSCF) approach, which provides an unbiased treatment of the different electronic configurations that play a role in the total electronic structure of the cluster.The multiconfigurational CASSCF wave function is obtained by distributing a limited number of valence electrons over a set of valence orbitals in all possible ways.In the present case, the so-called complete active space (CAS) contains 10 electrons and 11 orbitals, in short CAS (10,11).After the orbital optimization procedure, the active orbitals can be identified as five Ni-3d orbitals, five more diffuse Ni-d orbitals in the so-called 3d shell, 72 and one O-2p orbital.This active space offers a solid footing to describe the electronic ground state of NiO (Ref.100) and the excited states arising from internal Ni d-d excitations and from the lowest LMCT excitations.The so-obtained zero-order wave function is further improved by the aforementioned CASPT2 approach, a second-order perturbation treatment of the remaining dynamic electron correlation.We performed the CASPT2 calculations with a shift factor of 0.10 a.u., 101 which we previously tested to be the minimum value to eliminate possible intruder states.
The intensity of the considered electronic transitions was computed by state interaction calculations between the ground state and all the excited states at the CASPT2 level, replacing the Hamiltonian by the transition dipole operator. 102The electronic-structure calculations were carried out with the MOLCAS 7.2 package. 103,104The resulting data set of electronic transitions was statistically treated with the R package 105 to generate the corresponding spectra or DOS.For this purpose, we have employed a kernel density estimation with Gaussian kernel functions with a variable bandwidth (bw) depending on the length of the data set.

B. Self-consistent DRF
The electronic excitations such as the local d-d transitions do not change substantially the electronic density on the site of absorption.Therefore, the polarization on the rest of the crystal in response to such excitation is expected to be small.However, this is no longer the case for LMCT states and even less for ionization and electron addition states, which produce a larger electronic change in the atoms directly involved in the excitation and will induce a major response in the surroundings of that site.The adaptation of the extra-cluster region to the changes of electron density in the cluster will be first electronic and later structural, forming a polaron.The static embedding used so far in the model of NiO is not capable to deal with this situation.To overcome this limitation, we apply the DRF method, which offers a rigorous treatment of the electric response of the embedding through the inclusion of a set of dipole polarizabilities and which can be combined with the quantum description of the cluster.
Even though the DRF method was already applied to NiO successfully, 71 we present here an integration of DRF with CASSCF calculations in a self-consistent manner for highly excited states using the embedded cluster model.In the first step of the developed procedure, the electronic density of the cluster is computed in a static embedding and used to generate the corresponding induced extra-cluster polarization through DRF.The polarized embedding is formed by the same AIMPs and point charges of the static embedding, but adding induced electric dipoles on the lattice sites around the cluster.Hence, the response to changes in the cluster density is taken into account for all lattice sites that are in a box of ∼47 Bohr centered around the cluster.The CASSCF step is then rerun with the DRF polarized embedding to produce a new electronic density, and the process is repeated until the energy of the excited state of interest remains constant.However, there are two main issues with this basic design.The first one is that the electronic density used by the DRF method should be the one from the electronic state that causes the perturbation in the system.For example, the inclusion of the relaxation due to the LMCT process requires the calculation of the electronic density of the 11th excited state (the first 10 states are d 8 spin-triplet states).Unfortunately, state-specific CASSCF calculations for such high states do not converge and, thus, we are forced to work with state-average calculations from which the state-specific response can not be extracted directly.Nevertheless, it is possible to get closer to the state-specific situation by changing the weight of the states in the average.We observed that the induced dipoles vary linearly with respect to the weight of the LMCT state in the state-average CASSCF procedure (Fig. 3).This behavior is observed in several conformations and over a wide range of weights.Therefore, under the hypothesis that the linear behavior is maintained along the complete range of LMCT weights, it is possible to obtain the polarization in the embedding caused by the LMCT state through linear extrapolation.
The second issue finds its origin in the fact that the cluster model lacks the translational symmetry found in the crystalline material.The isolation of the active site and its surroundings in an infinite void causes the DRF embedding to be highly polarized even in the ground state.This behavior has no physical basis because, from a global view of the material in its ground state, the electronic density on one site is equivalent to its neighboring sites, and no polarization should occur.In order to study the physically meaningful electronic transition between the ground state and the LMCT state in its polarized embedding, we first generate an "artificial" DRF embedding for the ground state and subtract it from the DRF embedding induced by the LMCT.With this procedure, we obtain a polarized embedding with contributions exclusively from the charge-transfer process.
The last step consists of the state interaction calculation at the CASPT2 level between the ground state in the nonpolarized embedding and the LMCT state in its polarized embedding to obtain the corresponding transition energies and intensities.This step requires us to compute the LMCT state of the system including the DRF embedding previously generated.However, as shown before, it is in general not possible to perform the CASSCF state-specific calculations of these excited states.Therefore, we perform once again a state-average CASSCF calculation including the LMCT state and the ground state, but this time with the polarized embedding by the LMCT.Nevertheless, for one conformation, we were able to converge the wave function for the LMCT state alone.The final result is nearly identical to the one obtained using a state-average reference wave function, less than 0.2 eV difference in the energy of the LMCT state.
In the case of ionized (d 8 L and d 7 ) and electron addition (d 9 ) states, the procedure is exactly the same but without the difficulties inherent to the treatment of the LMCT states.All these states are the lowest states with different spin multiplicities and, hence, they can be calculated separately by state-specific CASSCF calculations.The DRF method requires us to define the atomic polarizability of each atom in the system.We set α O 2− = 1.988Å3 (Ref.106) for the O atoms and α Ni 2+ = 0.653 Å3 for the Ni atoms, based on the total polarizability of NiO (α NiO = 2.641 Å3 ). 8 The DRF calculations were performed with the DRF90 software. 107

Optical absorption
Figure 4 shows the computed UAS spectrum of NiO at 300 K.This spectrum is formed by electronic transitions of 100 conformations modeled with a NiO 6 cluster in a static nonpolarized embedding.For each structure, we performed a state-average CASSCF/CASPT2 calculation of the lowest 12 triplet electronic states and another state-average calculation of the lowest 17 singlet states.Although there is no symmetry in the analyzed distorted structures, we will use the term symbols of the O h symmetry group to facilitate the discussion and the comparison with previous studies.Due to the conformational distortions, the electronic transitions from the triplet 3 A 2g ground state are no longer forbidden by the dipole selection rule and, therefore, it is possible to compute their intensity.Moreover, by including the spin-orbit coupling, we could also take into account spin-forbidden transitions to singlet excited states, which are less intense but nevertheless do significantly contribute to the final spectrum.The combination of all the electronic transitions from the 100 conformations, 14 431 in total, was performed by means of a density estimation using Gaussian kernel functions with a minimal bandwidth of 0.038 eV.From this density estimate, the UAS spectrum bands of NiO appear naturally.
The energy distribution of the calculated triplet and singlet excited states is detailed in Table I.A comparison with the experimental spectrum reveals good agreement of the computed spectrum in the region below 4 eV, where   the overall agreement with previous data is good, the results suggest that some of the assignments in the spectrum should be reconsidered.Most notably are the three features found in the third d-d band (2.5-4.0 eV), where the first signal at 2.75 eV can indeed be attributed to the 3 A 2g → a 1 T 2g excitation, but the other two features at 2.95 and 3.25 eV both belong to transitions to the different nondegenerate components of the b 3 T 1g state.The excitations to the a 1 A 1g and a 1 T 1g states that also appear in this zone are completely masked by the triplet band.Figure 5 decomposes the d-d region of the global UAS NiO spectrum in the contributions of each single band and shows how the transitions to a 1 E g and a 1 T 2g states appear as weak shoulders on the a 3 T 1g and b 3 T 1g bands, respectively.][110] On the other hand, the high absorption plateau extending over 4 eV is clearly overestimated by ∼2 eV.This band is mainly formed by dipole and spin-allowed LMCT processes where an electron is transferred from an O(2p) orbital to a 3d(e g )-like orbital on Ni.Hidden in this band there are also less intense spin-forbidden LMCT and higher d-d transitions.The large overestimation of the LMCT band is produced by the shortcomings of the NiO 6 cluster employed.This small cluster combined with a static nonpolarized embedding is not well suited to properly describe the charge-transfer process, which will polarize the surroundings to some extent.Nevertheless, the relative intensity between the optically weak d-d transitions and the dipole-allowed LMCT processes is correctly reproduced.The increase in intensity among the d-d bands observed in the experimental spectrum, but not found in the computed spectrum, is due to a background signal of raising intensity seen below 4 eV.This signal is caused by impurities in the sample analyzed in the original experiment, 7 but not treated in our model.

Conduction band gap
We applied a similar approach to describe the band gap of NiO.We increased the model size in this case because the involved cationic and anionic states produce large electronic distortions in the system.As a first step, we used a NiO 14 cluster with a static nonpolarized embedding and the smaller set of 32 conformations due to the increased computational cost.For each conformation, we calculated the lowest anionic doublet ( 2 A n+1 ), neutral triplet ( 3 A n 2g ), and cationic quartet ( 4 A n−1 ) states, which are the ground states of their respective electronic configurations. 111Furthermore, these anionic, neutral, and cationic ground states were calculated using two different active spaces: one with Ni(3d,3d ) MOs and another one with Ni(3d,3d ) + O(2p) MOs.The difference in energy between the cationic/anionic states and the neutral triplet state is the n-electron equivalent of the binding energy of the corresponding monoelectronic levels.The smaller Ni(3d,3d ) active space is used to study the Ni(3d n ) states restricting the ionization to localize on the Ni atom and, thus, obtaining the energies of the 3d 7 , 3d 8 , and 3d 9 levels, which are required to calculate the Coulomb interaction U in NiO.On the other hand, the larger Ni(3d,3d )+O(2p) active space allows the ionization to delocalize on the ligand, enabling the description of ligand-hole levels, such as the 3d 8 L, and the calculation of the parameter.We have calculated the U and parameters of NiO by means of the relation given by Zaanen et al. 37 : 2g .Figure 6 shows the DOS nearby E f for NiO at 300 K computed at the CASSCF level.Both 3d 7 and 3d 8 L bands are highly overlapping, but the top 3d 8 L levels appear at slightly lower binding energies.This DOS is qualitatively in good agreement with the generally accepted behavior of NiO, with a conductivity gap defined by the charge-transfer process and a large mixing between the metal-hole and ligand-hole levels.Moreover, the calculation of NiO had to be done with data from only 24 out of the 32 conformations, because only those 24 structures have a 4 A n−1 electronic ground state of 3d 8 L character, with the hole localized in an O atom.The other 8 conformations have the hole mainly localized in the Ni atom independently of the active space used, being more similar to a 3d 7 level.This behavior confirms the large mixing between the 3d 7 and 3d 8 L levels in NiO.However, even though the hole-level character varies between different conformational structures, we have not observed any direct relation with simple distortions of the NiO 6 octahedron or the Madelung potential on the Ni site.
The mean computed values of U and parameters for NiO are U NiO = 13.8 eV and NiO = 13.0 eV.As expected, both parameters are largely overestimated in comparison to data derived from experiment (U ≈ 5 eV and ≈ 4 eV).However, looking at the computed DOS (Fig. 6), the actual picture is somewhat better.The inclusion of dynamical data and the associated phononic broadening reduce the gap by ∼3 eV.Since the transitions involved are intersite processes, it will be more probable that the two active sites are found in different conformational states.This circumstance in the extreme case reduces these parameters to U NiO = 10.9 eV and NiO = 10.4 eV, which is of course a significant improvement, but unfortunately these values are still too large.Nevertheless, the observed degree of broadening points toward a dominant role of the phonons over the hopping phenomena in determining the character of these ionic states, including the ligand-hole states, thus, confirming the local nature of this system.
Moreover, the calculation of the DOS for NiO at a higher level of theory, incorporating dynamical correlation effects through CASPT2, resulted in largely uneven reductions of the binding energy between the two cationic states.The 3d 7 levels shift to ∼2 eV lower energies, while there is practically no change for the 3d 8 L levels, resulting in U NiO = 8.4 eV and NiO = 11.0 eV.The difference in correlation between these two types of states is caused by a deficient description of the ligand-hole levels compared to the metal-hole levels by the employed NiO 14 cluster model because the hole in a O ligand is surrounded by five Ni AIMP, which do not contribute to correlation.
Therefore, we have obtained a good qualitative description of the DOS for NiO, but we also identified two main issues derived from the embedded cluster model employed.The static model employed is not able to correctly describe the large polarization that the addition/subtraction of an electron produces on the system and, on the other hand, the NiO 14 cluster does not offer a good treatment of ligand-hole states.Both effects lead to a deficient description of the states that take part in these phenomena.

Optical absorption
Previous results obtained with a cluster model of the material in a static nonpolarized embedding revealed the problems in treating the more polarizing excited states such as the LMCT or the onsite electron-removal and electron-addition states.These phenomena generate large electronic distortions in the system that go far beyond the active site modeled by the cluster.In order to correct this shortcoming of the model, we included extra-cluster polarization effects by means of the DRF method as described in Sec.II B.
We started with a NiO 6 cluster and computed the induced polarization in the embedding caused by the lowest LMCT process.Figure 7 shows the modulus of the induced dipoles on the lattice sites surrounding the cluster (left) and a qualitative topological representation of the induced polarization in a slab of the embedding (right) for one arbitrary conformation.We observe large polarization effects near the cluster, which gradually decrease and become very small for lattice sites at more than 20 Bohr from the central Ni atom.Although the chosen embedding size demonstrates to be large enough to account for the major part of the bulk polarization, it is preferable to increase the NiO 6 cluster by including the closest 32 oxygen ions with large induced dipoles and treat their polarization in the ab initio wave function.
Notwithstanding that such a NiO 38 cluster gives a quite accurate treatment of the response phenomena, this large cluster combined with the DRF embedding turned out to be computationally too demanding and therefore we decreased the size of the cluster to NiO 14 .This medium-sized cluster model incorporates into the wave function most of the system polarization with a reasonable computational cost, leaving in the embedding only five highly polarized ions (1 Ni 2+ + 4 O 2− ) found at ∼9 Bohr from the central Ni atom.In order to check the consistency of this model, we performed DRF calculations for the local d-d transitions because these kinds of excitations produce very small electronic changes in the crystal and, thus, the corresponding transition energies should not vary significantly by adding the polarization effects that they induce in the embedding.The resulting energies of all the d-d transitions increase by less than 0.03 eV on application of DRF, confirming the validity of the model.
Next, we repeated the calculations of the UAS spectrum of NiO but now with the bigger NiO 14 cluster and incorporating the bulk polarization induced by excited states.We performed two series of calculations on 32 conformations applying the DRF method separately to the lowest nine 3d excited states and to the lowest two LMCT states.For both series, we obtained the average polarization caused by these two types of transitions.In this case, we neglected the spin-orbit coupling and hence only considered spin-allowed transitions.Figure 8 shows the computed spectrum including both the d-d region and the charge-transfer band.The triplet d-d bands remain practically unchanged compared with the nonpolarized static embedding.The mean energy of each band changes by ∼0.1 eV, which is solely caused by the increase of the cluster model.On the other hand, the situation for the LMCT band is completely different.The relaxation produced by the induced bulk polarization shifts the LMCT band to lower energy by ∼2 eV bringing the computed spectrum in good agreement with the experimental one.
A detailed analysis of the origin of the broadening of the LMCT band reveals some interesting aspects.Figure 9 shows a comparison of the charge-transfer band with and without extra-cluster polarization.Once the bulk polarization is conveniently treated, the induced polarization by the LMCT becomes the main factor in defining the broadening of the band.The polarization caused by the LMCT in the outer side of the Ni-O bond subject to the charge-transfer process varies linearly along the band, experiencing a larger relaxation the LMCT states with a larger induced polarization.Since the LMCT excitations are highly localized in NiO, we quantified the degree of distortion of the NiO 6 octahedron along the LMCT band to check if there is a correlation between the induced polarization and the local geometry of the absorption The mean values of σ r and σ ∠ in intervals of 1.0 eV for each LMCT band are shown with overprinted filled circles in Fig. 9.Although we have no sufficiently large samples to go in deeper detail than 1.0 eV, and these are simple measures of the cluster distortion that only account for the immediate environment of the Ni central atom, they do offer a first glimpse to the internal structure of the LMCT band.In the first place, we concentrate on the LMCT band obtained before including the extra-cluster polarization.The center of the band arises mainly from conformations that are weakly distorted with respect to the ideal octahedral structure, with low values of both σ ∠ and σ r , whereas the distortion of the system increases to both extremes of the band.The situation changes significantly once the electronic response of the environment is taken into account.The broadening of the band is now primarily driven by the LMCT induced polarization.In this case, the conformations with a highly distorted NiO 6 cluster, mainly because of a larger σ ∠ , produce a higher polarization in the Ni-O bond direction where the LMCT occurs and experience a larger stabilization of their LMCT transition energies.As a consequence, these states become the front of the band.On the other hand, the more symmetric structures induce a lower polarization on their surroundings having higher LMCT states and appearing at the high-energy end of the band.Moreover, we performed calculations of the LMCT states for an idealized octahedral structure of NiO using the same model employed for the rest of the conformations.The LMCT energy of this ideal conformation of NiO resulted in 8.55 eV with a nonpolarized embedding and 6.83 eV with DRF, which are in good concordance with the analysis previously made.Nonetheless, this is indeed a first simple analysis of the broadening of the LMCT band, which indicates that the crystalline structure has some influence in the stabilization of the LMCT.Further analysis of this influence will require a more detailed study with a larger conformational set and considering long-range structural distortions.
These results demonstrate that the response phenomena have a crucial role in those excited states that cause large electronic distortions on the absorption site.Moreover, the embedded cluster model proves to be adequate to study these phenomena on highly correlated materials such as NiO, where the underlying phenomena of the UAS spectroscopy can be considered highly localized.The theoretical treatment of the electronic structure based on the ideal and highly symmetric geometry usually ascribes band broadening to hopping processes of the holes or electrons.3][114][115] Here, the broadening of the valence and conduction bands is caused by the inhomogeneity of the structure due to the movement of the nuclei, which introduces variations in the relative energies of the relevant n-electron states.Furthermore, the DRF method showed as an adequate approach to treat the polarization in crystalline environments under the condition that the largest polarization effects are treated in the wave function, and that the artificial polarization of the ground state is removed from the system.
It has been suggested that aside from the above-investigated local oxygen-to-nickel excitation as origin of the onset of the optical absorption edge, nickel-to-nickel charge-transfer excitations could also contribute to this onset. 70We have explored this hypothesis by performing some calculations on an embedded cluster with two Ni centers, similar to those used in Ref. 70.The lowest optically allowed excitations have a strongly mixed character making the identification as metalto-metal or ligand-to-metal charge-transfer state meaningless.This was actually also observed in the work of Minami and Fujimori, who described large screening effects of the metal-to-metal excitation by the neighboring oxygens.

Conduction band gap
We applied the DRF method to the calculation of the U and parameters of NiO.We performed this study on only five conformations due to the problems found in the DOS calculation without DRF (Sec.III A), but nonetheless we tried to address the issue regarding the lack of extra-cluster polarization offered by the static embedding.The resulting parameters at the CASPT2 level are U NiO = 3.0 eV and NiO = 4.6 eV, which are in much better agreement with the expected values of U ≈ 5 eV and ≈ 4 eV than the values obtained with the nonpolarized embedding.Hence, the inclusion of electronic response effects turns out to be very important for the description of the band gap, as it is expected due to the large electronic distortions caused by the ionized states involved.
However, these values for U and still suffer from the poor description of the ligand-hole states offered by the NiO 14 cluster.NiO should be lower than U NiO , which is the case at the CASSCF level, but U shifts from 7.2 to 3.0 eV upon application of the CASPT2 step, whereas practically does not change.Figure 10 shows a qualitative representation of the polarization induced by the cationic 3d 7 (left) and 3d 8 L (right) levels.There are large differences in the embedding polarization depending on where the hole is localized.In the case of the ligand-hole level involved in NiO (3d 8 L), there is a large polarization on atoms outside the cluster region which, preferentially, should be treated in the n-electron wave function by enlarging the cluster.
Furthermore, the calculated U value is slightly lower than previous theoretical estimations. 47The source of this deviation could be an overestimation of the extra-cluster polarization by the DRF method, which relies on atomic polarizability parameters.It must be noted that the calculation of the DOS with the DRF embedding is not possible due to a limitation of the employed model.The potential energy between the cationic, anionic, and neutral systems results in large nonphys- ical differences.The NiO 14 cluster with 42 Ni 2+ AIMPs is far from the stoichiometry of the NiO crystal and has a total formal charge of 56+.This positively charged cluster combined with a polarized embedding puts each state in a different potential background, preventing the calculation of the electron affinity and ionization potential and, consequently, the computation of the DOS of NiO with DRF.Notwithstanding this shortcoming of the material model, the calculation of U and parameters is possible because the involved cationic and anionic states produce opposed electric potential shifts, which are mutually compensated on the calculation of the band-gap parameters applying the expression of Zaanen et al. 37 A more detailed explanation of this limitation of the employed cluster model is available in the Appendix.

C. NiO 38
We compared the method proposed here for including the extra-cluster polarization through DRF against calculations of a big NiO 38 cluster in a static nonpolarized embedding.The NiO 38 cluster is expected to include the major part of the response due to the O 2− ions, the most polarizable atoms in the lattice.The static model potential representation of Ni 2+ neglects any polarization effect due to these cations.We have calculated the LMCT absorption edge and the U and parameters with this big cluster (Table II).In the case of the LMCT process, we observe an important lowering of the excitation energy from 8.53 to 6.11 eV.This behavior is expected, as we already showed that the ions with the largest induced dipoles are left outside the cluster in the NiO 14 case, whereas the NiO 38 incorporates them.However, the bulk polarization extends beyond the limits of the NiO 38 cluster (see Fig. 7) and, moreover, the induced dipoles on Ni 2+ ions in the cluster surrounding are significant.Therefore, it is not possible to obtain accurate results of the LMCT energies with a static nonpolarized model.We tried to improve this result by using the DRF embedding obtained with the NiO 14 cluster in the NiO 38 , but even though the resulting LMCT energy is better, it does not offer good accuracy (5.39 eV).Obviously, the induced dipoles and the cluster electronic density must be included in a self-consistent manner, which is computationally too expensive, at present.The situation is similar for the calculation of U : the NiO 38 cluster recovers ∼40% of the polarization caused by the 3d 7 and 3d 9 ionized levels.In the case of cationic/anionic states, the polarization spreads further away and induces weaker dipoles in the vicinity of the cluster than in the LMCT transition, reducing the amount of polarization incorporated in the NiO 38 cluster.On the other hand, in the case of the parameter, there is no gain by using this bigger cluster.An important reason for this situation is that the hole in the 3d 8 L state is not localized on the center of the cluster (see Fig. 10); the localization of the hole in the O atom induces large polarization on its neighboring region, which will need even a bigger cluster than NiO 38 to be treated effectively.
In terms of computational cost, the NiO 14 cluster with the DRF embedding calculations are the most expensive we performed.Due to the iterative process required to reach self-consistency, treating the bulk polarization has a cost ∼5 times larger than using a static nonpolarized embedding.However, the difference becomes rapidly smaller as the cluster size increases.Compared to NiO 38 , the difference is reduced to 1.5 times approximately, but still resulting in a better description of the system.

IV. CONCLUSIONS
We have combined successfully three different methodologies to provide a (nearly) ab initio description of the electronic structure of local excitations in ionic insulators taking NiO as a prototype.We used the Car-Parrinello MD simulations to incorporate the effect of the lattice vibrations on the local structure at one of the Ni sites.By means of the CASPT2 method, we were able to study the electronic excited states with an efficient treatment of the electron correlation.Moreover, we could include electronic response effects in the form of bulk polarization in the material model by coupling the electronic-structure calculations self-consistently with the DRF method.To do so, we assign polarizabilities to the ions taken from the literature.A complete ab initio treatment in the most strict meaning can be obtained by determining these polarizabilities from calculations.
As a first application of the procedure, we have studied the absorption spectrum of NiO and its DOS near the Fermi energy.These features are convenient to be treated by the proposed procedure as they imply localized changes in the electronic structure, which can be described using the embedded cluster model.The best results for both the absorption spectrum and the DOS have been obtained with the NiO 14 cluster incorporating the DRF embedding to take into account extracluster polarization effects.
The computed UAS spectrum is in good agreement with the experimental data (Fig. 8).The results obtained by incorporating the Car-Parrinello MD simulations allowed us to pinpoint the origin of some features observed in the spectrum, such as the d-d transitions to a 1 E g and a 1 T 2g singlet states, at 1.97 and 2.86 eV, respectively (Fig. 5).Moreover, by integrating the DRF into the model, we could explain the origin of the broadening of the LMCT band as caused by a linear decrease of the induced polarization on the different conformational structures that form the band (Fig. 9).Finally, the computed U and parameters are in good agreement with their expected values.Even though the phononic broadening reduces the conductivity gap by ∼3 eV, the main reduction factor with respect to a static embedding appears to be the bulk polarization, which lowers this gap to U NiO = 3.0 eV and NiO = 4.6 eV.
The major drawback of the method is the high computational cost, which makes cumbersome the use of larger conformational sets to analyze with deeper detail the resulting data for both the UAS spectra and the DOS.Furthermore, it has been shown that it might be also necessary to use a larger cluster than NiO 38 to treat ligand-hole states such as the 3d 8 L, where the hole is localized in a ligand O atom away from the center of the cluster.The effect of the extra-cluster polarization in cationic/anionic states is still an issue not completely resolved, which will need further investigation.
Nevertheless, the method presents an important advance in the ab initio description and analysis of the electronic structure of transition-metal oxides.We conclude that the local LMCT process is the origin of the optical gap starting at 4.1 eV and that the conductivity gap is slightly larger, between 4.5-5.0eV, being caused by U or depending on the conformational distortions.The computational issues are expected to become less and less important given the ever increasing computer power.potential upon electron addition will be Analogously, the electric potential at Q caused by the induced isotropic polarization from the ionization (electron removal) of the site (Fig. 12) is always negative and, hence, the total potential will be There is no such a change in potential for the LMCT process (Fig. 13).The highly localized charge-transfer process in NiO generates both a hole and a particle in the cluster, producing a closed electric field around the absorption site that does not affect the electric potential of the system.The induced dipoles parallel [V Q (p)] and perpendicular [V Q⊥ (p)] to the LMCT  direction become mutually opposed and cancel each other: This situation prevents the computation of the DOS including extra-cluster polarization effects by means of the standard clusters used so far.Due to the aforementioned electric potential shifts, the electron affinity and ionization potential energies obtained including a polarizable embedding will be larger, in absolute terms, than the nonpolarized counterpart.However, it is actually possible to calculate the U and gaps including the polarization response of the embedding.Applying the definition of these parameters done by Zaanen et al., 37 the contributions to the potential from the particle and the hole mutually counterbalance: . The calculation of the electron affinity and ionization potential requires a description of the cationic, anionic, and neutral states on the same potential energy background.Since the polarization caused by these three states is completely different, the only option to not drag the electric potential shift to the electronic-structure calculation is to work with a neutral cluster as reference system.The electric potential energy of the central charge in the former simplified bidimensional model is defined as Therefore, a embedded cluster model with a total cluster charge Q = 0 will be capable to describe the different ionic and neutral states of the system in the same potential energy background, regardless of the polarization of the embedding.We built various cluster models to explore and test this possibility.Table III shows the nuclear-external field potential energy of three different cluster models with varying global charge.For all cases, we used the polarization produced by a cationic state.The cluster charge is changed by including additional model potentials or point charges in the quantum region to compensate its charge.The results show how the interaction between the embedding dipoles and the cluster lowers to practically zero energy upon neutralizing the cluster.The lack of a complete suppression of this interaction is due to slight asymmetries in the model used, and the fact that the cluster has a determined volume.This approach opens new possible applications of DRF in future studies of ionic states in crystalline systems.

FIG. 1 .
FIG. 1. Thermal distributions of Ni-O bond lengths resulting from the Car-Parrinello MD simulation at 300 K for the 32 conformations set (gray) and the 100 conformations set (black).The dashed line indicates the Ni-O bond length (3.9343 Bohr) derived from x-ray diffraction (Ref.91).

FIG. 2 .
FIG. 2. (Color online) Structure of the three cluster models of NiO.From left to right they are NiO 6 , NiO 14 , and NiO 38 .The central Ni atom is represented in yellow, O atoms in the cluster are in red, and the smaller transparent gray spheres represent Ni-AIMP.

FIG. 3 .
FIG. 3. Variation of eight induced dipoles of an arbitrarily chosen conformation with respect to the increase of LMCT weight in a state average including the ground state.The dipoles shown are those which suffer the largest change with the increase of LMCT weight.The dashed lines represent the linear adjustment performed on each dipole with the following linear regression coefficients: R 2 • = 0.9991; R 2 = 0.9959; R 2 • = 0.9994; R 2 × = 0.9995; R 2 = 0.9994; R 2 = 0.9969; R 2 = 0.9993; R 2 + = 0.9969.

FIG. 4 .
FIG. 4. (Color online) Absorption spectrum of NiO at 300 K obtained from 100 conformations modeled by NiO 6 clusters with a static embedding (black) and the experimental spectrum (red) (Ref.7).The circles mark single electronic transitions from the ground state.The average energy of the absorption bands is noted on the upper part of the graph and in Table I.Density estimation with Gaussian kernel functions (bw = 0.038 eV).energetics, intensities, and bandwidths of the d-d bands are accurately reproduced.The three spin-allowed transitions from the ground state 3 A 2g to the excited states 3 T 2g , a 3 T 1g , and b 3 T 1g are the origin of the three main bands on this region of the spectrum.The remaining features below 4 eV are spin-forbidden transitions to singlet excited states.Although

FIG. 5 .
FIG. 5. (Color online) Breakdown of the absorption spectrum of NiO into its contributing bands at the d-d region.Density estimation with Gaussian kernel functions (bw = 0.038 eV).

FIG. 6 .
FIG. 6. Density of states of NiO for the states nearby E f at the CASSCF level.Obtained from a set of 32 conformations modeled with NiO 14 clusters with a static nonpolarized embedding.The energy of the 3 A 2g state is set as the reference value.Density estimation with Gaussian kernel functions (bw ≈ 0.3 eV).The mean values of U and are representative for the ideal cubic structure and the smaller values are taken with respect to the onset of the valence and conduction bands.

FIG. 7 .
FIG. 7. (Color online) Induced dipoles by the LMCT in a NiO 6 cluster.Left: Variation of the modulus of the dipoles with the distance to the Ni central atom.The empty region between 0-5 Bohr is the cluster's space.Right: Topological qualitative representation of the electronic polarization in a slab of the embedding (red) and the cluster (gray).

FIG. 8 .
FIG. 8. (Color online) Absorption spectrum of NiO at 300 K computed with 32 conformations modeled by NiO 14 clusters with a polarizable DRF embedding (black) and the experimental spectrum (red) (Ref.7).The circles mark single electronic transitions from the ground state.The average energy of the absorption bands is noted on the upper part of the graph.Density estimation with Gaussian kernel functions (bw = 0.092 eV).site.Therefore, we defined two distortion parameters for the NiO 6 cluster, the standard deviation of the 6 Ni-O bonds (σ r ) and the standard deviation of the 12 O-Ni-O angles (σ ∠ ) with respect to the ideal O h structure of NiO 6 :

FIG. 9 .
FIG. 9. (Color online) Comparison between the LMCT bands with polarization of the embedding (black) and without (gray).Filled circles show the NiO 6 octahedron distortion along the band, σ ∠ in green and σ r in blue.Empty orange circles show the total induced dipole in the five most polarized ions outside the cluster.Density estimation with Gaussian kernel functions (bw = 0.229 eV).

FIG. 10 .
FIG. 10. (Color online) Topological qualitative representation of the polarization in a slab of the embedding of a NiO 14 cluster by the 3d 7 level (left) and the 3d 8 L level (right).

FIG. 12 .
FIG. 12. Schematic representation of the induced polarization caused by a created hole (cationic state).

FIG. 13 .
FIG. 13.Schematic representation of the induced polarization caused by a charge-transfer process.

TABLE I .
8verage CASPT2 energies of the Ni-3d8and the lowest LMCT states in NiO obtained from 100 conformations modeled with a NiO 6 cluster in a static nonpolarized embedding.The results provided are the mean energy (E) and the standard deviation (σ ) of the absorption bands.Experimental data are also shown for comparison.

TABLE II .
Comparison of the LMCT energy and U and parameters obtained with the different models employed in this work.

TABLE III .
Nuclear-external field potential energy for various cluster models in a polarized embedding.The polarization corresponds to a cationic state.