Ogre: A Python package for molecular crystal surface generation with applications to surface energy and crystal habit prediction

The surface properties of molecular crystals determine critical performance parameters for pharmaceuticals,^1–3 organic electronics,^4–6 and energetic materials.^7–9 Polymorphism, which is the capability of a molecule to crystallize into multiple forms,^10,11 expands the design space and increases the opportunities for surface property optimization. Computer modeling of organic surfaces can promote rational molecular crystal habit design and solid form selection for pharmaceuticals, can lead to improved performance of organic electronic devices, and can help tune the sensitivity of energetic materials.

In the pharmaceutical industry, crystallized drug molecules are commonly compacted into tablets. The ease of production and the mechanical properties of tablets can be largely dependent on the surface properties of the molecular solid form.^1,3,12 It has been demonstrated that the sticking of an active pharmaceutical ingredient (API) to the tablet punch during compaction is a function of the amount of polar groups on the surface of the API crystal¹ and that the surface energy of the particles strongly affects the tablet tensile strength.¹³ In addition, the crystal habit of solid form APIs has been linked to differences in bioavailability.^2,14,15 Understanding the origin of these properties and how to control them is important to drug manufacturers and researchers. The incorporation of surface modeling into the solid form selection process could help save time and money by reducing manufacturing costs and improving the efficacy of the solid form API.³ 

Surfaces control the most important properties of organic solid-state devices, such as the efficiency of charge separation in organic solar cells and the exciton recombination in light-emitting diodes.^16–19 In addition, surface energies determine the growth behavior of organic electronic materials.^4,6 Growing high-quality crystalline organics requires choosing commensurate substrates and tuning growth parameters for specific systems. Modeling the surface properties of organic electronic materials could have a substantial impact on the design and growth of organic electronic devices.^20–23 Surface modeling can lead to rational choices for growth parameters and has been used for determining the equilibrium growth conditions with respect to the crystal shape^21,24 and understanding surface defects in thin film materials.²⁵ 

For energetic materials, lower mechanical and thermal sensitivities would result in improved safety during storage and transport.²⁶ The sensitivity of energetic materials is known to be dependent on their surface properties and characteristics.^7,8,27 It has been observed that crystals of the same energetic material with more angular morphologies have worse impact sensitivity than those that are spherical while controlling for the influences of particle size and void density.^7–9 Surface modeling and morphology prediction of energetic materials can aid in the design of safer solid forms.²⁸ For example, surface modeling could elucidate how the exposure of polar groups on different surfaces leads to differences in surface growth inhibition due the solvent used for crystallization.²⁹ 

Software that can generate molecular crystal surfaces from bulk crystal structures is generally applicable to many areas of research. Despite the importance of surface modeling across industries, the available computational tools for simulations of organic molecular crystal surfaces lag behind the tools available for modeling inorganic surfaces. Popular and widely used open-source Python packages for atomistic simulations, such as ASE³⁰ and pymatgen,³¹ have libraries capable of generating inorganic surfaces by cleaving bulk crystal structures. However, currently, these packages do not work correctly out-of-the-box for generating organic surfaces. For systems containing organic molecules, additional physical constraints must be respected because surface creation cannot cleave through the covalent bonds of any molecule in the crystal. Software for the creation of molecular surface slabs from bulk crystals must recognize molecular units and preserve all covalent bonds during surface creation. An additional problem that arises during surface slab generation is the presence of multiple unique terminations for a given Miller index.^32,33 Here, we present a methodology, based on a graphical representation of molecules, to identify and enumerate all possible surface terminations in molecular crystal surface slabs.

Surface structures may be used to calculate surface energies. Traditionally, the surface energy of molecular crystals has been evaluated by calculating the attachment energy, which is the energy released upon the addition of a new layer of molecules to the surface.^34–36 Attachment energies are typically calculated using inter-atomic force fields. Summation is performed over the atomic interactions within a surface slice and over the atomic interactions between the slice and the bulk crystal. The interactions within the surface slice are then subtracted, leaving the attachment energy. The attachment energy may be coupled to periodic bond chain (PBC) analysis.³⁷ PBC takes into account the impact of steps and kinks to model the nucleation of the bulk molecular crystal.³⁸ 

The surface energies of inorganic crystals are often calculated with ab initio methods, such as density functional theory (DFT), using surface slab geometries.^33,39 Surface slabs are a super-cell of the bulk crystal, oriented to expose a particular facet. In DFT codes that impose three-dimensional periodic boundary conditions, a large vacuum region is added above the surface to prevent spurious interactions between periodic replicas.^39,40 The surface energy, γ, is defined as the difference between the energy of the surface slab, E_Slab, and the bulk crystal, E_Bulk, multiplied by the number of layers in the slab, N, and divided by twice the surface area, A,

In practice, because the computational cost increases with N, the surface energy is converged as a function of the number of layers in the surface slab. The surface energies of inorganic materials obtained using Eq. (1) have been shown to converge poorly with increasing N. This has been attributed to differences in the numerical Brillouin zone integration when the slab energy and the bulk energy are calculated separately.^33,41 This can be corrected by using a very dense k-point grid, which comes at a significant computational cost,⁴¹ or by using a surface-oriented bulk structure.³³ Alternatively, the bulk energy may be evaluated based on the incremental energy of adding one layer to the slab, as suggested by Boettger,⁴² 

Equation (1) may also be rearranged into the following linear relation:

where E_Bulk is the slope and the surface energy, γ, is proportional to the intercept of the line.^43–45 For inorganic surfaces, the linear approximation of the surface energy given by Eq. (3) has been shown to yield the best convergence behavior.^33,46 Studies of the convergence of the surface energies of molecular crystals using ab initio methods remain scarce.²² Therefore, we assess the performance of the linear method and the Boettger method here.

An important difference between molecular crystals and typical inorganic crystals is that molecular crystals are bound by weak dispersion interactions. DFT with van der Waals corrections is a popular ab initio method for modeling molecular crystals. Dispersion-inclusive DFT often yields lattice energies in excellent agreement with the experimentally measured values^47,48 and is the method of choice for molecular crystal polymorph ranking.^49–51 Yet, few studies of the surface energies of molecular crystals have been conducted using ab initio methods.^20,22 The availability of an open-source software for molecular crystal surface slab creation may facilitate more detailed studies of molecular surface properties by ab initio techniques. Here, we compare the Tkatchenko–Scheffler (TS)⁵² pairwise dispersion method to the many-body dispersion (MBD)^53,54 method, which accounts for higher-order dispersion interactions and for the effect of long-range electrostatic screening on the atomic polarizabilities in an extended system.

The computed surface energies can be used to predict the crystal shape, which can be validated by comparison to experimentally grown crystallites.^37,55,56 Under equilibrium conditions, the crystal shape given by the Wulff construction minimizes the surface free energy for a constant volume.⁵⁷ Programs to compute and visualize the Wulff shape require the space group, Miller indices, and surface energies of the crystal structure.^39,58,59 We note that the experimentally observed crystal shapes may deviate from the predicted Wulff shape owing to non-equilibrium growth conditions, interactions with solvents⁶⁰ and additives, and kinetic effects. Accounting for solvent effects on the growth shape requires modeling how the solvent will impact the relative growth rate of each facet⁶¹ or the attachment energy.^28,55,62 DFT may be used in the future to calculate the interaction between molecular crystal facets and solvent molecules to account for solvent dependent crystal morphologies. Here, we restrict the discussion to equilibrium crystal shapes. To validate our results, we have chosen experiments known to produce crystal shapes that are close to equilibrium. This requires growing the crystals from low supersaturation, with a low driving force, from vapor, or from a solvent that does not exhibit significant interactions with the crystal facets.^34,38,63 Such examples may be found in the literature for the pharmaceutical aspirin,^15,64 for the organic semiconductor tetracene,⁶³ and for the energetic material HMX.^65–68 

Here, we present Ogre, an open-source code for generating molecular crystal surface slabs, which interfaces with both ASE and pymatgen. Ogre takes a bulk crystal structure as an input and generates all symmetrically unique surfaces up to a user-defined maximum Miller index. Ogre ensures that molecules remain intact by reconstructing the cleaved molecules during the surface creation process. This is done by matching graphical representations of molecules in the bulk structure with the molecules at the surface. In addition, Ogre considers all possible surface terminations for each Miller plane. Ogre streamlines the computation of surface energies and the prediction of crystal shapes by automatically generating slab models for all required surfaces. Ogre is interfaced with the FHI-aims⁷⁵ code for DFT calculations. It utilizes Pymatgen libraries and Wulffmaker⁵⁹ to create Wulff diagrams. In the following, we provide a full account of the implementation and features available in Ogre. We then use Ogre to calculate the surface energies and predict the crystal shapes of three benchmark systems: the pharmaceutical compound, aspirin; the organic semiconductor, tetracene; and the energetic material, HMX. The results are in good agreement with the experimentally observed crystal shapes.

Ogre is written in Python 3 and uses the ASE³⁰ and pymatgen³¹ libraries. The code is available for download from www.noamarom.com under a BSD-3 license. Ogre takes a bulk crystal structure as an input. Several common structure input formats are supported, including CIF, Vienna ab initio Simulations Package (VASP)⁶⁹ POSCAR, and FHI-aims⁷⁵ geometry.in. Ogre generates surface slabs by cleaving the bulk crystal along a user-specified Miller plane. Ogre outputs surface slab models with the number of layers requested by the user. Additionally, the user may specify the number of unit cells to include in the x and y directions and the amount of vacuum to add in the z direction. An algorithm for generating inorganic surface slabs is already available in ASE.³⁰ In Ogre, we have implemented a modified algorithm for organic surface slab generation that keeps molecules intact when cleaving molecular crystal surfaces.

The workflow of Ogre is shown in Fig. 1. Ogre first uses the ASE library for cleaving the surface. The cleaving algorithm in ASE works well for generating inorganic surface slabs. However, for organic molecular crystals, it may produce unphysical surfaces by cleaving through the molecules and breaking covalent bonds. Ogre identifies and repairs these broken molecules. The code first identifies the molecules comprising the bulk structure using a graphical representation, as described in Sec. II A. Subsequently, Ogre compares the molecules comprising the slab to the bulk and identifies the broken molecules, which are not isomorphic to any molecule in the bulk. The code then repairs the broken molecules, as described in Sec. II B. All possible surface terminations for given Miller indices are identified, as described in Sec. II C. Finally, Ogre outputs the slab structures. Ogre is accompanied by the OgreSWAMP (Smooth Workflow for Accurate Molecular habit Predictions) utilities to streamline surface energy calculations and predict the crystal shape by constructing Wulff diagrams, as described in Sec. II D.

The input structure information is usually stored as three lattice vectors, atomic coordinates, and species. Given a structure, Ogre first assigns atoms to molecules. If the distance between two atoms is within their bond length plus a tolerance value of 0.56 Å, then they are assigned to the same molecule. If two atoms have more than one possible bond length (e.g., a single bond or a double bond), the maximal possible bond length is considered. The bond distance values are taken from Open Babel.⁷⁰ To identify molecules in the bulk, a 3 × 3 × 3 super-cell is constructed (such that molecules that lie across the boundaries of the central unit cell are intact), and a list of molecules in the bulk is obtained by applying the aforementioned criterion to assign atoms to molecules in the central unit cell. A list of molecules in the slab is obtained similarly by constructing a 3 × 3 × 1 super-cell because there is no periodicity in the direction of the surface.

To compare the molecules in the slab to the molecules in the bulk, we use an undirected graphical representation. For each molecule, the nodes represent the atomic species and the edges represent the distance between two atoms.^71–73 For instance, the graphical representation of a C₂H₆ molecule is shown in Fig. 2. This representation facilitates the comparison between molecules because it is invariant to the translation and rotation of the atomic coordinates. Molecular graphs are considered isomorphic if their nodes are the same species, and the difference between corresponding edges is less than 10⁻⁵ Å.⁷⁴ 

Figure 3 shows an example of generating a (111) surface slab of pentacene. The first step is cleaving the structure. To this end, Ogre determines the lattice vectors of the bulk crystal, $a1→,a2→,a3→$⁠, and the Miller indices $h,k,l$ of the target surface. Then, the new basis vectors for the surface model, $v1→,v2→,v3→$⁠, are calculated by ASE. The resulting slab with broken molecules is shown in Fig. 3(a). The broken molecules are identified, as described in Sec. II A, and colored in yellow and red in Fig. 3.

The second step is to make a copy of each broken molecule and translate it by $v3→$ such that the broken molecules on the top surface are repaired. This results in some residual fragments above the top surface and at the bottom surface, as shown in Fig. 3(b). Therefore, the third step is to check for broken molecules again (now the molecules at the top surface are intact and no longer regarded as broken) and identify the residual fragments, which are deleted in the fourth step. This produces a surface slab with the requested number of layers, as shown in Fig. 3(d).

Molecular crystal surfaces can have more than one possible termination, depending on the symmetry of the crystal and its constituents. Here, different terminations are defined as slabs that have the same Miller indices and the same number of molecules but differ in their surface chemistry (e.g., because different chemical groups of the molecules are exposed on the surface). Ogre automatically identifies the possible terminations of a given surface. We do not consider surface reconstructions or defects, such as missing molecules. For example, the (100) surface of form I of aspirin has two possible terminations, as shown in Fig. 4. To find these, Ogre starts from the structure obtained after repairing the broken molecules in the slab, as shown in Fig. 4(a). Ogre stores the initial configuration and then moves the molecules in the top layer down to the bottom of slab by subtracting $v3→$ from the corresponding coordinates. The new configuration is compared with the original configuration. If determined to be different, Ogre stores the current configuration and moves the top molecules downwards again. Otherwise, Ogre outputs the structures of all stored configurations. For aspirin-I (100), Ogre performs one iteration and finds two terminations, as shown in Figs. 4(a) and 4(b). As shown in Fig. 3 and in panels (a) and (b) of Fig. 4, upon surface generation, there is no guarantee that $v3→$ would be perpendicular to $v1→$ and $v2→$⁠. A more convenient choice for adding the requested amount of vacuum is to define $v3̃→$ to be perpendicular to $v1→$ and $v2→$ while maintaining the total volume of the unit cell, as shown in Figs. 4(c) and 4(d) for the two configurations of aspirin-I (100). These are the final surface structures output by Ogre.

Ogre is accompanied by the OgreSWAMP (Smooth Workflow for Accurate Molecular habit Predictions) utilities to streamline the calculation of surface energies and Wulff shapes, as shown in Fig. 5. Given a molecular crystal structure, OgreSWAMP determines the Miller indices of all unique surfaces up to a user-defined maximal Miller index based on the space group symmetries. For each unique surface, Ogre generates slab models for all possible terminations with an increasing number of layers within a user-defined range. DFT calculations are performed to obtain the total energy of each structure. OgreSWAMP calculates surface energies by the Boettger or linear methods depending on user input options and outputs the surface energy convergence plots. The default convergence criterion is that the relative deviation of the surface energy should be within 0.35% with the addition of one layer. If the surface energy is not converged, the user needs to perform DFT calculations for slabs with a higher number of layers. The converged surface energies are then used as input to call WulffMaker or pymatgen to generate the Wulff crystal shape depending on user input options. To rigorously determine the maximal Miller index to be calculated, the user may increase the maximal Miller index until no new surfaces appear in the Wulff shape.

DFT calculations were performed with the all-electron electronic structure code FHI-aims using the light numerical settings and tier 1 basis sets.⁷⁵ The generalized gradient approximation (GGA) of Perdew, Burke, and Ernzerhof (PBE) was employed for the description of exchange–correlation interactions among electrons.^76,77 Structural relaxations were performed for the bulk crystal structures of aspirin, tetracene, and HMX, using PBE with the TS pairwise dispersion method,⁵² until the residual force per atom was less than 10⁻². The number of k-points in each direction was set to the closest integer greater than 24 divided by the norm of the respective lattice parameter. The relaxed structures are in close agreement with experiment, as shown in the supplementary material.

Surface slabs were constructed based on the relaxed bulk geometry with a varying number of layers and a vacuum region of 40 Å. Due to the large size of surface slabs (up to 1596 atoms), single point calculations were performed using PBE+TS and PBE+MBD.^53,54 It has previously been demonstrated that relaxation has a negligible effect on the surface energies of molecular crystals.^23,78 To ensure an approximately uniform k-point sampling, for each slab, we set the k-point grid in the xy plane such that k_i × v_i ≈ 24, ∀i ∈ 1,2. In the z direction, k₃ was set to 1. For instance, if the lattice parameters of a slab are v₁ = v₂ = 12 Å, then the k-point grid is 2 × 2 × 1. Dipole corrections⁷⁹ were included for the slab calculations. Surface energies were calculated using Eq. (1) with the Boettger method [Eq. (2)], and the linear method [Eq. (3)] used to evaluate the bulk energy. Within the localized linear method, for a slab with N layers, the slope was extracted from three points with N − 2, N − 1, and N layers. This leads to an improved convergence behavior. The surface energy convergence plots for all the surfaces studied here are provided in the supplementary material. Representative examples are shown in Sec. IV. Finally, based on the surface energies obtained with the linear method for all non-equivalent surfaces, the Wulff shapes of aspirin, tetracene, and HMX were calculated using pymatgen.⁸⁰ 

Interaction energies (IE) were calculated for the most important intermolecular interactions in the HMX crystal, using PBE+TS and PBE+MBD. First, the four different types of nearest neighbor interactions found in β-HMX were identified. Then, the periodic chains of nearest neighbor molecules were extracted from the relaxed crystal structure. A unit cell was constructed around each interaction such that the lattice vector in the direction of the interaction was the same as in the bulk crystal and a vacuum region of 40 Å was added in the other two directions. The k-grid was defined as above. The interaction energy was computed as the difference between the total energy of the IE simulation and the product of the energy of the isolated HMX molecule and the number of molecules in the IE simulation unit cell.

Aspirin is one of the most important pharmaceutical compounds used to reduce pain, fever, or inflammation. Aspirin has three known polymorphs.^81–83 Here, we focus on form I⁸¹ because of the availability of high-quality experimental micrographs of the crystal shape.¹⁵ Aspirin-I has the monoclinic space group P2₁/c, which leads to 13 unique surfaces with a maximal Miller index of 1: three of the {100} family, six of the {110} family, and four of the {111} family. Of these, the (100) surface has two possible terminations.

Figure 6 shows representative surface energy convergence plots for both terminations of the (100) surface and the (010) surface, obtained using PBE+TS and PBE+MBD with the Boettger and linear methods. Similar plots for all other unique surfaces of aspirin-I are provided in the supplementary material. The Boettger and linear methods converge to the same surface energy values. The linear method exhibits a smoother convergence behavior than the Boettger method, which exhibits some oscillations around the converged value. The oscillatory behavior of the Boettger method has also been reported for inorganic materials.³³ 

For aspirin, the hydrogen bonds between molecules play an important role in determining the surface energies. Cleaving through hydrogen bonds significantly increases the surface energy.⁶⁴ For the (100) surface of aspirin-I, one possible termination, denoted as type I in Fig. 6(a), cleaves through the strong OH⋯O bonds between the carboxylic acid cyclic dimers, exposing the carboxyl groups on the surface. This surface termination has a significantly higher surface energy than the type II termination in Fig. 6(b), which cleaves through weaker intermolecular interactions and exposes the acetyloxy groups on the surface. Similar to the (100) type I surface, the (010) surface, shown in Fig. 6(c), cleaves through several hydrogen bonds, exposing the hydroxy group on the surface, which leads to a higher surface energy. These high-energy surfaces are not expected to be present in the equilibrium crystal shape.

Accounting for the many-body dispersion interactions has been shown to be essential for the correct description of the relative stability and mechanical properties of form I and form II of aspirin.⁸⁴ For aspirin-I, we find that MBD reduces all surface energies by about 30 mJ/m². The TS dispersion method tends to overestimate the lattice energies of molecular crystals compared to the MBD method.⁴⁸ The many-body dispersion contributions and the electrostatic screening included in the MBD method tend to reduce cohesion.⁵⁴ The overbinding of the TS method leads to the overestimation of the strength of inter-layer dispersion interactions and thus of the surface energies.

The Wulff shapes obtained based on surface energies calculated with PBE+TS and PBE+MBD are compared to the experimentally observed morphology^15,64 in Fig. 7. Despite the significant differences in the surface energies produced by PBE+TS vs PBE+MBD, the resulting Wulff shapes are nearly identical in this case. We attribute this to the fact that the relative energies of different surfaces are more important for the Wulff construction than the actual values of the surface energy. The largest facets in the Wulff shape correspond to surfaces that do not cleave through hydrogen bonds, i.e., the (001) and (100) type II surfaces. The smaller (⁠$1¯$10) facets cleave half the number of hydrogen bonds and are therefore more stable than the (100) type I and the (010) surface. Of the facets present in the Wulff shape, (100) exposes the acetyloxy groups on the surface, while (⁠$1¯$10) and (001) both expose the benzene ring on the surface. The results obtained with both PBE+TS and PBE+MBD are in closer agreement with the experimental crystal habit than an earlier study using force-field based attachment energies.⁸⁵ 

Tetracene is a representative π-conjugated organic semiconductor used in organic field-effect transistors⁸⁶ and organic light-emitting diodes.^87,88 Recently, there has been increasing interest in tetracene as a singlet fission material.⁸⁹ Tetracene crystallizes in the space group P$1¯$⁠, which leads to 13 unique surfaces with a maximal Miller index of 1: three of the {100} family, six of the {110} family, and four of the {111} family. In addition, we considered one unique surface of the {210} family and two surfaces of the {211} family (the {200} surfaces are equivalent to {100}). These higher index families contain surfaces with an interplanar spacing above 3 Å, which may be of morphological importance because they are relatively smooth.^63,78

Figure 8 shows the representative surface energy convergence plots for the (100), (001), (⁠$1¯$01), and (20$1¯$⁠) surfaces obtained using PBE+TS and PBE+MBD with the Boettger and linear methods. Similar plots for all other unique surfaces of tetracene are provided in the supplementary material. The Boettger and linear methods converge to the same surface energy values. Similar to aspirin, the Boettger method exhibits oscillations around the converged value for some surfaces [such as the (011) surface, shown in the supplementary material].

In tetracene, the π system has a large effect on the surface energy. The (100) surface, shown in Fig. 8(a), cleaves through intermolecular π–π interactions. The molecules are oriented such that the π electrons are exposed on the surface. This leads to a relatively high surface energy. In comparison, the (001) surface, shown in Fig. 8(b), cleaves through weaker intermolecular interactions. The π-stacking direction is almost perpendicular to the plane, and the CH groups are exposed on the surface. This leads to a significantly lower surface energy.

The elongated shape of the tetracene molecule produces a significant corrugation for some surfaces, which is another decisive factor for the surface energy. The (⁠$1¯$01) surface of tetracene, shown in Fig. 8(c), is particularly jagged with hills and valleys between the molecules on the surface. The (20$1¯$⁠) surface with the molecules similarly oriented, as shown in Fig. 8(d), is much smoother in comparison. Therefore, although the (20$1¯$⁠) surface has a higher Miller index and a smaller interplanar spacing, which would generally lead to a decreased morphological importance, it is, in fact, more stable than the (⁠$1¯$01) surface.

It has been shown that the MBD method yields an improved description of the lattice energies and dielectric constants of acene crystals.⁹⁰ We find that PBE+MBD reduces all surface energies of tetracene by about 50 mJ/m², compared to PBE+TS. Of the materials studied here, this is the most significant difference. This may be attributed to depolarization, arising from electrostatic interactions between molecules in the tetracene crystal, which significantly reduces the van der Waals contribution to the lattice energy obtained with PBE+MBD.⁹⁰ 

The Wulff shapes obtained based on surface energies calculated with PBE+TS and PBE+MBD are compared to the experimentally observed morphology⁶³ in Fig. 9. For tetracene, the differences in the surface energy of different facets are not as large as for aspirin-I, which leads to a more complex crystal shape with more facets present. Overall, the crystal shapes obtained by PBE+TS and PBE+MBD are similar and in good agreement with the experiment. Because (001) is the only surface that is nearly perpendicular to the π-stacking direction, it is the most stable and largest facet. Other prominent facets present in both the PBE+TS and the PBE+MBD Wulff shapes are the (110), (20$1¯$⁠), and (⁠$1¯$10) surfaces. Smaller facets present in both cases are the (100) and (⁠$1¯$11). All of these have relatively smooth surfaces.

The differences between the PBE+TS and PBE+MBD Wulff shapes are due to the presence or absence of several very small facets. The (⁠$1¯1¯$2), (011), (1$1¯$1), and (010) surfaces are present in the PBE+TS Wulff shape and absent from the PBE+MBD Wulff shape. The (101) surface is present in the PBE+MBD Wulff shape and absent from the PBE+TS Wulff shape. These surfaces are relatively rough and lack strong in-plane intermolecular interactions. They may appear in the Wulff construction if their interplanar spacing and surface energy are comparable to other surfaces. However, they are unlikely to be present in experimentally grown crystals owing to kinetic effects that are not taken into account in the equilibrium Wulff model.⁹¹ Such rough surfaces without strong in-plane interactions tend to exhibit a greater density of kink sites.^92,93 These sites are more favorable for molecular attachment and incorporation than the smoother surfaces. Hence, these surfaces are expected to grow rapidly in solution and disappear from the experimental crystal shape. Such kinetic considerations are applied in the periodic bond chain (PBC) method.^37,38 Thanks to the absence of several small facets, the PBE+MBD shape is in better agreement with the experiment than the PBE+TS shape. We note that in Ref. 63, the tetracene crystal has a (11$2¯$⁠) surface that is not present in our predicted morphology. However, Ref. 63 states that under different growth conditions, the (11$2¯$⁠) surface is absent, whereas the (⁠$1¯$11) and (110) surfaces are present. In addition, the (11$2¯$⁠) surface has not been reported in any other studies of tetracene single crystals.^94–97 The crystal habit predicted with PBE+MBD is in better agreement with the experimental structure than predictions obtained previously, using other ab initio methods,²² as well as force-field based attachment energy calculations,⁶³ and comparable to the shape produced by periodic bond chain analysis.²³ 

HMX is an important and commonly used energetic ingredient in various high performance explosives and propellant formulations, thanks to its thermal stability and high detonation velocity relative to other explosives.^98,99 HMX has three polymorphic forms α, β, and σ.^100,101 Here, we focus on β-HMX because it is the most stable form at room temperature and atmospheric pressure and therefore industrially preferred.^102,103 β-HMX crystallizes in the monoclinic space group P2₁/c. The unit cell, illustrated in Fig. 10, contains two molecules, which occupy special Wyckoff positions with inversion center symmetry at the origin and the body center site. In the gas phase, HMX has two conformers, a boat form and a chair form.^104,105 In the solid state, HMX adopts the chair conformation with two axial nitro groups, lying in the plane of the HMX ring [colored in blue in Fig. 10(b)], and two pseudo-equatorial nitro groups, oriented perpendicular to the ring [colored in green in Fig. 10(a)].

The axial and equatorial nitro groups participate in distinct intermolecular interactions, as illustrated in Fig. 11. The corresponding interaction energies (IE) were calculated for the periodic chains of nearest neighbor molecules extracted from the relaxed crystal structure, using PBE+TS and PBE+MBD. In the c direction, there are two types of interactions. The stronger of these is the cooperative chain of axial and equatorial interactions, as shown in Fig. 11(a). The equatorial interactions are along the c direction, and the axial interactions are along the ⟨101⟩ direction. The weaker interaction, shown in Fig. 11(b), involves intermolecular contacts only between axial nitro groups. The equatorial groups participate in cyclic CH⋯O interactions along the a and b directions, as shown in Figs. 11(c) and 11(d), respectively. The interactions in the a direction are significantly stronger than the interactions in the b direction. For all the interactions involving the equatorial nitro groups along the a, b, and c directions, MBD yields a lower IE than TS. For the purely axial interactions along the c direction, both methods give the same IE.

β-HMX, upon relaxation with DFT, was identified as the P-1 space group, which has 13 unique surfaces with a maximal Miller index of 1: three of the {100} family, six of the {110} family, and four of the {111} family. All surfaces have one possible termination. We note that the (111) and (11$1¯$⁠) surfaces are polar such that they have two terminations with the top and bottom layers interchanged. These lead to identical results. Figure 12 shows the representative surface energy convergence plots for the (100), (001), (1$1¯$1), (01$1¯$⁠), (111), and (110) surfaces, obtained using PBE+TS and PBE+MBD with the Boettger and liner methods. Similar plots for all other unique surfaces of HMX are provided in the supplementary material. The Boettger and linear methods converge to the same surface energy values. Similar to aspirin and tetracene, the Boettger method exhibits oscillations around the converged value for some surfaces, such as the (01$1¯$⁠) and (110) surfaces, as shown in Figs. 12(b) and 12(d). The MBD method reduces the surface energies of HMX by about 30 mJ/m², compared to the TS method.

The nitro interactions are the dominant factor in the surface energies of β-HMX. The (100) surface, shown in Fig. 12(a), cleaves through the equatorial nitro interactions, exposing the equatorial nitro groups on the surface. Because these are the strongest interactions in the bulk crystal structure, this leads to the highest surface energy for the {100} family of planes. In comparison, the (001) surface, shown in Fig. 12(b), cleaves through the axial nitro interactions without disrupting the equatorial nitro interactions in the plane of the surface. This leads to a significantly lower surface energy. The (1$1¯$1) surface, shown in Fig. 12(c), cleaves both the equatorial and axial interactions, effectively removing all nearest neighbor interactions for the molecules on the surface. Moreover, the (1$1¯$1) surface is highly corrugated. This causes the (1$1¯$1) facet to have the highest surface energy of all surfaces considered here. The facet with the lowest surface energy is (01$1¯$⁠). This surface preserves the cyclic equatorial nitro group interactions in the ab plane while cleaving the cooperative axial and equatorial interactions and the axial interactions in the c direction. The (01$1¯$⁠) surface has a smooth termination and is the close packed plane of β-HMX. These factors contribute to lowering its surface energy. The (111) and (110) surfaces, shown in Figs. 12(e) and 12(f), respectively, have a similar zigzag motif. Although both the (111) and (110) surfaces cleave the equatorial and axial nitro interactions, the surface energy of (111) is about 20 mJ/m² higher than (110). This may be attributed to the corrugation and polarity of the (111) surface.

The Wulff shapes obtained based on surface energies calculated with PBE+TS and PBE+MBD are compared to the experimental crystal habit,^65,66 as shown in Fig. 13. The Wulff shapes obtained by PBE+TS and PBE+MBD are similar. The most prominent facet is (01$1¯$⁠), which has the lowest surface energy. The second largest facet is the (001) surface, which is relatively stable because it only cleaves through the weak axial nitro interactions. Smaller facets present in the Wulff shape are (110), (10$1¯$⁠), (1$1¯$0), and (010). All of these cleave through equatorial and axial nitro interactions but are relatively smooth. The Wulff shape obtained with PBE+MBD has a more elongated (001) surface, which is in somewhat better agreement with the experiment than the Wulff shape obtained with PBE+TS. The crystal habit predicted by PBE+MBD is in better agreement with experiment than a previous prediction obtained using a spiral growth model⁶⁸ and is comparable to a result obtained with the force-field based attachment energy calculations.⁶⁷ 

We have introduced Ogre, an open-source Python code for generating molecular crystal surface slabs and the accompanying OgreSWAMP utilities for streamlining the calculation of surface energies and Wulff shapes. Ogre can take as an input a bulk crystal structure in various formats and generate surface slab models with the user-specified Miller indices, number of layers, and vacuum space. When cleaving molecular crystal surfaces, Ogre uses a graphical representation to identify molecules and ensure that the molecules at the surface remain intact. Ogre identifies all the non-equivalent surfaces based on the crystal’s space group symmetry and outputs slab structures with all possible terminations for each unique surface. These can be used to calculate the surface energies and predict the crystal shape. Ogre has been successfully applied to three diverse systems: aspirin, tetracene, and HMX.

To converge the surface energy as a function of the number of layers in the slab, we evaluated the Boettger and linear methods. We find that the two methods consistently converge to the same values. The linear method displays a more stable convergence behavior, whereas the Boettger method exhibits oscillations around the converged value in some cases. Therefore, we recommend using the linear method for molecular crystals.

For surface energy calculations and crystal shape prediction with DFT, we assessed the performance of the pairwise TS dispersion method compared to the many-body dispersion method, when paired with the PBE functional. We find that PBE+MBD generally reduces the surface energies of molecular crystals, compared to PBE+TS. This may be attributed to the many-body dispersion contributions and the electrostatic screening included in the MBD method, which correct the overbinding of the TS method. Despite the differences in the values of the surface energies, PBE+TS and PBE+MBD produce similar Wulff shapes because the relative energies of different surfaces are more important for the Wulff construction. We find that the Wulff shapes produced using PBE+MBD are in better agreement with experiment. Therefore, we recommend using the MBD method for surface energy calculations.

The materials studied here are characterized by different intermolecular interactions, which govern the surface energies and the resulting crystal shape. Aspirin is bound by strong hydrogen bonds between carboxylic acids, tetracene is bound by π–π interactions between the aromatic systems along the stacking direction, and HMX is bound by nitro group interactions. In general, the surfaces that cleave through the strongest intermolecular interactions have higher surface energies and are not present in the Wulff shapes. Surfaces that preserve the strongest in-plane intermolecular interactions are relatively more stable. In addition to the strength of intermolecular interactions, surface roughness plays a significant role in the surface energy. Highly corrugated surfaces are less stable than smoother surfaces. Therefore, they are typically absent from the Wulff shape or present as very small facets, which are likely to disappear during crystal growth. Another factor that may contribute to the surface energy is polarity, with polar surfaces expected to be less stable. Owing to the interplay of these different factors, the surface energies of molecular crystals do not necessarily correlate well with the inverse of the interplanar spacing, as shown in the supplementary material.

In summary, we have demonstrated that Ogre is a useful tool for studying molecular crystal surfaces, calculating surface energies, and predicting crystal shape. Ogre is expected to be broadly applicable for designing molecular crystals with optimal shapes and surface properties for various purposes, including pharmaceuticals, organic semiconductors, and energetic materials.

See the supplementary material for a comparison of the relaxed crystal structures to experiment; all surface energy convergence plots for aspirin, tetracene, and HMX; an analysis of the surface energy vs interplanar spacing; information on the computational cost of Ogre; and an analysis of the effect of the surface energy convergence threshold on the resulting Wulff shape.

S.Y. and I.B. contributed equally to this work.

The data that support the findings of this study are available from the corresponding author upon reasonable request. The Ogre code is available for download.¹⁰⁶ 

Funding was provided by the U.S. Army Research Office (ARO) under Grant No. W911NF1810148. This research used resources of the National Energy Research Scientific Computing Center (NERSC), a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.