Computational Procedure to an accurate DFT simulation to solid state systems

Eduardo O. Gomes, Guilherme S. L. Fabris, Mateus M. Ferrer, Fabiana V. Motta, Mauricio R. D. Bomio, Juan Andres, Elson Longo and Julio R. Sambrano LSQM – Laboratory of Chemical Synthesis of Materials, Federal University of Rio Grande do Norte, Natal, RN, Brazil. Modeling and Molecular Simulation Group CDMF, São Paulo State University, 17033-360, Bauru, SP, Brazil. Department of Physics, Federal University of Pelotas, 96010-610, Pelotas, RS, Brazil Department of Analytical and Physical Chemistry, University Jaume I (UJI), Castelló 12071, Spain. Chemistry Institute CDMF, Federal University of São Carlos, P.O. Box 14801-907, São Carlos, SP, Brazil.


INTRODUCTION
In theoretical approaches, the elaboration of computational models that represent the reality is still a challenge, especially for material science. 1 However, the search for models with an accurate approximation of structural, vibrational, and electronic characteristics is very useful for trends, elucidations, and even in findings of new observable characteristics and properties. 2,3 The solid-state simulation still presents a set of hindrances for the correct description of the electronic structure by purely ab initio methods. As much as the approximation is more and more accurate, the calculations strategies constantly confront the results with experimental parameters, such as band gap energy, vibrational modes, and other target properties, in order to orientate the adopted parameters of the theoretical methods. The constant existence of distortions and stacking faults in real crystals can also generate structural characteristics very different from the expected theoretical global minimum.
Nowadays, DFT methods are the first choice for the computational treatment of solids, providing results comparable to multiconfigurational methods with a considerably low computational cost.
Usually, the choice of the density functional, that includes the exchange-correlation functional, and the use of an appropriate atomic basis set is the first procedure in DFT simulations. An appropriate choice can lead to accurate target properties being simulated reliably, especially when dealing with models for which there are no experimental information for comparison. 9 In particular, hybrid functional has been employed with success to simulate solid-state systems, such as the study of the bulk, 10,11 surfaces, 12 nanotubes 13 , adsorption, and doping processes. 14 As a popular example, the B3LYP 15-17 hybrid functional is a common choice to describe the properties of several systems.
As we know, the exact exchange-correlation functional is unknown; therefore, a considerable amount of research has been carried out in the past decades to find accurate approximations for the exact functional. There is no unique method to find these approximations, and, among these, the introduction and adjustments of empirical parameters are the most used methodology.
According to Bredow and Gerson, 18 the optimal percentage of each functional term can change depending on the system and may also be the target of study for a better approximation of the results.
When considering the atomic basis set, there are several available bases of each atom described in the literature that can be used for molecular or crystalline systems in order to expedite the study of materials. However, it is important to highlight that the development of a specific basis set for a molecule cannot always be successfully applied to solid-state systems, or a basis generated for a specific solid-state system can be appropriately used for others. Mostly, it is necessary to generate a new basis or adjust a predefined basis set, such as the optimization of the coefficients and/or exponents of the most external or more specific shells.
Based on that the above, the aim of this study was to investigate two different strategies to conceive computational models, taking as an example the particular case of the metallic molybdate, BaMoO4 (BMO). Both strategies consider the optimization of the atomic basis set and new parameterization of the mixing fraction of the HF exchange functional. Structural, electronic, and vibrational property data obtained by means of both strategies were compared with experimental and theoretical data reported in the literature to compare and evaluate their accuracy

PROBE SYSTEM
To verify the efficiency of the proposed strategy, it was chosen a system belonging to an important family of inorganic materials (AMoO4) that have a high potential for application in a great diversity of technological devices. [19][20][21] Under typical conditions, when solids of this class are composed of relatively large bivalent cations (A = Ba, Ca, Sr, and Pb), they assume a scheelite-type tetragonal structure. 22 In particular, (BMO) is cited as one of the most important members of this family, mainly because of its photoluminescent properties. [23][24][25] However, a considerable number of studies have reported several potential applications of BMO, such as its uses as anodes in lithium and sodium batteries, 26 catalyzers, 27 component in white LEDS 28 , and thermoelectric devices. 29 This versatility is the result of the structural and shape modifications of its nanoparticles, 30,31 as a function of the synthesis methodology, such as sonochemical, 32 electrochemical, 33 microwaveassisted citrate complex, 34 and microwave-assisted hydrothermal. 27 Looking for theoretical studies on BMO, it was noticed that there have been few reported in this field. Among these, Zhao et al. applied the DFT and fully relativistic self-consistent Dirac-Slater theory to explore the behavior of interstitial oxygen around the Mo 6+ ions on the structure. They concluded that there was an interstitial oxygen combined with the nearest-neighbor lattice oxygen ions to form molecular ions 3 4− and 4 4− . 35 Qin et al. used DFT combined with PBE functional to give support to experimental results about the dielectric improvement as a function of pressure. 36 Sczancoski and coworkers reported a combined experimental and DFT study with B3LYP hybrid functional, to explain the optical properties of BaMoO4 obtained by a coprecipitation method in the presence of polyethylene glycol. 37

COMPUTATIONAL SETUP AND MODELS SYSTEM
The structural and electronic properties of BMO bulk were simulated by means of periodic DFT using the CRYSTAL17 package. 38 The BMO has a tetragonal unit cell that belongs to the space group I41/a (symmetry C 6 4h.  Most of the available atomic basis sets are pre-optimized for a specific system that includes different structures and atomic interactions. This means an optimization process may be necessary to adjust, at the same time, all selected atomic bases to a target system. In this work, Powell's algorithm method 47 , combined with the CRYSTAL program, was used to perform the basis set optimization of the exponent of the valence shell functions to an energy convergence of 10 −6 Hartree. The level of accuracy on the convergence criteria for bi-electronic integrals was controlled by a set of five thresholds (10 -7 , 10 -7 , 10 -7 , 10 -7 , and 10 -14 ). These parameters represent the overlap and penetration for Coulomb integrals, the overlap for HF exchange integrals, and the pseudooverlap, respectively. The shrinking factor (Pack-Monkhorst and Gilat net) was set to 6, corresponding to 36-k-points in the irreducible Brillouin zone.
All stationary points were characterized as a minimum by diagonalizing the Hessian matrix with respect to atomic coordinates and unit cell parameters. The vibrational modes at the Γ point were evaluated using the numerical second derivatives of the total energies estimated with the coupled perturbed HF/Kohn-Sham algorithm. 48,49 The convergence was checked on gradient components and nuclear displacements with tolerances on their root mean square set to 0.0001 and 0.0004 a.u., respectively. The band structure and density of states (DOS) were analyzed with the same k-point sampling employed for the diagonalization of the Fock matrix in the optimization process.

Choice of functional, basis set, and strategies of calculations
As a starting point, the optimization of the structural parameters and band gap energy calculations was performed using the selected atomic basis and different hybrid functionals, as described in Table S1 (supplementary material) and shown in Figure 2. and band gap energy with the best average accuracy to the reference experimental system. The HSE06, B1WC, and B3LYP functionals also showed good and acceptable structural parameters but with discrepant band gap energy.
Nevertheless, the B3LYP functional was also used to perform the simulations for the sake of comparison with the WC1LYP, because it is one of the most used functional to study crystalline materials. 50 Although the results have been somewhat favorable, they can still be improved. The second strategy was the inverse of the first procedure. Structural optimization was performed in which the influence of the HF% variation was analyzed, and, subsequently, the atomic basis set was optimized (same most external orbitals of each atomic center) with the new hybrid functional that considered the new HF%.
The preliminary structural and electronic results of both paths are described below.

HF%)
The valence shell values of the optimized basis with standard WC1LYP and B3LYP functional can be found in Table 1. The index 1 and 2 refer to the position of coefficient in the respective basis set.
The calculated lattice parameters and band gap energy after this procedure were a = 5.62  Table S1). However, there were improvements in band gap energy on both functionals, with a correction of 3%. An interesting point to highlight in Table 1 is that the coefficients 1 αsp of Ba and ²αsp of Mo atoms are the ones that had the greatest value variation.
These new outputs were then used to elaborate the starting point of the gradual modification of the HF% of the WC1LYP and B3LYP functionals. The structural parameters and band gap energy of each HF% modulation are shown in Figure 3 and Table S2. These results indicated a better description with 12% of HF contribution, in relation to the experimental characteristic previously described.  Figure 4 and Table S3 depict the structural parameters and band gap values of all tested HF% using the original atomic basis. The values depicted in Figure 4 show that the 12% percentage had better proximity to the experimental results and, therefore, was selected to be used in the atomic basis optimization.

Second path: HF% analysis followed by basis set optimization (HF%/BasisOpt)
The valence shell values of the optimized atomic basis with modified WC1LYP(12%) and B3LYP(12%) can be found in Table 2. The index 1

Electronic properties
Each one of the final four models was well explored to understand the consequences of the calculation route choices. A better description of electronic characteristics was performed by evaluating the band structure and DOS, depicted in Figures 5 and 6, respectively. As can be observed in Figure 5, the band structure profiles of all models are very similar. In addition, the top of the valence band and the bottom of the conduction band are located at the Г point, which indicates a direct band gap for the four models.

Vibrational properties
Besides the study of the electronic structure, a detailed study of the vibrational models, which is an important way to evaluate the reliability of the calculations since it depends on structural and electronic characteristics, was performed.
Considering the group theory, the BaMoO4 structure has 26 zone-center phonon modes: Γ= 3Ag + 5Eg + 5Bg + 5Au+ 5Eu+ 3Bu, where 3Bu are silent modes and 1Au+ 1Eu are acoustic modes. The active optical modes are 3Ag + 5Eg + 5Bg Raman-active and 4Au + 4Eu as IR-active. Table 3 lists the expected number of the active modes indicated by the group theory, as well as the experimental modes and theoretical results obtained by means of lattice dynamics calculation (LDC) according to studies found in the literature for comparison.
According to the results shown in Table 3, the computational simulations describe the Raman modes with appreciable accuracy if compared with the experimental observations. In addition, the DFT calculations provide higher precision than the lattice dynamical calculations reported by Jindal et al. 51 An earlier study based on DFT calculations, 37 showed some considerable divergences, mainly for R6 (Eg), R7 (Ag), and R8 (Bg), with frequencies of 308, 329, and 423cm -1 , respectively.
However, the systematic study of the HF% modification and basis set optimization resulted in values very close to those of the experimental sample.
The theoretical IR modes for both functionals are also in accordance with the experimental results reported by Xiao et al. 53 The third mode (Eu) was not visually observed in the cited experimental study, but it is in accordance with the number of active modes in the group theory.
The percentage of variation of each functional depends on each mode, but all of them showed a low deviation of the experiment.  In addition, a prediction of Raman and IR spectra (absorbance) were also elaborated, as shown in Figures 8a and 8b, respectively. The simulated Raman spectra (Figure 8a) are similar to the experimental spectrum reported by Basiev. 52 The most apparent differences are the higher intensity of the R3 mode (~100cm -1 ) and the separation of the R9 and R10(~350cm -1 ) modes in two separated bands in comparison with the experimental results. The similarity also can be observed in the reflectivity spectrum ( Figure   9) when compared with the experimental spectra reported by Xiao. 53 The Raman and IR theoretical spectra are a robust indication of the ability of well-elaborated models to support experimental observations. In view of the result, all the models can describe the system with good precision, where the method shows slightly better results. According to the present study, a flowchart was elaborated containing systematic steps for the development of computational models for solid materials, as can be seen in Figure 10. It is believed that these steps can be used to elaborate models with a higher degree of reliability when compared directly with experimental systems or to extract unknown properties from an experimental point of view; thus, they can be the starting point for studies with dopant addition or structural defects.

Figure 10.
Steps for the elaboration of theoretical models with a higher level of accuracy The first path proposed starts with a basis set optimization changing coefficients or exponents in a function of the total energy, showing agreement with the experimental or literature data. If the structure is good enough, it can be the final structure; otherwise, it is advisable to do an HF exchange modification and find the best percentage that describes the structure.
The second path proposed starts with the inverse of the first one. Initially, the HF exchange modification is done, and it is verified whether the structural and electronic properties are well described. If they are good enough, it can be the final structure; otherwise, a basis set optimization should be done, and then the structure is final.

Conclusion
Periodic DFT simulations were done using eight different hybrid functionals together with an all-electron basis set. The development of the best approach to describe accurately the structural, electronic and vibrational properties of a target system was proposed in two ways: 1) basis set optimizations and, when necessary, the HF% adjustment and 2) HF% adjustment followed by a base set optimization, when necessary.
Both strategies were applied to a BMO probe system, and the description, regardless of the path chosen, of the results is quite the same, with a precision of structural parameters, band gap energy, and vibrational spectra never before presented in theoretical studies of the respective materials.
The adjustment of the HF exchange functional rate should be done with caution, because the abrupt change in this percentage can result in an anomalous description of the electronic structure or perhaps other properties.
Furthermore, both pathways represented by the flowchart may help other researchers in more accurate calculations. It is also believed that the improvements generated by this methodology consist of an important step to a better description of the 2D and 1D systems because these models are usually based on the bulk results. In addition, this procedure can be applied independently of the applied package program.