Naphthalene crystal shape prediction from molecular dynamics simulations

Abstract

We used molecular dynamics simulations to predict the steady state crystal shape of naphthalene grown from ethanol solution. The simulations were performed at constant supersaturation by utilizing a recently proposed algorithm [Perego et al., J. Chem. Phys., 2015, 142, 144113]. To bring the crystal growth within the timescale of a molecular dynamics simulation we applied well-tempered metadynamics with a spatially constrained collective variable, which focuses the sampling on the growing layer. We estimated that the resulting steady state crystal shape corresponds to a rhombic prism, which is in line with experiments. Further, we observed that at the investigated supersaturations, the {001̄} face grows in a two step two dimensional nucleation mechanism while the considerably faster growing faces {11̄0} and {201̄} grow new layers with a one step two dimensional nucleation mechanism.

---

1 Introduction

Controlling the crystal shape of organic molecules is crucial for process optimization in a wide range of industries, in particular in the pharmaceutical industry. Active pharmaceutical ingredients (APIs) are commonly produced as crystal powders purified by crystallization from solution in stirred suspension crystallizers. The resulting crystal shape affects among other properties the filterability, tabletability, and bio-availability of the drug product. Hence, control of the crystal shape is crucial in the manufacturing of pharmaceutical compounds.^1,2 Experimental methods aimed at determining the crystal shape are often expensive and cumbersome. For this reason optimizing the crystallization conditions to obtain the desired crystal shape can be a laborious trial-and-error endeavor. Obtaining insight into the growth process at the molecular level would be useful to make a rational choice for the crystallization conditions.

Molecular dynamics (MD) simulations have recently emerged as a useful tool to obtain atomistic insight into crystal growth.^3–10 In particular, the work of Salvalaglio et al.¹¹ showed the usefulness of MD simulations for the paradigmatic case of urea. In that work, the growth rate was calculated and the growth mechanisms, rough growth or two dimensional (2D) nucleation, were identified from simulations in which the crystal surface of interest was exposed to a supersaturated solution. Afterwards, the crystal shape was predicted based on the relative growth rates of the different crystal faces.¹²

The present work is dedicated to the study of naphthalene crystal growth in ethanol solution under different supersaturations.^13,14 After urea, naphthalene is a step further towards the simulation of systems as complex as real life APIs.

The study of crystal growth of naphthalene using MD simulations presents at least two major challenges. First, the growth rates of the slowest growing crystal faces of naphthalene are so small, that growth of new layers is not observable within the timescale of a standard MD simulation (∼1 μs). To overcome the simulation time limitation we employ well-tempered metadynamics¹⁵ (WTMetaD), which is a well established enhanced sampling method for the study of rare events.

Secondly, the strong finite size effects hamper direct comparison with experiments. For instance, typical experiments are performed at approximately constant supersaturation while standard MD simulations are unable to maintain such conditions.^8,11 The finite size effects are addressed using the constant chemical potential molecular dynamics (CμMD) method developed by Perego et al.,¹⁶ which maintains a constant supersaturation in the region adjacent to the crystal. With this simulation setup, relative growth rates of faces growing in the 2D nucleation regime can be obtained as a function of supersaturation.

2 Methods

2.1 Chemical potential control

In this section, we briefly summarize the CμMD method.¹⁶ The main features of this method are schematically depicted in Fig. 1: the system has periodic boundary conditions in all directions and is segmented along the z-axis into one reservoir, two control regions, as well as transition regions which surround the crystal slab. The slab is oriented such that the face of interest is exposed to the solution along the z-axis.

Fig. 1 Schematic of the CμMD method illustrating the segmentation of the simulation box. The two crystal surfaces are perpendicular to the z-axis and exposed to the solution at positions −z_S, and z_S respectively. The external force positions are marked at −z_F and z_F.

The method works as follows: two forces, F^μ_L and F^μ_R, are introduced in order to control the flux of molecules between the reservoir and the left and right control regions, respectively. Each force acts along the z-axis at a distance D_F from the crystal surface position z_S, namely at positions −z_F and z_F, where z_F = z_S + D_F. The functional forms of the forces are

F^μ_L(t,z) = k(C_CR,L(t) − C_0,L)G_ω(z, −z_F) (1)

F^μ_R(t,z) = k(C_CR,R(t) − C_0,R)G_ω(z, −z_F), (2)

where t is the time, C_CR,L(t) and C_CR,R(t) denote the instantaneous solute concentrations in the left and right control regions, while C_0,L and C_0,R are the predefined set-values of these concentrations. k is a force constant, and G_ω(z,z_F) corresponds to the bell-shaped function defined aswhose height and width are controlled by the parameter ω.

The idea behind the method is that depending on the concentration in the control regions at time t, F^μ will accelerate the molecules towards or away from the control regions. In this way, the concentration in the control regions can be kept constant. Further details on the CμMD method and the parameter optimization used in this work are reported in ref. 16.

2.2 Metadynamics

The crystal growth rate depends on the growing crystallographic plane and the solution supersaturation. In the case of naphthalene, even at high supersaturations, the growth of certain faces takes place in timescales much longer than those that can be afforded using standard MD simulations. A variety of different enhanced sampling techniques exist which can be applied to bring the growth event within the timescale of the simulation, like umbrella sampling,¹⁷ adaptive biasing force,¹⁸ or crystal-adiabatic free energy dynamics,¹⁹ to name a few. Here, we shall use WTMetaD. In this approach, a bias potential V(s,t) that evolves in time is constructed as a function of a small number of collective variables (CVs) s, which are functions of the atomic coordinates. The objective of V(s,t) is to discourage frequently visited configurations and aid the system to surmount free energy barriers. V(s,t) is defined as a sum of Gaussian functions of height W and width σ. V(s,t) converges asymptotically to²⁰where F(s) is the free energy surface (FES); γ is the so called bias factor, which takes a value greater than 1, and dictates how fast the height of the Gaussian functions decays with time. The time dependent energy offset c(t) is independent on s. The FES can be calculated as a function of a different set of CVs s′ using the reweighting procedure described in ref. 21. The interested reader is referred to ref. 15 and 20 for further details on WTMetaD.

2.3 Surface layer crystallinity

An important ingredient of WTMetaD is the choice of CVs. The CVs are typically non-linear functions of the atomic coordinates and describe the progress (state) of the process under study. In the case of crystallization, a reasonable choice is to use as CV the number of crystalline molecules in the growing layer. This CV would enhance both the growth and the dissolution of a single layer on the surface of the seed crystal. Here, we shall use a simple yet effective modification of the CV described in ref. 8 to localize the CV to the growing surface layer.

The CV should be able to account for the fact that molecules have a different environment in solution and in the crystal.²² Molecules in solution are characterized by a small coordination number and almost random relative orientation of the neighbors. On the other hand, molecules in the crystal have higher coordination numbers and are aligned with respect to their neighbors. With this in mind, we define the surface layer crystallinity of solute molecule i as

s_slc,i = ϕ_iρ_iΘ_i, (5)

where the function ϕ_i localizes the CV to the surface layer, ρ_i describes the local density, and Θ_i the local order.

The CV s_slc,i needs to be continuous and differentiable. Function ϕ_i is defined for each solute molecule i to have a non-zero value if solute molecule i is located within the surface layer. Defining the z-axis as perpendicular to the crystal surface, ϕ_i takes the form

where z_i is the position of molecule i along the z-axis, ζ_c and ζ_l determine the boundaries of the surface layer, while σ_c and σ_l define the steepness of ϕ_i individually on the crystal and on the liquid side of the surface layer. Fig. 2a shows a snapshot of a simulation box of the {001̄} face of naphthalene in which the surface layer is shaded in red. The ϕ_i as a function of z_i for this configuration is shown below.

Fig. 2 a) Snapshot of a simulation box of the naphthalene {001̄} face exposed to the solution (ethanol molecules are shown in faded colors). The red shaded region indicates the region of the growing surface layer. The surface layer crystallinity only takes molecules into account where function ϕ(z_i) is not zero (see graph below). b) Illustration of the vectors that characterize the orientation of naphthalene molecules i and j and the corresponding angle θ_ij. c) Probability density function of cos θ_ij for liquid (orange line) and crystalline (green line) naphthalene. d) Probability density functions of s_slc,i for each naphthalene molecule i in the case of a fully liquid surface layer (orange line) and fully crystalline surface layer (green line).

The local density ρ_i is defined with the following rational switching function

where n_i is the coordination number of molecule i, n_cut a predefined threshold value, and a a parameter that determines the steepness of the switching function. If n_i surpasses n_cut, ρ_i tends to 1, whereas otherwise it tends to 0; n_i is defined aswhere N is the number of solute molecules in the system and f_ij is defined for solute molecule pairs i and j as

Function ϕ_j, which was introduced in eqn (6), provides non-zero values only for neighboring molecules j within the surface layer. The distance between the molecule pair is given by r_ij and r_cut is the cutoff for the calculation of the coordination number.

Now we turn to describe the part of s_slc,i that characterizes the local order. The orientation in space of each molecule is characterized by a vector l. In the case of naphthalene, l is defined along the long axis of the molecule, see Fig. 2b. The relative orientation between solute molecules i and j is characterized through an angle θ_ij defined as cos θ_ij = l_il_j/(|l_i∥l_j|). In Fig. 2c we show the probability density function of cos θ_ij for neighbors ij inside the cutoff r_cut both for the liquid and the crystalline phase of naphthalene. The histogram shows four well defined peaks, each peak corresponding to one of the reference angles present in the crystal. We can now define the local order in eqn (5) as

where σ_κ is the width of the angular distribution associated to . K denotes the total number of reference angles occurring in the crystal. The purpose of f_ij in eqn (10) is to restrict the summation on j only to the neighbors within r_cut. Also, dividing by n_i normalizes Θ_i in the interval [0, 1]. Values of Θ_i close to 1 correspond to an environment where neighbors are aligned according to the reference angles , while values of Θ_i close to 0 characterize environments with random orientations relative to their neighbors.

The combined effect of ϕ_i, ρ_i, and Θ_i leads to a crystallinity measure s_slc,i that takes into account both density and order around the molecule i within the surface layer. Taking the sum over all N solute molecules of the system leads to the surface layer crystallinity

which quantifies the crystallinity of the growing layer and avoids distortions of the bulk liquid as well as the bulk crystal in biased simulations by the use of the function ϕ_i. In Fig. 2d we show probability density functions (PDFs) of s_slc,i for a fully liquefied layer and for a fully grown crystal layer. It is clear from the figure that s_slc,i distinguishes well these two states.

2.4 Biased collective variable

Although s_slc could have been used to construct the WTMetaD bias potential, we found it beneficial to use a slightly different CV. This is due to the fact that the two minima in the FES corresponding to the fully dissolved layer (A) and fully grown layer (B) have rather different widths (see the illustrative scheme in Fig. 3a). In particular basin A is rather narrow and close to the boundary of the CV. The first feature requires the use of very narrow Gaussians. These narrow Gaussians are totally unnecessary in the filling of basin B and thus convergence is slowed down. On the other hand, the closeness of basin A to the boundary creates artificial effects.^23,24

Fig. 3 a) Scheme of a FES, F, in dependence of the surface layer crystallinity s_slc. The local minimum A is the fully dissolved and local minimum B the fully grown surface layer. b) Corresponding scheme of the FES in dependence of the biased CV, s, in which states A and B feature similar widths.

In order to solve these problems, we could have used the variable Gaussian width approach of Branduardi et al.²⁵ Here instead we prefer to use a much simpler fixup. We use as CV not directly s_slc but the power law modification

s = s^χ_slc, (12)

where the exponent χ < 1 is chosen to make the two minima approximately equally broad. A qualitative example of the effect of this change of CV is shown in Fig. 3b. Not only are the two minima similarly broad but the minimum at A has been pushed fundamentally from the s = 0 border. Adapting the CV in this way leads to an improved sampling of the low crystallinity states, while not affecting the sampling of high crystallinity states in an undesired manner.

The parameters used for s_slc, s, and WTMetaD are tabulated in the ESI.†

2.5 Force fields

The general AMBER force field (GAFF)^26,27 was used for naphthalene and ethanol. Both types of molecules have a fully atomistic description. For naphthalene, the relaxed structure and the restrained electrostatic potential (RESP) charges²⁸ were calculated with Gaussian 09 (ref. 29) using the density functional theory structural optimization routine performed at the B3LYP/6-31G(d,p) level. The partial charges of the naphthalene atoms were fitted to the RESP charges. The lengths of the bonds involving hydrogen have been fixed at the equilibrium value. Details on the force field parameters of naphthalene and the force field validation can be found in the ESI.†^30–32 Force field parameters of ethanol were taken from van der Spoel et al.³³

2.6 Simulation runs

All simulations were performed with Gromacs 2016.5 (ref. 35–39) patched with a private version of Plumed 2.5.⁴⁰ The temperature was set to 280 K in line with the experiments reported in ref. 13 and was kept constant with the velocity rescaling thermostat.⁴¹ Periodic boundary conditions were used. The particle mesh Ewald approach⁴² was employed for the long-range electrostatic interactions, and the LINCS algorithm^38,43 was used to constrain the bonds involving hydrogens. The non-bonded cutoff for the van der Waals and electrostatic interactions was set to 1 nm.

For each of the three crystal faces, two concentrations of C = 1.2 nm⁻³ and C = 1.5 nm⁻³ were considered, which correspond to supersaturations of approximately S ≈ 0.5 and S ≈ 0.9 (see ESI† for details). The simulation boxes contained around 1300 molecules. The equilibration procedure⁴⁴ to obtain the initial configuration is reported in the ESI.† To compute the free energy barrier, simulation runs of 1 μs each were performed using WTMetaD with s. Application of the CμMD algorithm ensured the constancy of the chosen concentration. Simulation box visualizations of the three faces, {001̄}, {11̄0}, and {201̄}, are shown in Fig. 4. To keep the enhanced sampling protocol simple, only the right surface was biased at supersaturated conditions. The left control region concentration was kept approximately at solubility to prevent either growth or dissolution of the left crystal surface.

Fig. 4 Simulation box visualization of the biased simulation setups for the three faces, {001̄}, {11̄0}, and {201̄}. Ethanol molecules are shown in faded colors. The visual molecular dynamics (VMD) software was used to visualize the trajectories and prepare the images.³⁴

During the course of the simulation, the crystal slab can experience a small random displacement. This can affect the determination of the position of the growing layer. To avoid this, the position of the center of mass of the two inner crystal layers was constrained using a harmonic potential. Further details can be found in the ESI.†

3 Results

3.1 Steady state crystal shapes

The goal of this work is to predict the steady state crystal shape of naphthalene grown from an ethanol solution. The first step in this direction is to calculate the free energy barrier, , for the growth of each face.¹¹ We shall calculate from the WTMetaD simulations.

We describe a typical simulation taking as example the face {001̄} at concentration C = 1.5 nm⁻³. In Fig. 5a, we show the time evolution of the CV used in WTMetaD, namely s. During the simulation the system explores configurations in which the crystal layer is grown and dissolved. This should be compared to a standard MD simulation in which the system would remain for a much longer time in the dissolved state. From the simulation, the FES can be calculated using eqn (4) and is shown in Fig. 5b (blue curve). As expected, F(s) shows two minima that correspond to the grown and dissolved states. To improve sampling a wall potential was imposed to prevent the growth of the full layer. The wall potential is shown as a red curve in Fig. 5b (see ESI† for details). Representative simulation box snapshots for the two minima are shown in Fig. 5c.

Fig. 5 a) Time evolution of the biased CV, s. The orange shaded region relates to the dissolved crystal layer, while the green shaded region corresponds to the grown layer. b) Corresponding FES, F(s) (blue curve), and wall potential (red curve). c) Representative simulation box snapshots of states at each minimum of the FES.

Although F(s) exhibits a free energy barrier, we decided to calculate the height of the barrier using a CV with a clear physical interpretation. For this reason, we reweighted²¹ our results to calculate the FES as a function of s_slc, which approximates the number of crystalline molecules in the growing layer. In Fig. 6, we show F(s_slc) for all simulations, and the corresponding values of the are summarized in Table 1. For both concentrations, is substantially larger than and . Furthermore, is less sensitive to changes in concentration than and ΔF_{201̄}. For face {001̄}, the free energy barriers differ between the two concentrations by ∼2.4 kJ mol⁻¹, while for faces {110̄} and {201̄} the free energy barriers differ by ∼5.0 kJ mol⁻¹.

Fig. 6 Free energy, F(s_slc), in dependence of crystallinity of the crystal surface, s_slc. Green denotes the concentration at C = 1.2 nm⁻³, and purple denotes the concentration at C = 1.5 nm⁻³.

Table 1 Free energy barriers (in [kJ mol⁻¹])

C = 1.2 nm^−3	24.5	11.2	7.2
C = 1.5 nm^−3	22.1	6.2	2.1

---

It is important to note, that the sampled are strongly affected by the finite size effects.⁸ In a small system, such as the one used in this work, the periodic boundary conditions of the simulation box induce an artificial stabilization of the emerging 2D nucleus on the surface of the crystal. This leads to an underestimation of the critical 2D nucleus size and as a consequence an underestimation of . It is however possible to obtain trends which are consistent with experiments as shown in the following. Also one has to point out that here we are interested in the relative behaviors and we expect that finite size effects cancel out approximately.

Now since we have the free energy barriers, we can compute the relative growth rates for each face and concentration. For the subsequent analysis we shall assume that all faces grow at the chosen supersaturations in the 2D nucleation regime. The perpendicular growth rate of a face {hkl} growing by 2D nucleation can be written as^13,45

where d_{hkl} is the interplanar spacing of the face and τ_{hkl} denotes the time to completely cover the crystal surface with a new layer of naphthalene molecules. Numerical values of d_{hkl} are reported in the ESI.†τ_{hkl} depends on the nucleation rates and edge velocities of the emerging 2D nucleus. If we assume that nucleation is the rate limiting step, then we can approximate τ_{hkl} as the reciprocal of an Arrhenius type function^8,45,46where ξ is a pre-exponential coefficient which is assumed to be the same for all faces.

According to Chernov,^12,47 the crystal morphology in a typical crystallization environment is determined by the crystal growth kinetics. The ratio of the perpendicular distances of faces {hkl} and {h′k′l′} from the center of mass of the crystal, H_{hkl} and H_{{h′k′l′}}, is proportional to the ratio of the growth rates of the two faces, v_{hkl} and v_{{h′k′l′}}

With the ratios H_{hkl}/H_{{h′k′l′}}, the steady state crystal shape can be estimated for each supersaturation by solving a system of plane equations in the Cartesian coordinates.⁴⁸

For both concentrations, the predicted crystal shape is a wafer-thin rhombic prism, as shown in Fig. 7 for C = 1.2 nm⁻³. The blue colored large surface corresponds to the slow growing {001̄} face, while the red colored thin surfaces correspond to the fast growing {11̄0} face. The very fast growing {201̄} face is not present in the resulting steady state crystal shapes.

Fig. 7 Steady state crystal shape obtained from simulation at C = 1.2 nm⁻³. The blue colored large surface corresponds to face {001̄} and the red colored thin surfaces correspond to face {11̄0}.

The activation energy barrier changes more rapidly with supersaturation for face {11̄0} than for face {001̄}. Consequently, the ratio H_{11̄0}/H_{001̄} is about three times larger for C = 1.5 nm⁻³ than for C = 1.2 nm⁻³. Therefore the rhombic prism is about three times thinner for this higher concentration.

The experimental steady state crystal shape at high supersaturation reported by Lovette et al. in ref. 13 shows a wafer-thin rhombic prism. Only faces {001̄} and {11̄0} are present and face {11̄0} grows much faster in comparison to face {001̄}. These trends are captured well by the simulation result at C = 1.2 nm⁻³.

3.2 Growth mechanisms

The biased simulations also allow to gain insight into the details of the growth mechanism. One should keep in mind that the simulations were performed with small surface areas exhibiting finite size effects. However, here again we are interested into trends which are expected to be the same for larger system sizes. Fig. 8 shows the free energy surface, F, as a function of the CV s_slc, that approximates the number of crystalline naphthalene molecules in the surface layer, and CV N_e,sl, i.e. the number of ethanol molecules residing in the surface layer.

Fig. 8 Free energy surfaces, F, in dependence of s_slc, the approximate number of crystalline naphthalene molecules in the surface layer, and N_e,sl, the number of ethanol molecules in the surface layer. The white isocontours are drawn in 2.5 kJ mol⁻¹ energy steps. The red lines correspond to the minimum free energy paths. Fig. 8a)–f) show the free energy surfaces for different faces and supersaturations.

We shall focus on the growth mechanism of face {001̄} first (see for instance Fig. 8a). The dissolved free energy minimum is strongly oriented along the N_e,sl axis. Therefore the system will experience fluctuations in the number of ethanol molecules in the surface before starting to form a crystalline 2D nucleus. Of course, the number of naphthalene molecules in the surface is anticorrelated to the number of naphthalene molecules in the same region. Indeed, we observe in the simulation trajectories that an amorphous agglomerate of naphthalene molecules appears before the onset of crystallinity. Only at a subsequent stage the system surmounts the barrier in the s_slc axis thus forming a crystalline 2D cluster that grows spontaneously.

From these observations we conclude that this face grows through a two step mechanism. The first step is characterized by the formation of an amorphous agglomerate of naphthalene molecules while during the second step the molecules order into the crystalline lattice. This mechanism can be observed for both concentrations, see Fig. 8a and b. The minimum free energy path for the transformation is presented in red in the FES and is compatible with the description given above. The nudged elastic band method⁴⁹ was used for the calculation of the minimum free energy paths.

Now we turn to describe faces {11̄0} and {201̄} that exhibit a growth mechanism different from the one of face {001̄}. We analyze in particular Fig. 8c. In this case the decrease in the number of ethanol molecules in the surface occurs concomitantly with the increase in the number of crystalline naphthalene molecules. This is confirmed by the atomistic trajectories in which we also observe that as soon as naphthalene molecules reach the surface, they simultaneously tend to orient according to the crystal lattice. This description is compatible with a one step mechanism different from the scenario in face {001̄}. From the minimum free energy paths presented in red in Fig. 8 the two different growth mechanisms can be clearly distinguished.

4 Conclusions

In this work, we have used molecular dynamics to predict the steady state crystal shape of naphthalene grown from ethanol solution at two different supersaturations. We calculated the free energy barriers for the growth of faces {001̄}, {11̄0}, and {201̄} using well-tempered metadynamics. In order to mimic the experimental conditions the simulations were performed at constant supersaturation using the CμMD algorithm.¹⁶ From the free energy barriers we could estimate the relative growth rates and with them we predicted the steady state crystal shape. The resulting crystal shape is a wafer-thin rhombic prism, where only the {001̄} and {11̄0} faces are present and the former face grows considerably slower. These predictions are in line with experiments.

We identified two different growth mechanisms. For the {001̄} face, a two step nucleation was observed in which an amorphous agglomerate of naphthalene molecules has to form before the molecules arrange into the crystalline lattice. For the {11̄0} and {201̄} faces the growth occurs with a concerted agglomeration and orientation of naphthalene molecules in the crystalline 2D nucleus.

The framework developed to enhance the growth of a single naphthalene crystal layer under constant supersaturation can be adapted to study the growth of more complicated organic molecules like APIs.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

Z. B. and M. M. are thankful for the financial support provided by Novartis Pharma AG, Lichtstrasse 35, 4056 Basel, Switzerland. P. M. P., T. K. and M. P. acknowledge support from the NCCR MARVEL funded by the Swiss National Science Foundation and from European Union Grant No. ERC-2014-AdG-670227/VARMET. Z. B. would like to thank Matteo Salvalaglio, Dan Mendels, Claudio Perego, and Riccardo Capelli for valuable discussions. Z. B. thanks GiovanniMaria Piccini for supplying the code used to calculate the minimum free energy paths presented in Fig. 8. The computational resources were provided by ETH Zürich and the Swiss Center for Scientific Computing at the Euler Cluster.

References

1. D. Winn and M. F. Doherty, AIChE J., 2000, 46, 1348–1367.
2. J. Li, C. J. Tilbury, S. H. Kim and M. F. Doherty, Prog. Mater. Sci., 2016, 82, 1–38.
3. X. Y. Llu, E. S. Boekt, W. J. Brlels and P. Bennema, Nature, 1995, 374, 342–345.
4. D. Winn and M. F. Doherty, Chem. Eng. Sci., 2002, 57, 1805–1813.
5. S. Piana and J. D. Gale, J. Am. Chem. Soc., 2005, 127, 1975–1982.
6. S. Piana, F. Jones and J. D. Gale, J. Am. Chem. Soc., 2006, 128, 13568–13574.
7. J. Chen and B. L. Trout, Cryst. Growth Des., 2010, 10, 4379–4388.
8. M. Salvalaglio, T. Vetter, F. Giberti, M. Mazzotti and M. Parrinello, J. Am. Chem. Soc., 2012, 134, 17221–17233.
9. T. Karmakar, P. M. Piaggi, C. Perego and M. Parrinello, J. Chem. Theory Comput., 2018, 14, 2678–2683.
10. D. Han, T. Karmakar, Z. Bjelobrk, J. Gong and M. Parrinello, Chem. Eng. Sci., 2019 DOI:10.1016/j.ces.2018.10.022.
11. M. Salvalaglio, T. Vetter, M. Mazzotti and M. Parrinello, Angew. Chem., Int. Ed., 2013, 52, 13369–13372.
12. A. A. Chernov, Sov. Phys. Crystallogr., 1962, 7, 728–730.
13. M. A. Lovette and M. F. Doherty, Cryst. Growth Des., 2012, 12, 656–669.
14. R. F. P. Grimbergen, M. F. Reedijk, H. Meekes and P. Bennema, J. Phys. Chem. B, 1998, 102, 2646–2653.
15. A. Barducci, G. Bussi and M. Parrinello, Phys. Rev. Lett., 2008, 100, 020603.
16. C. Perego, M. Salvalaglio and M. Parrinello, J. Chem. Phys., 2015, 142, 144113.
17. G. M. Torrie and J. P. Valleau, J. Comput. Phys., 1977, 23, 187–199.
18. E. Darve and A. Pohorille, J. Chem. Phys., 2001, 115, 9169.
19. E. Schneider, L. Vogt and M. E. Tuckerman, Acta Crystallogr., Sect. B: Struct. Sci., Cryst. Eng. Mater., 2016, 72, 542–550.
20. O. Valsson, P. Tiwary and M. Parrinello, Annu. Rev. Phys. Chem., 2016, 67, 159–184.
21. P. Tiwary and M. Parrinello, J. Phys. Chem. B, 2015, 119, 736–742.
22. E. E. Santiso and B. L. Trout, J. Chem. Phys., 2011, 134, 064109.
23. C. Micheletti, A. Laio and M. Parrinello, Phys. Rev. Lett., 2004, 92, 170601.
24. Y. Crespo, F. Marinelli, F. Pietrucci and A. Laio, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2010, 81, 055701.
25. D. Branduardi, G. Bussi and M. Parrinello, J. Chem. Theory Comput., 2012, 8, 2247–2254.
26. J. Wang, R. M. Wolf, J. W. Caldwell, P. A. Kollman and D. A. Case, J. Comput. Chem., 2004, 25, 1157–1174.
27. J. Wang, W. Wang, P. A. Kollman and D. A. Case, J. Mol. Graphics Modell., 2006, 25, 247–260.
28. C. I. Bayly, P. Cieplak, W. Cornell and P. A. Kollman, J. Phys. Chem., 1993, 97, 10269–10280.
29. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G. A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino, G. Zheng, J. L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J. A. Montgomery, J. E. Peralta, F. Ogliaro, M. Bearpark, J. J. Heyd, E. Brothers, K. N. Kudin, V. N. Staroverov, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, N. Rega, J. M. Millam, M. Klene, J. E. Knox, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, R. L. Martin, K. Morokuma, V. G. Zakrzewski, G. A. Voth, P. Salvador, J. J. Dannenberg, S. Dapprich, A. D. Daniels, O. Farkas, J. B. Foresman, J. V. Ortiz, J. Cioslowski and D. J. Fox, Gaussian 09, Revision B.01, 2009.
30. A. W. Sousa da Silva and W. F. Vranken, BMC Res. Notes, 2012, 5, 367.
31. D. W. J. Cruickshank, Acta Crystallogr., 1957, 10, 504–508.
32. D. H. Andrews, G. Lynn and J. Johnston, J. Am. Chem. Soc., 1926, 48, 1274–1287.
33. D. van der Spoel, P. J. van Maaren and C. Caleman, Bioinformatics, 2012, 28, 752–753.
34. W. Humphrey, A. Dalke and K. Schulten, J. Mol. Graphics, 1996, 14, 33–38.
35. H. J. C. Berendsen, D. van der Spoel and R. van Drunen, Comput. Phys. Commun., 1995, 91, 43–56.
36. E. Lindahl, B. Hess and D. van der Spoel, J. Mol. Model., 2001, 7, 306–317.
37. D. van der Spoel, E. Lindahl, B. Hess, G. Groenhof, A. E. Mark and H. J. C. Berendsen, J. Comput. Chem., 2005, 26, 1701–1718.
38. B. Hess, C. Kutzner, D. van der Spoel and E. Lindahl, J. Chem. Theory Comput., 2008, 4, 435–447.
39. M. J. Abraham, T. Murtola, R. Schulz, S. Pall, J. C. Smith, B. Hess and E. Lindahl, SoftwareX, 2015, 1–2, 19–25.
40. G. A. Tribello, M. Bonomi, D. Branduardi, C. Camilloni and G. Bussi, Comput. Phys. Commun., 2014, 182, 604–613.
41. G. Bussi, T. Zykova-Timan and M. Parrinello, J. Chem. Phys., 2009, 130, 074101.
42. T. Darden, D. York and L. Pedersen, J. Chem. Phys., 1993, 98, 10089–10092.
43. B. Hess, J. Chem. Theory Comput., 2008, 4, 116–122.
44. M. Parrinello and A. Rahman, J. Appl. Phys., 1981, 52, 7182–7190.
45. X.-Y. Liu and P. Bennema, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 53, 2314–2325.
46. S. Arrhenius, Z. Phys. Chem., 1889, 4, 226–248.
47. 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.
48. Y. Zhang, J. P. Sizemore and M. F. Doherty, AIChE J., 2006, 52, 1906–1915.
49. G. Henkelman and H. Jonsson, J. Chem. Phys., 2000, 113, 9978–9985.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c9ce00380k

---