Density functional theory predictions of the mechanical properties of crystalline materials

The mechanical properties of crystalline materials are crucial knowledge for their screening, design, and exploitation. Density functional theory (DFT), remains one of the most effective computational tools for quantitatively predicting and rationalising the mechanical response of these materials. DFT predictions have been shown to quantitatively correlate to a number of experimental techniques, such as nanoindentation, high-pressure X-ray crystallography, impedance spectroscopy, and spectroscopic ellipsometry. Not only can bulk mechanical properties be derived from DFT calculations, this computational methodology allows for a full understanding of the elastic anisotropy in complex crystalline systems. Here we introduce the concepts behind DFT, and highlight a number of case studies and methodologies for predicting the elastic constants of materials that span ice, biomolecular crystals, polymer crystals, and metal – organic frameworks (MOFs). Key parameters that should be considered for theorists are discussed, including exchange – correlation functionals and dispersion corrections. The broad range of software packages and post-analysis tools are also brought to the attention of current and future DFT users. It is envisioned that the accuracy of DFT predictions of elastic constants will continue to improve with advances in high-performance computing power, as well as the incorporation of many-body interactions with quasi-harmonic approximations to overcome the negative effects of calculations carried out at absolute zero.


Introduction
8][9] The mechanical properties of a material tell us a lot about its longevity, its ability to withstand damage, and its potential applications. 10We are now in a period where shape-memory alloys, 11 flexible electronics, 12,13 ferroelastics, [14][15][16] and other mechanically-complex and exciting materials are becoming commonplace.As well this, modelling tools such as finite element analysis (FEA) often require prior knowledge of properties such as the Young's modulus to accurately predict the behaviour and performance of materials. 17,18More recently, expansions of Pugh mechanical analysis have been used in the derivation of hardness descriptors, 19 in the search for new materials for hard coating applications.This has also guided researchers in the search for an inorganic compound with a hardness greater than diamond. 202][23] By knowing the full anisotropic elastic tensor the elastic response of composite materials can be predicted when combined with mathematical homogenization theories.This has allowed for the design of such materials with engineered stiffness. 24,25Additionally, elastic properties are widely used in the field of geophysics, where acoustic velocities can be used to interpret seismic data. 26,275][36] Herein we discuss the various DFT methodologies that can be used to calculate the mechanical properties of crystals, as well as exploring other computational chemistry methods that are used today in this endeavour.The limitations of and challenges facing DFTbased predictions mechanical properties are discussed, as well as the exciting future that lies ahead if these challenges are overcome.

Why are mechanical properties important?
The elastic and mechanical properties of materials can be characterized by calculating second-order elastic constants.These are the measure of proportionality between stress (or energy) and strain within the region of Hooke's law.Elastic constants are broadly determined by applying a strain to a crystal, and measuring the resulting stress. 37Elastic constants of molecular crystals can be determined experimentally using ultrasonic methods that rely on elastodynamics. 38In recent years, theoretical studies have been used to estimate elastic constants with the help of firstprinciples methods in order to guide experiments and alleviate pressure on trial-and-error investigations.An increase in computational power has provided an avenue for DFT-led studies that can screen for materials with optimum properties prior to being fabricated or grown in the lab.
Research into the mechanical behaviour of a new class of solid-state materials is central to both the design and optimal performance of potential technological applications.Take for example metal-organic frameworks (MOFs) where theorists and experimentalists can examine the elasticity of these hybrid frameworks by examining their Young's modulus, Poisson's ratio, bulk modulus and shear modulus. 39Also crucial are discussions on their hardness, plasticity, yield strength and fracture behaviour.For these materials predicted elastic properties such as compressibility and bulk moduli can be compared to high-pressure X-ray crystallography. 40pectroscopic ellipsometry has also been used to estimate the elastic moduli of MOF nanoparticles and deposited films. 41anoindentation has emerged as a key technique for quantifying the mechanical properties of crystalline materials (Fig. 1), [42][43][44] and recently nanoindentation data has been used to train machine learning algorithms. 45n the pharmaceutical industry, it is of utmost importance to understand the elastic and mechanical properties of active pharmaceutical ingredients (APIs).For the API, mechanical properties govern physicochemical properties such as solubility, tabletability, stability and the bioavailability of a drug substance.It is especially important as roughly one in two APIs can exist in multiple solid forms, with each form markedly showing different physicochemical and mechanical properties. 47,480][51] In the case of polymorphism, the solid-solid transformations can also cause formulation problems.Since polymorphs of the same molecular crystal have differences in interaction energies, other polymorphs tend to transform into the polymorph with the least free energy, and therefore the most stable polymorph.The resulting polymorphs can have undesirable properties, such as in the case of ritonavir and rotigotine. 52,53Unwanted phase transformations can affect the drug stability during its lifespan and in the handling of the drug, particularly if a shearing stress is applied.The piroxicam-succinic acid co-crystal for example is formed via the application of mechanical stress to the two components, but undergoes decomposition when shearing occurs.These process-induced transformations are difficult to predict and control due to lack of understanding of the mechanochemical process at an atomistic level.DFT calculations can be a crucial tool to understand and quantify polymorph stability, and can be used to study the interactions between APIs, co-formers, and excipients in both amorphous and crystalline environments.

What is density functional theory?
For single electron systems (2D potential well, hydrogen atom) it is possible to obtain an exact solution to the Schrödinger equation. 54However, for any system more complicated than this it is impossible to precisely solve the Schrödinger equation and obtain a mathematical description of the system i.e. the wavefunction. 55he term DFT comes the fact that the functional used in DFT calculations is the electron density, 56 which is itself a function of space and time (mathematically, a functional takes a function and gives a resulting scalar value).The Hohenberg-Kohn theorem 57 tells us that the total ground state energy of a many-electron system is a functional of the density.The total energy of the system is written in terms of a number of individual energy contributions, 58 each of which are functionals of the charge density: • ion-electron potential energy • ion-ion potential energy • electron-electron energy • kinetic energy • exchange-correlation energy.The most computationally challenging energy contributions are the kinetic energy and the exchangecorrelation energy.The kinetic energy is calculated using the Kohn-Sham orbitals. 59Generally, these do not correspond to actual electron orbitalsthey are orthonormal orbitals.The Kohn-Sham orbitals map the system of interacting electrons on to a system of non-interacting electrons moving in an effective potential.The exchange-correlation energy accounts for the exchange interaction due to repulsion between electrons with parallel spins, and the correlation interaction, which is the correlated motion between electrons of antiparallel spins due to their mutual coulombic repulsion.In its simplest implementation, exchange-correlation effects are treated via the Perdew, Burke, and Ernzerhof (PBE) 60 implementation of the generalised gradient approximation (GGA). 61GGA builds on what is known as the local density approximation (LDA), by considering both the local electron density and its gradient, as the electron density can vary rapidly over a small region of space.

Long-range considerations: what is dispersion-corrected DFT?
Considering dispersion corrections in DFT calculations is important as the exchange functionals have difficulty in reproducing long-range behaviour, even those explicitly parameterised for long range behaviour. 62This is because approximations such as GGA within these functionals 63 result in potential energy surfaces that lack sufficient crystal packing effects (CPEs) without the addition of a dispersion correction. 64hese are influenced by van der Waals forces that are better estimated through dispersion corrections as their inclusion follows a simple formula where the displacement energy is added to the exchange correlation function so as not to be included in the initial DFT calculations.
A vital aspect to the dispersion correction is dampening, with methods that lack an adequate dampening failing to give consistent results for crystal structure and energies. 65The dampening function within the dispersion correction determines the range at which the dispersion correction acts 66 as well as the steepness of the cut off of the dispersion correction. 67The dampening function means that the dispersion effects approach 1 at long distances, this meaning that it is purely a dispersion interaction, but at as the distance between the dipoles shortens the dispersion correction gets dampened eventually going to zero. 68This means that there is no dispersion effect at shorter distances where the XC functionals perform better.The dampening effects are also intended to reduce double counting effects. 69Further studies are needed on which dampening function performs best with each dispersion correction and exchange functional. 68However with increasing numbers of atoms in the unit cell the difficulty arises in deciding where dampening should take place, as this needs to be symmetric. 70ispersion corrections can be calculated in several ways.This can be with a pairwise approach, a three-body approach or a many body approach.The pairwise approach which is chosen in dispersion corrections schemes such as Grimme-D2 (ref.71 and 72) and Tkatchenko-Scheffler (TS) 73 summates over the C 6 R −6 potential where R is the atomic distance and the C 6 is the dispersion coefficient.For the Grimme-D2 method these are multiplied by a global scaling factor. 71where the TS scheme calculates the pairwise dispersion energy using the formula: where f damp is the dampening function, R ij is the distance between atom i and j and R 0 ij is the sum of the van der Waals equilibrium radii. 74,75These methods diverge as the DFT-D2 method developed by Grimme utilises the atom-type model whereas the TS model utilises the atomic volume.This use of the atomic volume allows for the approximation of manybody effects. 70An example of the three bodied approach is used in the Grimme-D3 correction scheme.The simplest way to employ a three-body approach is to use a non-additive third order Axilrod-Teller-Muto (ATM) term: where θ represents the internal angles formed by R AB , R BC and R CA . 76This uses a C 9 triple dipole constant that can be approximated in two possible ways, either through integrating the average dipole polarizability at an imaginary frequency for all three atoms, or through calculating the square root of the three-atoms dispersion coefficient. 76This three-body treatment has been found to contribute minimally to the overall dispersion energy with one paper saying that the three body treatment contributes <5-10% of the dispersion energy, 77 and another paper stating that the three body effect contributes 7.2% of the lattice energy. 78he many body dispersion method, as used in the many body dispersion scheme (MBD), builds on the pairwise TS approach and addresses the fact that the nature of longrange energy is many-body in nature. 74The main drawback of the MBD scheme is the high computational cost 70 which is due to the fact that within this method it involves having to calculate both pairwise and three-body dispersion energy utilising the formula: This is where E (2) is a pairwise and E (3) is three body energies. 65This is extended up to N atoms up to the Nth order this is coupled with a self-consistent screening in order to reduce the error. 79

Elastic properties of electroactive materials
Our work to date studying biomolecular crystals using DFT relies on the Vienna ab initio simulation package (VASP) 80 which uses plane wave basis sets, 81 and the projector augmented-wave (PAW) method. 82However, many DFT packages, such as CASTEP, 83 ABINIT, 84 and CP2K 85 also have efficient schemes for calculating elastic properties.For VASP standard PAW pseudopotentials are used in all calculations, as supplied with the software.As our focus has been on the piezoelectric properties 86 of crystals, elastic constant calculation is part of a calculation workflow, as the elastic stiffness constants c kj , are required to calculate the piezoelectric strain constants d ik .The elastic compliance (s kj ) can easily be derived from the stiffness and measured alongside electromechanical properties using impedance spectroscopy, 87,88 as can the various elastic moduli.
Using the piezoelectric charge coefficients, e ij , which are also calculated directly by VASP, 89 and the elastic stiffness constants, c kj , we can calculate the more useful piezoelectric strain coefficient, d ik , using the relationship A finite difference method can be used to calculate the elastic stiffness tensor, with every atom in the unit cell being displaced along each Cartesian axis by a default value of ±0.01 Å.A Γ-centred k-point grid increases mechanical stability, with an observed negligible dependence of predicted elastic constants on the number of k-points (once N > 1). 90For organic materials high plane wave energy cutoffs of up to 1000 eV are recommended allow the stress tensor to fully converge due to the presence of oxygen and nitrogen atoms. 91,92Young's moduli can be derived from the stiffness matrix components.This can be done using the approximations of Nye, 93 or Voigt-Reuss-Hill (VRH). 94,95here the software does not automatically derive elastic moduli the ELATE online tool 96 is recommended for quick validation of mechanical stability, VRH derivation, and 3D visualisation of elastic properties.The ElAM program 97 can also be used to derive and visualise elastic properties and their anisotropy, with more thorough options for plot generation and customisation.
Using this methodology, we have calculated to a high accuracy versus experiment the elastic constants of amino acid 90,[98][99][100] and peptide crystals, 101,102 co-crystals, 103,104 and biominerals 105 (Fig. 2).We have also recently calculated the elastic properties of large protein crystals using classical force fields.Elastic constants of the transmembrane protein ba3 cytochrome c oxidase, as well as lysozyme, and aldehyde dehydrogenase were predicted using the classical CHARMM forcefield model for the protein, ions and water with structures calculated using the CP2K modelling software augmented with homemade subroutines to impose crystal symmetry.In any calculation of elastic properties if crystal symmetry is not preserved then stiffness and compliance tensors will be incorrect.To evaluate the mechanical stability of the crystal we can apply the Born-Huang elastic stability criteria for the appropriate crystal system.As an example for the cubic crystal system, it is required that c 11 − c 12 > 0, c 11 + 2c 12 > 0, c 44 > 0. If one or more of these criteria is violated, one or more of the elastic tensor eigenvalues is negative and the crystal is mechanically unstable.
Fig. 2 shows the high quantitative accuracy that can be obtained using the above PBE-only methodology for individual stiffness tensor components and derived Young's modulus values for a small sample of the different classes of crystal that we have studied.As the piezoelectric response is inversely proportional to the elastic stiffness if the material is predicted to be more flexible than it is this will lead to an overestimation of the electromechanical coupling.Full mechanical testing of piezoelectric materials is always recommended as properties such as fracture limit and hardness ultimately determine the specific applications and environments the material can be used in.
Pei & Zeng 106 computed structural and elastic properties of nine phases of piezoelectric polyvinylidene fluoride (PVDF) crystals using DFT with and without a variety of dispersion corrections.In addition to the four known crystalline forms the mechanical properties of five theoretically predicted crystalline forms of PVDF were also investigated.The DFT/ PBE calculations show that the cell parameters of four known crystalline forms are in good agreement with experiment.However, they identified that including empirical van der Waals corrections, specifically the Grimme-D2 method, led to a large error in the calculated unit cell lattice parameters.By comparing conventional PBE (without dispersion corrections) and DFT-D2 calculations the authors could highlight that the PBE method provides a better description of the structural and mechanic properties of PVDF crystals.

Complex hybrid systems: limitations of DFT
Kosa et al. recently demonstrated a combined computational and experimental approach in studying the elastic anisotropy of a zinc phosphate phosphonoacetate (ZnPA) 107 metalorganic framework (MOF).This 3D framework is composed of Zn-O-Zn layers that are connected by phosphate groups bridging the ZnO 4 tetrahedra.This in turn sections off small channels that run parallel to the crystallographic a-axis.DFT calculations rationalise that the pore morphology contributes towards the elastic anisotropy.An efficient computational scheme Again a simple and efficient GGA calculation scheme was used, in this case to predict the Young's modulus and the Poisson's ratio (n) along the primary crystallographic axes.Notably, they avoided a finite difference methodology in order to circumvent computing the full elastic stiffness tensor (c ij ), which can be a computationally heavy approach for low-symmetrical crystal systems. 38,108The relative stiffness of the different crystal facets was found to be in reasonably good agreement with experiments.In terms of their absolute values, however, the calculated moduli are consistently higher, by as much as 25% to 40%, than those determined experimentally.Given that the hybrid compound considered here was not moisture sensitive the Young's moduli measured by nanoindentation are expected to be reliable, 109 as in the authors' previous work.This study pointed out that the higher stiffness calculated by DFT could be due both from the derivation scheme and deficiencies associated with the electronic structure methods used.They discuss how theoretical calculations at 0 K are expected to overestimate the elastic stiffness when compared to experiments performed at ambient 300 K.In relation to the electronic structure method, the authors note that an incorrect description of the position of transition metal d-states can cause overhybridization of the Zn-O bonds within GGA because of a self-interaction error (SIE).This can lead to artificial stiffening in the zinc-oxygen cores of the ZnPA framework and thus an overestimation of mechanical strength.They emphasise that due to the current limitations of DFT, elastic property predictions for complex hybrid systems are potentially sensitive to the choice of, among others, the exchange-correlation functional and the pseudopotential.
This journal is © The Royal Society of Chemistry 2021 Choosing your functional for high accuracy elastic constants Dengetal et al. highlight the benefits of testing multiple functionals to obtain accurate bulk moduli predictions.The elastic moduli from all functionals (PBE, PBEsol, optB88-vdW) are significantly higher than the values reported experimentally, with PBEsol deviating the least from experimental measurements.The authors attribute this difference to many factors.For example, the samples for characterization in the reported experiments were polycrystalline with both finite grain size and porosity.They are modelled however, as infinite single crystals in planewave DFT.The elastic moduli of samples with different grain sizes or porosities can vary, which can also lead to discrepancies in experimental values.The authors also derive Pugh's ratio, G/B, which is commonly used to evaluate the brittleness of materials, with a larger G/B indicating that the material is more brittle.
Another excellent screening of DFT functionals for elastic constant prediction was carried out by Rego & de Koning in their recent study on hexagonal proton-disordered ice using the Quantum ESPRESSO software package. 110They evaluated nine different exchange-correlation functionals, four of which include long-range dispersion interactions through the non-local van der Waals (vdW) approach.While we have observed that dispersion corrections can over-bind crystals and induce polymorphic transitions in small crystals, 90,106 they are known to play an important role in the condensed phases of water.The authors utilise the well-established energy-strain approach to calculate the elastic constants, in which one exploits the relation between the energy of a crystal and its state of deformation, using the equation With the exception of the inferior metaGGA SCAN functional all functionals predict an excessively stiff response to tensile/ compressive distortions, as well as shear deformations along the basal plane.These overestimates correlate with lattice constants and other intermolecular distances that are underestimated, and intramolecular separations that are too large.The inclusion of non-local van der Waals corrections generally improves these structural parameters and softens the elastic response functions (Fig. 3).

DFT screening of mechanical properties and database generation
Efforts aimed at developing databases of elastic moduli from first-principles computational methods have been undertaken in numerous studies. 112,113These computational approaches are advantageous as all of the data can be derived in a consistent manner, allowing for unambiguous comparisons across many classes of material.De Jong and colleagues have expanded on this approach, producing the largest database of calculated elastic properties of crystalline inorganic compounds to date, ranging from metals and metallic compounds to semiconductors and insulators.The calculations formed part of the larger high-throughput effort, 114 undertaken by the Materials Project. 115Using DFT the calculated elastic constants were consistently shown to be within 15% of experimental values, which represents a smaller scatter than that observed in experimental measurements in some cases (Fig. 4).Pearson (r) and Spearman (ρ) coefficients indicate that the calculations performed in this work yield elastic properties that show an excellent correlation with experimental values, making the database useful for screening materials with properties based on elastic tensors.
For this study the elastic constants were calculated using a stress-strain methodology.Starting from a relaxed structure Fig. 3 Relative difference of DFT values with respect to experimental data for the structural parameters of ice. 110,111Reproduced with permission from AIP Publishing.Fig. 4 Comparison of experimental and calculated bulk moduli for a selected set of systems, with calculated Pearson correlation coefficient r and Spearman correlation coefficient ρ reported. 19or each compound, a set of distorted structures is generated.For each of the applied strains ∈ ij , the full stress tensor is obtained from a DFT calculation in which ionic positions are relaxed.One row (or equivalently, column) of the elastic matrix is obtained from a linear fit of the calculated stresses over the range of imposed strains.Repeating this procedure for each of the 6 independent strain components, all elements of the elastic modulus tensor can be calculated.The result is a calculated set of c ij values that can be used to calculate properties such as the bulk modulus K and the shear modulus G as in previously mentioned studies.

Mechanical properties of crystal polymorphs
DFT is widely used to study crystal polymorphism, most commonly through the calculation of the ground state energy difference between polymorphic crystal structures, 116 or a simpler analysis of intermolecular interactions (Fig. 5). 117rom this, the relative stability of polymorphs can be determined, though high accuracy functionals and dispersion corrections are required for this. 1189][120][121] As with any crystal, DFT-calculated physical properties, be they mechanical, 122 dielectric, 116 or physiochemical 123 can also be used to distinguish and rationalize polymorphic behaviour.Since both polymorphism and distinct mechanical properties arise from differences in crystal packing of molecules (in the case of packing polymorphism), polymorphs can be identified by their unique mechanical properties.In addition to colour or morphology, mechanical properties can also be used for the rapid screening and sorting of concomitant polymorphic forms. 47This is especially useful for structurally similar polymorphic forms, such as the aspirin polytypes.The structural differences in aspirin I and aspirin II are not obvious.They both constitute hydrogen-bonded dimers ordered in layers, with the layers arranged differently relative to each other. 124Reilly and Tkatchenko have shown that their Young's moduli exhibit distinct anisotropies. 125In their work they use dispersion-inclusive DFT, by incorporating manybody van der Waals (vdW) interactions to calculate the elastic tensor of the two polymorphic forms.Their elastic property analysis shows that both forms are expected to be mechanically stable under compression and shearing at room temperatureshedding light on the widely-discussed mechanical stability of form II. 38,126 The work of Reilly and Tkatchenko shows that extensive characterization of polymorphic forms can be achieved from the determination of elastic tensor using DFT and other in silico methods.A strong interconnect between molecular modelling and experiment is preferable, as even with structure-property correlations, a polycrystalline lens is needed to fully understand properties like plasticity and tabletability.The tabletability of a polycrystalline material is influenced by factors such as the relative movements of grains and grain size, which are difficult to extrapolate from single crystal DFT.Karki et al. have successfully used the eigenvalues of the compliance tensor to rationalize the differences in mechanical behavior of paracetamol polymorphic formswhich show significant differences in structural featuresas well as paracetamol form II cocrystals. 48They used the value of the highest compliance eigenvalue to draw conclusions about the compliance of a crystal and the relative strengths of shear planes.From that, they establish that form II shows a higher compliance eigenvalue compared to form I, because it is compliant to shearing, which leads to plastic deformation during tableting.This result is consistent with the layered structure of form II which grants it preferable compaction properties. 127Despite these milestones of linking mechanical behavior to structure using DFT, crystal engineering of pharmaceuticals with desired mechanical behavior is still to be achieved.With the maturity of DFT, crystal structure prediction, and ever-increasing computational power, there is continued opportunity for new insights.

The right direction: the importance of elastic anisotropy
The mechanical properties of crystals are naturally anisotropic, with different stress-strain responses depending on the crystallographic axis of interest.Azuri et al., in their investigation of the high mechanical strength of amino acid crystals, demonstrated that experimental measurements perpendicular to different crystal planes yield unique Young's modulus values. 128The mechanical properties of cellulose nanocrystals (CNCs) have been difficult to experimentally characterize owing largely to extreme anisotropy and uncertainties about the structure of these materials.0][131] Three dimensional surfaces, which are colour contours showing the crystallographic dependence of the Young's modulus and Poisson's ratio, were computed by Dri et al. to examine the extreme anisotropy of these important elastic properties using DFT. 132A clear correlation between the stiffness of the crystal and the different deformation mechanisms was noted.The largest Young's modulus (206 GPa) was found to be aligned with the crystallographic c-axis where covalent bonds drive the mechanical response of the crystal.Perpendicular to the cellulose chain axis, the b-direction shows the next greatest value for the Young modulus (98 GPa), explained by the presence of the hydrogen bond network linking the cellulose chains.Finally, a Young modulus value of only 19 GPa was computed along the direction perpendicular to the previous two, where weak vdW interactions play a dominant role in the mechanical response of the material.The 0 K calculations, carried out with dispersion-corrected DFT in VASP predicted a transverse Young modulus for crystalline cellulose in the range between 13 and 98 GPa, in good agreement with the reported range of experimental results.
The elastic stiffness tensor predicted by some DFT methods is calculated as a 6 × 6 matrix that naturally describes the elastic anisotropy of a crystalline system.However, it is only recently that elastic anisotropy has begun to be considered in experimental measurement of mechanical properties.Mishra and co-workers systematically examined the mechanical properties of dimorphic forms, forms I and II, of a 1 : 1 caffeine-glutaric acid cocrystal on multiple faces.Here nanoindentation was used to fully understand the co-crystal mechanical anisotropy and mechanical stability under an applied load. 133The higher hardness and elastic modulus of stable form II was rationalized on the basis of its corrugated layers, higher interlayer energy, lower interlayer separation, and the presence of more intermolecular interactions in the crystal structure compared to metastable form I. The results show that mechanical anisotropy in both polymorphs arises due to the difference in orientation of the identical 2D structural features, namely, the number of possible slip systems and the strength of the intermolecular interactions with respect to the indentation direction.It is hoped that studies like these will influence future experimental investigations, where directional elastic stiffness can be correlated with DFT-Fig.6 Three-dimensional plot in GPa showing the magnitude and anisotropy of the material's Young modulus a. experimental measurement 137 for deuterated ammonia (ND 3 ) b. DFT-calculated Young's modulus with no dispersion corrections (PBE), and many-body dispersion corrections (PBE + MBD).c. DFT-calculated Young's modulus using the same two methods but with the inclusion of the quasi-harmonic approximation (QHA) to simulate the Young's modulus at a temperature of 194 K. Adapted from ref. 134 with permission from Wiley.
predicted tensors to rationalise mechanical properties from the nanoscale up.

Looking forward: room temperature property predictions
Current standard ab initio methods describe ground-state properties at zero temperature and pressure.One implication of this that the theoretical structures do not expand with increasing temperature and the elastic constants do not account for zero-point or temperature effects.A quick fix for this is to compare the ground-state properties determined from theoretical structures to those of experimental XRD structures grown and characterized at the lowest-possible temperature, where temperature effects are dampened.However, even with this work-around, elastic constants can still be overestimated by up to 40%. 134hermal contributions to the elastic constants can be incorporated by treating the lattice dynamics of a crystal structure within the quasi-harmonic approximation (QHA), as opposed to the more commonly adopted harmonic approximation (HA) approach. 135In the QHA, the lattice dynamics of the structure are modelled within the HA at several unit-cell volumes, therefore incorporating the volume dependence missing in the pure HA.It follows then that with the volume dependence, elastic constants now depend on temperature.The incorporation of temperature effects into elastic constant predictions from QHA calculations accounts for approximately 30% of the disagreement observed 134,136 compared to elastic constants at 0 K (Fig. 6).
Treating lattice dynamics in this way is known to increase computational expense and effort when considering complex systems, and even more so with the inclusion of dispersion corrections and many-body effects.Lower-level densityfunctional based methods can be considered in efforts to reduce computational expense. 138Density-functional tight binding (DFT-B) is one such method, 139 being an approximate treatment of the Kohn-Sham DFT formalism with less empirical parameters compared to classical force fields.It therefore lies between ab initio methods and classical force fields in terms of time scales and attainable system sizes.Despite its computational efficiency, DFT-B is still to be fully explored as a tool for predicting the elastic constants of materials.

Conclusions
In this Highlight we have introduced the theory of DFT, and multiple methodologies and software packages that can be used to calculate the elastic properties of periodic systems.From these calculations further mechanical properties can be derived, and if desired fed into FEA models for structural and behavioural analysis.While all DFT methodologies overestimate elastic constants to some extent, many combinations of exchange-correlation functional and dispersion corrections match well with experimental values.
Extensive benchmarking is recommended for any DFT investigation into the mechanical response of materials, not only to ensure accuracy but to confirm mechanical stability and correct imposition of symmetry.It is envisioned that the accuracy of DFT predictions of elastic constants will continue to improve with advances in high-performance computing power, as well as the incorporation of more efficient manybody interaction schemes with quasi-harmonic approximations in order to overcome the negative effects of calculations carried out at absolute zero.

Fig. 1
Fig. 1 Schematic of the nanoindentation process and measurement (a) diagram showing the working principles of indenting the sample via loading and unloading, (b) a corresponding load-displacement curve showing the effect of the loading and unloading process.S is the contact stiffness of unloading.Reproduced with permission from Wiley.46

Fig. 5
Fig. 5 DFT-calculated intermolecular interactions that contribute to the Hirshfeld surfaces of polymorphic form-II (top) and form-I (bottom) of a 1D Ni(II) polymer with iminodiacetic acid (IDA), namely diaquaiminodiacetatonickel(II).Reproduced from ref. 117 with permission from the Royal Society of Chemistry.