Molecular Simulation of Chemical Reaction Equilibrium by Computationally Efficient Free Energy Minimization.

The molecular simulation of chemical reaction equilibrium (CRE) is a challenging and important problem of broad applicability in chemistry and chemical engineering. The primary molecular-based approach for solving this problem has been the reaction ensemble Monte Carlo (REMC) algorithm [Turner et al. Molec. Simulation2008, 34, (2), 119-146], based on classical force-field methodology. In spite of the vast improvements in computer hardware and software since its original development almost 25 years ago, its more widespread application is impeded by its computational inefficiency. A fundamental problem is that its MC basis inhibits the implementation of significant parallelization, and its successful implementation often requires system-specific tailoring and the incorporation of special MC approaches such as replica exchange, expanded ensemble, umbrella sampling, configurational bias, and continuous fractional component methodologies. We describe herein a novel CRE algorithm (reaction ensemble molecular dynamics, ReMD) that exploits modern computer hardware and software capabilities, and which can be straightforwardly implemented for systems of arbitrary size and complexity by exploiting the parallel computing methodology incorporated within many MD software packages (herein, we use GROMACS for illustrative purposes). The ReMD algorithm utilizes these features in the context of a macroscopically inspired and generally applicable free energy minimization approach based on the iterative approximation of the system Gibbs free energy function by a mathematically simple convex ideal solution model using the composition at each iteration as a reference state. Finally, we additionally describe a simple and computationally efficient a posteriori method to estimate the equilibrium concentrations of species present in very small amounts relative to others in the primary calculation. To demonstrate the algorithm, we show its application to two classic example systems considered previously in the literature: the N2-O2-NO system and the ammonia synthesis system.

Introduction

The chemical reaction equilibrium (CRE) problem, most typically at a specified absolute temperature T and pressure P, is of considerable importance in chemistry and chemical engineering. It can be posed as the solution of the nonlinear optimization problem of minimizing the system Gibbs function subject to element abundance conservation constraints, and macroscopically based thermodynamic models are widely used for this purpose (see, for example, Leal et al.[1]).

The CRE problem has been approached by means of molecular simulation in three different ways, each of which is based on a different molecular model of the system. When classical force fields (CFFs) are used to model the molecular environment, the reaction ensemble Monte Carlo (REMC) algorithm[2−4] provides a numerically exact solution of the problem. The algorithm is essentially a direct translation of the macroscopic thermodynamic approach to the corresponding molecular environment. It requires a specified list of species (which may be expanded to include species present in minor amounts by means of an a posteriori approach described in Section ), and implicitly incorporates the nonideal (residual) portions of the chemical potentials by molecular simulation involving the species force fields (FFs). Standard state ideal-gas quantities μ0(T, P) or ideal-gas reaction free energy changes, Δ0(T; P0) are also required input information, both of which may be separately obtained either from thermochemical tables, e.g., the JANAF Thermochemical Tables,[5] or from partition functions calculated using quantum mechanical software. We emphasize that such quantities are not required in the simulation of phase equilibria, since they cancel when the total interphase species chemical potentials are equated. REMC has been applied to many systems (see Turner et al.[4] for a review).

A second methodology describes the system using a quantum mechanical (QM) model involving its atoms. The recently introduced “first-principles” RxFPMC algorithm of Fetisov et al.[6] solves the CRE problem by means of an approximate QM simulation of the atoms of the system and an empirically based cluster analysis of its equilibrated atomic structure using a molecular identification criterion. Standard chemical potential data are not required. The original applications[6] have thus far been applied to a relatively simple system involving N and O atoms. It is very computationally expensive, and its extension to larger molecules and their structural isomers will require more detailed identification criteria for particular species, somewhat analogously to the use of a “species list”.

A third methodology describes the molecular environment by means of a “reactive force field” (ReaxFF) approach,[7−9] which uses force fields (whose parameters are typically obtained by matching to QM data), to model chemical reactions by means of the formation and dissociation of chemical bonds. van der Waals and Coulomb interaction terms are also included, and the ReaxFFs can model reaction kinetics in addition to reaction equilibria when implemented within a molecular dynamics (MD) simulation. This approach also does not require standard chemical potential data, and its underlying molecular description thus lies intermediate between the CFF and QM methods. FF parameter sets have thus far been developed primarily for atoms involved in combustion (C/B/N/H/O) systems and in aqueous (Ni/C/H/O) systems,[8] with some extensions to systems involving other atomic species (see Senftle et al.[8] for a review).

We note that, as for macroscopic CRE calculations using sets of chemical reactions to implement mass conservation and contrary to statements in the literature (e.g., Fetisov et al.[6]), REMC does not require a prescribed set of chemical reactions. As described in Smith and Missen,[10] an equilibrium algorithm based on the use of chemical reactions incorporates the atom conservation constraints by calculating an appropriate stoichiometric matrix from the formula matrix , and recalculating it as necessary during the course of the algorithm. This is a nontrivial task for multireaction systems, but for simple systems described by a single reaction (such as the examples considered herein), the stoichiometry is usually obvious, and it does not change during the course of the algorithm. In such cases, the reaction may be provided to the equilibrium algorithm at the outset.

It has also been stated[6] that “RxFPMC simulations are capable of generating new molecules, whereas the REMC (sic) simulations are limited to only those species appearing in the reaction set”. We show in this paper to the contrary that species not considered in the initial species list may be considered by means of a simple a posteriori procedure, which may be viewed as the molecular counterpart of an approach implemented by macroscopic algorithms.

A common disadvantage of all currently available molecular-simulation-based CRE algorithms for their routine implementation is their computational complexity. Although REMC[4] is arguably the currently most general and useful algorithm for the simulation of a broad range of CRE problems, its routine implementation for systems involving complex molecules, strong intermolecular interactions, high densities, and low temperatures is not easily formulated, since its computational efficiency typically depends on specific features of the molecules involved. REMC functions by implementing sequences of “reaction moves” interspersed with equilibration sequences; this process is carried out until numerical convergence is deemed to have been achieved. The major computational bottleneck lies in the efficient implementation of the reaction moves, which are equivalent to alchemical changes involving reactant and product species according to prescribed probabilities. For all but the simplest of systems, this requires the incorporation of special MC techniques, such as replica exchange,[11,12] configurational bias,[13,14] continuous fractional component,[15,16] and related[17,18] approaches.

The goal of this paper is to describe a novel classical FF-based CRE algorithm (ReMD), which is both computationally efficient and can be routinely applied to a wide range of systems, and whose implementation does not depend on the specific nature of the molecules involved. The algorithm is based on the core use of modern molecular dynamics (MD) simulation technology, which has undergone vast improvements in computational efficiency since the original development of the REMC algorithm, and which are currently implemented in widely available MD simulation packages such as GROMACS.[19] The relevant improvements for our purposes are due to both the widespread availability of increased hardware computing speeds (roughly 2 orders of magnitude since 1994[20]) and utilization of the software parallelization capabilities inherent in MD algorithms.[21−23] The key feature of the proposed ReMD algorithm is the coupling of these features with a computationally efficient free energy minimization strategy based on classical thermodynamic considerations.

In the next section of the paper, we describe the ReMD algorithm. The subsequent section describes the example systems studied, and the following section describes the simulation details and protocols. For illustrative purposes, we consider two classic single-phase CRE simulation problems: the N2–O2–NO system[2,6,24] and the ammonia synthesis system,[16,25−27] leaving application to more complex problems to future work. A Results and Discussion section follows, with a final section containing our conclusions.

ReMD Algorithm

We describe the algorithm for a reacting system at specified overall thermodynamic variables (T, P), involving a single fluid phase using mole fraction concentration variables; extensions to other thermodynamic and composition variables, and multiple phases, are straightforward in principle.

The ReMD algorithm is based on the translation to the molecular-simulation context of physical chemistry concepts underlying the macroscopic thermodynamic formulation of the CRE problem,[10] which we first briefly review. The CRE problem at specified (T, P) seeks to minimize the total system Gibbs free energy G(T, P; ) in a closed system, given bywhere is the species mole number vector with entries N giving the molar amount of species i, Ns is the number of species, μ is the chemical potential of species i, is the system mole fraction vector, with entries x = N/N, and N is the total number of moles.

The atom conservation constraints may be incorporated by means of a stoichiometric matrix that can be interpreted as describing a set of chemical reactions, and the mole number variables can be expressed aswhere R is the maximum number of linearly independent chemical reactions, and ν is the stoichiometric coefficient of species i in reaction j (we use the convention ν < 0 for a reactant species and ν > 0 for a product species). We emphasize that any reaction set may be used, whose sole purpose is to incorporate the element conservation equations of the system (see Smith and Missen,[10] Chapter 2). N0 is an arbitrary set of molar amounts satisfying the constraints, and ξ is a reaction extent variable. The necessary condition for the minimum of G in terms of the ξ variables is

The molar species chemical potentials calculated using molecular simulation can be expressed in the following form:where μ0(T; P0) is the ideal-gas chemical potential of species i at T and the standard state pressure P0 (typically 1 bar). μres,(T, P; ) is the residual chemical potential with respect to the underlying (T, P) ideal-gas model in the NPT ensemble. We have explicitly incorporated the fact that the reference state composition is x0 = 1, in view of the following development.

The basis of the ReMD algorithm is the iterative construction of ideal solution approximations to the system Gibbs function at a sequence of reference state compositions, which are obtained by the rapid solution of simple ideal solution CRE problems and which ultimately converge to the solution of eq .

We construct each ideal solution approximation by first using eq to write an exact expression for each chemical potential at an arbitrary composition in terms of its value at a composition , which will be used in the algorithm as a reference state composition:whereThe ideal solution approximation is constructed by neglecting the bracketed term in eq , to give

At each iteration, whose composition is denoted by , the ReMD algorithm solves the CRE problem for this ideal solution approximation to yield the composition at the next iteration. Solving this problem is a relatively simple task, even in the case of a multireaction system, for which many computationally efficient algorithms are available.[1] For the example problems considered in this paper, it involves finding the solution of a simple polynomial equation. The iterations are repeated until convergence is achieved.

We can examine the quality of the ideal solution approximation of eq by considering deviations of the values of G and the components of its gradient ∂G/∂ξ from the corresponding exact values at the reference state , which are given by

These equations indicate that the ReMD algorithm uses a convex approximation (the ideal solution model) that matches the exact values of both G(ξ) and its gradient at each iteration . In comparison, a Newton–Raphson strategy would also match the exact values of G and its gradient at each iteration, but would use a linear approximation to G(ξ) beyond . The general strategy of approximating a complex convex function by a sequence of simpler convex functions was first suggested by Folkman and Shapiro.[28] Its implementation in the context of a macroscopic CRE algorithm using eq was applied to an aqueous electrolyte system of 12 species and 5 elements (including the atomic charge) in Example 7.1 of Smith and Missen.[10]

Minor Species Concentrations

Provided that their concentrations are relatively small (“minor species”), the ReMD algorithm can accurately calculate the concentrations of species not included in the original list (thereby “creating new molecules”) using a relatively simple a posteriori procedure by means of a macroscopic thermodynamic device described in Section 9.5 of Smith and Missen.[10] This can also be viewed as a generalization of a somewhat related approach recently used by several research groups for solubility simulations of sparingly soluble solids.[29−32] The amounts and chemical potentials of the major species at the equilibrium composition calculated by omitting the minor species will be relatively unchanged from their values when the additional species are added to the system. In such cases, we calculate μres,(T, P; ) for a single molecule of the minor species m in the equilibrium mixture, and use eqs and 7 to approximate its chemical potential. We then write the chemical reaction forming 1 mol of the minor species from the major species, and setting its Δ value to zero yields the result

For each minor species, eq is used in conjunction with the reaction forming the species from the system’s major species to calculate its equilibrium composition. The assumption that the amounts of the hypothesized minor species are small may be tested following the implementation of the indicated procedure, and if a species amount is not small, it is added to the species list and the equilibrium composition quickly recalculated using the equilibrium composition of the major species as the initial estimate.

Example Systems

We consider two classic example reacting systems previously considered by others.

(1) The following reacting system at 3000 K and 30 GPawas one of the first systems studied by the REMC algorithm.[2] It was originally considered by Shaw,[24] and Fetisov et al.[6] recently considered it using their RxFPMC approach. We consider here calculations at several overall N:O ratios; Shaw,[24] and Smith and Tríska used EXP-6 FF models for the species which we employ here (see the original papers for the relevant force-field parameters). Ideal-gas standard chemical potentials for all species considered are taken from the JANAF thermochemical tables,[5] yieldingfor reaction . In addition to simulating the equilibrium of the N2–O2–NO system, we also calculate the compositions of the minor species {NO2, N2O, N, O, O3}. The FFs for these species are taken from Fried et al.[33]

(2) The ammonia synthesis reactionhas been considered by Turner et al.,[25] Lisal et al.,[27] and Poursaeidesfahani et al.[16] All three studies used the same Lennard-Jones-based FFs, which we also use herein (see the papers for the parameters). We consider calculations at several isotherms as a function of pressure for the stoichiometric N:H ratio, using ideal-gas data from the JANAF thermochemical tables.[5] Numerical values of Δr0 for reaction are given in the second row of Table .

Table 3

Ideal-Gas Standard Free Energy Changes Δr0(T; P0) in kJ mol–1 and Equilibrium Concentrations of Species for the Ammonia Synthesis Reaction of Eq at the Indicated Temperatures and Pressures Using the Methodology of This Work [Δr0(T) = 2μ0(NH3) – μ0(N2) – 3μ0(H2)]

T (K)	573	673	773	873
ΔGr0	25.756	48.084	70.982	94.348

Simulation Details

We used GROMACS 2016.3 for all simulations. For the EXP-6 FF, we used the following λ-scaling scheme[34] within GROMACS for implementation of the BAR method:where r is the interparticle distance, and n = 5.5, m = 1.0, and γ = 0.9. {ϵ, α, r0} are the FF parameters, whose values are given in the original papers.[24,33]

We used between 400 and 408 particles for the N2–O2–NO system, depending on the N:O ratio. For the ammonia synthesis system, the total number of particles used is equivalent to 250 particles of pure NH3.

The following general procedure was used to implement the ReMD algorithm:

The initial composition estimate was obtained by means of an ideal-gas equilibrium calculation (this is for convenience only; any other estimate can be used).

At each iteration, the residual chemical potentials are simultaneously calculated (in parallel) for all species using the BAR method,[35] using parallel calculations over 20 equally spaced λ intervals.

The ideal solution approximation of eqs and 7 is used in eq , which is solved for the extent of reaction; the particle numbers on the next iteration are then obtained from eq .

When the iterations are deemed to have converged, a final iteration is performed using 40 equally spaced λ intervals. The average of a set of 10 runs is used to calculate the residual chemical potentials and their uncertainties (one standard deviation), which are used to obtain a final equilibrium composition from the ideal solution approximation at the final state. The uncertainties (one standard deviation) in the mole fractions are calculated by the method described in the Supporting Information.

Results and Discussion

N2–O2–NO System

For the N2–O2–NO system, Step 3 of the algorithm solves the following equation for ξeq:The iterations of the algorithm for several N:O atomic ratios are shown numerically in Table . Convergence is seen to be very rapid, requiring at most two iterations, followed by a final iteration obtained from more precise values of the residual chemical potentials. The equilibrium compositions in the final rows for each ratio are in good agreement with previously obtained results for this system shown in parentheses[2,24] and in slightly poorer agreement with the RxFPMC results of Fetisov et al.[6]

Table 1

Convergence of the ReMD Algorithma for the N2–O2–NO System at T = 3000 K and P = 30 GPa Discussed in the Text at Several N:O Atomic Ratios Using the Indicated Total Number of Particles

(k)	N(N2)	N(O2)	N(NO)	x(NO)	μres,NPT(N2)	μres,NPT(O2)	μres,NPT(NO)	ΔGr	ξ
N:O = 1:1, 400 Particles
0	188	188	24	0.060	546.2631	498.0238	476.7631	–88.5779	42.731
1	145	145	110	0.275	545.8922	497.7521	475.6730	–1.2067	0.967
2	144	144	112	0.280	545.0258	497.7635	475.1624	–1.3683	0.096
final	144	144	112	0.280	545.8126	497.4717	475.6914	0.4441	–0.356
	144.3634	144.3634	111.2967	0.278217 (0.27554,[2] 0.276822,[24] 0.2699[6])
N:O = 2:1, 402 Particles
0	257	123	22	0.055	546.6131	498.7852	476.9928	–244.9882	40.156
1	217	83	102	0.254	546.0027	498.0744	476.5227	0.1675	–0.114
final	217	83	102	0.254	545.9911	498.2218	476.4518	–0.1242	0.087
	216.9130	82.9130	102.1760	0.25427 (0.25225,[2] 0.254216,[24] 0.2298[6])
N:O = 5:1, 402 Particles
0	327	59	16	0.040	547.2317	499.4415	477.6232	–91.4339	28.461
1	299	31	72	0.179	546.6726	498.7515	477.0020	–0.0136	0.456
final	298	31	72	0.179	546.3618	498.5419	477.2726	–1.0837	0.108
	298.8916	30.8916	72.2232	0.17968 (0.17854,[2] 0.177612,[24] 0.15710[6])
N:O = 11:1, 408 Particles
0	368	28	12	0.029	547.3037	499.1228	477.7640	–92.5492	17.226
1	351	11	46	0.113	546.7541	499.2541	477.4934	–1.1489	0.253
final	351	11	46	0.113	546.9614	498.7827	478.1027	0.3462	–0.075
	351.0714	11.0714	45.8528	0.11247 (0.11256[24])

The initial composition in each case indicated as (0) corresponds to rounded integer values of the indicated ideal-gas equilibrium composition. (k) is the iteration number, n is the number of particles of the indicated species, x is its mole fraction, μres, is its residual chemical potential in eq , and Δr = 2μ(NO) – μ(N2) – μ(O2) is the free energy change for the reaction. ξ is the change in the number of particles of N2 and O2 calculated from eq of the text on each iteration. The row indicated as “final” denotes the higher precision calculations described in the text at the approximately converged equilibrium composition. The final row gives the result of the ideal solution extrapolation using these results. Subscripts on the numerical results denote the uncertainties in the indicated number of final digits. Subscripts denote the uncertainties of the indicated values. The values in parentheses in the last row are previous REMC algorithm results[2,24] and the RxFPMC(BLYP) results of Fetisov et al.[6] using 192 atoms (see the indicated reference for details).

The progress of the iterations for the representative 1:1 N:O atomic ratio is shown graphically in Figures and 2. The solid curves are not needed for the algorithm implementation, but are shown to demonstrate the accuracy of the ideal solution extrapolation of eqs and 9 that uses the initial estimate as the reference state for the ideal solution chemical potential approximation. This approximation gives results for G and Δr that are almost indistinguishable from the exact curves. A significant source of this accuracy is that the approximating functions incorporate the logarithmic form of compositional dependence of the correct functions. For example, the S-shaped Δr curves shown in Figure are generic, in the sense that their logarithmic dependence on the species compositions results in limiting values of −∞ when product species approach zero amounts (ξ′ = 0) and +∞ when reactant species approach zero amounts (ξ′ = 1).

Figure 1

Exact and ideal solution approximation for the dimensionless relative Gibbs energy for the reacting system N2 + O2 = 2NO at T = 3000 K and P = 30 GPa with an atomic ratio of N:O = 1:1 and a total of 400 particles on the scaled extent of reaction parameter, ξ′ = n(NO) /400, where n(NO) is the number of particles of NO. G(ξ′) is the Gibbs energy of the system at ξ′, and G(0) is the Gibbs energy at ξ′ = 0, corresponding to 200 particles of N2 and 200 particles of O2. The points are simulation results for directly simulated values of G(ξ′) with their indicated uncertainties (one standard deviation), and the solid curve drawn through them is shown as an aid to the eye. The dashed curve shows G(ξ′) for the ideal solution approximation to G(ξ′) of eqs and 7 using the initial estimate as the reference state, which is the equilibrium composition of an ideal-gas model of the system (ξ′ = 0.06, corresponding to 24 particles of NO and 188 particles of each of N2 and O2). The vertical dashed lines show the progress of the ReMD algorithm from the initial estimate to the first iteration. The equilibrium solution for the system is indicated by the arrow.

Figure 2

Dependence of the dimensionless reaction free energy change Δr = 2 μ(NO) – μ(N2) – μ(O2) on the scaled extent of reaction parameter, ξ′, for the system described in Figure . The solid and dashed curves show the gradients with respect to ξ′ of the functions described by the corresponding curves shown in that figure. The simulation uncertainties lie within the symbol sizes.

For the N2–O2–NO system, the RxFPMC algorithm of Fetisov et al.[6] requires the selection of an approximate DFT method for implementing the calculation of the interatomic forces, and a method for analyzing the equilibrium atomic distributions to infer that particular configurations indicate separate molecular aggregates. Using four different DFT methods and an interatomic distance criterion indicating molecular aggregates based on the location of the first minimum of the radial distribution function, they obtained equilibrium compositions reasonably compatible with our results. Their indicated uncertainties for each DFT method are calculated from 32 independent simulation runs. The spread of these results arises from the different DFT approximations and shows much larger overall uncertainties in the computed compositions; the corresponding spread of our results reflects the use of different force fields within an essentially exact calculation. Overall, our results using the Fried FF are more consistent with those obtained using the Cheetah equation of state (EOS) than with those obtained using the Shaw force field. This is consistent with the fact that the Fried Exp-6 FF parameters were obtained by fitting to the Cheetah EOS.[33]

We consider the minor species {NO2, N2O, N, O, O3}, for which EXP-6 FF parameters are available[33] and the reactions areWe calculate the equilibrium composition for each reaction using eq for each minor species keeping the residual chemical potentials of the (major) reactant species fixed at their previously calculated equilibrium values. For example, for N2O, we solve the following equation for ξ, where denotes the previously calculated equilibrium composition of the major species:

The resulting mole fractions of N, O, and O3 were very small (less than 10–4 in all cases), but those of NO2 and N2O were both near 0.05; it was deemed that these species should be included in the equilibrium calculations along with N2, O2, and NO. The equilibrium composition of the resulting system with 3 reactions was then determined, requiring one or two additional iterations in all cases. At convergence, the residual chemical potentials arising from adding one particle each of N, O, and O3 were again calculated and their equilibrium mole fractions determined. The equilibrium compositions of all species are shown in Table , and compared with the results of Fetisov et al.,[6] and with those of the Cheetah thermochemical software program, quoted in Fetisov et al.[6]

Table 2

Equilibrium Concentrationsa of Species for the N2–O2–NO System at T = 3000 K and P = 30 GPa Discussed in the Text at Several N:O Atomic Ratios, Using Different Approaches and FFs

source	x(NO)	x(N2O)	x(NO2)	x(N)	x(O)	x(O3)
N:O = 1:1
REMC[24]	0.276822
REMC[2]	0.27554
ReMD with Shaw FF	0.278217
ReMD with Fried FF[33]	0.324920
ReMD with Fried FF[33]	0.320512	0.01051	0.02131	2.88131 × 10–6	8.760813 × 10–4	6.259314 × 10–5
RxFPMC(BLYP)-96	0.2792	0.0022	0.0042
RxFPMC(BLYP)-192	0.2699	0.0032	0.0043
RxFPMC(M06)	0.2753	0.0032	0.0042
RxFPMC(rVV10)	0.2853	0.0022	0.0062
Cheetah	0.3383	0.0080	0.0061
N:O = 2:1
REMC[24]	0.254216
REMC[2]	0.25225
ReMD with Shaw FF	0.25427
ReMD with Fried FF[33]	0.294122
ReMD with Fried FF[33]	0.28568	0.01131	0.01401	3.44781 × 10–6	5.870858 × 10–4	2.62823 × 10–5
RxFPMC(BLYP)-96	0.2323	0.0021	0.0032
RxFPMC(BLYP)-192	0.2298	0.0022	0.0022
Cheetah	0.2939	0.0069	0.0028
N:O = 5:1
REMC[24]	0.177612
REMC[2]	0.17854
ReMD with Shaw[24] FF	0.17968
ReMD with Fried FF[33]	0.20068
ReMD with Fried FF[33]	0.19242	0.00871	0.00531	3.81421 × 10–6	3.046515 × 10–4	5.49822 × 10–6
RxFPMC(BLYP)-96	0.1624	0.0062	0.0011
RxFPMC(BLYP)-192	0.15710	0.0001	0.0043
Cheetah	0.1892	0.0040	0.0008

The RxFPMC results are from Fetisov et al.,[6] and the Cheetah results are obtained from the Cheetah 8.0 thermochemical software package as quoted in Fetisov et al.,[6] which employs an equation of state whose pure fluids are tailored to the EXP-6 FF and dipolar interactions between molecular gas products.

We remark in passing that the ability of the RxFPMC algorithm to detect the presence of minor species in a CRE calculation is limited by the length of the simulation run. The RxFPMC simulations of Fetisov et al.[6] involved 96 or 192 atoms, and they averaged their results over 32 independent simulations. This only permitted concentrations down to 0.001 to be detected. The ReMD approximation method for minor species has no such limitations.

Ammonia Synthesis Reaction System

For the ammonia synthesis reaction of eq , the ideal-gas values Δr0(T) are given in the upper portion of Table . Our simulation results are based on data from the JANAF thermochemical tables.[5] At the considered isotherms, the Δ0(T) values were calculated using polynomial interpolation from values of Δf(NH3) given in these tables. Our values corresponding to those of Turner et al.[25] were calculated according to the description in their paper by adding 0.7 kJ mol–1 to our values, and those of Poursaeidesfahani et al.[16] were calculated from the data in Tables S2 and S3 of their paper. In comparison with our values, those of Turner et al. are slightly larger, and those of Poursaeidesfahani et al. are slightly smaller.

Step 3 of the algorithm solves the following equation for ξeq:

Our simulations converged in at most 3 iterations, similarly as in the case of the N2–O2–NO system. Numerical results are given in Table at several isotherms as a function of pressure in the case of the stoichiometric N:H ratio, and are compared in Figure with those of other research groups. All three sets of results are very similar. Our values are seen to lie between those of Turner et al.[25] and Poursaeidesfahani et al.,[16] consistent with the relative values of the ideal-gas quantities Δr0(T) shown in Table .

Figure 3

Dependence of the equilibrium composition of the ammonia synthesis reaction on pressure for the indicated isotherms and an overall N:H stoichiometric ratio. The filled symbols are the results of this work; the open circles are the results of Turner et al.,[25] and the open triangles are the results of Poursaeidesfahani et al.[16] The dashed lines are the results of an experimentally based equation of state.[36] The simulation uncertainties lie within the symbol sizes.

Conclusions

We have introduced a new molecular-simulation algorithm for calculating chemical reaction equilibrium, which we call the reaction ensemble molecular dynamics (ReMD) algorithm. We have illustrated its use in the case of specified (T, P) problems and provided examples. It minimizes the system Gibbs free energy based on iterations involving calculations of species residual chemical potentials at a sequence of reference state composition, and their subsequent extrapolation using a macroscopic ideal solution chemical potential model, whose equilibrium composition may be rapidly calculated to yield the next iteration.

The species residual chemical potentials are calculated by means of the BAR method within a standard MD simulation package; we have used GROMACS 2016.3 for illustrative purposes, which implements multicore parallelized calculations for all particle trajectories. Each iteration can also be performed rapidly by parallelizing both the set of species chemical potentials and their increments over each λ window within the BAR method.

The algorithm is general and can be readily implemented for any reacting system independently of the complexity of the molecules composing it; it requires only a list of chemical species and their ideal-gas free energy data, in addition to force-field models for the moleculars. The ReMD algorithm thus alleviates the difficulties encountered by the REMC algorithm in the calculation of chemical reaction equilibrium for challenging systems, which require the incorporation of special system-dependent Monte Carlo approaches[11−18] to ensure convergence.

We have also described a general procedure for calculating the equilibrium compositions of species that have not been included in the species list and are present in very small amounts (called “minor species”), which entails an a posteriori calculation of their chemical potentials due to the presence of a single molecule in the converged equilibrium system involving the major species.

Both main features of the algorithm (the iterative ideal solution extrapolation and the minor species calculation) are based on similar approaches previously used in the case of equilibrium algorithms employing macroscopic thermodynamic models.

We tested the algorithm for two classic benchmark reaction equilibrium problems at specified (T, P) considered by previous workers: the N2–O2–NO system and the ammonia synthesis system. Our results are in agreement with those of other workers for these systems, and convergence is achieved in most cases within two iterations.