Thermoelasticity of Flexible Organic Crystals from Quasi-harmonic Lattice Dynamics: The Case of Copper(II) Acetylacetonate.

A computationally affordable approach, based on quasi-harmonic lattice dynamics, is presented for the quantum-mechanical calculation of thermoelastic moduli of flexible, stimuli-responsive, organic crystals. The methodology relies on the simultaneous description of structural changes induced by thermal expansion and strain. The complete thermoelastic response of the mechanically flexible metal-organic copper(II) acetylacetonate crystal is determined and discussed in the temperature range 0-300 K. The elastic moduli do not just shrink with temperature but they do so anisotropically. The present results clearly indicate the need for an explicit account of thermal effects in the simulation of mechanical properties of elastically flexible organic materials. Indeed, predictions from standard static calculations on this flexible metal-organic crystal are off by up to 100%.

Thermoelasticity represents the dependence of elastic mechanical properties of materials on temperature. In particular, the thermoelastic response of crystalline materials is described by the thermal dependence of all the isothermal or adiabatic elastic constants defining the fourth-rank elastic tensor, which provides the formal description of the anisotropic mechanical properties of the material in the elastic regime.[1] The accurate description of thermoelasticity is relevant to many areas of research including (i) geophysics, where the elastic properties of minerals at temperatures of the Earth mantle determine the velocity of propagation of seismic waves;[2−5] (ii) refractory materials, whose mechanical stiffness must not be deteriorated at high temperature;[6−9] (iii) pharmacology, where most potential drugs are synthesized in the form of molecular crystals, whose mechanical stability at room temperature is crucial for an effective tableting process;[10−15] (iv) catalysis, where the mechanical instability of porous frameworks poses serious limitations to their effective use, as for metal–organic frameworks.[16−19]

There is great interest in the possibility of describing the thermoelasticity of soft, flexible organic materials by computational approaches, which would ensure access to the complete characterization of their elastic anisotropic response at operando conditions. This stems from experimental evidence that the anisotropy of the elastic response of such systems can be largely affected by temperature,[20] which makes a purely static description (either scaled or unscaled by volume expansion) inadequate. Indeed, while for certain organic crystals the thermomechanical response nicely correlates with the sole thermal expansion,[20,21] this cannot be expected in general. It has recently been shown on single crystals of elastically flexible copper(II) acetylacetonate that the molecular motions responsible for thermal expansion and bending are different.[22,23] Therefore, a reliable description of the thermoelastic response of such materials requires the simultaneous description of structural changes induced by thermal expansion and strain. A successful prediction relies on the balanced description of the interplay of these effects. For these reasons, we focus on a quantum-mechanical approach, which ensures an accurate account of all the relevant chemical interactions.

While calculations based on the density functional theory (DFT) have proved reliable in the prediction of many static properties of materials in the last decades, the effective inclusion of thermal effects still represents a challenge to state-of-the-art DFT methodologies, particularly so when one needs to go beyond the usual harmonic approximation in the description of the lattice dynamics.[24−26] For instance, this is the case when (i) cubic interatomic force constants are needed to compute phonon lifetimes and the thermal lattice conductivity and (ii) the free energy dependence on volume and strain is needed to determine thermal expansion and thermoelasticity. While the former class of physical properties above requires the explicit evaluation of anharmonic terms of the nuclear potential,[27−33] the latter has often been tackled within the so-called quasi-harmonic approximation (QHA), as long as inorganic solids were considered.[34−39]

Owing to the progress made in the inclusion of dispersive interactions into the DFT and in the exploitation of parallel computing, the QHA could finally be effectively applied to organic molecular crystals and mixed organic–inorganic materials such as metal–organic frameworks to determine their thermal expansion.[40−48] Let us stress that while the evaluation of the thermal expansion of the system requires the calculation of the dependence of the free-energy on volume, the evaluation of the thermoelasticity involves the extra dependence of the free-energy on strain, which makes the corresponding lattice dynamical calculations much more demanding and challenging, particularly so for soft, flexible materials with low-frequency phonon modes.

In this Letter, we introduce an affordable quasi-harmonic computational protocol for the description of the thermoelasticity of organic crystals within the framework of dispersion-corrected DFT calculations. We apply our scheme to the investigation of the thermoelasticity of the stimuli-responsive, metal–organic copper(II) acetylacetonate crystal: a system that has recently attracted a lot of attention because of its unusual high flexibility[22] and its different structural mechanisms induced by temperature and strain.[23] Present quantum-mechanical calculations allow to report the first complete characterization of the 3D anisotropic elastic response of this system, as well as its thermal evolution. Inclusion of thermal effects up to room temperature results in dramatic changes of the statically computed values. Comparison with the experimentally measured values of the Young modulus along two crystallographic directions confirms the reliability of the description to an unexpected degree.

The full thermoelastic characterization of a crystal requires the determination of the fourth-rank thermoelastic tensor at each desired temperature T; that is, the full set of isothermal elastic constants (second free-energy density derivatives with respect to pairs of strain types):where V(T) is the volume of the equilibrium structure at temperature T, F is the Helmholtz free energy, η is one of the six independent components of the strain tensor η, and v, u = 1, ..., 6 are Voigt indices.[49] The calculation of thermoelastic coefficients from eq represents a formidable computational task for systems with more than a few atoms per cell and/or belonging to low-symmetry space groups.[37] This is because it involves the calculation of the free energy dependence on all lattice parameters to get the anisotropic thermal expansion and the second derivatives with respect to strain.[50−53] As such, the use of eq is unpractical for soft (metal-)organic crystals. However, based on experience gained in previous investigations of the anisotropic thermal expansion of soft molecular crystals,[40,43,44] we know that, at variance with strongly bound inorganic crystals, the dependence on the lattice parameters of the free energy, F(η; T), can be safely substituted by that of the static energy, E(η; T0), for weakly bound organic crystals. Therefore, here, we introduce a different strategy where thermoelastic constants are obtained fromwhere T0 is the absolute zero. While eq still relies on the quasi-harmonic determination of V(T), now the second energy derivatives with respect to strain are evaluated from the static internal energy E and not from the free energy F, and this remarkably simplifies the corresponding calculations. See the Supporting Information for a more detailed description of the algorithm. We stress that, crucially, the scheme represented by eq , while neglecting vanishingly small terms in weakly bound crystals, still couples thermal expansion with strain.

Calculations are performed with the CRYSTAL program for quantum-mechanical simulations of the condensed matter.[54,55] The hybrid PBE0 exchange-correlation functional is used[56] as corrected for missing dispersive interactions with Grimme’s D3 scheme:[57,58] PBE0-D3. A basis set of triple-ζ quality plus polarization is used.[59] Reciprocal space is sampled on a Monkhorst–Pack mesh with a shrinking factor of 2.

The copper(II) acetylacetonate molecules are almost planar and crystallize in a monoclinic lattice (space group P21/n) that is shown in Figure . In particular, Figure b shows the structure of the crystal in the ac crystallographic plane (XZ Cartesian plane). Green dashed arrows identify the two crystallographic directions along which the Young modulus has been measured experimentally by Worthy et al.[22] The two directions are [101] and [101̅] and correspond to the two directions in the plane connecting the Cu centers of the molecules along the two molecular chains.

Figure 1

(a) Atomic structure of copper(II) acetylacetonate crystals (belonging to the P21/n space group); (b) atomic structure of the crystal in the ac crystallographic plane (XZ Cartesian plane). The two directions, [101] and [101̅], along which the elastic Young modulus of the system has been experimentally determined are identified by green arrows.

The first step of our methodology consists of a standard QHA calculation to determine the thermal expansion of the system. Harmonic vibration frequencies have been computed at four volumes, from a −2.5% compression to a +5% expansion relative to the static optimized equilibrium volume. Harmonic frequencies as a function of volume are then fitted to a quadratic polynomial function, and the fitting coefficients are used to set up the canonical vibrational partition function and compute the Helmholtz free energy F(V; T) on a dense grid of volumes at each desired temperature. The V(T) relation is obtained by minimizing the free energy with respect to volume at several temperatures.[60−63] The computed volumetric thermal expansion V(T) is reported in Figure (black line) where it is compared with available experimental data (red squares) from ref (22). The agreement between the two sets is remarkable, which proves that the QHA provides a reliable description of the thermal expansion in this case. The metal–organic molecular crystal of copper(II) acetylacetonate is found to expand by 1.7% just by inclusion of zero-point motion effects at the absolute zero, which is followed by a further thermal expansion of about 4% when temperature is raised from 0 to 310 K. The thermal expansion along individual lattice parameters is reported in the Supporting Information.

Figure 2

Volumetric thermal expansion of copper(II) acetylacetonate crystals in the 0–310 K temperature range. Available experimental data (red squares) are from ref (22). The black line is the expansion obtained from our QHA calculations.

Starting from the thermal expansion determined by the QHA, isothermal elastic constants have been computed through eq at four different temperatures: 0, 100, 200, and 300 K. At the absolute zero, two sets of elastic constants have been evaluated, one corresponding to the equilibrium volume obtained from the standard static geometry optimization and one corresponding to the equilibrium volume at 0 K upon inclusion of zero-point motion effects. The explicit dependence on temperature of the full set of 13 independent elastic constants can be found in the Supporting Information, where it is shown that very regular trends are obtained for all elastic coefficients, including small ones, Once the full set of elastic constants is obtained, a variety of mechanical features of the material can be derived such as the bulk modulus, the shear modulus, the Young modulus, Poisson’s ratio, the velocity of propagation of directional elastic waves, etc. In particular, the directional Young modulus can be determined, which measures the stiffness of the material along any direction in space.[49,64] The thermal evolution of the Young modulus of copper(II) acetylacetonate crystals is shown in Figure where it is reported in 2D maps in three Cartesian planes (XY, YZ, and XZ). The following can be observed: (i) as expected, the mechanical stiffness decreases with temperature (i.e., the Young modulus decreases); (ii) the Young modulus does not shrink isotropically with temperature and rather the anisotropy of the mechanical response significantly evolves as temperature changes. For instance, let us have a closer look at the spatial dependence of the Young modulus in the YZ plane. While the Young modulus almost does not change along the Y axis (with values close to 6 GPa at 0 and 300 K), its value passes from 27 to 11 GPa in the same temperature range in the diagonal YZ direction. A more quantitative analysis on the thermal evolution of the computed Young modulus is provided in the table in the bottom right corner of Figure , where the two crystallographic directions in the ab plane, [101] and [101̅], are considered, which were graphically marked by the green arrows in Figure b. These are the only two directions along which the Young modulus of copper(II) acetylacetonate crystals has been experimentally measured.[22] Also in this case, along one direction, [101̅], the Young modulus shows little thermal evolution, while along the second one, [101], it changes from 25.3 GPa at 0 K in the static limit to 11.3 GPa at 300 K. In the Supporting Information, we report the thermoelastic response of the system as obtained from a nonhybrid functional of the DFT, PBE-D3: while the absolute values of the computed mechanical properties are systematically slightly lower than those from PBE0-D3 (for reasons also discussed in the Supporting Information), thermal trends are found to be very consistent with those discussed above, which confirms the reliability of the approach. The comparison with the experimental values shows that it is absolutely mandatory to take into account thermal effects to accurately reproduce the room temperature mechanical features of this class of soft metal–organic molecular crystals and, at the same time, that the simplified quasi-harmonic scheme introduced here provides a reliable and yet feasible way to do so.

Figure 3

2D plots of the Young modulus (in GPa) of copper(II) acetylacetonate crystals in three planes (XY, YZ, and XZ) and at four temperatures (0, 100, 200, and 300 K). The directional Young modulus at 0 K is computed both without and with inclusion of the zero-point energy (ZPE). The XZ plane is the one displayed in Figure b. The table in the bottom-right corner of the figure reports the thermal evolution of the computed Young modulus along the two crystallographic directions, [101] and [101̅], probed in the experiments in ref (22).

The evolution with temperature of the mechanical stiffness of copper(II) acetylacetonate crystals is shown even more explicitly in Figure , where a 3D representation of the spatial dependence of the Young modulus is reported. The figure clearly shows how, still down to the absolute zero of temperature, the large expansion of about 1.7% due to the lattice-dynamical zero-point motions produces a sizable modification to the elastic mechanical response of the system. Both the magnitude of the stiffness and the anisotropy of its spatial distribution are significantly affected by the expansion. As temperature increases from 0 to 300 K a further volume expansion by about 4% is observed, which is again reflected in the thermoelastic features of the crystals.

Figure 4

3D plots of the spatial distribution of the Young modulus of copper(II) acetylacetonate crystals as a function of temperature. Data are in GPa.

In conclusion, we have presented a computationally affordable scheme to compute the anisotropic thermoelastic response of soft, flexible (metal-)organic crystals based on the quantum-mechanical description of their quasi-harmonic lattice dynamics. This scheme relies on the strain dependence of the static internal energy and, crucially, couples thermal expansion with strain. We have applied the scheme to stimuli-responsive, elastically flexible copper(II) acetylacetonate monoclinic crystals and reported the first characterization of its anisotropic mechanical behavior in the temperature range 0–300 K.