Unveiling oxygen vacancy impact on lizardite thermo and mechanical properties

Here, we performed a systematic DFT study assisted by the workflow framework SimStack for the mechanical and thermodynamic properties of the clay mineral lizardite in pristine and six different types of O vacancies configurations. In most cases, the defect caused a structural phase transition in the lizardite from the trigonal (pristine) to the triclinic phase. The results show that oxygen vacancies in lizardite significantly reduce the lattice thermal conductivity, accompanied by an elastic moduli reduction and an anisotropy index increase. Through the P–V relation, an increase in compressibility was evidenced for vacancy configurations. Except for the vacancy with the same crystalline structure as pristine lizardite, the sound velocities of the other vacancy configurations produce a decrease in these velocities, and it is essential to highlight high values for the Grüneisen parameter. We emphasize the great relevance of the punctual-defects introduction, such as O vacancies, in lizardite, since this microstructural design is responsible for the decrease of the lattice thermal conductivity in comparison with the pristine system by decreasing the heat transfer ability, turning lizardite into a promising candidate for thermoelectric materials


Methodology Atomic configurations
The Mg 3 (Si 2 O 5 )(OH) 4 lizardite structure is a layered trigonal crystal with a space group P31m 3,4 which has been studied in its pristine form and with six different types of oxygen vacancies, for which, we obtained the structural, mechanical, and thermodynamic properties.Following the Kroger-Vink notation 32 , we have organized the oxygen vacancies as follows: V × O1 -formed by removing one (basal) O atom that bonds to two Si atoms in the tetrahedra layer; V × O2 -by removing one (apical) O atom that bonds to a Si atom in the tetrahedra layer with three Mg atoms in the octahedra layer; V × O3 -by removing one OH from the top of the octahedra layer; the double vacancies: V × O1-O2 , V × O1-O3 , and V × O3-O3 are originated from the simultaneous formation of the vacancy pairs: , respectively.For clarity, Fig. 1 can be consulted to differentiate the six types of oxygen vacancies (top) and for a more precise notion of the layered lizardite structure (bottom).

Total energy calculations
Our ab initio calculations are based on DFT, as implemented in the Quantum ESPRESSO (QE) package 33 , using the ultrasoft method 34 within the generalized gradient approximation with the semilocal Perdew-Burke-Ernzerhof (PBE) 35 formulation for the exchange-correlation energy functional.As DFT-PBE cannot provide an accurate description of the nonlocal long-range van der Waals (vdW) interactions, we employ the vdW D2 correction proposed by Grimme 36,37 , which possess a good compromise between accuracy and computational coast 38,39 .The geometry optimization was carried out by the Broyden-Fletcher-Goldfarb-Shanno (BFGS) quasi-Newton algorithm 35 , with the convergence thresholds on forces and total energy for ionic minimization established in 10 −6 eV/Å and 10 −5 eV, respectively.The kinetic energy cutoffs related to the wave functions and charge density calculations in the summation were 0.8 keV and 8.2 keV, respectively.After convergence tests, a Monkhorst-Pack

SimStack workflow
SimStack 41 is a robust workflow framework that ensures the reproducibility and transferability of the simulation protocols 42 .Additionally, it simplifies the creation of custom-tailored simulation protocols using various computer simulation approaches.Here we use the Workflow Active Node (WaNo) DFT-QE 43 , which saves time by automating and reducing protocol complexity, permitting the monitoring of multiple sets of calculations for independent DFT protocols for different Lizardite configurations.The right side of Fig. 1 depicts this workflow, streamlining the optimization and elastic constants calculations using QE code.The output of DFT simulations for each Lizardite system automatically generates .ymlfile database, which we automatically may transfer to a Colab notebook to query data and calculate the lattice thermal conductivity κ L .To meet the FAIR principles 44 , we make all input data and the Colab notebook available on the following repository github.com/KIT-Workflows/Lizardite.

Property analyses
To deepen the understanding of lizardite, we have systematically analyzed its main mechanical and thermodynamic properties, including the lattice thermal conductivity.Further details on the properties are provided in the Supporting Information (SI) material.In terms of mechanical properties, we had considered the elastic stiffness ( c ij ) and compliance ( s ij ) constants (and their relations) as implemented in the Thermo_pw package 45 , checking the stability criterion when the eigenvalues of the elastic stiffness matrix were greater than zero [46][47][48][49] .Posteriorly, assuming the relaxation process under hydrostatic pressure 50 , the Voigt and Reuss bulk ( B V , B R ) and shear ( G V , G R ) moduli, respectively, were obtained from c ij and s ij .From which, the Hill bulk ( B H ) and shear ( G H ) moduli are given by averaging the corresponding values 50 : B(G) H = 1/2(B(G) V + B(G) R ) .Considering Hill's modulus and the mathematical theory of elasticity 51 , the Young's modulus ( E H ) and Poisson's ratio ( ν H ) were also considered.At the same time, the anisotropy of the material was obtained from the universal elastic anisotropy index ( A U ) 52 .Finally, we also have the possibility of estimating the ductility or brittleness of material through the criterion suggested by Pugh 53 , for which ratios B H /G H larger than 1.75 represent a ductile material.
Considering a P-V equations of state, the thermodynamic properties were obtained from the structural relaxation process under a pressure range from 0.0 to 8.0 GPa.The isochoric heat capacity ( C V ) was achieved by post-processing with Thermo_pw package 45 , considering the Debye temperature ( D ) within the Debye model.To get D , the elastic constant data was considered through the Voigt-Reuss-Hill average of the bulk and shear moduli to calculate the average sound velocities.Following Anderson 54 , where D depends on the average sound velocity ( v av ), which is given in terms of the longitudinal ( v ℓ ) and transverse ( v t ) sound velocities, ).On the right side: We display the workflow that we used to do all the simulations and calculate the properties of Lizardite.In the steps DFT-QE and DFT-QE_1, we optimized geometry and the elastic constants.These steps created a database that we then fed into a Colab notebook to calculate lattice thermal conductivity ( κ L ).
Vol:.( 1234567890 .Similarly, the Grüneisen acoustic constant ( γ ) was calculated as a v ℓ and v t function, following Belomestnykh 56 .See more details about thermodynamic properties in SI.
For the lattice thermal conductivity ( κ L ), we have considered two models: (i) for the pristine lizardite, the Slack model 57,58 has been adopted since it is widely appropriated for defect-free crystals and lattice thermal resistance resulting only from intrinsic phonon-phonon interactions [59][60][61][62] .The κ L is obtained from D and γ , which can be obtained from lattice dynamic calculations or experimental measurements 57,61 .However, in the present work, we have obtained κ L considering directly the D and γ calculation (see SI for mathematical expressions), in complete agreement with Xia et al. 63 approach.(ii) For lizardite with oxygen vacancies, κ L is calculated following the seminal works of Klemens 64,65 and corroborated by Callaway et al. 66 and Abeles 67 .Basically, κ L for lizardite with an oxygen vacancy type ( , it is given in terms of κ L for the pristine lizardite and u is written as a function of D , κ L , v av , and Ŵ (which is given in terms of the concentration and mass of the atom type i).Finally, the minimum lattice thermal conductivity of lizardite with oxygen vacancy ( κ V × α (min) ) is obtained using the minimum value of κ L for pristine lizardite ( κ L(min) ), according to Clarke 68 .More details are presented in the SI material.
In Table 1 we present the accuracy of the Belomestnykh approximation used in the present work to determine the Grüneisen parameters, in relation to the DFT calculation , for some examples from literature for structure simpler than the lizardite one.It is observed that the results from the Belomestnykh approach are in reasonable agreement with the DFT calculations.

Geometric optimization
We have performed the structural optimization of pristine and oxygen vacancies lizardite within the stacking layer context, each of which is composed of two sheets: a Mg-centered octahedra over a Si-centered tetrahedra sheet, as shown in the bottom of Fig. 1.The six types of O-vacancies considered are given by: V , and V × O2-O3 , as shown at the top of Fig. 1 (from a) to f), respectively).The structural relaxation procedure of these systems and the subsequent comparison between the pristine and defect systems are of great relevance for establishing the main structural changes resulting from vacancy formation.Thus, the main structural properties discussed here, e.g., lattice parameters (a, b, c), angles ( α , β , γ ), volume, density, and interlayer distance, are also presented in SI (Table S1).
For the pristine lizardite, we have obtained a trigonal lattice with parameters: a = b = 5.276 Å and c = 7.117 Å ; α = β = 90 • and γ = 120 • , and a layer distance of 1.783 Å , results that are in excellent agreement with previous works 5  is the number of O atoms removed from a tetrahedron vertex (the top of the octahedra sheet), and zero otherwise.The E o v values, along with the percentage differences in the volume ( Vol ) and interlayer distance ( id ) concerning the pristine system, are shown in Table 2.
From Table 2, it is observed that the lowest formation energy is obtained for V × O2 vacancy-type, followed by V × O1 and V × O3 .In contrast, the formation energy of a double vacancy is approximately the sum of the energies of the two single vacancies that originated it.This result is directly associated with the fact that only V × O2 preserved the trigonal lattice structure in agreement with the pristine lizardite.At the same time, the other vacancy systems changed the lattice structure type from trigonal to triclinic.For V × O2 , the removal of an oxygen atom between the octahedra and tetrahedra sheets causes a moderate volume contraction ( Vol = −2.67%),which is directly linked to an approximation between the sheets ( id = −5.89%)concerning pristine lizardite, all without chang- ing the lattice angles of the trigonal structure.On the other hand, all different vacancy types represent structural (1) Table 1.Comparision of the Grüneisen parameters obtained from the Belomestnykh approximation with the DFT calculation for some compounds.. (a) Ref. 69 (b) Ref. 70 (c) Ref. 71 (d) Ref. 63 Compound γ (DFT) γ(Belomestnykh) d) 1.65 (d) PbSe 2.66 (d) 1.69 (d) PbS 2.46 (d) 1.67 (d) Vol.:(0123456789) www.nature.com/scientificreports/changes in the pristine lizardite that, to a greater or lesser extent, lead to distortions that penalize the stabilization of the system with point defects, with volume variations ranging from a slight expansion ( Vol = 0.39%) for V × O3 , due to the small approximation between the sheets ( id = −0.17%), to significant volume contractions ( Vol = −6.45%)for V × O1-O2 , due to the largest approximation between sheets ( id = −9.48%),all of them within the structural context of the triclinic lattice.In short, the six vacancy types studied here represent, together with pristine lizardite, a very diversified set from the point of view of possible structural alterations.At the same time, they are energetically viable possibilities to happen in a practical application.

Mechanical properties
We started the study of the mechanical properties of pristine and vacancy-type lizardite from the calculation of the elastic stiffness constants c ij and the elastic compliances constants s ij , which are presented in the SI (Table S2).From this, we performed a two-level characterization of the main mechanical properties, a primary and complementary description of the mechanical response properties of the studied materials.Firstly, from the matrices [ c ij ] and [ s ij ], it is possible to check the mechanical stability criterion, which establishes that all eigenvalues of these matrices must be greater than zero to achieve the elastic stability [46][47][48][49] .All structures studied here meet this criterion, showing elastic stability in pristine and vacancy-type constitutions.Our results are in excellent agreement with literature for the case of pristine lizardite 8 .
Within the context of basic mechanical properties, we can check on (i) the in-plane and out-of-plane mechanical stiffness, (ii) the fracture resistance (or rigidity), and (iii) the stress anisotropy of the lizardite configurations with and without defect, which turns out to be interesting in the scenario of how pristine and vacancy-types lizardite resists concerning the stress application.For (i), an estimate can be made through the ratio c 11 /c 33 , whose values are presented in Table 3.We obtain a ratio c 11 /c 33 = 1.707 for the pristine configuration, which implies a greater mechanical stress resistance in-plane than out-of-plane.As expected, for configurations with vacancy, we have a stiffness dependence on the sites from which the O atom has been removed.Except for V × O1-O3 , for which the in-plane and out-of-plane stiffness are equivalent, all other cases of vacancy cases behave like the pristine lizardite, i.e., are weaker bonds along the [001] direction.For (ii), the estimate is made by the ratio c 66 /c 44 , which represents how much the basal plane is more resistant to fracture than the axis perpendicular to this plane when shear stress is applied.The rigidity values ( c 66 /c 44 ) are presented in Table 3.We observed that, with pristine liz- ardite, the basal plane is more (less) resistant to fractures than the c-axis for V × O1 , V × O3 , V × O1-O3 and V × O2-O3 (V × O2 and V × O1-O2 ).Finally, for (iii), we can estimate the stress anisotropy from the linear compressibility ( β ), which represents the relative variation in the length of a line when the body is subjected to unit hydrostatic pressure 73 and depends on the elastic compliance constants s ij .Following the Nye work 49 , we have estimated β for trigonal and triclinic crystal systems in each axial direction, as follows: β 1 = s 11 + s 12 + s 13 (both), β 2 = β 1 (trigonal), β 2 = s 12 + s 22 + s 23 (triclinic), β 3 = 2s 13 + s 33 (trigonal), and β 3 = s 13 + s 23 + s 33 (triclinic).The β values are also presented in Table 3, in which we observe that the two configurations with trigonal lattice (pristine and V × O2 ) present isotropic in-plane stiffness ( β 1 = β 2 ) as expected, due to the crystalline symmetry.In contrast, the other vacancy-types lizardite configurations show planar anisotropy.

Table 2. The formation energy of O vacancy per oxygen atom ( E o
v ), the percentage difference in the volume of structures concerning the pristine one ( Vol ), and the percentage difference in the interlayer distance ( id ) for the pristine system, for different types of vacancies in lizardite.Negative (positive) Vol and id values represent a decrease (increase) in the volume and interlayer distance of the structures with pristine lizardite.. www.nature.com/scientificreports/Complementarily, using the elastic constants tensor and the Voigt, Reuss, and Hill models, we have obtained the bulk ( B V , B R ), shear ( G V , G R ), and Young ( E H ) moduli, as well as the Poisson's ratio ( ν H ) and the universal elastic anisotropy index ( A U ) for a complete mechanical property characterization of the pristine lizardite and with different O vacancy-types.All these properties are shown in Table 4. First, we observe that the vacancy formation leads to a decrease in the B V , B R , G V , G R , and E H values, which is directly associated with the degree of the constitutional importance of the O atom(s) removed to form the vacancy.For example, the smallest elastic modulus values occur for a double vacancy, specifically for V × O1-O3 , when two O atoms are removed, one bonded to two Si atoms at the base of the tetrahedra sheet and the other at the top of the octahedra sheet.Even more specifically, the V × O3 vacancy type leads to elastic moduli smaller than V × O1 one, i.e., removing a O atom from the top of the octahedra sheet decreases the elastic moduli more than removing a O atom from the base of the tetrahedra sheet.All studied configurations prefer to be shear deformed since the shear moduli are smaller than the bulk moduli.In the case of the ν H values, we observed a slight variation ( ν H = 0.30 ± 0.01 ).Similarly to the small variation of the ratio B H /G H , i.e., between 2.01 and 2.24, which, according to the criterion suggested by Pugh 53 , characterizes all our configurations as ductile, since the values are more significant than 1.75.Finally, the smallest A U values occur for the trigonal structures (pristine and V × O2 ), which are isotropic in-plane in terms of stiffness, in contrast to triclinic (all other cases), in complete agreement with the linear compressibility results.

Thermodynamic properties
Once the mechanical properties of pristine and vacancy-types lizardite were established, there is a need for a thermodynamic property analysis taking into account pressure and temperature variations to trace the real implications of O vacancies in the lizardite constitution.First, in terms of the pressure variation, we obtained the P-V equations of state from the structural relaxation under pressure (from zero to 8.0 GPa), which are presented in Fig. 2. The normalized volume (relative to the structure at zero pressure) decreases with increasing pressure for all configurations, and the P-V curve shows a greater slope with O vacancies, i.e., the lizardite configurations with vacancies present a more significant decrease in volume with increasing pressure, which is directly correlated with the mechanical properties elucidated by the bulk, shear, and Young modulus values (see Table 4).Up to 4.5 GPa, the P-V curves of the pristine and V × O2 configurations have the same response to pressure increase, keeping in mind that both have the trigonal crystalline structure; for values larger than 4.5 GPa, we observe that the V × O2 configuration starts to have a more accentuated response to the pressure increase, which is associated with the structural change that is evident from the lattice parameters (a and c) variation in the insert of Fig. 2. All other configurations with vacancy, possessing the triclinic crystalline structure, are more susceptible to pressure variations, respond more accentuated, and are more compressible.The most significant decrease in volume with the pressure increase occurs for V × O1-O3 and V × O2-O3 (whose P-V diagrams are practically overlapping), which are the two configurations with double vacancies that have the highest (and identical) vacancy formation energies.These three lizardite configurations (V × O2 , V × O1-O3 , and V × O2-O3 ) had a change in the expected linear behavior between volume and pressure, showing an increase in the slope from 4.5 GPa, which is associated with the rise in compressibility that occurs mainly on the c axis, which implies a reduction in the interlayer space and significant changes in the O-H-Mg angle, positioning O-H outside the normal direction.
For the temperature variation, we performed the isochoric heat capacity ( C V ) calculation, which first requires obtaining D from the Debye model, whose expression depends on v av , given by the angular average of the sound velocities calculated for each propagation direction, i.e., v ℓ and v t , both being important for the Grüneisen acous- tic constant ( γ ) determination.Thus, before we get into C V , we briefly discuss the v ℓ , v t , v av , γ , and D results for lizardite configurations with and without vacancies, as presented in Table 5.Our sound velocity results are in excellent agreement with the experimental results obtained for pristine lizardite.For example, Kern et al. 74 have obtained v ℓ and v t in the serpentine, at low pressure, around 7.2 and 4.8 km/s, respectively.According to Christensen 75 , which presented v ℓ values in the range of 5.0-8.3km/s, depending on the sample serpentiniza- tion.Or the work by Seipold and Schilling 21 , for which v av varies from approximately 5.02 km/s under ambient conditions to 4.5 km/s at 900 K.
From Table 5, we observe that, in general, the different vacancy-types formation is associated with a sound velocity decrease, except for the pristine and V × O2 configurations, the only vacancy configuration that preserves Table 4.The bulk ( B V , B R ), shear ( G V , G R ), and Young ( E H ) moduli (in GPa); the Poisson's ratio ( ν H ) and the universal elastic anisotropy index ( A U ) of the lizardite with and without vacancies, where the subscript- labels: V, R, and H refer to the Voigt, Reuss, and Hill models, respectively.In parentheses, we have added the respective property values of the pristine lizardite from literature 8 .the trigonal crystalline structure.Unanimously, we found that v ℓ is a bit more than 1.8 times v t , which means that sound easily propagates in the plane of the layers than along the direction of the stacking layers ([001]).Our v ℓ /v t results agree with the experimental results of Christensen 75 , by which this ratio increases systematically from 1.78 to 2.21 with increasing sample serpentinization.We correlate the mechanical properties (Table 4) and the sound velocities.For example, V × O1-O3 has the lowest sound velocities in correspondence to the smallest volume, shear modulus, and Young's modulus.For γ , we observed a slight variation among configurations, leading to an average value of 1.76.On the other hand, D presents smaller values for the structures with vacancy, being the smallest values for the configurations with a triclinic crystalline lattice than with a trigonal one.
After these preliminary calculations, we performed the C V calculations concerning the temperature, from 0 to 800 K, for the lizardite configurations with and without O vacancies, which are presented in Fig. 3a.We observed the same behavior trend in the C V versus temperature curves for all lizardite configurations, with an exponential increase in C V up to approximately room temperature and subsequent logarithmic growth for C V  at higher temperatures, with a trend asymptotic to converge to the Dulong-Petit limit at high temperatures.The most significant C V differences among the configurations occur in the 100-300 K interval, where we observe that, with the V × O2 exception, the C V values are more prominent for the lizardite configurations with vacancy than that for pristine lizardite.

Lattice thermal conductivity
From γ and D , we can calculate the lattice thermal conductivity ( κ L ), as well as the minimum thermal con- ductivity ( κ min )), which are shown in Table 5, for the pristine and O vacancy-types lizardite.The behavior trend between κ L and κ min is the same across the configurations studied, i.e., for κ min (in W/(m.K)): pristine 0.40), in agreement with κ L as a temperature function, in Fig. 3b.As evidenced, the O vacancies formation in lizardite leads to a decrease in heat transfer ability compared to the pristine configuration.
Thus, from Fig. 3b, we observe that the O vacancies produce a significant κ L decrease, following the sound velocities (Table 5), which represents our main result here, since the κ L values, the material conducts less heat energy.Furthermore, we highlight that this important thermoelectric parameter, directly linked to determining the energy conversion efficiency of thermoelectric materials, has a very low variation with increasing temperature for systems with the vacancy in contrast to the pristine system.Consequently, the introduction of the defect (O vacancies) in lizardite is part of the successful strategy of intelligent microstructural design, which plays an essential role in the task of obtaining materials with high phonon-phonon scattering rates, i.e., materials with significant Grüneisen parameters.These behavior are in agreement with the results obtained by Shen et al. 76 and Jia et al. 77 , which studied the effect of vacancies on the lattice thermal conductivity of CuGaTe 2 and In 2 Te 3 -InSb, respectively.

Conclusions
In the present work, we performed a systematic first-principles study, based on DFT, for the main mechanical and thermodynamic properties of the clay mineral lizardite in pristine and six different types of O vacancies configurations.After the structural relaxation process, except for a vacancy configuration, we obtained a change in the crystalline structure of the vacancy configurations, going from trigonal (pristine) to triclinic.As the main result found, we highlight that oxygen vacancies in lizardite significantly reduce the lattice thermal conductivity, accompanied by an elastic moduli reduction and an anisotropy index increase.Through the P-V relation, an increase in compressibility was evidenced for vacancy configurations.Except for the vacancy with the same crystalline structure as pristine lizardite, the sound velocities of the other vacancy configurations produce a decrease in these velocities, and it is essential to highlight high values for the Grüneisen parameter.Finally, we emphasize the great relevance of the punctual-defects introduction, such as O vacancies, in lizardite, since this intelligent microstructural design procedure is responsible for the decrease of the lattice thermal conductivity in comparison with the pristine system, which decreases the heat transfer ability, giving rise to an essential candidate for thermoelectric materials.

Figure 1 .
Figure 1.Left bottom: Mg 3 (Si 2 O 5 )(OH) 4 one-layer representation, consisting of two sheets: an octahedra sheet, with Mg atom at the center of the octahedron; over a tetrahedra sheet, with Si atom at the center of the tetrahedron.Each Mg (Si) atom is surrounded by six (four) O atoms, and the atomic site on the octahedron top is occupied by OH.Top: the six types of oxygen vacancies: (a) removal of one (basal) O atom that bonds two Si atoms at the base of the tetrahedra sheet (V × O1 ); (b) elimination of one (apical) O atom that bonds one Si atom in the tetrahedra sheet with three Mg atoms in the octahedra sheet (V × O2 ); (c) elimination of one OH from the top of the octahedra sheet (V × O3 ); (d) double vacancy originated from the simultaneous V × O1 & V × O2 vacancies (V × O1-O2 ); (e) double vacancy originated from the simultaneous V × O1 & V × O3 vacancies (V × O1-O3 ); and (f) double vacancy originated from the simultaneous V× O2 & V × O3 vacancies (V × O2-O3).On the right side: We display the workflow that we used to do all the simulations and calculate the properties of Lizardite.In the steps DFT-QE and DFT-QE_1, we optimized geometry and the elastic constants.These steps created a database that we then fed into a Colab notebook to calculate lattice thermal conductivity ( κ L ).
. The formation of different types of vacancies in lizardite can lead to further structural parameter responses.In this context, we initially estimated the formation viability of the O-vacancies studied based on the formation energy ( E o v ) as established by Emery and Wolverton 72 : where E(Mg 3 Si 2 O (5−α) (OH) (4−β) ) is the total energy of Mg 3 Si 2 O 5 (OH) 4 with vacancy, µ O and µ H are the O and H chemical potentials, respectively, E(Mg 3 Si 2 O 5 (OH) 4 ) is the total energy of the pristine Mg 3 Si 2 O 5 (OH) 4 , and α ( β )

Figure 2 .
Figure 2. Variation of normalized volume in relation to pressure for pristine and O vacancy-types lizardite.Two inserts are shown for the normalized (scaled by a 0 , which consists of a at zero pressure) lattice parameters (a and c) variation in relation to pressure to explain the escape from a linear behavior for the P-V diagrams of the V × O2 , V × O1-O3 , and V × O2-O3 configurations.

Figure 3 .
Figure 3. (a) The isochoric heat capacity ( C V ) of pristine and O vacancy-types lizardite as a function of temperature.(b) The lattice thermal conductivity ( κ L ) as a temperature function of pristine and O vacancy-types lizardite.
51,55 turn, v ℓ and v t depends on the mechanical properties ( B H and G H )51,55