Molecular Dynamics Simulation on Nucleation of Ammonium

The free energy of nucleation obtained from the cluster size distribution definitely shows the energy barrier ..... This material is available free of...
1 downloads 0 Views 4MB Size
Download PDF
Article pubs.acs.org/crystal

Molecular Dynamics Simulation on Nucleation of Ammonium Perchlorate from an Aqueous Solution Hong-Min Shim,† Jae-Kyeong Kim,† Hyoun-Soo Kim,‡ and Kee-Kahb Koo*,† †

Department of Chemical and Biomolecular Engineering, Sogang University, Seoul 121-742, Korea Agency for Defense Development, Daejeon 305-600, Korea

‡

S Supporting Information *

ABSTRACT: The simulation system that yields reliable physicochemical properties of ammonium perchlorate (AP) such as lattice energy, solution density, diffusion coefficient, and radial distribution functions was constructed. Dissolution of AP in water was computationally confirmed to be an endothermic reaction, which corresponds to the fact that crystallization of AP from an aqueous solution is an enthalpically favored process. On the basis of the Yasuoka−Matsumoto method, the MD simulation on the homogeneous nucleation of AP from an aqueous solution was carried out, and the calculated nucleation rate was shown to be 1 order of magnitude different from that predicted by the classical nucleation theory (CNT). The critical sizes of nucleus (12 and 8 at S = 5.1 and 6.4, respectively) are also very close to those estimated by the CNT (12.4 and 8.2 at S = 5.1 and 6.4, respectively). The free energy of nucleation obtained from the cluster size distribution definitely shows the energy barrier for nucleation, which is consistent with the trends predicted by the CNT. The time-dependent concentration and maximum cluster size obtained by the present simulation clearly show that the crystallization of AP can be divided into two stages, nucleation and nuclei growth. At the earlier stage a lot of nuclei are spontaneously formed, and at the stage of nuclei growth AP molecules tend to be incorporated into the existing nuclei and continuously merged into a larger one.

1. INTRODUCTION Nucleation, the initial step of crystallization, is of undisputed importance because the crystal structure of nuclei and nucleation rate affect the polymorphism, morphology, crystal size, and their distribution. At the atomic level, nucleation is elucidated by aggregation of solute molecules, and hence, it can be characterized by the collective behavior of solute molecules as a function of time. In this respect, the molecular dynamics (MD) simulations are regarded as a powerful tool and expected to provide a deeper understanding of nucleation. In recent years, the nucleation mechanism of methane hydrates that is difficult to access experimentally was successfully explained through the MD simulations by Jiménez-Á n geles and Firoozabadi.1 Furthermore, Salvalaglio et al. revealed a key role of additives on the surface nucleation by the MD simulation.2 The pioneering study on nucleation by the MD simulation was conducted by Anwar et al.3,4 They described nucleation as aggregation of simple particles induced by the Lennard-Jones potential. The distinctive ordering of particles during clustering represents intuitively a nucleation process. Moreover, the importance of Coulomb-type contribution on nucleation that the partial charge of a molecule increased aggregation rate was emphasized by Gavezzotti.5 Cheong and Boon studied the relationship of force fields and atomic partial charges available for promoting 2-D nucleation.6 Those studies provide a good explanation for the nucleation process in a limited system. © 2014 American Chemical Society

However, they are devoid of nucleation kinetics such as nucleation rate. The nucleation rate is another key factor characterizing the nucleation phenomenon because it is profoundly related with the supersaturation ratio, critical size of nucleus, and free energy change of nucleation. Valeriani et al. tried to calculate the nucleation rate of sodium chloride from its melt, where the relationship of free energy barrier and supersaturation is well established.7 Yasuoka and Matsumoto proposed a method for the calculation of the nucleation rate by analyzing change of clusters with time. They gave a good explanation on homogeneous nucleation from a vapor phase.8−10 However, those methods remain a challenge when applied to the nucleation from a liquid solution because solute molecules strongly interact with solvent molecules by electrostatic force as well as van der Waals force in a very short molecular distance. Therefore, the MD simulations often produce unrealistic results where dissolution of solute molecules from the crystal faces is observed even at high supersaturation.11 In the present work, we aim to simulate the homogeneous nucleation of ammonium perchlorate (AP), which is a kind of solid oxidizers used for the solid propellants. Although heterogeneous nucleation frequently occurs, a study on Received: July 24, 2014 Revised: October 2, 2014 Published: October 3, 2014 5897

dx.doi.org/10.1021/cg501112a | Cryst. Growth Des. 2014, 14, 5897−5903

Crystal Growth & Design

Article

as the value of a concentration of the system divided by an equilibrium concentration at simulation temperature of 273.15 K, where solubility is XAP = 0.0184. The supersaturation ratios given at 353.15 and 373.15 K were S = 5.1 and 6.4, respectively. Initially, each system was geometrically optimized and followed by an NPT simulation for 1 ns, and it was subsequently quenched to 273.15 K with an NPT simulation for 10 ns, where the pressure was also coupled to 1.013 bar. In the present work, the time step was 2 fs, and a cluster was defined as an aggregate of molecules within a radius of 0.21 nm. The solvation free energy of AP was calculated by the coupling parameter approach. In this approach, each state is a function of a coupling parameter, λ, the extent to which the interaction potential U is perturbed. In the present work, the full solute−solvent interaction is done when λ = 0, and solute does not interact with the solvent when λ = 1. In NPT ensemble, the free energy change can be determined by

homogeneous nucleation is fundamentally important for understanding phase transition. This article reveals the potential of a computational approach to the kinetics of nucleation of ionic compounds from an aqueous solution. The correctness of the system we construct is confirmed by comparing the calculated physicochemical properties such as lattice energy, solution density, diffusion coefficient, and radial distribution functions (RDFs) with the reported data. Moreover, the related energy and molecular configuration as a function of time are analyzed to investigate the characteristics of nucleation at nonequilibrium state. The nucleation rate, critical size of nucleus, and free energy change of nucleation are calculated by the Yasuoka−Matsumoto (Y−M) method. However, there are some limitations when quantitatively being compared with experimental results owing to various experimental difficulties. Therefore, the simulation results are compared with the theoretical prediction. Lastly, the distinctive steps observed in the MD simulations such as nucleation and nuclei growth are specifically tackled.

ΔG =

The force field used in the present work was CHARMM (Chemistry at HARvard Macromolecular Mechanics).12,13 Atomic partial charges of the NH4+ and ClO4− were obtained by using AM1 basis set (Materials studio, ver. 6.1) and LCAO MO SCF calculation by Wagner, respectively.14,15 That of a water molecule was obtained from the SPC/E model.16 Simulation details are provided by Supporting Information 1. The intramolecular motions of ions were quenched by applying restraints. The MD simulations were performed under periodic boundary conditions using Gromacs (ver. 4.5.5).17−20 The cutoff distances were 9 Å for electrostatic interaction and 12 Å for LennardJones interaction, respectively. Typically, the Lennard-Jones potential is truncated at a cutoff distance of rc = 2.5σ to save computational time. In the present work, the largest value of σ is 4 Å of chlorine atom, and hence the cutoff distance of 12 Å is sufficient for the MD simulation. Long-range electrostatic interactions were treated by using the Particle-Mesh Ewald (PME) summation method.21 Temperature and pressure coupling were achieved by using V-rescale and Parrinello− Rahman, respectively.22,23 For the calculation of equilibrium properties of AP such as solution density and diffusion coefficient, the initial system was constructed by satisfying the solubility of AP in water at 293.15 K as shown in Table 1.

a

Nammonium

Nperchlorate

XAP

293.15 353.15 373.15

16000 16000 16000

576 1656 2160

576 1656 2160

0.0347 0.0938 0.1192

Table 2. Solvation Free Energy and Solvation Enthalpy of NH4+ and ClO4− in Water ΔGsolv,293.15 K (kJ/mol)

ΔHsolv (kJ/mol)

−279.1 −312.4

−277.4 −312.1

−310.8 −317.9

(1)

(2)

where Ecryst is the total energy of a crystal that is normalized as a value per mol (Ecryst = −668.71 kJ/mol), Ecation and Eanion are each the energy of isolated cation and anion (Ecation = −28.26 kJ/mol and Eanion = −28.18 kJ/mol). The calculated lattice energy was −612.27 kJ/mol, which corresponds to the experimental lattice energy of −601.66 kJ/mol with the deviation of 1.7%.25 When nucleation from a liquid solution is simulated, solution density and diffusion coefficient are important factors because the former represents how closely the molecules exist in a given solution and the latter influences the rate of collision of solute molecules ensued by nucleation. Through analyzing MD simulation results by an NPT ensemble at 293.15 K and 1.013 bar, the solution density was calculated to be 1.10 g/cm3, which is consistent with the experimental density of 1.09 g/cm3.24 The diffusion coefficient of AP molecules was calculated by an increase of mean squared displacement (MSD) with time. The diffusion coefficient, D, is expressed as

After geometry optimization, an NPT simulation was carried out for 1 ns to reach the equilibrium state of a solution at the temperature of 293.15 K and the pressure of 1.013 bar. However, the supersaturated solutions were constructed by quenching the saturated solutions at 353.15 and 373.15 K to the temperature of 273.15 K, respectively (Table 2). In the present work, the supersaturation ratio was defined

ΔGsolv,278.15 K (kJ/mol)

dλ NPT, λ

E latt = Ecryst − Ecation − Eanion

N is the number of molecules, and X is the mole fraction of AP.24

NH4+ ClO4−

∂U ∂λ

3. RESULTS AND DISCUSSION 3.1. Properties of Equilibrium State. For the validation of the system, where NH4+, ClO4−, and H2O molecules are interacted with each other, the physicochemical properties of AP including lattice energy, solution density, diffusion coefficient, and RDFs were estimated. The lattice energy is fundamentally defined as the intermolecular interactions between solute molecules contained in a crystal. The unit cell of AP (the space group, Pna21; lattice parameters, a = 9.22 Å, b = 5.81 Å, c = 7.45 Å, Z = 4) was expanded to 4, 5, and 6 times in the x-, y-, and z-directions, respectively, and followed by energy minimization. The lattice energy, Elatt, was calculated by eq 2.

Table 1. Simulated Concentration of AP at Equilibriuma Nwater

1

In the present work, the electrostatic and Lennard-Jones interactions were separately calculated. For the electrostatic interactions that depend linearly on λ, the decoupling was performed at values of λ = 0.0, 0.2, 0.4, 0.6, 0.8, and 1.0. However, for the LennardJones interactions whose derivative does not change monotonically, a narrow value of λ spacing value was employed for the calculation: λ = 0.00, 0.05, 0.10, 0.15, 0.20, 0.25, 0.30, 0.35, 0.40, 0.45, 0.50, 0.55, 0.60, 0.65, 0.70, 0.75, 0.80, 0.85, 0.90, 0.95, and 1.00. Simulations were carried out at constant pressure of 1.013 bar and temperatures of 278.15 and 293.15 K. The resulting solvation free energies were used for calculating the enthalpy change of a solution of AP in water.

2. COMPUTATIONAL METHODS

temperature (K)

∫0

5898

dx.doi.org/10.1021/cg501112a | Cryst. Growth Des. 2014, 14, 5897−5903

Crystal Growth & Design D = lim

t →∞

Article

2 ⟨| r (⃗ t ) − r (0) ⃗ |⟩ 6t

3.2. Characteristics of Nucleation at Nonequilibrium State. In this section, let us analyze the characteristics of nucleation at nonequilibrium state such as time-dependent molecular configuration and enthalpy change, which are usually detected from a process of nucleation. The focus will shift from validating solely the equilibrium properties of a system to how nucleation of AP takes place from an aqueous solution. Initial and final configurations of the system are shown in Figure 2,

(3) 2

The calculated diffusion coefficient was 1.48 m /s, which is consistent with the experimental value of 1.47 m2/s.26 Figure 1

Figure 1. Radial distribution functions of N(NH4+)−O(H2O) and Cl(ClO4−)−O(H2O).

represents the RDFs of NH4+ and ClO4− ions in an aqueous solution that provide the information on solution structure of ions and water molecules. The first peak of the N_O is 2.8 Å presenting a first hydration shell around NH4+ ions. It is generated by strong hydrogen bonding between NH4+ and water molecules and corresponds to the value of 2.8 Å obtained by Jorgensen and Gao.27 The first peak of ClO4− ions with water molecules (3.6 Å) was also found to correspond to that obtained by General et al. (3.7 Å) with a high accuracy.28 The simulated values of physicochemical properties are in remarkable agreement with experimental data, suggesting that the system constructed in the present work is sufficient to describe the AP/water solution. In cooling crystallization, the enthalpy change of solution is a critical factor to determine the quantity of heat gain or loss for the separation of solute molecules. Experimentally, ΔHsol of the AP in water is 26.4 kJ/mol.29 It means that the crystallization of AP from an aqueous solution is an exothermic reaction, and therefore, AP can be crystallized by cooling. The enthalpy change of solution can be calculated by the energy change of the solvation enthalpy and crystal energy.6 In the present work, the ΔHsol of AP from an aqueous solution is defined as ΔHsol = ΔHsolv,cation + ΔHsolv,anion − Ecryst

Figure 2. Snapshots of MD simulations on nucleation at supersaturation S = 5.1 (a) and 6.4 (b).

where water molecules are omitted for clarity. The initial configurations represent a diffuse structure of the NH4+ and ClO4− ions from an aqueous solution for each system, where the solute molecules are saturated in water at each temperature. After 10 ns, the systems generate numerous nuclei from a solution, and some of them are merged into a larger one. Naturally, the nuclei at S = 6.4 are shown to be more densely distributed in a given system because the number of solute molecules at S = 6.4 is larger than that at S = 5.1. Figure 3 shows the variation on the potential energy, where the values are normalized by the number of molecules. The magnitude of potential energy at S = 6.4 is larger than that at S

(4)

where the ΔHsolv,cation and ΔHsolv,anion are each solvation enthalpy of cation and anion, and Ecryst is the crystal energy that was used in eq 2. The solvation enthalpy is obtained from the Clausius−Clayperon equation given by ⎛1 ΔG2 ΔG1 1⎞ − = ΔH ⎜ − ⎟ T2 T1 T1 ⎠ ⎝ T2

(5)

where ΔG1 and ΔG2 are free energy changes at temperatures of T1 and T2, respectively, and ΔH is the enthalpy change. The simulation results on solvation free energy were given in Supporting Information 2. The solvation free energy and solvation enthalpy are summarized in Table 2. The negative value of solvation free energy means that the dissolution of NH4+ and ClO4− ions is thermodynamically favored in water. The enthalpy change of solution of AP estimated by eq 4 is 39.9 kJ/mol and underpins the exothermic nucleation process of AP from water in the MD simulation.

Figure 3. Profiles of potential energy at supersaturation S = 5.1 (solid) and 6.4 (open). 5899

dx.doi.org/10.1021/cg501112a | Cryst. Growth Des. 2014, 14, 5897−5903

Crystal Growth & Design

Article

Figure 4. Number of clusters larger than a given size with respect to time at supersaturation S = 5.1 (a) and 6.4 (b).

Figure 5. Free energy change of nucleation calculated by the MD simulation at supersaturation S = 5.1 (a) and 6.4 (b). Open circles, solid circles, and solid lines indicate steady state, equilibrium, and CNT results, respectively.

suggests that the critical size of nucleus decreases at relatively high supersaturation, where the nucleation rate is 3.02 × 1034/ m3 s. The number of nuclei begins to decrease after 2 ns, which is shorter than the required time at S = 5.1. It means that nuclei are rapidly formed at high supersaturation, and subsequently, homogeneous nucleation is terminated in a rather short period. The calculated results of the critical size of nucleus by the CNT (each ncrit for S = 5.1 and 6.4 was calculated to be 12.4 and 8.2, respectively) were found to be in a remarkable agreement with the MD results. Equations of the CNT used were given in Supporting Information 3. The nucleation rates were calculated to be 2.08 × 1032/m3 s and 3.16 × 1033/m3 s for S = 5.1 and 6.4, respectively, and found to be quite comparable with the MD results. The free energy change of nucleation can be evaluated from the cluster size distribution at equilibrium.

= 5.1 because the molecular interaction by electrostatic force with surrounding molecules becomes stronger as the solute concentration increases. The nucleus was found to keep growing as the energy continuously decreases. It implies that crystallization of AP from an aqueous solution is an exothermic reaction. Fundamentally, the reason behind the nucleation under supercooling is from the decrease on kinetic energy that is dependent only on temperature. At equilibrium, the solute molecules maintain a dispersive structure from a solution because there are competitive interactions between kinetic energy and potential energy. Because of the kinetic energy corresponding to a given equilibrium temperature, the dispersive structure does not collapse. However, when a solution was cooled down to a lower temperature, the decreased kinetic energy makes the solute molecules become more frequently trapped into the potential energy field generated by AP molecules. As a consequence, the solute molecules begin to assimilate into nuclei because the solute molecules do not have kinetic energy enough to be disintegrated from nuclei. 3.3. Nucleation Rate. The nucleation rate is defined as the number of nuclei that is larger than the critical size nucleus generated per unit volume per unit time. According to the Y− M method,8−10 the nucleation rate was determined from a slope, which does not change if size of a cluster reached a critical value. In the present work, the NH4ClO4 was assumed to be a monomer. Figure 4 displays simulation results for the supersaturation of 5.1 and 6.4. At S = 5.1, the slope slowly decreases until the cluster size exceeds the value of 12, above which the slopes are almost constant and the nucleation rate is 7.67 × 1033/m3 s. The number of clusters was found to continuously decrease after 3 ns. However, at S = 6.4, the slope does not change when the size of a cluster is larger than 8. It

ΔG = −kT ln(c(n)/c(1))

(6)

where c(n) is the number density of n-mer and c(1) is the number density of monomer. However, the number density of nuclei at equilibrium cannot be directly obtained from the MD simulation. With the help of the Y−M method the equilibrium number density can be predicted by employing the factor, B(i), which is a transition rate of an i-mer to an (i + k)-mer or (i − k)-mer for a given simulation volume (Supporting Information 4), as follows: n ⎡ ⎤ 1 ⎥ ceq(n) = css(n) exp⎢JMD ∑ ⎢⎣ B(i)css(i) ⎥⎦ i=1

(7)

where ceq(n) and css(n) are the equilibrium number density of nuclei and the steady state number density, respectively. JMD is the nucleation rate calculated by the MD simulation. In the present work, the transition of i-mer was sampled every 50 ps. 5900

dx.doi.org/10.1021/cg501112a | Cryst. Growth Des. 2014, 14, 5897−5903

Crystal Growth & Design

Article

Figure 6. Concentration changes with time at supersaturation S = 5.1 (a) and 6.4 (b). Solid and dashed lines denote initial and final constant slope, respectively.

Figure 7. Maximum cluster size at supersaturation S = 5.1 (a) and 6.4 (b).

Figure 8. RDFs of N (NH4+) and Cl (ClO4−) as a function of simulation time at supersaturation S = 5.1 (a) and 6.4 (b).

In the present work, the steady state number density, css(n), was obtained by averaging over the trajectories from the range, 1 ns < time < 3 ns at S = 5.1, and 0 ns < time < 2 ns at S = 6.4, where the constant nucleation rate was observed. Figure 5 displays the trend in free energy change of nucleation calculated by the MD simulation and the CNT. For the steady state number density, the free energy change continuously increases even after the size of a cluster reaches a critical value. The nuclei are in steady state and become larger by merging nuclei after the specific point of time. At the beginning of nucleation, solute molecules are aggregated and form a lot of nuclei in a constant rate, but not merged into a larger one. From this point of view, the concentration of nuclei smaller than ncrit can be assumed to be at quasi-equilibrium. For this reason, the free energy change of nucleation results in unrealistic values for the size of a nucleus larger than ncrit. However, the free energy change of nucleation estimated by the equilibrium number density definitely shows the critical point

where the free energy change starts to decrease. The peak of ΔG/kT at S = 6.4 was found to be smaller than that at S = 5.1, which implies relatively fast nucleation. 3.4. Growth of Nuclei. Figure 6 presents the concentration as a function of time that is calculated by the number of AP molecules less than the critical size of nucleus divided by the total number of molecules in a liquid phase. It is natural that the concentration of AP in solution decreases as crystallization continues in the solution. Hence, the change of concentration could be an indicator to deal with a crystallization process. The slopes at S = 5.1 and 6.4 begin to change after about 3 and 2 ns, respectively, which correspond to the points in which the number of clusters starts to decrease (Figure 4). On the basis of those points, the simulation results could be interpreted in two stages: nucleation and nuclei growth. The former is a step in which a number of nuclei are dramatically formed, whereas the latter is a step in which nuclei tend to grow and some nuclei become merged with each other. 5901

dx.doi.org/10.1021/cg501112a | Cryst. Growth Des. 2014, 14, 5897−5903

Crystal Growth & Design

Article

Figure 9. (a) Configuration of a cluster at supersaturation S = 5.1 and (b) its enlarged image. Blue and cyan beads visualize NH4+ and ClO4− molecules, respectively. Solid red and white lines indicate fully bonded connections of NH4+−NH4+ and ClO4−−ClO4−, whereas dashed lines present the virtual connections that lack at least one molecule.

According to the Y−M method, the nucleation rate was 7.67 × 1033/m3 s and the critical size of nucleus was 12 at S = 5.1, whereas those at S = 6.4 were 3.02 × 1034/m3 s and 8, respectively. Those simulations show that the high supersaturation ratio induces fast nucleation with a small-sized nucleus. The present simulation results were shown to be in good agreement with the CNT model, where the nucleation rates are 2.08 × 1032/m3 s at S = 5.1 and 3.16 × 1033/m3 s at S = 6.4, respectively. The critical size of nucleus (12.4 and 8.2 at S = 5.1 and 6.4, respectively) also considerably agree with those by the MD simulation. The free energy change of nucleation definitely shows the energy barrier that nuclei should be overcome to keep growing. In the present work, a process of nucleation was divided into two stages such as nucleation and nuclei growth. In the nucleation step, spontaneous nucleation is induced by initial large supersaturation and thus nucleation rate is approximately constant. However, in the nuclei growth step, AP molecules tend to be incorporated into existing nuclei rather than continuous homogeneous nucleation. The MD simulations performed herein were found to provide a mechanistic understanding of nucleation of AP from an aqueous solution as well as reliable kinetic properties. Nevertheless, a few cautions on the supersaturation ratio should be kept in mind in the MD simulations of nucleation because strong supersaturation may lead to spinodal decomposition. The energy barrier of nucleation vanishes at much higher supersaturation, and it brings on one cluster of critical size that directly triggers the phase transition.30 In that case, the metastability of the mother phase also depends on the system size. Therefore, a moderate supersaturation ratio and system size should be carefully chosen in the MD simulations of nucleation. However, there still remains the challenge of dealing with the two-step mechanism of nucleation of crystal in solution that the crystalline nuclei appear inside a dense liquid.31,32 The observations of dense liquid for protein, colloid, and organic molecules are reported,33−36 but the dense liquid for smallsized molecules is scarcely observable. In the present work, a cluster of AP molecules was assumed to be a nucleus, and then a comparable nucleation rate was obtained compared with the CNT model. According to the two-step mechanism of nucleation, the free energy of formation of a dense liquid phase and the free energy of formation of a crystalline nucleus inside the dense liquid should be overcome to form a crystalline nucleus in solution. In this respect, it can be carefully surmised that the high accordance of simulation results with the CNT

In the nucleation step, a high supersaturation ratio suitable for overcoming interfacial free energy for the formation of a nucleus induces homogeneous nucleation of AP from an aqueous solution. Once AP nucleation occurs, the supersaturation ratio continuously decreases and leads to a decrease of nucleation rate, after which growth of nuclei is observed. It was confirmed that the nucleus whose size is larger than critical value is not generated after the nucleation step. However, the concentration of AP still continuously decreases. It implies that the AP molecules do not tend to form new nucleus but rather to be incorporated into existing nuclei. However, compared to the maximum cluster size in Figure 7, nuclei of AP are shown to increase gradually during simulation, which means that the generated AP nuclei merged with each other, and the number of nuclei consistently decreases after the nucleation step. Large amplitude of oscillation is originated from the contact between large clusters. At the stage of nucleation, the maximum cluster size does not much vary because small clusters collide with most of the AP molecules. However, the maximum cluster size starts to oscillate dramatically in the stage of growth due to interaction with large clusters. Through analyzing RDFs of AP molecules in Figure 8, the structure of growing clusters was found to be a crystal-like phase because the distinctive RDF peaks of clusters correspond to those of AP crystal (Supporting Information 5). The characteristic peaks of RDF become clear as time elapsed. Figure 9 shows the largest nuclei of AP, which exhibits an ordered pattern of the NH4+ and ClO4− molecules.

4. CONCLUSIONS We simulated the homogeneous nucleation of AP from an aqueous solution to gain a fundamental understanding of phase transition from a liquid to a solid phase and the kinetics of nucleation. The simulation system constructed in the present work was found to reproduce the physicochemical properties such as lattice energy of AP crystal, solution density, diffusion coefficient, and RDFs. Furthermore, the thermodynamics of AP crystallization from an aqueous solution was verified by providing the enthalpy change of AP from a liquid to a solid, ΔHsol, whose positive value is consistent with the experimental fact that dissolution of AP in water is an endothermic reaction. The MD simulations show that the total energy of the system decreases as nuclei grow, which implies that the nucleation of AP from an aqueous solution is an enthalpically favored process. 5902

dx.doi.org/10.1021/cg501112a | Cryst. Growth Des. 2014, 14, 5897−5903

Crystal Growth & Design

Article

(23) Parrinello, M.; Rahman, A. J. Appl. Phys. 1981, 52, 7182−7190. (24) William, F. L. Solubilities; American Chemical Society: Washington, D.C., 1965. (25) Keenan, A. G.; Siegmund, R. F. Q. Rev. Chem. Soc. 1969, 23, 430−452. (26) Hiquily, N.; Clifton, M. J. J. Chem. Eng. Data 1984, 29, 371− 373. (27) Jorgensen, W. L.; Gao, J. J. Phys. Chem. 1986, 90, 2174−2182. (28) General, I. J.; Asciutto, E. K.; Madura, J. D. J. Phys. Chem. B 2008, 112, 15417−15425. (29) Mullin, J. W. Crystallization, 3rd ed.; Butterworth-Heinemann: Woburn, MA, 1993. (30) Wedekind, J.; Chkonia, G.; Wölk, J.; Strey, R.; Reguera, D. J. Chem. Phys. 2009, 131, 114506. (31) Erdemir, D.; Lee, A. Y.; Myerson, A. S. Acc. Chem. Res. 2009, 42, 621−629. (32) Vekilov, P. G. Cryst. Growth Des. 2010, 10, 5007−5019. (33) Brode, M. L.; Berland, C. R.; Pande, J.; Ogun, O. O.; Benedek, G. B. Proc. Natl. Acad. Sci. U.S.A. 1991, 88, 5660−5664. (34) Bonnett, P. E.; Carpenter, K. J.; Dawson, S.; Davey, R. J. Chem. Commun. 2003, 698−699. (35) Tan, P.; Xu, N.; Xu, L. Nat. Phys. 2014, 10, 73−79. (36) Vekilov, P. G. Nat. Mater. 2012, 11, 838−840.

model may be ascribed to the low energy barrier of transition of a dense liquid to a crystalline nucleus. The simple structure and high polarization of partial charge can make it easy to form a crystalline nucleus. In the future we further simulate the nucleation of AP on the crystal faces, where the growth behavior of a crystal is treated at the atomic level.

■

ASSOCIATED CONTENT

S Supporting Information *

Details of computational methods, solvation free energy, CNT, transition rate, and RDF of AP crystal. This material is available free of charge via the Internet at http://pubs.acs.org.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Tel: +82-2-705-8680. Fax: +82-2705-8630. Notes

The authors declare no competing financial interest.

■

■

ACKNOWLEDGMENTS This work was supported by Next-Generation Converged Energy Materials Research Center and by Hanwha and Agency for Defense Development (UC120019GD; Republic of Korea).

■

NOTE ADDED AFTER ASAP PUBLICATION This paper was published ASAP on October 21, 2014, with an error to equation 2. The corrected version was reposted on October 24, 2014.

REFERENCES

(1) Jiménez-Á ngeles, F.; Firoozabadi, A. J. Phys. Chem. C 2014, 118, 11310−11318. (2) Salvalaglio, M.; Vetter, T.; Giberti, F.; Mazzotti, M.; Parrinello, M. J. Am. Chem. Soc. 2012, 134, 17221−17233. (3) Anwar, J.; Boateng, P. K. J. Am. Chem. Soc. 1998, 120, 9600− 9604. (4) Anwar, J.; Boateng, P. K.; Tamaki, R.; Odedra, S. Angew. Chem., Int. Ed. 2009, 121, 1624−1628. (5) Gavezzotti, A. Chem.Eur. J. 1999, 5, 567−576. (6) Cheong, D. W.; Boon, Y. D. Cryst. Growth Des. 2010, 10, 5146− 5158. (7) Valeriani, C.; Sanz, E.; Frenkel, D. J. Chem. Phys. 2005, 122, 194501. (8) Yasuoka, K.; Matsumoto, M. J. Chem. Phys. 1998, 109, 8451− 8462. (9) Yasuoka, K.; Matsumoto, M. J. Chem. Phys. 1998, 109, 8463− 8470. (10) Suh, D.; Yasuoka, K. J. Phys. Chem. B 2011, 10631−10645. (11) Banerjee, S.; Briesen, H. J. Chem. Phys. 2009, 131, 184705. (12) Mackerell, A. D.; Feig, M.; Brooks, C. L. J. Comput. Chem. 2004, 25, 1400−1415. (13) Bjelkmar, P.; Larsson, P.; Cuendet, M. A.; Hess, B.; Lindahl, E. J. Chem. Theory Comput. 2010, 6, 459−466. (14) Dewar, M. J. S.; Zoebisch, E. G.; Healy, E. F.; Stewart, J. J. P. J. Am. Chem. Soc. 1985, 107, 3902−3909. (15) Wagner, E. L. J. Chem. Phys. 1962, 37, 751−759. (16) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269−6271. (17) Berendsen, H. J. C.; van der Spoel, D.; van Drunen, R. Comput. Phys. Commun. 1995, 91, 43−56. (18) Lindahl, E.; Hess, B.; van der Spoel, D. J. Mol. Model. 2001, 7, 306−317. (19) van der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. J. Comput. Chem. 2005, 26, 1701−1718. (20) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. J. Chem. Theory Comput. 2008, 4, 435−447. (21) Darden, T.; York, D.; Pedersen, L. J. Chem. Phys. 1993, 98, 10089−10092. (22) Bussi, G.; Donadio, D.; Parrinello, M. J. Chem. Phys. 2007, 126, 014101. 5903

dx.doi.org/10.1021/cg501112a | Cryst. Growth Des. 2014, 14, 5897−5903