J. Phys. Chem. C 2007, 111, 365-373
365
Simulating Polymorphic Phase Behavior Using Reaction Ensemble Monte Carlo John K. Brennan,*,† Betsy M. Rice,† and Martin Lı´sal‡,§ U. S. Army Research Laboratory, Weapons and Materials Research Directorate, Aberdeen ProVing Ground, Maryland 21005-5066, E. Ha´ la Laboratory of Thermodynamics, Institute of Chemical Process Fundamentals, Academy of Sciences of the Czech Republic, 165 02 Prague 6-Suchdol, Czech Republic, and Department of Physics, Faculty of Science, J.E. Purkinje UniVersity, 400 96 U Ä stı´ n. Lab., Czech Republic ReceiVed: July 20, 2006; In Final Form: October 10, 2006
In this work, we present an application of the reaction ensemble Monte Carlo (RxMC) method (Johnson, J. K.; Panagiotopoulos, A. Z.; Gubbins, K. E. Mol. Phys. 1994, 81, 717; Smith, W. R.; Trˇ´ıska, B. J. Chem. Phys. 1994, 100, 3019) for predicting solid-state structural transitions, namely, conformational polymorphism. The RxMC technique provides a means of overcoming high-energy transition barriers which may be caused, for example, by sterically hindered atoms or molecules. We demonstrate the method on nitromethane (NM), CH3NO2, whose behavior has been studied extensively. There are many morphologies of crystalline NM created by the rotation of the methyl group about the C-N bond. In applying the RxMC method to the phase transition behavior of crystalline NM, we treat the rotation of the methyl group as a conformational isomerization reaction. We then predict the equilibrium concentrations of the rotamers as a function of pressure and temperature. We find that the simpler, less computationally expensive rigid model of nitromethane used in the RxMC simulation can capture the same behavior as the fully flexible model. Investigations of such systems using the RxMC method allow for the identification of the key components of the mechanisms through the decoupling of the molecular degrees of freedom. In the scope of other similar solid-state transitions, the RxMC method can be a powerful tool in the interpretation of experimental measurements as well as in the confirmation of experimental findings. In this work, we also introduce an algorithm within the RxMC framework that ensures the number of accepted reactions is equivalent for each reaction type.
I. Introduction The macroscopic behavior of a solid is strongly dependent on its crystal structure. Any transformations in its crystal structure caused by, for example, temperature or pressure, will dramatically affect a wide range of its properties, including mechanical, thermal, electrical, optical, and magnetic. Pressureinduced changes in crystalline morphology are particularly interesting and are known to influence a wide range of phenomena, including deflagration and detonation processes,1 phase transitions,2 and various geophysical events.3 A fundamental understanding of pressure-induced changes in crystalline morphology could, in principle, be used in the design or modification of materials. Various crystallography methods, including Raman spectroscopy, X-ray diffraction, and neutron diffraction,4 along with density and calorimetry measurements can provide key data for solid-state transitions. Unfortunately, these measurements only provide information about the timeaveraged, space-averaged contents of a unit cell in a crystal; they cannot disclose the nature of the fluctuations from this averaged structure or the molecular-level rearrangements that occur during these transitions.5 It is well-established that molecular simulation is a valuable tool in predicting and understanding behavior that is not amenable to experimentation. Unfortunately, pressure-induced solid-state structural transitions pose considerable hurdles to not only experimental measurements but also molecular simulation approaches. Foremost, the development of an accurate model that describes these systems †
U. S. Army Research Laboratory. Academy of Sciences of the Czech Republic. § J.E. Purkinje University. ‡
over a large pressure range requires significant computational resources and human effort. Further, these molecules can exist in different molecular conformations, giving rise to polymorphism. In this work, we apply the reaction ensemble Monte Carlo (RxMC) method6,7 to the simulation of materials that exhibit pressure-induced solid-state structural transitions. The method does not explicitly trace a solid-solid coexistence boundary but rather determines the equilibrium concentrations of the conformers of the polymorph at any given temperature or pressure. A key feature of the RxMC method is that it is not required to follow naturally occurring trajectories through highenergy states found in many solid-state transitions. In this work, we also introduce an algorithm within the RxMC framework to efficiently select reactions so that the number of accepted reactions is equivalent for each reaction over a simulation run; details are given in the Appendix. We demonstrate the method on nitromethane (NM), CH3NO2, shown in Figure 1. This system was chosen because (1) its properties and behavior have been studied extensively by both experiment and other computational approaches (see, e.g., refs 8 and 9 and references therein), (2) it exhibits a pressureinduced structural transition, and (3) a molecular simulation interaction potential is available and has been shown in various molecular dynamics (MD) simulations10 to accurately predict numerous properties of condensed phase NM, including the pressure-induced structural transition. While our main objective in this work is to present a computational method that will allow for the prediction of solid-state structural transitions, we will also identify key features of the potential energy surface used in the simulation that affect the system’s molecular structure at
10.1021/jp0646170 CCC: $37.00 © 2007 American Chemical Society Published on Web 12/07/2006
366 J. Phys. Chem. C, Vol. 111, No. 1, 2007
Brennan et al. III. Methodology
Figure 1. Molecular structure of nitromethane corresponding to the crystalline state at 4.2 K and 1 atm. The atom labels are consistent with the atom indices discussed in the text.
various conditions, and which influence the pressure-induced structural transition. This study will demonstrate that over a wide pressure range the molecular structure of solid nitromethane can be modeled using a substantially reduced interaction potential compared to that used in previous theoretical studies.10 II. BackgroundsNitromethane The low-temperature, ambient-pressure crystalline structure of NM belongs to the P212121 space group and has four molecules in the unit cell.11 The internal rotation of the methyl group in the crystal has received significant attention, since experimental evidence indicates that this motion is almost completely governed by intermolecular interactions.12-14 In the gas phase, the barrier to rotation is ∼6 cal/mol,15 indicating the methyl group is essentially a free rotor under these conditions. In the crystalline state, the temperature dependence of the degree of libration of this motion at atmospheric pressure has been measured, as has the activation barrier to this rotation.12 Of interest in this study is the influence of pressure on the internal rotational motion of the methyl group. X-ray diffraction measurements at room temperature for pressures ranging from 0.3 to 6 GPa showed that, relative to the low-temperature, atmospheric-pressure structure, the space group and packing arrangement remain unchanged.13 However, the measurements indicated a change in the orientation of the methyl group with increasing pressure. For pressures below 0.6 GPa, the methyl group is believed to be freely rotating in the crystal;11,13 however, due to experimental difficulties, this has yet to be confirmed in the laboratory. At intermediate pressures between 0.6 and 3.0 GPa, the rotation of the methyl group is hindered, while, for pressures above 3.0 GPa, the orientation of the methyl group becomes fixed8,13 and is rotated by 45° relative to its lowtemperature, atmospheric-pressure structure.13 Recently, Sorescu and co-workers9,10 extended a rigidmolecule intermolecular interaction potential developed for CHNO materials16 to include intramolecular terms that would properly describe molecular vibrations and the deformation of NM.10 They demonstrated that this model interaction potential is well suited to predict properties of condensed phase nitromethane. Isothermal-isobaric MD simulations over a wide range of pressures and temperatures produced results for both static and dynamic properties that were in outstanding overall agreement with experimental measurements, including the prediction of the experimentally observed 45° change in methyl group orientation in the high-pressure regime relative to the lowtemperature crystal configuration.9,10
A. RxMC Method. RxMC is a classical-level computational method that simulates chemical reaction equilibria. RxMC determines the equilibrium states of reaction species by minimizing the Gibbs free energy, thereby determining the true chemical equilibrium state irrespective of the height of the activation energy barrier. Using the internal molecular partition functions for each of the reaction species, RxMC predicts shifts from the ideal-gas reaction equilibria due to nonideal effects, which are in the form of intermolecular interactions. Intermolecular contributions are included for the molecular species that are present in the reaction mixture, but other nonideal contributions can also be included.17-26 The internal molecular partition functions for each reaction species are specified which account for the internal modes of the isolated molecule. RxMC does not provide information about the kinetics of the reacting system, but in turn, it is not reaction rate limited and thus is not hindered by slow or fast occurring reactions. While RxMC does not require an interaction potential that mimics bond breaking and forming, it does require specifying the chemical reactions to be simulated. The RxMC method has been applied to a range of chemical reaction phenomena occurring in a variety of nonideal environments.17-31 All of these previous applications have been for typical reacting systems in which the reaction pathway from reactants to products involves the breaking or forming of chemical bonds. The reader is directed to these references along with the original papers for a more complete discussion of the method and the underlying physics.6,7,32 Only a brief summary of the RxMC method is given next. The isothermal-isobaric ensemble version of the RxMC method involves the following trial moves: (1) a change in the center-of-mass position or orientation of a molecule which is chosen at random, (2) a random change in the simulation box volume and/or shape, and (3) a randomly chosen reaction and reaction direction (forward or reverse), whereby reactant molecules are randomly chosen and deleted while product molecules either replace reactant molecules or are inserted randomly into the simulation box. Step 1 ensures thermal equilibrium is established for the specified temperature, step 2 ensures mechanical equilibrium is established for the specified pressure, and step 3 ensures chemical equilibrium is established for the specified reactions. The acceptance probabilities for steps 1 and 2 can be found elsewhere.7 For reactions in step 3, the acceptance probability of state k going to state l for reaction j is given by
{∑ cj
rxn Pj,kl ) min 1,
i)1 (Ni
Ni!
( ) qint,iV
+ νjiξj)! Λi3
νjiξj
}
exp(-β∆Ukl)
(1)
where cj is the total number of species in reaction j; Ni is the total number of molecules of species i; νji is the stoichiometric coefficient of species i in reaction j; ξj is the molecular extent of reaction for reaction j; qint,i is the internal molecular partition function of species i, which includes vibrational, rotational, and electronic modes; Λi is the thermal de Broglie wavelength of species i; V is the total volume of the system; ∆Ukl ) Ul - Uk is the change in the configurational energy; and β ) 1/kBT, where kB is the Boltzmann constant and T is the temperature. Equation 1 is appropriate for both forward and reverse reaction steps (ξj ) 1 for a forward step and ξj ) -1 for a reverse step), where the stoichiometric coefficients are taken to be positive for product species and negative for reactant species. B. The Reaction Set. The reaction pathway for the conversion of chemically distinct species does not always involve the
Reaction Ensemble Monte Carlo
J. Phys. Chem. C, Vol. 111, No. 1, 2007 367
breaking or forming of chemical bonds. For example, the antibutane S gauche-butane conformational isomerization produces chemically distinct species through a rotation around the central C2-C3 bond.34 Although no bonds are broken or formed as in the usual chemical reaction, the process is controlled by the same physical criteria. For some solid materials, pressure governs similar chemical transformations. In this work, we treat the rotational behavior of the methyl group about the C-N bond as a conformational isomerization reaction and predict the concentration of the rotamers as a function of the pressure and temperature. An atomistically detailed model of nitromethane is used, where internal rotations of the methyl group are modeled every 3° of rotation over a range of 120°. The simulated reaction species set is therefore 40 conformational isomers for this rotational range, each of which is identical in molecular structure except for the orientation of the methyl group relative to the heavy-atom frame. Rotations of the methyl group every 1° and 2° were also considered but require larger system sizes to minimize statistical uncertainty. A schematic of the reaction set used to simulate conformers of NM is given in Figure 2, where a total of 78 reactions is simulated. As illustrated in Figure 2, the reactant and product rotamers for each reaction differ by a 3° rotation. The chemical equilibria condition for this reaction set is then
µrotamer 1 - µrotamer 2 ) 0 µrotamer 2 - µrotamer 3 ) 0 µrotamer 3 - µrotamer 4 ) 0 etc.
(2)
where µi is the chemical potential of rotamer i. An alternative reaction set was also considered where the reactant is the same rotamer in all reactions, with the choice of this rotamer being inconsequential. The chemical equilibria condition for this reaction set is then
µrotamer 1 - µrotamer 2 ) 0 µrotamer 1 - µrotamer 3 ) 0 µrotamer 1 - µrotamer 4 ) 0 etc.
(3)
The results for these two reaction sets were in excellent agreement; however, the reaction set shown in eq 3 resulted in an overall lower reaction step acceptance ratio and higher statistical noise. The results reported in this work are for the reaction set shown in eq 2. IV. Model and Simulation Details A. Intermolecular Potential Models. The model potentials for the nitromethane crystal used in this work are derived from a previously determined fully flexible model potential which contains both intra- and intermolecular interactions.10 Although various physical and chemical processes cannot be simulated without an adequate description of the interatomic interactions, some processes do not require a complete description of all interactions. For example, melting simulations using only the intermolecular portion of the force field for NM produced results that were similar to those using the fully flexible potential, indicating that intramolecular motions do not strongly affect the process of melting NM.35 However, other studies have
Figure 2. Set of conformational polymorphic reactions used in the RxMC simulations to mimic the methyl group rotation in crystalline nitromethane.
shown the limitations of the models at high pressures, demonstrating that a more complete description of the interatomic interactions is required for these regimes.36 Therefore, in this work, we will attempt to identify the key features of the interatomic interactions that influence the pressure-induced methyl group reorientation by performing simulations with and without all or some of the molecular vibrations. Results using a fully flexible model, a model with reduced flexibility, and a rigid model in which the molecular vibrations are fixed will be compared. The reduced-flexible model will assume that all intramolecular degrees of freedom are fixed except for the orientation of the methyl group relative to the heavy-atom frame of NM. In other words, the only intramolecular term that will be included in this model corresponds to the methyl group rotation. Hereafter, we will term the three models, respectively, as the flexible-SRT model, the reduced-flexible-SRT model, and the rigid-SRT model. In the rigid-SRT model, the intermolecular portion of the potential uses a Buckingham exponential-6 form
{
r < rcore
∞
Uexp-6 (r) ij
)
Aij exp(-Bijr) -
Cij r
6
r g rcore
(4)
plus electrostatic interactions accounted for by a charge-charge Coulombic potential
(r) ) UCoulombic ij
qiqj 4π0r
(5)
where r is the interatomic distance between atoms i and j; Aij, Bij, and Cij are the interaction parameters for the i-j pairs; qi and qj are the electrostatic charges on the atoms; and 0 is the dielectric permittivity constant of free space. The cutoff distance (rcore) is included to avoid the unphysical singularity in the potential function as r f 0. The atom-centered partial charges used in these calculations were determined through fitting to
368 J. Phys. Chem. C, Vol. 111, No. 1, 2007
Brennan et al.
the quantum-mechanically derived electrostatic potential for an isolated molecule whose atoms are arranged in the experimental crystallographic arrangement; details are given elsewhere.10 The unlike interactions between atoms i and j were approximated by the Lorentz-Berthelot mixing rules.37 A spherical cutoff of 11.55 Å was applied with long-range corrections added to the exponential-6 terms to account for interactions beyond this distance.38 Ewald summations were employed to account for the long-range electrostatic interactions using the same spherical cutoff.33,38 The flexible-SRT model is an extension of the rigid-SRT model through the inclusion of intramolecular terms that allow for the description of molecular deformation and vibration. The intramolecular portion is taken as a superposition of bond stretching, bond bending, and torsional contributions. Force constants for the intramolecular terms were determined by adjusting the constants to reproduce eigenvectors and eigenvalues of the normal modes of vibration of the isolated molecule calculated at the B3LYP/6-31G* level.10 The reduced-flexible-SRT model excludes all intramolecular terms in the flexible-SRT model except for the cosine-type torsional potential that describes the orientations of the hydrogen atoms relative to the heavy-atom frame of the NM molecule. The form of the torsional potential is 6
U
torsional
)
UΦ [1 + cos(3Φi - δ)] ∑ i)1 i
(6)
where UΦi is the torsional barrier, Φi denotes one of the six H-C-N-O dihedral angles in NM, and δ is either 90.0 or -90.0°. Thus, for the reduced-flexible-SRT model, all bond distances and angles are fixed; the only intramolecular structural parameters that vary are the six HCNO dihedral angles upon a single reaction step in the RxMC scheme. The potential parameters for the rigid-SRT and reduced-flexible-SRT models used in this study are given in Table 1; parameters for the fully flexible-SRT model can be found elsewhere.10 The rigid structures of the rigid-SRT and reduced-rigid-SRT models are taken to be the molecular structure consistent with that of the experimental crystal at 4.2 K, 1 atm.11 For all models, rotation of the methyl group relative to the heavy-atom frame of NM will occur. However, for the reducedflexible-SRT model and the rigid-SRT model, rotations of the methyl group will occur in the form of a set of conformational isomers, each of which is identical in molecular structure except for the orientation of the methyl group relative to the heavyatom frame. Since all intramolecular contributions are ignored for the rigid-SRT model, contributions due to the orientation of the methyl group will only arise in the intermolecular interactions. B. Molecular Partition Functions. The RxMC simulation method requires inputting the internal molecular partition functions for each of the reacting species present. For polyatomic molecules, the total molecular partition function can be written as
q ) qtransqelecqvibqrotqint rot
(7)
where the contributions in eq 7 are, respectively, translational, electronic, vibrational, rotational, and internal rotational.39 The rotational contribution (qrot) refers to rotations of the molecular center of mass, while the internal rotational contribution (qint rot) refers to rotations of a group of atoms within the molecule, such as the methyl group in NM. For the set of reacting species considered in this work, qtrans, qelec, qvib, and qrot will be identical
TABLE 1: Potential Parameters for the Crystalline Nitromethane Rigid-SRT Modela Atom-Atom Exponetial-6 Potential Parameters i-j pair
Aij (kJ/mol)
Bij (Å-1)
Cij [kJ/(mol Å6)]
H-H C-C N-N O-O
9213.510 369 726.330 264 795.246 290 437.820
3.74 3.60 3.78 3.96
136.3800 2439.3459 1668.3316 1453.3114
Electrostatic Charges atom
charge (e)
C1 N2 O3 O4 H5 H6 H7
-0.305 342 0.820 603 -0.470 445 -0.484 277 0.143 365 0.155 443 0.140 652 Torsional Potential Parameters
angle
UΦi (kJ/mol)
δ (deg)
Hi-C1-N2-O3 (i ) 5,6,7) Hi-C1-N2-O4 (i ) 5,6,7)
0.270 00 0.270 00
-90.0 90.0
a
The atom indices are defined in Figure 1.
for each rotamer. This allows for a simplification of the reaction step acceptance criteria given in eq 1. For the rigid-SRT model, we ignore internal rotational contributions; therefore, the (qint,iV/ Λi3)νjiξj term in eq 1 is unity. However, for the reduced-flexibleSRT model, internal rotational contributions are included for each of the rotamers and further analysis is required. Consider the first reaction given in Figure 2, rotamer #1 f rotamer #2. The acceptance criterion is then
{
rxn P1,kl ) min 1,
( ) }
N1
qint rot,2 -β∆Ukl e (N2 + 1) qint rot,1
(8)
In the classical limit, the internal rotational partition function can be expressed as
qint rot(T) )
(2πkBTIred)1/2 nrh
∫02πe-U
tors(Φ)/k
BT
dΦ
(9)
where Ired is the reduced moment of the rotating group, nr is the number of maxima in Utors(Φ), and h is the Planck constant.40 (By considering a system in the limit of zero density, Frenkel and Smit derived a similar expression relating the intramolecular contributions of an isolated molecule and the internal partition function; see Appendix G in ref 33.) The internal rotation energy, Utors(Φ), is described by the torsional potential given in eq 6. As merely a simulation device to explore the effect of the methyl group rotation in crystalline NM, we have discretized these torsional contributions among the various rotamers based on the dihedral angle. In this work, we have defined these intervals to be the same for each rotamer, where the interval (Φinterval) is taken as the ratio of the rotational range (Θ) to the number of rotamers in that rotational range (Nrot). By additivity then, the definite integral in eq 9 can be written as the sum of the integrals for these intervals over the range of rotation
qint rot(T) )
(2πkBTIred)1/2
∑ Θ
3h
∫Φ
e-U
tors(Φ)/k
BT
dΦ
(10)
interval
where nr ) 3 in eq 9 for the methyl group in NM. In our study,
Reaction Ensemble Monte Carlo
J. Phys. Chem. C, Vol. 111, No. 1, 2007 369
for each rotamer i, we have assumed that Utors(Φi) is constant for each interval. From eq 10, we can then determine an internal rotation partition function for each rotamer i
(2πkBTIred)1/2 qint rot,i(T) ) 3h
∫Φ
)
e-U
tors(Φ )/k T i B
interval
dΦi
(2πkBTIred)1/2 -Utors(Φi)/kBT e dΦi ) Φinterval 3h (2πkBTIred)1/2 tors Φintervale-U (Φi)/kBT (11) 3h
∫
where Utors(Φi) is evaluated from eq 6 for each rotamer i. The terms outside the exponential term in eq 11 will be identical for each rotamer so that our acceptance criterion (eq 8) simplifies to
{
rxn ) min 1, P1,kl
N1
e-β(U (N2 + 1)
tors(Φ )-Utors(Φ )) 2 1
}
e-β∆Ukl
(12)
Analogous expressions for the acceptance criteria for all other reactions in the reaction set can readily be derived from eq 12. C. Simulation Details. Constant-pressure RxMC simulations were initiated from a simulation box consisting of 5 × 4 × 3 unit cells of NM (240 molecules), with the contents of each unit cell arranged in the low-temperature, ambient-pressure configuration. Standard periodic boundary conditions and the minimum image convention were used.38 Simulations were equilibrated for 7 × 106 steps after which averages of the quantities were taken over 8 × 106 steps, with the exception for simulations at 4.2 K and 1 atm where quantities were averaged over an additional 18 × 106 steps. The probability of selecting each type of RxMC step, displacement and rotation of the molecular center of mass, reaction, and simulation cell shape change was approximately 10, 99, and 1%, respectively. The maximum range for the attempted displacements, rotations, and volume changes were adjusted to achieve an acceptance fraction of approximately 0.4. The number of accepted reaction steps varied over the conditions considered but ranged from approximately 0.1 × 104 to 7 × 104 accepted reaction steps per simulation run. A scheme to adjust the reaction-type selection was implemented to ensure that the number of accepted reactions is approximately the same for each reaction type over the simulation run. This scheme is briefly explained next with further details found in the Appendix. In an RxMC simulation, each reaction type need not be selected with an equal probability. Rather, it can be selected with an arbitrary probability which, however, cannot change during the simulation. During an equilibration period, we determined these arbitrary probabilities based on the requirement that the acceptance ratios for all reaction types be approximately equal. The development of this scheme was essential, since a large number of reaction types were considered in the reaction set and the acceptance of these reaction types varied considerably. By implementing the scheme, phase space near the torsional barriers is sampled with the same precision as phase space corresponding to torsional minima. V. Discussion of Results Lattice parameters and the orientation of the methyl group relative to the heavy-atom frame calculated using the RxMC method for both the rigid-SRT and reduced-flexible-SRT-models are compared with results generated using NPT-MD simulations
and the flexible-SRT model10 at various temperatures and pressures. These results are also compared with experimental data, where available. The orientation of the methyl group can be described by three of the six HCNO dihedral angles in this system; in order to compare with the Sorescu et al.10 results, we will use the HiC1-N2-O3 (i ) 5-7) dihedral angles. The experimentally determined values of these three dihedral angles at 4.2 K and 1 atm are 151, -88, and 32°, respectively.11 Normalized cumulative distributions of the Hi-C1-N2-O3 (i ) 5-7) dihedral angles determined from molecular configurations recorded during NPT-MD trajectories under these conditions have maxima at 161, -78, and 42°, respectively.10 The measured lattice dimensions for the a, b, and c axes under the same conditions are 5.1832, 6.2357, and 8.5181 Å, respectively. The values of the lattice dimensions for the a, b, and c axes at 4.2 K and 1 atm calculated using NPT-MD and the flexible-SRT model are 5.2116, 6.3416, and 8.6464 Å, respectively. Note that, for all RxMC simulations in this study, the cell shape was allowed to change. However, the simulation cell maintained its orthorhombic shape for the duration of the simulations. Within statistical uncertainty, the interaxial cell angles remained at 90° in agreement with the NPT-MD simulations of Sorescu et al.10 A. Rigid-SRT Model. RxMC simulations were performed using the rigid-SRT model at 0.3, 3.5, 5.4, and 7.0 GPa and at T ) 293 K. Lattice dimensions for the a, b, and c axes are compared with experimental values in Table 2. The calculated lattice dimensions are in good agreement with experiment and the NPT-MD simulations under all conditions; however, at the low pressures, the orientation of the methyl group is not consistent with those obtained from NPT-MD simulations and the flexible-SRT model (Figure 3a). At 0.3 GPa, the normalized cumulative distribution of the Hi-C1-N2-O3 (i ) 5-7) dihedral angles are peaked at 191, -49, and 71°, indicating that the methyl groups are rotated 30° relative to those in the NPTMD crystal under the same conditions. At the higher pressures, however, the maxima of the distributions are in much closer agreement with those of the NPT-MD simulations. The most distinguishing features of these distributions are the intensities. Those corresponding to the RxMC simulations using the rigidSRT model suggest that the methyl groups are only partially hindered at 0.3 GPa, whereas the methyl groups in the NPTMD/flexible-SRT simulations are fixed at approximately the low-temperature, ambient-pressure configuration. This behavior is more consistent with the experimental measurements than those produced by the flexible-SRT results, since experimental evidence suggests that the methyl group is freely rotating under these conditions.13 A freely rotating methyl group would correspond to a uniform (flat) distribution. While neither the rigid-SRT model nor the flexible-SRT models are flat in Figure 3a, a comparison of the curves shows that the rigid-SRT model is flatter than the flexible-SRT model. Neither of these models reproducing a uniform distribution is attributed to a deficiency in the original model. A similar effect was observed in NPT-MD simulations at 4.2 K and 1 atm using a flexible interaction potential identical to the flexible-SRT potential described herein (denoted in Sorescu et al.10 as Mgas) except for the description of the internal methyl group rotation. In Mgas, the methyl group rotational motion was described using a six-fold torsional potential with an extremely small barrier (6 cal/mol, in agreement with the experimental value15). For comparison, the flexible-SRT and reduced-flexible-SRT models have a torsional barrier of 774
370 J. Phys. Chem. C, Vol. 111, No. 1, 2007
Brennan et al.
TABLE 2: Experimental and Predicted Lattice Parameters of Nitromethane for Pressures Ranging from 0.3 to 7.0 GPa at 293 K a (Å) P (GPa)
exptl
0.3 3.5
5.182a 4.890a 4.884a 4.94b 4.769a 4.8480f 4.7844f
5.4e 7.0 a
b (Å)
flexibleSRTc
rigidSRTd
reducedflexible-SRTd
5.2070
5.1474
5.1779
4.8603
4.8420
4.8605
4.7426
4.7410
4.7633
4.6666
4.6779
4.6969
b
c
exptl 6.259a 5.885a 5.878a 5.95b 5.777a 5.8257f 5.7392f
d
c (Å)
flexibleSRTc
rigidSRTd
reducedflexible-SRTd
6.4829
6.4659
6.4643
6.0448
6.0299
6.0266
5.9412
5.9266
5.9192
5.8733
5.8591
5.8517
e
exptl 8.645a 8.099a 8.116a 8.15b 7.958a 8.1565f 8.0868f
flexibleSRTc
rigidSRTd
reducedflexible-SRTd
8.8594
8.8017
8.8008
8.3392
8.3042
8.3017
8.1666
8.1771
8.1747
8.1162
8.0977
8.0949
a,f
Reference 13. Reference 14. Reference 10. This work. Experimental measurements were made at 5.45 GPa; likewise, the value of the imposed pressure in the NPT-MD simulations of the flexible-SRT model10 was set at 5.45 GPa. f Experimental value from high-pressure data fits.14
TABLE 3: Experimental and Predicted Lattice Parameters of Nitromethane for Temperatures Ranging from 4.2 to 250 K at 1 atm a (Å)
b (Å)
c (Å)
T (K)
exptla
flexibleSRTb
rigidSRTc
reducedflexible-SRTc
exptla
flexibleSRTb
rigidSRTc
reducedflexible-SRTc
exptla
flexibleSRTb
rigidSRTc
reducedflexible-SRTc
4.2 78 150 250
5.1832 5.1983 n/ad n/a
5.2116 5.2195 5.2298 5.2620
5.0337 5.0831 5.1305 5.1910
5.0886 5.1814 5.1876 5.2296
6.2357 6.2457 n/a n/a
6.3416 6.4032 6.4605 6.5691
6.3234 6.3763 6.4315 6.5373
6.3282 6.3688 6.4301 6.5301
8.5181 8.5640 n/a n/a
8.6464 8.7543 8.8204 8.9561
8.6068 8.6658 8.7440 8.8713
8.5760 8.6688 8.7479 8.8778
a
Reference 12. b Reference 10. c This work. d Experimental data is not available.
cal/mol. While the resulting lattice dimensions of the simulated crystal using Mgas were in reasonable agreement with experiment, the methyl groups were rotated by approximately 30° relative to the experimental orientation (consistent with the results described heretofore assuming the rigid-SRT potential). Subsequently, Sorescu et al.10 modified the description of the internal methyl group rotation to a three-fold form and increased the barrier (denoted therein as model M1 and as the flexibleSRT model in this work). NPT-MD simulations using this model produced distributions of the Hi-C1-N2-O3 (i ) 5-7) dihedral angles within 10° of the experimental values. However, this is most likely the reason behind the undesired features of the dihedral angle distributions calculated using the flexible-SRT model at 0.3 GPa. Clearly, the distinctions between the rigid-SRT and flexibleSRT models are in the intramolecular motions. The use of the rigid-SRT model in RxMC calculations is equivalent to having a model with a barrierless torsional motion. However, both models produce similar distributions at increased pressures, indicating that the torsional and other intramolecular terms do not significantly contribute to the behavior of the methyl group orientation at higher pressures. Therefore, at the higher pressures, the main influence of the methyl group rotation must be due to intermolecular terms. Implementation of the RxMC technique in such a manner has allowed us to decouple the intramolecular contributions in the models. B. Reduced-Flexible-SRT Model. In order to assess whether the main distinction between results using the rigid-SRT and flexible-SRT models is due to that single intramolecular torsional term, we have performed simulations using the reduced-flexibleSRT model, in which a single intramolecular term describing the internal rotation of the methyl group has been included. The effect of pressure on the rotamer concentration for crystalline NM determined from the RxMC method using the reduced-flexible-SRT model at 293 K is shown in Figure 3b. As done previously, the rotamer concentrations are plotted in terms of the Hi-C1-N2-O3 (i ) 5-7) dihedral angles. Also shown in Figure 3b are the cumulative distributions of the dihedral angles determined from NPT-MD simulations using
the flexible-SRT model. The agreement at all pressures between the two methods is excellent. Figure 3c illustrates the effect of temperature on the methyl group rotation at 1 atm. While the RxMC and MD distributions are peaked at the same dihedral angle, the overall distributions are slightly different. Moreover, as the temperature decreases, the disagreement between the two methods increases. This could be a reflection of the neglect of other intramolecular motions that might enhance or influence methyl group libration under these conditions or perhaps as the temperature decreases longer equilibration is required in the molecular dynamics simulations. Notwithstanding, the calculated lattice dimensions for the a, b, and c axes are in reasonable agreement with experiment and the NPT-MD simulations under all conditions and are given in Table 3. These results demonstrate that the reduced-flexible-SRT model is capable of predicting structures of crystalline NM over a large temperature and pressure range. In comparison to simulations using the flexible-SRT model, this is significant, since the reduced-flexible-SRT model is considerably less expensive computationally. Furthermore, application of the RxMC approach has allowed us to decouple the intramolecular contributions in the various models, providing insight into melting and decomposition behavior by isolating the key components of the pressureinduced behavior of NM. Finally, a note regarding the notion of implementing a standard Monte Carlo (MC) approach to simulate this type of polymorphic behavior. Consider an isothermal-isobaric ensemble MC simulation of the rigid-SRT or reduced-flexibleSRT model based on the following. In addition to the standard MC moves that displace and rotate the molecular center of mass, internal rotations of the methyl group are attempted. To mimic the RxMC algorithm, attempted methyl group rotations would be constrained to a set of dihedral angles ranging from 0 to 120°. The configurational energy term (β∆U) used in the standard MC acceptance criteria33 would then include these torsional contributions. In practice, this MC move would be performed the same as a reaction step in the RxMC algorithm
Reaction Ensemble Monte Carlo
J. Phys. Chem. C, Vol. 111, No. 1, 2007 371
Figure 4. Comparisons of the dihedral angle distributions determined from NPT-MC simulations (dashed) and RxMC simulations (solid) at T ) 293 K for the reduced-flexible-SRT model at P ) 0.3 GPa.
features. The discrepancies arise because the RxMC algorithm rigorously ensures that the rotamers are in chemical equilibrium while the NPT-MC simulation does not. Note too the additional peaks from the NPT-MC simulation which are the result of the system becoming trapped in a local energy minimum during the simulation. VI. Conclusion
Figure 3. Comparison of the dihedral angle distributions at various pressures (a and b) and various temperatures (c) for the rigid-SRT model (r) using RxMC, the flexible-SRT model (f) using MD, and the reducedflexible-SRT model (rf) using RxMC. The data shown in parts a and b are at 293 K, while those in part c are at 1 atm.
used in this work; however, the acceptance criteria would differ for each approach, interestingly, only by the factor N1/(N2 + 1) given in eq 12. A comparison of the dihedral angle distributions from such an MC simulation and an RxMC simulation for the reduced-flexible-SRT model is shown in Figure 4. The distributions, while peaked at the same angle, have distinctly different
In this work, we demonstrated that the reaction ensemble Monte Carlo method can be used to efficiently predict polymorphic equilibrium. Also in this work, we presented a case in which a simple, atomistic model with reduced flexibility is capable of capturing the same behavior as a fully flexible atomistic model. This is significant, since this reduced-flexible model requires less computational resources to simulate. This is also significant because we can avoid the often daunting endeavor of generating a fully flexible model that accurately depicts the complete set of intramolecular vibrational motions. Finally, replacing a flexible model with a reduced-flexible version allows one to decouple crucial information and identify the key components of the mechanisms controlling the dynamic event being simulated. The RxMC method allows for the isolation of various intramolecular degrees of freedom in the molecule, for example, the rotational motion of a group of atoms, from which a detailed and systematic study can be carried out to identify the governing transition behavior. The advantages of using a rigid or reduced-flexible model as opposed to a fully flexible model increase substantially as the complexity of the molecule increases. A caveat to using a reduced-flexible model is that since the system is modified it can introduce nonphysical behavior. Furthermore, the RxMC method can assist in the interpretation of experimental measurements as well as in the confirmation of experimental findings. RxMC studies can determine the simplest model possible that reproduces the experimentally observed behavior, helping experimentalists to better resolve their measurements and identify key mechanisms. By applying the RxMC approach to the pressure-induced polymorphic behavior of NM, we verified the key phenomena by isolating it through the use of the reduced-flexible-SRT model. Using the RxMC method, we calculated dihedral angle distributions, which are the experimental analogues of an intensity-wave number plot determined from Raman spectra. This allowed for a direct
372 J. Phys. Chem. C, Vol. 111, No. 1, 2007
Brennan et al.
comparison with experimental measurements, a notable achievement since verification of the hypotheses arising from oftenambiguous crystallographic measurements can prove elusive. By generating a model that isolates the molecular motion that is hypothesized to govern the transition, we have shown that RxMC simulations can prove or disprove hypotheses determined from the crystallographic data. Appendix Adjusting the Selection of the Reaction Type. To increase the sampling efficiency of phase space near torsional barriers, the following generalized algorithm could be used for selecting the reaction types: (1) Loop over total # of reaction types, i ) 1 to # of reaction types (2) attempt_rxn ) attempt_rxn + rxn_try(i) (3) If (rand_rxn_step e attempt_rxn) Then (4) {attempt reaction type i} (5) End loop In the above algorithm, attempt_rxn is used as a counter and rand_rxn_step is a random number ranging from 0 and 1. The quantity rxn_try(i) controls the selection of the reaction type, where the lower the acceptance ratio for a reaction step is the larger its value should be. The value of rxn_try(i) for reaction type i is determined by the following formula:
rxn_try(i) )
1 rxn_acc_rate(rxn_type_arb) normal rxn_acc_rate(i)
(A.1)
where rxn_acc_rate(i) is the acceptance ratio of reaction type i and a normalizing factor, normal, is based on an arbitrarily chosen reaction type, rxn_type_arb, and is determined by the following formula:
normal ) # of reaction types
∑i
[
normal +
]
rxn_acc_rate(rxn_type_arb) rxn_acc_rate(i)
(A.2)
Equation A.1 ensures that reaction types with lower acceptance ratios have a higher probability of being chosen at each reaction attempt. At the start of the simulation, the acceptance ratios for each reaction type are unknown, so an equilibrium simulation is performed where for all reaction types rxn_try(i) ) (1/# of reaction types); that is, all reaction types are selected with equal probability. The equilibrium simulation should be sufficiently long so that, for all reaction types, rxn_acc_rate(i) * 0 to ensure that eqs A.1 and A.2 are defined. Note that the values of rxn_try(i) are real numbers between 0 and 1 with ∑#i of reaction types rxn_try(i) ) 1. The forward and reverse reaction directions for each reaction type should be set to the same value of rxn_try(i) to maintain symmetry in the Markov chain. The forward and reverse reaction acceptance ratios determined from the equilibrium simulation should be equivalent but will not be exactly equivalent. A simple method to ensure their exactness is to set them to a linear average of the equilibrium rxn_try(i) values. Provided that the value of rxn_try remains the same throughout the simulation run, there are no violations to microscopic reversibility. The sampling efficiency gained using the adjusted approach is illustrated in Figure 5, where the total number of accepted reaction steps for both schemes are shown using the set of conformational reactions implemented in our NM study at P )
Figure 5. Comparison of the number of accepted reactions for the adjusted and unadjusted approaches.
7.0 GPa and T ) 293 K. Data are shown for a simulation of 1.0 × 106 total attempted steps. The total number of attempted reaction steps for the adjusted and unadjusted approaches are essentially equivalent, that is, 890 081 and 890 153, respectively. For the adjusted approach, the number of accepted reaction steps is approximately the same for all reaction types, while, for the unadjusted approach, the number of accepted reactions steps varies considerably. Acknowledgment. The calculations reported in this work were performed at the Army Research Laboratory Major Shared Resource Center (Aberdeen Proving Ground, MD) and the Engineer Research and Development Center (Vicksburg, MS) Major Shared Resource Center. M.L. acknowledges the Grant Agency of the Czech Republic (Grant No. 203/05/0725) and the National Research Programme “Information Society” (Project Nos. 1ET400720507 and 1ET400720409). References and Notes (1) Cooper, P. W. ExplosiVes Engineering; Wiley-VCH: New York, 1996. (2) Callister, W. D. Materials Science and Engineering, An Introduction; Wiley: New York, 1997. (3) Poirier, J. P. Introduction to the Physics of the Earth’s Interior; Cambridge University Press: Cambridge, U.K., 1991. (4) Massa, W. Crystal Structure Determination; Springer-Verlag: Berlin, 2004. (5) Dunitz, J. D. Phase Transitions in Molecular Crystals from a Chemical Viewpoint. Pure Appl. Chem. 1991, 63, 177-185. (6) Johnson, J. K.; Panagiotopoulos, A. Z.; Gubbins, K. E. Reactive Canonical Monte Carlo - A New Simulation Technique for Reacting or Associating Fluids. Mol. Phys. 1994, 81, 717. (7) Smith, W. R.; Trˇ´ıska, B. The Reaction Ensemble Method for the Computer Simulation of Chemical and Phase Equilibria. I. Theory and Basic Examples. J. Chem. Phys. 1994, 100, 3019. (8) Courtecuisse, S.; Cansell, F.; Fabre, D.; Petitet, J. P. Phase Transitions and Chemical Transformations of Nitromethane up to 350 °C and 35 GPa. J. Chem. Phys. 1995, 102, 968. (9) Sorescu, D. C.; Rice, B. M.; Thompson, D. L. Molecular Dynamics Simulations of Liquid Nitromethane. J. Phys. Chem. A 2001, 105, 9336. (10) Sorescu, D. C.; Rice, B. M.; Thompson, D. L. Theoretical Studies of Solid Nitromethane. J. Phys. Chem. B 2000, 104, 8406. (11) Trevino, S. F.; Prince, E.; Hubbard, C. R. Refinement of the Structure of Solid Nitromethane. J. Chem. Phys. 1980, 73, 2996. (12) Trevino, S. F.; Rymes, W. H. A Study of Methyl Reorientation in Solid Nitromethane by Neutron Scattering. J. Chem. Phys. 1980, 73, 3001. (13) Cromer, D. T.; Ryan, R. R.; Schiferl, D. The Structure of Nitromethane at Pressures of 0.3 to 6.0 GPa. J. Phys. Chem. 1985, 89, 2315. (14) Yarger, F. L.; Olinger, B. Compression of Solid Nitromethane to 15 GPa at 298 K. J. Chem. Phys. 1986, 85, 1534.
Reaction Ensemble Monte Carlo (15) Tannenbaum, E.; Myers, R. J.; Gwinn, W. D. Microwave Spectra, Dipole Moment, and Barrier to Internal Rotation of CH3NO2 and CD3NO2. J. Chem. Phys. 1956, 25, 42. (16) Sorescu, D. C.; Rice, B. M.; Thompson, D. L. Intermolecular Potential for the Hexahydro-1,3,5-trinitro-1,3,5-s-triazine Crystal (RDX): A Crystal Packing, Monte Carlo, and Molecular Dynamics Study. J. Phys. Chem. B 1997, 101, 798. (17) Turner, C. H.; Johnson, J. K.; Gubbins, K. E. Effect of Confinement on Chemical Reaction Equilibria: The Reactions 2NO S (NO)2 and N2 + 3H2 S 2NH3 in Carbon Micropores. J. Chem. Phys. 2001, 114, 1851. (18) Boro´wko, M.; Zago´rski, R. Chemical Equilibria in Slit-like Pores. J. Chem. Phys. 2001, 114, 5397. (19) Turner, C. H.; Pikunic, J.; Gubbins, K. E. Influence of Chemical and Physical Surface Heterogeneity on Chemical Reaction Equilibria in Carbon Micropores. Mol. Phys. 2001, 99, 1991. (20) Turner, C. H.; Brennan, J. K.; Johnson, J. K.; Gubbins, K. E. Effect of Confinement by Porous Materials on Chemical Reaction Kinetics. J. Chem. Phys. 2002, 116, 2138. (21) Turner, C. H.; Brennan, J. K.; Pikunic, J.; Gubbins, K. E. Simulation of Chemical Reaction Equilibria and Kinetics in Heterogeneous Carbon Micropores. Appl. Surf. Sci. 2002, 296, 1. (22) Turner, C. H.; Gubbins, K. E. Effects of Supercritical Clustering and Selective Confinement on Reaction Equilibrium: A Molecular Simulation Study of the Esterification Reaction. J. Chem. Phys. 2003, 119, 6057. (23) Lı´sal, M.; Brennan, J. K.; Smith, W. R.; Siperstein, F. R. Dual Control Cell Reaction Ensemble Molecular Dynamics: A Method for Simulations of Reactions and Adsorption in Porous Materials. J. Chem. Phys. 2004, 121, 4901. (24) Peng, X.; Wang, W.; Huang, S. Monte Carlo Simulation for Chemical Reaction Equilibrium of Ammonia Synthesis in MCM-41 Pores and Pillared Clays. Fluid Phase Equilib. 2005, 231, 138. (25) Hansen, N.; Jakobtorweihen, S.; Keil, F. J. Reactive Monte Carlo and Grand-Canonical Monte Carlo Simulations of the Propene Metathesis Reaction System. J. Chem. Phys. 2005, 122, 164705.
J. Phys. Chem. C, Vol. 111, No. 1, 2007 373 (26) Lı´sal, M.; Brennan, J. K.; Smith, W. R. Chemical Reaction Equilibrium in Nanoporous Materials: NO Dimerization Reaction in Carbon Slit Nanopores. J. Chem. Phys. 2006, 124, 064712. (27) Lı´sal, M.; Smith, W. R.; Nezbeda, I. Computer Simulation of the Thermodynamic Properties of High-Temperature Chemically-Reacting Plasmas. J. Chem. Phys. 2000, 113, 4885. (28) Lı´sal, M.; Smith, W. R.; Buresˇ, M.; Vacek, V.; Navra´til, J. REMC Computer Simulation of the Thermodynamic Properties of Argon and Air Plasmas. Mol. Phys. 2002, 100, 2487. (29) Brennan, J. K.; Rice, B. M. Molecular Simulation of Shocked Materials Using the Reactive Monte Carlo Method. Phys. ReV. E 2002, 66, 021105. (30) Brennan, J. K.; Rice, B. M. Efficient Determination of Hugoniot States Using Classical Molecular Simulation Techniques. Mol. Phys. 2003, 101, 3309. (31) Turner, C. H. Monte Carlo Simulation of Equilibrium Reactions at Vapor-Liquid Interfaces. J. Phys. Chem. B 2005, 109, 23588. (32) Johnson, J. K. Reactive Canonical Monte Carlo. AdV. Chem. Phys. 1999, 105, 461. (33) Frenkel, D.; Smit, B. Understanding Molecular Simulation; Academic Press: San Diego, CA, 2002. (34) Vollhardt, K. P. C. Organic Chemistry; Freeman: New York, 1987. (35) Agrawal, P. M.; Rice, B. M.; Thompson, D. L. Molecular Dynamics Study of the Melting of Nitromethane. J. Chem. Phys. 2003, 119, 9617. (36) Sorescu, D. C.; Rice, B. M.; Thompson, D. L. Theoretical Studies of the Hydrostatic Compression of RDX, HMX, HNIW, and PETN Crystals. J. Phys. Chem. B 1999, 103, 6783. (37) Reed, T. M.; Gubbins, K. E. Applied Statistical Mechanics; McGraw-Hill: New York, 1973. (38) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: Oxford, U.K., 1987. (39) McQuarrie, D. A. Statistical Mechanics; Harper: New York, 1976. (40) Frenkel, M. Thermochemistry and Equilibria of Organic Compounds; VCH Publishers: New York, 1993.
