Automated all-functionals infrared and Raman spectra

Infrared and Raman spectroscopies are ubiquitous techniques employed in many experimental laboratories, thanks to their fast and non-destructive nature able to capture materials' features as spectroscopic fingerprints. Nevertheless, these measurements frequently need theoretical support in order to unambiguously decipher and assign complex spectra. Linear-response theory provides an effective way to obtain the higher-order derivatives needed, but its applicability to modern exchange-correlation functionals remains limited. Here, we devise an automated, open-source, user-friendly approach based on ground-state density-functional theory and the electric enthalpy functional to allow seamless calculations of first-principles infrared and Raman spectra. By employing a finite-displacement and finite-field approach, we allow for the use of any functional, as well as an efficient treatment of large low-symmetry structures. Additionally, we propose a simple scheme for efficiently sampling the Brillouin zone with different electric fields. To demonstrate the capabilities of our approach, we provide illustrations using the ferroelectric LiNbO$_3$ crystal as a paradigmatic example. We predict infrared and Raman spectra using various (semi)local, Hubbard corrected, and hybrid functionals. Our results also show how PBE0 and extended Hubbard functionals yield in this case the best match in term of peak positions and intensities, respectively.


INTRODUCTION
Vibrational spectroscopies are one of the most powerful, fast and reliable methods for materials' structure identification and characterization.Among the different techniques probing atomic motions, infrared (IR) absorption and Raman experiments are widespread in most experimental set-ups, thanks to their low cost and precision in identifying microscopic features.Despite these advantages, data interpretation for many cases calls for theoretical support.Density-functional theory (DFT) is most often employed as a reliable framework for accurate predictions of such spectra, providing in principle exact estimates of the secondand third-order derivatives of the energy functional that are needed to predict vibrational spectra.Over the past decade, many techniques have been developed to efficiently obtain non-resonant spectra of insulating crystals, such as linear response theory [1,2] for infrared absorption, second-order response of the density matrix [3] and the 2n + 1 theorem [4] for Raman couplings.While being elegant and attractive, these methods come at the cost of rather complex code implementations, limiting in practice applications to selected semilocal functionals and often only to norm-conserving pseudopotentials.This is especially the case for the third-order derivatives of the total energy, required for the Raman calculations.
One common expedient for the latter is to resort to finite differences of the dielectric tensor upon atomic displacements, computed using density-functional perturbation theory (DFPT) [5,6].However, this approach scales in general as 6N at with the number of atoms N at , a burden in low-symmetry systems, and becomes out of reach for the non-linear optical susceptibility tensor, which has been shown to be fundamental to achieve good agreement both in single crystals [7] and powder spectra [8].An alternative framework for carrying out such derivatives is the finite numerical differentiation of forces and polarization upon the application of a homogeneous electric field [4,[9][10][11][12].By this technique, along with finite displacements [13], one can compute straightforwardly all the quantities needed for IR and Raman spectra, as well as allowing to use any functional (Hubbard corrected, hybrids, ...), and pseudopotential formalism (norm-conserving, ultrasoft, PAW) with little effort.Furthermore, the use of finite fields overcomes the expensive computation and poor scaling of the variational dielectric tensor approach, allowing for the study of complex materials such as amorphous and defective structures [12,14].In fact, within this scheme, one would need to perform typically only 12 ground-state calculations in the presence of a homogeneous electric field (see Eq. 13 and related discussion).While appealing, this approach has seen scarce applicability, possibly due to the many steps and pitfalls related to the underlying algorithm, which can be highly non-trivial even for expert users.
A key solution is to streamline all of these operations through a modern workflow manager, adhering in passing to the FAIR [15] principles for data sharing, while being able of automatically submitting, parsing, and processing the outputs generated by multiple calculations.In this study, we present a comprehensive formulation for a unified finitedisplacement and finite-field approach, along with efficient Brillouin zone sampling, and encode this formalism in a self-contained Python package that operates within the AiiDA framework [16,17], offering a scalable computational infrastructure for automated and reproducible workflows and data provenance.
To showcase the capabilities of the approach, we demonstrate its effectiveness in predicting the infrared and Raman spectra of LiNbO 3 using seven different functionals: LDA, PBE, PBEsol, PBEsol+U, PBEsol+U+V, PBE0, and HSE06.

RESULTS
The finite-displacement and finite-field method Infrared and Raman spectroscopy are investigation methods probing the atomic vibrations of materials.In both techniques, a source of light is shone on the sample, and the incident photons will be either absorbed or scattered by the material.In IR experiments, peaks of the absorption spectra correspond to collective atomic motions in resonance with that frequency, corresponding to long-wavelength polar optical phonons in crystals.Raman spectroscopy records instead a scattered frequency ω S ; in this case, the excited lattice vibration corresponds to the shift from the incoming laser frequency ω L .Far from resonance the two phenomena are described as the change in polarization P and the change in polarizability χ, for IR and Raman respectively, leading to different selection rules for the activated phonon modes which ultimately provide a comprehensive fingerprint of the material.The former mechanism is coupled to the Born effective charge Z * I [1,18] of an atom I in the unit cell, and the latter instead to the Raman tensor defined as dχ/dτ I in the Placzek approximation [19][20][21], where χ is the electronic susceptibility tensor and τ I the atomic displacement.
In terms of these quantities, the scattering intensities of a normal mode ν are for infrared spectra, and for Raman Stokes processes (i.e. when the scattered light frequency is lower than the incident one).The tensors defining these amplitudes are defined via the polarization vector and the Raman susceptibility tensor where ω ν and e ν are the frequency and the eigenvector of the phonon ν, n ν (T ) is the Bose-Einstein occupation function, m I the mass of the atom I, and Ω the volume of the unit cell In the harmonic approximation, the phonon normal modes at zone center are obtained from the interatomic force constants (IFCs) where E is the potential energy surface.The IFCs can be written in terms of first-order derivative of atomic forces F I with respect to atomic displacements, and thus it can be obtained using a set of small displacement configurations, as discussed in the literature [1,13].The diagonalization of the mass-scaled IFCs, describing the harmonic ion dynamics, gives the phonon frequencies and eigenvectors: While the IFCs extend over an arbitrary (but decaying) number of unit cells, these photon spectroscopies utilize photons which carry a negligible momentum q ≈ 0, thus only probing phonons in the vicinity of the zone center, and justifying the evaluation via finite displacements in the primitive cell only.The Born effective charges and Raman tensors can also be defined as derivatives of forces with respect to a macroscopic electric field E [4,12].For the former we have: while for the latter, the Raman tensor reads: In a finite-difference scheme, this means that a small homogeneous electric field must be applied in order to carry out the numerical differentiation.In periodic-boundary conditions, this can be accomplished by extending the energy functional E[ψ] to include a term [9,10,22] accounting for the coupling of the electric field and the total polarization of the system (often referred as electric-enthalpy functional): where P ion and P el are the ionic and electronic polarizations, respectively; the latter being described by the modern theory of polarization [23,24].Such functional, while not bounded from below, admits local minima representing self-consistent stationary solutions [25,26], if E is properly chosen; this guarantees the simultaneous evaluation of forces and electronic polarization at convergence.As a consequence, the electronic polarization (that comes for free) can be exploited to extract its first and second-order derivatives with respect to an electric field [4,10].These correspond respectively to the high-frequency dielectric tensor ϵ ∞ , used to account for the non-analytical contribution to the IFCs [2], and the non-linear optical susceptibility χ (2) , key in non-centrosymmetric crystals where the Frölich contribution to the Raman tensor could be sizable [7,8].Their relationship with the electronic polarization derivatives is as follows:

Numerical differentiation
The possibility of consistently describing homogenoues electric fields allows to access any derivative in atomic displacements and electric field, while being easily extensible to advanced functionals and pseudopotential formulations [12].We calculate here the tensors for infrared and Raman cross-sections using a central derivative formulation.These derivatives are performed around E = 0, which resembles the common experimental setup where an external electromagnetic source is negligible, either natural or artificial.Indicating (by dropping the indices) a computed quantity as A, either a force or the electronic polarization, we calculate its m-th derivatives with a numerical accuracy of n-th order via the discretization formula: where ∆E i > 0 is the small electric field, which we refer to as the finite difference step, along the i-th Cartesian direction, and c m l the central derivative coefficients that can be automatically found using the Fornberg's algorithm [27].The discretization is performed using uniform steps such that l = − n 2 , . . ., 0, . . ., n 2 is always an integer number and the applied fields are evenly spaced (see Fig. 1(a)).For mixed second-order derivatives of A, i.e.
for the Raman and non-linear susceptibility tensors, we exploit a formula devised in [12] to reduce the number of calculations.We assign E i = E j = λ and express the mixed derivative as: where each second-order derivative is evaluated through Eq. 12; this can be visualized schematically in Fig. 1(a).For example using a 2 nd order numerical accuracy, we would need only 12 self-consistent field (SCF) calculations with a non-zero electric (∆E-SCF) field to carry out all the needed tensors both for IR and Raman spectra, even in low-symmetry systems such as amorphous materials [12,14]; this number can be reduced if one exploits the symmetries of the material.In fact, there could exist symmetries belonging to the point group of the crystal that transform an electric field direction into an other.For example, if the crystal satisfy inversion symmetry, this is sufficient to reduce by a half the number of ∆E-SCF calculations needed, since forces and electronic polarization produced by ∆E i are connected by a similarity transformation to the one produced by −∆E i ; keeping the previous numerical example, we would need then just 6 ∆E-SCFs.For highly-symmetric crystals, such as cubic Si, only two electric field directions must be evaluated; in this case, the number of calculations reduces further to only 2. In general, the number of ∆E-SCF calculations needed scales with the numerical accuracy n of the central formula as where While this approach is advantageous, we also highlight that the choice of the field magnitude is restricted.As already stated before the electric enthalpy functional 9 has no global minima, and long-lived meta-stable states exist only below a certain critical field E c , where the Zener tunneling is suppressed [9,22].An estimate of such critical value can be performed using [9]: where e is the electric charge, a i the lattice vectors defining the unitcell, E gap the electronic band gap, and N i the number of k-points sampling the Brillouin zone in the i-th crystal direction.For very small band gap semiconductors this critical field can be very small, limiting in practice the convergence due to numerical noise.In these circumstances, a balance must be found; nevertheless, for the majority of cases this does not represent an actual limitation.

Directional sampling
A subtle issue related to the Berry phase formalism used in the electric-enthalpy formulation is the dense sampling of the Brillouin zone required to achieve well converged values of the polarization, needed to carry out the numerical derivatives.In fact, one can write the electronic polarization using the string-averaged discretized Berry phase formulation [23] along the direction of a primitive reciprocal lattice vector b i as: Here S nm (k, k ′ ) = ⟨u nk |u mk ′ ⟩ represents the overlap matrix between Bloch's states u lk belonging to occupied bands l, and ⊥ is the number of strings (i.e. the number of k-points orthogonal to b i ), each having N (i) ∥ k-points.Thus, when we apply a static electric field along lattice vectors, an electronic polarization is induced according to the field direction, properly sampled by increasing the number of k-points within each string N (i) ∥ .While the derivatives can be performed in the crystal reference system [28], we describe in the following some possible implementations using a Cartesian reference system (the one employed in the code).In this coordinates system one can follow two approaches, namely either performing an orthogonal transformation of the lattice to have the chosen field direction along one of the new basis vectors, or generalise Eq. 15 to a generic reciprocal vector.It is clear it is more convenient to choose the latter, as it requires only to modify Eq. 15 to account for different ⊥ , while the former would require in general large supercells to accomodate the rotated crystal.Nevertheless, many codes use automatically generated Monkhorst-Pack (MP) meshes, and it is not intuitive how to properly choose a regular grid that samples correctly an arbitrary direction.To automate the choice of such grids, we propose here the concept of directional sampling targeting the increase of the number of k-points in any desired direction and crystal lattice, while still retaining a regular MP mesh.
The main idea is to define the BZ mesh by two quantities: an orthogonal and a parallel k-point distance, d ⊥ and d ∥ respectively.The orthogonal distance represents a minimum distance between k-points in the BZ that are commonly used to achieve convergence of total energies and forces, thus usually producing rather sparse meshes -this quantity is then used as a reference minimum for the generation of the grid.The parallel distance instead is used to define a denser number of k-points in the desired direction, which in our case is referring to the applied electric field.To do so, we transform the electric field direction in crystal coordinates.This step allows us to assign weights for the sampling in reciprocal space.Let Ẽ = ( Ẽ1 , Ẽ2 , Ẽ3 ) be the electric field in the lattice reference system, we then define the weights as: aiming at maximizing the sampling for lattice directions where the electric dipole is expected to be more relevant.Finally, we can assign the number of k-points N i of an MP mesh as: This methodology is easy to implement and generates meshes that better sample the direction where the induced polarization is expected, significantly reducing the number of k-points needed instead when using a constant distance to define the entire MP grid.
We point out that the parallel distance defined may not resemble the true distance among k-points within strings along a defined direction.This is easily understood with a trivial example reported in Fig. 1 substituted with generalised MP grids [29] that can be ad hoc fine tuned for effective and inexpensive samplings.However, this is far from being a trivial task and thus out of scope for the present work.The black circles on top of boxes mean that the calculation starts using charge density and wavefunctions from a previous calculation.

Computational workflow
In this section, we finally present the computational workflows that allow to simulate the desired spectra, which can be infrared, Raman, or both, in an automated fashion.We designed two main independent workflows as AiiDA WorkChains [16]: one carrying out solely the finite displacement calculations for the phonon modes to evaluate the IFCs, which we called PhononWorkChain, and the other for the mixed total-energy derivatives via the homogeneous electric fields, the DielectricWorkChain.When joined together, the above two provide the IRamanSpectraWorkChain, as schematically outlined in Fig. 2.
The IRamanSpectraWorkChain takes as inputs a structure, a code and other details on how to run the workflows (e.g.SCF parameters, wall-time, parallelization options, and so on).These information are then split and/or shared between the two sub-workflows.Both start running an initial SCF ground-state calculation, meant to produce the charge density ρ 0 and wavefunctions {ψ i } that will be used as a starting point for the next runs.We generally represent this restart throughout the flowchart of Fig. 2 by a curved white arrow in a black circle, and we use the dashed arrows to indicate that the files have been copied over to a fresh new folder.When these are finished, the PhononWorkChain generates the structures with irreducible atomic displacements and computes their forces.The latter are then used to produce the matrix of IFCs, which will be employed in a post-processing step to obtain the phonon frequencies and eigenvectors.In parallel, the DielectricWorkChain evaluates E c through a non-SCF calculation and extracts the numerical accuracy n and the electric field step ∆E, as follows: when the critical value is larger than 10 −4 (Ry)a.u.
(Rydberg atomic units; 1 (Ry)a.u.≈ 36.3609V/ Å), a 4 th order numerical expansion is chosen, otherwise a 2 nd order is used.The 4 th order allows to achieve a greater numerical accuracy for the finite differences, thus it is utilized whenever it is possible to apply an electric field above threshold.In cases the electric field that can be applied is smaller than 10 −4 (Ry)a.u., a 2 nd order expansion is sufficient, thus avoiding calculations that would be affected by too much numerical noise [21]; the amplitude of the electric field step in fact may result too small to produce significant digit changes in forces and polarization, hence possibly worsening the numerical derivative even if adding extra points to the formula.The step ∆E is instead selected according to: where the electric field units are again expressed in (Ry)a.u. and round(x,y) rounds the number x to the y digit after the decimal point.These choices should guarantee larger values of the finite step to reduce the numerical noise on forces and polarization, and at the same time higher accuracy to remove the step size dependence of the finite difference formula, which can be verified a posteriori when n > 2 (see also Fig. 1(a)).We point out here that eventually the user can still enforce the accuracy n and/or the ∆E value that will be used in Eqs. 12 and 18.Once an electric field is chosen and the ground-state SCF calculation is finished, a symmetry analysis is performed to find the inequivalent directions Ê of the electric fields to apply, where Ê is a versor.For each such inequivalent direction an MP mesh is determined through Eq. 17 and a series of SCF calculations are then launched at E = 0 with the corresponding mesh restarting from the ground-state density of the previous SCF.The reason for computing the SCF at E = 0 with these meshes is twofold: first, it is crucial to evaluate the polarization operator, set up in following finite electric-field calculations at the corresponding mesh reference, and, second, it guarantees the same level of convergence on forces and polarization for the finite differentiation.In fact, Ê and − Ê produce the same k-point meshes (as can be verified via Eqs.[16][17], hence all the calculated values, at E = 0 included, employed in the finite formula 12 have the same level of convergence with respect to the particular k-point grid.Once finished, electric-enthalpy calculations start for each inequivalent direction Ê at the smallest step value ∆E = ∆E • Ê.If the numerical accuracy n is greater than 2, further SCF calculations with homogeneous electric fields are run, restarting from the previous steps till the maximum value (n/2)∆E; this loop corresponds to increasing/decresaing the value of l in Eq. 12.At this point, the computed forces and polarizations are gathered together to perform the numerical differentiation according to Eqs. 12 and 13 for achieving second and third rank tensors 7-8-10.As a final step, the IRamanSpectraWorkChain joins all the information in a unique AiiDA data type, named VibrationalData, able to efficiently store large arrays and information in the AiiDA repository and database, and with provenance graph for full reproducibility [16].This data type is a self-contained Python class, and it is provided with a number of post-processing methods and features that allow the user to compute powder and single crystal polarized vibrational intensities and frequencies, along with the symmetry labels of the respective phonon modes.We also implemented an IntensitiesAverageWorkChain which performs the spherical average spectra, important in Raman for non-centrosymmetric crystals [8].
As mentioned earlier, the inputs of these workflows can be fully customized by the user.
We do not report here the full list of these inputs, which can be explored through the usual command line interface of AiiDA [16] and through the on-line documentation made available.However, we would like to emphasize the possibility of fully defining the entire set of inputs via a limited set of them, as already been done in other AiiDA packages [17].
In similar fashion to Ref. [17], we identified the fundamental inputs to be the code, i.e.
the Quantum ESPRESSO [30][31][32] pw.x binary, associated with the computer, that runs the calculations, the relaxed structure and the protocol, a string defining in a general and user-friendly fashion the accuracy of the simulation.This minimal set, used through the get builder from protocol method [17] designed for all the workchains, provides the user with a pre-filled builder (the main unit in AiiDA used to run workflows) ready to be submitted.This is adequate for non experts in the field that might be interested in simulating their structure's spectra with minimal knowledge, but also in the context of high-throughput searches or for teaching or testing.
All-functionals LiNbO 3 spectra We demonstrate the power of this approach studying the ferroelectric phase of LiNbO 3 (space group R3c, no.161) for which the vibrational spectra have been long debated in the literature [33][34][35][36] and exhibit large LO-TO phonon splitting at Γ [33] that are ideal for exploring different functional flavours.We selected three commonly used local and semilocal funcationals, LDA, PBE and PBEsol, then Hubbard-corrected PBEsol in the standard and extended version, namely DFT+U [37,38] and DFT+U+V [39][40][41], and two of the most used hybrids functionals for inorganic crystalline materials, HSE06 [42,43] and PBE0 [44].
Through the use of the IRamanSpectraWorkChain, we carried out the simulation of the vibrational spectra for all the functionals above and we compare our results to experimental measurements from the most recent literature [8,35,36].In particular, we compare our simulations in Fig. 3 to (a) different Raman single-crystal polarization setups [35] and (b) two polarized infrared absorption spectra [36], and in Fig. 4 to polycristalline Raman powder spectra [8] for which we use the spherical average procedure accounting for the χ (2)   contribution, as outlined in Ref. [8].The single-crystal data allowed us to compare directly the frequencies and, in the case of Raman, the intensities.In intensities for Fig. 3(a).These spectra setups allowed to accurately estimate the error of the intensities evaluated through the minimization of the relative error: where f represents the unknown scaling factor between theoretical I theo and experimental FIG. 4. Raman powder spectra using the spherical average formula [8] and accounting for the Frölich term to reproduce the characteristic asymmetric line shapes.Theoretical intensities are smeared with a constant 14 cm −1 broadened Lorentzian.The theoretical spectra are reported in order of decreasing error |∆ω|, and Hubbard-corrected PBEsol are abbreviated only with their corrective parameters, as in Fig. 3  I exp intensities, due to the use of arbitrary units in the experimental spectra.The smallest absolute error on frequencies |∆ω| is provided by the PBE0 functional with 8.9 cm −1 , and on intensities |∆I| by the PBEsol+U sc +V sc functional with 15%.It is to note that the latter also shows a remarkable |∆ω| of 9.7 cm −1 , less than 1 cm −1 apart from the best performing functional for this material, whereas the former has 25% error for |∆I|, 10% more with respect to PBEsol+U sc +V sc .We also reported the renormalized absolute error |∆ω| * , for which the frequencies are shifted by the bare mean error ∆ω before the modulus average; this gives some insights whether the discrepancy arises from a systematic error, e.g.due to the geometry used.Interestingly, after renormalization, most of the functionals get rather close error values.While this might be due to the renormalization of the phonon frequencies due to the thermal expansion which tends to soften the vibrational modes, all the theoretical spectra were performed at constrained experimental lattice (see the Method section), meaning the systematic error is unrelated to this effect.While other thermal contributions should in principle be accounted for, their inclusion would require more sophisticated level of theory [47][48][49] which goes beyond the scope of the present work.Moreover, at the experimental conditions of our references, such contributions are expected to produce only minor phonon shifts.Given these points we might conclude in this case that such errors are intrinsic to the functionals.

DISCUSSION
In this paper, we have presented a comprehensive and efficient approach that utilizes finite displacements and finite fields to calculate infrared and Raman spectra with any modern exchange-correlation functional and pseudopotential formulation (norm-conserving, ultrasoft, PAW).Our approach overcomes the current expensive alternative of the variational dielectric tensor technique [6] for Raman coupling calculations, which remains impractical in large or low symmetry systems.This is also very relevant when generating coupling tensors in low-symmetric supercells to incorporate temperature and quantum effects [48][49][50], or to train a tensorial machine-learning potential useful for molecular dynamics simulations.
The formalism has been implemented as a highly optimized, automatic workflow package for the prediction of vibrational spectra.The technical and practical challenges have been successfully addressed with the help of the AiiDA infrastructure [16], which empowered the design of fully reproducible, reusable, user-friendly workflows for one of the most important class of spectroscopic techniques for materials characterisation.As such, we believe it will be broadly useful to both the computational community and the experimental community at large.Furthermore, we have shown that the workflows can effectively operate with state-of-theart functionals and can be employed for a comprehensive analysis of both single crystals and powder spectra.In conclusion, we believe that the present approach will pave the way for accelerated materials characterization of materials displaying complex, challenging chemistry.

METHODS
The infrared and Raman calculations were carried out using the workflows developed, relying on aiida-phonopy [51] for pre-and post-processing and aiida-quantumespresso [16,17], having Quantum ESPRESSO [30][31][32] as DFT quantum engine.Before performing the vibrational spectra calculations, the positions of the atoms in the LiNbO 3 primitive cell are optimized at fixed experimental lattice geometry [52] for all functionals until the total energy and forces acting on atoms are less than 10 −4 Ry/atom and 5 • 10 −5 Ry/Bohr, respectively.
Norm-conserving pseudopotentials from the Pesudo-Dojo library [53] are used for the LDA and hybrid functionals, employing PBE pseudopotentials for the latter, and pseudopotentials from the precision SSSP (version 1.1) [53][54][55][56][57][58][59] library for the others.A cutoff of 90 Ry is used for the wavefunctions; 4 and 12 times greater cutoffs are used for the charge density (respectively for Pseudo-Dojo and SSSP based pseudopotentials).The Brillouin zone is sampled uniformly with a 0.3 Å−1 k-points distance.The Hubbard parameters are computed following the self-consistent procedure exploiting density-functional perturbation theory, as explained in [45,46,60], using a q-point distance of 0.4 Å−1 , orthogonalized atomic orbitals [46,61] for the occupation matrices, and with onsite (intersite) parameters converged within 0.01 (0.005) eV (a very conservative threshold).The hybrid calculations were performed using a 90 Ry cutoff and a single q-point mesh for the exact exchange operator, and the adaptively compressed exchange (ACE) technique [62] was used to speed up the convergence of the calculations.Finally, finite differences were carried out using a 0.01 Å displacement distance and a denser k-point distance of 0.15 Å−1 for the phonon calculations (in the

FIG. 1 .
FIG. 1.(a) Schematic illustration of the central discretization formula 12 and 13 for derivatives with respect to the electric field in the directions i = {x, y, λ = (x, y)}.The black dots represent the applied electric fields l∆E i used in the SCF calculations to evaluate the electric-field-dependent forces and polarization to employ in 12 and 13, while the different coloured circles are associated to the numerical accuracy n achieved in the numerical differentiation (here shown up to n = 6).It can be noticed that the evaluation of the numerical derivative at n > 2 provides automatically derivatives at lower accuracy n; this is helpful to check the convergence a posteriori in respect to the finite step size ∆E.(b) Square reciprocal lattice showing a shifted regular k-points sampling (black dots) of the BZ generated via Eqs.16-17 for the shown electric field Ẽ in reciprocal coordinates.The strings of k-points are indicated by straight black lines, and are only shown on part of the BZ for clarity.In this particular case Ẽ1 = Ẽ2 , hence w 1 = w 2 = 1 from Eq. 16, meaning the generated mesh via Eq.17 have k-points along b 1 and b 2 equally spaced by d ∥ , indicated on the figure by an horizontal arrow.It can be seen that this introduced distance differs from the actual distance d true ∥

− 1 )FIG. 3 .
FIG. 3. Raman and IR spectra of ferroelectric LiNbO 3 as obtained with seven different functional, and compared to experimental measurements.(a) Raman single-crystals spectra for different polarization setups (reported in paranthesis) for pure TO modes corresponding to E and A 1 symmetries.(b) Infrared absorption using polarized light for E (upper panel) and A 1 (lower panel) symmetry modes.The theoretical spectra are ordered according to the error |∆ω| (shown in Tab.I), from smallest (pink, PBEsol+U sc +V sc ) to largest (blue, LDA).Hubbard-corrected PBEsol+U sc and PBEsol+U sc + V sc are labeled just U sc and U sc + V sc , the subscript sc refers to the self-consistent calculation of U and V [45, 46].Theoretical intensities are smeared with a constant 8 cm −1 broadened Lorentzian.

cm − 1
for frequencies ω and in percentage for the intensities I.The estimators are: |∆ω| = N ν=1 |∆ω ν |/N is the mean absolute error; ∆ω = N ν=1 ∆ω ν /N is the mean error; |∆ω max | is the absolute max difference; |∆ω| * = N ν=1 |∆ω ν − ∆ω|/N is the absolute error for frequencies shifted by their mean difference; |∆I| = 100 • min f {∆I(f )} represents the percentage error of the optimal scaled intensities; N is the number of considered modes and ∆ω ν is the difference between experimental and theoretical functional mode frequency.

TABLE I .
.LDA PBEsol+U sc PBEsol PBE PBEsol+U sc +V sc HSE06 PBE0 Statistical error estimators between functional predictions and experiments, given in