An atomistic mechanism for elasto-plastic bending in molecular crystals

Mechanically flexible single crystals of molecular materials offer potential for a multitude of new directions in advanced materials design. Before the full potential of such materials can be exploited, insight into their mechanisms of action must be better understood. Such insight can be only obtained through synergistic use of advanced experimentation and simulation. We herein report the first detailed mechanistic study of elasto-plastic flexibility in a molecular solid. An atomistic origin for this mechanical behaviour is proposed through a combination of atomic force microscopy, μ-focus synchrotron X-ray diffraction, Raman spectroscopy, ab initio simulation, and computed elastic tensors. Our findings suggest that elastic and plastic bending are intimately linked and result from extensions of the same molecular deformations. The proposed mechanism bridges the gap between contested mechanisms, suggesting its applicability as a general mechanism for elastic and plastic bending in organic molecular crystals.

Electronic Supplementary Material (ESI) for Chemical Science. This journal is © The Royal Society of Chemistry 2023 S1 EXPERIMENTAL DETAILS S1.1 Materials| 2-amino-6-methyl pyridine (98%; CAS No:1824-81-3) and 3.5-dichlorosalisaldehyde (98%; CAS No: 90-60-8) were purchased from JK Chemicals and abcr respectively and used as received. S1.2 Crystal growth| The compound, 2,4-dichloro-6-[(6-methylpyridin-2-ylimino)methyl]phenol (DMP) was synthesized by liquid assisted mechanochemical grinding of equimolecular mixture of 2-amino-6methyl pyridine (108 mg, 1 mmol) and 3.5-dichlorosalisaldehyde (175 mg, 1 mmol) in a mortar and pestle using 100 L of methanol. Slow evaporation of a saturated CH 2 Cl 2 solution of ground powder at room temperature in beaker yielded long needle-shaped (size: 5 cm x 1 mm x 2 mm) orange coloured crystals. S1.3 Single Crystal X-ray Diffraction| Single crystals of 2,4-dichloro-6-[(6-methylpyridin-2ylimino)methyl]phenol were individually mounted on a glass tip. Single crystal XRD data were, collected on a Bruker D8 Venture system with graphite-monochromatized Mo Kα radiation (λ = 0.71073 Å) at 300 K and 100 K. Data reduction was performed with Bruker AXS SAINT [1] and SADABS [2] packages. The structure was solved by SHELXS 2018 [3] using direct method and followed by successive Fourier and difference Fourier synthesis. Full matrix least-squares refinements were performed on F 2 using SHELXL-2018 [3] with anisotropic displacement parameters for all non-hydrogen atoms. All other calculations were carried out using SHELXS 2018, [3] SHELXL 2018, [3] and WinGX (Ver-1.80), [4] whereas Mercury v3.6 [5] was used to visualize and to draw some of the figures for the structures. Detailed data collection strategy and structure refinement parameters along with crystallographic data are presented in Table S3.1. S1.4 -Focus Raman spectroscopy| All Raman spectra were collected using a Horiba Jobin Yvon Labram HR800 Raman microscopy system with an Olympus BX41 microscope. A continuous-wave diode-pumped solid-state laser ( nm) was used for Raman excitation. The laser spot was focused onto the sample = 633 by a 50×/N.A. = 0.75 objective (N.A. = numerical aperture). The scattered light was collected by the same lens, with the reflected and Rayleigh scattered light filtered by a bandpass filter. The Stokes-shifted Ramanscattered light was dispersed by an 1800 mm -1 grating. The signal was detected by a Peltier-cooled (-60 o C) charge-coupled device (CCD) camera. The spectrometer entrance slit was 100 μm wide with a confocal pinhole 1000 μm in diameter. With this configuration, the spectral resolution varies from 0.5 cm -1 per CCD pixel (at 100 cm -1 Raman shift) to 0.3 cm -1 per CCD pixel (at 3800 cm -1 Raman shift). The laser spot resolution is approx. 1 μm x 20 μm. 1 S1.5 Atomic Force Microscopy| Atomic force microscopy (AFM) measurements were performed with an MFP-3D Asylum (Oxford Instruments Asylum Research Inc., Santa Barbara, CA). AFM Force distance curves (FDC) were recorded with a frequency of 1 Hz in a force volume with an area of 10 x 10 µm and 10 x 10 measuring points yielding 100 single curves. All FDCs were corrected for the point of contact and subsequently averaged with SOFA. As AFM probe a biosphere B150-NCH (Nanotools GmbH, Munich, Germany) was used with a spring constant kc = 36.5 N/m, determined by a non-invasive thermal noise method. The tip is made from electron beam deposited high density carbon, with a tip radius R = 150nm ± 10nm, a Young's modulus E tip = 350 GPa and a Poisson ratio ν tip = 0.3.

S1.6 Microfocus Synchrotron X-ray Diffraction|
Microfocus synchrotron X-ray diffraction data were collected on the precision diffractometer equipped with a hybrid photon counting detector, EIGER X 1M detector (DECTRIS) in the SPring-8 BL40XU beamline. The X-ray beam (λ = 0.81042 Å) was focused to 0.922 (vertical) x 3.67 (horizontal) μm using a zone plate. To irradiate only a focused beam, a 30-μm diameter sorting aperture and a 40-μm-diameter centre stop were used. One end of the crystal was first glued on a pin and the crystal was bent to form a loop, attaching the second end to the pin by glue. This pin was mounted on the goniometer head for data collection. The beam was focussed on the middle of the bent position. After measurement of diffraction from the centre of the elastically bent crystal, one end of the crystal is released to confirm reversibility of the bending phenomenon. The same diffraction measurements were carried out on the re-straightened crystal, confirming that the crystals returned to their original configuration after removal of the deforming stress. The focused X-ray was incident normal to the bending plane (ω = 0°). The ω range, oscillation angle (Δω) and exposure time were ±10°, 0.1° and 0.5 s, respectively. Data were collected at room temperature. The measured position was moved from the convex to the concave side of the bend in steps of 50 μm, and 500 diffraction images were collected from each position.

S2 COMPUTATIONAL DETAILS
S2.1 Geometry optimisation and phonon simulations| Solid state simulations were performed within the framework of plane wave Density Functional Theory (DFT), as implemented in CASTEP v19.11. 2 Input structures were taken from experimentally determined geometries (See Section S3, S11, and S14). The electronic structure was expanded in plane waves to a kinetic energy cutoff of 1000 eV, and sampled on a Monkhorst-Pack 3 grid of spacing 0.05 Å -1 . Simulations were performed using the exchange-correlation functional of Perdew-Burke-Ernzerhof (PBE), 4 along with the D2 dispersion correction. 5 Norm-conserving pseudopotentials were generated 'on the fly' as implemented in CASTEP. Convergence on the geometry optimization was accepted with energy convergence < 3 x 10 -8 eV/atom, electronic eigenvalue convergence < 1 x 10 -12 eV, residual atomic forces < 1 x 10 -3 eV/Å, and maximum stress of 1 x 10 -2 GPa. The charge density was integrated over a fine grid of G max =64.8 A -1 , to facilitate convergence of vibrational frequencies.
The vibrational frequencies were calculated within the framework of density functional perturbation theory (DFPT) at the Brillouin zone centre (k= ). 6 LO-TO splitting was not explicitly considered.

S2.2 Elastic Tensors|
Elastic tensors were calculated at PBE-D3/def2-mSVP level of theory as implemented in CRYSTAL17. 7 The D3 correction 8 was applied with BJ damping. Following recent indications that computational elastic tensors can be qualitatively produced using semi-empirical methods, we additional calculated the elastic tensors using the HF3c 9 and sHF-3c methods, as implemented in the same software suite. In all cases, the electronic structure was sampled using SHRINK factors 4 4, with convergence of energy gradients of 4E-5 eV. The remaining parameters were held at the default values.

S2.4 Energy Framework Calculations|
The energy frameworks calculations 10 relating to intermolecular interactions of DMP were performed using the software suite Crystal-Explorer 11 based on B3LYP/DGDZVP molecular wavefunctions calculated using the CIF files. For calculations, the hydrogen atoms were normalized to standard neutron diffraction values. The energy frameworks constructed were based on the crystal symmetry and total interaction energy components of which included electrostatic, polarization, dispersion and exchange repulsion components scaled by 1.057, 0.740, 0.871 and 0.618, respectively. The interaction energies below 5 kJ.mol -1 are omitted for clarity and the cylinder thickness is proportional to the intermolecular interaction energies along the parallel vector passing through the cylinder.

S3.1 Face indexing of DMP single crystals|
To correlate the macroscopic bending against molecular packing, the orientation of DMP single crystals were indexed by means of single crystal X-ray diffraction,  The experimentally obtained unit cell of DMP (see Table S3.1) was fully relaxed at PBE-D2 level of theory (see Section S2). The relaxed unit cell dimensions are in good agreement with the experimentally obtained values, Table S3.2, with the unit cell volume underestimating the 100 K crystal structure by < 3 %. Variable temperature experimental data suggest significant thermal expansion of DMP single crystals, with the 300 K and 100 K structures differing by ca. 3%. Hence, this slight underestimation between the DFT (0 K) and experimental (100 K) geometries is expected and reasonable. The single crystal structures suggest that no major structural changes occur upon cooling to 100 K, Table  S3. 1 and Fig. S3.3, beyond the thermochromic transition (see Section S5). Residual electron density maps at 300 K and 100 K display an obvious change of electron density volume of the enol hydrogen (-O1-H1) with decreasing temperature (Fig. S3.4). Residual density suggests proton transfer towards the imine nitrogen (N1) was not observed. This is likely due to strong scattering of Cl atoms which precludes the detections of presumably partially transfer proton.   Table S3.3. The interaction energies due to π-π stacking interactions along the b-axis are the most dominating interactions (−52.5 kJ.mol -1 ) that stabilizes molecular columns running along the needle axis. In contrast, the interaction energies along the herringbone chain are −21.9 kJ.mol -1 , and the interchain interaction (along the a axis) energies are −82.3 kJ.mol -1 . Based on the previous model, calculation suggests plasticity should dominate this system. Despite the experimental finding shows marked elasticity before goes to the plastic deformation. This could be attributed due to the zig-zag interlock crystal packing.  deformation was fully reversible, but became irreversible at higher loadings (e.g. Fig. S4.1o). In contrast, three-point bending over the (100) face, Fig. S4.2, shows a small degree of macroscopic elasticity at low deforming loads, followed by partial fracture of the convex surface at higher loadings. Interestingly, the unfractured crystal remains flexible (see Fig. S4.2 e-f). Following from the apparent selfhealing of DMP crystals during AFM deformation over this face, we suspect the apparent elasticity stems from self-healing of microfracture on the (100) face (see discussion of the atomic force microscopy results in the main manuscript) along the crystallographic (001)/(001−) face: a-c) Elastic bending: the crystal regains its pristine shape after withdrawal the stress. i-k) Brittle fracture: crystal broke into pieces by increasing stress.

S4.2 Calculation of Elastic Strain|
Using microscopic images of bending, Fig. 4.3, the elastic strain of elastically bent single crystals of DMP were calculated according to the Euler-Bernoulli equation, which relates the beam thickness, t and radius of curvature, r=d/2, to elastic strain, Hence, the elastic strain observed for DMP is 2.45 %

S5 Thermochromic Behaviour of DMP
o-Hydroxy Schiff base molecules are well-known to exhibit thermochromism, resulting from intramolecular proton transfer. Small amounts of DMP crystals were put into glass test tubes which were then submerged into liquid nitrogen. DMP crystals changed colour from orange to yellow, and upon removal from liquid nitrogen the crystal return to their original orange colour, Fig. S5.1. To further study the reversible thermochromism of DMP crystal, we controlled the temperature of the crystal by using a Linkam LST 420 hot-stage. The videos were collected between 20 C and -170 C under a NikonAZ100 MULTIZOOM microscope with a DS-Ri1 12.7 megapixel camera. The crystal started to change colour from orange to light orange to yellow ( Fig. S5.2; Video S3). Upon re-heating the crystal to 20 C, the crystal returned to its original orange colour ( Fig. S5.2; Video S3). A plastically bent crystal also showed the same reversible thermochromism over similar controlled temperature conditions ( Fig. S5.3

; Video S4).
To verify if the temperature affects the elasticity, we first cooled a metal plate by adding excess liquid nitrogen. The crystal was placed on the cooled plate (became visibly yellow) and bent around a metal rod using wooden sticks. The elasticity of DMP crystals was also observed at low temperature with distinct colour change from orange to yellow. The crystal could be bent repeatedly under the low temperature condition (Video S3). Nikon AZ100 MULTIZOOM microscope with a DS-Ri1 12.7 megapixel camera. Nikon AZ100 MULTIZOOM microscope with a DS-Ri1 12.7 megapixel camera.  Importantly, crystal structure analysis revealed that no structural transition occurred in DMP across the thermochromic transition (see Section S3), suggesting the mechanical properties could be largely preserved at low temperature. Moreover, the thermochromic behaviour of these crystals appeared to be unaffected by bending, with plastically bent crystals transitioning continuously and reversibly from orange to yellow across the same temperature range as the straight crystals, Fig. S5.3. Our results therefore suggest that flexible o-hydroxy Schiff base crystals may be suitable for applications as flexible thermal-sensors or thermo-optical devices and will be explored in a detailed follow up investigation.

S6 Force-Displacement Curves -Atomic Force Microscopy
Force-displacement curves (FDCs) for a crystal of DMP were measured by AFM with different maximum loads as shown in Fig. S6.1 (red markers = 1.5 µN and brown markers = 3.5 µN). Up to an applied load of 1.5µN both experimental curves are in very good agreement with a fit by the Hertz theory (black dashed line) which describes the elastic deformation of a homogeneous material. From this fit the Young's modulus of the sample can be estimated as E =5.7 GPa (with Poisson ratio ν = 0.3). However, for higher applied forces > 1.5µN (brown markers) a characteristic deviation from the fit is observed, which is due to yield, i.e. plastic deformations. This can also be seen in the hysteresis of the approach (full markers) and retraction (empty markers) part of the FDC, as well as in the offset of the displacement D at F = 0. For applied forces < 1.5 µN the hysteresis is small ( Nm) and the deformation recovers nearly completely (dD < 1 nm), which is consistent with elastic deformations. However, for applied forces > 1.5 µN the hysteresis is significantly higher (( Nm) and deformations do not recover ℎ = 2.5 × 10 -14 (dD = 9nm).

S7 Synchrotron Based Microfocus X-Ray Diffraction
Typical laboratory X-ray diffraction uses a beam size on the order of 100-1000 m in diameter. Correspondingly, resolving unit cell geometry across the bent positions of a single crystal is impeded by averaging over a large spatial distribution. To overcome this challenge we have used a micro-focus beam (0.922 × 3.67 μm beam size) available at the SPring-8 synchrotron facility (BL40XU beamline) to capture the spatially resolved evolution of unit cell geometry across each of the elastically and plastically flexible crystals, Fig.s 7.1 and 7.2, respectively.   imaginary (negative) frequencies confirms that our simulated structure corresponds to an energetic minimum on the potential energy surface.

S8.2 Experimental Raman Spectra|
The effects of mechanical bending on DMP single crystals were analysed by -focus Raman spectroscopy to explore the potential influence of bending on molecular deformation, Fig. 8.1. For the same elastically bent single crystal, Fig. 8.1b and 8.1d, the Raman spectra at the inner (IA), middle (MA), and outer arc (OA) remain effectively identical, with no detectable shift in the positions of the vibrational bands, Table S8.2. This strongly suggests that no significant perturbations on the molecular structure occur during elastic bending. We do note minor changes in the relative intensities of the Raman bands as compared with the straight portion (SP) of the elastically bend crystal and a separate single crystal (S), which presumably result from a change in the orientation of the crystal in the polarization direction of the laser, or due to formation of defects during bending. The same effect is observed in plastically bent single crystals of DMP, Fig. 8.1c, Fig. 8.1e, and Table S8.3, again highlighting that no major molecular level deformations occur during bending. We note that comparison of intensities between different crystals is precluded by different concentration of defects and microstructure, leading to difficulties in extracting meaning when comparing the intensities of data collected from straight crystals, elastically bent crystals, and plastically bent crystals. The only notable difference in the Raman spectra of the plastically bent crystals is the apparent splitting of the Raman band at ca. 190 cm -1 (M4 in Table S8.2). This effect is presumably due to the anisotropic strain causing resolution of closely spaced vibrational bands (see Table S8.1) in DMP.

S9 Anisotropic Compression of DMP
To assess the effects of anisotropic mechanical deformation on the structure of DMP we imposed anisotropic compression/tension over the unit cell. To mimic the direction of compression / tension experienced during elastic bending, expansion and contraction of the unit cell was conducted along the baxis (i.e. the long axis of the crystal). The b-axis was held constant at the given strain value ( = ) and the rest of the unit cell was allowed to relax in response to the perturbation, Table   100 × ( -0 )/ 0 S9.1. We note that the expansion / contraction of the b-axis was exaggerated beyond the stresses typically observable in molecular crystals. This was done to emphasis the structural response to bending. It is worth highlighting that despite applied strains ranging nearly 10%, the crystal volume varies by only ca. 3.5%. This demonstrates the relaxation of the crystal unit cell in directions orthogonal to the applied strain as being critical to the material response. This is consistent with experimental observation from -focus XRD (e.g. Fig. 3 in the main manuscript). Accompanying the adjustment of unit cell dimensions, anisotropic compression leads to an increase on the internal energy of the unit cell, Fig. S9.1. This strain energy represents the excess energy held by the stressed crystal, i.e. the energy the crystal will aim to alleviate through deformation such as delamination (see Section S). From -focus XRD, a strain of ca. 2% is expected in elastically bent crystals of DMP ( Fig.  3 in main text), which is consistent with the calculated elastic strain calculation suggestion of % (ESI ≈ 2.5 S4). According to our DFT models (PBE-D2), this degree of strain is accompanied by a ca. 3-kJ.mol -1 increase in the internal energy of system. These energies are well within reach of standard experimental procedures and thermal energy at room temperature (ca. 2.5 kJ.mol -1 ), and are on par with the energies that separate normal solid-state processes and weak intermolecular interactions. It is therefore very reasonable to assume that these energies are achieved upon bending without causing significant damage to the system.  When crystals of CMP were bent over the (010) face, exceptional and distinct plasticity was observed, Fig.  S10.1. Plasticity over the (010) face corresponds to slip in parallel planes that run along the needle axis (the ac plane). Consistent with DMP, this corresponds to delamination of the planes perpendicular to the herringboned layers. However, owing to the lack of strong interactions across these plans in CMP, plasticity is observed, rather than the elasto-plastic behaviour exhibited by DMP.

S10.2 Energy Framework Calculations of CMP|
The intermolecular interactions in CMP were calculated using the energy framework calculation (Fig. S10.2 and Tables S10.2). These calculations showed that the π⋯π stacking along the a-axis are the strongest interactions in CMP structure (-50.7 kJ.mol −1 ), stabilizing the molecular columns that run along the needle axis. We note that the interactions in CMP are ca. 2 ⋅⋅⋅ kJ.mol -1 weaker than in the DMP structure (see ESI S3.3), suggesting that compression / expansion of the a-axis in CMP should occur more readily. In contrast, the total interaction energies along the herringbone chain (b axis) are -19.6 kJ.mol -1 , and the interchain interaction energies (along the c axis) are -73.5 kJ.mol -1 . Again we find that the interchain interaction energies in CMP are significantly lower than those in DMP (73.5 vs 82.3 kJ.mol -1 ), indicating that delamination can occur much more readily in CMP than in DMP, presumably responsible for the markedly different mechanical behaviour of the materials.
Overall, the structure of CMP has anisotropic non-covalent interactions, consistent with qualitative models of plastic deformation 13

S11 Anisotropic Compression of CMP
To better capture the effects of bending on the structure of CMP we conducted a series of ab initio simulations at the PBE-D2 level of theory, shown to reproduce the experimental crystal structure well, Table S11.1. Computational details are given in ESI S2. Our simulation methodology leads to an overall good agreement with the experimental unit cell geometry, slightly underestimating the experimental geometry, presumably due in large to the thermal expansion of the material. To observe the effects of uniaxial compression on CMP we systematically expanded / contracted the unit cell along the crystallographic a axis, and allowed the remaining crystallographic parameters to relax, Table S11.2. The a axis was selected as it corresponds to the strained axis during bending (see Fig. 4 in the main text). Upon compression, the herringboned layers clearly flatten, leading to elongation of the crystallographic b axis, and shortening of the crystallographic c axis, Fig. S11.1. In contrast, during expansion, the opposite trend is observed. This is consistent with the structural effects observed in DMP (see Fig. 4 in the main text).  The internal energy of CMP increases slowly with anisotropic compression, Fig. S11.4, at a slower rate than for DMP, Fig. S9.1 and Fig. 5 in the main text.

S12 SLIP PLANE ENERGIES IN DMP
The general model for bending suggests that the slip plane energy changes as a function of anisotropic compression along the bending axis. To explore this, we imposed systematic distortions of molecules along the respective slip planes by a factor , where describes = / the Cartesian displacement in unit cell direction a (or b, c). A value of thus corresponds to a = 1 full translational unit of the unit cell dimension. In the perturbed geometry we performed single point energy calculations to obtain an estimated slip plane energy. At selected perturbations, the ionic positions and unit cell were left to relax. Relaxation was performed whilst freezing uniaxial compression and ionic positions in the slip direction. In this way, we ensure that structural relaxation occurs only in the dimensions that are not directly affected by the structural deformation of the model.
When considering the slip plane energies of DMP as a function of uniaxial compression/expansion ( ), Table S12.1, it is clear that small amounts of distortion have marked impact on the subsequent ability of DMP to undergo plastic deformation. This is reflected in the Poisson's ratio of DMP (see Section S13), which may act as an indicator for the relative elasto-plastic threshold for flexible crystals. However, a larger family of crystals with this behaviour are first required. Structural analysis of the simulated structures suggests that the reduction in slip plane energy with anisotropic compression results from the elongation of the herringbone chains (c axis in DMP), as well as of the slip planes, Fig. S12.1.

S13.1 Elastic Tensors of DMP|
Prior to simulating the elastic tensors, the geometry was relaxed at the corresponding level of theory, with the PBE-D3(BJ) functional instead leads to a ca. 5% underestimation of the 100 K unit cell geometry, dominated by an underestimation of the unit cell b-axis (along the herringboned chains). Following from the relative unit cell volumes in Table 13.1, it is consistent to find that the elastic tensor ( ) and its inverse the stiffness tensor ( ) are softer for HF-3c as compared with sHF-c3, Eqns S13.  14 provide reasonable estimate of the elastic tensors for molecular materials. Consistent with the underestimated unit cell volume, PBE-D3(BJ) predicts a somewhat more rigid material as compared with both HF-3c approaches. We note our AFM measurements (Fig. 2 in the main text) indicate a room temperature Young's modulus of ca. 6 GPa. This is generally consistent with our simulations when the effects of thermal expansion are considered (see for example Table S13.1 the softening of elastic tensors when a change in volume of ca 4% is considered).

S13.2 Elastic Tensors of CMP|
The geometry of CMP was relaxed at the corresponding level of theory, Table S13.3. The HF-3c method 9 underestimates the low temperature unit cell of CMP by 6.06%, providing marginally poorer agreement than for DMP. We note however, that the experimental geometry for CMP was measured at a higher temperature than for DMP, presumably accounting in part for this difference. By applying the scaling =0.7, the sHF-3c method instead 8 leads to a better agreement of the unit cell, underestimating it by 3.08%. The cell obtained using PBE-D3(BJ)/def2-mSVP is similar to that obtained by HF-3c, underestimating the 150 K unit cell by ca. 6%.

General Discussion of Mechanical Properties|
The elastic tensor describes the deformability of the crystal structure and therefore likely holds significant information regarding the mechanical flexibility of single crystal materials. In our proposed model we consider the effects of uniaxial compression along the bending axis (b for DMP and a for CMP). Analysis of the Young's moduli, Fig. S13.1, along the respective directions shows that both materials respond similarly to uniaxial compression. It follows that the amount of stress accumulated by each material in response to bending-induced compression should be similar and thus does not account for the different mechanical behaviour of the material. The Poisson's ratio for each material do however suggest that the CMP unit cell should exhibit a more significant structural response from uniaxial compression. This is consistent with the finding that CMP slip planes are more affected by uniaxial compression than are those found in DMP. The most significant difference between the two flexible crystals is the magnitude of their shear moduli, Fig. S13.2. For the elasto-plastic crystal DMP, the shear modulus along the slip plane vector is ca. 2.7 GPa. In contrast, the modulus along the slip plane of CMP is only ca. 1 GPa. This strongly supports the relative ease of delamination in the latter, and its lack of elasto-plastic response.