Predictive Calculation of the Crystallization Tendency of Model

Article pubs.acs.org/JPCB

Predictive Calculation of the Crystallization Tendency of Model Pharmaceuticals in the Supercooled State from Molecular Dynamics Simulations J. Gerges and F. Affouard* Unité Matériaux et Transformations (UMET), UMR CNRS 8207, UFR de Physique, BAT P5, Université de Lille 1, 59655 Villeneuve d’ascq, France ABSTRACT: Molecular dynamics (MD) simulations were used to perform a comparative study of the crystallization tendency from the melt of two model pharmaceutical compounds: felodipine and nifedipine. Two crystalline polymorphs of nifedipine (Nα, Nβ) and felodipine (FI, FII) have been studied. Calculations were performed on liquid and crystal systems separately in order to determine their main physical properties (diffusivity, density, and enthalpy). A fair agreement was found between the simulation and the known experimental data confirming the ability of the force field GAFF to reproduce accurately the experimental data for both compounds. Simulations of the crystal−liquid interface enabled the determination of the melting temperature and the interfacial free energy of the different polymorphs. Guided by the classical nucleation theory (CNT) predictions and different growth mechanism models (normal, twodimensional, and screw dislocation), the nucleation and growth rates have been determined. The present investigation particularly raises the very important role of the solid−liquid interfacial free energy and its interplay with the driving force during the crystallization. The origin of the higher crystallization tendency of nifedipine with respect to felodipine is discussed from the present computed kinetic and thermodynamical factors.

I. INTRODUCTION Upon cooling, there are two basic ways for a liquid to solidify: it may form a crystalline solid in which atoms are arranged in an orderly manner on a lattice; or it may form an amorphous solid, i.e., a glass, characterized by an atomic disorder and the absence of a long-range order. The latter can only be formed if the liquid is cooled rapidly enough with respect to the kinetics associated with crystallization. Although thermodynamically unstable, a solid can remain amorphous for a long time until its transformation into a stable crystalline state having a lower Gibbs free energy because the atoms do not have enough mobility to rearrange.1,2 Despite considerable attention and great importance to both basic science and applications for processing materials, the tendency for a given material to crystallize, vitrify, or revert to the more stable crystalline state from the amorphous state remains poorly understood.3−5 This fundamental issue has initiated intense activity. Recent developments and approaches in many different fields particularly focus on predicting and possibly controlling the conditions that will favor/unfavor the crystallization/vitrification depending on if it is desirable or not. For example, these transformations have been extensively studied and reviewed throughout the metallurgic literature.6,7 Indeed, understanding the glass-forming ability of metallic glasses is a long-standing fundamental problem.6−9 The possibility to obtain excellent glass-forming ability at a very low cooling rate for complex multicomponent alloys has clearly aroused much interest. In © XXXX American Chemical Society

pharmaceutical industries, many substances are preferably developed in the crystalline state for evident reasons of stability. However, either by accident or design, they may also exist in a total or partially amorphous state. This situation is encountered more and more frequently due to the increasing complexity of synthesized molecules. The amorphous state has thus become increasingly important in pharmaceutical applications owing to its higher apparent solubility and faster dissolution rate. There has been a considerable amount of work on the so-called amorphous solid dispersions aiming to predict the persistence of the active pharmaceutical ingredients (API) in the amorphous state as a function of temperature.10−12 The main motivation is to prevent recrystallization upon long-term storage since this transformation would obviously negate advantages of using the amorphous state. 13−18 Most pharmaceuticals are composed of small organic molecules showing a high contrast between intra- and intermolecular forces with respect to other classes of materials. Motivated by the need to predict if crystallization/vitrification will occur or not for a large number of materials including pharmaceuticals, many attempts have been made to derive general rules from easily accessible structural, dynamical, or thermodynamical experimental properties and thus define some quantitative Received: June 10, 2015 Revised: July 20, 2015

A

DOI: 10.1021/acs.jpcb.5b05557 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

the energy barrier associated with atomic motion. Both probabilities behave in opposite ways as a function of temperature and thus give rise to the bell-shaped dependence of the nucleation rate. Diffusion thus allows the nucleation phenomenon to occur but also modulates its kinetics.40,41 Growth is the second stage in which the supercritical nuclei formed will grow within the liquid phase.1,2,36−38,42 It occurs by aggregation of molecules of the phase undergoing a transformation. Crystal growth velocity, G (unit m s−1), is a measure of the growing of the crystal surface of the supercritical nuclei. Similarly to the nucleation rate, G depends on the product of a thermodynamic term that describes the probabilities of attachment/detachment of the molecules to/from the nuclei [1 − exp(−ΔG/kBT)] and the probability of the diffusion of the molecules across the interface which is also roughly proportional to the diffusivity, D. The net effect of combining these two terms is that the growth rate must go through a maximum and thus also give rise to a bell-shaped dependence as a function of temperature. The maximum of the growth rate can be different from the maximum associated with the nucleation rate. In other words, the curves representing their temperature dependencies can either overlap or not. When the curves representing the temperature dependencies of nucleation and growth rate do not significantly overlap, good glassforming ability is to be expected. Poor glass-formers would have a large temperature overlap in which both nucleation and growth rates are high. The knowledge of the two nucleation and growth rates as a function of temperature thus provides a guide to predict the crystallization tendency. It could be noted that from the combination of both nucleation and growth rates (assuming a constant variation as a function of time), it would be possible to determine the so-called TTT curve for a given crystallized volume fraction as it is frequently done in metallurgy.29,43 This relation links the degree of supercooling and the crystallization time. The glass-forming ability is directly determined by the position of the “nose” of the TTT curve. However, it is a difficult task experimentally. Indeed the direct study of the early stages of nucleation presents significant experimental difficulties since the embryos are small and have a transient nature. It is thus challenging to observe this phenomenon with enough precision by visual methods, and only a few studies have been performed especially for molecular compounds.16,37 Overall the phase transformation induced by cooling is the interplay of N and G which are basically determined to some extent by three main physical ingredients: (i) the driving force, ΔG; (ii) the diffusivity, D; and (iii) the interfacial free energy, γ.34,39 The first two parameters can be estimated without major difficulty from experimental data. ΔG can be usually obtained either from the temperature dependencies of the heat capacities of the liquid and crystalline phases if experimental data are available with enough accuracy or from analytical expressions derived from the melting enthalpy and temperature such as Turnbull or Hoffman equations.44,45 Experimental diffusivity data are often scarce, but assuming the validity of the Stokes− Einstein relation, D can be replaced by the much more experimentally accessible shear viscosity, although this replacement is questionable due to possible decoupling of diffusion and relaxation.2,46 The determination of the last parameter, i.e., the crystal−liquid interfacial free energy, which is a fundamental parameter for nucleation and is involved in several mechanisms of growth,6,7,39 is problematic. Despite its importance to theoretical models and to practical applications,

predictors. For example, it was suggested that good glassforming ability is linked to the complex composition of molecules (such as metallic alloys19), low molecular symmetry, and high flexibility due to a high number of internal degrees of freedom.20 The role of the molecular mobility and its temperature dependence has been also pointed out. The mobility versus the steepness of the temperature curve is usually measured by the so-called fragility parameter defined as the slope at the glass transition temperature Tg of the viscosity or the primary relaxation time obtained from dielectric relaxation experiments.21 Recent results suggest a weak correlation between the glass-forming ability and the fragility which could be also highly dependent on the considered class of materials. While this correlation seems accurate in metallic materials, it seems more questionable in other molecular glassformers.18,22−26 It can also be noted that some pharmaceuticals may show same molecular mobility but completely different crystallization tendencies.24,27 Thermodynamically, a low melting entropy, ΔSm, and a limited range between the melting temperature, Tm, and the glass transition temperature, Tg, appear to be an advantage for glass formation.20,22,28,29 For a long time, it has been wellknown that the ratio (Tg/Tm) > 2/3 may provide a fair prediction of the glass-forming ability,30,31 but there are some counterexamples.32,33 No predictors can be, unfortunately, considered as fully reliable14,29,34,35 since the structure, the dynamics, and the thermodynamics seem interrelated. Moreover, in order to really investigate the influence of one specific parameter, the others should be the same which is not the case in general. From a fundamental point of view, crystallization in the undercooled melt is a phenomenon involving complex nucleation and growth stages, both having their own specific properties.34 Nucleation is the first stage. It can be described in the framework of the classical nucleation theory (CNT) which provides a reliable basic description despite significant shortcomings amply reviewed in the literature.1,2,36−38 Nucleation is related to the stochastic spontaneous formation of small crystalline nuclei made of few molecules in the liquid. The appearance of a nucleus of the new phase implies the creation of an interface between liquid and solid phases. This transformation is thermodynamically controlled by a subtle balance between a favorable term, the so-called the “driving force” for the phase transformation ΔG which is the difference between the Gibbs free energies of the liquid and the solid and an unfavorable term, the crystal−liquid interfacial free energy, γ, that is the resistance to the creation of the solid−liquid interface. This latter term mainly originates from an entropy loss due to the increased ordering of the liquid as the crystal nucleus surface is approached.39 This balance gives rise to an activation energy barrier (ΔG*) that needs to be overcome to produce a nucleus having a size larger than a minimum critical size needed for the nucleus to expand; otherwise it will shrink. In addition, diffusion is also required to make this process possible since molecules have to aggregate to form the nucleus. Nucleation is often characterized by the nucleation rate, N (unit m−3 s−1), which is the number of formed critical nuclei per unit of volume and time. It depends on the product of two activated processes: (i) the number of critical size nuclei ∼ exp[−ΔG*/ kBT] and (ii) the probability of a molecule to migrate across the interface separating the critical size nucleus and the liquid (thus forming a supercritical size nucleus) which is roughly proportional to the diffusivity D ∼ exp[−A/kBT] , where A is B

DOI: 10.1021/acs.jpcb.5b05557 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B accurate data for γ are in general not known, even for simple cases.6,7 Only a few interfacial free energy values are available in the literature for pharmaceuticals and more generally for molecular compounds.37,47−49 A further difficulty is that, unlike a liquid−gas interface, the crystal−liquid interfacial free energy may show anisotropy and dependence on the orientation of the crystalline face in contact with the liquid. There are few experimental techniques capable of measuring this fundamental quantity.2,44,50−56 It is therefore not surprising that the values determined experimentally differ considerably depending on the specific technique used.54,57 The determination of the parameter γ from simulation thus presents a clear interest but also remains challenging.6,7,58,59 In recent years, several approaches based on molecular dynamics (MD) simulations have been proposed and validated mostly on simple systems composed of hard spheres, Lennard-Jones particles, metal atoms, or small molecules such as water. Broughton and Gilmer60 proposed a method called “adiabatic cleaving” based on the fact that the free energy of the interface is a thermodynamic state function. Using this approach, they were able to calculate the reversible work required to create a solid−liquid interface, i.e., the interfacial free energy for a set of particles interacting via a Lennard-Jones potential. This approach has been extended by Davidchack and Laird61,62 to take into account the anisotropy of the crystal interface of the compounds mainly made either of hard spheres61 or LennardJones particles.62 Handel et al.63 extended the cleaving method to molecular systems where they performed direct calculations of the ice−water interfacial free energy for the TIP4P model. Moreover, the critical nucleus method (CNM)64 is based on the classical nucleation theory. γ is estimated from the Gibbs− Thomson effect and from the determination of the critical undercooling corresponding to a given critical radius. The thermodynamic integration scheme is a more recent method used to determine γ for Lennard-Jones systems65 and for hard spheres.66 Six consecutive steps were used in this method to confine the interface by short-ranged Gaussian walls enabling them to overcome the hysteresis problems that generally occur in the reported previous integration schemes.62,67 An alternative approach proposed by Hoyt et al.,68,69 called “the capillary fluctuation method”, has been successfully validated in recent years for monatomic70 and binary atomic simple systems71 and more realistic models such as metallic compounds,58,72 alloys,73 and molecular materials such as succinonitrile.48 In this approach the interface is intrinsically treated as “rough”. The amplitude of the fluctuations of the position of the crystal−liquid interface is used as a measure of its roughness in order to determine its stiffnesssoft interfaces have more fluctuationswhich can be directly correlated to the interfacial free energy. For a few years, due the importance of the determination of γ, lots of new theoretical approaches have been developed to tackle the issue of the interfacial free energy via computer simulations.74 The most recent approaches are based on metadynamics,75 superheating−supercooling method,76 Gibbs−Cahn integration for single component systems,77 extraction of the interfacial free energy from entropy/energy changes across the interface,74 the Gibbs−Thomson approach,59,78,79 and the mold integration method.80 All of these methods seem promising, but they have only been applied on model systems so far. It should be also mentioned that a major issue in most methods is the definition of a local parameter at the molecular scale that can distinguish

very clearly between the liquid and the solid phases in order to obtain the position of the interface.59,69,72,81 By combining guidance from the available theories of nucleation and growth with calculations of fundamental parameters (ΔG, D, and γ) obtained from MD simulations there is thus an interesting opportunity to investigate the conditions that will favor crystallization/vitrification of low molecular weight molecules. Indeed, in addition to γ above all, MD simulations provide at the same time a relatively small computational effort to determine the diffusivity and a slightly larger effort to obtain ΔG.82−85 This approach avoids the difficulty of simulating the crystallization process itself. Crystallization is an extremely improbable event in the current accessible MD time scales, i.e., nanoseconds to microseconds, and in the MD typical system size containing a few thousand atoms, i.e., system volume of about 10−18 cm3. Direct approaches to crystallization from the melt have thus been restricted to model systems such as monatomic LennardJones.35,86−88 However, it should be mentioned that several computer simulation studies were also performed on nucleation89,90 and growth91,92 for molecular systems but from solution. In this work, by means of MD simulations, we performed a comparative study of two well-known pharmaceutical systems (see Figure 1): nifedipine (C17H18N2O6) and felodipine

Figure 1. Chemical structure of felodipine (a) and nifedipine (b).

(C18H19Cl2NO4), known to have different remarkable crystallization tendencies.15,16,93,94 Both drugs usually serve as model systems for studying the crystallization of organic glasses. It has been reported that nifedipine crystallizes more readily than felodipine either from the melt alone or in the presence of poly(vinylpyrrolidone) (PVP).15,16,93,94 The significant easier crystallization capability of nifedipine is difficult to clearly understand owing to the observed similarities of both materials in molecular structure and molecular mobility (fragility and glass transition temperature (Tg ≈ 43−49 °C15,16,94). Both drugs also exhibit a rich crystalline polymorphism, and polymorphic selectivity of crystallization is observed depending on temperature.95 In nifedipine, three crystalline polymorphs can be obtained from the melt: a metastable form, noted β that is preferentially obtained from the glassy state, and another metastable polymorph that crystallizes concomitantly with the β form from the supercooled melt and transforms to the β form at room temperature. The most stable form noted as α only grows preferentially above 120 °C.5,16,93,94 In felodipine, two polymorphs are currently well established and have been reported to be produced from the melt: the metastable form II that grows preferentially at all temperatures and the most stable form I that may also grow close to heterogeneities.96 The most stable form I grows faster than the metastable form II at all temperatures.13,15,97 To the authors’ best knowledge, it can be noted that additional forms may be also obtained from coC

DOI: 10.1021/acs.jpcb.5b05557 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Table 1. Comparison between Experimental Data and Simulation Data Obtained Using the GAFF100 Force-Field for Nifedipine α and β and Felodipine I and II Crystal Polymorphsa polymorph nifedipine nifedipine felodipine felodipine

α (Nα) β (Nβ) I (FI) II (FII)

ρcx(exp)4,5/ρcx(sim) (at T (K)) (g/cm3)

Tm(exp)4,5/Tm(sim) (K)

ΔHm(exp)4,5,16/ΔHm(sim) (kJ/mol)

1.38/1.37 (283−303) (1.39 to 1.36)/1.37 (296−298) 1.45/1.44 (123) 1.42/1.40 (150)

444/448 −/445 416/426 416/410

39.9/43.9 −/43.5 31.5/31.8 29.1/29.3

γm(sim) (mJ/m2) 21.5 14.4 28.7 15.5

± ± ± ±

1.6 1.0 2.5 1.2

The exp and sim between parentheses indicate the experimental and the simulated values, respectively. ρcx, Tm, ΔHm, and γm represent the density, the melting temperature, the enthalpy, and the interfacial free energy at the melting temperature, respectively. Uncertainties of the interfacial free energy, γm, have been estimated from taking into consideration the uncertainties on the melting temperature (10 K) and the fitting.

a

crystallization4 in solution but not from the melt. A wellcharacterized form III, almost isoenergetic to form I, has been also thus obtained. Recent studies pointed out the possible important role of the often overlooked interfacial free energy in addition to molecular mobility and thermodynamic driving force.6,29,34,35 In this work, the capillary fluctuation method is thus used to determine the crystal−liquid interfacial free energies of nifedipine (forms α and β) and felodipine (forms I and II) polymorphs. In addition, D and the Gibbs free energy difference, ΔG, between the crystalline and the liquid phase have been also determined. It allows us to discuss the relative importance of these different parameters during nucleation and growth and thus gain a better understanding of the factors that govern the crystallization of these pharmaceuticals as a function of temperature This work is organized in the following manner. Section II presents the simulation details that enabled us to determine the main parameters of the nucleation and growth processes (ΔG, D, and γ). Section III is dedicated to the results and the comparison of the estimated values to the experimental data with discussions related to previous work. Conclusions and perspectives are given in section IV based upon our results.

free energy (ΔG) between the liquid and the crystalline state), the self-diffusion coefficient (D) and the interfacial free energy (γm) between the crystal and its melt. A. Calculation of Some Structural, Dynamical, and Thermodynamical Properties of the Liquid and Crystalline States via MD Simulations. The crystallographic data for nifedipine α (Nα), nifedipine β (Nβ), felodipine I (FI), and felodipine II (FII) polymorphs were extracted from refs 4 and5 and from information taken from the Cambridge Structural Database.101,102 The crystalline structures of FI, FII, and Nα are monoclinic with respective space groups P21/c, C2/c, and P21/c while the structure of Nβ is triclinic with a space group P1̅. Small simulation boxes of felodipine and nifedipine crystalline polymorphs were generated with about 100 molecules. For each system, a series of short MD simulations of 200 ps in the NPT statistical ensemble were performed from 100 to 700 K to obtain the evolution of the density and the enthalpy of the crystal as a function of temperature. A good agreement is found between the numerical and the simulated data as reported in Table 1. At high temperature the melting of the crystalline structure is observed. The resulting liquid system is then cooled to a low temperature. At each temperature upon cooling, the system was first equilibrated in the NPT ensemble for 1 ns and then a production run was performed in the NVT ensemble for another 3 ns at the equilibrated volume. These longer runs allowed us to obtain the self-diffusion coefficient D from the calculation of the mean square displacement ⟨r2(t)⟩. For each temperature the long-time behavior of ⟨r2(t)⟩ was fitted using the linear relation ⟨r2(t)⟩ = 6Dt in order to obtain D. Since no experimental values were available for diffusivity, we have also calculated the shear viscosity from the stress−stress autocorrelation functions (Green−Kubo relations)103 even if this parameter will not be used in the following to discuss crystallization. This calculation enabled us to check the fair ability of the present model to reproduce transport properties of both felodipine and nifedipine33 (see Figure 2). These simulations also allowed us to determine the density and the enthalpy of the liquid as a function of temperature upon cooling. B. Calculation of the Melting Temperature via MD Simulations. The numerical melting temperature was estimated using solid−liquid coexistence simulations104−108 requiring MD simulations of a biphasic crystal−liquid structure. For each polymorph, a large crystalline simulation box has been generated from the experimental crystallographic data4,5 but with a larger number of molecules than the previous section ranging from about 1200 to 1700. Each system was first equilibrated at the known experimental melting temperature for 100 ps in the statistical ensemble NPT and then for another 100 ps in the NVT ensemble using the equilibrated volume obtained from the NPT simulation. This crystalline box was

II. SIMULATION DETAILS Molecular dynamics simulations have been performed using the DLPOLY package98 and the force field GAFF99 (general amber force field). This force field was chosen due to its capability of reproducing successfully a large number of experimental data for molecules with low molecular weight.99,100 Simulations have been conducted either in the NPT or the NVT statistical ensemble where N is the number of molecules, P the pressure, T the temperature, and V the volume. The number of molecules, N, was fixed in the simulations, and pressure and temperature were controlled with a Berendsen barostat and thermostat, respectively. All NPT simulations were realized at atmospheric pressure. The equilibrated volume at a given P and T was used for subsequent NVT simulations. The time step to integrate Newton’s equation of motion was chosen as 0.001 ps. A cutoff radius of 10 Å has been used to calculate short-range van der Waals interactions. Ewald summation was employed in order to calculate the long-range electrostatic interactions with the same cutoff radius. Periodic boundary conditions were applied in all directions. In order to investigate nucleation and growth tendencies of nifedipine and felodipine, different systems composed of either molecules in the crystalline or the amorphous state have been simulated. Then, pertinent involved parameters were then extracted from simulation data such as the crystal and the liquid densities as a function of temperature (ρcx(T) and ρlq(T)), the enthalpy of melting (ΔHm that further allows the determination of the driving force of crystallization, i.e., the difference in Gibbs D

DOI: 10.1021/acs.jpcb.5b05557 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

shows the evolution of the total density of a biphasic system of the crystalline Nα polymorph in contact with its melt at different temperatures. Tm was estimated to be the temperature at which the density remains roughly stable as a function of the simulation time, in other words at the temperature where the liquid and crystal coexist. In the case of Figure 3 the melting of Nα was estimated to be 450 K. Table 1 shows the comparison of the experimental and the numerical melting temperatures for the different polymorphs. A very fair agreement is obtained. The method in use predicted a slight overestimation of the melting temperatures, but the correct tendencies of the melting temperature for the different polymorphs are maintained. The determination of the melting temperature allowed us to estimate the enthalpy of melting from the difference of the enthalpy of the crystal and the liquid according to the MD simulations data. A fair agreement is also seen with respect to the experimental data confirming the reliability of the GAFF force field (see Table 1). C. Calculation of the Interfacial Free Energy, γm, via MD Simulations. In spite of its great importance as a key parameter in crystallization processes,39,44 the experimental and numerical determination of the solid−liquid γ remains very challenging.60,69,70,109 In the following, we used the capillary fluctuation method. This method is precisely described in refs 48, 68−70, 110, and 111. So we give here only essential details. The typical geometry of the required system needed to perform such calculations is shown in Figure 4a. The longest direction is perpendicular to the average interfacial plane. The other two directions parallel to the plane are the thickness b and the width w, where b ≪ w in order to obtain a quasi-one-dimensional height function h(x) measuring the position of the interface (see Figure 4a). The interfacial free energy is determined from the power spectrum of the fluctuations of the interfacial position. The amplitudes of the capillary fluctuation modes in a quasi-one-dimensional interface are related to the interfacial stiffness γ̃ by the following relation:

Figure 2. Evolution of the shear viscosity and diffusivity (in inset) as a function of the inverse of temperature for nifedipine and felodipine in the liquid state. Experimental data have been extracted from ref 33. The diffusion coefficients are determined from the long-term evolution of the mean square displacement (⟨r2(t)⟩ ∼ 6Dt)). Shear viscosity, η, has been computed from the stress−stress autocorrelation functions. In the inset, solid lines indicate fits using eq 4. A1 = 18.17, A2 = 491.37 K, and A3 = 874.81 K are obtained for nifedipine while A1 = 18.03, A2 = 445.84 K, and A3 = 874.81 K are obtained for felodipine.

then melted at high temperature (∼800 K), subsequently cooled at the experimental melting temperature, and then equilibrated at this temperature using first NPT then NVT simulation for 100 ps each. Both crystal and liquid boxes are combined together in order to form only one biphasic system. For all investigated systems, the orientation of the crystal in contact with the liquid is [100] (010) [along the interface] (normal to the interface). A small gap is initially set between the crystalline and liquid parts in order to avoid overlapping between the atoms. The whole system is then equilibrated for 10 ps in order to remove this gap. Due to the periodic boundary conditions, a solid−liquid− solid structure was finally obtained. Finally, this biphasic system was simulated at different temperatures above and below the experimental melting temperature. As an example, Figure 3

⟨|h(q)|2 ⟩ =

kBTm bwγ q̃ 2

(1)

where Tm is the melting temperature, h(q) is the onedimensional Fourier transform of h(x) with q as the wavenumber, the symbol represents the time average, and kB is the Boltzmann constant. The value of the function h(x) is measured at discrete set of points xn = nΔ, where n = 1, ..., M and Δ = w/M. In the present study we found that M = 32 is well-adapted to the simulation box dimensions. Hence the Fourier modes hq in eq 1 are defined as hq = (1/M)∑M n=1 h(xn)eiqxn, where q is the wavenumber equal to 2πk/w with k = 1,..., M. The value of the interfacial stiffness can be determined by fitting ⟨|h(q)|2⟩ to a function with a quadratic form in q. In order to reduce the computational effort, in the following we will identify the interfacial stiffness γ̃ determined using one single orientation to an isotropic γ. It is usually reported in many model systems69,71,72,109 and realistic molecular systems48,63,112,113 that the anisotropic effect is usually relatively weak of the order of a few percent and smaller than the uncertainties of the employed method. For example, for the cubic molecular solid (succinonitrile),48 the interfacial free energy has been estimated to be 7.06 ± 0.4, 7.0 ± 0.4, and 7.01 ± 0.4 mJ/m2 for the (100), (110), and (111) orientations, respectively. For the hexagonal ice Ih−water interface system in

Figure 3. Evolution of the density of the biphasic crystal/liquid Nα system as a function of time at three different temperatures. The increase (at T = 350 K) or decrease (at T = 500 K) of density indicates crystallization and melting, respectively. The temperature (T ∼ 450 K) at which the density remains roughly stable as a function of simulation time, in other words at the temperature where the liquid and crystals coexist, was estimated to be the numerical melting temperature, Tm. E

DOI: 10.1021/acs.jpcb.5b05557 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 4. (a) Snapshot of a simulation box used to determine the interfacial stiffness of felodipine FII. The interface is parallel to the (xz) plane. The system size is as follows: w = Lx = 224.31 Å, b = Lz = 26.05 Å, and Ly = 283.20 Å (see text). Lx, Ly, and Lz are the dimensions of the whole simulation box along the x, y, an z directions. (b) Evolution of the order parameter ⟨q6⟩ along the direction orthogonal to the interface (y). The line indicates the fitting using eq 2.

two crystal−melt interfaces due to boundary periodic conditions for which the profile is described by a function h(x). In order to determine the interface position, i.e., interface profile, a local order parameter must be used to distinguish the liquid-like molecules from the solid-like molecules. Due to the complexity of the present molecular systems and their crystalline symmetries, a simple rotational-invariant order parameter ⟨q6⟩58,81,114−117 is used. It is based on a spherical harmonics development of vectors joining the center of mass of neighboring molecules. Generally this order parameter is used for systems having a cubic symmetry and allows distinguishing between different crystal cubic systems, but it can be used as well for noncubic symmetries.116,117 In this work, this parameter allows us to distinguish the liquid and the solid phases in order to obtain the position of the interface. Figure 4b shows an example of the evolution of ⟨q6⟩ along the direction perpendicular to the interface (y direction in Figure.4b) for one instantaneous configuration of the FII system for a given layer [x, x + Δ]. This figure shows clearly the existence of two domains allowing the identification of a crystal part (high values roughly above ⟨q6⟩ ≈ 0.15) and a liquid part (low values below ⟨q6⟩ ≈ 0.06) separated by an interface. In order to precisely determine the position of the interface h(x), the evolution of ⟨q6⟩ has been fitted to an hyperbolic tangent function:58

which molecules may form a high number of hydrogen bonds, the anisotropy was also found to be relatively weak: 23.3 ± 0.8, 23.6 ± 1, and 24.7 ± 0.8 mJ/m2 for the basal, prism, and (112̅0) orientations, respectively.63 Moreover, it can be mentioned that in the present work the chosen crystal orientation in contact with the liquid, [100](010), has approximately the same planar density for all polymorphs (1.2 × 10−2, 1.4 × 10−2, 1.4 × 10−2, and 1.3 × 10−2 molecules/ Å2 for Nα, Nβ, FI, and FII, respectively) which made possible the comparison of the different interfacial free energies. It should be also noted that this property is not completely verified by some of the other faces for which PD may reach (1.2−1.8) × 10−2 molecules/Å2 depending on the considered face and polymorph. These other faces may possibly possess a different interfacial free energy. It is supposed by the authors that the difference remains within the uncertainties of the method as shown for succinonitrile and water. Note that the systems’ anisotropy is a very important issue and it would clearly deserve further studies. In the present study, the crystal-melt biphasic system was constructed similarly to the simulation box previously used to estimate the melting temperature. First the crystal and liquid phases were equilibrated separately at Tm with similar crosssection b × w (see Figure 4a). Then the two systems were joined together along the y direction to form an interface where a small gap region is created to avoid overlapping. The total system is then equilibrated in a NPT ensemble for 50 ps while atomic positions of solid molecules were fixed leaving the melt molecules to fill in the gap. Runs of 100 ps each were then conducted to equilibrate the whole system freeing the solid molecules. These simulations were followed by an NVT simulation of 100 ps for additional equilibration and a production run of 500 ps in which 5000 configurations were stored for analysis. This procedure leads to the formation of

⟨q6⟩(y) =

qs + ql 2

+

⎛ y − h1(x) ⎞ ⎛ y − h2(x) ⎞⎞ qs − ql ⎛ ⎜⎜tanh⎜ ⎟ + tanh⎜ ⎟⎟⎟ 2 ⎝ δ1 δ2 ⎝ ⎠ ⎝ ⎠⎠

(2)

This function describes the shape of a flat interface with effective width δ; qs and ql are the average values of the ⟨q6⟩ parameter in the solid and the liquid domains, respectively; and h(x) is the mean position of the interface. Note that the subscripts 1 and 2 represent the two interfaces created due to the boundary conditions. Figure 5 shows an example of the F

DOI: 10.1021/acs.jpcb.5b05557 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

h(q)|2⟩ versus q. The interfacial free energies obtained by the fitting curves for different polymorphs are summarized in Table 1. The error bars are calculated taking into consideration the uncertainties on the estimated melting temperature (∼10 K) and the fitting procedure. As compared to previous studies on molecular compounds using the capillary fluctuation method,48,112 the linear domain is more limited and the data are much noisier. It might indicate that some interfaces are not completely rough or a possible size effect. This latter was investigated by simulating a Nα system that contains twice more atoms than the system reported in Table 1 corresponding to a 20% larger interface width w. The fit was not really improved in this case. The interfacial free energy was slightly higher than the value calculated for the smaller box (23.7 ± 3.8 mJ/m2), but the difference remains within the error bars.

evolution of the height function h(x) represented for two interfaces of the FII system.

III. RESULTS AND DISCUSSION In the following, we will first present comparative results of the basic factors, i.e., thermodynamic driving force, mobility, and crystal−liquid interfacial free energy, of the investigated compounds based on the present MD simulations. Then, we will employ the general mathematical forms developed to describe homogeneous nucleation and crystal growth rates in order to discuss crystallization tendencies. A. Driving Force, Diffusivity, and Crystal−Liquid Interfacial Free Energy. Here, we used the Hoffman equation45 to estimate ΔG between the crystal and the liquid states as a function of temperature for the different investigated systems. It is given by

Figure 5. Evolution of the interface position h(x) corresponding to the same interface as described in the caption of Figure 4.

Parts a and b of Figure 6 show the behavior of ln(⟨|h(q)|2⟩), the power spectrum of h(x), as a function of ln(q) for nifedipine and felodipine polymorphs, respectively. A small linear range (with a slope of −2) is seen corresponding to the domain where eq 1 is valid. The value of the interfacial free energy can now be obtained by the intercepts of plots of ⟨|

ΔG = ΔHm

(Tm − T )T Tm 2

(3)

where ΔHm and Tm are the simulated melting enthalpy and melting temperature obtained as described previously. Note that in the following ΔGv is the Gibbs free energy difference per unit volume between crystal and liquid states. This equation was already shown to reproduce with good agreement the experimental data determined by specific heat thermodynamic integration for nifedipine and felodipine, respectively.16 Figure 7 shows the evolution of the driving force ΔG as a function of temperature. The expected tendency is obtained for the

Figure 6. Fluctuation spectrum ⟨|h(q)|2⟩ of the interface height h(x) for felodipine FI and FII (a) and nifedipine Nα and Nβ (b). The solid lines represent fits to the simulation results using eq 1. The fit is limited to the small q range. It corresponds to the domain where eq 1 is valid and thus in which the plot ln(⟨|h(q)|2⟩) vs ln(q) has a slope with the value −2. For larger values of q, this equation is no longer valid since the wavelength becomes commensurate with the crystal lattice spacing. For larger values of q, the straight line thus deviates significantly with data points.

Figure 7. Difference of the free energy between the liquid state and the crystal state as predicted by the Hoffman equation using Tm and ΔHm determined from MD simulations for the different investigated nifedipine and felodipine polymorphs. G

DOI: 10.1021/acs.jpcb.5b05557 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B different systems and their polymorphs.4,5,16 It is a direct consequence of the fair agreement obtained for ΔHm and Tm compared to the experimental data as reported in Table 1. In Figure 7, one particularly may notice that we were able to reproduce the experimental tendencies of all of the systems from the most to the less stable following this order: Nα, Nβ, FI, and FII. The greater this number, the greater the tendency for crystallization confirms that the thermodynamic driving force for crystallization of nifedipine is greater than that for felodipine. However, this information would suggest that the most stable polymorph will crystallize which does not seem to be the case following the well-known Ostwald rule of stage.118 This issue will be discussed further in subsequent text. D and η of nifedipine and felodipine in the equilibrated liquid state as a function of temperature T are shown in Figure 2. Upon decreasing temperature, both transport properties show the classical deviation from an Arrhenius behavior which is a typical behavior of the so-called fragile glass-formers.17 A fair agreement is found with experimental data for viscosity.33 Very similar mobility for both nifedipine and felodipine liquids, as measured by values of D and η, are found at all temperatures covered by MD simulations, which is in line with data reported in the literature.16 This result is consistent with the same experimental glass transition temperature Tg ≈ 304−318 K that has been particularly found for both pharmaceuticals.33 However, in the high-temperature range investigated from MD simulations, molecular mobility of nifedipine is found to be a bit slower than that of felodipine. It is opposite to the trend observed experimentally at lower temperatures either from viscosity measurements or dielectric relaxation spectroscopy data for which felodipine relaxation times are about 1 order of magnitude longer than nifedipine ones.33,119 This disagreement between the simulation and the experiments can be a consequence of the inability of the GAFF force field to reproduce dynamics with sufficient accuracy. However, from the trend reported from experimental viscosity measurements (see Figure 2), one may reasonably well also speculate about a possible crossover between felodipine and nifedipine dynamics that may occur from high to low temperatures. In the simulated high-temperature range, both diffusion coefficient and shear viscosity follow similar temperature dependence (see Figure 2). This link between these two transport properties is welldescribed by the Stokes−Einstein relation D ∼ η/T,120 which was shown to be perfectly valid (data not show). Inspired by this similarity and in order to fit the diffusion coefficient as a function of temperature, an approach similar to that described in ref 121 for viscosity was used. D has been thus fitted using the following equation: ⎡A ⎤ A ⎛1⎞ ln⎜ ⎟ = A1 + 2 exp⎢ 3 ⎥ ⎝D⎠ ⎣T ⎦ T

the nucleation/crystallization rates, the interplay between the driving force and the crystal−liquid interfacial free energy should be taken into consideration as it will be described in the following. In contrast to the case of the liquid−vapor interface, direct experimental measurements of the crystal−liquid interfacial free energy for organic materials are generally difficult.57,122 It is commonly determined at the melting temperature and taken as a fit parameter assuming that the homogeneous N, η, and ΔGv are known.50 It is also treated, to a first approximation, as a temperature-independent quantity.44,50 Therefore, only a few data of the crystal−liquid interfacial free energies are available in the literature. Most reported values are for metallic systems.50 Data for molecular compounds are very scarce, and only a few values are known.37,47−49 The paucity of the experimental data especially for pharmaceuticals has been an additional motivation to the present computational work. Using the fluctuation method, the crystal−liquid interfacial free energies at the melting temperature γm = 21.5, 14.4, 28.7, and 15.5 mJ/m2 have been determined for Nα, Nβ, FI, and FII, respectively. Uncertainties are of the order of 7−8%. Results including uncertainties are also reported in Table 1. For both nifedipine and felodipine compounds, the interfacial free energies are clearly the highest for the more stable polymorph (see phase diagram displayed in Figure 7). Indeed, since, at the crystal−liquid interface, the surface of a less stable phase is likely to be more disordered than the surface of a most stable one, hence the interfacial free energy is likely to be smaller. This trend is also well in line with the well-known Ostwald rule of stage118 suggesting that the crystal phase that nucleates is not the most thermodynamically stable phase but rather another metastable phase that is closest in Gibbs free energy to the parent phase. Moreover, in Table 1 one may notice that the interfacial free energy of the metastable Nβ is significantly smaller than the values obtained for all felodipine polymorphs. Such a result appears very reasonable in explaining why nifedipine may crystallize more easily than felodipine as reported in many investigations.16 We were unfortunately unable to compare our simulated values to experimental ones since no data are available for nifedipine and felodipine. However, the present obtained values are of the same order of reported data on other similar molecular compounds.49,123 For example, for indomethacin polymorphs, Andronis et al.37 estimated using nucleation rates and the CNT predictions, a crystal−liquid interfacial free energy of 27 and 17 mJ/m2 for the stable γ form and the metastable α form, respectively. They also found that the crystal−liquid interfacial free energy follows the Ostwald rule of stage: γ is higher for the stable polymorph than the metastable one. It can be also noted in the present numerical study that the trend of the crystal−liquid interfacial free energies between the different polymorphs is consistent with the melting entropies ΔSm = ΔHm/Tm. In other words, for the metastable polymorph, one finds a smaller value of the interfacial free energy and a smaller value of the melting entropy. The same trend is observed for felodipine (FI, FII) and indomethacin (α, β) in experiments.4,37 To the best knowledge of the authors, no other values were determined for pharmaceuticals. Few values have been reported for molecular compounds such as 7.0 mJ/m2 for succinonitrile,48 7.9 mJ/m2 for neopentylglycol,49 2.8 mJ/m2 for pivalic acid,124 4.6 mJ/m2 for cyclohexane,125 and 10 mJ/m2 for carbon tetrabromide.57 B. Prediction of the Nucleation Tendencies. The main parameters involved in nucleation as described by the CNT, i.e.,

(4)

where A1, A2, and A3 are adjustable parameters. It has been suggested that this alternative relation to the usual Vogel− Tamman−Fulcher (VFT) equation was particularly able to reproduce very well both high- and low-temperature ranges.120 It thus provides a possibility to estimate diffusivity in the lowtemperature range below the melting temperature required to estimate the kinetic prefactor involved in the CNT. Results of the fitting procedure are shown in the inset of Figure 2 for which very similar fitting parameters are obtained for both compounds (see caption Figure 2). Hence, since the molecular mobility of both compounds is very similar, in order to predict H

DOI: 10.1021/acs.jpcb.5b05557 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B ΔGv, γ, D, and liquid and crystal densities (ρlq and ρcx), have been calculated from the present MD simulations. The steadystate nucleation rate N can be thus estimated from the expression: ⎡ ΔG* ⎤ N = AN (T ) exp⎢ − ⎥ ⎣ kBT ⎦

proposed.87,88 For a hard-sphere model, a proportionality with the melting temperature γm ∝ kBTmρcx2/3 has been suggested which is consistent with the Turnbull relation since for all fcc metals the ratio ΔHm/kBTm is roughly the same.130 In Table 1, both the melting enthalpy and the crystal−liquid interfacial free energy of nifedipine and felodipine polymorphs are reported. At first glance, by comparing the melting enthalpies and the interfacial free energies of both compounds, it seems difficult to reconcile the higher values of the melting enthalpy found for nifedipine with their lower interfacial free energies. Polymorphs of the same compound show similar melting enthalpies but quite different interfacial free energies although the trend seems in agreement with Turnbull law considering a similar density for both polymorphs. In other words, it thus seems that the Turnbull relation may not be applicable in our case with the same CT. In our calculations, we have estimated CT values for each compound to be 0.28, 0.40, 0.44, and 0.75 for Nβ, Nα, FII, and FI, respectively. A similar calculation made for the indomethacin polymorphs using data taken from ref 37 would give us CT values of about 0.43 for the α polymorph and 0.59 for the γ polymorph close to the values obtained in the present study. Several works have mentioned the possible temperature dependence131 of γ, and the use of an effective temperaturedependent interfacial free energy has been shown from MD simulations to reproduce nucleation rates in the framework of the CNT.87,88,132−136 There are good reasons, related to the entropy loss due to the ordering of the liquid near the interface, suggesting that the interfacial free energy could increase with T.7 In the following, we have assumed the arbitrary form suggested in ref 88 which has been inspired from the Turnbull law for the temperature dependence of the interfacial free energy:

(5)

where AN(T) is the kinetic prefactor and ΔG* = (16π/3)[γ3/ ΔGv2] is the nucleation barrier. The prefactor is usually expressed as AN(T) = Znlq(24D(n*)2/3/λ2), where λ is the atomic jump distance approximated to 1/nlq1/3 in the following, n* is the number of atoms in the critical nucleus, nlq is the number density of liquid, Z = (ΔGv/6πkBTn*ncx)1/2 is the Zeldovich factor, and ncx is the number density of solid. We could write as well n* = (4/3)(π(r*)3ncx) and r* = 2γ/ΔGv, where r* is the critical radius of the nucleus.126 It should be noted that the validity of CNT is frequently questioned because of a lack of agreement between these predictions with measured nucleation rates either from experiments or from direct determination of these rates from simulation. Many checks of the CNT have been reported in the literature, ranging from fair agreement to severe disagreement.50 Nevertheless, because of the simplicity of its formalism, it is still commonly used to analyze crystal nucleation in experiments and simulations. The origin of the discrepancies has been mostly attributed to the kinetic prefactor rather than the nucleation barrier which is approximately correct in most cases.7 For certain glass-formers such as sodium metasilicate glasses127 differences by more than 10 orders of magnitude have been reported. The direct determination of the kinetic prefactor from MD simulations in the nucleation rate is found to be some 2 orders of magnitude larger than predicted by classical nucleation theory for a simple atomic system.128 It should be mentioned that disparities between predicted and fitted prefactors can be greatly reduced by using the crystal− liquid interfacial free energy as a fit parameter, i.e., a temperature-dependent parameter, or by changing the assumptions in using the viscosity instead of the diffusivity which is usually treated via the questionable validity of the Stokes−Einstein relation.41,46 A thorough determination of N requires a precise knowledge of the different physical parameters involved. In the present study, a direct determination of the liquid and crystal densities and D without the need to use the Stokes−Einstein relation has been obtained. ΔGv has been indirectly estimated from the Hoffman relation, and Tm and ΔHm have been directly determined. However, this latter relation has been generally shown to reproduce rather fairly the experimental data, especially in for felodipine and nifedipine.16 The remaining parameter is γ which has been directly determined at the melting temperature Tm. For a long time, it has been well-known that the interfacial free energy is correlated with both ΔHm and ρcx. This correlation was first established by Turnbull44 on the basis of experimental data of metallic systems. It is usually expressed as γm = CTΔHmρcx2/3, where CT is the Turnbull coefficient which was found to be about 0.45 for most metals and 0.32 for semimetals and water.48 Other studies showed that this coefficient, for metallic alloys, ranges from 0.13 to 1129 or from 0.21 to 0.7749 for a variety of molecular materials (camphene, benzene, lauric acid, stearic acid, and dibromobenzene, etc.).49 Therefore, the value of CT seems to change with the type of materials. Alternative rules to the Turnbull relation have been also

⎡ ρ (T ) ⎤2/3⎡ ΔH(T ) ⎤ ⎥ ⎢ γ(T ) = γm⎢ cx ⎥ ⎢⎣ ρcx (Tm) ⎥⎦ ⎣ ΔHm ⎦

(6)

where ΔH(T) is the difference of the enthalpy between the crystal and the liquid as a function of temperature. Using this form, the parameter CT is obviously different for each investigated compound. Figure 8 shows the steady-state N obtained from the CNT as a function of temperature for the different investigated polymorphs. Nucleation rates are represented as a function of the undercooling ΔT = Tm − T in order to facilitate the comparison between the different polymorphs. This representation also allows the correction from the slight disagreement between the experimental and the numerical melting temperatures (see Table 1). All nucleation curves show the same expected behavior. N can be first described by a sharp increase upon cooling due to the increasing driving force and then a decrease upon further cooling to low temperatures due to the decreasing molecular mobility. All curves thus exhibit a clear maximum. The value of this maximum significantly increases in the following order: FI, FII, Nα, and Nβ. The present results thus express the more favorable nucleation tendency of nifedipine with respect to felopidine as observed experimentally.16 Figure 8 also clearly suggests the predominant nucleation of the metastable Nβ polymorph as observed in spontaneous crystallization5 and well in line with the Ostwald rule of stage. The fact that nifedipine crystallizes more readily than felodipine was mainly attributed to the higher driving force for crystallization for nifedipine than felodipine since both I

DOI: 10.1021/acs.jpcb.5b05557 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 8. Predicted evolution of the nucleation N (a) and the growth G (b) rates as a function of temperature of the different investigated nifedipine and felodipine polymorphs. A comparison between experimental and numerical data is also given in the inset in b. Constant parameters 1/k = 60, 570, and 39 have been used to rescale the numerical growth rate in order to fit with experimental values using a growth model with f = 1 (see eq 7). The choice of the value of k does not change the order of the growth rate between the different polymorphs. Experimental growth rates for nifedipine Nβ and felodipine FI and FII were extracted respectively from refs 94 and 15.

molecules have a very similar mobility and glass transition temperature.16 The same argument cannot be applied when comparing polymorphs of the same molecule since, for monotropic systems, the driving force of the more stable phase is higher than the metastable phase for all temperatures. Our results show clearly an inverse correlation between the increasing value of γ and the decreasing maximum nucleation rate (see Figure 8 and Table 1) for each polymorph. The order observed for the maxima of nucleation (FI, FII, Nα, and Nβ) can be mainly explained by the influence of the interfacial free energy that acts as a barrier for the nucleation. However, the nucleation tendencies are not only controlled by the interfacial free energy since γm of Nα is higher than γm of FII but the nucleation rate of Nα is much higher. In that case, the higher driving force of Nα should explain its higher nucleation rate. For other compounds having approximately the same mobility but different crystallization rates18,137 as felodipine and nifedipine, it could be suggested that crystallization is controlled by the interplay between ΔGv and γ. This result is well in line with the nucleation barrier (ΔG* ∼ (γ3/ΔGv2) as predicted by the CNT. Despite a reasonable trend in the nucleation order, it can be also noted that the nucleation rate for the different polymorphs predicted from simulations deviates by several orders of magnitude from the experimental ones16 (see Figure 9). Moreover, the predicted nucleation rates for the different polymorphs are also separated by several orders of magnitude which does not seem to be the case experimentally as reported in ref 16. The present absolute value of the predicted nucleation rates should be thus considered with caution. An ideal situation certainly far from the reality is considered in the framework of the CNT where homogeneous bulk nucleation is assumed. Moreover, a direct comparison between numerical and experimental values is also difficult since it is not clear which polymorph is nucleated in the experimental study.16 As mentioned previously, significant deviations between the CNT

Figure 9. Comparison of the predicted evolution of the nucleation rate, N, as a function of temperature for two temperature-dependent models of the interfacial free energy, γ, of the different investigated nifedipine and felodipine polymorphs: constant value γ(T) = γm (a) and γ(T) = γm[ρcx(T)/ρcx(Tm)]2/3[ΔH(T)/ΔHm] (b). Experimental data extracted from ref 16 are also given for comparison. For each polymorph, either a dotted or dotted−dashed line indicates the prediction of N using values of the interfacial free energy (γ + Δγ) and (γ − Δγ). This calculation takes into account the uncertainties of the calculated γ of the order of Δγ/γm ≈ 7−8% (see Table 1). Using the γ(T) = γm model, the nucleation rate of felopidine FI polymorph is so low that it is not represented on panel a.

predictions and the experimentally measured nucleation rate are expected.50,127 In the present investigation, the nucleation rate is mostly influenced by the value of γ. In contrast, we checked that a small change in D or ΔG has much less influence on the nucleation rate than γ. As reported in ref 136, a small change of a few percent in the value of γ can significantly alter by several orders of magnitude the nucleation rates. In order to evaluate the sensitivity of γ on N, we have recalculated the nucleation rate by taking into consideration the Δγ/γm ≈ 7−8% uncertainties obtained on the value of γm determined in this work. In other words, we have computed the nucleation rates using either an estimated maximum or minimum value of γ using (γ + Δγ) and (γ − Δγ), respectively. We have also compared the nucleation rate previously determined with the mathematical form given in eq 6 with a behavior using a constant value γ(T) = γm owing to the same uncertainties. Parts a and b of Figure 9 show N obtained using the two different models of γ(T). The general trend of the nucleation rates is found to be similar using both models following the order of decreasing maximum of the nucleation rates: Nβ, Nα, FII, and FI. One may particularly notice the overall higher nucleation rate of the metastable β polymorph of J

DOI: 10.1021/acs.jpcb.5b05557 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

point along a direction normal to the surface. This behavior could be expected from a rough interface where on average all sites are equivalent.2,39,139 In that case, f ≈ 1 and growth at deep undercooling is mainly controlled by dynamics. In contrast, for perfectly smooth surfaces, growth requires the formation of two-dimensional nuclei on the interface. Such surface nucleation mediated growth proceeds via lateral aggregation. Models that are based on the two-dimensional nucleation are treated similarly to the three-dimensional theory.2,39,139 The function is given by f = exp[−(πaγ2/kBTΔGv)]. However, a crystal often possesses defects such as dislocations that may induce continuous steps on their surface. The dislocation thus also provides a perpetual ledge and there is no longer need for the nucleation of new layers as in the previous two-dimensional model. Steps form a spiral on the surface, and f is given by the number of growth sites where the spiral was developed on the crystal surface which is inversely proportional to the spacing between spiral steps.2,39,139 For dislocation-controlled growth, it can be shown that f = aΔGv/γ. It can be mentioned that it is believed that specific treatments to eliminate lattice defects are required in order to observe surface nucleation mediated growth as shown for molecular compounds such as salol.2,39 Experimentally, an adapted plot representing the reduced growth G/[AG(T)(1 − exp(−ΔG/kBT))] against undercooling ΔT = Tm − T may allow discrimination between these three models.2,37,39,96,138 In most cases, the inverse of the liquid viscosity 1/η(T) instead of diffusivity is supposed to provide a reliable description of the molecular mobility AG(T). For indomethacin, growth by two-dimensional nucleation has been found to be applicable37 from such an approach. However, recent works by Wu and Yu138 have also shown that viscosity may not accurately reproduce molecular mobility upon deep cooling owing to the breakdown of the Stokes−Einstein relation. Use of a fractional Stokes−Einstein relation46 where ηξ(ξ