Probing anisotropic mechanical behaviour in carbamazepine form III

Nanoindentation measurements in single crystals of carbamazepine form III show that the (020) face is stiffer and harder than the (002) and (101) faces. AFM imaging and molecular simulations reveal that the (020) plane is the most likely slip plane.


Introduction
Beyond their required biological activities, active pharmaceutical ingredients (APIs) need to have appropriate further physicochemical properties that enable their development into effective drug products. For example, APIs need to have suitable aqueous solubilities and kinetics of dissolution, 1 be able to grow from solution into adequately sized crystals with good filterability 2 and have the required mechanical properties for tabletting. If the physical properties of an API are poor, attempts to improve them are often sought through the design of alternative solid forms such as cocrystals, 3 hydrates, [4][5][6] or polymorphs. 7 This approach requires the development of new solid forms often involving different chemical compositions, which may require new and costly solid form screenings together with corresponding stability studies. Another route to modify properties is through the design, addition and fine tuning of various processes used for the production of the API. For example, sonocrystallisation, 8 temperature cycling 9 or appropriate use of milling 10 are processes shown to enable modifications of particle size and shape distributions in crystalline powders which in turn can impact their filtration behaviour 11,12 or mechanical properties. 8 Tablets contain a mixture of crystalline powders of the API and excipients. This mixture is referred to as the formulation and the exact proportion of the powder components is usually established on a trial-and-error basis working around the mechanical behaviour of the API powders so as to achieve a powder mixture that can be tabletted. Tablets formulated directly by compacting the API in isolation (without excipients) have the advantage of being smaller, not needing the formulation step and thus being cheaper to produce. Ultimately, the mechanical properties of the material API in powder form are primarily dictated by the mechanical properties of its single crystals. In that context, the study of the mechanical properties of single crystals directly using nanoindentation 13,14 offers the opportunity to determine an API's mechanical behaviour at early stages of drug development (when the material is only available in small quantities) and to probe a number of other solid-state phenomena upon exposing the material to tabletting (i.e. stressed induce phase transformations, 15

dehydration 4,16 …).
The early acquisition of this information would provide extra tools to influence process design and control, allowing improvements to the drug's mechanical behaviour prior to its transition to formulation and tabletting.
Molecular crystals, particularly those with low symmetry forms, are well known to be highly anisotropic, with their properties varying along different crystal axes and faces. Growth rates, 2 dissolution 17 and mechanical properties 18 are, amongst others, properties known to be face and direction dependent. Of particular interest, nanoindentation enables the evaluation of anisotropy of mechanical properties in molecular crystals, when combined with basic face indexing pre-experiments. It can probe specific crystallographic surfaces over penetration depths <1 μm and requires an exposed area of a few μm 2 . For this, however, large single crystals (hundreds of μm for at least two dimensions) are required so that repeated measurements on different {hkl} faces are possible. Most molecular crystals only expose one perhaps two {hkl} families of faces, hence it remains very challenging to perform nanoindentation experiments along different crystal directions. Despite these difficulties, the number of reports on anisotropic mechanical behaviour in molecular crystals measured with nanoindentation has increased considerably in recent years. Examples of these include works on saccharin, 19 L-alanine, 20 sucrose, 21 difluoroboron avobenzone, 22 piroxicam 15 or a series of amino-acids 23 amongst others. Nanoindentation allows for the determination of a reduced modulus E R (related to the modulus of elasticity E (hkl) ≈ E R,(hkl) (1 − ν), [24][25][26] with ν the Poisson's ratio) and hardness H of materials.
Here, we have studied the mechanical properties of carbamazepine (CBZ) form III (Fig. 1) along a number of crystallographic directions. CBZ, a drug used to treat epilepsy, has an interesting solid form landscape with five polymorphs [27][28][29][30][31] and a wide range of cocrystals. 32 We have chosen the stable form III 28 for our nanoindentation studies since it readily crystallises as equant large crystals under various crystallisation conditions. Growth of appropriately sized crystals of CBZ-III from methanol allows for nanoindentation of three distinct faces: (020), (101) and (002). Nanoindentation combined with imaging analysis reveals distinct mechanical behaviour in the three directions. These experiments are supported by computer simulations to aid our understanding of plastic deformation of CBZ-III along the three directions studied. Overall, the study provides an enhanced understanding of the anisotropic mechanical behaviour in form III carbamazepine from a combined experimental and computational perspective.

Sample preparation
Slurries of carbamazepine in methanol (anhydrous, 99.8%, Sigma Aldrich) were heated to 40°C. Solutions were then filtered and transferred into petri dishes which were left to cool (naturally from 40°C to room temperature, ∼21°C) and slowly evaporate under ambient conditions over 48 to 96 hours. The harvesting of the CBZ-III single crystals from their mother liquor was performed as follows. Due to the very high solubility of CBZ in methanol, a layer of hexane (HPLC, >95%, Fisher Chemical) was poured gently on top of the methanolic solution. The two solvents being immiscible, it allowed for the removal of the saturated solution from the crystals' surface while harvesting them with tweezers from solution, preventing secondary nucleation. Handling of crystals was always done with the aid of PEEK tweezers. Crystals retained for the nanoindentation experiments had dimensions of 0.5 × 0.5 × 0.25 mm, or larger. After face identifications, crystals were glued in the correct orientation using cyanoacrylate glue.

Face indexing
Representative single crystals of CBZ-III grown from methanol were face indexed by means of a goniometer mounted on a single crystal X-ray diffractometer (XCalibur-2, Rigaku Oxford Diffraction) with a MoKα X-ray source (λ = 0.71073 Å). Photographs of the single crystals were taken before each X-ray diffraction measurement (whilst rotating the crystals) and were used for subsequent face indexing. The unit cell parameters and the space group were determined and verified to match those of CBZ-III as reported in the Cambridge Structural Database 33 (CSD) refcode CBMZPN02. 28 Thus, the Miller indices used throughout the manuscript refer specifically to the unit cell parameters defined in this refcode. Indexing of the faces was carried out using the software CrysAlisPro (v. 171.40.14d, Rigaku Oxford Diffraction). 34 Nanoindentation Experiments were carried out using a TI 950 Triboindenter (Bruker Nanomechanical Test Instruments, Minneapolis, USA) equipped with a diamond Berkovich probe. Indents were carried out on selected single crystals whose crystal shape clearly matched those identified during face indexing. Indentations were conducted under load-control mode using a unique trapezoid load function which is presented in the ESI. † Crystals of carbamazepine were large enough to allow for them to be firmly glued on a glass slide prior to the indentation experiments. 10 crystals were indented per crystal face tested, 9 indents were carried out per crystal and the distance between individual indents was set to 20 μm. Reduced modulus and hardness were extracted from the indentation curves as stated in ISO 14577. 35 Calibration of the tip was made using a PS-1 sheet sample (E R = 2.91 GPa, VPG Micro-measurements, Malvern, USA), which was indented in between each single crystal tested to verify for possible tip contamination and calibration errors. As previously described in detail, 14 prior to indentation we investigated the optimal maximum displacement, load rate and holding times which were fixed to 5 second load, 2 second hold and 5 seconds unload segments. The maximum load P Max used was fixed to 3 mN as no further changes in the values of E R and H were observed for higher loads. 14

Imaging
Crossed polarised light was used to identify single crystals of good quality using an Axioplan2 Zeiss Microscope in transmission mode. Imaging of the indents was carried out using an atomic force microscope (AFM, Multimode8, Bruker), using PeakForce QNM in Air mode, with a TESP-V2 probe (k = 42 N m −1 ).

Computation
All simulations were done with Materials Studio 2019 (ref. 36) using the COMPASS-II 37 forcefield with its own atomic charges, known to be a reliable forcefield for the simulation of molecular crystals. 38 Energy values of all simulations were converted to kJ mol −1 and reported per molecule of CBZ.
Structure retrieval, geometry optimisation and attachment energy calculations. The crystal structure of CBZ form III (CBMZPN02) was retrieved from the CSD. Using the Forcite module, the crystal structure was geometry optimised allowing all atomic positions to relax while keeping the unit cell parameters constant to the observed room temperature unit cell parameters. The COMPASS-II optimised structure served as the initial model for the rest of the computations. Attachment energy calculations were then computed with the same model.
Simulation of compression. CBZ unit-cell parameters were compressed progressively in three different crystallographic directions. CBZ-III compressed cells were then geometry optimised with fixed unit cell parameters in two sets of simulations: treating the molecules as flexible and rigid. For the compressions, the required unit cell axis (or axes) was progressively reduced by 0.05 Å from its original experimental value whilst keeping the other cell axes constant, until the corresponding lattice energy became positive. The axes compressed were the b-axis, the c-axis and both the a and c-axis together to simulate compression upon indentation on the (020), (002) and (101) faces respectively.
Molecular dynamics (MD). Slip planes and slip directions of CBZ-III were studied by means of shearing MD simulations. Supercells of 96 molecules of CBZ-III were generated from the relaxed unit cell. A MD run at 298 K was first performed for 50 ps in the NVT ensemble to allow for the system to equilibrate. The system was then run for a further 50 ps in the NPT ensemble. Shearing simulations were then carried out in different directions (see ESI † and Results section) at 298 K using constant shearing rates of 5.10 −4 ps −1 for 1500 ps (NVT). Two simulations were done per shearing direction whereby the shear force was applied along one of the supercell axes, maintaining one of the two planes containing the sheared axis constant.
Mechanical properties. After the 50 ps-NVT and 50 ps-NPT MD equilibration runs, the resulting supercell was run for an additional 500 ps-NPT at 298 K. Every 5 ps, a series of strains were applied to the simulation cell. The stiffness and compliance matrices were then computed using the constant strain method. The components of the final stiffness and compliance matrices are averages of the entire simulation run. The Young's modulus was then computed in 3D and projected onto the (100), (010) and (001) planes using the ElAM computer program 39 and the values along the desired directions retrieved and presented in 2D plots in the Results section. Poisson ratios were also extracted from such simulations.
Surface rugosities. Surface rugosities were calculated for a number of (hkl) faces using the Bryant et al. method 40 as implemented in the CSD Python API. 33 Positive values indicate separation between planes in Å whilst negative values indicate layer intercalation.

Experimental results
3.1.1. CBZ-III crystal habit, face indexing and crystal structure. Thick crystals of CBZ form III were obtained from methanol. Fig. 2 shows the face indexing of a representative single crystal. From methanol, CBZ-III exhibits a large dominant (020) face and, however less significant, still large (002) and (101) lateral faces. When viewed under the optical microscope, the shape of these two latter faces were identical. Hence to identify each face with certainty prior to sample mounting for indentation, the angle between the inclined top section of the crystal and the (002) and (101) faces were measured under the microscope and compared with the crystallographic data (Fig. 2, face indexing). The angles between the faces measured from the micrographs of 23 different crystals were 116.4°± 0.8 and 123.4°± 1.2, in excellent agreement with the theoretical crystallographic values. The dominance of the (020) face was further confirmed by preferential orientation powder-X-ray diffraction (PXRD, ESI †).
The crystal packing of CBZ-III along the three directions indented is presented in Fig. 3. Hydrogen bonded dimers (R 2 2 (8) amide⋯amide) and aromatic interactions are highlighted in Fig. 3a (black dashed circles and black arrows respectively). The lattice energy of CBZ-III U LATT was calculated to be −145.6 kJ mol −1 with contributions of around 20% of electrostatic interactions and 80% of van der Waals interactions (see ESI †). Dimers are oriented very closely along [010] (Fig. 3a) whilst van der Waals interactions are dominated mostly from the aromatic interactions between neighbouring aromatic rings of CBZ.
3.1.2. Nanoindentation experiments. Crystals were mounted and indented on their three dominant faces (020), (002) and (101). Thirty crystals were mounted, ten on each of the faces. Fig. 4 shows one typical load-depth curve for each face indented. In all curves, multiple small bursts occurred upon loading (commonly referred to as "pop-ins"). The intensity of pop-ins observed on the loading section of nanoindentation experiments of molecular single crystals have been correlated with relevant crystallographic parameters such as specific d hkl values or particular unit cell lengths by other groups. 19,41,42 We attempted to establish a similar correlation for carbamazepine, with no success (analysis presented in ESI †). The maximum depth of penetration h Max reached is lower for the (020) faces compared to the (002) and (101) faces (Fig. 4).
The associated averages and standard deviations of E R and H obtained from these nanoindentation experiments are given in Table 1, whilst associated histograms of these datasets are presented in ESI. † The reduced modulus measured on the (020) faces is higher than those measured on the (002) and (101) by 26% and 29% respectively. Similarly, the hardness measured on the (020) face is higher Fig. 2 Micrographs of a single crystal of CBZ form III grown from methanol, with its associated face indexing.   than that of the two other faces, with an increase of approximately 20%. Differences between the (002) and (101) faces are less prominent and overall averages of both E R and H on these faces are extremely similar. 3.1.3. Imaging of indents. Atomic force microscopy (AFM) was carried out on a number of indents in order to assess the degree of plastic deformation around the indents. Images revealed that the three surfaces responded very differently upon indentation. A typical image of an indent on the hardest (020) surface in Fig. 5a shows the indent is well defined as a near-perfect equilateral triangle. No pile-ups, sink-ins or cracks were observed around indent edges. However, in Fig. 4, the load-depth trace of an indent carried out on the (020) face clearly contains pop-ins, which cannot be visually related to the formation of radial cracks because none were observed in Fig. 5a. These pop-ins are probably related to plastic deformation through dislocation nucleation, these dislocations propagate into the material to accommodate the volume displaced to generate the indent. The second and third AFM images ( Fig. 5b and c) show indents carried out on the softer (002) and (101) faces respectively. Plastic deformations of several nature are observed in both cases with evidence of slips, radial cracks, small pile-ups and sink-ins. On the (002) face, slip lines run along [100] and individual cracks propagate along a direction close to (but not quite) [010] with cracks having a zig-zag shape. On the (101) face, similar deformation features are observed, however it was not possible to unequivocally determine the direction of the slip bands and radial cracks. We note that cracks on the (101) face also propagated alternatively along two directions, as seen in Fig. 5c. The presence of small pile-ups and sink-ins along some indent edges is visible in the cross-sections shown in Fig. 5b′ and c′. Such dramatic differences in the plastic response of CBZ crystals across different crystallographic faces are consistent with the differences in the values of E R and H in Table 1.
We note that whilst indentation of the (002) and (101) faces resulted in cracks at the working P Max of 3000 μN (which may result in an artificial reduction of the derived stiffness or hardness), 43,44 we did not observe a significant impact of these cracks on the resulting E R or H values calculated at different maximum loads (Fig. 6). If cracks were to have a significant effect on these derived values, an increase in load should result in more cracks and thus a significant lowering of the calculated E R and H values. As discussed previously, 14 the maximum load P Max used in nanoindentation measurements needs to be above a minimum value below which the data becomes inconsistent due to calibration issues. In this case, the plateau for acquisition of consistent data is reached in all three faces at around 2000 μN. At loads above 2000 μN, the derived values of E R and H as a function of maximum load remain constant within the expected experimental variability for these measurements. For example at P Max of 3000 μN, the derived mean E R value for the (101) face on 90 independent indents on 9 different crystals is 10.74 GPa with a measured  (Fig. 6) lie within these ranges. This suggests that cracking, in this case, has little impact on the derived stiffness and hardness values. We also note that since the pile-ups and sink-ins observed in the (002) and (101) indents were small, no corrections were applied to the indentation area calculations.

Computational results
3.2.1. Elasticity. We calculated the elastic properties of CBZ-III using the COMPASS-II forcefield as implemented in Materials Studio. Calculations were made twice: 1) with the unit cell parameters set to room temperature values (embedded in the refcode CBMZPN02): RT-computed and 2) optimising them (≈ optimisation at 0 K): 0 K-computed. These results were obtained from the computed matrix of elastic compliance constants (presented in ESI †). 39 An average of both was also calculated and results along the indentation directions are shown in Table 2. We note that the RT-computed data are ≈20% below the experiment values whilst the 0 K-computed data are ≈20% above. The average computed values of Young's modulus E, however, match very well those measured experimentally. These were calculated from experimental reduced moduli E R using the computed Poisson's ratio ν for each associated compressed direction (averaged in a similar way over both "0 K" and "RT" computed data). All values are given in Table 2. Fig. 7 shows the 2D projection of the "averaged-E" matrix on the (100), (010) and (001) planes, in which the distance of each data point from the centre of the graph corresponds to the absolute value of E computed for a given crystallographic direction. The highest computed value read is 19.0 GPa along [−302] and the lowest is 9.1 GPa along the [301] as shown in Fig. 7b. If we define here the anisotropy ratio (AR) as the ratio of the largest over the smallest value of a property, hence, for the Young's modulus, AR COMP-AVG = 2.1. From Fig. 7a, we can extract E Z = E (002) = 9.9 GPa and E Y = E (020) = 12.8 GPa, thus the calculations confirm that E (002) < E (020) as observed experimentally. Moreover, the computed value for the direction normal to (101) (i.e. direction of indentation, ≈[301] as explained in ESI †) can be read in Fig. 7b as E (101) = 9.1 GPa. Again, E (101) < E (020) whilst E (101) is slightly lower in value than E (002) as observed experimentally. We note that the computed values for E are close to the measured values and predict the correct anisotropy of mechanical behaviour in CBZ-III, which is encouraging given the simplicity of the model (a classic COMPASS-II forcefield).
3.2.2. Molecular flexibility. As seen in Fig. 1b, CBZ molecules have a butterfly like shape with the two aromatic rings (butterfly wings) found in its form III at an angle θ 0 = 130.7°. Simulations of the deformation of the CBZ crystal structure along the three nanoindentation directions were carried out to study the influence of molecular shape on mechanical properties (Fig. 8a). Upon compression along [010], [001] or [301] directions (directions of indentation), intramolecular and intermolecular bonds will deform elastically prior to plastic deformation. Fig. 8b shows how the lattice energy U LATT increases upon such compressions as a function of relative volume reduction (V/V 0 where V 0 is the volume of the uncompressed cell). The compression simulations were carried out in two ways: treating   Fig. 8b) or flexible (solid lines Fig. 8b). For any value of relative volume V/V 0 , U LATT is seen to be lower for the compressions with the flexible molecule assumption, hence the system is more stable. In Fig. 8c, the deformation of the "wings" in the molecule is shown by displaying the change in angle between these two rings upon compression (using V/V 0 as the measure) along three indentation directions. The molecular deformation along the three directions is very different with compression along the [010] axis resulting in the opening of the butterfly wings (θ 1 > θ 0 ), whilst compression along the [001] and [301] directions resulting in their closing (θ 1 < θ 0 ). In Fig. 8b, the lattice energy functions can be conveniently fitted to a second order polynomial function. The modulus of elasticity along a crystal direction is theoretically proportional to the second derivative of the corresponding function, when  assessed at U LATT = U LATT,MIN . Full details of the equation linking E and U LATT , obtained for these derivatives, including numerical values and statistical measures of fit, are reported in ESI. † Here, we draw one simple conclusion: the second derivative is always lower for the potential energy functions of flexible molecules as compared to those for rigid molecules, which indicates that flexibility makes the material less stiff in CBZ-III.
3.2.3. Shear deformation. Given the difficulty in determining the (hkl) indices of the slip planes from analysis of the specimen surfaces after indentation of the faces, simulations of shear deformation of the crystals were performed to identify the influence of crystal structure on slip. The slip behaviour of planes in CBZ-III was computed by carrying out shearing simulations on supercells of 96 molecules. Extended details of the method are presented in the ESI. † For each indentation plane, two shearing simulations were performed, and for each simulation, pure shear was applied along the indentation direction, whilst keeping one plane fixed (which contains that direction) at one end of the simulation box whilst shearing it at the other end of the box. Thus, for example, for the indentation of the (002) plane, shearing along the [001] was studied in two simulations: i) fixing the (010) plane at one end of the box whilst shearing it at the other, and ii) fixing the (100) plane at one end of the box whilst shearing it at the other. The first of such simulations is shown in Fig. 9 which results in the shear of the (020) plane. Thus, the simulation reveals which plane will shear upon applying a shear force on the studied nanoindentation directions. An interesting result here is the sequential slip of the (020) plane, with slips occurring by a third of the c-unit cell length at a time (∼4.6 Å, this sequence is illustrated in Fig. 9). The first slip results in the breaking of CBZ⋯CBZ hydrogen bonded dimers (amide⋯amide R 2 2 (8) hydrogen bonds). Two further slipping bursts result in reformation of the dimers. This could possibly be interpreted as indicating that the slip mechanism occurs by a sequence of partial dislocations.
The results from the shearing simulations are presented in Table 3. At least one slip event occurred over the 1500 ps simulation time along the (020) plane in four simulations, which include all three indented faces. A single simulation showed slip along the (002) plane. One of the simulations showed no slip events prior to structural collapse (loss of long-range order in the supercell). Slip events are abrupt, occurring over a short time <10 ps. It is notable that only two simulations revealed repeating slip events (simulations 3 and 6, Table 3), indicating (020) is the most likely slip plane. These simulations only show repetitive crystal slips when indentation is carried out on the (002) or (101) face. We further note that the intensity of the shear (distance made by the (020)-molecules per slip event) is not an integer of a unit cell axis (1/3 of c-axis ∼4. 6   crystals, and in particular to identify potential slip planes. 45 Initially, the energy of attachment U ATT was used to predict weak planes along which slip systems can be most easily activated as it is proportional to the energy required to slice one (hkl) plane, but this hypothesis had limited success. 46 A second approach, developed by Bryant and co-workers, 40 is to compute the rugosity and the degree of interpenetration between the main crystallographic planes, assuming that slip is easier along low rugosity planes. Here the rugosity parameter is positive if the planes are not interlocked, thus indicating separation between the planes, and negative if the planes are interlocked. A third approach is to compute the energy between molecules in the crystals to reveal periodic bond chains (PBCs) also known as energy frameworks. 47,48 Bonds of different strengths are represented as tubes between the centroid of molecules within a supercell, revealing the strength and dimensionality of intermolecular interactions. Materials Studio with the forcefield COMPASS-II was used here to compute the three predictors in CBZ-III, the results are shown in Table 4 and Fig. 10.
One may state that the (011) plane is the most likely slip plane based on attachment energy computations or that the (020) and the (101) planes are the most likely slip planes based on the surface rugosities (these are the smoothest surfaces). Table 4 shows thus that the rugosity predictor seems to be in agreement with the experimental observations and the MD-shear simulations for the (020) face, however it shows that the attachment energy model fails.
We note that these are simple predictive tools and should be used with care. In Fig. 10, PBCs in CBZ-III are shown in a 3 × 2 × 2 supercell (only interactions that are more stabilising than −20 kJ mol −1 shown). This reveals a 3D network in CBZ-III, made up of hydrogen bonds (blue tubes) corrugated with aromatic interactions of various strengths (purple and red tubes). Such 3D-network materials have been reported to i) exhibit brittle mechanical responses 47 and ii) possess higher mechanical properties than their 2D counterparts. 49 Concerning this latter point, values reported in Table 1 are significantly higher than those reported for smaller 2Dnetwork molecules such as in benzoic acid derivatives, even if systems with higher mean molecular volumes like carbamazepine should have reduced mechanical properties. 49

Discussion
We have studied CBZ-III experimentally with nanoindentation and have interpreted its mechanical behaviour using a number of simulation methods using the COMPASS-II forcefield. For the experimentation, we note that appropriate crystal growth is needed to obtain morphologies allowing nanoindentation of multiple crystal faces and thus assessment of the mechanical properties along independent directions for this system. Here, we found that growth from methanol provided a suitable chunky block morphology allowing indentation on three distinct faces. Table 3 Results of molecular dynamics shearing simulations, including three shearing directions and two planes sheared per shearing direction, on the crystallographic planes indented. When slip was observed, the plane and the direction along which it was observed is reported   For the computations, since calculations of mechanical properties (i.e. stiffness matrix) require the computation of the second derivatives of the system forces, the model describing the force/potential needs to be very accurate. A very accurate method based on DFT-d, however, would require of significant computational time and would not allow for MD or shearing simulations involving supercell models. It is thus challenging to find an accurate method that is also affordable computationally for the prediction of mechanical properties in molecular crystals. 45 Here the COMPASS-II model seemed to perform relatively well when both resultswith and without unit cell parameter optimisationof E were averaged. Further to this, the model allowed for the PBC calculations which revealed that CBZ-III is made of strong 3D continuous interactions (absolute values higher than 20 kJ mol −1 ), so the material was expected to be harder than other comparable 2D or 1D interacting systems, and the identification of soft/slip planes was thus more complex. We also used this model to gain insights as to the mechanisms of inter and intra molecular deformations upon compression. Whilst there are a number of studies looking at the impact of intermolecular interactions on the mechanical properties of molecular crystals, to the best of our knowledge this is the first study which has studied the contribution of molecular flexibility.
Combining the simulation results and experimental data, we can assess the extent of anisotropy of the mechanical properties in CBZ-III in various ways. We extend here the definition of the anisotropy ratio (AR) as the ratio of the largest value of a property over the smallest for values measured/computed in a number of directions (i.e. E MAX /E MIN or H MAX /H MIN ). A large AR would indicate a highly anisotropic behaviour whilst an AR = 1 would indicate a perfectly isotropic material. CBZ-III crystals had experimental AR values for E and H equal to 1.3 and 1.2 respectively (using data from Table 1). These AR values are in good agreement with typical AR results measured experimentally and reported in the literature. In molecular crystals, these ARs range from 1 to ∼2 for Young's modulus and from 1 to ∼1.6 for hardness. 50 From this, CBZ-III appears to be at the middle-lower end range of ARs for molecular crystals. For comparison, in crystalline metals ARs of the Young's modulus were reported to range from 1.12 (magnesium) to 3.42 (zinc). 51 The anisotropy in mechanical behaviour was also manifest through the imaging of the indents on the three different crystal faces studied. This revealed that, upon indentation of the (020) face, neither significant surface displacement around the indent nor cracking was seen. However, indents on both the (002) and the (101) faces clearly displayed small material pile-up, slip traces and cracks. Our E R values calculated at different P Max and depth of penetrations (Fig. 6) show that for this system, the cracking seems to have no or little impact on the derived mechanical properties. This is also corroborated in the computations of E: the predicted stiffness for the studied faces match well the experimentally derived values. Interestingly, the faces that showed some degree of cracking, (002) and (101), were softer than the face that did not show cracking, (020). Normally, brittle faces are expected to be harder but we observed the opposite behaviour for CBZ-III.
It was difficult to accurately identify the slip planes from the images of the slip traces but, with the aid of simulations, these could unequivocally be assigned to the (020) planes. This result is surprising because i) the slip of (020) (Fig. 3a) requires of the breakage of hydrogen bonds and ii) (020) does not have the lowest energy of attachment (three other planes have a lower U ATT ). However, the associated rugosity of (020) predicts the highest probability of slip amongst the potential slip planes computed (Fig. 9). Hence, in this case it seems that interplanar rugosity is the most significant predictor of the primary slip plane.
To enable the slip of planes with negative rugosities, one can think of two alternatives. First, a large and local increase of the interplanar spacing d hkl (repulsive intermolecular interactions) would make the rugosity positive, enabling a slip or possibly a crack. A second option would require the modification of molecular flexibility (angles between intramolecular bonds, for example θ 1 in Fig. 8) upon external stress in such a manner that a negative rugosity becomes positive. Both alternatives would probably require significantly more energy than that required to break the hydrogen bonds between two (020) planes. We note that the three options discussed (slicing and slip of planes with positive rugosity, local enlargement of d hkl or modification of angles between intramolecular bonds) all require energy. To correctly identify slip planes in molecular systems, the challenge would be hence to quantify which one is the most favourable energetically. The shearing experiments do indeed predict easiest slip on the (020) plane even upon indentation of the (020) plane. This again confirms this plane contains the easiest/lowest energy deformation mechanism. However, upon indentation of (020), the (020) slip must be partly hindered because (020) is orthogonal to the indentation (shearing) direction, leading to higher values of hardness in this direction than in the other two. We note that the shearing simulations only afforded one slip event on the (020) simulations whilst several events occurred on the other faces for simulations of the same length of time.
Coming back to the rugosity calculations, planes with similar rugosities can indeed behave differently. For example, the rugosity calculations show that whilst (020) and (−101) faces have similar rugosities (+0.09 and +0.08 Å repectively), only the (020) plane serves as a slip plane as shown by the MD simulations. Here, the hydrogen bonds that need breaking for this slip to occur are important. In previous works, the breaking of hydrogen bonds has been acknowledged as playing an inhibitting role in slippage but it has not been quantified. 40 Here, we propose the calculation of what we define as the (hkl) plane hydrogen bond density (HB-density-(hkl)). For this, a 2D-unit cell on the plausible slip plane under consideration is calculated by identifying the shortest translations along two perpendicular directions to the plane direction. This is, crystallographically speaking, the 2D-unit cell on the plane of interest. The number of HBs that need breaking in this area for the plane to slip are then counted and the HB-density (number of bonds per unit of area) is calculated. Table 5 presents the rugosity data together with the HB-density for the most important planes. When seen together then, the combination of low rugosity (flat surface) and low HB-density-{hkl} (less energy required to break the HBs during slippage) does indeed unambigously reveal the (020) plane as the most likely slip plane in agreement with the MD simulations.
The simulation of unit cell compression also demonstrated that molecular deformation of CBZ is markedly different, depending on the deformation orientation. Compression of the hardest face, (020), leads to an opening of the CBZ butterfly wings whilst compression of the softer faces results in their closing. In both cases, this flexibility results in a softer material than that of an equivalent rigid molecule ( Fig. 6b and ESI †).
Finally, we note that whilst this is the first study of mechanical behaviour in CBZ-III single crystals upon nanoindentation, previous works by Rowe and Roberts reported the mechanical properties of CBZ-III powder compacts with different levels of porosities. However measured on powder compacts, they reported a value of 13.2 GPa for the Young's modulus of CBZ-III. 52,53 Their reported value is close to those we obtained using nanoindentation on different faces.

Concluding remarks
The large size of crystals obtained by crystallising carbamazepine form III from methanol enabled nanoindentation testing along more than one crystal orientation. Nanoindentation showed that the (020) face of CBZ-III is stiffer and harder than the two other faces tested, (002) and (101). Imaging the indents using AFM revealed different plastic deformation processes between the (020) face and the (002) and (101) faces, with slip lines observed parallel to (020). On the surface around the indents performed on (020), there is no resolvable features or cracking that could be attributed to the nanoindentation experiment, whilst the latter two faces show indents around which nucleated slip bands and cracks as well as small pile up of the material.
Overall, whilst CBZ has a strong 3D network of interactions, the system deforms plastically through slip on the (020) plane, which requires the breaking (and reforming) of hydrogen bonds. The material also shows some brittleness and surface cracking was observed on indentation. MD-Shearing simulations showed results consistent with our experimental observations of the slip plane, whilst energy of attachment failed at predicting such a slip plane. Rugosity data combined with a new parameter we defined as the HBdensity of the slip plane of interest, were able to reveal the (020) as the most likely slip plane in agreement with the experiments and the MD simulations. Computed Young's moduli gave a good picture of the anisotropy of mechanical behaviour in CBZ-III. Of particular interest, MD-shearing simulation add to a set of computable parameters such as attachment energy, planes rugosity, whose efficiency at predicting weak planes has often been challenged.

Conflicts of interest
The authors declare no conflicts of interest.