Study of the solvent-dependent crystal shape of theophylline using constant chemical potential molecular dynamics simulations

Abstract

The crystal habit of an organic crystalline material impacts its properties, processing, and performance, especially in pharmaceuticals. In solution crystallization, solvents play a crucial role in modulating the crystal habits by interacting with the different growing facets – either enhancing or inhibiting the growth of specific facets. Thus, an in-depth understanding of the role of the solvent in crystal shape selection is of paramount importance for the design and growth of crystals. In this work, we used constant chemical potential molecular dynamics (CμMD) to simulate the growth of theophylline crystals in solvents with decreasing polarity, i.e. water, isopropyl alcohol, and dimethylformamide. Our results indicate that as the polarity of the solvent increases, theophylline crystallizes into rod-shaped crystals; the aspect ratio is modulated by the relative growth of the (001) and (010) facets. Furthermore, solvent profile analyses revealed that the desolvation of the crystal facets plays a major role in the growth process.

---

1 Introduction

The large diversity observed in crystal habits,¹ ranging from the simple cubic to hexagonal prism shape of snowflakes,² is due to the complex interplay of various microscopic events such as local density fluctuations, molecular ordering, nucleation and dissolution that take place during the crystallization process. Understanding the phenomenon of crystal growth is not only important in academics but also critical in various industrial processes, especially in pharmaceuticals.^3–5 Most of the active pharmaceutical ingredients (APIs) in pharmaceutical industries are manufactured via a solution crystallization route. The resulting crystal shape impacts various properties like filterability, tabletability, wettability, and flow properties during various production processes. Afterward, it can also impact the physicochemical properties of the drug like bioavailability.^6–8 Therefore, controlling crystal habits during manufacturing processes is of critical importance to improve the properties of the material for desired applications.

The notion of equilibrium crystal shape has been given by Gibbs⁹ and is made practical by the Wulff construction. According to Gibbs and Wulff, the shape of the crystal is determined by thermodynamic factors, and the equilibrium shape is the one that minimizes the surface energy of the crystal for a given volume. But in reality crystal growth is determined by kinetic factors also; especially in the case of solution crystallization, the morphology dictated by kinetic factors is known as growth morphology.¹⁰ The growth morphology depends on the relative growth rates of different crystal facets in a given solvent. It is well-entrenched in the literature that the aspect ratio of crystals depends on the relative polarity of the solvent.¹¹ The steady-state habit of a crystal in the solution can be predicted by studying the solid-state structure of molecules and the solvent effects.^12–14

Experiments can be used to determine the shape of the resulting crystal in solution; however, they are often limited to unveiling the crystal growth mechanism and the role of a solvent in mediating the preferential growth of certain facets. In the literature, there are well-known computational methods based on the attachment energy (AE) model,¹⁵ which has been used extensively to predict crystal morphology.^16,17 According to the AE model, the relative growth rate of a face always increases with the increase in AE. Although the functional form of the relative growth rate depends on the growth mechanism and various factors such as supersaturation, temperature, and solid–fluid interaction, the AE method provides the crystal morphology but cannot provide the crystal growth mechanism. Molecular dynamics (MD) simulations decipher crystal growth at the molecular level, offering insights into different crystal shapes and even assisting in achieving desired crystal shapes by understanding the underlying reasons.^18–21 Salvalaglio et al.²² used MD simulations to understand the growth of urea crystal from aqueous solution. They unveiled the growth mechanism of the (001) and (110) facets of urea and predicted the shape of urea crystals in different solvents by calculating the relative growth rates between different facets.²³ Bjelobrk et al.²⁴ predicted the shape of naphthalene crystal in ethanol using constant chemical potential molecular dynamics (CμMD) simulations in combination with well-tempered metadynamics (WTMetaD) techniques. In 2019, Han et al.²⁵ predicted the solvent-dependent morphology of an anti-tuberculosis drug, isoniazid. These works manifested the effective use of CμMD simulations in carrying out crystallization simulations under constant supersaturation conditions and opened the possibility of studying crystallization of APIs under realistic experiment-like conditions.

In this work, we studied the solvent-dependent shape evolution of theophylline using CμMD simulations. Theophylline is a pharmaceutically active compound used to treat asthma and cardiopulmonary diseases. It is known to crystallize into lath- or rod-shaped crystals in polar solvents such as alcohols and solvents having H-bond donor and acceptor groups, and triangular-shaped crystals in the case of non-polar or H-bond acceptor solvents like toluene and ethyl acetate.²⁶ We have selected three solvents with their relative polarity (R)²⁷ (a scale that refers to how polar a solvent is with respect to other solvents) given in parentheses: water (R = 1.0), DMF (R = 0.4) and IPA (R = 0.5) to study the solvent-dependent shape of theophylline. Furthermore, the role of solvents in crystal growth processes has been deciphered by analyzing the desolvation of various crystal facets.

2 Computational methods

2.1 Constant chemical potential molecular dynamics (CμMD)

Unlike an experimental setup, standard MD simulations cannot maintain the condition of constant supersaturation and suffer from finite size effects.^23,28 Therefore, the results cannot be directly compared with the experiments. The method developed by Perego et al.²⁹ commonly known as CμMD takes care of the finite size effect and maintains a constant supersaturation in the region adjacent to the growing crystal surface.

In a CμMD simulation, a crystal slab is placed in the middle of an orthorhombic box filled with the solution on both sides of the solid surface. On each side, we define a control region (CR) from the solid–solution interface at which the solution concentration is controlled at a target value using an external force (force region, FR) applied to solutes belonging to the CR.²⁹ The region between the solid–solution interface (Z_I) and the CR is called the transition region (TR). The remaining volume beyond the FR acts as a reservoir that supplies molecules to the CR and thereby helps in maintaining the constant supersaturation condition. The schematic of CμMD is shown in Fig. 1, where n_R and n_B are the number density of solute in the CR and the buffer region (bulk), respectively.

Fig. 1 Schematic of CμMD. Theophylline molecules belonging to the crystal are shown in blue color, those belonging to the solution are shown in red colour, and the solvent molecules are shown in grey colour. Z_I is the interface between the crystal and solution, ε is the width of TR and Z_I + D_F represents the position of the force region. F^μ is the force that acts on the molecules that enter the FR.

We start by calculating on-the-fly the solid–solution interface, Z_I, from which the TR, CR, and FR are defined. The solution concentration (n^CR_i, where the index i refers to solute or solvent species) is defined as the number of solutes (N_i) per unit volume of the CR (V^CR) using eqn (1):

where
where z_in and z_out are the inner and outer CR boundaries, respectively, and α is a parameter that controls the steepness of the switching functions.

The force that restraints the instantaneous solution concentration n^CR_i at a target value (n_i0) is given by

F^μ_i(z) = k(n^CR_i − n_i0)G(z,Z_F) (2)

where k is the force constant and G(z,Z_F) is a bell-shaped function that is non-zero in the vicinity of the force center Z_F = Z_I + D_F. The force defined above acts on the solute molecules to compensate for the deviations in the instantaneous CR density n^CR_i from the target value, n_i0. The function G(z,Z_F) localizing the function action of force is defined aswhere ω is the width of the function. The parameters used in CμMD simulations in the case of all three solvents are listed in Table 1 and the details related to choosing TR, CR and FR are provided in the ESI† (Section 1). During all simulations, a constant number density of the solute is maintained in the CR that corresponds to the supersaturation (S) range between 4 and 7.

Table 1 The parameters used for all CμMD simulations

System	ε (nm)	CR (nm)	D F (nm)	σ F (nm)	k (nm^3 kJ mol^−1)
Water (001)	2.0	3.5	5.7	0.5	50 000
Water (010)	2.7	4.0	6.9	0.5	80 000
DMF (001)	2.0	3.0	5.2	0.5	50 000
DMF (010)	2.0	3.0	5.2	0.5	80 000
IPA (001)	2.0	3.0	5.7	0.5	50 000
IPA (010)	2.0	3.0	5.7	0.5	80 000

---

In our simulations, the supersaturation (S) is defined in the following way. At first, the molar concentration is calculated from the formula:

where C_s and C_l are the concentrations of the solute and the solvent molecules in the CR, respectively. The solution supersaturation (S) is then defined aswhere x₀ is the experimental solubility value.

2.2 Definition of the number of crystalline molecules, N

The number of crystalline molecules as a function of time is calculated throughout the trajectory to provide a quantitative description of the crystallization process. A function defined by Salvalaglio et al.²² is used to calculate the number of crystalline molecules, which can distinguish between molecules in the crystal and solution.³⁰ In the crystalline state, the ith molecule has a well-defined orientation and a fixed local density, while in the liquid state, molecules are disordered. A function called the degree of crystallinity is defined by Γ_i(t) that measures the product of two functions describing the local density ρ(n_i(t)) and local order around a given molecule. If a molecule has at least as many neighbors within a cutoff as a molecule present in the interface has and its orientation matches with the crystal angles, it is assigned a local density value of 1; otherwise, it is assigned a value of less than 1. The expressions of the functions are as follows:where n_i(t) is the number of neighbours of the ith molecule, σ_n is the standard deviation in the solid state, and is a vector that represents mutual orientations θ_i,j between the vectors associated with the ith and jth molecules. Since the molecules in the crystalline solid are not completely static, their orientation with respect to each other fluctuates over time and σ_θ expresses the extent of these fluctuations, whereas [small theta, Greek, macron]₀ and [small theta, Greek, macron]₁ are the mean values corresponding to the mutual orientations of the internal vector defined for the ith molecule with respect to its neighbors in the crystalline state. Once Γ_i(t) values are computed for all the theophylline molecules in the system, the number of molecules belonging to the crystalline phase at each time step t can then be obtained as

A summary of the parameters used in the calculation of N_c(t) is given in Table 2 and the ESI† (Fig. S2).

Table 2 Parameters used for the calculation of N_c(t)

Parameter	Value	Parameter	Value
r 0	0.5	n 0	3
θ 1	0.10	σ 1	0.20
θ 2	0.60	σ 2	0.20
m	12	n	24

---

2.3 Simulation details

All simulations were performed using the Gromacs software package version 2021.4^31–34 patched with a private version of Plumed-2.8.0.^35,36 The CHARMM general force field (CGenFF)³⁷ with full atomistic description was used to model theophylline, IPA, and DMF used in this work. All simulations were performed at a constant temperature and pressure of 350 K and 1 bar, respectively. The steepest descent algorithm was used to minimize the initial configurations. A cut-off of 0.7 nm was used for both short range Coulomb and van der Waals interactions. The particle mesh Ewald algorithm was used to calculate long-range electrostatics. The velocity rescale thermostat³⁸ was used to keep the temperature at 350 K during the simulations.³⁹ The Berendsen barostat⁴⁰ was used to keep the pressure constant during NPT equilibration, and the Parrinello Rehman barostat⁴¹ was used to keep the pressure constant at 1 bar in all NPT production runs. The covalent bonds connected to hydrogen atoms were constrained by the LINCS method. An integration time step of 0.5 fs was used during NPT equilibration and 2 fs was used in all production run simulations.

2.4 System preparation

The crystallographic structure of theophylline was obtained from the CCDC database.⁴² It is found that theophylline belongs to an orthorhombic crystal system with the space group Pna2₁ and unit cell dimensions of a = 24.61 Å, b = 3.83 Å, and c = 8.50 Å. We considered the growth of the two facets (001) and (010) of theophylline that predominantly determine the crystal morphology. Initial coordinates required for the simulation were taken from the experimentally determined crystal structure. A 2 × 2 × 2 supercell was constructed from the experimentally determined crystal structure using the AVOGADRO software⁴³ and then subjected to an NPT equilibrium simulation at a temperature of 350 K and a pressure of 1 bar to check its stability. The insignificant deviation of the equilibrated lattice parameters from that of the experimental ones indicates the robustness of the force field used in our simulations. The genconf utility of GROMACS was used to prepare larger crystal structures for the simulations.

The generated crystal slabs that were used to study (001) and (010) surface growth have an approximate size of 4.72 × 1.47 × 6.53 nm³ and 4.72 × 2.94 × 3.26 nm³, respectively. Each slab has 256 theophylline molecules. In the subsequent step, the crystal slabs generated for (001) and (010) surfaces were placed in the center of the simulation box with these facets exposed to the solution along the z-axis and having dimensions of 4.72 × 1.47 × 25.00 nm³ and 4.72 × 3.26 × 20.00 nm³, respectively, except for the water (010) facet having dimensions of 4.72 × 3.26 × 25.00 nm³. More details related to the preparation of the systems are provided in the ESI† and in Table 3. We have performed two sets of simulations for each system mentioned in Table 1. We have obtained similar results for both sets of simulations. The results of one set of simulations are discussed in this manuscript (Fig. 3 and 5), while the results of the other set of simulations are discussed in the ESI† (Fig. S6 and S7).

Table 3 Details of the systems simulated with the value of solution supersaturation^a

Surface	System	Solute	Solvent	C s (nm^−3)	S
a According to Joswiak et al.^44 and Zimmermann et al.,^45 the simulated solubility should be used to calculate the supersaturation values. However, due to the lack of simulated solubility data, the experimental solubility^46 observed at 300 K has been considered in our simulations. The experimental solubility of THE in IPA is not known.
0 0 1	Water	127	4382	0.1	4.72
	IPA	127	944	0.3	—
	DMF	127	899	0.3	6.76
0 1 0	Water	208	11 948	0.1	4.03
	IPA	208	1850	0.3	—
	DMF	258	2159	0.3	6.11

---

3 Results and discussion

In 2015, Mark et al. studied the solvent-dependent shape of theophylline crystal using transmission electron microscopy.²⁶ They found that theophylline crystallizes into unusual triangular plate-like crystals when crystallized in solvents that possess hydrogen bond acceptor groups, but no donor groups like nitromethane, ethyl acetate, and dioxane. Additionally, triangular crystals were also obtained using highly non-polar solvents such as toluene. On the other hand, solvents with both hydrogen bond donor and acceptor groups or polar solvents, such as methanol, ethanol, and IPA, yield rod- or lath-shaped crystals. They have explained the crystal habit by the relative growth rate of (001) and (010) facets. Based on this study, we selected these two facets to investigate the solvent-dependent crystal habit of theophylline in solvents of different polarities. The CμMD method has been used to study theophylline growth in water, DMF, and IPA with the desired solute concentration (number density) in the CR (shown in Section 2, ESI†).

3.1 Crystal growth of the (001) surface

Let us focus on the crystal growth simulations along the 〈001〉 direction in water, DMF, and IPA. Along this direction, the (001) and (001̄) facets of theophylline are exposed to the solution. It is observed that these crystal facets are not equivalent at the molecular level as hydrogen bond donating –NH groups are exposed in the (001) facet, whereas hydrogen bond accepting nitrogen atoms (possessing a lone pair of electrons) are exposed in (001̄) as shown in Fig. 2A.

Fig. 2 (A) The N (blue arrows) and NH (yellow arrows) groups exposed to the solution in (001̄) and (001) crystal facets, respectively. (B) (01̄0) and (010) crystal facets. The interfacial layers of both surfaces are well decorated with kink (A) and step sites (B) (shown by red arrows).

We carried out 1 μs of simulation of each system and monitored the crystal growth. In Fig. 3, the number of growing crystalline molecules (N) is plotted as a function of simulation time. In the case of IPA, one layer of theophylline grew on the (001) surface within ∼300 ns. The other, i.e. (001̄) facet, starts growing after the growth of one layer along the (001) facet. On the other hand, in the case of DMF, one partially grown layer on both sides of the crystal slab was observed. In the case of water, theophylline grows in a similar way to that in the case of IPA but the growth is less as compared to other solvents.

Fig. 3 Growth profiles of (A) (001) and (B) (010) surfaces in water (w), DMF (d), and IPA (i) obtained from simulations; (C) and (D) the growth of theophylline in IPA along 〈001〉 and in water along 〈010〉, respectively. Snapshots were taken at ∼1 μs. Theophylline molecules belonging to the crystal are shown in blue and those belonging to the solution are shown in red, while the IPA and water molecules are shown in ice blue and grey colour, respectively.

3.2 Crystal growth of the (010) surface

Now, we switch our attention to the crystal growth along the 〈010〉 direction. Along this direction, the (010) and (01̄0) facets of theophylline are exposed to the solution. On the molecular level, these two facets are similar and have non-polar methyl groups exposed to the solution. Both the facets (010) and (01̄0) grow to an almost equal extent in all three solvents. The maximum growth that we observe along (01̄0) is for water, and that is approximately one layer of growth on both facets in the simulation time of 1 μs (Fig. 3). In the case of DMF and IPA, the growth is relatively less than one layer on both sides as compared to water.

3.3 Details of the growth mechanism

Crystal growth theories are based on the consideration of the crystal surface structure. One of the most accepted models of crystal growth was given by Kossel.^47,48 This model envisions the crystal structure as being made of cubical units and explains layer, and birth and spread crystal growth mechanisms. The event of crystal growth consists of a series of steps through which the atoms or molecules get incorporated into the crystal. These processes can be summarised in four steps: (1) transport of atoms through the solution, (2) attachment of atoms to the surface, (3) movement of atoms on the surface, and (4) attachment of atoms to edges or kinks. Step (1) is called the transport process, and steps (2)–(4) are the surface processes. These processes naturally occur in a series, and the overall growth rate depends upon the slowest step of these processes. So the growth can be either surface-controlled or transport-controlled depending upon the slowest step.^14,49

3.4 Number of solvent and solute molecules near the (001) surface interface

To understand the growth mechanism, we have plotted the number of solute and solvent molecules (N) within 0.5 nm near the interface for both crystal surfaces. First, let us discuss the growth along the 〈001〉 axis along which (001) and (001̄) facets are exposed to the solution. These facets are not smooth and are already well decorated with kink and step sites as can be seen from Fig. 2A and B. From the visualization of the simulation trajectories in all three solvents, we observe that when the solute molecules come close to the interface they quickly start occupying these sites one by one. When all the sites are occupied, the same type of step or kink sites are generated in the new interfacial layer.

The number of solute molecules around the interface for the (001) surface depends on the nature of the solvent. The value of N is the highest for IPA compared to DMF and water, as shown in Fig. 4A. The highest growth in the case of IPA can be validated by the least number of solvents around the interface (Fig. 4B). One of the possible reasons for the highest growth in the case of IPA can be its weak interactions with the polar (001) surface as can be seen from Fig. 5(c) due to which solute molecules can easily access the interface. In the case of water being most polar among all the solvents, it shows stronger interaction with the polar interface as it forms extended H-bonded chains with the nitrogen and carbonyl of the interfacial theophylline due to which the growth of the solute is hindered as compared to IPA (Fig. 5(a)). One unexpected result that we observed was the case of DMF which being the least polar should show the least interaction with the polar interface and should show the highest growth among all the solvents. One reason might be the slow desolvation of DMF (can be seen from the slope of DMF in Fig. 4B) from the kink site which might compensate for the polarity, and it can be the rate-limiting step in the case of DMF due to which the solute grows less than IPA. Fig. 4B of IPA shows the fastest and maximum decrease, indicating the highest growth rate in IPA. The order of growth is IPA > DMF ≈ water.

Fig. 4 TThe number of solute (A and C) and solvent (B and D) molecules (N) within a 0.5 nm region around the (001) (A and B) and (010) (C and D) crystal interfaces as a function of simulation time in water, DMF, and IPA simulation systems.

Fig. 5 Snapshots taken during simulation showing the interaction of interfacial (001) theophylline with (a) water, (b) DMF and (c) IPA.

3.5 Number of solvent and solute molecules near the (010) surface interface

From visual inspection, we observed that the (010) surface is decorated with step sites as can be seen from Fig. 2B. From the simulation trajectories, we have observed that when solute molecules come closer to the surface they randomly get attached to the interfacial step sites. The growth on this surface is controlled by diffusion, which is a feature of the rough growth mechanism. In the case of water, the maximum number of solutes get attached to the surface compared to DMF and IPA (Fig. 4C). This might be due to the polarity of the solvent. As the surface is non-polar, more polar solvents like water do not interact much with the non-polar surface, and hence the crystal grows faster, and less polar and bulky solvents like IPA and DMF hinder solute attachment to the growing surface due to non-polar type interactions. Therefore, for the (010) surface the order of growth is water > DMF ≈ IPA.

3.6 Calculation of the growth rate

The growth rate of theophylline in all solvents is calculated from simulation trajectories. The growth rate was calculated from the slope of the N vs. time plot shown in Fig. 4A and C. The slope and intercept of the plots along with the aspect ratio are given in Table 4.

Table 4 Slope, intercept and aspect ratio (ratio of the slopes of 010 and 100) of the plots shown in Fig. 4

System	Slope	Intercept	Aspect ratio
Water (001)	12.79	265.38	2.59
Water (010)	33.17	311.53
DMF (001)	24.05	267.53	1.10
DMF (010)	26.63	278.69
IPA (001)	55.73	264.18	0.49
IPA (010)	27.60	273.27

---

4 Conclusion

In this work, we have investigated the solvent-dependent crystal growth of theophylline and provided a rationale for the experimentally observed crystal shapes obtained from solution crystallization using solvents like water, DMF, and IPA. Experimentally, it has been observed that theophylline crystallizes into lath- or rod-shaped crystals in the case of polar solvents having H-bond donor or acceptor groups. From the long timescale CμMD simulations, we found that the growth rate of the (010) facet is higher than (001) with the exception of IPA. In the case of the (010) facet, the relative growth rate decreases from water through DMF to IPA, while for the surface (001), the relative growth rate increases from water through DMF to IPA. We concluded that the aspect ratio of crystal decreases from water through DMF to IPA. Unfortunately, the unavailability of experimental growth rates does not allow us to make a direct comparison between the experimental and computational results. Our study helps understand the solvent-dependent crystal growth of theophylline and delineates the role of the solvent in the growth mechanism. This work further epitomizes how one can choose the right set of solvents in pharmaceutical crystallization to synthesize the crystals of the desired shape with computer assistance.

Data availability

Input files for simulation and the code for the calculation of the probability distribution of angle can be found at https://github.com/Neha-Mehta98.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

Neha thanks IIT Delhi for the institute PhD fellowship. T. K. acknowledges the Science and Engineering Research Board (SERB), New Delhi, India for the Start-up Research Grant (File No. SRG/2022/000969). We also acknowledge IIT Delhi for the Seed Grant and the institute's HPC facility for computational resources.

References

1. I. Sunagawa, Morphology of Crystals, Parts A and B, 1989.
2. W. A. Bentley and W. J. Humphreys, Snow crystals, Courier Corporation, 2013.
3. J. Swarbrick, Encyclopedia of Pharmaceutical Technology: Volume 6, CRC press, 2013.
4. A. Bongioanni, M. S. Bueno, B. A. Mezzano, M. R. Longhi and C. Garnero, Crystal Growth and Chirality-Technologies and Applications, IntechOpen, 2022.
5. L. S. Taylor, D. E. Braun and J. W. Steed, Cryst. Growth Des., 2021, 21(3), 1375–1377.
6. N. Variankaval, A. S. Cote and M. F. Doherty, AIChE J., 2008, 54, 1682–1688.
7. S. R. Modi, A. K. Dantuluri, V. Puri, Y. B. Pawar, P. Nandekar, A. T. Sangamwar, S. R. Perumalla, C. C. Sun and A. K. Bansal, Cryst. Growth Des., 2013, 13, 2824–2832.
8. P. Bukovec, P. Benkič, M. Smrkolj and F. Vrečzer, Die Pharmazie-An Int. J. Pharm. Sci., 2016, 71, 263–268.
9. J. W. Gibbs, Thermodynamics, Longmans, Green and Company, 1906, vol. 1.
10. A. Chernov, Crystallogr. Rep., 1963, 7, 728–730.
11. Y. Wang and Z. Liang, CrystEngComm, 2017, 19, 3198–3205.
12. X. Liu, E. Boek, W. Briels and P. Bennema, Nature, 1995, 374, 342–345.
13. D. Winn and M. F. Doherty, AIChE J., 1998, 44, 2501–2514.
14. P. Cubillas and M. W. Anderson, Zeolites and catalysis: synthesis, reactions and applications, 2010, pp. 1–55.
15. P. Hartman and P. Bennema, J. Cryst. Grow., 1980, 49, 145–156.
16. T. Beyer, G. M. Day and S. L. Price, J. Am. Chem. Soc., 2001, 123, 5086–5094.
17. R. A. Sullivan and R. J. Davey, CrystEngComm, 2015, 17, 1015–1023.
18. D. Winn and M. F. Doherty, Chem. Eng. Sci., 2002, 57, 1805–1813.
19. S. Piana and J. D. Gale, J. Am. Chem. Soc., 2005, 127, 1975–1982.
20. S. Piana, F. Jones and J. D. Gale, J. Am. Chem. Soc., 2006, 128, 13568–13574.
21. J. Chen and B. L. Trout, Cryst. Growth Des., 2010, 10, 4379–4388.
22. M. Salvalaglio, T. Vetter, F. Giberti, M. Mazzotti and M. Parrinello, J. Am. Chem. Soc., 2012, 134, 17221–17233.
23. M. Salvalaglio, T. Vetter, M. Mazzotti and M. Parrinello, Angew. Chem., 2013, 125, 13611–13614.
24. Z. Bjelobrk, P. M. Piaggi, T. Weber, T. Karmakar, M. Mazzotti and M. Parrinello, CrystEngComm, 2019, 21, 3280–3288.
25. D. Han, T. Karmakar, Z. Bjelobrk, J. Gong and M. Parrinello, Chem. Eng. Sci., 2019, 204, 320–328.
26. M. D. Eddleston, K. E. Hejczyk, A. M. Cassidy, H. P. Thompson, G. M. Day and W. Jones, Cryst. Growth Des., 2015, 15, 2514–2523.
27. C. Reichardt, Chem. Rev., 1994, 94, 2319–2358.
28. T. Karmakar, P. M. Piaggi, C. Perego and M. Parrinello, J. Chem. Theory Comput., 2018, 14, 2678–2683.
29. C. Perego, M. Salvalaglio and M. Parrinello, J. Chem. Phys., 2015, 142, 144113.
30. C. Micheletti, A. Laio and M. Parrinello, Phys. Rev. Lett., 2004, 92, 170601.
31. B. Hess, C. Kutzner, D. Van Der Spoel and E. Lindahl, J. Chem. Theory Comput., 2008, 4, 435–447.
32. H. J. Berendsen, D. van der Spoel and R. van Drunen, Comput. Phys. Commun., 1995, 91, 43–56.
33. D. Van Der Spoel, E. Lindahl, B. Hess, G. Groenhof, A. E. Mark and H. J. Berendsen, J. Comput. Chem., 2005, 26, 1701–1718.
34. S. Pronk, S. Páll, R. Schulz, P. Larsson, P. Bjelkmar, R. Apostolov, M. R. Shirts, J. C. Smith, P. M. Kasson and D. Van Der Spoel, et al. , Bioinformatics, 2013, 29, 845–854.
35. G. A. Tribello, M. Bonomi, D. Branduardi, C. Camilloni and G. Bussi, Comput. Phys. Commun., 2014, 185, 604–613.
36. Nat. Methods, 2019, 16, 670–673.
37. K. Vanommeslaeghe and A. D. MacKerell Jr, J. Chem. Inf. Model., 2012, 52, 3144–3154.
38. G. Bussi, D. Donadio and M. Parrinello, J. Chem. Phys., 2007, 126, 014101.
39. B. Hess, H. Bekker, H. J. Berendsen and J. G. Fraaije, J. Comput. Chem., 1997, 18, 1463–1472.
40. H. J. Berendsen, J. v Postma, W. F. Van Gunsteren, A. DiNola and J. R. Haak, J. Chem. Phys., 1984, 81, 3684–3690.
41. M. Parrinello and A. Rahman, J. Appl. Phys., 1981, 52, 7182–7190.
42. A. A. Naqvi and G. C. Bhattacharyya, J. Appl. Crystallogr., 1981, 14, 464.
43. M. D. Hanwell, D. E. Curtis, D. C. Lonie, T. Vandermeersch, E. Zurek and G. R. Hutchison, J. Cheminf., 2012, 4, 1–17.
44. M. N. Joswiak, B. Peters and M. F. Doherty, Cryst. Growth Des., 2018, 18, 6302–6306.
45. N. Zimmermann, B. Vorselaars, J. R. Espinosa, D. Quigley, W. R. Smith, E. Sanz, C. Vega and B. Peters, et al. , J. Chem. Phys., 2018, 148, 222838.
46. P. Cysewski, T. Jeliński, P. Cymerman and M. Przybyłek, Int. J. Mol. Sci., 2021, 22, 7347.
47. V.-W. Kossel, Ann. Phys., 1934, 413, 457–480.
48. W. R. Wilcox, Prog. Cryst. Growth Charact. Mater., 1993, 26, 153–194.
49. 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.

---

Footnote

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

---