Analytical approach to structural chemistry origins of mechanical, acoustical and thermal properties

ABSTRACT Crystalline matters with periodically arranged atoms found wide applications in modern science and technology. To facilitate the design of new materials and the advancement of existing ones, accurate and efficient models without relying too much on known inputs for predicting the functionalities are essential. Here, we propose an analytical approach for such a purpose, with only the knowledge of the structural chemistry of crystals. Based on the electrostatic interaction between periodically arranged atoms, the 1st, 2nd and 3rd derivatives of interatomic potential, respectively, enable a prediction of ten kinds in total of mechanical, acoustical and thermal properties. Over a thousand measurements are collected from ∼500 literatures, this results in the symmetric mean percentage error (SMPE) within ±25% and the symmetric mean absolute percentage error (SMAPE) ranging from 22%∼74% across all properties predicted, which further enables a revelation of bond characteristics as the most important but implicit origin for functionalities.


INTRODUCTION
Electrostatic interactions are ubiquitous in nature, having a vital significance in the formation of condensed matters.They provide robust interaction force among various charged species, such as elementary particles, molecules and proteins [1 ].These microscopic forces largely determine the macroscopic properties of matter.For instance, sound waves (mechanical waves) in a one-dimensional osci l lator are generated and transmitted by forces based on Hooke's law, while thermal resistance is the result of anharmonicity which is the high-order derivative of the interatomic forces.Therefore, a deeper understanding of electrostatic forces could shed new light on materials' properties ranging from structure formation to lattice dynamics, help improve existing functional materials and design new ones.
Macroscopic functionality of materials mainly originates from the microscopic functional units and the structure [2 ].At a microscopic level, lattice, charge, orbital and spin are the fundamental units [3 ], and each of them dominates the corresponding electrical and magnetic properties.Moreover, interactions between different units give birth to many unique properties, e.g.superconductivity [4 ] and thermoelectricity [5 ] through electron-phonon coupling.In fact, interatomic force is the key functional unit determining the mechanical, acoustical and thermal properties of materials.
The development of materials can be categorized into two perspectives.One focuses on structural geometry based on point group and space group, leading to the knowledge of the structural chemistry of ∼10 5 inorganics and 10 6 organics according to the Inorganic Crystal Structure Database (ICSD) and the Cambridge Crystallographic Data Centre (CCDC) databases.The other is a chemical algebraic approach based on quantum [6 ] or classical mechanics [7 ,8 ] for understanding the functionality.However, knowledge of functionality is far behind that of structural chemistry.Specifically, < 5% of known inorganics are documented with their thermal conductivities, which of course limits the exploration of thermal functional materials.
In order to facilitate the manipulation and design of functionalities at equilibrium, one has to know the energy of the system as the first step.Surely, atoms standing at their equilibrium positions in space Figure 1.Measurements vs. model predictions for mechanical, acoustical and thermal properties in crystals (detailed units in Fig. 3 ), along with the corresponding statistic error analyses.
corresponds to the lowest-energy state, while the difficulty lies in the explicit and accessible expression of microscopic interatomic energy.Quantitatively, first-principle calculations enable an estimation of the relationship between structure, chemistry, energy and property.However, it is sti l l chal lenging in accuracy, cost-efficiency and productivity particularly for complex structures [9 -11 ].
To balance the efficiency and accuracy of predictions, many approximations and models have been developed.For example, the L-J [12 ], EAM [13 ,14 ], Tersoff [15 ] and machine leaning potentials [16 ] give the effective potential based on the known binding energy; the Griffith model [17 ] gives the theoretical fracture strength based on the known surface energy; the Slack [18 ] model and Synder model [19 ] give the lattice thermal conductivity based on the known sound velocity and anharmonicity.However, these approaches largely depend on the amount and quality of known inputs and involve numerical fittings [20 ].Otherwise, the predictions could tend to be much less-accurate or even physically lessmeaningful.Therefore, there is always a desire for an efficient and accurate approach for predictions without many material-dependent inputs and fittings.
In this work, we develop an analytical model for predicting mechanical, acoustical and thermal functionalities of crystals, using MS Excel as the calculation tool and using the knowledge of structural chemistry from ICSD as the only input.Over 20 0 0 measurements involving ten kinds of properties (mostly at room temperature), are collected from ∼500 literatures, the symmetric mean percentage error (SMPE) of −7.8% and the symmetric mean absolute percentage error (SMAPE) of 32% validate the prediction accuracy (Fig. 1 ).
For each of the predicted properties, the resultant SMPE are within ±25% and SMAPE range from 22% ∼74%.Based on the electrostatic interaction between ionic cores, this model approximates the Coulomb, exchange and overlap integrals by analytical functions (Fig. 2 ).Holding the atomarrangement periodicity in crystalline matters for simplicity in interatomic potential, the first-, second-and third derivatives of which enable many lattice-dynamic related properties including elastic constants, sound velocities, Grüneisen parameter, linear thermal expansion coefficient, etc. (Fig. 3 ).Importantly, this approach reveals the key structural chemistry parameters for manipulating each individual functionality (Fig. 4 ).

The analytical model
To establish an analytical model (details in Supplementary Section 1), the following approximations are taken in this work: first, a chemical bond is divided into ionic ( U i ), covalent ( U c ) and metallic ( U m ) components; second, the ionic cores are treated as stationary point charges under the Born-Oppenheimer approximation.In this way, the electrodynamic interaction reduces to electrostatic interaction, and the potential between ionic cores has an analytical form of Coulomb potential, as shown in the first terms in Equations 1 a, 1 b and 1 c: where k is the Coulomb constant, CN is the coordination number, r is the interatomic distance and r 0 is the equilibrium bonding length.The Madelung constant ( α i ) is a geometry factor, which is used for determining the superposition of the long-range electrostatic interactions in the lattice.Similarly, the other geometry factors ( α c for covalent bond and α m for metallic bond) are defined in this work for determining short-range electrostatic interactions.Due to the short range, one can approximately consider only the repulsion between the nearest neighbors in which each atom has a CN nearest to the neighboring atoms.Thus, the repulsive potential components are proportional to the CN .It should be noted that the chemical bonds are localized in covalent compounds  while they are delocalized in metallic compounds.This difference leads the geometry parameter for metallic bonds ( α m ) to be further weighted by the CN .Eventually, α c = 10 CN 2 (100 + 19 CN ) and α m = α c / CN .The Q i , Q c and Q m are ionic, covalent and metallic components of ionic charge for ionic cores i and j , respectively ( Equation S2).
To avoid the collapse or decomposition of the lattice, the total potential requires a short-range repulsive potential of anti-bonding state ( E − ) for ionic component or a long-range attractive potential of bonding state ( E + ) for covalent and metallic components.Learning from the Coulomb, exchange and overlap integrals of hydrogen molecule ions [21 ], E ± are approximated to have analytical forms, as given by the second terms in Equations 1 a-1 c.Based on the exact definitions that alkalis are purely metallic, alkali salts are purely ionic and diamond-structured IVA elements are purely covalent, the constants A, B and C in Equation 1 a-1 c are respectively determined to be 1, 9 and 10, after fitting potential energy at equilibrium position ( Fig. S1).
Meanwhile, the structural factors are normalized using a pre-factor β.The total potential energy ( U ) becomes the sum of βE ± and Coulomb potential between ionic cores (Fig. 2 ).The pre-factor β can then be solved according to the equilibrium conditions ( Equation S1): the first and second derivatives of U at the equilibrium have, respectively, to be 0 and > 0. Since the three types of strong chemical bonds (ionic, covalent and metallic) are considered in this work, an arbitrary chemical bond is approximated here to be composed only of ionic ( p i ), covalent ( p c ) and metallic ( p m ) components, then the total potential is the summation of the individual contributions.It should be noted here that there is no strict boundary between chemical bonds.According to Pauling and Phi l lips [22 -24 ], the difference of electronegativity ( χ ) measures the polarity of the chemical bond, while the average of electronegativity ( χ ) measures the delocalization of the bonding electrons.The former determines p i and p c , and the latter determines p c and p m .Noted here that the bonding electrons are shared among coordination number bonds in crystals.Consequently, p i , p c and p m could be obtained according to Equation S3.
With basic structural chemistry parameters including atomic radius, equilibrium bond length ( r 0 ) and coordination number ( CN ) and so on, potential function and the charge Q of ionic cores at equilibrium can be determined ( Fig. S2).To simplify the derivatives process of the potential energy in this work ( Fig. S3), the lattice structure for modeling is approximated to be simple cubic, which is similar to the real structures for IV elements (C, Si and Ge), III-V semiconductors (Ga A s, GaSb, etc.), alkali salts (NaCl, KCl, etc.) and alkali metals (Li, Na, K, etc.).

Measurements vs. predictions
The relationship between interatomic force (stress) and distance (strain) can be obtained by taking the first-order derivatives of the potential energy (Fig. 3 a), of which the slope is identified as elastic modulus or elastic force constant around the equilibrium position r 0 .In this work, elastic ( E ), shear ( G ) and bulk ( B ) modulus are predicted.For simplicity, we use fourth-order stiffness tensor C ijkl [25 ] to quantify these moduli by considering the nearest and the next nearest neighbor interactions with a cubic lattice approximation.According to the Voigt and Reuss methods [26 ], E , G and B can be derived from three independent components ( C 11 , C 12 and C 44 ) of the stiffness tensor ( Section S2).It should be noted that the observed deviations between the predicted values and the experimental results primarily stem from the discrepancy between the approximated cubic structure of the current model and the actual material structures.For example, alumina and its analogs (e.g.Fe 2 O 3 , Sc 2 O 3 , etc.) crystallize in R3-CH structure.This structure induces extra elements of C 13 and C 44 for the stiffness tensor in addition to C 11 , C 12 and C 44 as given by Equation S5.In this case, it is reasonable to expect deviations when estimating the mechanical properties of a hexagonal structure using the stiffness tensor of a cubic structure.
The second-order derivative of the potential energy is the elastic restoring force constant of a vibrating object, which describes the propagation of sound waves through a continuous medium (Fig. 3 b).The propagating speed of sound energy (sound velocity, v ) [27 ] is highly anisotropic even in a cubic structure, because the bonding strength is highly directional.For comparison ( Section S2), predictions and measurements for v in polycrystals are included in Fig. 3 b.The transverse ( v T ) and longitudinal ( v L ) sound velocities depend on both the nearest-neighbor ( f 1 ) and next-nearest-neighbor ( f 2 ) force constant.For a cubic approximation, the v T is roughly half of the v L , because f 1 is 8 0.5 times f 2 based on Equations S4 and S10.The harmonic mean sound velocity [28 ] , which is formulated as Equation 2. In addition, acoustic impedance Table 1.Error analysis of symmetric mean percentage error (SMPE), symmetric mean absolute percentage error (SMAPE) and mean absolute percentage error (MAPE) for each of the properties, namely shear modulus ( G ), bulk modulus ( B ), elastic modulus ( E ), transverse sound velocity ( v T ), longitudinal sound velocity ( v L ), average sound velocity ( v s ), acoustic impendence ( Z ), Grüneisen parameter ( γ ), linear thermal expansion coefficient ( α) and lattice thermal conductivity ( κ L ).  [29 ] ( Z ) that is analogous to electrical impedance, can be predicted according to the product of mass density and v L .

Error analysis
Anharmonicity is the non-zero third-order and higher-order derivatives of potential energy.This results in a non-linear combination of amplitude during phonon collisions, which determines the phonon scattering and thermal properties of crystals [30 ] (Fig. 3 c).Grüneisen parameter [ γ = −(r 0 / 2 ) ( ... U / Ü )] at r = r 0 measures the anharmonicity, which is also crucial for thermal expansion ( α T ) and lattice thermal conductivity ( κ L ).The linear thermal expansion coefficient ( α T ) is determined by α T = (γ /B ) (C V /V ) / 3 according to the Grüneisen law [31 ], where C V is the heat capacity at constant volume and V is the average atomic volume.With the isotropic and cubic approximations, γ is reduced to a scalar from a tensor due to the rotational invariance, which is mainly determined by the fractions of bond characters as well as the coordination number: For lattice thermal conductivity ( κ L ) of a crystalline matter, the phonon gas model is utilized in this work to calculate the lattice thermal conductivity via κ L = 1/3 C V v s 2 τ , in which the phonon relaxation time ( τ ) is inversely proportional to the square of the Grüneisen parameter ( γ 2 ).

Error analysis
Since symmetric mean percentage error (SMPE), symmetric mean absolute percentage error (SMAPE) and the mean absolute percentage error (MAPE) can better i l lustrate the deviation between predicted values and measured values, we estimate the error analysis for each of the properties based on the definitions of SMPE, SMAPE and MAPE ( Equations S12a, S12b and S12c), as shown in Table 1 .All of the collected and the predicted data are listed in Tables S1-S3.
It is important to point out here that the relatively large deviation in κ L is related to dynamic factors rather than to the thermoelastic parameters predicted here.In regard to the dynamic factors, we learned from the Born-von Karman boundary condition [32 ] that the dispersion of acoustic phonons is a sine-type ( Equation S11a).This is a good approximation for those materials with large number of atoms in the primitive cell, but might not be applicable for materials with simple structures, such as IV elements and III-V semiconductors.
Consequently, although the sound velocity ( v s ) and Grüneisen parameter can be predicted quite well ( Fig. S4), the lattice thermal conductivity is statistically underestimated.This is attributed to the underestimation of the contribution of high-frequency acoustic phonons to κ L .This underestimation can be partially eliminated if a linear-type dispersion is used ( Equation S11b).The physical scenario is that the group velocity of high-frequency acoustic phonons is similar to that of low-frequency acoustic phonons, which is consistent with the actual phonon dispersion of IV elements, III-V semiconductors, alkali salts and alkali metals.In this way, the SMPE decreases from 88.2% for sine-type dispersion to 18.3% for linear-type dispersion ( Fig. S5).

Key structural chemistry parameters for functionalities
To better reveal the key structural chemistry parameters for manipulating the material's properties, we quantitatively analyze the implications that affect some mechanical, acoustical and thermal properties.As shown in Fig. 4 , chemical bond length ( r ) and characteristics along with charge quantity ( Q ) are important.
Common wisdom of electrostatic interaction is that a low charge quantity and a long equilibrium distance result in a weak interatomic force thus a low elastic module.This is reproduced in Fig. 4 a, showing that shear modulus ( G ) is clearly proportional to r 0 −4 and Q 2 for exemplary purely covalent IVA elements (C, Si and Ge), purely ionic alkali salts (NaCl,

Figure 2 .
Figure 2. Development of analytical model.Schematic of the analytical potential based on the electrostatic interactions.The first step is to approximate Coulomb/exchange/overlap integrals by analytical function E according to the molecular orbital theory.The pre-factor β can then be solved according to the equilibrium condition and the electrostatic charge Q can be parameterized using basic structural chemistry parameters.

Figure 3 .
Figure 3. Prediction of properties.Schematic of the differentiations of potential energy, which determine the mechanical (a), acoustical (b) and thermal properties (c), along with a comparison of measurements (detailed ∼500 references are given in Supplementary Materials).

Figure 4 .
Figure 4. Key structural chemistry parameter for functionalities.Structural-parameter-dependent measured (scatters) and predicted (curves) shear modulus (a and b), sound velocity (c and d) and Grüneisen parameter (e and f).