Elastic and piezoelectric properties of β-glycine – a quantum crystallography view on intermolecular interactions and a high-pressure phase transition

The effect of hydro­static compression on the elastic and electronic properties of β-glycine was studied using a quantum crystallography approach. The interrelations between the changes in the microscopic quantum pressure in the electronic continuum, macroscopic compressibility and piezoelectricity were considered. The geometries and energies of hydrogen bonds in the crystal structure of β-glycine were considered as functions of pressure before and after a phase transition into the β′-phase in relation to the mechanism of this phase transition.

1. Introduction

The combination of quantum chemistry calculations and diffraction experiments provides a synergistic effect. It allows one to refine and improve structural data, to study electron density distribution, to interpret, or even to predict physical and chemical properties. This approach is known as quantum crystallography (Grabowsky et al., 2017; Tsirelson, 2018; Genoni & Macchi, 2020; Grabowsky et al., 2020; Tsirelson & Stash, 2020; Macchi, 2022; Matta et al., 2023). One of the important aims of quantum crystallography is to elucidate the mechanical characteristics of crystals, relating microscopic properties (distortions of chemical bonds, changes in the intermolecular contacts under high pressure) and macroscopic properties, such as structural response to hydro­static pressure, compressibility and piezoelectric properties (Zhurova et al., 2006; Coudert & Fuchs, 2016; Riffet et al., 2017; Tsirelson et al., 2019; Evarestov & Kuzmin, 2020; Bartashevich et al., 2020; Korabel'nikov & Zhuravlev, 2020; Mishra & Tewari, 2020; Bartashevich et al., 2021; Feng et al., 2021; Matveychuk et al., 2021; Bogdanov et al., 2022; Gajda et al., 2022; Zhuravlev & Korabel'nikov, 2022; Stachowicz et al., 2023). Revealing relations between structural changes and mechanical properties of crystals helps to understand the nature of the piezoelectric effect and to design new piezoelectric materials (Guerin et al., 2019; Vijayakanth et al., 2022; Ivanov et al., 2023), as well as to rationalize thermo- and photomechanical effects with applications for design of materials and devices (Naumov et al., 2013, 2015; Koshima et al., 2021; Karothu et al., 2022; Awad et al., 2023).

Different types of noncovalent intermolecular interactions, especially hydrogen bonds (Dougherty, 1998; Desiraju & Steiner, 1997), play an important role in the formation of molecular crystal structures and their evolution under temperature and pressure variations (Katrusiak & Dauter, 1996; Katrusiak & Szafrański, 1996; Katrusiak, 2003; Resnati et al., 2015). Hydrogen bonds in organic crystals have been extensively studied under external hydro­static pressure (Katrusiak & Szafrański, 1996; Boldyreva, 2003, 2004; Gatti & Macchi, 2012; Isono et al., 2013; Zakharov & Boldyreva, 2019; Bartashevich et al., 2020). At the same time, attempts to relate macroscopic mechanical and piezoelectric properties to the distortions of crystal structure on hydro­static compression are rare (Bartashevich et al., 2020; Boldyreva, 2021; Zhu et al., 2022; Qiao et al., 2023; Khainovsky et al., 2023).

Orbital-free quantum crystallography approaches provide new opportunities to characterize crystals under compression (Zhurova et al., 2006; Casati et al., 2016, 2017; Coudert & Fuchs, 2016; Tsirelson et al., 2016; Riffet et al., 2017; Tsirelson et al., 2019; Tsirelson & Stash, 2020; Evarestov & Kuzmin, 2020; Bartashevich et al., 2020; Korabel'nikov & Zhuravlev, 2020; Mishra & Tewari, 2020; Bartashevich et al., 2021; Gajda et al., 2022; Zhuravlev & Korabel'nikov, 2022; Stachowicz et al., 2023). As in the classical theory of elasticity (Landau et al., 2009), the compression and stretching of the electronic continuum in molecules and crystals can be described by the density of the stress tensor, σ(r) (Rogers & Rappe, 2002; Tao et al., 2008; Tsirelson et al., 2016). In the framework of density functional theory (DFT) (Dreizler & Gross, 1990; Tao et al., 2008)

where σ_kin(r) is the quantum contribution of electron kinetic energy, σ_xc(r) is the exchange-correlation term, and σ_es(r) is the classical part resulting from electrostatic interactions of all electrons and nuclei [Maxwell's tensor (Craggs, 1962)]. In the equilibrium state of a many-electron multinuclear system, the forces of internal stresses, determined by the divergence of the density of the tensor σ(r), balance the external forces. Details of determining the components of the stress tensor of the electronic continuum are given in Appendix A.

According to Bader (1990), kinetic and potential electron energy densities characterize the internal quantum pressure of the inhomogeneous electron cloud spanning the nuclei positions, . The changes in the electron continuum in the intermolecular space under pressure characterize the changes of mechanical properties of solids. A local increase in p(r) corresponds to an increase in internal electron energy per unit volume, while the decrease in electron energy density corresponds to the compensation of internal compressions in a system. The effects of external pressure are additive to the internal changes. It is appropriate to consider the local pressure at the critical points in the electron density to receive a concise picture of the changes in the electronic structure and bonding in a crystal under pressure (Zhurova et al., 2006; Dorbane et al., 2019; Bendjemai et al., 2020; Tsirelson & Ozerov, 2020).

Modeling of crystal structures under pressure by DFT methods has been performed for many compounds and functional materials: organic aromatic compounds such as syn-1,6:8,13-bis­carbonyl­[14]annulene (Casati et al., 2016, 2017); perovskites (Tariq et al., 2015; Dar et al., 2017; Coduri et al., 2019; Ciupa-Litwa et al., 2020), kaolinites (Richard & Rendtorff, 2022; Richard & Rendtorff, 2022), carbonates (Zhuravlev & Atuchin, 2021; Zhuravlev & Korabel'nikov, 2022) and other inorganic compounds (Faridi et al., 2018; Nazir et al., 2018; Yaseen et al., 2021); grossular (Gajda et al., 2020), ice (Chodkiewicz et al., 2022), zeolites (Stachowicz et al., 2023) and also molecular crystals (Schatschneider et al., 2013; Liu et al., 2014; Moellmann & Grimme, 2014; Matveychuk et al., 2021). Such studies make it possible to predict mechanical, thermodynamic and optoelectronic properties, as well as phase transitions of various types. Piezoelectric properties of crystals under pressure were studied for some inorganic crystals, such as wurtzite (Marana et al., 2017; Almaghbash & Arbouche, 2021), potassium niobate (Faridi et al., 2018; Stash et al., 2021) and others (Daoud, 2019; Mubarak & Tariq, 2021). However, to the best of our knowledge, molecular crystals were not considered from this point of view yet. In this context, crystalline glycine, ⁺NH₃-CH₂-COO⁻ (Boldyreva, 2021), is one of most interesting compounds to be studied using a quantum crystallography approach. To date, it is the only amino acid that can form several polymorphs on crystallization and that different polymorphs can coexist and be stored under ambient conditions for a very long time. Two of the three polymorphs (β- and γ-forms) are piezoelectric (Iitaka, 1960, 1961) and are used in various applications, including biomedical ones (Chen et al., 2003; Heredia et al., 2012; Guerin et al., 2018; Sui et al., 2022; Ivanov et al., 2023). Glycine was the first amino acid, for which the effect of hydro­static pressure on the crystal structure was studied. It was shown that the effect of pressure depends on the starting polymorph. The α-polymorph does not undergo phase transitions at least up to 50 GPa on hydro­static compression (Murli et al., 2003; Sharma et al., 2012; Shinozaki et al., 2018; Hinton et al., 2019). At the same time, hydro­static compression of two other polymorphs with piezoelectric properties, β- and γ-forms, gives different high-pressure phases, depending on which polymorph was taken initially. Crystals of the γ-polymorph are transformed into another high-pressure δ-form (Boldyreva et al., 2004, 2005; Dawson et al., 2005; Goryainov et al., 2006). On decompression, the high-pressure δ-form first transforms into ζ-form at 0.62 GPa (Goryainov et al., 2006; Bull et al., 2017) and into the α-polymorph only at ambient pressure on contact with the atmosphere over time after opening the cell (Goryainov et al., 2005). In this case, the structural rearrangement is very large and a single crystal is fragmented into powder. Crystals of β-glycine are transformed reversibly into a structurally related β′-form at 0.76 GPa and single crystals are not destroyed in this case (Dawson et al., 2005; Goryainov et al., 2005; Tumanov et al., 2008). The electronic structure of crystalline glycine has been repeatedly analyzed using quantum chemical calculations at ambient pressure for the α-modification and, much less, for the β- and γ-modifications (Behzadi et al., 2007; Flores et al., 2008; Stievano et al., 2010; Lund et al., 2015; Marom et al., 2013; Moggach et al., 2015; Seyedhosseini et al., 2015; Rodríguez et al., 2019; Guerra et al., 2020; Xavier et al., 2020). Attempts to describe the structure of high-pressure phases of glycine using quantum chemical calculations were also undertaken (Chisholm et al., 2005; Guerin et al., 2018; Szeleszczuk et al., 2018; Hinton et al., 2019; Mei & Luo, 2019). The piezoelectric properties of mixed crystals containing glycine as one of the components were also reported (Guerin et al., 2021). The calculations were aimed at estimating the differences in the relative thermodynamic stability of the polymorphs, to interpret their optoelectronic properties and to reproduce (or predict) the structure of the high-pressure equilibrium phase. To the best of our knowledge, the mechanical properties of glycine polymorphs have not been analyzed yet.

X-ray diffraction data collected at high pressure from a powder sample of a low-symmetry organic crystal with weak scattering atoms can provide the atomic coordinates, but are not sufficient to make unambiguous conclusions on the fine effects of the redistribution of the electron density, which are of utmost importance for phase transitions or piezoelectric response. In the present contribution we report, how DFT calculations, which are interpreted from a viewpoint of orbital-free quantum crystallography (Tsirelson & Stash, 2020, 2021) make it possible to follow the electron continuum evolution based on the experimental powder diffraction data collected from a relatively weak scatterer, glycine, years ago (Tumanov et al., 2008).

1.1. Calculation details

The DFT modeling of the electronic structure of the glycine β-polymorph and its changes on hydro­static compression was carried out using the Kohn–Sham method in the CRYSTAL17 software package (Dovesi et al., 2018). We took into account the periodicity of the electronic wavefunctions of the ground electronic state. Since the elastic properties of organic piezoelectrics are largely determined by hydrogen bonds, we chose the exchange-correlation functional PBE0 (Adamo & Barone, 1999) because makes it possible to model crystals and their elastic properties more precisely (Erba et al., 2013). The 6-31G(d,p) basis set from Pritchard et al. (2019) was used. External pressure corresponded to the points: 0, 0.2, 0.4, 0.7, 0.9 and 1.7 GPa. As starting unit-cell parameters and atomic coordinates, we used those obtained in an earlier variable-pressure powder diffraction experiment performed at room temperature at the ESRF synchrotron source (Tumanov et al., 2008). The experimental values of unit-cell parameters and atomic coordinates were important to have to establish a starting model corresponding to the real structure. At the same time, the experimental data obtained from a weak scattering powder sample in a diamond-anvil cell (DAC) years ago, before new detectors became available, needed to be additionally refined by DFT. Both unit-cell parameters and atomic coordinates were therefore optimized to avoid the risk of obtaining the structural model corresponding to a non-equilibrium point of the potential energy surface (PES). The hydrogen-atom positions were derived from theoretical calculations by CRYSTAL17.

The energy convergence criterion for geometry optimization was 10⁻¹⁰, gradient RMS < 0.0003, displacement RMS < 0.00012. The SHRINK factor, which determines the number k of points in reciprocal space in the Pack–Monkhorst scheme, in which the Kohn–Sham matrix is diagonalized, was equal to 8 8. The TOLINTEG parameter, which is responsible for the values of the overlap integrals, was set to 10 10 10 10 20, for ensuring a sufficiently high calculation accuracy (Pascale et al., 2022). All calculations were carried out either applying the dispersion correction D3 (Grimme et al., 2007) or without it. The IR vibrational frequencies were calculated and none of them were imaginary, confirming that energy minima were reached. Quantum-topological analysis of electron densities was performed using TOPOND and MultiWFN (Lu & Chen, 2012). Piezoelectric coefficients (both strain and stress ones) were obtained through the Berry phase approach and visualized used the MTEX Matlab package (Bachmann et al., 2010).

The internal quantum electron pressure at each point of the system was calculated. The pressure has the dimension of J m⁻³ and describes the self-organization of the quantum part of the electron energy density of the system in the physical space (rather than the electron density itself). In the absence of external fields, regions of positive pressure (compressed) are characteristic of atomic cores and intra­molecular interactions, while stretched (soft) regions of negative pressure are observed for intermolecular inter­actions. To find out the behavior of the internal electron pressure under external action, we approximated the electron kinetic energy density g(r) according to Kirzhnits (1957) and moved over to the field of orbital-free quantum crystallography (Tsirelson & Stash, 2020, 2021) using formula (10) in Appendix A.

Elastic macroscopic properties of crystals were characterized by elasticity descriptors, such as Young's modulus (E), shear modulus (G), hydro­static (H) and volumetric (K) compressibility (Reuss, 1929; Hill, 1952). To obtain the values and to determine the spatial anisotropy of the hydro­static compressibility H (the compliance of a crystal under an isotropic external stress), we calculated the elastic tensors for the low- and high-pressure phases of β-glycine (Erba, 2016) and analyzed them in the ELATE software package (Gaillac et al., 2016).

2. Results and discussion

2.1. Unit-cell parameters, volume and the geometry of hydrogen bonds

The values of the unit-cell parameters optimized during the quantum chemical calculations are given in Table S1 in supporting information. Since the calculations were performed for the structural model at 0 K, for a comparison with the calculated values at atmospheric pressure we used the values of the unit-cell parameters and volume measured in the temperature range 100–300 K (Boldyreva et al., 2003) extrapolated to 0 K. The crystal setting used by Boldyreva et al. (2003) was preliminary changed to the setting used by Tumanov et al. (2008) – the same as was used also in this work.

The best agreement between the computed and the experimental unit-cell parameters and volume was achieved without using the dispersion correction (Fig. 1). We note that application of dispersion corrections (Grimme et al., 2007) has already been doubted when studying the mechanical properties of crystals of various amino acids and their cocrystals (Kiely et al., 2021). We conclude that the PBE family of functionals without dispersion correction may more accurately describe the structural (and elastic properties) of glycine, since this correction tends to overestimate van der Waals inter­actions in the structures, where the long-range electrostatic interactions and various hydrogen bonds are important and may play a major role in determining structures and their response to pressure. This is in agreement with earlier findings (Azuri et al., 2015).

Figure 1 Relative changes in (a) a/a0, (b) b/b0, (c) c/c0, (d) β/β0 and (e) V/V0 as a function of external compression according to experimental (Tumanov et al., 2008) and calculated values with and without dispersion correction D3. Parameters with index 0 correspond to normal conditions. The gap corresponds to a 0.76 GPa phase transition and separates the low-pressure phase existence region from the high-pressure phase existence region, in accordance with experimental data (Tumanov et al., 2008). Connecting lines are calculated using B-spline fitting to facilitate the visualization of trends and are guides to the eye.

The largest linear strain in the structure of β-glycine is observed in the b direction and the smallest one in the a direction (Table S1 and Fig. 1). This agrees with the fact that molecular chains of zwitterions in the (ac) layer are parallel to the a axis and orthogonal to the b direction. During the phase transition, the c parameter doubles due to a change in the symmetry of the structure, although the changes in intermolecular distances in the c direction are small.

The values of covalent bond lengths and bond angles optimized by quantum chemical calculations agree with the previously obtained experimental values (Tumanov et al., 2008) (Table S2). As pressure increases to the point of polymorphic transition, the intramolecular geometry changes insignificantly, which is typical for hydro­static compression of molecular organic crystals. At a pressure of 0.76 GPa, a reversible first-order phase transition to the high-pressure phase β′ occurs (Dawson et al., 2005; Goryainov et al., 2005; Tumanov et al., 2008). Every second zwitterion rotates slightly, and, as a result, the asymmetric unit of the high-pressure β′-polymorph contains two independent zwitterions, so that the unit-cell c parameter doubles (Figs. 1 and 2). The intramolecular geometry of symmetrically independent zwitterions differs insignificantly, except for the angles characterizing the position of the amino group relative to the backbone of the zwitterion.

Figure 2 Fragments of β-glycine (on the left side) and β′-glycine (on the right side) crystal structures in two orientations. The blue and yellow colors correspond to two symmetrically independent zwitterions. Hydrogen bonds numbered I, II, III and IV correspond to those in the text, tables and other figures. Hydrogen bonds are shown by cyan dotted lines. Symmetry codes of acceptors in hydrogen bonds: (I) −x, y + ½, −z; (II) x, y, z − 1; (III) −x, y + ½, −z; (IV) −x, y + ½, −z + 1.

Based on interatomic distances, one could conclude that zwitterions in the crystal structure of β-glycine are linked by shorter hydrogen bonds N1—H3⋯O1 (I) (symmetry code −x, y + ½, −z) and N1—H4⋯O2 (II) (symmetry code x, y, z − 1) into two-dimensional layers in the (ac) plane, similar to those present in the α-polymorph, as well as in the high-pressure δ- and ζ-phases (Boldyreva, 2021). Interlayer hydrogen bonds N1—H5⋯O1 (III) and N1—H5⋯O2 (IV) are longer (symmetry codes are −x, y + ½, −z and −x, y + ½, −z + 1, respectively). On hydro­static compression these hydrogen bonds shorten and the distances between molecules in the layers decrease (Table S4).

To estimate the energy values of hydrogen bonds E_{H bond}, we used the local electron kinetic energy density g_b(r) at the critical points of hydrogen bonds, namely, the empirical relation E_{H bond} = 0.429g_b(r) (Espinosa et al., 1998; Vener et al., 2012, 2015). These values are listed in Table S3.

The analysis showed that the hydrogen bonds become stronger with increasing external hydro­static pressure: E_{H bond} increases and the distances between the proton and the acceptor of the energy value of hydrogen bonds decrease (Table S4).

The agreement between the calculated and the experimental values of distances between the non-hydrogen atoms was reasonable (Table S4); the main trends in the relative values of distances between non-hydrogen atoms were reproduced. Exact O⋯H distances could not be found by X-ray diffraction, so the calculated values are the only information we can rely upon in this case (Fig. 3).

Figure 3 Dependence of the energy value of hydrogen bonds, EH bond (a), quantum electronic pressure (QEP) (b) and bond distances (H⋯O) (c) at the hydrogen bond critical points in β- and β′-glycine on external compression. An increase in EH bond corresponds to a decrease in QEP. The gap corresponds to a 0.76 GPa phase transition and separates the region of existence of the low-pressure phase from the region of existence of the high-pressure phase. Symmetry codes of acceptors in hydrogen bonds: (I) −x, y + ½, −z; (II) x, y, z − 1; (III) −x, y + ½, −z; (IV) −x, y + ½, −z + 1. Connecting lines are made using B-spline fitting to facilitate the visualization of trends and are guides for the eye.

From the geometric and energy characteristics, it became clear that the relatively strong non-covalent bonds I and II are less compressible, than the bonds III and IV. With increasing pressure, the resistance to deformation along the molecular layers of glycine in the crystallographic direction a increases. Strain is the largest in the direction b, normal to the molecular layers. The trends in the changes in the intermolecular distances for hydrogen bonds with pressure are in agreement with experiment data (Fig. 4, Table S4).

Figure 4 Intermolecular distances of weak hydrogen bonds in β- and β′-glycine as a function of external compression. Connecting lines are calculated using B-spline fitting to facilitate the visualization of trends and are guides for the eye. Symmetry codes of acceptors in hydrogen bonds: (III) −x, y + ½, −z; (IV) −x, y + ½, −z + 1.

Due to rigidity of strong bonds I and II, the elastic and piezoelectric properties of β-glycine depend on the features of weak hydrogen bonds (III and IV). Therefore, we considered these interactions in more detail. For this purpose, various binding descriptors based on the electron density and its topology can be used; that makes it possible to study the space distribution of electrons and associate the changes in the electronic continuum with changes in the crystal structure (Bader, 1990).

2.2. Microscopic features of weak hydrogen bonds on hydro­static compression

We estimated the curvature of the electron density (ED) at the bond critical point (BCP) using the eigenvalues of the Hessian of the electron density, λ_i (Table S3), and the parameter associated with them – the ellipticity of the bond ɛ = (λ₁/ λ₂) − 1. The index ɛ at 0.0001 GPa for bonds III and IV is larger than for bonds I and II (Table S1). Weaker bonds are characterized by a more uneven distribution of ED normal to the bond lines. The values of the ellipticities of bonds I and II are close to each other.

Non-covalent interactions are characterized by positive values of the Laplacian at the BCP (Bader, 1990), see Table S3. The H5 atom forms a bifurcated hydrogen bond N1—H5⋯(O1/O2), which, for convenience of presentation, we divided into two components – bonds III and IV. This bifurcated bond is asymmetric (Table S3); with increasing pressure after the phase transition the values of the electron density in the two branches (III and IV) change (see Section 2.3).

The distribution of the Laplacian in the hydrogen bonds of β-glycine (Fig. S1) shows that in all cases the lone pair of the oxygen atom involved in the hydrogen bond points towards the H atom. After the phase transition this pattern does not change.

Evaluating the properties of the electronic continuum from the point of view of its response to external action, we found (Table S3) that the greatest compliance of the electronic continuum of both β- and β′-glycine phases (and hence of their structural micro-characteristics) takes place in the intermolecular space. At the same time, regions of relatively high positive pressure (atomic cores, covalent bonds, lone electron pairs) show the greatest resistance to external action. This is an expected behavior for a molecular crystal.

Less predictive a priori is the response to pressure of different hydrogen bonds in the structure. All hydrogen bonds of glycine show negative values of quantum pressure, QEP(r), at BCPs [Fig. 3(b) and Table S3], revealing soft regions of the electron continuum. The absolute values of QEP(r) for weak bonds III and IV are smaller than those for bonds I and II. As the external pressure increases, the internal pressure on these bonds slightly increases. The anisotropy of the resistance to external compression is related to the elastic properties of β-glycine. Strain along the molecular layers is smaller due to relatively strong hydrogen bonds and is facilitated normal to them due to weaker hydrogen bonds.

2.3. Changes in the asymmetry of weak hydrogen bonds

With increasing pressure, the ellipticity increases for bonds I, II and III, and decreases for bonds IV (Table S4). This indicates an increase in the asymmetry of the electron density distribution as a result of hydro­static compression for all bonds, except for IV. Geometry characteristics of weaker bonds III and IV (which are branches of a bifurcated hydrogen bond) show that one can observe a change in the asymmetry of this bifurcated hydrogen bond with increasing pressure: two unequal branches first interchange their lengths (the longer becomes shorter and vice versa), after that, at an even higher pressure, the branches become practically equivalent. During the phase transition, a jump in the ellipticity index for weaker bonds III and IV was observed, indicating an increase in the mutual influence of bonds III and IV. The same hydrogen atom H5 is involved in forming these bonds, and the calculated change in the bond length III during the phase transition Δr_III^{(0.7→0.9 GPa)} = 0.041 Å is comparable with the change in the bond length IV Δr_IV^{(0.7→0.9 GPa)} = 0.049 Å (Table S3). Thus, the calculated O1—H5—O2 angle increases from 90° to 92°. This indicates that the inter-electronic repulsion between bonds III and IV increases slightly. This may be related to an increase in the electron concentration normal to the connecting line and indicates a redistribution of electron density between bonds III and IV during the phase transition.

Thus, during the phase transition at 0.76 GPa, the branches III and IV of a bifurcated hydrogen bond `switch over', although the rest of the hydrogen bond network does not change significantly (Table S3). This phenomenon leads to a denser molecular packing. At the same time, more rigid hydrogen bonds I and II keep the structure from a complete rearrangement: the H⋯O distances in these bonds decrease continuously and monotonically even when passing through a phase transition. This keeps the single crystal from destruction and provides an easy return to the initial structure on decompression. Similar phenomena, when most of the hydrogen bonds are continuously compressed with increasing pressure, and only a few selected bonds switch over in a jump-wise manner providing a structural reorganization, were described in the literature earlier for reversible single-crystal-to-single-crystal high-pressure phase transitions (Zakharov & Boldyreva, 2013). With a further increase in pressure, branches III and IV of the bifurcation hydrogen bond become equal (see Section 2.4).

2.4. Equalization of branches of the bifurcated hydrogen bonds

As noted above, after the phase transition, with a further increase in the external pressure to 1.7 GPa, the quantum pressure in the BCP of weaker bonds III and IV becomes similar. This confirms the hypothesis that the weak bonds III and IV tend to equalize their characteristics under external hydro­static compression, which corresponds to the symmetrization of the N1—H5⋯(O1/O2) bifurcated hydrogen bond. This leads to a significant increase in the resistance of the structure to strain in the b direction after the phase transition, and may be one of the reasons for the stability of the β′-glycine phase with respect to new phase transitions, at least up to 7.6 GPa (Goryainov et al., 2005). This is much higher than the pressure at which the γ-polymorph of glycine, which is thermodynamically stable under ambient conditions, undergoes a phase transition to the δ-form (Boldyreva et al., 2005).

Quantum electronic pressure maps at 0.9 GPa and 1.7 GPa (Fig. 5) reflect clearly this change in the system of weaker bonds III and IV in the high-pressure phase: during the phase transition, the quantum pressure in the intermolecular space decreases, which corresponds to an increase in the concentration of electrons in both bonds. At the same time, for relatively strong hydrogen bonds I and II the intermolecular space is compressed only slightly as external pressure increases.

Figure 5 Maps of quantum electronic pressure in the structures of β-glycine and β′-glycine for weak hydrogen bonds in the corresponding planes at different pressures: (a) 0 GPa, (b) 0.9 GPa and (c) 1.7 GPa. Isoline interval is ± (2, 4, 8) × 10^n a.u. (−3 ≤ n ≤ 3).

All these observations agree with the geometry and energy analysis: under compression, the hydrogen bonds in β-glycine shorten, and this is accompanied by an increase in the bond energy E_b and a decrease in the quantum pressure.

Thus, in a β-glycine crystal, the main changes in the elastic properties under pressure are determined by the reorganization of the intermolecular space corresponding to van der Waals interactions and hydrogen bonds, and not by intramolecular bonds. This is what one can expect from a molecular crystal with low conformational flexibility at relatively low pressures not affecting the covalent bonds and angles.

The anisotropy of compressibility decreases clearly with increasing pressure after the phase transition. We associate this with the equalization of branches III and IV of the bifurcated hydrogen bonds, which are responsible for the maximum compressibility (Fig. S5).

The main components of the elasticity tensor (Table 1), its eigenvalues, are positive at all points of external pressure, which indicates that the crystal of β-glycine is mechanically stable. These values increase with increasing pressure, reflecting an increase in the resistance of the crystal to compression. The moduli of elasticity calculated according to Voigt (1887), Reuss (1929) and Hill (1952) are given in Table 1.

Table 1
Calculated components of the main diagonal of the elasticity tensor C of β-glycine

P (GPa)	0	0.2	0.4	0.7	0.9	1.7
C11	98.11	101.12	103.70	107.45	102.22	109.12
C22	28.77	30.14	31.53	33.78	25.96	30.55
C33	49.87	50.86	51.93	53.72	53.95	62.25
C44	6.64	7.06	7.54	8.13	8.97	9.67
C55	15.72	16.03	16.30	16.75	17.25	19.09
C66	8.63	8.93	9.28	9.67	11.64	12.26

The linear compressibility anisotropy index shows (Table S5) that crystals of β-glycine are rather poorly compressible on applying isotropic pressure. The compressibility of a crystal under external pressure can be conveniently illustrated by spatial visualization of the maximum and minimum values of the individual components of the linear (hydro­static) compressibility (Fig. 6). In β-glycine, the maximum hydro­static compressibility at 0.0001 GPa was observed normal to the molecular layers in the III and IV bonds directions. With an increase in pressure, the maximum hydro­static compressibility decreases, which is accompanied by an increase in the values of E_{H bond} when the corresponding bonds are compressed. Weaker hydrogen bonds contribute little to the resistance to strain in the unit-cell direction b (Table 1, Fig. 1). The direction of minimum hydro­static compressibility corresponds to the compression of the crystal in the plane of molecular layers (ac) at an angle of 25° relative to the C1—C2 bond. This can be explained by the fact that stronger hydrogen bonds I and II make a significant contribution to the resistance to strain along the molecular layers. It should be noted that the angle between H_min and bond II is 21°, whereas between H_min and bond I it is 86°. It can be assumed, that bond II makes a greater contribution to the resistance to these deformations. The contribution of individual bonds (I and II) to H_min is difficult to differentiate.

Figure 6 Visualization of the anisotropy of the compressibility under hydro­static compression for β-glycine. The yellow axis corresponds to the direction of minimum hydro­static compressibility, the cyan axis corresponds to maximum hydro­static compressibility. Pink, red, green and blue colors of dotted lines correspond to hydrogen bonds I, II, III and IV, respectively. Symmetrically non-equivalent glycine molecules are colored green and red. Symmetry codes of acceptors in hydrogen bonds: (I) −x, y + ½, −z; (II) x, y, z − 1, (III) −x, y + ½, −z; (IV) −x, y + ½, −z + 1.

After the phase transition, the hydro­static compressibility does not change its anisotropy, but changes in value; H_max also decreases, which is accompanied by gradual equalization of the branches III and IV of the bifurcated hydrogen bonds on further hydro­static compression of the β′ form above the phase transition pressure.

As long as macroscopic elastic properties are closely related to hydrogen-bond distances, it is natural to focus on the microscopic evolution of the electronic continuum.

2.5. A `new bond problem' and topology analysis in detail

A comparison of the distances obtained experimentally at 300 K (R_O1⋯H4 = 2.459 Å) and theoretically (R_O1⋯H4 = 2.57 Å) with van der Waals radii [R_O + R_H = 2.720 Å (Batsanov, 2001)] suggests that a new II′ bond forms in the high-pressure phase, in accordance with the geometrical criterion for the existence of hydrogen bonds (Hamilton & Ibers, 1968; Koch & Popelier, 1995) (Fig. 1). Since changes in a system of hydrogen bonds can lead to changes in the elastic and piezoelectric macroscopic properties, we studied this issue in detail. A search was made for critical points of non-covalent bonds at various values of external hydro­static pressure (Table S4). The interatomic bond paths and the corresponding critical points of bonds made it possible to establish a molecular graph (Fig. S1) that satisfies the Poincaré–Hopf relationship (Bader, 1990), Unexpectedly, we found that the bond path corresponding to bond II′ does not exist. The electron density distribution (Fig. 7) also did not reveal the formation of any immediate bonding contacts between the O1 and H4 atoms.

Figure 7 Electron density of glycine. (a) At critical points of hydrogen bonds for β and β′ phases as a function of pressure. The gap corresponds to a 0.76 GPa phase transition and separates the region where the low-pressure phase exists from the high-pressure phase existence region, in accordance with experimental data (Tumanov et al., 2008); (b) in β′-glycine in the plane of molecular layers. The black dotted line corresponds to bond II, the red dotted line corresponds to the hypothetical bond II′ (not confirmed). Isolines interval is 0.02 a.u. Symmetry codes of acceptors in hydrogen bonds: (I) −x, y + ½, −z; (II) x, y, z − 1; (III) −x, y + ½, −z; (IV) −x, y + ½, −z + 1.

This result illustrates that conclusions related to the formation of hydrogen bonds based on the geometrical criteria (distances between atoms) may need correction after a more sophisticated topological analysis of the electron density distribution. Thus, the topological analysis did not confirm the formation of bond II′ during the phase transition: in contrast to what the geometric crystallographic criterion yields (Hamilton & Ibers, 1968; Koch & Popelier, 1995). As such, the number of hydrogen bonds in the β-glycine/β′-glycine system does not change during the high-pressure phase transition.

2.6. Piezoelectric properties of β-glycine

The piezoelectric coefficients for β-glycine are closely related to elastic characteristics, especially compressibility. Piezoelectric tensor axes were chosen as described in Appendix B. In the Voigt notation (Nye, 1985) d₂₁ = d_yxx, d₂₂ = d_yyy, d₂₃ = d_yzz, d₂₅ = d_yxz, d₁₄ = d_xyz, d₁₆ = d_xxy, d₃₄ = d_zyz and d₃₆ = d_zxy. Only a few (d₂₁, d₂₂ and d₂₃) are responsible for hydro­static compression in β-glycine. The piezoelectric coefficients calculated using Berry phase approach (Resta, 1994) are given in Table 2.

Table 2
Calculated piezoelectric coefficients of β-glycine (pm V⁻¹) at different external pressure values

P (GPa)	0	0.2	0.4	0.7
d21	−2.1	−2.0	−1.9	−1.8
d22	4.9	4.9	4.8	4.6
d23	−1.2	−1.2	−1.14	−0.9
d14	20.5	18.0	15.7	13.4
d16	−15.1	−13.8	−12.5	−11.2
d25	0.6	0.3	−0.1	−0.3
d34	−29.6	−25.1	−21.1	−17.1
d36	28.5	25.2	22.1	19.4

The absolute values of the coefficients decrease with increasing external pressure. One can expect a decrease in the polarization in the y direction due to the maximum hydro­static compressibility H_max (see coefficient d₂₂). The piezoelectricity tensor was visualized using the MTEX program (Bachmann et al., 2010). This allowed us to determine the directions of maximum hydro­static polarization due to the compression of β-glycine via the piezoelectric coefficients (Fig. 8). The compressibility of hydrogen bonds III and IV can be related to quantum pressure in corresponding BCPs. An increase in external pressure leads to the decrease of maximum hydro­static compressibility (which means a decrease in the compressibility of weaker bonds III and IV). The latter is related to the E_H bond value when the corresponding bonds are compressed (quantum electronic pressure in BCPs decreases). We assume that quantum electronic pressure is directly proportional to the polarization along these hydrogen bonds.

Figure 8 (a) Visualization of hydro­static components of the piezoelectricity tensor in relation to the crystal structure of β-glycine. (b) Polarization vector P is oriented along the axis of maximum hydro­static compressibility, which is parallel to the y axis (and b lattice vector). Red, green and blue colors of dotted lines correspond to II, III and IV hydrogen bonds, respectively. Symmetrically non-equivalent glycine molecules are colored in green and red. Symmetry codes of acceptors in hydrogen bonds: (II) x, y, z − 1; (III) −x, y + ½, −z; (IV) −x, y + ½, −z + 1.

Only parameter d₁₆ = 178 pm V⁻¹ for β-glycine, measured by resonance piezometry, has been reported in the literature (Guerin et al., 2018). In the same work, a calculated value of this coefficient d₁₆ = 195 pm V⁻¹ was also reported. The possible reasons for such a discrepancy are discussed here in the supporting information.

Thus, increasing external pressure resulted in a decrease in the polarization in the direction of maximum hydro­static compressibility (H_max), which, in turn, is closely related to the microscopic characteristics [E_{H bond} and QEP(r)] in BCPs of the hydrogen-bond system of β-glycine. It is natural to associate the deterioration of the piezoelectric response when external pressure is increased with the upcoming change in the asymmetry of the weak bifurcated hydrogen bonds with `switching-over' of its branches III and IV pointing in the y direction.

The piezoelectric constants for β-glycine have not been measured experimentally yet, excluding coefficient d₁₆ (Guerin et al., 2018). The possible reason may be in the difficulty of growing single crystals that would be large enough for such a measurement. Thus, the experimental determination of these values remains a challenge for experimentalists. This is important because the piezoelectric properties of β-glycine are used in nanodevices, for example, for biomedical applications (Boldyreva, 2021). The relation between the networks of hydrogen bonds, the mechanical properties and electron polarization can stimulate further studies.

3. Conclusions

In this work, relations between the changes in microscopic quantum pressure of the electronic continuum and the macroscopic compressibility, piezoelectric properties, geometries and energies of hydrogen bonds in the crystal structure of β-glycine were established. Compression of β-glycine crystals under the external hydrostatic pressure is accompanied by a change in the internal quantum pressure in the intermolecular space and along the weak noncovalent hydrogen bonds III and IV, which are the two branches of a bifurcated hydrogen bond. The quantum pressure in the intermolecular space decreases, which corresponds to an increase in the electron concentration in both bonds. At the same time, on compression further to 1.7 GPa, the compression regions of these bonds are equalized (relative to each other). On the one hand, this leads to the strengthening of hydrogen bonds (decrease in the H⋯O distance and increase in the bond energy). On the other hand, this results in a change in macroscopic properties – hydro­static compressibility in directions corresponding to the directions of hydrogen bonds and bulk compressibility of the crystal as a whole.

The reversible phase transition in β-glycine into the high-pressure phase, β′-glycine, at 0.76 GPa is an example of a first-order phase transition that does not break the single crystal. This is of particular interest for molecular crystals and solids, which can potentially be used as functional materials due to their structure-correlated properties. It has been found that as the external pressure increases, weak hydrogen bonds are `switched over', which is accompanied by a rotation of every second zwitterion. Due to this rotation, the parameter c and the unit-cell volume are doubled. In this case, relatively strong hydrogen bonds forming molecular layers prevent a crystal structure from a significant rearrangement and the crystal from disintegration. Equalization of weak hydrogen bonds is one of the reasons for the high stability of β′-glycine on further hydro­static compression.

The predicted piezoelectric properties of β-glycine indicate that when hydro­static pressure is applied, the piezoelectric response decreases. This prediction is especially important because of the complexity of experimental measurements of piezoelectric coefficients. This indicates how the piezoelectric properties can be controlled via pressure variation. The polarization in the direction of maximum compressibility, y (which is parallel to b), of this polymorph decreases as well. Piezoelectric properties are related to the weaker hydrogen bonds (III and IV), and the change of the piezoelectric properties under external pressure correlates with the changes in the asymmetry of the bifurcated hydrogen bonds.

Using the concept of quantum electronic pressure makes it possible to comprehensively describe the pressure-induced phase transition. The piezoelectric and elastic properties of β-glycine are considered in relation to the changes in the 3D architecture of the electronic continuum. As long as microscopic and macroscopic properties are related, quantum electronic pressure serves as a universal electron density descriptor for solids/crystals under external compression.

4. Related literature

The following references are cited in the supporting information: Baima et al. (2016), Wu et al. (2005).

APPENDIX A

Quantum electronic pressure

The local stress tensor of the second rank, σ(r) (Rogers & Rappe, 2002; Tsirelson et al., 2016) in the framework of density functional theory (Dreizler & Gross, 1990) reads as (Tao et al., 2008):

Here σ_kin(r) is the quantum contribution of the electron kinetic energy, term σ_xc(r) takes into account the exchange-correlation electronic effects and σ_es(r) is the classical component, which includes the electrostatic interactions of electrons and nuclei. Treating the latter as a contribution of the external field, one can write the density of the quantum stress tensor as

Here σ_kin(r) and σ_x(r) are kinetic and exchange-correlation components, respectively. The Coulomb electron correlation is small compared with the exchange correlation, so it was neglected (Burke et al., 1998).

We define the kinetic part of the quantum stress tensor according to Tao et al. (2008) and adopt the minimal geometric coupling of electrons to the deformations, which originates from the description of quantum dynamics in local non-inertial frames (Tokatly, 2005). This definition allows for a physical interpretation of the local stress in terms of the internal energy response to a local deformation of small elements of the electron fluid. Explicitly, the Kohn–Sham kinetic stress tensor in the atomic units reads as

where ψ_lk(r) are Kohn–Sham spin orbitals and ρ(r) is electron density.

By definition, local pressure is p(r) = − . Therefore, the kinetic and exchange-correlation contributions to the quantum electron pressure at each space point is obtained by taking the trace of the corresponding parts of the stress tensor,

Here v_x(r) is electron exchange potential, which should be specified; s(r) = |∇ρ(r)|/2k_F(r)ρ(r), k_F(r) = 3π²ρ(r)^1/3. The electron exchange energy density, e_x(r), is defined according to (Dreizler & Gross, 1990).

In the framework of orbital-free DFT, the kinetic part of the internal quantum hydro­static pressure at each point of the electronic continuum has the form

where t_s(r) is positive-definite kinetic energy density of non-interacting electrons according to Kohn–Sham (1965), which may be represented as an explicit functional of the density. We adopted the simple approximation based on the second-order gradient expansion (Kirzhnits, 1957):

This function shows incorrect asymptotic behavior near the nuclear positions; therefore, small areas around nuclei should be excluded from consideration.

The exchange component of the pressure, p^x can be taken into account using the exchange-correlation functionals of DFT. Within the local electron density approximation (LDA), the exchange contribution to the electronic internal pressure reads as

The final form of the approximate quantum pressure function is

The accepted definition of pressure is the closest to the macroscopic definition of pressure in classical thermodynamics as the response of the internal Helmholtz energy to a change in volume for a fixed number of particles. In the same way, the quantum stress tensor (Rogers & Rappe, 2002) describes the change in internal energy (virtual work) due to the virtual local deformation of an infinitesimal element of the continuum, provided that the number of particles inside the deformable element at the point r is preserved, while its shape is unchanged.

APPENDIX B

Piezoelectric tensors

Piezoelectric tensor for β-glycine (space group P2₁, axis 2₁ is parallel to axis y) reads as

Chosen tensor axes for visualization in the MTEX program (Bachmann et al., 2010) are X||(bc), Y||(ac), Y||c. In the Voigt notation: d₂₁ = d_yxx, d₂₂ = d_yyy, d₂₃ = d_yzz, d₁₄ = d_xyz, d₁₆ = d_xxy, d₃₄ = d_zyz, d₃₆ = d_zxy. Coefficients d_yzy(d₂₅), d_zxy(d₃₆), d_zyz(d₃₄), d_xyz(d₁₄), d_xxy(d₁₆) are associated with shear stresses in the corresponding planes, whereas d_yxx(d₂₁), d_yyy(d₂₂), d_yzz(d₂₃) correspond to hydro­static stress.