Ab initio simulation of molecular crystal regrowth of paracetamol from solution

Abstract

The morphology of molecular crystals depends strongly on both thermodynamic stability and the growth kinetics which are themselves dependent on the fine details of intermolecular interactions and challenging to model with ab initio methods. Here, the combination of density functional theory with the effective screen medium-reference interaction site model (DFT/ESM-RISM) is used to study the fast regrowth of a form I paracetamol crystal post-breakage, recently reported by [Bade et al., Mater. Horiz., 2023, 10, 1425–1430]. It is demonstrated that both the thermodynamic and the kinetic factors affecting regrowth are successfully captured by DFT/ESM-RISM with relatively low computational costs. With inclusion of all the externally observed facets, the morphology predicted from thermodynamic considerations alone is found to agree well with observation. Deviation from this morphology is predicated upon inclusion into the model of the fast-growing internal (010) plane, indicating the strong influence of kinetic effects on morphology. The paracetamol molecules at the surface are characterised by unsaturated hydrogen bonds; the resultant strong interaction with the solutes and the solvent significantly altering surface thermodynamics and the structure of the near-surface solvent. For example, the competition between ethanol and solvated paracetamol molecules for the formation of hydrogen bonds is found to reduce the growth rate due to steric hindrance. This effect becomes less prominent for the (010) surface, which presents no broken hydrogen bonds, resulting in a more uniform near-surface solvent structure that facilitates surface growth. As the first attempt to investigate the complicated solid–liquid interface of molecular crystals, this study broadens the applicability of DFT/ESM-RISM. The kinetic mechanisms underpinning the fast regrowth of form I paracetamol post-breakage are qualitatively elucidated, suggesting new strategies for efficient morphology control in molecular crystals.

---

New concepts

As an important purification technique and phase transition process, crystallisation is highlighted in the surface science and fine chemical industries. In this study, we pioneer the use of quantum mechanics/molecular mechanics (QM/MM) to model the crystallisation of paracetamol from solution at the atomic scale, which successfully combines the advantages of both methods to capture the critical thermodynamic and kinetic factors underpinning the delicate and complex interactions at solid–liquid interfaces. Mechanisms of the unique post-breakage growth kinetics recently observed in paracetamol crystals [Bade et al., Mater. Horiz., 2023, 10, 1425–1430] are firstly elucidated, revealing the non-negligible influence of solvent structures during the crystallisation of molecular crystals. The unique role of hydrogen bonds in the selection of adsorbates and the tuning of steric hindrances of the near-surface solvent revealed in this paper offers novel strategies for granting morphological and process controls over crystallisation. This study is also the first attempt to apply the density functional theory with the effective screen medium-reference interaction site model (DFT/ESM-RISM) to investigate surfaces of molecular crystals, which successfully broadens its applicability and offers a prototype for future atomic-scale simulation studies for the finer tuning of crystal morphologies.

Introduction

Crystallisation is an important purification technique in the pharmaceutical industry, where randomly distributed molecules within the fluid accumulate and pack into structurally ordered and compositionally pure crystals. This delicate process can be influenced by thermodynamic and kinetic factors, both of which could potentially dominate the morphology of products and alter their properties. Practically, downstream processes, such as milling and tableting, can largely benefit if crystals with the desired shape and surface properties can be reliably reproduced by crystallisation. Therefore, a thorough understanding of the thermodynamics and kinetics of crystallisation is needed for a precisely controlled manufacturing, and any unpredicted growth phenomenon might lead to costly mitigation.^1,2 It is widely observed in the growth of pharmaceuticals that the morphology of molecular crystals can be sensitive to multiple environmental factors,^3–6 which can be attributed to the influence of the subtle nature of the interactions in molecular crystals. A delicate balance is maintained between the weak intermolecular interactions (e.g., van der Waals dispersions and hydrogen bonds) and the strong intramolecular covalent bonds, providing substantial challenges to both theoretical and experimental investigations.^7,8

Paracetamol is a widely used and important pharmaceutical and has become a prototype system for both process engineering and theoretical chemistry research. Multiple factors influencing the morphology of paracetamol crystals are reported in the literature, including supersaturation,⁶ surface chemistry,^9,10 milling⁹ and impurities.⁵ An unconventional growth kinetic of paracetamol has recently been reported by Bade et al.¹¹ The cleaved form I paracetamol crystal was observed to restore its original shape in saturated ethanol solution with a growth rate approximately three times faster than that of non-cleaved crystals. The dominating factors and microscopic mechanisms of this so-called ‘self-healing’ phenomenon remain unknown. Similar post-breakage regrowth phenomena of aceclofenac has been observed by Jiang et al.¹² in acetone and methyl acetate, suggesting that this is a common phenomenon of both fundamental and commercial interest. It is also a particularly vivid illustration of the peculiar nature of molecular crystals and their crystallisation processes. Simulations at the atomic scale provide a probe for examining the fundamental mechanisms governing this behaviour. In recent years, ab initio simulations based on density functional theory (DFT) have been adopted to study a broad range of properties of molecular crystals. Successful applications include lattice energy,^13–15 vibrational spectra^16,17 and mechanical properties.^18,19 However, DFT-based morphology studies have so far been limited to vacuum morphologies of inorganic materials such as alloys and minerals,^20–23 where stronger interactions exist, enhancing the robustness of predictions of crystallisation morphologies and weakening the influence of solvent effects. For molecular crystals, simulations based on pairwise interactions and molecular dynamics have illustrated the non-negligible influences of solvent effects during crystallisation.^4,24,25 The specific mechanisms, however, are likely to depend upon the interaction model adopted and may not be transferred to configurations for which the models have not been validated. Therefore, a predictive theory, which captures the critical thermodynamic and kinetic processes in the growth of molecular crystals and combines the predictive power of DFT with an appropriate solvent model, would be of great value.

Explicit atomistic models of solvents such as water have been used successfully for interfacial reactions on structurally simple surfaces,^23,26 but the computational cost makes this approach unfeasible for molecular crystals, where, in typical scenarios, much larger simulation cells with lower symmetry must be modelled. Common implicit solvent models, such as the widely adopted polarisable continuum model,²⁷ on the other hand, approximate the solvent as the averaged potential field, making it difficult to simulate the kinetic phenomena arising from non-uniform solvent structures. By iteratively solving the DFT Kohn–Sham equation and a classical solvent correlation function, the DFT/effective screen medium-reference interaction site model (DFT/ESM-RISM)^28–30 achieves a balance between cost and accuracy. In DFT/ESM-RISM, the solvent structure and orientation of solvent molecules are implicitly described by particle distribution functions. Therefore, a non-uniform potential commensurate with the solvent structure arises and is applied to the electronic structure of the solute. As a hybrid quantum mechanics/molecular mechanics (QM/MM) method, DFT/ESM-RISM has been successfully adopted to study highly heterogeneous and dynamic electrode systems,^31–33 where fields originating from solvation effects and electrical potentials lead to non-negligible influences on the thermodynamics and kinetics of solid–liquid interfaces.

In this work, DFT/ESM-RISM is adopted to elucidate the microscopic mechanisms underpinning the recently reported ‘self-healing’ phenomenon of paracetamol crystals.¹¹ The monoclinic form I paracetamol crystal studied is the thermodynamically most stable phase under ambient conditions and the phase grown from saturated ethanol solution. Multiple facets of form I paracetamol are modelled as slabs, including the (001), (011), (110) and (201̄) facets exposed in the experimental morphology^4,9,10,18,34 and the fast-growing internal facet (010).¹¹ Surface formation energies with and without implicit ethanol are computed and used for Wulff constructions³⁵ to investigate the influence of the solvent on the thermodynamically equilibrium morphology. The structural and the electronic properties of the solvated systems are analysed to understand the kinetic process during crystallisation. The model predicts distinct particle and charge distributions at surface–solvent interfaces driven by different distributions of the unsaturated hydrogen bonds exposed on surfaces. The competition between solvent and solute molecules at different adsorption sites is analysed to qualitatively elucidate the dominant factor in the ‘self-healing’ phenomenon post-breakage. This work reports on the first extensive application of the DFT/ESM-RISM for the simulation of crystallisation in molecular crystals and thus broadens its application scenario significantly.

Computational details

DFT/ESM-RISM calculation

All calculations reported in this work are based on DFT in the planewave, projected augmented wave (PAW) method combined with ESM-RISM as implemented in Quantum Espresso 7.1.^30,36,37 As discussed above, the ESM-RISM method provides a suitable balance between an unfeasible fully atomistic model of the solvent and an oversimplified continuum model which neglects the important molecular-scale effects near surface. The B86bPBE functional is used for electron exchange–correlation, and the exchange-dipole moment (XDM) method is adopted for van der Waals dispersion. B86bPBE-XDM has been successfully applied in calculations of both structures and lattice energies of molecular crystals, as reported by several benchmarking studies.^15,38–40 The B86bPBE functional is also used to generate PAW pseudopotentials for H, C, N and O (Section S1.1, ESI†). The kinetic energy cutoffs for wavefunctions and charge density are set to 100 Ry and 600 Ry, respectively. A 2 × 2 × 2 k grid generated by the Monkhorst–Park method is used to sample the Brillouin zone of the bulk crystal. Convergence tests are performed on both the bulk crystal and the slabs to ensure that the adopted parameters limit the error of total energies within 1.2 × 10⁻⁴ Ry per molecule. The k grids of 2D periodic slabs employed to model the surfaces are adapted to best match the sampling of the bulk and are reported in Table S1 (ESI†). For self-consistent field (SCF) calculations, the projected preconditioned conjugate gradient (PPCG) diagonalisation algorithm is used due to its robustness. The total energy per cell is converged to 10⁻⁸ Ry during geometry optimisation and 10⁻⁹ Ry during SCF total energy calculations. A full geometry optimisation is performed on the bulk crystal with energy change between two optimisation steps smaller than 10⁻⁵ Ry per cell and interatomic forces smaller than 5 × 10⁻⁴ Ry Bohr⁻¹. For slab models, convergence thresholds commensurate with those of bulk crystal are used; the details are summarised in Table S1 (ESI†).

The ‘Vacuum-Slab-Vacuum’ ESM method is used with the 2D periodic ‘Right-Hand’ Laue-RISM model to simulate the solvated surfaces.³⁰ The length of the extensive solvent region is set to 50 Å to ensure convergence of the solvent correlation function (Fig. 1(a)). Parameters of the OPLS-AA force field⁴¹ are used to model the dispersion interactions between solvent and surfaces while the implicit solvent molecules are parameterised with the same force field.⁴² The temperature of the solvent is 300 K. The Kovalenko and Hiratas (KH) closure equation²⁸ is applied. The kinetic energy cutoff of the reciprocal space expansion of the solvent correlation function is 400 Ry, as recommended in Quantum Espresso.³⁰ To simulate the competing adsorption of ethanol and paracetamol in solution, the implicit paracetamol model is parameterised with the OPLS-UA force field to reduce computational load.⁴³ As required by the model, a dilute solution is used, with a paracetamol density of 0.01 g cm⁻³. The MOL file used in Quantum Espresso and other ESM-RISM parameters are documented in Section S1.2 (ESI†).

Fig. 1 (a) The ‘Vacuum-Slab-Solvent’ model of DFT/ESM-RISM. The plane-averaged electrostatic potential is plotted to indicate 2D periodicity. The green slashed area is the buffer region of Laue-RISM. As an example, (010) slab is used. (b) Form I paracetamol crystal grown from ethanol solution by slow solvent evaporation.¹¹ The internal (010) facet is indicated by the red dashed line. (c) and (d) Wulff constructions (c) without and (d) with the (010) facet.

Modelling

For the monoclinic form I paracetamol crystal, various space group settings have been reported in the literature. In this work, P12₁/a1 is adopted to represent the primitive cell.⁴⁴ The neutron diffraction resolved structure HXACAN13 from the Cambridge Crystallographic Data Centre (CCDC) is referenced.⁴⁵ Approximately 1% error of the optimised bulk lattice parameters against the experimental data is achieved (Table 1).

Table 1 Experimental and optimised bulk lattice parameters

Source	a/Å	b/Å	c/Å	β/Å	Volume/Å^3
Experiment	12.667(4)	9.166(3)	7.073(3)	115.51(2)	741.2
This work	12.7681	9.0304	7.0025	114.160	736.7
% Error	0.80	−1.48	−1.00	−1.17	−0.61

---

Slabs of (001), (010), (011), (110) and (201̄) facets of form I paracetamol are created as bulk cuts with minimal surface unit cell and finite thickness using the optimised bulk lattice parameters. The GDIS software⁴⁶ provides a slab cutting algorithm that preserves the integrity of molecules crossing Miller planes when generating slabs representing different surface terminations. Convergence tests are performed on pre-optimised slabs to obtain the optimal number of layers and favourable terminations based on the total energy per formula. For all slabs, a vacuum layer of 25 Å is added along the non-periodic direction. The top two layers of molecules are relaxed, while the bottom layers are fixed to represent the bulk phase (Fig. 1(a)). Other structural parameters of slabs are summarised in Table S2 (ESI†). VESTA⁴⁷ and XCrySDen⁴⁸ are used to visualise atomic structures and 3D grid data.

To characterise the influences of solvation on dielectric properties of surfaces, the change of surface dipole moment is computed by eqn (1):⁴⁹

The product of z and the plane-averaged charge density difference (Δ[small rho, Greek, macron]) between ethanol and vacuum is integrated over the surface region, whose boundaries z₁ and z₂ are defined by Δ[small rho, Greek, macron]. In this work, the surface region is defined as |Δ[small rho, Greek, macron]| ≥ 2.5 × 10⁻⁵ e Å⁻³. z₁ and z₂ of different surfaces are listed in Table S2 (ESI†).

Wulff construction

Wulff construction predicts the thermodynamic morphology of a crystal by assuming that the length of a facet norm from the geometric centre of a crystal is proportional to its surface free energy.³⁵ As an approximation of surface free energy (γ), the surface formation energy (E_surf) derived from total energy is popularly used in DFT-based Wulff constructions,^20–23 which consists of the cleavage energy (E_cle) and relaxation energy (E_rlx):

γ ≈ E_surf = E_cle + E_rlx (2)

Cleavage energy is the energy needed to expose a facet of the bulk. Since two surfaces are exposed simultaneously during cleavage, E_cle is calculated as:

E_cle = (E^unrlx_slab − nE_bulk)/2A (3)

E^unrlx_slab is the total energy of the unrelaxed slab. n is the number of molecules and A is the cross-sectional area of the slab along its non-periodic direction. E_bulk is the total energy per formula of the bulk crystal.

Relaxation energy characterises surface relaxation due to exposure. Since the bottom layers are fixed for all slabs, E_rlx is calculated by:

E_rlx = (E^rlx_slab − E^unrlx_slab)/A (4)

E^rlx_slab is the total energy of the relaxed slab.

In this work, the surface formation energies of slabs in vacuum are used to approximate the surface free energies between solid and saturated vapor (γ_SV). Similarly, the surface formation energies of slabs in ethanol are used to approximate the surface free energies between solid and liquid (γ_SL). It has been widely reported that the solvation free energy (Δμ_solv) by ESM-RISM highly depends on parameterisation of force fields and might largely deviate from experimental data.^30,50,51 Therefore, for Wulff construction, Δμ_solv of slabs solvated in ethanol are subtracted to get their DFT total energies (E^rlx_slab and E^unrlx_slab), which also keeps the level of theory consistent when calculating γ since entropic term is not included in E_bulk either.

E_rlx = F − Δμ_solv (5)

Results and discussions

Morphology prediction

Surface formation energies and their components are computed and listed in Table 2. As expected for the observed cleavage plane, the (010) surface has the lowest surface energy amongst those computed, in agreement with the experimental results.¹¹ This stability originates from the packing motif of form I paracetamol, where a corrugated network of hydrogen bonds parallel to the (010) surface is formed. Interlayer interactions perpendicular to the (010) surface are therefore dominated by weak van der Waals dispersion, leading to the lower E_cle. Differences in E_rlx should be attributed to the varied surface geometry and chemistry before and after solvation, which is discussed in detail in the following sections.

Table 2 Surface energies (unit: mJ m⁻²)

(hkl)	Vacuum			Ethanol
	Ecle	γSV	Ref. γSV^9	Erlx	γSL	Ref. γSL^a
a Derived from Owens–Wendt method^52 and experimental data published by Heng et al.^9 and Panzer.^53 Details are reported in Section S2 (ESI).
(001)	267.36	214.71	72.4	−44.31	223.04	16.78
(010)	112.59	109.00	52.1	−0.66	111.92	7.53
(011)	228.86	191.35	66.5	−29.08	199.78	13.48
(110)	197.20	178.45	54.4	−15.44	181.76	7.37
(201̄)	148.97	138.49	62.4	−7.15	141.82	11.10

---

Multiple factors may contribute to the disagreement between the experimental and theoretical energies, including the computationally prohibitive entropy and thermal expansion (eqn (2)).^54,55 In addition, the corrugated surface morphologies increase the contact area between the cleaved surfaces, which is larger than the conventionally adopted cross-sectional area. The ratio between solvent accessible area and the cross-sectional area was adopted to correct solvent attachment energies,^24,25,56 which, however, might not be applicable here since RISM is based on implicit solvents. Nevertheless, one may expect these factors only exert marginal influences and might systematically affect all surfaces considered, which is reflected in the fact that the experimental and theoretical rankings of γ_SV and γ_SL agree well; this is important for morphology prediction. For γ_SV:

DFT/ESM-RISM: (001) > (011) > (110) > (201̄) > (010)

Experiment: (001) > (011) > (201̄) > (110) > (010)

For γ_SL:

DFT/ESM-RISM: (001) > (011) > (110) > (201̄) > (010)

Experiment: (001) > (011) > (201̄) > (010) > (110)

In both cases, the (110) facet exhibits the greatest discrepancy between model and observation. This reconfirms the result and the explanations by Ristic et al.,⁶ that the discrepancy originates from the unique growth mechanism of the (110) surface, where surficial microstructures like voids dominate the growth. Influences of thermodynamic factors such as surface energy could be marginal in this scenario considering the large deviation from the idealised ‘layer-by-layer’ growth, stressing the necessity of studying growth kinetics at atomic scale. In addition, origins of the morphological instabilities of the (110) surface are discussed in detail in Section S3 (ESI†), in a similar way as presented in the main text for the (010) surface.

The thermodynamic morphology of form I paracetamol is predicted using Wulff construction and the computed surface formation energies. By including only the externally observed surfaces, i.e., (001), (011), (110) and (201̄), the morphology predicted (Fig. 1(c)) agrees well with the experimental morphology (Fig. 1(b)). The predicted morphology is significantly altered when including the low energy internal (010) facet (Fig. 1(d)). However, this facet has never been experimentally observed and rapid regeneration is reported when exposed to saturated ethanol solution.¹¹ This strongly suggests that the growth of the (010) surface is dominated by kinetic factors.^35,57 Further support for this proposal is provided by using the measured γ_SL to produce a Wulff construction. This produces a similar morphology in which a high surface area ratio of the (010) facet is present, while the morphology constructed without the (010) facet reproduces the observed morphology (Fig. S1, ESI†). This leads directly to the conclusion that to reliably predict the growth morphology, one must analyse the effect of structural and chemical factors on different paracetamol–ethanol interfaces and probe the kinetics of crystallisation.

Solvation of the paracetamol molecule

Molecules are the ‘building blocks’ of molecular crystals, whose integrity remains unaffected during crystallisation. Solvation of the paracetamol molecule in ethanol is discussed first to understand the dominating ethanol–paracetamol interactions during solvation. 3D-RISM, the three-dimensional periodic version of RISM,⁵⁸ is used for a paracetamol molecule solvated in the implicit ethanol solvent. A cubic cell of 25 × 25 × 25 ų is used to prevent interactions with periodic images of the solute. The mole ratio between paracetamol and ethanol is 0.0062, significantly lower than the solubility of paracetamol in ethanol at 298 K.⁵⁹

In Fig. 2(a) the Löwdin atomic charge difference between solvated and gaseous paracetamol molecules is presented. Due to solvation effects, significant charge redistribution within the molecule is observed between atoms on the ‘upper’ and ‘lower’ sides. The amine H, which forms the NH⋯OH hydrogen bond in form I crystal, has the largest charge depletion of −0.0124 e. Comparatively, the acetyl O, which is the acceptor of the OH⋯O=C hydrogen bond, has the highest charge accumulation of 0.0478 e. The spatial charge density difference (Fig. 2(b)) agrees with the atomic charge difference. A relatively uniform distribution of charge depletion regions is formed around the charge depleted atoms in Fig. 2(a). The distribution of charge accumulation region is more localised. A significant charge accumulation region is formed around the acetyl O and its neighbouring aromatic H, indicating the enhanced interaction between two atoms. The distance between two atoms is reduced from 2.24 Å to 2.22 Å after solvation, still longer than the typical length of hydrogen bonds. Smaller charge accumulation regions are formed around the hydroxyl O, N and methyl C. The charge accumulation and depletion regions at the C=O bond and O–H bond show that both bonds become increasingly polar due to solvation. The clear division between the charge accumulation and depletion regions throughout the molecule increases its dipole moment by 1.34 Debye along almost the same direction as in vacuum (Fig. S2a and b, ESI†), indicating that solvation does not significantly change the distribution patterns of charge and potential across the molecule, but enhances their magnitudes. Besides, the overall electrostatic potential of the paracetamol–ethanol solution clearly shows the solvent screening effect, where the induced solvent field significantly counterbalances the solute field, reducing the overall polarity (Fig. S2c–f, ESI†).

Fig. 2 Difference maps of (a) Löwdin atomic charge and (b) spatical charge density distribution during the solvation of paracetamol. Yellow and blue isosurfaces respectively denote charge accumulation and depletion regions in (b), whose values are ±1.5 × 10⁻⁴ e Bohr⁻³. Particle density distributions of (c) O and (d) H^oh of ethanol. The isosurface value is 1.7 ρ/ρ_bulk. Implicit particles plotted are labelled in the inset.

Particle density distributions of ethanol are analysed to understand the origin of charge redistribution and solvent screening. In Fig. 2(c), the ethanol O atom forms an extensive solvation shell around almost all the H atoms of paracetamol, showing that interactions between two types of atoms induce the widespread charge depletion near the H atoms of paracetamol. A moderately increased concentration is observed around H atoms of the –NH and –OH groups. In Fig. 2(d), the ethanol H^oh atom forms discrete accumulation regions, with the highest concentration around the acetyl and hydroxyl O atoms. Furthermore, on the side of charge depletion, the H^oh rich regions are ∼1 Å outside the O rich region, which is consistent with the typical O–H bond length and can be explained by the orientation of ethanol molecules. Among all the particles plotted in Fig. S3 (ESI†), the highest accumulation of 8.4 ρ/ρ_bulk is reported for H^oh around the acetyl O, which significantly enhances local charge density and induces a strong solvent potential. Comparatively, H¹, H² and C² form almost uniformly distributed solvation shells around the molecule. A moderately increased density of C¹ is observed, whose pattern is similar to that of H^oh because of the lower electronegativity compared to O. Particle distributions indicate that the ethanol solvation shell around the paracetamol is non-uniform and oriented, making dimers with a certain configuration preferred and leading to a counterbalancing solvent field. Among various competing intermolecular interactions, the interaction between ethanol H^oh and acetyl O is the most prominent one. Other important interactions include those between ethanol H^oh and hydroxyl O, ethanol O and amine H, ethanol O and hydroxyl H.

Solvation of paracetamol surfaces

According to the analysis above, solvation in ethanol induces a significant redistribution of the charge within the paracetamol molecule and alters its properties. Since such effect is non-uniform and specific to certain functional groups, it is necessary to investigate and compare solvation effects on different surfaces, where various functional groups are exposed, leading to distinct physicochemical properties. In Fig. 3(a), the plane-averaged charge density difference between solvated and non-solvated surfaces (Δ[small rho, Greek, macron]) reveals a significant charge redistribution on the (001) surface. Δ[small rho, Greek, macron] of (110) is slightly higher, but comparable to those of (011) and (201̄) surfaces, while Δ[small rho, Greek, macron] of the (010) surface is the least prominent. On the (001) surface, the solvent electrostatic potential (V_solv) decreases, inducing a surface dipole pointing from the surface to the solution. For other surfaces, however, backward surface dipoles are induced. Large oscillations of V_solv on (001), (011) and (201̄) surfaces indicate a stronger surface–solvent interaction and the ordering of solvent structure. Comparatively, surface–solvent interactions of (010) and (110) surfaces are weaker, probably leading to more uniform solvent structures.

Fig. 3 (a) Plane-averaged charge density difference (Δ[small rho, Greek, macron]) and electrostatic potentials of solvent (V_solv) and solute (V_solu) along z axis. The grey dotted line denotes the boundary of DFT/ESM-RISM cell. Slabs in the background denote atomic positions. Charge density differences between solvated and non-solvated (b) (001) and (c) (010) surfaces (colour coded as Fig. 2(b)). The isosurface value is ±3 × 10⁻⁴ e Bohr⁻³.

Spatial distributions of charge density difference (Δρ) are illustrated in Fig. 3(b), (c) and Fig. S4 (ESI†). In general, strong and extensive charge accumulation regions form around O atoms near surfaces, which have minimal steric hindrances and strong interactions with solvent. Charge depletion regions are more scattered across the exposed hydrogen bond donor atoms and H atoms attached. The largest Δρ is noticed around atoms with unsaturated hydrogen bonds, indicating the probable formation of hydrogen bonds between paracetamol and ethanol. Interestingly, electrostatic properties of surfaces are largely influenced by surface geometries. Δρ of the (001) surface can be divided into accumulation regions on the top and depletion regions beneath. Therefore, a large surface dipole change (Δp) perpendicular to the surface is induced after solvation (Table 3). For (011), (110) and (201̄) surfaces, such divisions are ambiguous, so Δ[small rho, Greek, macron] and Δp are reduced due to plane-averaging. The negative Δp indicates that, overall, H atoms dominate surface properties, which becomes electron deficient through the interaction with the hydroxyl O of ethanol. The (010) surface has minimal interactions with solvent since its hydrogen bonds remain intact. Δρ shows that surface–solvent interactions are limited around the acetyl group on the top with much smaller amplitudes. Different from other surfaces, the (010) surface might have a more uniform solvent structure across its surface due to the smaller Δρ, while in other cases, localised solvent structures parallel to surfaces might still exist.

Table 3 Change of surface dipole moments

(hkl)	(001)	(010)	(011)	(110)	(201̄)
a Oriented along z axis.
Δp^a (×10^−2 Debye)	8.90	−1.58	−3.33	−5.81	−2.53

---

To learn how surface–solvent interactions affect the electronic structures at the interface, in Fig. 4, density of states (DOS) and integrated local density of states (ILDOS) of the (010) surface are plotted. Projected DOS of the molecule M1 shows that both the valence band maximum (VBM) and the conduction band maximum (CBM) of M1 move to slightly lower energies. This increases the electron affinity of M1 and induces near-surface charge depletion regions and the dipole pointing towards the surface as previously analysed. ILDOS plots in Fig. 4 shows that both VBM and CBM consist of π orbitals and are similar to VBM/CBM of the form I crystal (Fig. S5, ESI†), which is also true for other surfaces (Fig. S6–S9, ESI†), indicating the limited interactions between surfaces and solvents in all the cases.⁶⁰ ILDOS of VBM shows that M1 has negligible contributions due to its interaction with the solvent, and contributions from M2 of the top layer are also reduced. ILDOS of CBM is almost limited to M1 only, with the largest isosurfaces around acetyl and amine groups. Due to the accumulated charge around the top-layer acetyl O, on the (001) surface, the VBM and CBM of the top molecule M1 are shifted to higher energies (Fig. S6, ESI†), enhancing the peak in the valence band and quenching peaks in the conduction band. ILDOS plots show that the surface VBM lies in the top layer molecules, while the CBM lies in the molecules beneath. Therefore, the solvent-induced surface dipole pointing from the (001) surface to solvent can be attributed to the increased valence electron states of M1. VBM states of the top layer acetyl O merit particular attention, as broken OH⋯O=C hydrogen bond and weak steric hindrance might further promote the interaction between O and ethanol H^oh, which significantly influences the solvent structure around the paracetamol molecule (Fig. 2(d)). CBM states increase at M2 due to the exposed hydroxyl group, where interactions between H and ethanol O are enhanced due to the broken OH⋯O=C hydrogen bond. (011), (110) and (201̄) surfaces are illustrated in Fig. S7–S9 (ESI†), with detailed analysis in Section S3 (ESI†). In all cases, moderate solvation effects are observed around molecules near surfaces, while DOS and ILDOS of molecules below quickly converge to bulk states. Surface structures also have significant influences. In particular, the exposed acetyl O brings significant local effects which might differ from the averaged properties of surfaces.

Fig. 4 (left) Total density of states (DOS) and its projections onto the top two layers of molecules (M1–4) on the (010) surface. The corresponding molecules are indicated by the grey dashed lines. VBM and CBM peaks integrated for the integrated local DOS (ILDOS) are marked in purple and red, respectively. Energy is aligned to vacuum. (right) ILDOS of VBM and CBM. Energy ranges of integration are annotated below the figures.

Implicit adsorption

Besides the electronic properties of surfaces, the surface–solvent interaction also alters the solvent structure. A highly localised particle density distribution is expected considering the non-uniform surface structure and chemistry. The solvation-induced ordering of the solvent structure probably plays a crucial role in the activation or passivation of certain adsorption sites, leading to distinct growth rates. Therefore, particle distributions of the implicit solvent (ethanol) and solute (paracetamol) are analysed in this section.

To illustrate the response of the entire system, the plane-averaged distributions of ethanol and paracetamol on different surfaces are plotted in Fig. 5(a) by summing the densities of their constituent atoms and dividing the sum by its counterpart in the bulk solution. As indicated by the completely overlapped plots of pure ethanol (denoted by (p)) and ethanol in the paracetamol–ethanol solution (denoted by (s)), the distribution of ethanol is robust at the low concentration of paracetamol discussed here. Surprisingly, paracetamol has higher average concentrations than ethanol on all surfaces considered. The highest concentration is observed on the (001) surface, apparently in conflict with the observed fast growth of the (010) surface, suggesting that macroscopic supersaturation might not be sufficient to understand ‘self-healing’, and the molecular-level variations of solvent structures need to be elucidated.

Fig. 5 (a) Plane-averaged particle distributions of ethanol (p, pure), ethanol (s, solution) and paracetamol (s, solution). The grey dotted line denotes the boundary of DFT/ESM-RISM cell. Slabs in the background denote atomic positions. (b) Particle distributions of ethanol H^oh, O and paracetamol H^o and O^h near the (010) surface. Implicit particles plotted are labelled in the inset. Semi-transparent molecules on the right column indicate the approximate positions of the adsorbed paracetamol. Atoms forming hydrogen bonds in adsorbates are marked by the blue dotted circle. The corresponding adsorption site on the surface is marked by the red dotted circle.

In Fig. 5(b), distributions of ethanol H^oh, O and paracetamol H^o, O^h near the (010) surface are illustrated. These particles are chosen due to their strong interactions with the acetyl and amine groups of M1, where prominent solvation effects are characterised by Δρ and DOS/ILDOS (Fig. 3 and 4). For both ethanol and paracetamol, strong interactions between surface O and solution H give rise to significantly localised particle accumulations near the surface. The similar chemical properties of ethanol H^oh and paracetamol H^o lead to competitions around acetyl O on the top of the surface, where the larger ρ(H^o) is observed. It is noteworthy that localised H^o accumulation does not necessarily promote the growth rate, as the (010) surface might have a nearly ideal ‘layer-by-layer’ growth due to the parallel hydrogen bond network of form I crystals.^8,61 Strong localised interactions with surface atoms might disrupt the ordering of paracetamol molecules in solution. However, interactions between surface H and solvent O are much weaker, as indicated by the nearly uniform distributions of ethanol O and paracetamol O^h across the solvent accessible surface. The higher ρ(O^h) indicates that the uniform paracetamol growth layer also has higher affinity than ethanol molecules on the (010) surface.

For (001), (011), (110) and (201̄) surfaces, the sequence and orientations of adsorbed paracetamol can be inferred by the broken hydrogen bond and their positions in the bulk, since the adoption process is largely dominated by hydrogen bonds, as revealed by previous analysis. Therefore, it is convenient for the ‘growth layer’ to be considered as split into individual molecules with labels to indicate their orientations (Fig. S10–S12, ESI†). For the (001) surface (Fig. S10a, ESI†), M2 should be adsorbed first to form the lower OH⋯O=C bond. However, more significant accumulations of ethanol H^oh and M1 H^o are observed near the M2 adsorption site. The accumulation of ‘M1-oriented’ molecules causes steric hindrances for ‘M2-oriented’ molecules near the adsorption site. In addition, ethanol O, rather than M2 O^c, has the higher concentration around the hydroxyl H of (001) surface, which is probably due to the weaker steric hindrance of smaller ethanol molecules. The case of the (001) surface reveals a dilemma of surface growth dominated by multiple hydrogen bonds, that the adsorption site close to surfaces is probable to induce a stronger accumulation of adsorbates due to its larger solvent accessible area, but it is the lower lying adsorption site that corresponds to the lower ‘correct’ adsorbate, if lengths and orientations of the two hydrogen bonds are similar. To follow the idealised sequence of adsorption, the adsorbate has to firstly deposit at the lower site with smaller solvent accessible area. That probably brings strong steric hindrances to the adsorbate and limits its diffusion. This mechanism is prominent for the (001) surface due to its relatively simple surface geometry and the large distance between the hydrogen bond acceptor (acetyl O) and donor (OH) sites. Detailed analysis of (011), (110) and (201̄) surfaces is included in Section S3 (ESI†). In all cases, the adsorption is influenced by steric hindrances whose origin depends on the surface chemistry. When the stronger OH⋯O=C hydrogen bond is broken, paracetamol near the top prevents adsorption to the lower sites. When the weaker NH⋯OH hydrogen bond is broken, the interaction between ethanol and surface is stronger, making ethanol the major contributor of steric hindrances.

Explicit adsorption

By analysing the averaged implicit particle distributions, the adsorbate is approximated as the field near surfaces, indicating its average orientation after equilibrium. To further inspect the kinetic process, the energy profile of an adsorption event can be scanned by varying the distance between the explicit paracetamol molecule and the surface. Due to the limited computational resource, the structural and computing parameters of the explicit adsorption models are slightly adjusted, which are believed to be sufficient to capture the key features (Section S1.2, ESI†).

The scanned adsorption energy profiles of paracetamol onto the (010) surface are plotted in Fig. 6(a). The adsorption energy (ΔE) is the total energy difference between the current system and the one where the adsorbate is the farthest from the surface. The minimum of ΔE corresponds to the equilibrium adsorption position inferred from the periodicity of form I crystal. This suggests that adsorption and recrystallisation are thermodynamically favourable in ethanol solution. Near the clean (010) surface, the larger ΔE of the lower M2 agrees with the ‘idealised sequence’, reducing steric hindrances during adsorption. The minima of ΔE of M1 and M2 differ by 23.40 kJ mol⁻¹, increasing selectivity and reducing competition between differently orientated molecules (i.e., M1 and M2). The plane-averaged solvent structures during the adsorption of M1 and M2 are visualised in Fig. 6(b) and (c), where scenes are compared between the equilibrium adsorption positions and the top of energy barriers. An adsorption process with negligible influences in the solvent structure is identified for the lower M2, making supersaturation the only controlling factor. Only the density of ethanol H^oh increases vaguely when M2 is on the top of barrier due to its strong interactions with acetyl O of the adsorbate, which is, however, inadequate to explain the high energy barrier of M2. It is probable that the kinetic diffusion of M2 cannot be fully depicted with a simple energy barrier obtained by DFT/ESM-RISM, which describes the binding between the solute and the surface. Since solvent molecules are not explicitly modelled with ESM-RISM, the energy of removing the solvent between solute and surface can only be indirectly captured with total energies perturbed by the variations in solvent structure and, consequently, in potential field. Although an accurate calculation of energy barriers requires more advanced techniques, such as metadynamics,²⁶ the particle density distribution probed with ESM-RISM provides a qualitative description.

Fig. 6 (a) Adsorption energy profiles onto the (010) surface. The z coordinate of molecular mass centre is used to indicate the position of adsorbate. Plan-averaged solvent particle density distributions of (b) M2 and (c) M1 adsorbed onto the clean surface, and (d) M1 adsorbed at the presence of M2 (M1@M2). Labels of molecules are consistent with Fig. 4. Adsorbates at the equilibrium position (Eq.) and the top of energy barrier (Top) are analysed. Implicit particles plotted are labelled in the inset.

Compared to M2, a non-negligible ethanol accumulation layer between M1 and the surface is observed when M1 is located at the equilibrium position or the top of barrier, indicating a strong steric hindrance of near-surface ethanol. The steric hindrance significantly limits the diffusion of M1, while the diffusion of molecules adopting the M2 orientation is free of it. This selective effect promotes adsorption in the ‘idealised sequence’ on (010) surface, leading to its fast growth behaviour. In comparison, the intermediate ethanol layers are present on other surfaces owing to the unsaturated superficial hydrogen bonds and the resultant increase in surface–solvent interactions, limiting their growth rates. Detailed analysis of other surfaces is available in Section S3 (ESI†).

To understand how hydrogen bond promotes the growth of (010) surface, the ‘step 2’ adsorption energy profile of M1 at the presence of a pre-adsorbed M2 (M1@M2) is plotted in Fig. 6(a). A substantial drop in ΔE is reported when M1 approaches its equilibrium position, while the energy barrier is slightly pushed outwards. The energy drop originates from the formation of hydrogen bonds, suggesting that the strong interactions between paracetamol overtake the paracetamol–ethanol interactions. It also reveals that the rate-limiting step during the growth of the (010) surface is the adsorption of M2, where the weak van der Waals dispersions dominate and only marginal steric hindrance is imposed as previously analysed. The most remarkable difference in particle density distributions of M1@M2 (Fig. 6(d)) is the absence of the near-surface ethanol layer when it is adsorbed to the equilibrium position, which can be attributed to the completion of the hydrogen bond network, leading to uniform surface chemistry. When M1@M2 is located on the top of energy barrier, the ethanol layer is pushed away from the surface with slightly reduced density and spatial range, probably due to the pre-adsorbed M2. This observation is consistent with the variations in the energy barrier of M1@M2, indicating that the barrier might originate from repulsions against the ethanol layers that exist between the adsorbate and the surface.

Conclusions

Based on Wulff construction and surface formation energy, the ESM-RISM method is successfully applied to predict the equilibrium morphology of crystals grown from solution. This method is proven to be sufficient to capture the critical thermodynamics underpinning the growth of form I paracetamol crystal in ethanol. ESM-RISM also exhibits its potential in studying solvent structure and growth kinetics by solving both the charge and particle distributions in the solvent. The structured ethanol shell is formed around the paracetamol molecule, which enhances the dipole of paracetamol by induction and screens the total potential. Distributions of charge and particle densities reveal that the interaction between acetyl O of paracetamol and hydroxyl H of ethanol dominates solvent structures and the electronic properties of surfaces. The complex surface–solvent interactions are probed with slab models, revealing the important influence of surface structure and broken hydrogen bonds. Ethanol strongly interact with the acetyl O atoms near the (001) surface, tuning its properties. Since no hydrogen bond is broken on the (010) surface, the solvation effect is smaller and the distribution of ethanol is more uniform. Depending on the dominant interactions between paracetamol and solvent, electron states of surfaces and near-surface molecules are shifted differently. The density of occupied states increases on the (001) surface, while the density of unoccupied states increases on other surfaces, altering surface reactivities.

The growth of different surfaces is investigated by using two approaches: with implicit or explicit paracetamol molecules. Counterintuitively, neither high accumulations (∼6.18 ρ/ρ_bulk) of paracetamol nor large adsorption energy (∼30.08 kJ mol⁻¹) have been observed for the fast growing (010) surface. Analysis elucidates that competitions among solute molecules might exist near surfaces. In particular, near surfaces with multiple unsaturated hydrogen bonds, paracetamol molecules in solution tend to strongly interact with the higher and more accessible adsorption sites, causing steric hindrances for molecules to bond to the lower and less accessible adsorption sites. Similarly, the solvent, ethanol, might also increase the steric hindrance and screen the adsorption site. This is more prominent when the relatively week NH⋯OH hydrogen bond is broken and ethanol interacts more strongly with the surface. Comparatively, the (010) surface benefits from the absence of unsaturated hydrogen bond and the relatively uniform surface chemistry. A relatively uniform distribution of the adsorbate is observed, while the adsorption energy shows good selectivity that promotes adsorptions in the idealised sequence, minimising steric hindrances among adsorbates. The ‘step 2’ adsorption onto the (010) surface reveals the substantial stabilising effect of hydrogen bond formation and the rate-dominating step.

In summary, this study shows the delicate and complex interplay of thermodynamic and kinetic factors during the growth of molecular crystals. The critical role that broken hydrogen bonds play in interactions with both solvent and solute is highlighted by the comprehensive analysis of the structural and electronic properties of molecule and slab models. Based on that, the microscopic mechanism of the experimentally observed ‘self-healing’ phenomenon is investigated and analysed in detail. With the DFT/ESM-RISM method, growth kinetics can be qualitatively assessed. To quantitively incorporate growth kinetics into morphology prediction, DFT/ESM-RISM based correction terms, such as attachment energy^24,56 and relative occupancy,²⁵ can be proposed in the future to correct surface formation energy. Methods explored here offer a prototype for studying the evolution of crystal morphology, which would be helpful to develop new strategies that control the shapes of crystals during manufacturing.

Author contributions

Huanyu Zhou: conceptualization, methodology, investigation, writing – original draft preparation, visualization. Giuseppe Mallia: supervision, writing – review and editing. Nicholas M. Harrison: supervision, writing – review and editing, funding acquisition, resources.

Data availability

The supplementary data and imagines supporting this article have been included as part of the ESI.†

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work used the ARCHER2 UK National Supercomputing Service (https://www.archer2.ac.uk) via our membership of the UK Car–Parrinello HEC Consortium, which is funded by EPSRC (EP/X035891/1), and Imperial College London's Research Computing Service (https://doi.org/10.14469/hpc/2232). HZ thanks BASF for PhD sponsorship. The authors are grateful for the experimental details shared by Professor Jerry Y. Y. Heng, Dr Isha Bade, Mr Deniz Etit and Ms Allison Arber.

References

1. M. A. Lovette, A. R. Browning, D. W. Griffin, J. P. Sizemore, R. C. Snyder and M. F. Doherty, Ind. Eng. Chem. Res., 2008, 47, 9812–9833.
2. P. Dandekar, Z. B. Kuvadia and M. F. Doherty, Annu. Rev. Mater. Res., 2013, 43, 359–386.
3. S. Datta and D. J. W. Grant, Nat. Rev. Drug Discovery, 2004, 3, 42–57.
4. J. Li and M. F. Doherty, Cryst. Growth Des., 2017, 17, 659–670.
5. Y. Liu, B. Gabriele, R. J. Davey and A. J. Cruz-Cabeza, J. Am. Chem. Soc., 2020, 142, 6682–6689.
6. R. I. Ristic, S. Finnie, D. B. Sheen and J. N. Sherwood, J. Phys. Chem. B, 2001, 105, 9057–9066.
7. G. J. O. Beran, Chem. Rev., 2016, 116, 5567–5613.
8. C. B. Aakeröy and K. R. Seddon, Chem. Soc. Rev., 1993, 22, 397–407.
9. J. Y. Y. Heng, A. Bismarck, A. F. Lee, K. Wilson and D. R. Williams, Langmuir, 2006, 22, 2760–2769.
10. J. Y. Y. Heng and D. R. Williams, Langmuir, 2006, 22, 6905–6909.
11. I. Bade, V. Verma, I. Rosbottom and J. Y. Y. Heng, Mater. Horiz., 2023, 10, 1425–1430.
12. L. Jiang, W. Chen, L. Zhou, L. Xu, F. He, J. Y. Heng, H. Shehzad and J. Ouyang, Sep. Purif. Technol., 2024, 347, 127577.
13. M. Rossi, P. Gasparotto and M. Ceriotti, Phys. Rev. Lett., 2016, 117, 115702.
14. J. L. McKinley and G. J. O. Beran, Faraday Discuss., 2018, 211, 181–207.
15. L. M. LeBlanc, A. Otero-de-la Roza and E. R. Johnson, J. Chem. Theory Comput., 2018, 14, 2265–2276.
16. N. Raimbault, V. Athavale and M. Rossi, Phys. Rev. Mater., 2019, 3, 053605.
17. M. De La Pierre and C. Pouchan, Theor. Chem. Acc., 2018, 137, 25.
18. K. Adhikari, K. M. Flurchick and L. Valenzano, Chem. Phys. Lett., 2015, 630, 44–50.
19. T. Beyer, G. M. Day and S. L. Price, J. Am. Chem. Soc., 2001, 123, 5086–5094.
20. L. Wang, K. C. Lai, L. Huang, J. W. Evans and Y. Han, Surf. Sci., 2020, 695, 121532.
21. X. Tian, T. Wang, L. Fan, Y. Wang, H. Lu and Y. Mu, Appl. Surf. Sci., 2018, 427, 357–362.
22. J. C. M. Silva, H. A. De Abreu and H. A. Duarte, RSC Adv., 2015, 5, 2013–2023.
23. T. Debroise, T. Sedzik, J. Vekeman, Y. Su, C. Bonhomme and F. Tielens, Cryst. Growth Des., 2020, 20, 3807–3815.
24. Q. Zhao, N. Liu, B. Wang and W. Wang, RSC Adv., 2016, 6, 59784–59793.
25. K. Wang, S. Huang and W. Zhu, J. Cryst. Growth, 2022, 585, 126605.
26. J. Rey, P. Clabaut, R. Réocreux, S. N. Steinmann and C. Michel, J. Phys. Chem. C, 2022, 126, 7446–7455.
27. J. Tomasi and M. Persico, Chem. Rev., 1994, 94, 2027–2094.
28. A. Kovalenko and F. Hirata, J. Chem. Phys., 1999, 110, 10095–10112.
29. M. Otani and O. Sugino, Phys. Rev. B: Condens. Matter Mater. Phys., 2006, 73, 115407.
30. S. Nishihara and M. Otani, Phys. Rev. B, 2017, 96, 115429.
31. J. J. Karnes, S. E. Weitzner, S. A. Akhade, S. E. Baker, E. B. Duoss and J. B. Varley, J. Phys. Chem. C, 2022, 126, 12413–12423.
32. Y. Onabuta, M. Kunimoto, S. Wang, Y. Fukunaka, H. Nakai and T. Homma, J. Phys. Chem. C, 2022, 126, 5224–5232.
33. H. Yu, S. E. Weitzner, J. B. Varley, B. C. Wood and S. A. Akhade, J. Phys. Chem. C, 2023, 127, 1789–1797.
34. J. Y. Y. Heng, F. Thielmann and D. R. Williams, Pharm. Res., 2006, 23, 1918–1927.
35. C. Boukouvala, J. Daniel and E. Ringe, Nano Converg., 2021, 8, 26.
36. P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. D. Corso, S. d Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari and R. M. Wentzcovitch, J. Phys.: Condens. Matter, 2009, 21, 395502.
37. P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. B. Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, N. Colonna, I. Carnimeo, A. D. Corso, S. d Gironcoli, P. Delugas, R. A. DiStasio, A. Ferretti, A. Floris, G. Fratesi, G. Fugallo, R. Gebauer, U. Gerstmann, F. Giustino, T. Gorni, J. Jia, M. Kawamura, H.-Y. Ko, A. Kokalj, E. Küçükbenli, M. Lazzeri, M. Marsili, N. Marzari, F. Mauri, N. L. Nguyen, H.-V. Nguyen, A. Otero-de-la Roza, L. Paulatto, S. Poncé, D. Rocca, R. Sabatini, B. Santra, M. Schlipf, A. P. Seitsonen, A. Smogunov, I. Timrov, T. Thonhauser, P. Umari, N. Vast, X. Wu and S. Baroni, J. Phys.: Condens. Matter, 2017, 29, 465901.
38. A. Otero-de-la Roza and E. R. Johnson, J. Chem. Phys., 2012, 137, 054103.
39. A. J. A. Price, R. A. Mayo, A. Otero-de-la Roza and E. R. Johnson, CrystEngComm, 2023, 25, 953–960.
40. A. J. A. Price, A. Otero-de-la Roza and E. R. Johnson, Chem. Sci., 2023, 14, 1252–1262.
41. W. L. Jorgensen, D. S. Maxwell and J. Tirado-Rives, J. Am. Chem. Soc., 1996, 118, 11225–11236.
42. S. Nishihara, nisihara1/MOLs, 2023, https://github.com/nisihara1/MOLs.
43. W. L. Jorgensen and J. Tirado-Rives, J. Am. Chem. Soc., 1988, 110, 1657–1666.
44. I. U. of Crystallography, Symbols for space groups, International Tables for Crystallography Volume A: Space-Group Symmetry, Springer, Dordrecht, Netherlands, 5th edn, 2006, ch. 4.3, pp. 62–76.
45. C. Wilson, Z. Kristallogr. - Cryst. Mater., 2000, 215, 693–701.
46. S. Fleming and A. Rohl, Z. Kristallogr. - Cryst. Mater., 2005, 220, 580–584.
47. K. Momma and F. Izumi, J. Appl. Crystallogr., 2011, 44, 1272–1276.
48. A. Kokalj, J. Mol. Graphics Modell., 1999, 17, 176–179.
49. J. Junquera, M. H. Cohen and K. M. Rabe, J. Phys.: Condens. Matter, 2007, 19, 213203.
50. D. Roy and A. Kovalenko, J. Phys. Chem. A, 2019, 123, 4087–4093.
51. D. Roy and A. Kovalenko, J. Phys. Chem. B, 2020, 124, 4590–4597.
52. D. K. Owens and R. C. Wendt, J. Appl. Polym. Sci., 1969, 13, 1741–1747.
53. J. Panzer, J. Colloid Interface Sci., 1973, 44, 142–161.
54. S. Schönecker, X. Li, B. Johansson, S. K. Kwon and L. Vitos, Sci. Rep., 2015, 5, 14860.
55. A. Forslund, X. Zhang, B. Grabowski, A. V. Shapeev and A. V. Ruban, Phys. Rev. B, 2021, 103, 195428.
56. Y. Liu, S. Niu, W. Lai, T. Yu, Y. Ma, H. Gao, F. Zhao and Z. Ge, CrystEngComm, 2019, 21, 4910–4917.
57. T. J. Frankcombe and O. M. Løvvik, J. Phys. Chem. B, 2006, 110, 622–630.
58. H. Sato, A. Kovalenko and F. Hirata, J. Chem. Phys., 2000, 112, 9463–9468.
59. M. A. Bellucci, G. Gobbo, T. K. Wijethunga, G. Ciccotti and B. L. Trout, J. Chem. Phys., 2019, 150, 094107.
60. N. Kovačević and A. Kokalj, J. Phys. Chem. C, 2011, 115, 24189–24197.
61. P. W. Carter, A. C. Hillier and M. D. Ward, J. Am. Chem. Soc., 1994, 116, 944–953.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d5mh00170f

---