ChemTraYzer2: Reactive MD Analysis¶

Building a periodic box of gas¶

To build our periodic box that contains our gas mixture, we will use the Builder tool of AMS. The Builder helps us to set up a initial stoichiometric composition and distribute our molecules in the box. We will simulate in quite a small box, so that the runtime of our simulation will be short.

Edit → Builder

Enter 15 as the length of the first lattice vector, the other vectors will adapt automatically

Enter “O2” in the field next to copies of, then select O2 (ADF) from the dropdown

Change the number of oxygen molecules from 100 to 20

Press to add a line for another species

Like before, enter “H2” and select H2: Hydrogen (ADF) from the dropdown

Change the hydrogen amount to 40

Click Generate Molecules

Close the Builder

Choosing simulation conditions¶

With the Builder tool, we already set the density of our system to \(60 / 15^3 \mathrm{molecules/{A^3}}\) which translates to \(29.5 \mathrm{kmol/{m^3}}\). So together with a very high temperature, we will create supercritical conditions, to make the combustion reactions happen faster. This will definitely not be a realistic setup, but it helps to see some reactions within a short calculation time.

Furthermore, we will use the default integration time step of 0.25 fs. With extreme conditions like in our case, the chemistry will be fast and usually smaller integration time steps are recommended to capture the molecular vibrations correctly. However, this would also increase the computation time, which is why we stick to the default value.

Model → MD or next to Task: Molecular Dynamics

Set Number of steps to 200000

Set Sample frequency to 10

Set Initial temperature to 3500 K

Setting the sample frequency to 10 means every 10th calculated step is written to disk and can be analyzed later, this translates to a time difference of 2.5 fs between frames. Decreasing this value will increase the accuracy of the analysis later, but also the disk usage.

Setting up a thermostat¶

ChemTraYzer 2 can calculate rate constants for all observed reactions, by evaluating the number of events and the concentration of reactants for each reaction. The computation assumes constant temperature (an NVT ensemble), which is why we need to set up a thermostat that balances out temperature changes once reactions start to happen.

Model → MD… → Thermostat or next to Thermostat

Add a thermostat by clicking and select NHC from the dropdown next to Thermostat:

Set the same temperature as before, which was 3500 K

Set Damping constant to 10 fs

Get reaction rate constants¶

For each reaction, ChemTraYzer 2 calculates a temperature-dependent rate constant based on the number of events and concentration of reactants during the simulation. You can see the calculated values in the upper list.

The units used here are cubic centimeters, mole, and second. So, for a bimolecular reaction the reported rate constant has a unit of \(\mathrm{cm^3*mol^{-1}*s^{-1}}\), a termolecular rate constant has a unit of \(\mathrm{cm^6*mol^{-2}*s^{-1}}\) and so on.

Let’s now compare two of the calculated rate constant to values from literature. For our temperature of 3500 K, there are no direct measurements, but we can still extrapolate a given Arrhenius equation to our temperature to get at least an idea of the right magnitude. For the reaction

we find a value of \(4.2*10^{10}\) from a 1986 review 1 and \(1.4*10^{11}\) from a theoretical high-pressure-limit rate equation 2. All rate constant units mentioned here are \(\mathrm{cm^3*mol^{-1}*s^{-1}}\), and the Arrhenius parameters are found via the NIST kinetics database. To find the respective value calculated by ChemTraYzer 2, do the following steps.

Click the search field from the upper list

Enter @Reaction Comp:O2 + H2 => to filter for reactions that consume \(\mathrm{O}_2\) and \(\mathrm{H}_2\)

In column Reaction Composition, look for the entry O2 + H2 => H + HO2

Look for the Rate Constant of this reaction

In this example, we found a value of \(1.6*10^{11}\), which is slightly higher than the theoretical high pressure value. Let’s check a second reaction

Enter @Reaction Comp:H2 + HO => in the search field

Look for the rate constant of the entry H2 + HO => H + H2O

An extrapolated literature value for this rate constant is \(1.1*10^{13}\) 3, and in this example we have calculated \(4*10^{13}\), which is again higher.

These values differ from experiment for a few reasons:

• We have hot and dense conditions, so the molecules did not have enough time to equilibrate in the beginning of the simulation and between reactions.

• The classical dynamics model ReaxFF cannot cover all quantum mechanical effects, and brings in some uncertainty.

• In this tutorial, the simulation size is small and the duration is short, so the rate constants are not representative. Rate constants become more accurate (and associated upper and lower bounds tighten) when there are more occurrences of each reaction. Our first reaction, for example, only occurred twice in the entire trajectory, which is certainly too small a population to accurately estimate its probability.

• Our literature values were extrapolated.

However, these quick calculations using ChemTraYzer2 come with very little computational costs and can be used as a first estimate and for comparison. The results will of course improve significantly if the user invests more time in a longer, larger simulation.

Visualize reactive events¶

We used the upper list, now let’s make use of the single reaction events from the lower list. Each entry there corresponds to a change in the molecular composition that happened at some point during the simulation (an event). The upper list is a condensed version of the lower list which is made by aggregating equivalent reactions.

When a reaction entry is selected, the lower list automatically shows all events of that reaction. The search bar below the events list then shows the search query that is used as a filter, e.g. (“@Index@3 45@”)

One of the reactions we expected to happen was \(\mathrm{HO_2 \rightarrow O + OH}\). Let’s use the search function to find its first occurrence.

Press the ESC key to deselect the reaction

Click the search field from the lower list

Enter @Products:[OH]

Now, only reactions with a hydroxyl radical \(\mathrm{OH}\) as a product are shown. Note here, that we have used the radical’s SMILES format for filtering (“[OH]”).

By clicking an entry of the events list, AMSmovie automatically hides all atoms that are not involved in this particular event. Now, we can take a look in detail what led to the first production of \(\mathrm{OH}\).

Switch to the ChemTraYzer 2 results window

Press the ESC key to deselect any reactions and make all atoms reappear in AMSmovie

Press the ESC key again to clear the search field

Using the “T_stable” setting¶

As mentioned in the introduction, the engine we used, ReaxFF, is a classical force field and therefore it does not include electrons or orbitals. Especially in our dense simulation, some unintuitive or unexpected molecules may appear, like \(\mathrm{HOHOH}\) or \(\mathrm{OH_3}\). These might be artifacts of imperfect force field parameters or actual short-lived, unstable intermediates.

In the latter case, we can make the criteria for stability more strict to filter out some of these short-lived species. This criteria is called T_stable and can be accessed in the CT2 menu as an analysis setting. T_stable is essentially the minimum amount of time a species must exist in order to be considered stable and appear in the results.

In this section of the tutorial, we will look for the unstable intermediate \(\mathrm{HOHOH}\) and increase \(\mathrm{t_{stable}}\) to see the effects on the reaction detection involving this intermediate.

Enter @Products:O[H]O in the search field of the lower list

Mind the number of found items

If we see \(\mathrm{HOHOH}\) in one of lists, it means that the intermediate was present in the simulation for longer than the time specified by T stable.

Wait for the analysis to finish, and look look again for \(\mathrm{HOHOH}\) in the reactions list.

Enter @Products:O[H]O in the search field of the lower list

Mind the number of found items

By increasing the time threshold of stable molecules, the number of unstable intermediates detected as reaction products should have decreased.

Investigating population properties of the trajectory¶

CT2 enables the user to access information about the dynamics of species counts, when reactions occur over the course of the trajectory, as well as many other general quantities that describe the evolution of the trajectory over time. This information is not calculated by default but can be requested on the ChemTraYzer2 settings window. We’ll go through an example illustrating how to access these properties.

Now, our output contains a third table, as shown below.

Now, click on a few of the molecules in the third table and take a look at the generated graphs. One of the graphs you’ll notice is the ChemTraYzer2 Events Per Time graph. This is a histogram of the number of reactive events and bond change events that occur in a certain bin of frames. Adjust the display style (lower left corner) to Bar and change the bin width to 2000 to generate an image like the one below.

Remember, bond change events occur any time a bond breaks or forms, even if this quickly re-forms the original molecule. Reaction events are sequences of bond change events that have not been filtered out (e.g., by the T_stable criterion). This plot provides an overview of when reactions and bond changes are occurring in a trajectory. This can be useful to determine if your reactive trajectory has reached equilibrium or if it terminated before all reactions completed. In the example above, it’s clear from the high bars near the right side of the graph that this trajectory should have been run for longer. Also, in general we expect there to be more bond change events than reactive events. Some bond change events are acceptable near the end of the trajectory, but significant numbers of reactive events in this region are an indication that the final frame of the trajectory has not reached reactive equilibrium.

Finally, we’ll take a look at some common species in this MD simulation. To do that, we’ll go to the third table and sort (decreasing) by Average Count. Clicking on the top 4 molecules in the sorted table (you may have to CTRL+click to select multiple) will produce a table like the following:

In this image, we can get a clear overview of what’s being consumed as a reactant and what’s being produced as a product. Clearly, \(\mathrm{H_2}\) and \(\mathrm{O_2}\) are reacting to form \(\mathrm{H_2O}\). The fourth species, called \(\mathrm{[OH]}\) above is an important intermediate in this reaction. Its rising population count is another indication that our trajectory is not in reactive equilibrium.

References¶

1

W. Tsang and R.F. Hampson, Chemical kinetic data base for combustion chemistry. Part I. Methane and related compounds, J. Phys. Chem. Ref. Data 15, (1986)

2

J.V. Michael et al., Initiation in H2/O2: Rate constants for H2+O2 -> H+HO2 at high temperature, Proc. Combust. Inst. 28, 1471 - 1478 (2000)

3

P. Roth and Th. Just, Kinetics of the high temperature, low concentration CH4 oxidation verified by H and O atom measurements, Symp. Int. Combust. Proc. 20 (1985)