Computing the Relative Nucleation Rate of Phenylbutazone and Sulfamerazine in Various Solvents Sharmistha Datta and David J. W. Grant* Department of Pharmaceutics, College of Pharmacy, University of Minnesota, Weaver-Densford Hall, 308 Harvard Street SE, Minneapolis, Minnesota 55455-0343
CRYSTAL GROWTH & DESIGN 2005 VOL. 5, NO. 4 1351-1357
Received December 4, 2003
ABSTRACT: This work simulates the nucleation of sulfamerazine in acetone, methanol, or water, and phenylbutazone in diethyl ether, methanol, or acetone. A new method is developed to estimate the time of onset of nucleation from supersaturated solutions using molecular dynamics simulation and knowledge of the solubility of the solute and the density of the saturated solution. Calculations are based on the dynamics of pairs of solute molecules at extreme supersaturations instead of dynamics of molecular aggregates at practical supersaturations. First, the characteristic radial distance (CRD) specific to each solute-solvent system is identified. Next, the time evolution of the radial distribution function (RDF), g(r), at CRD is evaluated. The onset of the nucleation is taken as the time after which g(r) at the CRD remains constant. Finally, the estimated relative nucleation times are compared with those measured experimentally. The calculated and experimental nucleation times in various solvents are related linearly and by rank order. However, the ratio of the calculated to the experimental nucleation time is on the order of 10-16 for sulfamerazine and 10-14 for phenylbutazone. This discrepancy arises from the underlying assumptions that are necessarily imposed by the limitation in the number of atoms the computation can treat. Introduction Crystallization of compounds from solvents is a common purification process. The physical and chemical properties of the solvent profoundly influence the crystallization kinetics of the solute compound and hence affect the physical properties, such as particle size, morphology, and polymorphic nature, of the crystallized material. The overall time of crystallization of a solute compound is essentially the sum of the nucleation time and the crystal growth time. Although the experimental measurements of the overall crystallization rate and the crystal growth rate are relatively simple, that of the nucleation rate may be quite difficult, primarily because the small size of a nucleus is not measurable by current experimental techniques. A widely used method to calculate the induction time for nucleation is to observe the change in the intensity of scattered laser light as nucleation begins. Other methods that have been utilized to measure nucleation rates are time-resolved simultaneous SAXS-WAXS experiments,1 laser polarized turbidimetry,2 refractive index measurements,3 and in situ atomic force microscopy.4 These experimental techniques indirectly indicate the beginning of nucleation. Currently, there is increasing interest in the prediction of nucleation rates and in the effect of additives on nucleation rates for better control of the final crystal properties.5 For example, Esselink et al.6 used molecular dynamics to study the nucleation and melting of n-alkanes, while Smith et al.7 used molecular mechanics to study the effect of acetals in inducing nucleation in polypropylene. * Corresponding author: David J. W. Grant, Department of Pharmaceutics, College of Pharmacy, University of Minnesota, WeaverDensford Hall, 308 Harvard Street SE, Minneapolis, MN 55455-0343, USA. Phone: (612) 624-3956. Fax: (612) 625-0609. E-mail: grant001@ umn.edu.
Molecular dynamics simulation is an important computational technique to understand the molecular mechanism of crystal nucleation and growth.8 However, the successful application of this method has primarily been limited to simulating crystallization from the melt because of the complexity of selecting suitable potential energy functions and because of the length of these calculations.9-15 Recently, Anwar and Boateng16 demonstrated that simulation of crystallization from solvents is feasible by performing molecular dynamics on a binary solution of model Lennard-Jones particles. Mucha and Jungwirth17 extended this idea by describing a method to simulate nucleation and crystal growth of sodium chloride from an aqueous solution using realistic interaction potentials. The selection of appropriate force field parameters for the particle system was perceived to be essential for successful simulation. In both of the above studies, the intermolecular interaction between the particles was of the Lennard-Jones type. The classical methods to identify crystallization using simulations are to note the absence of diffusion as well as the emergence of characteristic peaks in the radial distribution function (RDF), g(r), and the segmental orientation correlation, S(k).18 The RDF function, g(r), is used to determine how atoms organize themselves around one another and to reveal overall structural properties, such as packing, ordering, compressibility, and phase transitions. At the beginning of a simulation, when the solute and solvent molecules are randomly oriented, only short-range order is observed, as evidenced by peaks at short radial distances in the RDF plot. The onset of crystallization is indicated by the emergence of peaks at longer radial distances in the RDF plot, indicating long-range order. Molecular dynamics simulations have been utilized to simulate and to compare the effects of physical factors, such as temperature and supersaturation, on the nucleation rate.16 In this paper, we describe a
10.1021/cg0342462 CCC: $30.25 © 2005 American Chemical Society Published on Web 05/28/2005
1352
Crystal Growth & Design, Vol. 5, No. 4, 2005
method to estimate the relative nucleation rate of two pharmaceutically relevant compounds, phenylbutazone and sulfamerazine, in selected solvents using molecular dynamics simulations and knowledge of the solubility of the solute and the density of the saturated solution. We further compare the calculated nucleation rates with the experimental nucleation rates of these compounds in the specified solvents. Although the nucleation rates obtained from experimental and computational approaches differ by many orders of magnitude, the computational method may be used to predict the rank order of nucleation rates of solutes from solutions. Both sulfamerazine and phenylbutazone exhibit polymorphism, which is defined as the ability of a substance to crystallize in more than one crystalline form resulting from different arrangements and/or conformations of the solute molecules in the crystal lattice. For practical convenience, the solubility data necessary for calculating the solvent/solute ratios required for the molecular dynamics simulation refer to that of the most stable polymorph. For the calculations in the present work, molecular pairs of solute, which cannot represent a crystal lattice, are assumed instead of molecular aggregates. Therefore, this method cannot distinguish between polymorphs in the nucleation process. Finally, the extreme supersaturations and high pressure, which are necessary to increase the speed of the simulation, result in nucleation times that are many orders of magnitude shorter than the actual nucleation times. Sulfamerazine may crystallize in polymorphic forms I and II. Polymorph I, which is metastable at 25 °C, crystallizes before the stable polymorph II in all the solvents studied in this paper. During the solventmediated polymorphic transformation of I to II, the nucleation rate of II depends on the chemical properties of the solvent. Hence, the relative nucleation rates of sulfamerazine in the different solvents, predicted using molecular dynamics simulation, are compared with the nucleation rate of sulfamerazine II. Although phenylbutazone can crystallize in as many as five polymorphic forms depending upon the crystallization technique, isothermal crystallization of a relatively high supersaturated solution in any solvent at 25 °C yields only the stable δ polymorph. Hence, the measured nucleation rate of phenylbutazone is assumed to equal the nucleation rate of the δ polymorph. Materials and Methods Materials. Sulfamerazine (4-amino-N-[4-methyl-2-pyrimidinyl]benzenesulfonamide, Figure 1a, purity > 99.9%) and phenylbutazone (1,2-diphenyl-4-n-butyl-3,5-pyrazolidinedione, Figure 1b, purity > 99.9%) were purchased from Sigma Co., St. Louis, MO. Methanol, acetone, and diethyl ether were obtained from Fisher Scientific, Fairlawn, NJ. Residual water in the organic solvent was minimized by adding activated molecular sieves (Arcos Organics, Morris Plains, NJ). Methods. Measurement of Solubility and Density. The solubility of sulfamerazine polymorph II was determined in acetonitrile, methanol, and water by the equilibrium analytical method (Table 1). Excess solid sulfamerazine (form II) was suspended in 5 mL of each solvent, and the suspension was equilibrated by shaking at 100 strokes/min in a water bath (BT-47, Yamato Scientific Co. Ltd., Tokyo, Japan) maintained at 25 °C. After 5 days, the concentration of sulfamerazine in each of the solvents was determined by UV spectrophotometry (DU-7400,
Datta and Grant
Figure 1. Molecular structure of (a) sulfamerazine and (b) phenylbutazone. Table 1. Solubility and Density of Sulfamerazine in Acetonitrile, Methanol, or Water, and of Phenylbutazone in Diethyl Ether, Methanol, or Acetone at 25 °C
solvent acetonitrile methanol water diethyl ether methanol acetone
solubility (mM)
density of saturated solution (g/mL)
Sulfamerazine 16.3 15.2 1.23
0.790 0.800 0.997
Phenylbutazone 0.259 0.174 1.11
0.703 0.844 0.900
Beckman Instruments Inc., Fullerton, CA) at λ ) 307 nm, 307 ) 10 758 M-1 cm-1. The solubility data agree closely with those reported by Gu et al.19 Similarly, the solubility of phenylbutazone was determined in diethyl ether, methanol, and acetone by suspending excess solid phenylbutazone in 5 mL of each solvent and equilibrating the solution by shaking the suspension at 30 strokes/min in the water bath maintained at 25 °C. After 72 h, the concentration of phenylbutazone in each of the solvents was determined by UV spectrophotometry at λ ) 243 nm, 243 ) 15 084 M-1 cm-1. In both cases, equilibrium was judged to have been achieved when two successive absorbance readings did not differ by more than 2%. The density of the saturated solution of sulfamerazine (form II) or of phenylbutazone (form δ) in the various solvents (Table 1) was determined by weighing a known volume of the saturated solution. Experimental Measurement of Nucleation Rates. The nucleation times of sulfamerazine (form II) are described by Gu et al.19 In this case, supersaturation of polymorph II was achieved by preparing a saturated solution of polymorph I, which undergoes solvent-mediated polymorphic transformation to polymorph II. The nucleation times of phenylbutazone (form δ) were determined as follows. A defined weight of phenylbutazone was added to 5 mL each of methanol, acetone, or diethyl ether to prepare supersaturated solutions in which the ratio of concentration of phenylbutazone in the solvent, cs, to that in the saturated solvent, ceq, is 2.5 at 25 °C. The solutions were heated to 50 °C to dissolve the excess solid. After the last
Relative Nucleation Rate of Phenylbutazone and Sulfamerazine particle of phenylbutazone had dissolved during heating, the solution was maintained at 50 °C for 2 min to dissolve any small particles, including nuclei, which could not be seen by the naked eye. Furthermore, the isothermal heating for 2 min ensures that the different saturated solutions have the same thermal history. The supersaturated solutions were then placed in a water bath at 25 °C, and the time at which the first crystals appeared to the naked eye was taken as the induction time for nucleation. The heat capacity and thermal conductivity of the solvent influence the cooling rate of the supersaturated solution from 50 to 25 °C. However, in all cases, the induction time for nucleation was greater than the time required for the temperature to decrease to 25 °C. Hence, the initial change in temperature due to cooling from 50 to 25 °C was assumed not to affect the relative nucleation rates of phenylbutazone in different solvents. Any crystal that can be seen by the naked eye is considerably larger than a nucleus. Therefore, the measured induction time is actually the sum of the nucleation time and the crystal growth time to reach visibility. The crystal growth rate of phenylbutazone was determined by measuring the linear growth rate of the crystals along the needle axis. Because the crystal growth rate of phenylbutazone determined in the different solvents is similar, any difference in the induction time results from the difference in nucleation time in different solvents. Computational Estimation of Nucleation Rates. All the simulations and calculations were conducted using Cerius2 software (Accelrys Inc., San Diego, CA). The relative nucleation rate of the model compounds, phenylbutazone and sulfamerazine, in the various solvents was calculated in several steps using molecular dynamics, as exemplified by sulfamerazine in acetonitrile. (1) The detailed molecular structure of sulfamerazine was constructed using the 3-D Builder. The hybridization state of each of the atoms was checked and corrected, if necessary. The charge of each atom was calculated using the charge equilibration method developed by Rappe and Goddard,20 with convergence criterion of < 5 × 10-4 Qeq. The energy expression was set up using Drieding 2.21 force field.21 The atom types and bond types were assigned by following the rules of the selected force field. The energy of the molecule was minimized by Smart Minimizer, which uses the iterative Newton-Raphson (NR) method to minimize the energy function. The maximum atomic displacement per step was 2 Å, with an adopted-basis NR termination criterion of 50 kcal/mol and quasi-NR termination criterion of 0.1 kcal/mol. (2) The possible torsion angles of the sulfamerazine molecule selected for variation were C4-S1-N2-C7 and C3-C4-S1N2 (Figure 1a). For phenylbutazone, the torsion angles to be varied were N2-N1-C8-C9, C1-C2-C4-C5 and C2-C4C5-C6 (Figure 1b). In each case, the other torsion angles of the solute molecules and all those of the solvent molecule were held constant. (3) Using Amorphous Builder, a three-dimensional box was constructed, which contained two molecules of sulfamerazine and n molecules of acetonitrile. The number n was calculated from the solubility. For example, the solubility of sulfamerazine in acetonitrile is 0.016 M. Therefore, in a saturated solution, two molecules of sulfamerazine require 2400 molecules of acetonitrile. However, the time required for a dynamics calculation on more than 500 molecules is excessive. Furthermore, primary nucleation occurs only in a supersaturated solution. To address these two issues, the solution concentration for the dynamics calculation was 40 times the solubility. Consequently, for the dynamics calculation, the concentration of sulfamerazine was selected such that two molecules of sulfamerazine were present with 60 molecules of acetonitrile. Hence, for the sulfamerazine-acetonitrile system, the value of n was 60. Figure 2 shows the simulation box containing two molecules of sulfamerazine and 60 molecules of acetonitrile. The above method could not be applied to sulfamerazine in water because the solubility is so low that the number of water
Crystal Growth & Design, Vol. 5, No. 4, 2005 1353
Figure 2. The three-dimensional box containing 60 molecules of acetonitrile and two molecules of sulfamerazine. The dimension of the box is determined from the density of the saturated solution of sulfamerazine in acetonitrile, as described in the text.
Figure 3. The three-dimensional box containing 100 molecules of methanol and two molecules of phenylbutazone. The dimension of the box is determined from the density of the saturated solution of phenylbutazone in methanol, as described in the text. molecules necessary for the dynamics calculation, described above, is 2600. Because this number is too large for practical computation, we selected 100 molecules of water with two molecules of sulfamerazine. When the solute was phenylbutazone, the solution concentration for dynamics calculation was eight times the solubility. For the dynamics simulation, the lower supersaturation of phenylbutazone solutions than that of sulfamerazine solutions was possible because the solubilities of phenylbutazone in the studied solvents are higher than those of sulfamerazine in the same solvents. Figure 3 shows the simulation box containing two molecules of phenylbutazone and 100 molecules of methanol. The rationale for selecting two molecules of sulfamerazine is described in step 5 of this section. (4) Molecular dynamics was performed on the aforementioned box. The canonical function (NPT) was defined by constant values of the total number of molecules, N, the pressure, P, and the temperature, T, of the system. P was set at the default value, 1 GPa, while T was set at 298 K. The T-damping thermostat function was used to maintain temperature equilibrium during the dynamics calculation with 0.1
1354
Crystal Growth & Design, Vol. 5, No. 4, 2005
Datta and Grant
Figure 4. Representative plot of RDF against distance, r, obtained for the sulfamerazine-acetonitrile system after 11 ps of dynamics time. picosecond (ps) relaxation time. A total of 10 000 dynamic steps were performed in which each step lasted for 1 femtosecond (fs). The trajectory of the system upon dynamics simulations was determined by saving a trajectory frame for every 10 dynamic steps, i.e., after every 0.01 ps of dynamics time. Thus, a total of 1000 trajectory frames were saved. (5) The two sulfur atoms, S1 and S2, belonging to two different sulfamerazine molecules were placed into two separate groups, termed A and B, by individual selection. Then, the unnormalized distribution function DAB(r) was defined by eq 1, where rA and rB are the positions of sulfur atom 1, S1, and sulfur atom 2, S2.
DAB(r) ) 〈
∑ δ(r - |r
A
- rB|)〉
(1)
Figure 5. (a) Time evolution of the RDF, (a) averaged over 0.5 ps, g(r)0.5, and (b) over the total time period, g(r)full, as a function of distance, r, for the sulfamerazine-acetonitrile system.
The RDF between groups A and B, gAB(r), were calculated from DAB(r) using eq 2, where NA and NB are the number of atoms in group A and group B respectively, NAB is the number of atoms common to both groups A and B, and V is the volume of the three-dimensional box. The term ∆r was evaluated by dividing the cutoff radial distance, rcutoff, by the number of concentric bins, Nbins.
gAB(r) )
DAB(r) V (NANB - NAB)4πr2∆r
(2)
The function, gAB(r), is the spherically averaged distribution of interatomic vector lengths, r. For sulfamerazine or phenylbutazone, A and B are atoms of the same element, so gAB(r) was replaced by g(r) for simplicity. Furthermore, because only one atom is present in groups A and B, NA and NB are both unity. Moreover, because no atom is present in both groups, NAB is zero. Therefore, eq 2 can be replaced by eq 3.
g(r) )
D(r)V 4πr2∆r
(3)
Figure 4 shows a representative plot of g(r) after 11 ps of dynamics time against the radial distance (r, Å). The RDF can be calculated for every single time frame or can be averaged over multiple time frames. Statistically, the averaged RDF is more reliable than that obtained from a single frame. Hence, for preliminary examination, the RDF averaged over 0.5 ps, g(r)0.5, was used for further analysis. Figure 5a shows the time evolution of g(r)0.5. The radial distance at which g(r)0.5 is maximal is termed the characteristic radial distance (CRD), which initially increases with time and then becomes constant. For the sulfamerazine-acetonitrile system, the CRD at the beginning of the dynamics simulation is 12.55 Å and at the end of the simulation is 12.83 Å. Although the time evolution of the CRD is determinable from g(r)0.5, the evolution is easier to visualize when the RDF is averaged over the entire time frame, g(r)full, as shown in Figure 5b. However, for comparing g(r) at different simulation times, the number of frames over which RDF is averaged to calculate g(r)full is not constant with time. Because the above advantage out-
Figure 6. RDF, g(r)full, at a distance of CRDinitial ) 12.55 Å and CRDfinal ) 12.85 Å as a function of time for the sulfamerazine-acetonitrile system. weighs the disadvantage, g(r)full was used for determining the time when displacement of CRD is maximized in subsequent analysis. Similarly, to calculate g(r), and hence g(r)full, of phenylbutazone, the two nitrogen atoms, N1 and N2, belonging to two different phenylbutazone molecules were selected. (6) Figure 6 shows time plots of the numerical value of g(r)full at the initial CRD (12.55 Å) and at the final CRD (12.83 Å). The time at which g(r)full becomes constant was taken as the estimated nucleation time, which is 5 ps for sulfamerazine in acetonitrile. The nucleation time of phenylbutazone and sulfamerazine in other solvents was estimated by the same procedure using the parameters stated in column 3 and 4 of Table 2. Columns 5 and 6 of Table 2 show the estimated nucleation time, τcalc, and the experimental nucleation time, τexp, for phenylbutazone and sulfamerazine, respectively, in the various solvents.
Results and Discussion The total energy of a solute-solvent system is the sum of the various intramolecular and intermolecular interactions among the various molecules constituting the system. Hence, to calculate the potential energy of such a system, suitable parameters for both the solute and the solvent molecules are required. In the present
Relative Nucleation Rate of Phenylbutazone and Sulfamerazine
Crystal Growth & Design, Vol. 5, No. 4, 2005 1355
Figure 7. (a) RDF, g(r), at the respective characteristic radial distance, CRDfinal, plotted against time, when sulfamerazine is dissolved in (a) acetonitrile (CRDfinal ) 12.83 Å), (b) methanol (CRDfinal ) 14.83 Å), or (c) water (CRDfinal ) 6.31 Å), and when phenylbutazone is dissolved in (d) diethyl ether (CRDfinal ) 13.95 Å), (e) methanol (CRDfinal ) 10.51 Å), or (f) acetone (CRDfinal ) 7.53 Å). Table 2. Number of Solvent Molecules in the Dynamic Box, n, Final Characteristic Radial Distance, CRDfinal, Estimated Nucleation Time, τcalc, and Experimental Nucleation Time, τexp, Described in the Text for Sulfamerazine and Phenylbutazone Molecules in Various Solvents solute sulfamerazine phenylbutazone
solvent
n
CRDfinal (Å)
τcalc (ps)
τexp (min)
τcalc/τexp
acetonitrile methanol water diethyl ether methanol acetone
60 80 100 60 100 20
12.83 14.89 6.31 13.95 10.51 7.53
5 40 ∼120 7.5 13 20
120 7200 >21600 10 45 75
7.00 × 10-16 9.26 × 10-17 ∼9.26 × 10-17 1.25 × 10-14 4.83 × 10-15 4.44 × 10-15
solute-solvent systems, both the solute and the solvent molecules consist of atoms of different types and hybridizations. Force fields were used to generate explicit parameters for each molecule using selected atomic parameters coupled with rules derived theoretically and empirically. Hence, the need to calculate the parameters of individual atoms was avoided, facilitating quicker molecular dynamics calculations. Furthermore, the results of the molecular dynamics calculations for different solvents can be compared when a single force field is used for their parametrization. In this paper, the parameters used to calculate the intramolecular and intermolecular interactions of the solute and solvent molecules were derived from Drieding 2.21 force field,20 which is a purely diagonal force field with harmonic valence terms and a cosine-Fourier expansion torsion term. The umbrella function form was used for inversions, which are defined according to the Wilson outof-plane definition. The van der Waals interactions were described by the Lennard-Jones potential. Electrostatic interactions were described by the atomic monopoles and a screened Coulombic term. Hydrogen bonding was described by an explicit Lennard-Jones 12-10 potential. It is recognized that the force constants and the geometric parameters for the Drieding 2.21 force field are based on simple hybridization rules rather than on specific combinations of atoms. Hence, while the use of Drieding 2.21 force field simplifies the dynamics calculation, more realistic simulations could be performed, if the parameters of each system could be obtained separately. Estimation of the absolute values of the nucleation rates is not practically feasible. Molecular dynamics
simulations are performed on systems that have submicroscopic dimensions. Hence, these calculations indicate the changes at the submicroscopic levels that may exist for only a short time. Extending the calculations to larger systems with microscopic or greater dimensions would simulate a real system. However, the time required for simulations on such systems would be enormous. The real system is assumed to be an ensemble of the system on which molecular dynamics simulation is performed. Furthermore, the extreme supersaturations and high pressure assumed for the purposes of calculation are not practically feasible. Three assumptions are used to increase the speed of simulation. Hence, results from a molecular dynamics calculation serve only as indicators of a real system. The absolute time scales in a simulation study may not be comparable with real time scales. The short time scale required for molecular dynamics calculations also suggests that the method does not require excessive computational time and resources. This method therefore rapidly predicts the relative nucleation rates in different solvents. Figure 7a-f shows the time evolution of the radial distribution function, g(r), at final CRD (CRDfinal) of phenylbutazone and sulfamerazine in the different solvents. Table 2 lists the estimated nucleation times of phenylbutazone and sulfamerazine obtained from these plots. For sulfamerazine, the estimated nucleation time bears a similar ratio (∼10-16) to the experimental nucleation time in the three solvents. The different ratio for acetonitrile may be a computational artifact. When molecular dynamics simulation is begun, the contents of the system are oriented in random order. At least a
1356
Crystal Growth & Design, Vol. 5, No. 4, 2005
Figure 8. Correlation between the estimated nucleation time, τcalc, by the method described in the text and the experimental nucleation time, τexp, of (a) sulfamerazine in acetonitrile, methanol, and water (y ) 0.323x + 3.104, R2 ) 0.9992); (b) phenylbutazone in diethyl ether, methanol, and acetone (y ) 0.1913x + 5.2087, R2 ) 0.9872).
few molecular dynamics steps are required to bring the energy of the system to a perceivable equilibrium value. Hence, it is advantageous to disregard at least the first 100 dynamic steps. However, when the nucleation time is relatively short, as occurs in acetonitrile, nucleation proceeds during the time of energy equilibration, which results in a higher ratio, τcalc/τexp. On the other hand, the system in which water is the solvent differs from the systems with other solvents, because the number of water molecules used for the dynamics calculation is much smaller than that required by the method specified. Furthermore, the Dreiding 2.21 force field may be less suitable for simulating water dynamics. The use of other specialized force fields or an entirely separate parametrization of the water molecule may improve the dynamics result. Despite these irregularities in the calculation, Figure 8a shows that the estimated nucleation times and the experimental values, for sulfamerazine in the three solvents, are linearly related with only a small intercept. For phenylbutazone (form δ), the ratio of the calculated nucleation time to the experimental nucleation time is similar for methanol and acetone (∼4.6 × 10-15 in Table 2). The different ratio for diethyl ether results from the same experimental artifact that affects the τcalc/ τexp ratio in the sulfamerazine-acetonitrile system. Because the estimated nucleation time with diethyl ether as the solvent is short, it is difficult to subtract the time required for energy equilibration. Despite this anomaly, Figure 8b shows that the estimated nucleation times and the experimental values for phenylbutazone in the three solvents are linearly related. To summarize, Table 2 and Figure 8 show that the relative nucleation rates estimated by the above method in different
Datta and Grant
solvents are related linearly and by rank order to the actual nucleation rates in these solvents. This paper essentially explains the methodology of the computational technique to estimate the relative nucleation rate of a solute in different solvents. Here, we do not explore all the possible applications of this method. At the present time, the method can be used to indicate the differences in nucleation rates in different solvents without actually suggesting a reason for the differences. A qualitative reason for the difference in nucleation rate may be obtained from the individual trajectory frames. For example, the phenylbutazone molecules in acetone are more constrained than those in methanol, as evidenced by the reduced movement of the phenylbutazone molecules in the acetone system in consecutive trajectory frames. This constraint may result from weak hydrogen bonding between the carbonyl oxygens of acetone and the methine hydrogen of phenylbutazone. Careful analysis of the trajectory frames may be used to determine the reason for the differences in the nucleation rate in quantitative terms. While this present method can provide preliminary predictions of nucleation rates, some limitations need to be considered. The present method assumes that the nucleation rate of a solute is faster when the time required to reach equilibrium between the two solute molecules is shorter. While increasing the number of solute molecules may enhance the credibility of the calculations, the enormous increase in the computation time prohibits such practice. In fact, the computational resources required for the calculations increase exponentially as the number of molecules, on which dynamics in performed, increases. Hence, if the solubility of the solute is relatively low in a solvent and the density of the solution is relatively high, much time may be required for the completion of the computation, as explained for the sulfamerazine-water system. Furthermore, although the estimated nucleation time is in picoseconds in the examples presented in this paper, the estimated nucleation time can be much longer in other systems. Approximately 12 h were required for molecular dynamics simulation of each system. The use of Drieding 2.21 force field was found to be adequate for energy calculations in the present examples because the parametrization of the atom types seemed reasonable. However, for other systems the suitability of the force field should be checked before using the method. Finally, this method assumes that the solute and solvent molecules do not undergo any chemical reaction, which would need to be verified experimentally in each case. Conclusions We present a computation method to estimate the relative nucleation rates of phenylbutazone (form δ) and sulfamerazine (form II), each in three selected solvents. To increase the speed of calculation, extreme supersaturations, high pressure, and equilibrium between molecular pairs of solutes are assumed. For both the estimated and experimental nucleation times, the rank order for sulfamerazine in the studied solvents is acetonitrile < methanol < water, while that for phenylbutazone is diethyl ether < methanol < acetone. The ratio of the calculated nucleation time of sulfamerazine to the experimental nucleation time is on the order of
Relative Nucleation Rate of Phenylbutazone and Sulfamerazine
10-15 to 10-16, while that of phenylbutazone is on the order of 10-14 to 10-15. This discrepancy arises from a necessary limit on the number of atoms for the computation, which require restrictive underlying assumptions. The work illustrates a rapid method to predict the relative nucleation rates of solutes in solvents. Acknowledgment. We thank the United States Pharmacoepial (USP) Convention for the award of a USP Fellowship in Drug Standards Research to S.D. We also thank the Supercomputing Institute of the University of Minnesota for financially supporting our use of the Medicinal Chemistry/Supercomputing Institute Visualization - Workstation Laboratory. References (1) Ueno, S.; Ristic, R. I.; Higaki, K.; Sato, K. J. Phys. Chem. B. 2003, 107, 4927. (2) Martini, S.; Herrera, M. L.; Hartel, R. W. J. Agric. Food Chem. 2001, 49, 3223. (3) Devarakonda, S.; Evans, J. M. B.; Myerson, A. S. Cryst. Growth Des. 2003, 3, 741. (4) Beekmans, L. G. M.; Vallee, R.; Vancso, G. J. Macromolecules 2002, 35, 9383. (5) Nagarajan, K.; Myerson, A. S. Cryst. Growth Des. 2001, 1, 131.
Crystal Growth & Design, Vol. 5, No. 4, 2005 1357 (6) Esselink, K.; Hilbers, P. A. J.; van Beest, B. W. H. J. Chem. Phys. 1994, 101, 9033. (7) Smith, T. L.; Masilamani, D.; Bui, L. K.; Brambilla, R.; Khanna, Y. P.; Garbriel, K. A. J. Appl. Polym. Sci. 1994, 52, 591. (8) Allen, M.; Tildesley, D. Computer Simulation of Liquids; Oxford Science Publications: Oxford, 1987. (9) Boek, E. S.; Briels, W. J.; Feil, D. J. Phys. Chem. 1994, 98, 1674. (10) Mandell, M. J.; McTague, J. P.; Rahman, A. J. Chem. Phys., 1976, 64, 3699. (11) Mandell, M. J.; McTague, J. P.; Rahman, A. J. Chem. Phys. 1977, 66, 3070. (12) Swope, W. C.; Andersen, H. C.; Phys. Rev. B 1990, 41, 7042. (13) van Duijneveldt, J. S.; Frenkel, D. J. Chem. Phys. 1992, 96, 4655. (14) ten Wolde, P. R.; Ruiz-Montero, M. J.; Frenkel, D. J. Chem. Phys. 1996, 104, 9932. (15) Nose´, S.; Yonezawa, F. J. Chem. Phys. 1986, 84,1803. (16) Anwar, J.; Boateng, P. K. J. Am. Chem. Soc. 1998, 120, 9600. (17) Mucha, M.; Jungwirth, P. J. Phys. Chem. B 2003, 107, 8271. (18) Rigby, D.; Roe, R. J. J. Chem. Phys., 1988, 89, 5280. (19) Gu, C. H.; Young, V., Jr.; Grant, D. J. W. J. Pharm. Sci. 2001, 90, 1878. (20) Rappe A. K.; Goddard, W. A. J. Phys. Chem. 1991, 95, 3358. (21) Mayo, S. L.; Olafson, B. D.; Goddard, W. A. J. Phys. Chem. 1990, 94, 8897.
CG0342462
