Computational Techniques for Predicting ... - ACS Publications

Mar 5, 2019 - ... Laboratory, Department of Pharmaceutics, College of Pharmacy, University of .... using various force fields (Compass II, Compass, Dr...
0 downloads 0 Views 1MB Size
Subscriber access provided by UNIV OF TEXAS DALLAS

Article

Computational techniques for predicting mechanical properties of organic crystals - A systematic evaluation Chenguang Wang, and Changquan Calvin Sun Mol. Pharmaceutics, Just Accepted Manuscript • DOI: 10.1021/acs.molpharmaceut.9b00082 • Publication Date (Web): 05 Mar 2019 Downloaded from http://pubs.acs.org on March 13, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Molecular Pharmaceutics

1

Computational techniques for predicting mechanical properties of organic

2

crystals - A systematic evaluation

3 4 5

Chenguang Wang and Changquan Calvin Sun *

6 7

Pharmaceutical Materials Science and Engineering Laboratory, Department of Pharmaceutics, College of

8

Pharmacy, University of Minnesota, Minneapolis, MN 55455, USA

9 10 11 12 13 14 15 16 17

*Corresponding author

18

Changquan Calvin Sun, Ph.D.

19

9-127B Weaver-Densford Hall

20

308 Harvard Street S.E.

21

Minneapolis, MN 55455

22

Email: [email protected]

23

Tel: 612-624-3722

24

Fax: 612-626-2125

1 ACS Paragon Plus Environment

Molecular Pharmaceutics 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

25

Page 2 of 25

ABSTRACT

26

Understanding of the structure – mechanical properties relationship in organic crystals can

27

potentially facilitate the design of crystals with desired mechanical properties through crystal

28

engineering. To understand and predict crystal mechanical properties, including tableting behavior, a

29

number of computational methods have been developed to analyze crystal structure. These include

30

visualization, attachment energy calculations, topological analysis, energy framework, and elasticity

31

tensor calculation. However, different methods often lead to conflicting predictions. There is a need

32

for a computational tool kit for predicting crystal mechanical properties from crystal structures. Using

33

α-oxalic acid anhydrous (OAA) and dihydrate (OAD) as a model system, we have systematically

34

compared the predictive accuracy of the experimentally determined mechanical properties using

35

powder compaction and nanoindentation of several methods. We have found that crystal plasticity

36

can be accurately predicted based on energy framework combined with topological analysis and DFT

37

calculated elasticity tensor. Although very useful in characterizing crystal packing features, structure

38

visualization, topology analysis, and attachment energy calculations alone are insufficient for

39

accurately identifying the slip planes and predicting mechanical properties and tableting behavior of

40

organic crystals.

41

KEYWORDS: computational prediction, organic crystal, mechanical property, structural analysis,

42

nanoindentataion, intermolecular interaction

2 ACS Paragon Plus Environment

Page 3 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

43

Molecular Pharmaceutics

Introduction

44

A clear understanding of the relationship between structure and properties can guide the design

45

of new materials or solve material related problems in various disciplines.1 The understanding of the

46

relationship between crystal structure and mechanical properties for organic crystals has been

47

significantly improved in recent years due to the advances in crystal engineering, materials

48

characterization and computational techniques. In contrast to the traditional view that organic crystals

49

are

50

plastically/elastically deforming, brittle fracturing, and twisting, have been observed and explained

51

from molecular packing features in crystal.2-4 Experimental methods for quantitative measurements

52

of mechanical properties of crystals, including nano-, micro-, macro- indentation,5-7 crystal thermal

53

expansion, high pressure crystallography, Brillouin light scattering,8 and powder compaction,9 have

54

been applied. Computational tools for characterizing crystal properties, such as attachment energy,10

55

energy framework,11 energy-vector model,12,

56

potential,15 and elastic constants calculation,16 have also been made available.

57

brittle,

diverse

responses

to

an

external

13

mechanical

stress,

including

shearing,

topological analysis,14 molecular electrostatic

For pharmaceutical tablet manufacturing, adequate crystal plasticity is a prerequisite for forming

58

tablets with sufficient mechanical strength by powder compression.

A readily accessible

59

computational tool kit that can be used to accurately predict tabletability of an active pharmaceutical

60

ingredient from its molecular structure would be extremely useful to the development of high quality

61

tablet products. Attaining this goal requires reliable predictions of crystal structure from molecular

62

structure, crystal mechanical properties from crystal structure, and tabletability from crystal

63

mechanical properties and particle properties.

64

One of the major challenges for linking crystal structure with mechanical properties is the

65

accurate prediction of intermolecular interaction strength in a given crystal structure.17 The energy

66

framework is useful for quantifying and visualizing intermolecular interaction energies, which assist

67

in slip plane identification.18, 19 Crystal plasticity was accurately ranked for three pairs of polymorph,

68

when the energy framework was combined with topological analysis of the slip plane.18 However, 3 ACS Paragon Plus Environment

Molecular Pharmaceutics 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 25

69

further development of this tool is needed as some organic crystals with similar three-dimensional

70

packing features, in term of hydrogen bonds network and energy framework, still exhibit very

71

different mechanical properties.20,

72

system, Forms II and III showed similar tape-like structures and energy frameworks, but Form II is

73

elastic and exhibits a notably higher elastic modulus and hardness than the plastic Form III.21

74

Additionally, isostructural crystals of hexachlorobenzene (HCB) and hexabromobenzene (HBB)

75

showed similar columnar packing and energy framework, but only HCB is bendable at room

76

temperature, while HBB is bendable only at high temperature.19 Given the outstanding problems in

77

predicting mechanical properties from crystal structures of single component polymorphs and

78

isostructural crystals of structurally similar molecules, the structure – mechanical properties

79

relationship among multi-component crystals (hydrate, salt and cocrystal) is expected to be even more

80

challenging.

21

In the trimorphic 3-((4-chlorophenyl)imino)indolin-2-one

81

The significant challenges in computationally predicting mechanical properties of multi-

82

component crystals are likely responsible for the few studies that address the effects of crystal

83

hydration on mechanical properties of organic crystals. Studies of p-hydroxybenzoic acid (HBA),22

84

sodium naproxen (Na-NAP),23 sodium saccharin (Na-SAC),24 galactose derivatives (GD),25

85

theophylline (TH),26 and uric acid (UA) systems revealed that incorporation of water molecules into

86

the crystal lattice increases crystal plasticity, i.e., water molecules always plasticize these organic

87

crystals.27 The incorporation of water molecules into the crystal structure was thought to either cause

88

loosening of the crystal packing or the water functions as a ‘lubricant’ in the cases of HBA

89

monohydrate, Na-NAP hydrates, and GD dihydrate.22, 23, 25 For TH, the superior plasticity of the

90

monohydrate is rationalized by forming a ladder-like packing feature, which enables the propagation

91

of dislocations and facile movement of TH molecules under stress.26 The anhydrate and dihydrate

92

crystals of UA have very similar crystal packing patterns, but the anhydrate is much stiffer (E = 17.73

93

± 2.32 GPa on the (100) facet) than the dihydrate (E = 3.49 ± 0.54 GPa on the (001) facet) despite the

94

fact that the dihydrate has slightly higher packing efficiency.27 In addition, the UA dihydrate also

95

exhibited significantly higher viscoelasticity compared to anhydrous UA. For Na-SAC, the 15/8 4 ACS Paragon Plus Environment

Page 5 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Molecular Pharmaceutics

96

hydrate exhibited anisotropic mechanical properties, where the elastic modulus on (011) facet is 14%

97

higher than that on the (101) facet. Upon dehydration, the E values increased by 15% and 3% on the

98

(0 1 1) and (1 0 1) facets, respectively, suggesting a facet dependent plasticizing effect by water.24

99

When different computational tools, e.g., structural visualization, topological analysis,

100

attachment energy calculation, energy framework, and elasticity tensor calculation, were employed

101

in previous studies,2, 5, 8, 10-12, 14-16, 20-22, 26-32 conflicting predictions on the rank order of plasticity and

102

crystal anisotropy as well as slip plane assignment were often encountered.18 To develop a reliable

103

tool kit for predicting crystal mechanical properties from crystal structures, it is imperative to compare

104

the performance of all available computational methods in one study, which has not been done before.

105

We have applied currently available computational techniques to predict the mechanical properties

106

of α-oxalic acid anhydrate (OAA) and dihydrate (OAD). The accuracy of the predictions was

107

assessed by a comparison to experimentally determined mechanical properties of these two crystals.

108 109

Materials and Methods

110

Materials

111

Both α-oxalic acid anhydrous (purity > 98%) and dihydrate (purity > 99%) forms were

112

purchased from Sigma–Aldrich (St. Louis, MO).

113

Powder X-ray diffraction pattern

114

The α-oxalic acid anhydrous and dihydrate forms were characterized using a powder X-ray

115

diffractometer (PANalytical X’pert pro, Westborough, MA) with Cu Kα radiation. PXRD patterns

116

were obtained by scanning samples from 5° to 35° with a step size of 0.017° and 2 s dwell time at

117

room temperature and 3% RH for OAA and 35% RH for OAD. The tube voltage and amperage were

118

set at 45 kV and 40 mA, respectively. X-ray was generated with a copper source with wavelength of

5 ACS Paragon Plus Environment

Molecular Pharmaceutics 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 25

119

1.5418 Å. Mercury (V. 3.10, CCDC, Cambridge, UK) was used to calculate PXRD patterns with

120

step size of 0.02°.

121

Structure Visualization and Topology Analysis

122

OAA and OAD crystal structures were visually examined to identify crystallographic planes

123

using Mercury CSD (V. 3.10, CCDC, Cambridge, UK). To minimize the thermal effect on structures,

124

crystal structures determined at 298 K (CSD refcodes: OXALAC06; OXACDH11) were used for

125

comparison. This effort was aided by showing hydrogen bonds that meet the default criteria given in

126

the software; that is, when the distance between donor atom and acceptor atom is shorter than the

127

sum of their corresponding van der Waals radii.

128

The crystal structures of OAA and OAD were simplified using the ToposPro (V 5.3.0.4,

129

Samara, Russia) to show packing feature, where hydrogen bonds were indicated with a line and each

130

molecule was represented as a sphere.28

131

The quantitative layer topology analysis in crystals was obtained using CSD Python

132

program.14 This geometric analysis approach identified the most likely slip plane based on the lowest

133

degree of interpenetration and highest distance between the separated molecular layers. If interlayer

134

distance was positive and multiple orthogonal slip planes were identified, the crystal was predicted as

135

“bendy” or “flexible”.14

136

Energy Framework

137

The dispersion-corrected density functional theory (DFT-D) model used for estimating

138

intermolecular interactions was shown to be accurate.33 The pairwise intermolecular interaction

139

energy was estimated using CrystalExplorer and Gaussian09 with experimental crystal geometry.34,

140

35

141

calculation. For each molecule in the asymmetric unit of a crystal, the total intermolecular interaction

142

energy with another molecule, calculated using the B3LYP-D2/6-31G(d,p) electron densities model,

The hydrogen positions normalized to standard neutron diffraction values were used during the

6 ACS Paragon Plus Environment

Page 7 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Molecular Pharmaceutics

143

is the sum of electrostatic, polarization, dispersion, and exchange-repulsion components with scaling

144

factors of 1.057, 0.740, 0.871, and 0.618, respectively.36 The intermolecular interaction is neglected

145

with molecule - molecule distance more than 3.8 Å.11 The OAD asymmetric unit contains a half

146

oxalic acid and one water. The interaction energies of oxalic acid and water with neighboring

147

molecules were separately calculated. The interlayer or intralayer interaction energies were calculated

148

by adding the interaction energies between a given molecule in one layer and all interacting molecules

149

in a neighboring layer or within the same layer, respectively. The interaction energies below a certain

150

energy threshold (10 kJ/mol) were omitted for clarity, and the cylinder thickness was taken to be

151

proportional to the intermolecular interaction energies in the energy framework.

152

Attachment Energy Calculation

153

The attachment energy was calculated using the crystal growth module in Materials Studio

154

2016 (BIOVIA Inc., San Diego, CA) using all available force fields (Compass II, Compass, Dreiding,

155

cvff, and pcff) and charge assigments (Forcefield assigned, Qeq, and Gasteiger) at fine quality. The

156

“Ewald” electrostatic summation method and “atom based” van der Waals summation were chosen.

157

In addition, a minimum dhkl was set at 1.0 Å.37

158

Elastic Constants Calculation

159

The elastic constants were calculated with the forcite module in Materials Studio 2016 using

160

various force fields (Compass II, Compass, Dreiding, cvff, and pcff) and charges (Forcefield assigned,

161

Qeq, and Gasteiger) at ultra fine quality. The “Ewald” electrostatic summation method and “atom

162

based” van der Waals summation were chosen during all calculations. The elastic stiffness (Cij) and

163

compliance (Sij) constants were also calculated using dispersion corrected - density functional theory

164

(DFT-D), employing generalized gradient approximation (GGA) exchange‐correlation functional of

165

Perdew, Burke, and Ernzerhof (PBE) and corrected by Tkatchenko–Scheffler dispersive pairwise

166

schemes.38-40 The convergence thresholds during atomic coordinate geometry optimization were set

167

at 10-5 Hartree (Energy), 0.002 Hartree Å-1 (Max. force) and 0.005 Å (Max. displacement), while the 7 ACS Paragon Plus Environment

Molecular Pharmaceutics 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 25

168

unit cell dimensions were fixed. A 6×6 symmetric elastic constants matrix was calculated to describe

169

fully the stress-strain relationship for both crystals. Due to the symmetry of the crystal structure,

170

independent elastic constants were reduced to 9 values for OAA (orthorhombic) and 13 values for

171

OAD (monoclinic). The bulk moduli, shear moduli, Young’s moduli, and Poisson’s ratio were

172

calculated from the matrix based on the elasticity theory. The crystal anisotropy index was calculated

173

as the ratio of the largest to the smallest Young’s modulus.41

174

Dynamic Water Vapor Sorption Isotherm

175

Water sorption isotherms of the OAA and OAD powder (~15 mg) were obtained using an

176

automated vapor sorption analyzer (Intrinsic DVS, Surface Measurement Systems Ltd., Allentown,

177

PA) at 25 °C. The nitrogen flow rate was 50 mL/min. As appropriate, each sample was first purged

178

with either dry or humidified nitrogen until a constant weight was obtained. The sample was then

179

exposed to a series of relative humidities (RH) between 0% and 95% with the step size of 5% RH.

180

At each specific RH, the sample was assumed to have reached equilibration when either the relative

181

rate of mass change (dm/dt) was less than 0.002% per minute with a minimum equilibration time of

182

0.5 h or a maximum equilibration time of 6 h was met. The RH was then changed to the next target

183

value.

184

Thermal Expansion and High Pressure Crystallography

185

The available crystallographic parameters of OAA (OXALAC05-06) and OAD

186

(OXACDH15-24) as a function of temperature were used to calculate crystal thermal expansivity.42,

187

43

188

Nanoindentation

189

Triboindenter I980 (Hysitron, Minneapolis, MN) equipped with a three-sided pyramidal Berkovich

190

tip was employed to determine the OAD mechanical properties at ambient conditions (23 °C, 32%

191

RH). Large OAD crystals were mounted onto a glass slide using superglue, and the slide was fixed

High pressure crystallographic parameters of OAD were extracted from OXACDH35-38.44

8 ACS Paragon Plus Environment

Page 9 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Molecular Pharmaceutics

192

on the stage by vacuum. The indented crystal face was indexed as (1 0 -1) by XRD. A smooth area

193

on the crystal surface was chosen under optical microscopy and indented after the roughness was

194

verified by in situ imaging to be less than 30 nm. Indentations were conducted under displacement

195

control with a target maximum depth of penetration of 1000 nm. Loading was performed using a 200

196

nm/s loading rate, followed by 10 s holding at peak displacement, and then unloading at the same

197

rate. The initial unloading part of the force-displacement curve was used to extract the elastic contact

198

stiffness for calculating reduced elastic modulus (Er). The Er and indentation hardness (H) were

199

calculated using appropriated equations according to the Oliver–Pharr method.45 The Poisson's ratio

200

of 0.3 was used for calculating the Young’s modulus (E) from Er.16, 46, 47 Before each set of runs, the

201

tip area function was calibrated by indenting fused silica with a series of indents at various contact

202

depths.

203

Powder compaction

204

Sieved size fractions, 125–250, 250-355, and 355-500 μm, (USA standard sieves, W.S. Tyler

205

Industrial Group, Mentor, OH) of OAA and OAD were used for tableting. A compaction simulator

206

(Presster, Metropolitan Computing Corp., NJ), simulating a 10 stations Korsch XL100 press with a

207

10.0 mm diameter, flat-faced round tooling, was used to evaluate powder compaction properties.

208

Tablets were compressed without lubricant over a compaction pressure range of 25–250 MPa at a

209

speed of 100 ms dwell time. Tablets were broken diametrically on a texture analyzer (Texture

210

Technologies Corp., Surrey, UK) immediately after tablet ejection. Tablet tensile strength was

211

calculated from breaking force and tablet dimensions, following a standard procedure.48

212 213

Results and Discussion

214

Solid-state stability

215

In order to accurately characterize mechanical properties and tableting behavior, it is

216

necessary to ensure the solid-state stability of OAA and OAD. Considering that the hydrate form is 9 ACS Paragon Plus Environment

Molecular Pharmaceutics 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 25

217

stable at high RHs, the OAD dehydration experiment was started at 95% RH (Figure 1a). Dehydration

218

from OAD to OAA was observed at a RH of 5%. The 39.6% weight loss at 0% RH matched well

219

with its theoretical OAD water content (39.99 %) indicating complete dehydration. Hysteresis during

220

the dehydration - hydration process was observed, where the onset of OAA hydration occurred at

221

15% RH and continued through 70% RH. Therefore at 25 °C, OAA is kinetically stable at RHs ≤

222

5%, and OAD is kinetically stable at RHs ≥15% RH. Phase purity of OAA and OAD powders,

223

equilibrated at 3% RH and 35% RH for 6 h, were confirmed by PXRD. All peaks predicted from the

224

crystal structures were observed in the corresponding experimental PXRD patterns, and no extraneous

225

peaks were present (Figure 1b). Thus, PXRD data suggested that both OAA and OAD were phase

226

pure. The different peak ratios and peak shapes between experimental and predicted PXRD patterns

227

are attributed to preferred orientation of crystals during data collection. To maintain the phase purity,

228

OAA and OAD powders were equilibrated at 0% RH (over phosphorus pentoxide) and 32% (MgCl2

229

saturated aqueous solution), respectively, for more than two weeks before the compaction study. The

230

compaction behavior, evaluated at ambient temperature and 3% RH for OAA and 35% RH for OAD

231

powders, was used to evaluate the accuracy of the predicted crystal mechanical properties and

232

tabletability using different computation methods. a)

233 234 235

b)

Figure 1. a) Effect of relative humidity (RH) on the kinetic stability of OAD and OAA; b) verification of the phase purity of OAA and OAD powders by PXRD.

236 10 ACS Paragon Plus Environment

Page 11 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

237

Molecular Pharmaceutics

Structure analysis and mechanical properties prediction

238

The key crystallographic parameters of OAA and OAD are summarized in Table S1. The

239

orthorhombic α-OAA crystal consists of O–H···O (2.708 Å) hydrogen bonded layers parallel to the

240

bc plane and stacking along the a axis (Figure 2a). These layers interact through weak C–H···O

241

hydrogen bonds (2.884, 2.895 Å). Thus, (1 0 0) is identified as slip plane for OAA by structure

242

visualization. High plasticity of OAA is expected from such a structure. For monoclinic crystal

243

system OAD, which showed a dense three-dimensional (3D) hydrogen bonded network, limited

244

crystal plasticity is expected. Within the 3D network, each oxalic acid is connected with six water

245

molecules through different O–H···O (2.513 Å, 2.883 Å, 2.995 Å) hydrogen bonds. Each water

246

molecule interacts with three oxalic acid molecules through four O–H···O hydrogen bonds (2.513 Å,

247

2.863 Å, 2.883 Å, 2.995 Å) to form an isolated hydrate (Figure 2b). The qualitative structure analysis

248

predicts that the layered OAA is more plastic than OAD, which has a three dimensional hydrogen

249

bonded network. An inspection of the ToposPro crystal packing pattern analysis suggested a (2 0 0)

250

slip plane in layered OAA, but no clear slip planes could be identified for OAD due to the extensive

251

three dimensional hydrogen bonded OAD crystal (Figure 3). These results suggest that OAA is more

252

plastic than OAD. a)

b)

253 254 255

Figure 2. Crystal packing patterns of a) OAA and b) OAD. A likely slip layer of OAA by visualization is shaded in pink. No slip layer can be readily identified in OAD. 11 ACS Paragon Plus Environment

Molecular Pharmaceutics 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 25

256

A quantitative layer topological analysis predicts the most likely slip planes based on the ease of

257

molecular layers sliding, judged from a set of calculated geometric descriptors, such as layer separation,

258

presence of orthogonal planes, presence of interlayer hydrogen bonds, and dimensionality of hydrogen

259

bonding network in the structure.14 The relative plasticity between two crystals is predicted based on the

260

same set of descriptors, where larger separation distance, lack of interlayer hydrogen bonds, and

261

presence of orthogonal slip planes favor easier slip.14 Such analysis of OAA and OAD (Table S2)

262

predicted the primary slip planes of (2 0 0) for OAA and (1 0 1) for OAD, and that OAA was more

263

plastic than OAD. This is consistent with the predicted plasticity order between OAD and OAA by

264

qualitative ToposPro crystal packing pattern analysis (Figure 3). a)

b) (200)

(101)

265 266 267 268

Figure 3. Simplified crystal packing structure by ToposPro analysis of a) OAA and b) OAD. OA and water molecules are represented by black and yellow spheres, respectively. Hydrogen bonds are represented by lines. The likely slip layer in each crystal is shaded in purple.

269

It is common practice to identify possible slip planes as the Miller planes with the lowest

270

attachment energy in a crystal.10 The three planes in OAA with the lowest attachment energy are (2

271

0 0) (-48.18 kcal/mol), (1 1 1) (-56.13 kcal/mol), and (0 0 2) (-57.99 kcal/mol), and those in OAD are

272

(1 0 -1) (-16.34 kcal/mol), (1 0 1) (-21.39 kcal/mol), and (0 1 1) (-31.92 kcal/mol) by using the

273

classical molecular mechanics Forcite method with force fields of compass II and Qeq charge. Thus,

274

the Eatt of the primary slip plane (2 0 0) in OAA is about three-fold larger than that of the (1 0 -1) in

275

OAD. Therefore, slip between (1 0 -1) planes in OAD is energetically more favorable than that 12 ACS Paragon Plus Environment

Page 13 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Molecular Pharmaceutics

276

between (2 0 0) planes in OAA. However, slip planes identified based on attachment energy for both

277

crystals depend on the force field and charge applied for calculation (Table 1), as previously shown.10

278

In addition, the identified slip planes with the lowest Eatt may not be the easiest to slip, if they assume

279

a rough topology.10, 14 This explains why slip planes in the OAD crystal predicted by the attachment

280

energy (1 0 -1) and topological analysis (1 0 1) are different.

281

Table 1. Predicted crystal morphology and attachment energy for OAA and OAD.

Force field

COMPASSII

COMPASS

Dreiding pcff

cvff

Charge

OAA

OAD

Slip Plane

Eatt(kcal/mol)

Slip Plane

Eatt(kcal/mol)

Qeq

(2 0 0)

-48.18

(1 0 -1)

-16.34

Gasteiger

(1 1 1)

-34.47

(0 0 2)

-1.39

Forcefield assigned

(2 0 0)

-38.33

(1 0 -1)

-15.89

Qeq

(0 0 2)

-58.31

(1 0 -1)

-16.34

Gasteiger

(1 1 1)

-30.47

(0 0 2)

-1.39

Forcefield assigned

(2 0 0)

-38.33

(1 0 -1)

-15.89

Qeq

(2 0 0)

-34.63

unstable

Gasteiger

(2 0 0)

-17.31

unstable

Forcefield assigned

(2 0 0)

-68.97

(1 0 -1)

-25

Qeq

(2 0 0)

-68.54

(1 0 -1)

-25.39

Gasteiger

(1 1 1)

-47.26

Forcefield assigned

(2 0 0)

-34.47

(1 0 -1)

-17.93

Qeq

(2 0 0)

-48.26

(1 0 -1)

-19.8

Gasteiger

(2 0 0)

-30.95

(0 0 2)

-4.37

unstable

282

The energy frameworks of OAA and OAD are distinctly different (Figure 4). It is also

283

noteworthy that the major component of the intermolecular interactions in both OAA and OAD are

284

electrostatic, while non-directional dispersive components are minimal (Table S3).

285

polarization energy was observed in the water-oxalic acid interaction, which was attributed to the O–

286

H···O hydrogen bond (2.513 Å). In OAA, (1 0 0) has the highest intralayer bonding energy (-131.2

287

kJ/mol), which is more than two-fold greater than the interlayer bonding energy (-64.4 kJ/mol)

Strong

13 ACS Paragon Plus Environment

Molecular Pharmaceutics 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 25

288

(Figure 4a), suggesting that the (1 0 0) is the primary slip plane. In OAD, the intralayer bonding

289

energy for (1 0 -1) is higher (-207.4 kJ/mol) than (1 0 1) (-171.4 kJ/mol), while the interlayer bonding

290

energy between (1 0 -1) planes (-48.2 kJ/mol) is lower than that between than (1 0 1) planes (-102.4

291

kJ/mol). Such energetic differences imply that (1 0 -1) is more likely the primary slip plane for OAD.

292

Sliding between (1 0 -1) planes is also unobstructed based on the quantitative layer topological

293

analysis. Although both planes are topological favorable, the more energetically favored (1 0 -1)

294

plane is the primary slip plane of OAD, which has an approximately 25% lower interlayer bonding

295

energy than that between (1 0 0) in OAA. Therefore, OAD is expected to be more plastic.

296

The energy framework combined with topological analysis of crystal structures (Figure 4)

297

predicts a plasticity order that contradicts that deduced from the simple structure visualization (Figure

298

2) and ToposPro analysis (Figure 3). However, the former energy based method is expected to be

299

more reliable.18 To verify this, the full elasticity matrices were also calculated for OAA and OAD.

300

The matrices allow the calculation of crystal facet specific E as well as key mechanical properties,

301

including bulk modulus, Young's modulus, shear modulus, and Poisson's ratio. The calculation of

302

elasticity matrices was performed using both forcite method and DFT.

303

method, none of the combinations of forcefield and charge could generate E values for both OAA

304

and OAD (Table 2). Consequently, the order of plasticity between OAA and OAD could not be

305

predicted by the forcite method. a)

(100)

b)

(101̅)

However, in the forcite

(101)

306 307 308

Figure 4. Energy frameworks of a) OAA and b) OAD with a likely slip layer shaded in blue. The thickness of each cylinder (in blue) represents the relative strength of intermolecular interaction. The 14 ACS Paragon Plus Environment

Page 15 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Molecular Pharmaceutics

309 310

energy threshold for the energy framework is set at −10 kJ/mol. A (1 0 1) layer in OAD is indicated in red box.

311

Table 2. Predicted mechanical properties of crystals by the forcite method.

Force field

Charge

OAA Bulk modulus (GPa)

COMPASSII

COMPASS

Dreiding pcff

cvff

Qeq

OAD Anisotropy

Bulk modulus (GPa)

unstable

unstable

Gasteiger

12.517

4.573

unstable

Forcefield assigned

14.548

3.08

unstable

Qeq

Anisotropy

unstable

unstable

Gasteiger

14.82

4.765

unstable

Forcefield assigned

15.95

3.106

unstable

Qeq

22.392

20.65

unstable

Gasteiger

unstable

unstable

Qeq

unstable

Gasteiger

unstable

Forcefield assigned

unstable

Qeq

unstable

unstable

Gasteiger

unstable

unstable

Forcefield assigned

unstable

18.293

2.683 unstable

16.973

13.575

4.222

21.54

In contrast to the poor performance of forcite method, the DFT method successfully predicted elastic constants for both OAA and OAD. For OAA, the relatively high and comparable principal elasticity constants (C11 = 29.01, C22 = 31.34, and C33 = 38.45 GPa) suggest OAA is relatively isotropic in term of E. The Voigt (upper bound) and Reuss (lower bound) averaging schemes of the E (25.74 to 35.01 GPa) and shear modulus (10.06-14.48 GPa) of OAA are much higher than those of OAD (E: 11.86-21.29 GPa and shear modulus:4.37-8.20 GPa). Figure 5 displays 3D Emaps of OAA and OAD. The highest E for OAA (54.36 GPa) and OAD (44.67 GPa) both correspond to strong intralayer interactions within (1 0 0) and (1 0 -1). In OAA, the highest E corresponds to the compression direction along the O–H···O hydrogen bonded chains and the 15 ACS Paragon Plus Environment

Molecular Pharmaceutics 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 25

lowest E is along the direction that only weak C–H···O hydrogen bonds could be identified. The higher bulk modulus of OAA (19.49 GPa) compared to OAD (13.89 GPa) suggests that OAA is less compressible than OAD. The higher plasticity of OAD is also supported by its larger anisotropy index (7.86) relative to OAA (3.06).49

a)







b)



312

OAA

OAD

313

Figure 5. Three-dimensional distribution of Young's modulus for a) OAA and b) OAD.

314

Verification of Mechanical Properties Prediction

315

The predicted order of crystal plasticity and slip plane by various computational methods are not

316

in agreement. To determine unambiguously the accuracy of predictions by these methods, we

317

experimentally assessed plasticity of OAA and OAD and tabletability by powder compaction. For

318

the same particle size fraction (125-250 m), the tabletability of OAD was significantly better than

319

OAA (Figure 6). To assess a possible size effect,50 tabletability of two larger particle size fractions

320

(250-355 m and 355-500 m) of OAD were determined, but they only exhibited slightly lower tablet

321

tensile strength in the high pressure (150-250 MPa) range (Figure 6). Thus, the significantly better

322

tabletability of OAD over OAA was not due to differences in particle size distribution. Rather, it 16 ACS Paragon Plus Environment

Page 17 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Molecular Pharmaceutics

323

reflected the different mechanical properties between the two crystals. The higher plasticity of OAD

324

is supported by its better compressibility, where lower tablet porosity was obtained for OAD than

325

OAA at the same pressure. The lower tablet porosity is expected to correspond to a larger bonding

326

area between OAD crystals when compressed. When normalized by porosity, the tablet tensile

327

strength of OAD is also higher than that of OAA (Figure 6c). The higher compactibility of OAD

328

suggests it has higher bonding strength. Therefore, the significantly better tabletability of OAD is a

329

result of both larger bonding area and higher bonding strength.51 a)

b)

c)

330 331 332

Figure 6. Mechanical properties of OAA and OAD assessed by powder compaction. a) tabletability, (b) compressibility, and c) compactibility.

333

Having experimentally verified the higher plasticity of OAD than OAA, it was of interest to also

334

experimentally verify the predicted crystal anisotropy and slip planes. This was achieved using 17 ACS Paragon Plus Environment

Molecular Pharmaceutics 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 25

335

experimental anisotropic thermal expansion behavior, where larger thermal expansion is expected

336

along the direction of weaker intermolecular interactions.32 Unit cell parameters of OAD determined

337

at 100 and 300K and those of OAA determined at 130 and 298 K were used for this purpose.42, 43 At

338

room temperature, OAA density (1.905 gcm−3) and packing efficiency (79.2%) were both higher than

339

those of OAD (density: 1.639 gcm−3; packing efficiency: 73.4%). The larger volume thermal

340

gradient of OAD (0.047 Å3/K) compared to OAA (0.041 Å3/K) corresponded to the higher thermal

341

expansivity of OAD. For OAA, the axial thermal expansivity followed the order: a (0.00039 Å/K) >

342

c (0.00029 Å/K) > b (0.00020 Å/K). The largest expansivity along the a axis is consistent with the

343

absence of hydrogen bonds and the largest interlayer distance between (1 0 0) planes. Therefore, the

344

thermal expansivity data suggested (1 0 0) as the primary slip plane of OAA. For OAD, the axial

345

thermal expansivity follows the order: b (0.00061 Å/K) > c (0.00055 Å/K) > a (0.00011 Å/K),

346

although the β angle also slightly increased upon heating or decompression. This order is the same

347

as the observed order of axial pressure sensitivity b (-0.10 Å/GPa) > c (-0.083 Å/GPa) > a (-0.075

348

Å/GPa). Both confirmed the interactions along b were the weakest in OAD.44 However, the

349

molecular layers slip along the (0 1 0) is topologically unfavorable. Hence, the thermal and pressure

350

expansivity data could not correctly identify the slip plane for OAD. In summary, OAD is more

351

anisotropic than OAA based on the thermal expansivity, and the primary slip plane for OAA is (100). a)

b) Holding

352 353 354 355

Figure 7. Nanoindentation data (displacement control) obtained on (1 0 -1) facet of OAD. a) A representative loading force–displacement curve (holding stage is indicated by a dashed line), and b) three stages of loading force and tip penetration depth as a function of time. 18 ACS Paragon Plus Environment

Page 19 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Molecular Pharmaceutics

356

To verify the accuracy of the predicted elastic moduli by the DFT-D method, the mechanical

357

properties of OAD were also determined using nanoindentation. When the (1 0 -1) facet of OAD was

358

indented, large ‘pop-in’ events were absent during loading (Figure 7). This is consistent with a plastic

359

crystal having multiple slip systems.21 During holding at the maximum tip penetration, the force

360

continuously decreased by approximately 30% in 10 s. The significant relaxation behavior also

361

indicates high viscoplasticity of OAD. The closeness between measured E (18.810.96 GPa) and

362

calculated E (19-20 GPa) of OAD (1 0 -1) facet confirmed the reliability of the DFT-D method.

363

The E of OAD was much higher than the mean (12.97 GPa) and median (11.15 GPa) of the

364

reported values for organic crystals by nanoindentation (n=220), indicting high stiffness of OAD.

365

This is consistent with the presence of the three-dimensional hydrogen bond network in OAD. On

366

the other hand, the measured H, 0.350.03 GPa, fell in the range of the previously reported bendable

367

crystals of diphenhydramine-acesulfame salt (0.300.04 GPa) and -succinic acid (0.460.01 GPa),

368

which are plastic due to the presence of multiple slip planes in their structures.9, 29, 52 Thus, high

369

plasticity of OAD is also expected. The H value is also smaller than both the mean (0.55 GPa) and

370

median (0.43 GPa) of the reported H values of organic crystals by nanoindentation.29, 53 As a result

371

of the low H and high E values, the E/H value of OAD (53.7) is much higher than those of known

372

brittle pharmaceutical crystals (ranging 15 – 25), such as metformin HCl (22), metformin-acsulfmate

373

(18), piroxicam (25), famotidine Form A (15) and Form B (23), α- nifedipine (19), and olanzapine

374

(19).9

375

The visualization of hydrogen bonds network is the most commonly used approach to

376

characterize the crystal structures and predict mechanical properties. However, the hydrogen bonds

377

and other intermolecular interactions are geometrically sensitive. Thus, variations in interaction

378

strengths cannot be easily identified by qualitative visualization of structures.17, 30, 54 Consequently,

379

misleading predictions can be made. The topology analysis describes the major packing feature and

380

infers the ease of molecular layers movement. This is useful information, but prediction of mechanical

381

properties may be erroneous, because the interaction energy is neglected,14 which also plays an 19 ACS Paragon Plus Environment

Molecular Pharmaceutics 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 25

382

important role in the slip between molecular layers. The attachment energy calculation only accounts

383

for the interlayer interaction energy without considering the important factors of layer roughness and

384

layer interlocking on the ease of slip between layers.18 Hence, predicted slip planes based on

385

attachement energy alone are often incorrect. Accurate prediction of mechanical properties requires

386

systematic use of complementary computational methods.

387

collected using material-sparing techniques, such as nanoindentation, should be used to validate

388

predictions.

Where possible, experimental data

389 390

Conclusion

391

We have systematically examined the performance of a number computational techniques,

392

i.e., the visualization method, qualitative and quantitative topology analysis, attachment energy,

393

energy framework, and elastic constant calculation, for assembling a useful computational toolbox

394

for predicting mechanical properties from crystal structures. Using OAA and OAD as a model

395

system, we show that the crystal structure visualization, topology analysis, and attachment energy

396

calculation alone failed to accurately predict relative order of mechanical properties and tableting

397

behavior of OAA and OAD. However, the quantitative analyses of intermolecular interaction

398

strength by DFT-D model, graphically represented in the form of energy framework and 3D elasticity

399

map, successfully predicted the different mechanical properties of the two crystals, which were also

400

experimentally confirmed. Although the incorporation of water generates an extensive 3D H-bond

401

network, the OAD crystal exhibits higher interaction energy anisotropy and higher plasticity, which

402

explains the significantly better tabletability than the OAA crystal. Thus, the computational toolbox

403

should consist of principally the energy framework, 3D elasticity map, and layer topology analysis,

404

which are complemented by other computational methods, including attachement energy,

405

visualization, and qualitative ToposPro analysis, in order to make correct predictions of the

406

mechanical properties. The tool-kit described in this wrok is applicable to predicting mechanical

407

properties of other organic crystal systems. 20 ACS Paragon Plus Environment

Page 21 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Molecular Pharmaceutics

408 409

Supporting Information

410

Crystallographic parameters, topology analysis, and intermolecular interaction energies of OAA and

411

OAD.

412 413

Acknowledgments

414

Parts of this work were carried out in the Characterization Facility, University of Minnesota, a

415

member of the NSF-funded Materials Research Facilities Network (http://www.mrfn.org) via the

416

MRSEC program. We thank the Minnesota Supercomputing Institute (MSI) at the University of

417

Minnesota for providing the resources that contributed to the research results reported within this

418

paper. URL: http://www.msi.umn.edu. C.W. thanks Prof. N.A. Mara at Department of Chemical

419

Engineering and Materials Science, University of Minnesota for training on the triboindentor and Dr.

420

M. J. Bryant at the Cambridge Crystallographic Data Centre for his support in crystal topology

421

analysis.

422

References 1. Sun, C. C. Materials science tetrahedron—a useful tool for pharmaceutical research and development. J. Pharm. Sci. 2009, 98, 1671-1687. 2. Reddy, C. M.; Padmanabhan, K. A.; Desiraju, G. R. Structure−property correlations in bending and brittle organic crystals. Cryst. Growth Des. 2006, 6, 2720-2731. 3. Saha, S.; Desiraju, G. R. Crystal engineering of hand-twisted helical crystals. J. Am. Chem. Soc. 2017, 139, 1975-1983. 4. Saha, S.; Desiraju, G. R. Trimorphs of 4-bromophenyl 4-bromobenzoate. Elastic, brittle, plastic. Chem. Comm. 2018, 54, 6348-6351. 5. Kiran, M. S. R. N.; Varughese, S.; Reddy, C. M.; Ramamurty, U.; Desiraju, G. R. Mechanical anisotropy in crystalline saccharin: Nanoindentation studies. Cryst. Growth Des. 2010, 10, 46504655.

21 ACS Paragon Plus Environment

Molecular Pharmaceutics 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 25

6. Duncan-Hewitt, W. C.; Weatherly, G. C. Evaluating the hardness, young's modulus and fracture toughness of some pharmaceutical crystals using microindentation techniques. J. Mater. Sci. 1989, 8, 1350-1352. 7. Patel, S.; Sun, C. C. Macroindentation hardness measurement—modernization and applications. Int. J. Pharm 2016, 506, 262-267. 8. Singaraju, A. B.; Nguyen, K.; Gawedzki, P.; Herald, F.; Meyer, G.; Wentworth, D.; Swenson, D. C.; Stevens, L. L. Combining crystal structure and interaction topology for interpreting functional molecular solids: A study of theophylline cocrystals. Cryst. Growth Des. 2017, 17, 6741-6751. 9. Wang, C.; Paul, S.; Wang, K.; Hu, S.; Sun, C. C. Relationships among crystal structures, mechanical properties, and tableting performance probed using four salts of diphenhydramine. Cryst. Growth Des. 2017, 17, 6030-6040. 10. Sun, C. C.; Kiang, Y. H. On the identification of slip planes in organic crystals based on attachment energy calculation. J. Pharm. Sci. 2008, 97, 3456-3461. 11. Turner, M. J.; Thomas, S. P.; Shi, M. W.; Jayatilaka, D.; Spackman, M. A. Energy frameworks: Insights into interaction anisotropy and the mechanical properties of molecular crystals. Chem. Comm. 2015, 51, 3735-3738. 12. Upadhyay, P. P.; Sun, C. C.; Bond, A. D. Relating the tableting behavior of piroxicam polytypes to their crystal structures using energy-vector models. Int. J. Pharm 2018, 543, 46-51. 13. Zolotarev, P. N.; Moret, M.; Rizzato, S.; Proserpio, D. M. Searching new crystalline substrates for ombe: Topological and energetic aspects of cleavable organic crystals. Cryst. Growth Des. 2016, 16, 1572-1582. 14. Bryant, M. J.; Maloney, A. G. P.; Sykes, R. A. Predicting mechanical properties of crystalline materials through topological analysis. CrystEngComm 2018, 20, 2698-2704. 15. Thomas, S. P.; Shi, M. W.; Koutsantonis, G. A.; Jayatilaka, D.; Edwards, A. J.; Spackman, M. A. The elusive structural origin of plastic bending in dimethyl sulfone crystals with quasi‐isotropic crystal packing. Angew. Chem. Int. Ed. 2017, 129, 8588-8592. 16. Azuri, I.; Meirzadeh, E.; Ehre, D.; Cohen, S. R.; Rappe, A. M.; Lahav, M.; Lubomirsky, I.; Kronik, L. Unusually large young’s moduli of amino acid molecular crystals. Angew. Chem. Int. Ed. 2015, 54, 13566-13570. 17. Gavezzotti, A. The "sceptical chymist": Intermolecular doubts and paradoxes. CrystEngComm 2013, 15, 4027-4035. 18. Wang, C.; Sun, C. C. Identifying slip planes in organic polymorphs by combined energy framework calculations and topology analysis. Cryst. Growth Des. 2018, 18, 1909-1916. 19. Edwards, A. J.; Mackenzie, C. F.; Spackman, P. R.; Jayatilaka, D.; Spackman, M. A. Intermolecular interactions in molecular crystals: What's in a name? Faraday Discuss. 2017, 203, 93112. 20. Ghosh, S.; Mondal, A.; Kiran, M. S. R. N.; Ramamurty, U.; Reddy, C. M. The role of weak interactions in the phase transition and distinct mechanical behavior of two structurally similar caffeine co-crystal polymorphs studied by nanoindentation. Cryst. Growth Des. 2013, 13, 4435-4441. 21. Raju, K. B.; Ranjan, S.; Vishnu, V. S.; Bhattacharya, M.; Bhattacharya, B.; Mukhopadhyay, A. K.; Reddy, C. M. Rationalizing distinct mechanical properties of three polymorphs of a drug adduct by nanoindentation and energy frameworks analysis: Role of slip layer topology and weak interactions. Cryst. Growth Des. 2018. 22. Sun, C. C.; Grant, D. J. Improved tableting properties of p-hydroxybenzoic acid by water of crystallization: A molecular insight. Pharm. Res. 2004, 21, 382-386.

22 ACS Paragon Plus Environment

Page 23 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Molecular Pharmaceutics

23. Joiris, E.; Di Martino, P.; Malaj, L.; Censi, R.; Barthélémy, C.; Odou, P. Influence of crystal hydration on the mechanical properties of sodium naproxen. Eur. J. Pharm. Biopharm. 2008, 70, 345356. 24. Kiran, M.; Varughese, S.; Ramamurty, U.; Desiraju, G. R. Effect of dehydration on the mechanical properties of sodium saccharin dihydrate probed with nanoindentation. CrystEngComm 2012, 14, 2489-2493. 25. Panda, M. K.; Bhaskar Pal, K.; Raj, G.; Jana, R.; Moriwaki, T.; Mukherjee, G. D.; Mukhopadhyay, B.; Naumov, P. e. Flexibility in a molecular crystal accomplished by structural modulation of carbohydrate epimers. Cryst. Growth Des. 2017, 17, 1759-1765. 26. Chang, S.-Y.; Sun, C. C. Superior plasticity and tabletability of theophylline monohydrate. Mol. Pharm. 2017, 14, 2047-2055. 27. Liu, F.; Hooks, D. E.; Li, N.; Mara, N. A.; Swift, J. A. Mechanical properties of anhydrous and hydrated uric acid crystals. Chem. Mater. 2018, 30, 3798-3805. 28. Blatov, V. A.; Shevchenko, A. P.; Proserpio, D. M. Applied topological analysis of crystal structures with the program package topospro. Cryst. Growth Des. 2014, 14, 3576-3586. 29. Jing, Y.; Zhang, Y.; Blendell, J.; Koslowski, M.; Carvajal, M. T. Nanoindentation method to study slip planes in molecular crystals in a systematic manner. Cryst. Growth Des. 2011, 11, 52605267. 30. Khomane, K. S.; Bansal, A. K. Weak hydrogen bonding interactions influence slip system activity and compaction behavior of pharmaceutical powders. J. Pharm. Sci. 2013, 102, 4242-4245. 31. Krishna, G. R.; Devarapalli, R.; Lal, G.; Reddy, C. M. Mechanically flexible organic crystals achieved by introducing weak interactions in structure: Supramolecular shape synthons. J. Am. Chem. Soc. 2016, 138, 13561-13567. 32. Rather, S. A.; Saha, B. K. Thermal expansion study as a tool to understand the bending mechanism in a crystal. Cryst. Growth Des. 2018, 18, 2712-2716. 33. Thomas, S. P.; Spackman, P. R.; Jayatilaka, D.; Spackman, M. A. Accurate lattice energies for molecular crystals from experimental crystal structures. J. Chem. Theory Comput. 2018, 14, 16141623. 34. M. J. Turner, J. J. M., S. K. Wolff, D. J. Grimwood, P. R. Spackman, D. Jayatilaka and M. A. Spackman, Crystalexplorer17. University of western australia. Http://hirshfeldsurface.Net. 2017. 35. Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J., Gaussian 09. Wallingford CT, 2009. 36. Turner, M. J.; Grabowsky, S.; Jayatilaka, D.; Spackman, M. A. Accurate and efficient model energies for exploring intermolecular interactions in molecular crystals. J Phys Chem Lett. 2014, 5, 4249-4255. 37. Hartman, P.; Bennema, P. The attachment energy as a habit controlling factor: I. Theoretical considerations. J. Cryst. Growth 1980, 49, 145-156. 23 ACS Paragon Plus Environment

Molecular Pharmaceutics 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 25

38. Delley, B. An all ‐ electron numerical method for solving the local density functional for polyatomic molecules. J. Chem. Phys. 1990, 92, 508-517. 39. Delley, B. From molecules to solids with the dmol 3 approach. J. Chem. Phys. 2000, 113, 7756-7764. 40. McNellis, E. R.; Meyer, J.; Reuter, K. Azobenzene at coinage metal surfaces: Role of dispersive van der waals interactions. Phys. Rev. B 2009, 80, 205414. 41. Ranganathan, S. I.; Ostoja-Starzewski, M. Universal elastic anisotropy index. Phys. Rev. Lett. 2008, 101, 055504. 42. Wang, Y.; Tsai, C.; Liu, W.; Calvert, L. Temperature‐dependence studies of α‐oxalic acid dihydrate. Acta Crystallogr. B 1985, 41, 131-135. 43. Thalladi, V. R.; Nüsse, M.; Boese, R. The melting point alternation in α,ω-alkanedicarboxylic acids. J. Am. Chem. Soc. 2000, 122, 9227-9236. 44. Casati, N.; Macchi, P.; Sironi, A. Hydrogen migration in oxalic acid di-hydrate at high pressure? Chem. Comm. 2009, 2679-2681. 45. Oliver, W. C.; Pharr, G. M. Measurement of hardness and elastic modulus by instrumented indentation: Advances in understanding and refinements to methodology. J. Mater. Res. 2004, 19, 320. 46. Mazel, V.; Busignies, V.; Diarra, H.; Tchoreloff, P. On the links between elastic constants and effective elastic behavior of pharmaceutical compacts: Importance of poissons ratio and use of bulk modulus. J. Pharm. Sci. 2013, 102, 4009-4014. 47. Wildfong, P. L.; Hancock, B. C.; Moore, M. D.; Morris, K. R. Towards an understanding of the structurally based potential for mechanically activated disordering of small molecule organic crystals. J. Pharm. Sci. 2006, 95, 2645-2656. 48. Fell, J. T.; Newton, J. M. Determination of tablet strength by the diametral-compression test. J. Pharm. Sci. 1970, 59, 688-691. 49. Bahl, D.; Singaraju, A. B.; Stevens, L. L. Aggregate elasticity and tabletability of molecular solids: A validation and application of powder brillouin light scattering. AAPS PharmSciTech 2018, 19, 3430-3439. 50. Sun, C.; Grant, D. J. W. Effects of initial particle size on the tableting properties of l-lysine monohydrochloride dihydrate powder. Int. J. Pharm 2001, 215, 221-228. 51. Osei-Yeboah, F.; Chang, S.-Y.; Sun, C. C. A critical examination of the phenomenon of bonding area - bonding strength interplay in powder tableting. Pharm. Res. 2016, 33, 1126-1132. 52. Mishra, M. K.; Ramamurty, U.; Desiraju, G. R. Hardness alternation in α,ωalkanedicarboxylic acids. Chem Asian J 2015, 10, 2176-2181. 53. Wang, C.; Hu, S.; Sun, C. C. Expedited development of diphenhydramine orally disintegrating tablet through integrated crystal and particle engineering. Mol. Pharm. 2017, 14, 33993408. 54. Desiraju, G. R. Hydrogen bridges in crystal engineering:  Interactions without borders. Acc. Chem. Res. 2002, 35, 565-573.

24 ACS Paragon Plus Environment

Page 25 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Molecular Pharmaceutics

Information for Table of Content only Title: Computational techniques for predicting mechanical properties of organic crystals - A systematic evaluation Author: Chenguang Wang and Changquan Calvin Sun

Synopsis Performance of currently available computational tools for predicting mechanical properties from crystal structure was systematically evaluated using experimentally determined mechanical properties as a reference.

25 ACS Paragon Plus Environment