Monte Carlo Simulations of Binary Mixtures of Nitrotoluene Isomers

13752

J. Phys. Chem. B 2009, 113, 13752–13760

Monte Carlo Simulations of Binary Mixtures of Nitrotoluene Isomers with n-Decane† Katie A. Maerzke and J. Ilja Siepmann* Departments of Chemistry and of Chemical Engineering and Materials Science, UniVersity of Minnesota, 207 Pleasant Street SE, Minneapolis, Minnesota 55455 ReceiVed: March 24, 2009; ReVised Manuscript ReceiVed: June 2, 2009

Previous experimental studies [Sliwinska-Bartkowiak et al., Chem. Phys. Lett. 1983, 94, 609] indicate that o-nitrotoluene/n-decane and m-nitrotoluene/n-decane mixtures are of particular interest because they exhibit anomalous behavior in the nonlinear dielectric effect near the upper critical solution point, and it has been surmised that these effects are due to an increase in cluster formation as each mixture approaches its critical point. In this work, Monte Carlo simulations are performed using the TraPPE force field at multiple temperatures, concentrations, and electric field strengths to examine the structural microheterogeneity of these mixtures and determine the dielectric constants and nonlinear dielectric effect. At low field strength, the simulations indicate substantial structural microheterogeneities, but these persist over a wide range of conditions. At high field strength (an order of magnitude higher than in the experimental work), the simulations indicate a field-induced segregation. 1. Introduction Noncrystalline mixtures forming a homogeneous bulk phase, such as supercritical, liquid, and solid solutions, often exhibit significant microheterogeneity, that is, composition fluctuations on the microscopic (molecular) scale. Microheterogeneity can be inferred through deviations from ideal solution behavior (e.g., the entropy of mixing may be smaller than that for an ideal solution), but the structural details of these mixtures are challenging to probe by experimental means. One experimental observable that provides indirect information on the microheterogeneity of liquids is the nonlinear dielectric effect (NDE). The NDE is the change in electric permittivity (∆ε) induced by a strong electric field (E), given by

εE - ε0 ∆ε ) 2 E E2

(1)

where εE is the permittivity of the medium with the applied field E, and ε0 is the permittivity without the applied field. Significant insights into the NDE in liquids are due primarily to Piekara and co-workers in Poland, who started to investigate this topic in the 1930s.1 In fact, the nonlinear dielectric effect is sometimes referred to as the Piekara effect. Debye theory predicts that in the presence of a strong electric field, the electric permittivity will decrease (i.e., ∆ε < 0).2 Weakly interacting molecules, such as diethyl ether and chloroform, do, in fact, exhibit the behavior predicted by Debye theory; that is, the NDE is always negative. However, molecules with relatively strong intra- or intermolecular interactions can exhibit positive values of ∆ε under certain conditions. Pure liquids of some polar molecules, such as nitrobenzene and nitrotoluene, exhibit positive NDE values; however, upon dilution by a nonpolar solvent, the NDE decreases, becoming negative at sufficiently dilute concentrations.1–3 This implies that the positive NDE value is due to the interactions between the †

Part of the “H. Ted Davis Special Section.” * Corresponding author. E-mail: [email protected].

polar molecules at high concentration. The mechanism proposed by Piekara4 as explained by Davies2 and Kielich5 involves the transitory pairwise association of the polar molecules. In the absence of a strong electric field, the molecules will align in an antiparallel arrangement, which reduces the effective dipole moment and, hence, the electric permittivity. However, in the presence of a strong electric field, the molecules align in a parallel configuration, which increases the effective dipole moment and, hence, also increases the electric permittivity.2 At higher concentrations, the probability for forming dimers and larger-order aggregates increases, causing the NDE to increase. At low concentrations, the solution consists mostly of isolated polar molecules; hence, the electric permittivity decreases when an electric field is applied, in agreement with Debye’s theory.5 Of particular interest have been solutions of o- and mnitrotoluene in nonpolar hydrocarbons. Measurements performed by Sliwinska-Bartkowiak and co-workers indicate that mixtures of o-nitrotoluene with linear alkanes exhibit a positive NDE near the upper critical solution point.6–8 A slight increase in NDE is observed at temperatures as far as 30 K above the critical point, with the NDE becoming larger as the temperature approaches Tc. Positive NDE values near the critical point can be explained in terms of molecular clusters. Far away from the critical point, the mixture is relatively homogeneous and exhibits negative NDE values. As the mixture approaches the critical point, clusters of nitrotoluene molecules begin to form. As stated above, in the presence of a strong electric field, the molecules in the clusters will align in such a way as to increase the dipole moment, thus increasing the electric permittivity and giving rise to a positive NDE value. The closer the mixture gets to the critical point, the larger the critical fluctuations become, which creates an even larger dipole moment and NDE value. However, the reason these systems are particularly interesting is that mixtures of m-nitrotoluene with the same n-alkanes, although exhibiting similar NDE behavior, do not phase separate.6,7,9,10 This raises the possibility of a hypothetical (and experimentally inaccessible) critical point which exists below the melting point, but nevertheless affects the properties of the stable liquid phase. For instance, as the m-nitrotoluene/n-decane mixture approaches

10.1021/jp902631f CCC: $40.75  2009 American Chemical Society Published on Web 07/23/2009

MC Simulations of Nitrotoluene Isomers, n-Decane

J. Phys. Chem. B, Vol. 113, No. 42, 2009 13753

its hypothetical critical point, critical fluctuations begin to appear, which causes an increase in the NDE. However, before the mixture is able to phase separate, it reaches the freezing point and solidifies. The reasons for the difference in behavior between the two isomers are unclear, and the difference may be limited to the solid state. Moving the nitro group from the ortho to the meta position leads to an increases in the dipole moment by a small amount (from 3.69 to 4.17 D).6 This small change is evidently not sufficient to significantly effect some properties, such as the liquid density (e.g., the liquid density of o-nitrotoluene at 293 K is 1.1611 g/cm3, but the liquid density of m-nitrotoluene is 1.1581 g/cm3),11 whereas the change in intramolecular structure (and, maybe, dipole moment) significantly influences the free energy balance between the solid and liquid forms (e.g., the melting point of o-nitrotoluene is 262.6 K, but the melting point of m-nitrotoluene is 288.5 K).11 Sliwinska-Bartkowiak, et al. have determined that the critical point of the mixture o-nitrotoluene/n-decane occurs at Tc ) 276 K and xc ) 0.53 mol fraction o-nitrotoluene6,8 and have estimated that the hypothetical critical point of the mixture m-nitrotoluene/n-decane occurs at Tc ) 273 K and xc ) 0.54 mol fraction m-nitrotoluene (i.e., very close to that of the other isomer),12,13 but this mixture solidifies at 281 K.10,13 In a first attempt to understand these phase transitions, Ratajczak et al. have performed some Monte Carlo simulations for the system m-nitrotoluene/n-decane.13 However, due to their choice of a very simplistic single-site Lennard-Jones (LJ) model, which does not reflect the dipolar nature of nitrotoluene or the articulated nature of n-decane, these simulations offer only qualitative insight into the microscopic nature of these systems. In the present work, the TraPPE-UA force field is used to examine the structural microheterogeneity and determine the dielectric constants of the binary mixtures o-nitrotoluene/ndecane and m-nitrotoluene/n-decane at multiple temperatures, concentrations, and electric field strengths. By applying an electric field, we are able to calculate the nonlinear dielectric effect. To our knowledge, no one else has used a simulation to determine the NDE of a mixture. The remainder of this paper is organized as follows: the details of the electric field and dielectric constant calculations are described in Section 2. The details of the force field can be found in Section 3, and the simulation details can be found in Section 4. The results of the dielectric constant calculations and nonlinear dielectric effect can be found in Section 5.1. The microhetereogeneity of these systems is examined in Section 5.2. The effects of the electric field on the orientation of the nitrotoluene molecules and the liquid density is examined in Section 5.3. The conclusions of this work can be found in Section 6.

which eliminates the need to calculate an angle and, hence, reduces the computational complexity. The calculation of the dielectric constant depends on the method of treating first-order electrostatics and the boundary conditions. For simulations that use an Ewald sum with tinfoil boundary conditions, the dielectric constant in the absence of an applied electric field is given by14–16

ε0 ) 1 +

4π [〈M2〉0 - 〈M〉02] 3VkBT

(4)

where M is the total system dipole moment, the subscript 0 indicates a zero-field ensemble average. In principle, 〈M〉 should be zero for an isotropic system, but due to system size effects and the finite length of the simulation, it generally has a small but nonzero value. In the presence of a small applied electric field, E, the dielectric constant εE is given by16

εE ) ε0 +

( )[

4π 1 90V kBT

3

3〈M4〉0 - 5〈M2〉02]E2

(5)

where V is the volume of the system with the applied electric field, ε0 is the zero-field dielectric constant, and the subscript 0 again indicates a zero-field ensemble average. 3. Force Field For this work, we have used the TraPPE-UA force field for n-decane,17 the flexible nitro group,18 and the rigid nine-site benzene model.19 The TraPPE-UA force field employs pseudoatoms for all CHx groups (where 0 e x e 4), which are located at the position of the carbon atom. The nonbonded interactions are described by the simple pairwise-additive Lennard-Jones (LJ) and Coulomb potentials:

[( ) ( ) ]

unonbonded(rij) ) 4εij

σij rij

12

-

σij rij

6

+

qiqj 4πε0rij

(6)

where rij, εij, σij, qi, and qj are the separation, LJ potential well depth, LJ diameter, and partial charges, respectively. The unlike LJ interactions are determined using the Lorentz-Berthelot combining rules:20

1 σij ) (σii + σjj) and εij ) (εiiεjj)1/2 2

(7)

2. Methodology To probe the effects of an applied electric field on the structure of these liquid mixtures, the simulations must be performed with and without an applied electric field. The energy from an applied electric field is given by

ufield(θ) ) -µ · E ) -|µ||E|cos θ

(2) ubend(θ) )

where µ is the molecular dipole moment and E is the applied electric field. To simplify the calculation, we chose to apply the field in the z-direction only,

ufield(ri,z) ) µzEz ) -qiri,zEz

Nonbonded interactions are calculated for intermolecular interactions and intramolecular interactions of atoms separated by four or more bonds. Rigid bond lengths are used for bonded interactions. Bond angles are treated using a harmonic potential,

(3)

kθ (θ - θeq)2 2

(8)

where θ is a bond angle, θeq is the equilibrium value for that bond angle, and kθ is the force constant. For beads separated by three bonds, torsional interactions are computed using a cosine series, with

13754

J. Phys. Chem. B, Vol. 113, No. 42, 2009

Maerzke and Siepmann TABLE 1: Dielectric Constants Computed with the TraPPE-UA Force Field at Zero Applied Fielda o-nitrotoluene xNT 1.0 0.8 0.6 0.5 0.4 0.2 a

298 K 9.25 6.03 4.054 3.276 2.612 1.672

273 K 9.65 6.24 4.23 3.327 2.744 1.681

m-nitrotoluene 298 K 7.74 4.71 3.475 2.978 2.354 1.582

273 K 6.52 5.02 3.62 2.91 2.373 1.622

Subscripts indicate uncertainties in the final digit.

Figure 1. Dielectric constant computed during the production period for o-nitrotoluene (solid lines) and m-nitrotoluene (dashed lines) at xNT ) 1.0 (top) and 0.5 (bottom) and T ) 298 K (left) or 273 K (right). The different colors indicate the results for independent simulations. The initial value of ε0 computed from the fluctuation formula (eq 4) is equal to unity irrespective of the configuration. Note the different scales on the y-axis. 6

utor ) a0 +

∑ ai cos(iφ)

(9)

i)1

for the O-N-C-C torsion and

utor ) a1[1 + cos φ] + a2[1 - cos 2φ] + a3[1 + cos 3φ] (10)

4. Simulation Details

Figure 2. Nonlinear dielectric effect for o-nitrotoluene (top) and m-nitrotoluene (bottom). The symbols indicate the simulation data computed using eq 5 at 298 K and an electric field strength of E ) 0.009 V/Å. The numerical values for pure o- and m-nitrotoluene are 7100 ( 15800 and 2800 ( 14300 Å2/V2, respectively. The black, red, and green lines show the experimental data6 at 298, 288, and 278 K and an electric field strength of E ) 0.004-0.009 V/Å.

Monte Carlo (MC) simulations in the isobaric-isothermal (NpT) ensemble21 were carried out for the binary mixtures o-nitrotoluene/n-decane and m-nitrotoluene/n-decane to examine the structural microheterogeneity, dipole ordering, and densification effects. The program MCCCS (Monte Carlo for Complex Chemical Systems) developed by the Siepmann group was used for these simulations. This code has been extensively verified for simulations of rigid and flexible molecules in various ensembles through extensive comparison to analytical results for simple model systems and previous simulation data obtained with independent simulation codes. In particular, the addition of the external electric field was tested against the analytical orientational distribution of an isolated molecule with fixed dipole moment. The simulations were performed for 256 molecules at two temperatures (T ) 273 and 298 K), six concentrations (1.0, 0.8, 0.6, 0.5, 0.4, and 0.2 mol fraction nitrotoluene), five applied electric field strengths (E ) 0.0, 0.009, 0.1, 0.5, and 1.0 V/Å), and a pressure of 1 bar to examine the structural microheterogeneity and effects of the electric field. Simulations were not performed for pure n-decane because the TraPPE-UA decane model does not contain partial charges; hence, the applied electric field would have no effect. The applied electric field strength of 0.009 V/Å corresponds to the largest experimental field strength.6 Translational and rotational moves along with conformational changes were performed to equilibrate the system. The confor-

mational changes for n-decane and the flexible nitro group were performed using coupled-decoupled configurational-bias Monte Carlo (CBMC).22–25 Volume moves were performed to equilibrate the system with an external pressure bath. The fraction of volume moves and the corresponding maximum displacement were adjusted so that on average a volume move was accepted once every 10 Monte Carlo (MC) cycles (1 MC cycle ) N moves, where N is the number of molecules in the system). Of the remaining moves, approximately 1/3 were CBMC regrowths, with the rest divided evenly between translations and rotations. In addition, MC simulations in the canonical (NVT) ensemble were performed for these same systems at zero field strength to calculate the dielectric constant. For all systems, a spherical cutoff of rcut ) 14 Å with analytic tail corrections26 was used for the LJ interactions. The Ewald summation technique with tin foil boundary conditions was used to compute the Coulombic interactions.14,27,28 The real-space cutoff was equal to rcut, and the Ewald sum convergence parameter was set to 3.5/rcut. For every state point, four independent simulations were performed. Systems were equilibrated for at least 300 000 MC cycles before starting production runs. In the NpT and NVT ensemble, production runs consisted of 450 000 and of approximately 1 800 000 MC cycles, respectively. Calculated uncertainties are given as the standard error of the mean.

for the alkane torsions.



Figure 3. Snapshots of the mixtures of o-nitrotoluene (top) and m-nitrotoluene (bottom) with decane at zero field; T ) 298 K; and xNT ) 0.8 (left), 0.5 (middle), and 0.2 (right). The decane chains, aromatic rings, nitrogen, and oxygen atoms are shown in gray, green, blue, and red, respectively.

Figure 4. Nitrotoluene-nitrotoluene local mole fraction enhancement (L) for o-nitrotoluene (left) and m-nitrotoluene (right) at zero electric field (T ) 298 K, solid lines; T ) 273 K, dashed lines). Note the varying scale on the y-axis.

Figure 5. Decane-decane local mole fraction enhancement (L) for mixtures with o-nitrotoluene (left) and m-nitrotoluene (right) at zero electric field (T ) 298 K, solid lines; T ) 273 K, dashed lines). Note the varying scale on the y-axis.

Here, it should be noted that the lower temperature investigated in this study (273 K) lies about 5% below the experimental melting point for m-nitrotoluene (288.5 K).11 However, our simulations do not yield any indication of crystallization nor any dramatic difference between the two temperatures. It is wellknown that simulations for a finite system size in a cubic simulation box (i.e., one that is usually not commensurate with the crystal unit cell) allow for significant undercooling of a liquid phase. Hence, we can not tell whether the liquid phase is stable or only metastable for the TraPPE-UA force field at 273 K.

5. Results and Discussion 5.1. Dielectric Properties. Since the dielectric constant is determined via a fluctuation formula (see eqs 4 and 5), long simulations are required to obtain a converged value. To assess the convergence of our simulations, we have plotted the cumulative average of the dielectric constant as a function of run length for select systems (see Figure 1). As can be seen, the dielectric constants for the independent simulations are significantly better converged at 298 K than at 273 K. Thus, the statistical uncertainties at the lower temperature remain quite large.

13756


Figure 6. Nitrotoluene-nitrotoluene center-of-mass RDF for onitrotoluene (left) and m-nitrotoluene (right) at zero electric (T ) 298 K, solid lines; T ) 273 K, dashed lines). Note the varying scale on the y-axis.

The numerical values of the dielectric constant are listed in Table 1. At 298 K, the TraPPE-UA model underpredicts the dielectric constants of the isomers by a factor of 3 compared to the experimental values11 of 26.26 and 24.95 for o- and m-nitrotoluene, respectively (older experimental values29 are 27.4 and 23.0 ( 0.8, respectively). Lowering the temperature by 25 K does not have a substantial effect on the computed dielectric constant, but the uncertainties are larger at the lower temperature. The dielectric constant decreases smoothly as the concentration of nitrotoluene is decreased. Experimental values are not available for the mixtures. The reason for the underestimation of the dielectric constant is our choice of model. The TraPPE-UA force field does not reproduce the dipole moment of either nitrotoluene model very well because it does not allow for a polarization of the π cloud of the aromatic ring30,31 (and the nitro group carries the same charges as for nitroalkanes).18 The average liquid-phase dipole moments at 298 K and zero field of the TraPPE-UA model for o- and m-nitrotoluene are 3.45 and 3.46 D (note that the TraPPEUA model is nonpolarizable but flexible and that this small difference is due to slightly different average conformations), compared to an experimental gas-phase results of 3.69 and 4.17 D,6 respectively. Furthermore, we computed the dipole moments for both isomers using the CM4 charge model of Kelly et al.32 with 1-octanol as the continuum solvent, as done in the development of the TraPPE-EH model for aromatic heterocycles and substituted rings.30,31 To determine the dipole moments, the structure of each isomer was first optimized using Kohn-Sham density functional theory at the B3PW91/6-31++G(d,p) level of theory in the Gaussian 03 program.33 The optimized gasphase geometry was then used to optimize the charge distribution (and, hence, determine the dipole moment) at the B3P86/ 6-31+G(d,p) level of theory using the Minnesota Gaussian solvation model (MN-GSM) implemented in Gaussian 03.34 The resulting gas-phase dipole moments are 4.237 D for o-

Maerzke and Siepmann nitrotoluene and 5.109 D for m-nitrotoluene, which are significantly larger than the experimental values. Furthermore, the CM4 calculation predicts solution dipole moments of 5.881 and 6.935 D for o- and m-nitrotoluene, respectively. Hence, not only does our model severely underestimate the dipole moment of both nitrotoluene isomers, but the difference between the dipole moments of the two isomers is much smaller than is observed experimentally or calculated with the CM4 charge model. Nevertheless, despite the very similar dipole moments of the TraPPE-UA models for the two isomers, it correctly predicts that the dielectric constant for m-nitrotoluene is somewhat lower than that for o-nitrotoluene. Figure 2 shows a graphical comparison of calculated NDE at 298 K with the experimental data.6 First of all, it should be noted that the statistical uncertainties for the simulation data are rather large, particularly for high nitrotoluene concentrations. The total system dipole moment appearing in eqs 4 and 5 is a collective property that may not be sampled very well in a Monte Carlo trajectory utilizing rotational moves of a single molecule, and, hence, the statistical uncertainty increases with increasing nitrotoluene concentration. Furthermore, the term appearing in square brackets in eq 5 requires significant sampling of the tail of the total system dipole moment distribution. Since simulations in the canonical ensemble were not carried out for conditions with an external electric field, we cannot judge whether such calculations would offer more precise estimates of the NDE. However, it appears that the simulations yield a positive value and the correct order of magnitude of the NDE for the two nitrotoluene isomers. For both mixtures, there is also a hint of a secondary maximum at nitrotoluene concentrations of 0.4-0.5, followed by a minimum (possibly even negative values of the NDE) before rising to larger positive values for xNT g 0.8. Thus, we are cautiously stating that the calculated NDE shows a qualitatively similar concentration dependence as observed experimentally,6 which gives credence to the following structural analysis. As one would expect from eq 5, changes in the electric field strength do not lead to substantial changes in the calculated NDE because the effect on the molar volume is small for the applied fields studied here (see below). 5.2. Microheterogeneity. Snapshots of two mixtures at select state points (see Figure 3) give a visual impression of their microheterogeneous structures, that is, nonuniform spatial composition distribution. At high concentration of nitrotoluene, xNT ) 0.8, the polar molecules form a continuous aggregate interspersed by decane molecules. At xNT ) 0.2, individual nitrotoluene aggregates are clearly visible. The pattern of spatial distributions appears to be similar for both nitrotoluene isomers. To provide a quantitative measure of the microheterogeneity of these mixtures, the local mole fraction enhancement (LMFE) was determined as a function of distance from a chosen molecule.35 The local mole fraction is the average number of molecules of a given type up to a distance r from the chosen molecule divided by the total number of molecules of any type, where the number of molecules of each type is determined from the center of mass number integral (NI). To obtain the local mole fraction enhancement, the local mole fraction is divided by the bulk mole fraction; that is, for a binary system with molecule types a and b.

Laa(r) )

[

NIaa(r) 1 xa NIaa(r) + NIab(r)

]

(11)



Figure 7. Snapshots of the mixtures of o-nitrotoluene (top) and m-nitrotoluene (bottom) with n-decane at T ) 298 K, xNT ) 0.5, and E ) 0.0 (left), 0.009 (middle), and 0.1 V/Å (right). Colors as in Figure 3.

Figure 8. Distribution of the dipole vector orientation for nitrotoluene molecules with respect to the z-direction (direction of applied field) for o-nitrotoluene (left) and m-nitrotoluene (right) at 298 K and varying electric field strengths: E ) 0.0 (black), 0.009 (red), 0.1 (green), 0.5 (blue), and 1.0 V/Å (magenta).

Figures 4 and 5 show a comparison of the nitrotoluene-nitrotoluene and decane-decane LMFE for both isomers at all concentrations and at both temperatures. First of all, the LMFE shows aggregation at all conditions; that is, the solvation shell around a nitrotoluene molecule is enriched by other nitrotoluene molecules, and the same holds true for decane around decane. Second, there is relatively little difference between the mixtures containing o- or m-nitrotoluene, as also inferred from the experimental NDE measurements. Third, the nitrotoluene-nitrotoluene LMFEs are largest at the lowest nitrotoluene concentra-

Figure 9. Nitrotoluene-nitrotoluene local mole fraction enhancement (L) for o-nitrotoluene (left) and m-nitrotoluene (right) at 298 K and varying electric field strengths: E ) 0.0 (black), 0.009 (red), 0.1 (green), 0.5 (blue), and 1.0 V/Å (magenta). Note the varying scale on the y-axis.

tion (as one should expect from any preferential like-like aggregation); that is, there is no indication of a special enhancement near concentrations of xNT ≈ 0.5, where the secondary maximum in the NDE is observed. Fourth, the temperature dependence is relatively small, with the lower temperature consistently yielding slightly larger LMFEs (with the exception of the highest concentration for the m-nitrotoluene mixture). Upon examining the nitrotoluene-nitrotoluene radial distribution function, we again see that temperature has relatively little effect on the structure of the systems (see Figure 6),

13758


Maerzke and Siepmann

Figure 12. Temperature dependence of the nitrogen-nitrogen RDF (top) and of the distribution of dipole vector orientations for nitrotoluene molecules within the first N-N solvation shell (bottom) for onitrotoluene (left) and m-nitrotoluene (right) at xNT ) 0.5 and E ) 0.009 V/Å: 298 K (solid lines) and 273 K (dashed lines). Figure 10. Nitrogen-nitrogen RDF for o-nitrotoluene (left) and m-nitrotoluene (right) at 298 K and varying electric field strengths. E ) 0.0 (black), 0.009 (red), 0.1 (green), 0.5 (blue), and 1.0 V/Å (magenta).

TABLE 2: Effects of the Electric Field on Liquid Densities (in kg/m3) at 298 K and Different Mol Fractions of the Nitrotoluene Isomersa E

although the systems do become slightly more structured at low concentrations of nitrotoluene, which is where we see the largest LMFE. For all state points, one can observe two distinct shells with peak positions at r ≈ 4.2 and 8 Å. The increase in the peak heights with decreasing nitrotoluene concentration can again be attributed to favorable dipolar interactions that have the largest effect on the LMFE at the lowest concentration. The fact that the height of the first minimum exceeds unity for xNT g 0.5 is another reflection of the LMFE. Irrespective of concentration (and temperature), the height of the first peak is always higher for m-nitrotoluene than for o-nitrotoluene, which we attribute to slightly more favorable packing.

0.8

0.6

0.5

0.4

0.2

8534 8591 8601

7883 7911 7921

8511 8531 8541

7873 7881 7891

0.0 0.1 1.0

11735 11672 11722

10405 10431 10431

o-Nitrotoluene 9431 9001 9431 9001 9461 9041

0.0 0.1 1.0

11502 11451 11434

10274 10251 10241

m-Nitrotoluene 9312 8924 9321 8901 9331 8921

a

Figure 11. Distribution of dipole vector orientations for nitrotoluene molecules within the first N-N solvation shell for o-nitrotoluene (left) and m-nitrotoluene (right) at 298 K and varying electric field strengths. E ) 0.0 (black), 0.009 (red), 0.1 (green), 0.5 (blue), and 1.0 V/Å (magenta).

1.0

Subscripts indicate uncertainties in the final digit.

Overall, these structural observations agree with the microheterogeneity and their similarity for the two isomers inferred from the NDE measurements.6,10,13 However, our simulations at zero field do not support the hypothesis of SliwinskaBartkowiak et al.6 that the mixtures form larger and larger clusters as the upper critical solution point is approached, but seem to indicate that the distribution of cluster sizes shifts to larger clusters with increasing xNT. However, we do not know the upper critical solution point (nor the melting line) of our model systems. Determining liquid-liquid and solid-liquid phase diagrams of relatively complex molecules via molecular simulation is difficult, and the locations of the critical point and the melting line are very sensitive to the detail of the underlying force field.36,37 Given the limitations of our model discussed above in Section 5.1, it is unlikely that we would be able to accurately reproduce the liquid-liquid coexistence curve for these systems. 5.3. Effects of the Electric Field. Although the effects of temperature and composition observed in this work are somewhat gradual, it is found that the structure changes dramatically once the applied field exceeds a certain threshold. This is most evident for the equimolar mixture. Snapshots for these mixtures at field strengths of 0.0, 0.009 (the largest value used in the experimental studies),6 and 0.1 V/Å (approximately one order of magnitude larger) are shown in Figure 7. Visually, the extent of microheterogeneity appears to be quite similar at the two lower field strengths, but there is a clear indication of segregation into essentially two aggregates (either spherical or lamellar) at the higher field strength.

MC Simulations of Nitrotoluene Isomers, n-Decane Let us first examine the orientation of the dipole vectors with respect to the electric field. For simplicity, we have chosen to use the C-N bond vector as an approximation to the dipole vector. Given the charge distribution in our model, the C-N vector is, indeed, extremely close to the actual dipole vector. As can be seen in Figure 8, the dipole vectors show a uniform distribution for all mixtures at zero field. At E ) 0.009 V/Å, there is a small (but beyond the statistical uncertainties) preference for alignment with the field (e.g., for xNT ) 0.5, orientations with a parallel alignment are about 1.5 times more likely than orientations with an antiparallel alignment). Upon an increase in the field strength to 0.1 V/Å, there is a dramatic change in the distribution, and more than 90% of the molecules possess an alignment angle smaller than 90°, and this percentage climbs to more than 99% for the two strongest fields. We have also examined the dependence of the LMFE on electric field strength (see Figure 9). Again, there is little difference between zero field and E ) 0.009 V/Å. For the mixtures with o-nitrotoluene, a significant increase in LMFE for r > 5 Å with increasing field strength is observed for xoNT e 0.6, whereas a similar increase is found only for 0.4 g xmNT e 0.5 for m-nitrotoluene. For both isomers, the increases are particularly large for the equimolar mixture, that is, near the upper critical solution composition. With the exception of xoNT ) 0.2, these structural changes emerge only for field strengths higher than those used in the experimental studies. Nevertheless, we would like to argue that the increase in LMFE is related to the anomalous NDE behavior of these mixtures; that is, that the application of an external electric field is driving the system toward (at least) local phase separation. In recent theoretical work,38 significant changes in vapor-liquid and liquid-liquid equilibria upon application of an electric field were predicted for the Stockmayer model. Molecular dynamics simulation of nonpolarizable water models demonstrate electrofreezing.39 Furthermore, it should be noted that the experimental NDE measurements apply microsecond electric pulses, a duration that may be too short to allow for field-induced phase separation of the nitrotoluene-decane mixtures. Beyond the increase in LMFE, we also observe that the shape of the LMFE changes with increasing field strength. In particular, LMFE shows a maximum at r ≈ 5.3 Å and a steep decrease for shorter separations for both nitrotoluene isomers. Since the nitrotoluene molecules show significant alignment at the higher field strengths (see Figure 8), dipole-dipole repulsion disfavors the very close approach of two molecules. Considering the effect of the applied electric field at a constant temperature on the liquid structure (see Figure 10), we see that as the field strength increases, the system becomes more and more structured. However, the increase in the peak heights from low to high field does not exceed the increase in the LMFE, with the exception of xNT ) 0.8, where increased peak heights are observed for the RDF without a comparable change in the LMFE. As the proposed mechanism for the NDE in solutions of nitrotoluene relies on the relative orientation of the dipole moments (see Section 1), the relative orientations of the dipole vectors of neighboring nitrotoluene molecules is examined here. For nitrotoluene molecules whose nitrogen atoms are within the first N-N solvation shell (as determined from the N-N RDF; see Figure 10), we have calculated the angle between the dipole vectors. As the position of the first solvation shell does not vary greatly with either temperature or electric field strength, we have chosen 5.6 Å for all systems containing o-nitrotoluene and 6.2 Å for all systems containing m-nitrotoluene as the outer limit

J. Phys. Chem. B, Vol. 113, No. 42, 2009 13759 of the first solvation shell. The proposed NDE mechanism predicts that in the absence of an applied electric field, the molecules will align in an antiparallel arrangement (180° angle), whereas in the presence of a strong electric field, the molecules will align in a parallel arrangement. In the absence of an electric field, we observe a modest preference for antiparallel arrangements (see Figure 11). This preference changes only very slightly upon application of the weakest field investigated here. For the intermediate field (E ) 0.1 V/Å), the preference has reversed, and parallel arrangements are now modestly favored. For the two highest fields, we observe a strong preference for parallel alignment of neighboring molecules (see Figure 11). This should come as no surprise because most nitrotoluene molecules are aligned with field under these conditions (see Figure 8). In Figure 12, we depict the temperature dependence of the nitrogen-nitrogen RDF and of the dipole vector orientations at low field for the equimolar mixtures (i.e., near the experimental critical composition and for the experimental field strength). Decreasing the temperature by 25 K has a minor effect on the RDF that is mostly explained by the increase in the LMFE (very similar to that observed without the field, see Figure 6) but no discernible effect on the orientational distribution of neighboring nitrotoluene molecules. The application of a strong electric field can cause solid materials or liquid crystals to deform (a phenomenon known as electrostriction), affecting the density. Although the nitrotoluene molecules were found to orientationally align with the field at the higher field strengths, our simulations allowed for only isotropic volume changes (i.e., a safer choice because it is not known a priori whether the liquid will be isotropic or anisotropic). Kusalik has observed a slight increase in the density for simpler models at higher applied fields than used here.16 The specific densities obtained from our mixture simulations are listed in Table 2. Averaging over all mixtures, the density increase is 0.2% upon increasing from zero applied field to E ) 1.0 V/AA, which falls within the statistical uncertainties. Furthermore, it should be noted that the TraPPE-UA force field reproduces the experimental density11 at 298 K for both isomers to within 1%, but overestimates the density difference. 6. Conclusions Although the TraPPE-UA model significantly underestimates the dielectric constants for o- and m-nitrotoluene, the simulations yield a positive NDE value for pure liquids of both isomers. Albeit the statistical uncertainties are very large, it appears that the simulations show a composition dependence of the NDE for mixtures with n-decane that is in qualitative agreement with the experimental data.6 Structural analysis does not yield significant differences between the o-nitrotoluene/n-decane and m-nitrotoluene/n-decane mixtures, as also inferred from the experimental NDE measurements.6 Mixtures of both isomers with n-decane exhibit structural microheterogeneity at all concentrations and both temperatures investigated here. Most interestingly, the simulations indicate that the application of an electric field increases the local mole fraction enhancements of like species (i.e., like-like aggregation) and induces, at least, local segregation of the species. This effect is most pronounced for the equimolar mixtures and may offer an explanation of the anomalous NDE for these mixtures. In the future, it may be worthwhile to confirm these qualitative conclusions using a nitrotoluene model that yields a more accurate dielectric constant and to explore whether macroscopic phase separation under an applied electric field would be

13760


observed using much larger system sizes or simulations in an open ensemble, such as the Gibbs ensemble. In addition, one may want to carry out simulations for related systems that do not show anomalous NDE behavior, such as the mixture of nitrobenzene and benzene.4 Acknowledgment. Financial support from the National Science Foundation (CBET-0553911 and CBET-756641) is gratefully acknowledged. Part of the computer resources were provided by the Minnesota Supercomputing Institute. References and Notes (1) Piekara, A. Phys. ReV. 1932, 42, 448–449. (2) Hill, N.; Vaughn, W. E.; Price, A. H.; Davies, M. Dielectric Properties and Molecular BehaViour; Van Nostrand Reinhold: London, 1969. (3) Piekara, A. Proceedings of the 31st International Meeting of the Societe de Chimie Physique: Nonlinear BehaViour of Molecules, Atoms, and Ions in Electric, Magnetic, or Electromagnetic Fields; Elsevier: Amsterdam, 1979. (4) Piekara, A.; Chelkowski, A. J. Chem. Phys. 1956, 25, 794–795. (5) Kielich, S. Dielectric and Related Molecular Processes; The Chemical Society: London, 1972; Vol. 1. (6) Sliwinska-Bartkowiak, M.; Szurkowski, B.; Hilczer, T. Chem. Phys. Lett. 1983, 94, 609–613. (7) Sliwinska-Bartkowiak, M. Ber. Bunsenges. Phys. Chem. 1990, 94, 64–70. (8) Sliwinska-Bartkowiak, M.; Szurkowski, B.; Hilczer, T. Phase Transitions 1989, 18, 77–86. (9) Sliwinska-Bartkowiak, M.; Szurkowski, B.; Hilczer, T. Phys. Lett. 1981, 81A, 411–412. (10) Sliwinska-Bartkowiak, M. Phys. Lett. A 1988, 128, 84–88. (11) Lide, D. R., Baysinger, G., Berger, L. I., Goldberg, R. N., Kehiaian, H. V., Kuchitsu, K., Rosenblatt, G., Roth, D. L., Zwillinger, D., Eds.; CRC Handbook of Chemistry and Physics, 88th ed.; CRC Press: Boca Raton, FL, 2007. (12) Sliwinska-Bartkowiak, M. J. Phys.: Condens. Matter 1993, 5, 407– 422. (13) Ratajczak, B.; Sliwinska-Bartkowiak, M.; Coasne, B.; Gubbins, K. E. J. Non-Cryst. Solids 2007, 353, 4565–4569. (14) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: New York, 1987. (15) de Leeuw, S. W.; Perram, J. W.; Smith, E. R. Annu. ReV. Phys. Chem. 1986, 37, 245–270. (16) Kusalik, P. G. Mol. Phys. 1994, 81, 199–216. (17) Martin, M. G.; Siepmann, J. I. J. Phy. Chem. 1998, 102, 2569– 2577. (18) Wick, C. D.; Stubbs, J. M.; Rai, N.; Siepmann, J. I. J. Phys. Chem. B 2005, 109, 18974–18982.

Maerzke and Siepmann (19) Wick, C. D.; Siepmann, J. I.; Klotz, W. L.; Schure, M. R. J. Chromatogr., A 2002, 954, 181–190. (20) Maitland, G. C.; Rigby, M.; Smith, E. B.; Wakeham, A. Intermolecular Forces: Their Origin and Determination; Clarendon Press: Oxford, 1987. (21) McDonald, I. R. Mol. Phys. 1972, 23, 41–58. (22) Siepmann, J. I.; Frenkel, D. Mol. Phys. 1992, 75, 59–70. (23) de Pablo, J. J.; Laso, M.; Siepmann, J. I.; Suter, U. W. Mol. Phys. 1993, 80, 55–63. (24) Vlugt, T. J. H.; Martin, M. G.; Smit, B.; Siepmann, J. I.; Krishna, R. Mol. Phys. 1998, 94, 727–733. (25) Martin, M. G.; Siepmann, J. I. J. Phys. Chem. B 1999, 103, 4508– 4517. (26) Wood, W. W.; Parker, F. R. J. Chem. Phys. 1957, 27, 720–733. (27) Ewald, P. Ann. Phys. 1921, 64, 253–287. (28) Madelung, E. Phys. Z. 1918, 19, 524–532. (29) Lide, D. R., Ed. CRC Handbook of Chemistry and Physics, 72nd ed.; CRC Press: Boca Raton, FL, 1991. (30) Rai, N.; Siepmann, J. I. J. Phys. Chem. B 2007, 111, 10790–10799. (31) Rai, N.; Bhatt, D.; Siepmann, J. I.; Fried, L. E. J. Chem. Phys. 2008, 129, 194510. (32) Kelly, C. P.; Cramer, C. J.; Truhlar, D. G. J. Theor. Comput. Chem. 2005, 1, 1133–1152. (33) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, Jr., J. A.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A. Gaussian 03, ReVision C.01; Gaussian, Inc.: Wallingford, CT, 2004. (34) Chamberlin, A. C.; Kelly, C. P.; Thompson, J. D.; Xidos, J. D.; Li, J.; Hawkins, P. D.; Winget, G. D.; Zhu, T.; Rinaldi, D.; Liotard, D. A.; Cramer, C. J.; Truhlar, D. G.; Frisch, M. J. MN-GSM, Version 6.0; University of Minnesota: Minneapolis, MN, 2006, 55455-0431. (35) Chen, B.; Potoff, J. J.; Siepmann, J. I. J. Phys. Chem. B 2001, 105, 3093–3104. (36) Zhang, L.; Siepmann, J. I. J. Phys. Chem. B 2005, 109, 2911– 2919. (37) Zhao, X. S.; Chen, B.; Karaborni, S.; Siepmann, J. I. J. Phys. Chem. B 2005, 109, 5368–5374. (38) Ga´bor, A.; Szalai, I. Mol. Phys. 2008, 106, 801–812. (39) Svishchev, I. M.; Kusalik, P. G. J. Am. Chem. Soc. 1996, 118, 649–654.

JP902631F

Monte Carlo Simulations of Binary Mixtures of Nitrotoluene Isomers

Recommend Documents