Pritam Kumar Panda

Since I'm still new to LAMMPS, I just want to make sure of the whole concept: What should my LAMMPS input file contain so that its results yield a reaction? And how do I link these results to ReaxFF? Everything sounds too confusing at first glance.

Christina Algemayel

I have the same question as you did. But I can no longer see the answer Pritam wrote. So could you please share the references with me on how to calculate rate constants?

Thank you.

Xiao Zhao

Unfortunately, I am unable to see the links anymore, and it has been such a long time. However, some notable papers I referred to in my research:

1)

https://www.researchgate.net/publication/322328310_Molecular_dynamics_simulation_of_the_high-temperature_pyrolysis_of_methylcyclohexane 2) ReaxFF molecular dynamics study on the pyrolysis process of cyclohexanone by Arvelos et al. 3) Reactive molecular dynamics simulation and chemical kinetic modeling of pyrolysis and combustion of n-dodecane by Wang et al.

These explain how to navigate with ReaxFF. In addition, they elaborate on how to use the number of molecules versus time (like the concentration profile) to obtain k-values for reactions.

Plotting the natural log of these k-values versus the inverse of temperature, ln(k) vs 1/T, will provide you with the Arrhenius parameters.

All the best.

Thanks a lot! You are my angel!

I read the papers you recommended and they are really helpful.

By the way, I found the rate constant is the overall rate constant in the first order reaction mechanism. What if I want to calculate rate constant for each free radical reaction ?

Xiao Zhao

Most welcome :)

For free radical reactions, you have two options:

1) the free, open software ChemTraYzer (Chemical Trajectory Analyzer) available at https://sourceforge.net/projects/chemtrayzer/

It provides you with 'events' or elementary reaction steps, in addition to their rate constants, by applying transition-state theory. More on: https://www.ltt.rwth-aachen.de/cms/LTT/Forschung/Forschung-am-LTT/Model-Based-Fuel-Design/Aktuelle-Projekte/~kqbf/ChemTraYzer/lidx/1

2) alternatively, you can employ DFT to find your transition states and elementary reaction kinetics via the Gibbs free energy and transition states. There are numerous papers that use DFT so the literature here is more abundant. Note that DFT works on a smaller scale than molecular dynamics as it is based on quantum chemistry.

Christina Algemayel

Thanks for your recommedation!

I spent some time with ChemTraYzer. Since it's originally built with python2, it costs me several days to get it working on the latest linux platform.

Finally, I followed the example of wiki website of ChemTraYzer. The output file are Default.rate.tab, Default.reac.tab and Default.spec.tab.

What should I do with those files as a beginner?

May you have a nice day!🙂

Xiao Zhao

Glad to hear ChemTraYzer is up and running. For some reason, I am unable to visualize the files you attached, but I'll just explain from the files I already have as I'm assuming all the formats are the same.

First of all, export the .tab to excel or notepad++ to visualize them better, then:

1. for the reac.tab file:

• the first row is like a header, it gives you the first column as the time steps and the other columns as the indices of your reactions (R1, R2...)
• the second row defines your reactions based on species indices (something like S1-->S2)
• the third row does the same but with the actual species or 'sum formulas' (rather than S1--> S2 it'll tell you A-->B)

2. for the spec.tab:

• the first row gives you time steps and species indices
• the second row gives you the 'SMILES' of the species. SMILES, or Simplified Molecular Input Line Entry System, is like a representation to let the computer understand the chemical species. you can read more on the notations here: https://archive.epa.gov/med/med_archive_03/web/html/smiles.html#:~:text=SMILES%20(Simplified%20Molecular%20Input%20Line,learn%20a%20handful%20of%20rules

or convert from SMILES notation to chemical structure here: https://pubchem.ncbi.nlm.nih.gov//edit3/index.html

• the third row and beyond give you the timestep with its counts of the species. You know the count from the index: if you have at time step 120s: 0 0 0 3 0 0 1, then S1 is 0, S2 is 0...S4 is 3 and S7 is 1. the row will look like: 120 0 0 0 3 0 0 1. Yes, a bit confusing!

3. rate.tab (what you need)

• first three rows like reac.tab. note: '*' = back reaction
• fourth and beyond: timestep followed by rate constant of the reaction R1; rate constant of the back reaction R1*; and so on for R2, R3....

I believe spec and reac tab files are the most relevant and what you'll rely most on.

Kindly let me know if you need additional information.

Thanks for your detailed explanation Christina.

The attatched file is my output from reaxff run. When I use CTY to analyzer it, my command is like,

python2 processing.py H2O2+Ag/reaxff_out.bonds -i 1:1 -i 2:8 -i 3:47 -norec 5

I already added ".Val" and ".Cor" for silver in processing.py. And I got the following error,

Processing files ...

H2O2+Ag/reaxff_out.bonds reading all steps, staring at timestep 0,

using 10000 bytes memory, a static cutoff of 0.5, considering 5

steps for recrossing and reading each 1th step of connectivity.

Identifying 1: 1, 2: 8, 3: 47

/home/xiao/Softwares/anaconda3/envs/py2/lib/python2.7/site-packages/openbabel/__init__.py:14: UserWarning: "import openbabel" is deprecated, instead use "from openbabel import openbabel"

warnings.warn('"import openbabel" is deprecated, instead use "from openbabel import openbabel"')

Traceback (most recent call last):

File "/home/xiao/Softwares/ChemTraYzer_2.1/processing.py", line 1210, in

proc = Processing(Reader=reader, Identify=identify[f], Start=worktime[f], Storage=storage[f], Static=static[f])

File "/home/xiao/Softwares/ChemTraYzer_2.1/processing.py", line 211, in __init__

self.reaction, self.atom, self.molecule, self.bond[not self.switch] = self.initSystem()

File "/home/xiao/Softwares/ChemTraYzer_2.1/processing.py", line 443, in initSystem

mol[1] = self.canonical(mol[0],atoms)

File "/home/xiao/Softwares/ChemTraYzer_2.1/processing.py", line 791, in canonical

bond.SetBO(1)

AttributeError: 'OBBond' object has no attribute 'SetBO'

I searched online but got no clue on how to solve this problem. Could you give a try on my out file in your computer? Thanks in advance!

Never mind. I replaced the default 3.0.0 openbabel with 2.4.1 and this problem was fixed.

Xiao Zhao I'm glad to hear that! Keep me posted.

Hello Christina and Xiao,

I currently interested in try this software but I am not simulation a molecular reaction, instead I am working in reactions of nanoparticles with surrounding molecules. Do you think this can be applied for my system? Also Xiao, may I ask how do you solve the issues of Python2 in order to use Chemtrazer? Thank you both :)

Hi Diego Chaparro ,

I suppose it would work as ReaxFF has a wide range of applications, such as reactions occurring on surfaces (slabs), but would you care to elaborate more on your system?

Hi Christina,

Thanks for your reply, my system consists in a metal nanoparticle which reacts with gas molecules (nanoparticle in the center surrounded by the molcules inside a cube). I would like to try Chemtrazer and see if it possible to quantify the reaction rate. Do you think it would work or this software is only for molecular reaction rather than surfaces?

Diego Chaparro

I believe it would work in ReaxFF, as you can build such a system: a simulation box of a slab (like Nickel or any metal), with your molecules randomly placed around it yet confined inside the cube whose dimensions (and density) you can control.

However, I would suggest you check the literature first to see (1) if indeed ReaxFF has been used before for a system similar to yours and (2) what are the common methodologies of building such a system

Diego Chaparro

Your system sounds very common in this field. If you can find a proper reaxff potential file, then you can definitely simulate it. Although the results cannot be guaranteed. Because now I am having troubles with the simulatioin structure.

For python2 with chemtrayzer, you can build a virtual environment with anaconda in linux (I think it is the easiest way for newer linux computer).

Good luck!