7980
J. Phys. Chem. C 2007, 111, 7980-7985
Melting Point Determination from Solid-Liquid Coexistence Initiated by Surface Melting Ali Siavosh-Haghighi and Donald L. Thompson* Department of Chemistry, UniVersity of Missouri, Columbia, Missouri 65211 ReceiVed: January 11, 2007; In Final Form: March 22, 2007
A coexisting solid-liquid (s-l) system of nitromethane is created by surface-induced melting. A nitromethane crystal with a free surface is simulated by molecular dynamics (MD) in the constant-volume and -energy (NVE) ensemble for initial conditions generated by short MD simulations of the constant-volume and -temperature (NVT) ensemble at temperatures slightly above the melting point. Melting starts at the surface, initiating a solid-liquid interface, and the temperature drops as the system moves toward a state of equilibrium in which the solid and liquid phases coexist. The temperature at which the coexisting solid and liquid reach equilibrium is taken to be the melting point. The melting points of crystals with exposed (100), (010), and (001) crystallographic faces are predicted to be 238, 245, and 242 K, respectively. The predicted melting points are in good agreement with experiment (244.7 K) and previous simulations. The approach to equilibrium during the NVE simulation is monitored by calculating the orientational order parameter, diffusion coefficient, and density, which provide insights into the melting mechanism. The Sorescu-Rice-Thompson [J. Phys. Chem. B 2000, 104, 8406] force field, which accurately describes the inter- and intramolecular motions, was used.
1. Introduction The melting point of a perfect crystal is higher than the thermodynamic melting point as result of the free energy barrier to the formation of a solid-liquid (s-l) interface. A straightforward way to model thermodynamic melting is to initiate a constant volume and energy (NVE) simulation with the solid in contact with the liquid. The solid and the liquid portions of the system are separately prepared by performing simulations near an estimate of the melting temperature; then the two are joined in a single simulation with periodic boundary conditions.1,2 If the energy is too low there will be a net l f s change and if it is too high then s f l will dominate. The temperature and pressure change as the coexisting solid-liquid system comes to dynamic equilibrium. While practical, this method requires some effort to determine the conditions for equilibrium of the coexisting solid-liquid system. In the present study we have applied a method originally developed by Ercoslessi et al.3,4 This approach is more straightforward than the standard coexistence method1,2 because it avoids the separate preparations and then matching of solid and liquid simulation supercells. The basic idea is to initiate molecular dynamics (MD) simulations for the constant-volume-energy (NVE) ensemble for a solid in contact with vacuum at a temperature slightly above the melting point. Melting will begin at the free surface, forming a solid-liquid interface which progresses inward. The temperature decreases as a function of time as the melting progresses. If the initial temperature is relatively close to the melting temperature, the temperature of the system will asymptotically approach the melting point Tmp.The temperature T, computed over the time * Corresponding author.
of the simulation when both solid and liquid are present, can be fit to
T ) T0 + ae-bt
(1)
where t is the simulation time, T0 is the asymptotic value of the temperature, and a is a constant equal to 1 K. If a true s-l interface is created T0 is an estimate of Tmp. We have used this method to determine the melting point of crystalline nitromethane with exposed surfaces (100), (010), and (001), and to investigate the molecular-level mechanism of melting. The NVE simulations provide more realistic results for the time evolution of a melting system than do simulations in which constraints must be used to maintain the ensemble conditions, because applying a thermostat or barostat in a simulation of the NPT ensemble affects the dynamics of the system. We have monitored the molecular orientations and mobility as melting occurs in the NVE simulations to determine the details of the melting mechanism for a molecular solid. 2. Methods The initial configuration of the crystal structure was based on the neutron diffraction scattering data reported by Trevino and Rymes.5 The simulations were performed using the DL_POLY (v. 2.15) code.6 Periodic boundary conditions were applied in all directions, thus two exposed surfaces are necessary as illustrated in Figure 1. The integration step size was 0.75 fs, the cutoff range for the potential 11 Å, and the verlet shell 0.5 Å. Ewald sum method was used to calculate electrostatic interactions in periodic system. Equilibration in the NVE runs was determined by monitoring the density, orientational order parameter, diffusion coefficient, and temperature as functions of time; these quantities were computed every 500th integration
10.1021/jp070242m CCC: $37.00 © 2007 American Chemical Society Published on Web 05/15/2007
Melting Point Determination from s-l Coexistence
J. Phys. Chem. C, Vol. 111, No. 22, 2007 7981
Figure 1. Top panel: initial configuration of the simulation cell with the (100) surface exposed to vacuum. Bottom panel: illustration of the rectangular fixed-volume slabs (defined by the dashed vertical lines) approximately one-molecular size in width (2.56 Å) for which quantities were computed to monitor the progression of the melting through the crystal. Because of the symmetry of the simulation cell, the properties of the molecules in slabs n and n′ were combined to compute averages of quantities as functions of the distance from the exposed surface or the center of the simulated crystal.
TABLE 1: Simulations NVT simulation run 1 2 3 4 5 6 7 8 9 10 11
simulated crystal supercell (Å)a
exposed surface
10 × 4 × 3 157.32 × 25.28 × 26.19 10 × 4 × 3 157.32 × 25.28 × 26.19 10 × 4 × 3 157.32 × 25.28 × 26.19 10 × 4 × 3 157.32 × 25.28 × 26.19 10 × 4 × 3 157.32× 25.28 × 26.19 10 × 4 × 3 157.32 × 25.28 × 26.19 5×8×3 26.22 × 151.68 × 26.19 5×4×6 26.22 × 25.28 × 157.14 5×4×3 78.66 × 25.25 × 26.19 5×4×3 26.22 × 75.84 × 26.19 5×4×3 26.22 × 25.25 × 78.75
NVE simulation
time (ps)
T (K)
Eb
time (ps)
T0 (K)c
100
15
255
-23.248
1700
233 ( 6
100
15
260
-23.235
1700
238 ( 6
100
15
265
-23.229
1700
237 ( 8
100
15
270
-23.221
1700
238 ( 8
100
15
276
-23.214
1700
240 ( 8
100
15
280
-23.202
1700
242 ( 10
010
15
270
-23.206
1700
245 ( 8
001
15
270
-23.210
1700
242 ( 10
100
3.75
260
-23.174
900
237 ( 15
010
3.75
260
-23.195
900
245 ( 10
001
3.75
260
-23.203
900
241 ( 11
a Crystal dimensions (unit cells) and simulation box size. b The energy zero is for isolated molecules at 0 K. Units: kJ mol-1 molecule-1. c The uncertainties are the root-mean-square fluctuations.
step. The density, orientational order parameter, and diffusion coefficient were calculated within rectangular fixed-volume slices approximately one-molecular size in width (2.56 Å) as illustrated in Figure 1. Because of the symmetry of the simulation cell, the properties of the molecules in slabs n and n′ were combined to compute averages of quantities as functions
of the distance from the exposed surface or the center of the simulated crystal. The simulations that were performed are summarized in Table 1. In one series of simulations a large supercell (10 × 4 × 3 unit cells) was used to ensure sufficient time for the formation of a true solid-liquid interface. These simulations
7982 J. Phys. Chem. C, Vol. 111, No. 22, 2007
Siavosh-Haghighi and Thompson
Figure 2. Plot of the temperature as a function of time calculated for NVE simulations of supercells of dimensions 10 × 4 × 3 unit cells for initial conditions prepared by running short NVT simulations at 260 and 280 K. The temperature approaches the melting point as the system approaches equilibrium.
were done for the (100) crystallographic face exposed to vacuum. Another series of simulations were done using smaller supercells. Melting initiated at the (010) face was simulated using a supercell of dimensions 5 × 8 × 3 and for the (001) exposed surface a 5 × 4 × 6 supercell was used. The predicted melting points are statistically the same for the various sizes of simulation supercells. The larger supercell allowed for a more detailed analysis of the mechanism. The simulation supercell was constructed by placing the crystal in an empty larger orthorhombic box; see column 2 of Table 1. Conditions corresponding to 1 atm and temperatures 255, 260, 265, 270, 276, and 280 K were assigned to the atoms, and short NVT simulations were run to obtain initial conditions for the NVE simulations; see column 4 and 5 of Table 1. The NVE simulations were run until equilibration was achieved, generally in less than 1700 ps; see column 6 of Table 1. The total energies for the NVE runs are given in column 7 of Table 1. The temperature was calculated as a function of the simulation time and fit to eq 1. The values obtained for T0 are given in the last column of Table 1 with root-mean-square fluctuations in the temperature in the NVE simulations; these results will be discussed in the next section. The density, orientational order parameter, and diffusion coefficient were computed to determine the state of the system during the simulations. The orientational order parameter is defined as the ensemble average:
〈|cos θ|〉 )
1 N
∑i |ei‚k|
(2)
where ei is the unit vector along C-N bond and k is the unit vector along (i.e., along the [001] direction) of the unit cell. The angle between the C-N bond and edge c is close to 0° and 180° in the crystal; thus, 〈|cos θ|〉 ) 0.93 in the crystal and decreases as the crystalline order diminishes. The diffusion coefficient D was calculated using
〈∆r2〉 ) 6Dt
(3)
where r is the position vector of a molecule relative to the center of mass of the system and t is the time. The density was calculated in each of the slices as defined in Figure 1, except for slices 10 and 10', by considering a molecule to be within the slice in which its center of mass was
Figure 3. The (a) orientational order parameter and (b) density as functions of time computed for all the molecules in NVE simulations of supercells of dimensions 10 × 4 × 3 unit cells for initial conditions corresponding to 260 K (blue curve), and 280 K (red curve). (c) Extent of melting in simulations with initial temperatures 260 K (blue curve) and 280 K (red curve) illustrated by the percent solid as a function of time.
located. The average density of whole slab was obtained by averaging the densities of all slices. The force field is that developed by Sorescu et al.7 It has been shown to accurately predict the solid- and liquid-state properties of nitromethane.8,9 It was used with minor modifications by Agrawal et al. 10 to study the melting of nitromethane, giving results in good agreement with experiment. The intermolecular forces are described by the Buckingham (exp -6) potential plus Coulombic interactions with fixed partial charges. The intramolecular interactions are represented by a sum of Morse potentials for the bond stretches, harmonic oscillators for the bond angles, and cosine sum for the torsion angles. A complete description of the force field7 and values of the constants10 are available elsewhere. 3. Results and Discussion The density, order parameter, extent of melting, diffusion coefficient, and temperature as functions of time were computed in NVE simulations of crystalline nitromethane for energies corresponding to initial temperatures over the range 255-
Melting Point Determination from s-l Coexistence
J. Phys. Chem. C, Vol. 111, No. 22, 2007 7983
Figure 4. Fits of NVE simulation results to eq 1 for initial conditions corresponding to 255, 260, 265, 270, 276, and 280 K (see column 5 of Table 1 for corresponding E values) for the 10 × 4 × 3 supercell. The values of T0, which we take to be the melting point, for the range of temperature between 260 and 280 K are in the range 237.1-241.8 K; the uncertainties are in the range 5.7-10.1 K for NVE; thus, the T0 values are statistically the same. The initial temperature 255 K is too close to actual melting point for a true solid-liquid interface to be established; thus, it results in a quasi-liquid layer that does not have the properties of liquid nitromethane. The quasi-liquid surface is influenced by the immediate solid phase.
280 K; see columns 6 and 7 of Table 1. These quantities were used to determine the state of the system. The temperature at which the equilibrium between the solid and liquid was achieved is taken to be the melting point. Analyses of the behavior of the order parameter and diffusion coefficient were used to study the mechanism of melting. The temperature as a function of time during NVE simulations of supercells of dimension 10 × 4 × 3 with initial conditions prepared by running short NVT simulations at 260 and 280 K are shown in Figure 2. The density, orientational order parameter, and percent solid as functions of time for NVE simulations initially equilibrated at 260 and 280 K are shown in frames a, b, and c, respectively, of Figure 3. The order parameter, given by eq 2, is 0.93 for the crystal, and decreases as the melting progresses (see Figure 3a). It levels off at an equilibrium value, which is determined by the extent of melting. The magnitude of the order parameter for liquid nitromethane, calculated in a separate simulation, is 0.78. The state of the system can also be determined by monitoring the density. The average density as a function of time is shown in Figure 3(b) for initial temperatures 260 K (blue curve) and 280 K (red curve). To ascertain that a true liquid phase is formed we calculated the density of the melted portion; at 239 K the density of the liquid portion is 1.13 ( 0.12 g/cm3, which is essentially the value (1.17 g/cm3) calculated by Sorescu et al.9 These results show the same behavior as the order parameter. Because there is incomplete melting at equilibrium, the final equilibrium values shown in Figure 3b are more than that of liquid nitromethane, which for this force field at 240 K is 0.0118 Å-3.9 The extent of melting is shown in Figure 3c. The percent solid is determined by characterizing the solid layers (as defined in Figure 1) as either solid or liquid based on the order parameter. The fraction of the system that remains solid at the end of the simulations range from about 37% to 79%. The temperature as a function of time in the NVE simulations for initial temperatures over the range 255-280 K (see column 5 of Table 1) were fit to eq 1; the resulting curves are shown in Figure 4 and the asymptotes T0 are given in the last column of Table 1. If the initial temperature is too low the result is the formation of an amorphous surface layer, not a solid-liquid
Figure 5. (a) The temperature as a function of time for the NVE simulations of supercells of dimensions 10 × 4 × 3 unit cells for initial conditions corresponding to 270 K for exposed face (100), 5 × 8 × 3 unit cells for exposed surface (010), and 5 × 4 × 6 unit cells for exposed face (001). (b) Fits of the results in panel a to eq 1. The temperature decreases faster for the crystal with the (100) free surface. The rates of decrease in temperature for the (010) and (001) exposed surfaces are similar and result in higher asymptotic values, 244.7 and 242.5 K, respectively, compared to that (237.7 K) for the exposed (100) face.
interface. For instance, the diffusion coefficient for the outermost layer of the slab initially equilibrated at 255 K is equal to 0.03 Å2/ps, much lower than that of liquid (∼0.1 Å2/ps).14 The NVE simulation for initial temperature 255 K does not result in a solid-liquid equilibrium, rather in a solid with a disorder surface. The NVE simulations for the other initial temperatures reach equilibrium temperatures in the range 237-242 K, all within the range of uncertainties thus statistically the same. The predicted melting point is thus in the range 237-242 K. We did not observe any evaporation in the simulations. The pressure remained close to zero; always less than the fluctuation. For example, the pressure of the equilibrated system for the highest initial temperature studied, 280 K, was calculated to be 0.22 ( 0.37 kbar. The values of the melting point predicted by this method are 14-19 K lower than the value, 256 K, obtained by the standard solid-liquid coexistence method used by Agrawal et al.10 The pressure was 1 atm in the coexisting solid-liquid simulations performed by Agrawal et al., while in current study there is no external pressure applied to the simulation cell. We believe that the solid and liquid phases are not very sensitive to variations of 0 to 1 atm of pressure; thus, we would assume the differences in the computed melting points obtained by the different approaches should be comparable except for the differences resulting from the types of ensembles. Thus, we ascribe the difference to the differences in the thermodynamic averaging. Agrawal et al. determined the melting point by direct simulation of the solid-liquid coexistence system in the NVT ensemble at various temperatures. The calculation of a thermodynamic quantity requires long simulations (ergodicity theorem).11 Furthermore, the slow rate of energy transfer from the liquid to
7984 J. Phys. Chem. C, Vol. 111, No. 22, 2007
Siavosh-Haghighi and Thompson
Figure 6. The order parameter, density, and diffusion coefficient as functions of time for the core (left column) and the surface (right column) for NVE simulations of 5 × 4 × 3 crystals of nitromethane with exposed (100) surfaces for initial conditions corresponding to 260 K. The core is approximately one unit cell in width, slices 1 and 1′ in Figure 1, and surface layer is defined by 5 and 5′ in Figure 1; thus, the average of the surface layer is effectively over one-half the width of a unit cell. The dashed lines indicate the values of the quantities calculated in separate liquid-state simulations. Orientational disorder of the molecules precedes translational freedom.
the solid phase in NVE simulations requires long simulation times (on the order of nanoseconds).12,13 The length of the solid-liquid coexistence simulation performed by Agrawal et al. was 225 ps., while the trajectories in the present study were integrated for 1700 ps, and the temperature versus time curves were extrapolated to infinity to determine the temperatures at equilibrium. Also, we note that in another study14 of nitromethane using the same force field we calculated the melting point to be 251 K using 1200 ps simulations of the NVT ensemble, which is in close agreement with current results. In an earlier study,14. we calculated 239 ( 5 K for the melting point by fitting the thickness of the melted layer of a slab of nitromethane computed in an NVT simulations at different temperatures to the equation l(T) ) K(Tmp - T)-1/3, which was suggested by Mori et al.15 and Pluis et al.16 for van der Waals solids. The present results are consistent with that result. Simulations were performed for exposed (100), (010), and (001) crystallographic surfaces using supercells of dimensions 10 × 4 × 3, 5 × 8 × 3, and 5 × 4 × 6 unit cells, respectively, for initial conditions corresponding to 270 K. Figure 5a shows the changes in the temperature as functions of time for the three cases and the curves obtained by fitting the results to eq 1 are shown in Figure 5b. The rate of decrease in the temperature is faster for the melting initiated at the (100) free surface. Taking the values of T0 in eq 1 to be the melting points predicts 238, 245, and 242 K, respectively, for the (100), (010), and (001) cases. The present results are in accord with results of an earlier
study,14 using a different approach, in which we calculated values of 251, 266, and 266 K, respectively, for the melting points of (100), (010), and (001) surface-initiated melting. Figure 6 shows the variation of the order parameter, density, and diffusion coefficient of the core (left column) and the surface (right column) calculated in a NVE simulation for the total energy corresponding to 260 K of a 5 × 4 × 3 crystal with the (100) face exposed. There was complete melting in the case of the exposed (100) surface, however, the predicted melting point is still a good estimate. The dashed lines in each plot indicate the value of the quantity for the liquid state, which we calculated in separate single-phase simulations of liquid nitromethane. Orientational disordering of the molecules precedes mobility. The translational freedom of molecules was determined by calculating the diffusion coefficient using eq 3 within the volume elements (illustrated in Figure 1); the results are shown in Figure 6. The diffusion coefficient was averaged over 18.75 ps. While orientational disorder of molecules at the core reaches the liquidstate value within about 550 ps, the values of the density and diffusion coefficient at the core do not reach the liquid-state values until about 620 ps. The time gap between onsets of molecular orientational disorder and translational freedom at the surface is even larger. The orientational disorder reaches that of the liquid state at about 70 ps, while translational freedom is not reached until about 160 ps. Since the surface molecules are coordinated to fewer neighboring molecules than those in the core, the energy barriers to reorientation are lower.
Melting Point Determination from s-l Coexistence
J. Phys. Chem. C, Vol. 111, No. 22, 2007 7985
4. Conclusions
References and Notes
We have tested a simple and straightforward method suggested by Ercoslessi et al.3,4 for calculating melting points for molecular solids. This approach is more straightforward than the standard coexistence method1,2 because it avoids the separate preparations and then matching of solid and liquid simulation supercells. A coexisting nitromethane solid-liquid (s-l) equilibrium is created by surface-induced melting. A crystal with a free surface is simulated by MD in the constant-volume-energy (NVE) ensemble for initial conditions generated by short MD simulations of the constant-volume-temperature (NVT) ensemble at temperatures slightly above the melting point. In the NVE simulation melting begins at the surface and the temperature decreases as the system moves toward the solid-liquid coexistence equilibrium. The temperature at which the coexisting solid and liquid are at equilibrium is taken to be the melting point. The computed melting points are 238, 245, and 243 K, respectively, for the crystals with exposed (100), (010), and (001) crystallographic faces. The predicted melting points are in good agreement with experiment17 (244.7 K) and previous simulations.10,14,18 By monitoring the orientational order parameter, diffusion coefficient, and density as the system approaches equilibrium during the NVE simulation insight into the melting mechanism is obtained. The molecules first gain rotational freedom followed by translational freedom.
(1) Morris, J. R.; Wang, C. Z.; Ho, K. M.; Chan, C. T. J. Phys. ReV. B 1994, 49, 3109. (2) Morris, J. R.; Song, X. J. Chem. Phys. 2002, 116, 9352. (3) Ercolessi, F.; Tomagnini, O.; Iarlori, S.; Tosatti, E. In Nanosources and Manipulation of Atoms Under High Fields and Temperature: Applications; Binh, V. T., Garcia, N., Dransfeld, D., Eds.; Nation-ASI Series E; Kluwer: Dordrecht, The Netherlands, 1993; Vol. 235, p 185. (4) For a review of the method and some applications, see: Di Tolla, F. D.; Tosatti, E.; Ercolessi, F. In Monte Carlo and Molecular Dynamics of Condensed Matter Systems, Binder, K., Ciccotti, G., Eds.; SIF: Bologna, 1996; Conference Proceedings Vol. 49, Chapter 14, p 346. (5) Trevino, S. F.; Rymes, W. H. J. Chem. Phys. 1980, 73, 3001. (6) Smith, W. F.; Leslie, M.; Forester, T. R. DL POLY 2.15; : Daresbury Laboratory at Daresbury: Warrington, U.K., 2003. (7) Sorescu, D. C.; Rice, B. M.; Thompson, D. L. J. Phys. Chem. B 1997, 101, 798. (8) Sorescu, D. C.; Rice, B. M.; Thompson, D. L. J. Phys. Chem. B 2000, 104, 8406. (9) Sorescu, D. C.; Rice, B. M.; Thompson, D. L. J. Phys. Chem. 2001, 105, 9336. (10) Agrawal, P. M.; Rice, B. M.; Thompson, D. L. J. Chem. Phys. 2003, 119, 9617. (11) Reichl, L. E. A Modern Course in Statistical Physics, 2nd ed.; John Wiley & Sons: New York, 1997; p 296. (12) Nada, H.; van der Eerden, J. P.; Furukawa, Y. J. Cryst. Growth 2004, 266, 297. (13) Ferna´ndez, R. G.; Abascal, J. L. F.; Vega, C. J. Chem. Phys. 2004, 124, 144506. (14) Siavosh-Haghighi, A.; Thompson, D. L. J. Chem. Phys. 2006, 125, 184711. (15) Mori, H.; Okamoto, H.; Isa, S. Physica 1974, 73, 237. (16) Pluis, B.; Denier, van der Gon, A. W.; Frenken, J. W. M.; van der Veen, J. F. Phys. ReV. Lett. 1987, 59, 2678. (17) Jones, W. M. and Giauque, W. F. J. Am. Chem. Soc. 1947, 69, 983. (18) Zheng, L.; Luo, S.-N.; and Thompson, D. L. J. Chem. Phys. 2006, 124, 154504.
Acknowledgment. We thank Dr. Thomas D. Sewell for reading and commenting on the manuscript. This work was supported by a Multidisciplinary University Research Initiative (MURI) grant managed by the U.S. Army Research Office.
