Atomistic modeling of martensitic phase transition in hexamethylbenzene

Abstract

Materials exhibiting martensitic phase transitions are essential for applications in shape memory alloys, actuators, and sensors. Hexamethylbenzene (HMB) has long been considered a classical example of ferroelastic organic crystals since Mnyukh's pioneering work in the 1970s. However, the atomistic mechanism underlying this phase transition has not been clarified. In this work, we present a molecular dynamics (MD) simulation to directly model the phase transition mechanism in HMB. For the first time, we report simulation results that accurately reproduce both the transition temperature and hysteresis loop observed in previous experimental studies. By analyzing the MD trajectories and potential energy surface, we identified that a low-barrier atomic sliding mode along the close-packed (111̄) plane of the low-temperature phase is key to triggering the phase transition at the critical temperature window. This is further confirmed by the observed continuous softening of shear modulus around the transition window. Our results demonstrate that the integration of various atomistic modeling techniques can provide invaluable insights into martensitic phase transition mechanisms in organic crystals and guide the development of new organic martensites.

---

I. Introduction

Phase transitions are critical in understanding the relationship between material structure and properties. Among various types of transitions, the diffusionless phase transition plays a pivotal role in numerous industrial applications, particularly shape memory alloys (SMAs).^1,2 When heat is applied to SMAs, they undergo a phase transition from low-temperature martensite to high-temperature austenite, a process known as the martensitic phase transition (MPT).³ Many applications of SMAs (e.g. Ni–Ti alloy) are often incorporated into thin films, where they function as micro-actuators, making them highly valuable in various technological applications.^4–6

Despite extensive studies of polymorphic phase transition in metals, alloys or polymers,^7–9 MPT in molecular crystals have remained relatively underexplored. Molecular crystals have received growing attention in materials sciences due to their diverse group of properties ranging from thermal and mechanical to optoelectrical functionalities.^10–12 Unlike the inorganic counterparts, there exist more design choices to develop new molecular crystals with target properties via modern crystal engineering techniques.^10–12

Recently, reconfigurable molecular crystals have spurred the development of a burgeoning research field known as crystal adaptronics.¹³ Several organic crystals were experimentally reported to exhibit martensitic phase transitions.^14–26 The martensitic phase transitions in organic crystals have been found to couple with a range of functional properties, including thermoelasticity, thermosalience, superelasticity, and shape memory effects,²³ which make molecular crystals highly attractive for modern smart material applications, as they offer unique responsiveness to external stimuli and are well-suited for a variety of technological uses in areas like adaptive materials, sensors, and actuators.

Among molecular materials, hexamethylbenzene (HMB) is one of the very early examples that has been reported to exhibit martensitic phase transition in the history of small molecule organic crystallography. At low temperatures, HMB molecules crystallize in a triclinic P1̄ symmetry (referred to as HMB-II).²⁷ In 1952, Saito studied the phase transition of HMB, and found that HMB-II undergoes a phase transition at high temperatures, forming what is now commonly referred to as the orthorhombic Fmmm form I in the literature.²⁸ Although the HMB II-I transition had been observed for some time, it was not initially recognized as an example of MPT at that time. In the 1970s, Mnyukh performed a series of pioneering studies on the polymophic phase transitions of molecular crystals.^29–33 Among these studies, Mnyukh proposed that HMB, together with several other systems, possesses a likely martensitic phase transition behavior.³¹ For HMB, he determined that form II transforms to I at 382–385 K upon heating and reverse transition occurs at 380–382 K upon cooling. Based on X-ray data analysis, Mnyukh also suggested a molecular sliding mechanism to explain the crystallographic relation between two forms. Recently, Li and coworkers revisited the study of HMB's phase transition and further confirmed its martensitic nature.³⁴ More interestingly, they found that the rapid and complete structural switching behavior of HMB enables it to perform work comparable to many actuator materials. This exciting discovery suggests that organic crystals like HMB hold potential for applications in soft devices, such as actuators.

To our knowledge, identifying a new organic martensite in the laboratory often requires extensive experimental efforts. Computational techniques, however, may offer potential to accelerate this process by enabling high throughput in silico screening. To facilitate computational discovery, it is essential to first verify whether current modeling techniques can reliably reproduce experimental phenomena. If successful, one may take steps further to develop a systematic computational pipeline to screen potential organic martensites from existing databases.

In this work, we selected HMB as a model system to conduct molecular dynamics (MD) based atomistic simulation to elucidate the phase transition mechanism of HMB. In the following sections, we begin by providing the necessary computational details of our modeling approach. We then demonstrate how MD simulations can reproduce the experimentally observed MPT through simulated heating and cooling cycles of the HMB crystal. This is followed by an energy surface analysis to quantify low-energy barrier pathways during the phase transition. Additionally, we explore the use of elastic constants and stress–strain relations as indicators to evaluate the possibility of observing an MPT by starting from the known room-temperature phase. Finally, we conclude by discussing how our findings can facilitate the automated screening of new organic martensites.

II. Computational methods

A. HMB crystal and model setup

As shown in Fig. 1a, the low temperature HMB form-II features a single molecule in the unit cell with the P1̄ symmetry. Due to the molecular symmetry, the asymmetric unit contains a half molecule. In the standard crystallographic setting, it has a triclinic box with a = 5.260 Å, b = 6.199 Å, c = 8.004 Å, α = 103.8°, β = 98.7°, and γ = 100.2°. Fig. 1b displays a 2 × 2 × 4 supercell model based on the standard setting. Recently, we have implemented a molecular layer finding function into the open source Python library PyXtal³⁵ that can search for the molecular layer spacing in any arbitrary crystallographic plane. Using this module, we find that the HMB-II crystal has the largest separation of 3.595 Å at the 111̄ plane (highlighted in blue) of the unit cell. This plane is also parallel to the benzene ring of the HMB molecule. As discussed in the previous section, it has been found that the shearing and molecular sliding activities on the 111̄ plane are the key to triggering phase transition. To simplify the analysis, we reoriented the unit cell model using a rotation matrix,^1,2 [1, −1, 0], [4, 3, −1]. As depicted in Fig. 1c, this transformation maps the 111̄ plane of the original unit cell onto the (001) plane of the rotated supercell. In addition, the resulting triclinic box has α = 96.0°, β = 91.0°, and γ = 90.0°, which are closer to an orthorhombic box, thus allowing more robust neighbor finding in the MD simulation. Based on the rotated crystal, we created a 6 × 5 × 3 supercell model consisting of 345 600 atoms in a triclinic box (see Table 1) for the subsequent MD simulations.

Fig. 1 The crystal structure of HMB. (a) The HMB unit cell with only one molecule; (b) the HMB supercell highlighting the closest packing (111̄) plane (colored in blue) of the unit cell; (c) the rotated HMB supercell by placing the closest packing plane to the (001) plane (colored in blue) in the new setting.

Table 1 The cell parameters for the P1̄ HMB-II phase at 293 K

Source	a (Å)	b (Å)	c (Å)	α (°)	β (°)	γ (°)
Unit cell-expt.^31	5.260	6.199	8.004	103.8	98.7	100.2
Unit cell-MD	5.372	6.442	8.257	105.1	100.1	98.1
Supercell-expt.	185.083	88.112	173.514	96.0	91.0	90.0
Supercell-MD	188.640	89.492	183.119	92.2	96.4	90.2

---

B. Molecular dynamics simulation

It is well known that the accuracy of MD simulation results heavily depends on the quality of the interatomic force field. In this work, we employed our previously developed pipeline³⁶ to automate the generation of general amber force field (GAFF) interatomic potential^37,38 using AMBERTOOLS.³⁷ Specifically, we used the semiempirical AM1-BCC method³⁹ to compute the atomic partial charges for each atom in the molecule. To validate the force field, we performed NPT calculations of HMB form II at 300 K and obtained cell parameters that closely matched the experimental values (see Table 1) based on the standard settings. For the rotated supercell, the a, b and γ values match very well. However, the c value is about 9.4% mismatch. This is not unexpected since the layer spacing is more subtle due to the weak interlayer interactions. Due to the error in c, the variations on α and β values are also notable. Despite the discrepancy in c, the overall agreement ensures the reliability of our subsequent simulation results. We used the LAMMPS⁴⁰ package to carry out the MD simulations, mimicking the experimental heating and cooling processes of the HMB crystal over a range of temperatures. In these simulations, we employed the Langevin thermostat⁴¹ to control the temperature and the Nose–Hoover barostat⁴² to maintain atmospheric pressure, with a 1 fs timestep. After testing several heating and cooling rates ranging from 3.33 × 10⁹ to 1.25 × 10¹⁰ K s⁻¹, we selected 4 × 10⁹ K s⁻¹, which provided the optimal balance between accuracy and computational efficiency.

C. Evaluation of properties

To evaluate the mechanical and thermodynamic properties associated with MPT in our system, we conducted a comprehensive set of calculations, including γ-surface energy, elastic constants, and stress–strain relationships under mechanical loading. These analyses help us understand the structural and energetic factors influencing phase stability, deformation mechanisms, and the mechanical properties of HMB in this study.

In the context of stacking faults and plastic deformation in metals, the concept of the γ-surface energy has been widely used to evaluate the spatial distribution of penalty energy required for a crystal to glide along a particular plane. To compute the γ-surface energy, we employed a systematic approach where a chosen slip plane was displaced incrementally while relaxing the surrounding lattice. For each displacement, the total energy of the system was calculated using LAMMPS, capturing the energy landscape associated with partial or complete shifts of the crystal along the slip plane. This method allows us to map the energy variations across different slip vectors, thus building a comprehensive γ-surface profile. As we will discuss in section III C, the resulting energy-surface profile is instrumental in predicting the ease of slip in the material, which is crucial to the martensitic transformation process.

For the computation of elastic constants, we start by applying Hooke's law:

σ_ij = C_ijklε_kl, (1)

where σ is the stress tensor, ε is the strain tensor, C is the second order elastic constant, and i, j, k, and l are direction indices. Here, we employed LAMMPS to calculate the elastic constant of HMB. Specifically, we started by preparing a set of equilibrated structures at the desired temperatures by running the NPT MD simulation. For each equilibrated configuration, we applied a range of uniaxial and shear deformations in both positive and negative directions to measure the averaged stress values after an NPT equilibration with the constraints of strain deformation. The elastic constants were then derived from these stress values according to eqn (1).

With the elastic constants determined, we calculated mechanical properties such as Young's modulus and shear modulus and analyzed their spatial distribution in both 2D and 3D representations. For this analysis, we used the open-source Python library ELATE⁴³ to explore the anisotropic behavior of the shear modulus as a function of temperature. Since a shear mode depends on two orthogonal directions (the direction of the applied shear stress and the direction in which the shear strain is measured), we followed the notation system (θ, ϕ, χ) as defined in ELATE to denote these shearing directions. For clarity, this notation is illustrated schematically in Fig. 2, where vector a is perpendicular to the shear plane that measures the strain, while b indicates the direction in which the stress is applied.

Fig. 2 The angular notation systems θ, φ, and χ used to define the spatial shear modulus, where a is perpendicular to the shear plane (highlighted in gray eclipse), while b indicates the direction in which the strain is measured.

To calculate the stress–strain relationship under mechanical loading, we adopted a similar approach to that used for elastic constants, but with a controlled, incremental application of strain. The deformation strain was applied at a specified rate (2 × 10⁻⁷ fs⁻¹), and during the simulation, stress and strain values were recorded for further analysis. The simulation was continued until a maximum strain of 20% was reached, allowing us to capture the stress response and analyze the material's deformation behavior under various loading conditions.

III. Results and discussion

In the section, we discuss how to computationally reproduce the experimentally observed phase transition, and extract the atomistic insight by analyzing direct MD simulations. Furthermore, we conducted a series of energy and property calculations to understand the origin of the martensitic phase transition from the perspective of atomistic interactions and molecular motions. Using HMB as a model system, we aim to explore the potential of integrating computational techniques to provide quantitative predictions of MPT in organic martensites.

A. Heat cycling simulation

We first attempted to set up the direct MD simulation to mimic the experimental heat/cool cycling as reported in a recent study.³⁴ Specifically, we performed (i) an initial equilibration of 1 ns at 350 K under NPT ensemble; (ii) step-wise heating from 350 to 390 K over a period of 8 ns; (iii) a subsequent equilibration of 1 ns at 390 K; and (iv) reverse cooling from 390 K back to 350 K at the same rate used during the heating phase.

Despite the significantly faster heating/cooling rate in our MD simulation (4 × 10⁹ K s⁻¹) compared to the practical experiments, our simulation successfully captures the phase transition behavior analogous to experimental observations. In Fig. 3, we track the evolution of two representative variables, average interlayer spacing (d) and potential energy per molecule (PE), as functions of temperature throughout the entire cycle. Both d and PE, taken as order parameters, clearly exhibit hysteresis loop behavior. During heating, we first observed a steady increase of d up to 373 K, suggesting a normal thermal expansion behavior. Between 373 and 375 K, a sudden increase in d within a short temperature window indicates the onset of a phase transition. Beyond 376 K, there is another steady expansion of d, indicating a thermal expansion of the new phase. A similar trend is observed during the cooling phase, except that the phase transition temperature window shifts to 366–368 K, resulting in a hysteresis loop in the overall curve. A similar hysteresis behavior is also observed in the evolution of PE values, as shown in Fig. 3b, further confirming the phase transition dynamics during the heating and cooling cycles.

Fig. 3 The evolution of two representative system variables upon the heating and cooling cycle simulation of HMB between 350 and 390 K. (a) The change of the averaged interlayer separation d values along [001] against the temperature. (b) The averaged potential energy as a function of temperature.

In previous experiments, Mnyukh³¹ and Li³⁴ reported the phase transition that happened around 382–385 K upon heating and 380–382 K upon cooling. Our simulation results are qualitatively consistent with the reported experimental values. We also attempted to model this process with different parameter settings. Indeed, we found that a slower cooling rate and a larger model size result in slightly higher phase transition temperatures. However, the overall physical behavior remains unchanged. Therefore, we conclude that our MD simulation can qualitatively reproduce the experimentally observed MPT phenomena, despite the limitations of force field accuracy and the faster cooling rates used in the simulations.

B. MD trajectory analysis

After confirming the reliability of our modeling approach, we proceed to study the atomistic mechanism of this phase transition by analyzing the molecular dynamics trajectories. Fig. 4a and b plot two MD snapshots at 350 K and 390 K to represent the low temperature phase form II and high temperature phase form I, respectively. Clearly, we observe a rapid shearing on the xz component during the transition window, corresponding to molecular sliding along the (001)[100] direction. Such molecular sliding is expected to result in a change of stacking sequence between the adjacent (001) molecular layers.

Fig. 4 Comparison of molecular arrangements in HMB before and after heat-induced phase transition. The upper panel depicts the low-temperature phase (HMB-II): (a) the side view in the xz plane, (c) a zoomed top view of the rectangle region in (a) at the xy plane, focusing on bilayer molecular centers colored in blue and red, and (e) the density distribution of each molecule relative to the nearest four molecular centers in the adjacent layer. The lower panel provides a similar analysis for the high-temperature phase (HMB-I) obtained from our MD simulation: (b) the side view, (d) the zoomed top view, and (f) the density distribution of molecular centers.

To better understand the sliding process, we select two adjacent (001) molecular layers and track the evolution of their molecular centers on the projected xy plane as shown in Fig. 4c and d. Notably, we found a common pattern in both phases. That is, each molecule (i.e., an individual red or blue sphere) is surrounded by four molecules from the neighboring layer that forms a rhombus like shape. However, subtle differences in molecular positioning were observed between the phases. In the high-temperature phase, the shape is closer to an ideal rhombus, with the central molecule nearly located at the geometric center. In contrast, in the low-temperature phase, the shape deviates more from a perfect rhombus, and the central molecule is slightly displaced from the geometric center.

We further calculated the distribution of red molecules relative to the blue rhombus by considering all MD trajectories during the phase transition. This results in density plots as shown in Fig. 4e and f. From the density plots, it can be seen that each HMB molecule in the high temperature phase is located at the center of the rhombus of the ajacent layer, while the HMB molecule in the low temperature phase is slightly off from the rhombus' center. These results suggest that there is a systematic shift between the adjacent molecular layers upon heating. This process also results in a crystal symmetry change from the triclinic P1̄ at low temperature to the orthorhombic Fmmm symmetry at high temperature.^28,31

Overall, the underlying molecular mechanism presented here closely aligns with Mnyukh's interpretation, based on his experimental data.³¹ Additionally, we note that similar mechanisms have been observed in several other systems, as discussed in a recent review.⁴⁴

C. γ-Surface energy calculation

Although molecular sliding has been recognized as the key for the martensitic phase transition, the relationship between the sliding motion, crystal packing, and intermolecular interactions remains unclear. Following this approach, we performed the γ surface calculation on HMB at 300 K. Specifically, we evenly split the supercell model into two regions along the z-axis, and then systematically displaced the upper region molecule in the entire xy plane, while keeping the lower region fixed. For each displacement, the system was relaxed to the energy minimum under the rigid body assumption. Fig. 5 shows the resulting energy surface contour plot that represents the difficulty of molecular sliding in terms of energy.

Fig. 5 The calculated γ surface energy contour plot for the rotated HMB supercell at the (001)-xy plane.

Starting from the origin O at (0, 0), we can see that there mainly exist two adjacent low energy basins (A and B). In particular, basin A has a slightly lower energy minimum value than B. In addition, the cost energy barrier from O to A is only 0.082 meV per atom. Upon heating, the barrier is expected to become even lower. Such a low energy barrier value is comparable to magnitude of thermal fluctuations (23.6 meV at 300 K), and thus the sliding from O to A can be easily triggered upon heating.

We also found that the O → A direction (approximately the [1̄10] direction of the rotated cell) agrees with the molecular displacement in Fig. 4e and f. The excellent agreements suggest that we may use the γ-surface energy analysis to predict the likely slip mode during the phase transition of molecular crystals without running long-time MD simulation.

D. Evolution of elastic properties

In addition to molecular sliding, the phase transition between HMB II and I can be alternatively understood from the mechanical shearing on the xy plane. Hence, we investigated HMB's elastic properties as a function of temperature, following the method described in section II C.

We began with analyzing the three-dimensional plot of the minimum shear modulus for HMB-II at 300 K. As shown in Fig. 6, the shearing plane (a direction) can be visualized imagining a vector connecting from the origin to any point on the plane. The distance from the origin to any point on a line represents the value of minimum shear modulus measurement among all b directions that is orthogonal to a. Clearly, Fig. 6 suggests that the weakest shear modulus on the xy plane is generally no more than 0.5 GPa and systematically smaller than other directions. In particular, the minimum shear modulus value of 0.37 GPa can be found at θ = 184°, ϕ = 128°, and χ = 103°, which is approximately consistent with the (001) xy plane and [1̄10] direction as found in our previous γ-surface energy analysis.

Fig. 6 The calculated spatial distribution of minimum shear modulus values in HMB II at 300 K.

Furthermore, we plotted the temperature-dependent evolution of the shear modulus at the yz and xz planes in Fig. 7. In both plots, each contour line represents the shear modulus at a specific temperature. From Fig. 7, it is evident that the shearing on xy is much weaker than that on the yz and xz planes. In addition, we can see that the shear modulus on all directions continuously decreases as the temperature increases, which is expected given that the thermal expansion should generally weaken the shear modulus.

Fig. 7 The evolution of HMB's shear modulus as a function of temperature at the (a) yz and (b) xz planes.

However, the continuous decreasing trend in shear modulus diminishes when we examine the evolution of shear modulus contour lines in the xy plane. To illustrate this trend, we plot the contour line of the G_xy with an interval of 10 K from 320 to 400 K, as shown in Fig. 8. Between 320 and 380 K, there is clearly a consistent decrease in the shear modulus of HMB-II. To aid the analysis, a red line is drawn along the direction of minimum shear modulus on this plane, corresponding to the path from basin O to basin A as identified in our previous γ-surface calculation. Along this direction, the shear modulus at 380 K drops to approximately 0.1 GPa, indicating substantial softening that suggests potential mechanical instability under shear. Indeed, further heating to 390 K produces a new contour shape, implying the onset of a new phase (HMB-I) due to shear instability. In HMB-I, the rearranged molecular stacking enhances mechanical stability by increasing the weakest shear modulus. These findings align well with our γ-surface analysis from the previous section, reinforcing the relationship between thermal effects and mechanical stability.

Fig. 8 The evolution of HMB's shear modulus as a function of temperature at the xy plane. The arrow denotes the direction leading to shear instability.

E. The impact of mechanical loads

In addition to heat-induced phase transition, there has been a wealth of studies on molecular crystal phase transitions that can be triggered by mechanical work.^44–46 Hence, we investigated the impacts of mechanical loading on HMB. In particular, we focus on the shear stress for a few representative temperatures and directions.

Fig. 9a plots the simulated shear stress–strain relation at three main directions ε_xy, ε_yz, and ε_xz, as well as the direction around [1̄10](001) (will be called ε_min from now on)⁴⁷ for HMB I at 300 K. The results reveal that the peak shear stress in the σ_xy direction (approximately 2500 MPa) is significantly higher than in the σ_xz and σ_yz directions (both around 400 MPa). Meanwhile, σ_min exhibits a fluctuating stress response of approximately 50 MPa. These findings are consistent with our expectations, as the weakest interactions in the crystal are predominantly found between the xy planes.

Fig. 9 The simulation of HMB-II under various shear strains. (a) The stress–strain relation with the applied strain on ε_xy, ε_xz, and ε_yz and the minimal stress ε_min directions at 300 K; (b) the corresponding energy change for ε_xz, ε_yz and ε_min. For clarity, all curves were smoothed by applying the Savitzky–Golay filter as implemented in SciPy with 20 coefficients.

When the peak σ_xz is reached at ε_xz = 0.15, it corresponds to an approximate displacement of 0.54 Å along the x-axis (assuming molecular spacing d = 3.595 Å). This displacement drives the structure away from basin O towards to B as shown in our γ surface energy plot (see Fig. 4), thus leading to a phase transition and a subsequent decrease in the stress value σ_xz, as shown in Fig. 9a. Under ε_yz, the structure undergoes a phase transition from basin O to A. However, ε_yz does not align well with the OA direction. Hence, it reaches a peak stress with a less strain value (∼0.1). In the case of ε_min, the stress is even smaller, primarily due to the small shear modulus, as discussed in the previous section.

Fig. 9b further compares the energy–strain relation for ε_xz, ε_yz and ε_min. Since ε_yz does not align well with the OA direction, it requires a high penalty energy barrier in general. In contrast, both ε_xz and ε_min correspond to the direct deformation path along OA and OB, and hence they require a lower energy barrier to cross. Indeed, we found that ε_min is the preferred deformation path upon heating. On the other hand, the activation barrier of ε_xz is about the same as that of ε_min, suggesting that the mechanical work of ε_xz could potentially trigger a phase transition. This may be interesting for future experimental studies.

Finally, the barrier values are ∼0.30 meV per atom for ε_yz, and ∼0.20 meV per atom for ε_xz and ε_min. These values are about 3 times larger than the barriers found in γ-surface calculation. This discrepancy is expected, as the γ-surface setup assumes a collective slip of one entire plane, while shear deformation involves the continuous slipping of multiple planes.

IV. Concluding remarks and outlook

In this work, we present a direct molecular dynamics simulation to investigate the phase transition mechanism in the model system HMB. Our results demonstrate that direct MD simulation of heating and cooling can accurately reproduce both the transition temperature and hysteresis loop observed in previous experimental studies.

Through analysis of the MD trajectories, we found that a molecular slip mode along (001)[1̄10] in the rotated supercell is the key mechanism triggering the phase transition within a critical temperature range. Further γ-surface energy analysis confirms that this mode has a low energy barrier that can be activated by thermal vibration around 380 K. This result is also supported by the observed continuous softening of shear modulus at the (001) plane, in particular along the [1̄10] direction. Collectively, these results highlight that the integration of different atomistic modeling techniques can provide invaluable insights into the mechanisms underlying martensitic phase transitions in organic crystals.

For the case of HMB, it is evident that the crystal packing with large molecular layer separation allowing for slip or shear at low activation barriers is fundamental to enabling this martensitic phase transition. Indeed, this class of molecular crystal has been widely recognized due to its exceptional mechanical properties.⁴⁸

Through this work, we have developed an in-house pipeline integrating several key modules: (1) force field generation; (2) crystallographic model set up; (3) scheduled heating/cooling or mechanical loadings in the context of MD simulation. Additionally, we have developed multiple scripts performing molecular plane analysis and property calculation, most of which are entirely general for the study of other systems which may exhibit similar structural behaviors.

Our studies on HMB, along with recent studies,^36,46 demonstrate that generic force fields such as GAFF are sufficient for qualitative simulations of phase cycling. Additionally, atomistic simulations offer quantitative insights – such as energy barriers and elastic properties – that help predict the feasibility of phase transitions driven by external stimuli, such as heat or mechanical work.

Looking ahead, we will further automate this pipeline and make it publicly available, facilitating streamlined data collection and enabling comparative analyses of experimentally reported mechanically flexible organic crystals. Ultimately, our goal is to enable large-scale screening prior to experimental validation, accelerating the discovery of novel flexible organic materials.

Data availability

The source code, instructions, as well as scripts used in this study are available in https://github.com/zfahim57/2025-HMB.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This research was sponsored by the NSF (DMR-2410178). The computation were performed at the Anvil cluster of Purdue University through allocation TG-MAT230046 from the ACCESS program, which is supported by U.S. NSF grants (2138259, 2138286, 2138307, 2137603 and 2138296).

References

1. K. Otsuka and C. M. Wayman, Shape memory materials, Cambridge University Press, 1999.
2. K. Otsuka and X. Ren, Physical metallurgy of Ti–Ni-based shape memory alloys, Prog. Mater. Sci., 2005, 50(5), 511–678.
3. I. Mihálcz, Fundamental characteristics and design method for nickel-titanium shape memory alloy, Period. Polytech. Mech. Eng., 2001, 45(1), 75–86.
4. B. Winzek, S. Schmitz, H. Rumpf, T. Sterzl, R. Hassdorf and S. Thienhaus, et al., Recent developments in shape memory thin film technology, Mater. Sci. Eng., A, 2004, 378(1–2), 40–46.
5. S. Miyazaki and A. Ishida, Martensitic transformation and shape memory behavior in sputter-deposited TiNi-base thin films, Mater. Sci. Eng., A, 1999, 273, 106–133.
6. Y. Fu, H. Du, W. Huang, S. Zhang and M. Hu, TiNi-based thin films in MEMS applications: a review, Sens. Actuators, A, 2004, 112(2–3), 395–408.
7. K. Oikawa, T. Ota, T. Ohmori, Y. Tanaka, H. Morito and A. Fujita, et al., Magnetic and martensitic phase transitions in ferromagnetic Ni–Ga–Fe shape memory alloys, Appl. Phys. Lett., 2002, 81(27), 5201–5203.
8. P. A. Olsson, P. Hyldgaard, E. Schröder, E. P. Jutemar, E. Andreasson and M. Kroon, Ab initio investigation of monoclinic phase stability and martensitic transformation in crystalline polyethylene, Phys. Rev. Mater., 2018, 2(7), 075602.
9. E. Brown, D. Dattelbaum, D. Brown, P. Rae and B. Clausen, A new strain path to inducing phase transitions in semi-crystalline polymers, Polymer, 2007, 48(9), 2531–2536.
10. G. R. Desiraju, Supramolecular synthons in crystal engineering—a new organic synthesis, Angew. Chem., Int. Ed. Engl., 1995, 34(21), 2311–2327.
11. C. Sutton, C. Risko and J. L. Bredas, Noncovalent intermolecular interactions in organic electronic materials: implications for the molecular packing vs electronic properties of acenes, Chem. Mater., 2016, 28(1), 3–16.
12. S. K. Park, J. H. Kim and S. Y. Park, Organic 2D optoelectronic crystals: charge transport, emerging functions, and their design perspective, Adv. Mater., 2018, 30(42), 1704759.
13. E. Ahmed, D. P. Karothu and P. Naumov, Crystal adaptronics: mechanically reconfigurable elastic and superelastic molecular crystals, Angew. Chem., Int. Ed., 2018, 57(29), 8837–8846.
14. S. Takamizawa and Y. Takasaki, Versatile shape recoverability of odd numbered saturated long-chain fatty acid crystals, Cryst. Growth Des., 2019, 19(3), 1912–1920.
15. T. Sasaki, S. Sakamoto and S. Takamizawa, Twinning organosuperelasticity of a fluorinated cyclophane single crystal, Cryst. Growth Des., 2019, 19(10), 5491–5493.
16. T. Mutai, T. Sasaki, S. Sakamoto, I. Yoshikawa, H. Houjou and S. Takamizawa, A superelastochromic crystal, Nat. Commun., 2020, 11(1), 1824.
17. P. Naumov, S. Chizhik, M. K. Panda, N. K. Nath and E. Boldyreva, Mechanically responsive molecular crystals, Chem. Rev., 2015, 115(22), 12440–12490.
18. S. C. Sahoo, S. B. Sinha, M. Kiran, U. Ramamurty, A. F. Dericioglu and C. M. Reddy, et al., Kinematic and mechanical profile of the self-actuation of thermosalient crystal twins of 1, 2, 4,5-tetrabromobenzene: A molecular crystalline analogue of a bimetallic strip, J. Am. Chem. Soc., 2013, 135(37), 13843–13850.
19. M. K. Panda, T. Runčevski, S. Chandra Sahoo, A. A. Belik, N. K. Nath and R. E. Dinnebier, et al., Colossal positive and negative thermal expansion and thermosalient effect in a pentamorphic organometallic martensite, Nat. Commun., 2014, 5(1), 4811.
20. A. Colin-Molina, D. P. Karothu, M. J. Jellen, R. A. Toscano, M. A. Garcia-Garibay and P. Naumov, et al., Thermosalient amphidynamic molecular machines: Motion at the molecular and macroscopic scales, Matter, 2019, 1(4), 1033–1046.
21. A. Khalil, D. P. Karothu and P. Naumov, Direct quantification of rapid and efficient single-stroke actuation by a martensitic transition in a thermosalient crystal, J. Am. Chem. Soc., 2019, 141(8), 3371–3375.
22. Y. Diao, K. M. Lenn, W. Y. Lee, M. A. Blood-Forsythe, J. Xu and Y. Mao, et al., Understanding polymorphism in organic semiconductor thin films through nanoconfinement, J. Am. Chem. Soc., 2014, 136(49), 17046–17057.
23. H. Chung, S. Chen, B. Patel, G. Garbay, Y. H. Geerts and Y. Diao, Understanding the role of bulky side chains on polymorphism of btbt-based organic semiconductors, Cryst. Growth Des., 2020, 20(3), 1646–1654.
24. S. Q. Su, T. Kamachi, Z. S. Yao, Y. G. Huang, Y. Shiota and K. Yoshizawa, et al., Assembling an alkyl rotor to access abrupt and reversible crystalline deformation of a cobalt (II) complex, Nat. Commun., 2015, 6(1), 8810.
25. F. Borbone, A. Tuzi, A. Carella, D. Marabello, S. L. Oscurato and S. Lettieri, et al., High-temperature reversible martensitic transition in an excited-state intramolecular proton transfer fluorophore, Cryst. Growth Des., 2019, 19(11), 6519–6526.
26. D. Das, E. Engel and L. J. Barbour, Reversible single-crystal to singlecrystal polymorphic phase transformation of an organic crystal, Chem. Commun., 2010, 46(10), 1676–1678.
27. K. Lonsdale, The structure of the benzene ring, Nature, 1928, 122(3082), 810.
28. Y. Saito, X-ray Investigation of the Thermal Transitions in Some Molecular Crystals, X-RAYS, 1952, 7(1–4), 9–24.
29. Y. V. Mnyukh and N. Petropavlov, Polymorphic transitions in molecular crystals-I. Orientations of lattices and interfaces, J. Phys. Chem. Solids, 1972, 33(11), 2079-IN3.
30. Y. V. Mnyukh and N. Panfilova, Polymorphic transitions in molecular crystals—II. Mechanism of molecular rearrangement at ‘contact’ interface, J. Phys. Chem. Solids, 1973, 34(2), 159–170.
31. Y. V. Mnyukh, N. Panfilova, N. Petropavlov and N. Uchvatova, Polymorphic transitions in molecular crystals—III.: Transitions exhibiting unusual behavior, J. Phys. Chem. Solids, 1975, 36(3), 127–144.
32. Y. V. Mnyukh, Polymorphic transitions in crystals: Nucleation, J. Cryst. Growth, 1976, 32(3), 371–377.
33. Y. V. Mnyukh, Polymorphic transitions in crystals. V. Kinetics of crystal growth in crystal medium, J. Cryst. Growth, 1977, 38(2), 284–284.
34. L. Li, P. Commins, M. B. Al-Handawi, D. P. Karothu, J. M. Halabi and S. Schramm, et al., Martensitic organic crystals as soft actuators, Chem. Sci., 2019, 10(31), 7327–7332.
35. S. Fredericks, K. Parrish, D. Sayre and Q. Zhu, PyXtal: A Python library for crystal structure generation and symmetry analysis, Comput. Phys. Commun., 2021, 261, 107810.
36. P. A. Santos-Florez, S. Hattori and Q. Zhu, Bending deformation driven by molecular rotation, Phys. Rev. Res., 2023, 5, 033185.
37. D. A. Case, H. M. Aktulga, K. Belfon, I. Ben-Shalom, S. R. Brozell and D. S. Cerutti, et al., Amber 2021, University of California, San Francisco, 2021.
38. J. Wang, R. M. Wolf, J. W. Caldwell, P. A. Kollman and D. A. Case, Development and testing of a general amber force field, J. Comput. Chem., 2004, 25(9), 1157–1174.
39. A. Jakalian, B. L. Bush, D. B. Jack and C. I. Bayly, Fast, efficient generation of high-quality atomic charges. AM1-BCC model: I, J. Comput. Chem., 2000, 21(2), 132–146.
40. A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown and P. S. Crozier, et al., LAMMPS - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales, Comput. Phys. Commun., 2022, 271, 108171.
41. S. Adelman and J. Doll, Generalized Langevin equation approach for atom/solid-surface scattering: General formulation for classical scattering off harmonic solids, J. Chem. Phys., 1976, 64(6), 2375–2388.
42. W. G. Hoover, Canonical dynamics: Equilibrium phase-space distributions, Phys. Rev. A: At., Mol., Opt. Phys., 1985, 31(3), 1695.
43. R. Gaillac, P. Pullumbi and F. X. Coudert, ELATE: an open-source online application for analysis and visualization of elastic tensors, J. Phys.: Condens. Matter, 2016, 28(27), 275201.
44. S. K. Park and Y. Diao, Martensitic transition in molecular crystals for dynamic functional materials, Chem. Soc. Rev., 2020, 49(22), 8287–8314.
45. S. Takamizawa and Y. Miyamoto, Superelastic organic crystals, Angew Chem., 2014, 126(27), 7090–7093.
46. M. B. Al-Handawi, P. Commins, A. S. Dalaq, P. A. Santos-Florez, S. Polavaram and P. Didier, et al., Ferroelastic ionic organic crystals that self-heal to 95%, Nat. Commun., 2024, 15(1), 8095.
47. While it was mentioned that the direction is [1̄10], in LAMMPS simulation, we applied ”fix 1 all deform 1 xz erate −1.25 × ${srate} yz erate ${srate} units box remap x flip no” for the rotated cell as described in Table 1 .
48. C. M. Reddy, R. C. Gundakaram, S. Basavoju, M. T. Kirchner, K. A. Padmanabhan and G. R. Desiraju, Structural basis for bending of organic crystals, Chem. Commun., 2005, 3945–3947.

---