Home

ChemTraYzer -- Reaction Models from Molecular Dynamics Simulations

ChemTraYzer (CTY) is an integrated simulation and analysis tool to study reactive molecular dynamics of a gas phase.

---

How to run the ChemTraYzer:

The main script is simulation.py which takes a configuration file as input:

• python simulation.py < example/sim.conf

Available settings are desribed in the configuration file (example/sim.conf). CTY then builds the box with packmol and runs the dynamics with LAMMPS. Reactions and species are registered and sent to an quantum mechanics engine for optimization (Gaussian by default).
Once the high level calculations are done, CTY can build the reaction mechanism, and fit Arrhenius and NASA parameters:

• python master.py . -reac default.reac -source QM

where default.reac is the reaction list output from simulation.py and QM is the folder with the QM engine output. Final output is the reaction list default.reac, the observed rate constants reac.rate.tab, time-wise species and reactive event counts reac.spec.tab and reac.reac.tab, the mechanism in ChemKin format Default.therm and Default.chem, and the mechanism in a graph format Default.gml.

It is possible to run CTY without the dynamics part, if there is a LAMMPS-style bond table as in example/bonds.tatb:

• python processing.py example/bonds.tatb -i 1:6 -i 2:1 -i 3:8 -norec 5

where the repeated -i provides the atom/element mapping and -norec 5 enables the recrossing detection for an interval of 5 data points in the bond table. CTY then creates the reaction list output (bonds.tatb.reac), which in turn can be used to calculate rate constants from the observed events and concentrations:

• python analyzing.py examples/bonds.tatb.reac

Since in this case, CTY lacks the geometry information from the dynamics, a QM optimization and thermochemistry parameter fitting is not possible.

---

How to change the QM engine or the QM submission scripts:

The method submitQM() in simulation.py creates and submits a batch job containing the input for the QM engine (Gaussian09 in our case). It uses some static variables from the Simulation class to fill in e.g. the QM method. If you change the QM engine to something else than Gaussian, you also need to change the QM output read routine harvestGaussian() in harvesting.py. Keep in mind that the hindered rotor contribution to the RRHO partition functions is adapted to Gaussian output only.

---

Extending the supported elements of ChemTraYzer

CTY can handle the elements C, H O, N, F, S, Cl, He, Ne, and Ar by default. If you have other elements in your simulation, you need to add the maximum possible valency and the maximum possible coordination for those elements. Search for ".Val" and ".Cor" in processing.py and add the new numbers accordingly. CTY uses both values to distingush between long range interactions and actual bonds.

---

Python packages that are required to run ChemTraYzer:

• lammps (see https://lammps.sandia.gov/doc/Python_shlib.html)
• openbabel (version 2.4.0)
• numpy
• scipy
• sqlite3
• optional packages for plotting concentration profiles:
• pygraphviz
• matplotlib (pyplot)

---

Anharmonic correction to rate constants from RRHO-TST

For a short tutorial on how to correct your rate constants according to ACS Omega 2020, 5, 5, 2242-2253, please read the README.doc located in the folder small_tools/ACSOmega2020. You will need the python library TAMkin, which is not required to run ChemTraYzer.

---

Bug Reports:

2016-01-07: OpenBabel ignores the aromaticity of carbon atoms in SMILES in some cases

     - OpenBabel adds too many implicit hydrogen atoms to carbon atoms which should be aromatic
     -> analyzing.py returns a warning that the balance of elements failed, although the SMILES are balanced correctly

Bug Closed:

2015-07-30: Recrossing filter fails if final reaction of recrossing event is of type A -> A

    - Reaction A -> A is removed before tracing reaction history
    -> If Reaction A -> A is part of the recrossing process (only as one of the final steps), the filter tries to remove that reaction, but it cannot find it (due to previous point)

Version 1.1:
2015-07-15: (not directly related to the code) Installation of pygraphviz for windows failed

    - a user reported problems with installing pygraphviz for windows 8
    -> 'analyzing.py' cannot be used without pygraphviz
    **SOLUTION:** pygraphviz is optional now -> 'analyzing.py' runs witout

Bug Fixes:

Version 1.1:
2015-07-08: Disturbance of atom balance in 'analyzing.py' if the initial chemical composition contains single-atom molecules.

    - 'processing.py' does not recognize initializing reactions (atoms -> molecules) for single-atom molecules
    - 'analyzing.py' computes initial chemical composition from initializing reactions
    -> single-atom molecules are missed and concentration will become negative when these molecules are consumed. This in turn leads to wrong concentration integrals and therefore to wrong rate constants.
    **SOLUTION:** Initial reactions changed so that single atoms are recognized

2015-08-31: When being initialized, Open Babel displays H and H2 molecules as empty strings

    - Problems with displaying pure hydogen molecules
    -> String representations of Reactions may have commas enclosing an empty string 
    **SOLUTION:** Identify H and H2 separately

2015-09-29: Systems without any initial bond will end with an error in the 'splitBonds' function

     - Initialization of bond handler yields empty array
     -> 'splitBonds' trys to access empty array and fails
     **SOLUTION:** 'splitBonds' checks for empty bond array to avoid abort

2016-01-29: Recorssing filter period too long if first timestep in Bond-Order output non-Zero

     - Filter period is the the product of the number of timesteps between two entries in the Bond-Order file and a user-defined number of Bond-Order steps
     - the code assumes that the first timestep is zero and simply takes the second timestep as number of steps between two steps
     -> the computed number of steps between two timesteps is way to large, thus the filter period is way too long and actual reactions may be filtered
     **SOLUTION:** Initial timestep extraction method added

2016-01-29: Depending on the ReaxFF force field the Bond-Order between the two oxygens of a peroxide group (R-O-O.) is oscillating around 1.5

     - O-O bond is oscillating between single and double bond, thus the 'open' valence of the central O is oscillating between one and zero. As a consequence, the R-O bond is periodically forming and breaking
     -> the Recrossing filter is removing these oscillator events from the output until the R-O-O. undergoes further reactions (in which case R-O-O. would appear as an extremely short-lived intermediate).
     **SOLUTION:** half-step bond orders allowed -> issue shifted to 0.25 BO before and after integer bond orders, which are presumably less important

---

Project Members:

• Enia Mudimu (admin)
• Felix Schmalz (admin)
• Florian Solbach (admin)
• Kai Leonhard (admin)
• Malte Döntgen (admin)

[Sam Afari]

Hello,
Great software. I plan to use it but do get a warning and an error below for the example file. I am using the current version 2.1. Can you help me?

WARNING
attempt to read from bonds.tatb failed.Please check whether the file is broken or does not exist. File will be ignore...

   ERROR
     No input file found

[Felix Schmalz]

Hi Sam,

thanks for your interest, I hope I can help you with a general answer as a start.
You can use ChemTraYzer's parts separately like this:
python processing.py example/bonds.tatb -i 1:6 -i 2:1 -i 3:8 -norec 5
which creates a reaction list in the same place. The -i specifies elements, e.g. -i 1:6 defines carbon as first element. With -norec 5 you can filter out recrossing reactions, where 5 is the number of steps (corresponds to 5fs in the example case). After that,
python analyzing.py examples/bonds.tatb.reac
creates a second list with rate constants.

You can also use ChemTraYzer's integrated simulation and analysis via
python simulation.py < example/sim.conf
You will need the LAMMPS python package installed, packmol and a ReaxFF parameter file (or use the one from the example). Have a look at example/sim.conf to see all the options. By default, no Gaussian QM calculations are submitted, only input files are created, because the environment might be different on other machines. If you have Gaussian output (.log),
python master.py . -reac default.reac -source QM
creates the mechanism (.therm and .chem).

You might want to re-download the source code, because I fixed a file name error there, and I added packmol exe and a ReaxFF parameter file to the example folder. Hopefully that helped you, feel free to report back.
Best, Felix

[Raj maddipati]

Excellant software! Thank you so much!

[Ning Gu]

hello, what's version of openbabel? I meet the error about the ChemTraYzer-2.1 below.

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

I thought it's version problem about openbabel，and where I can find the corresponding dependency of ChemTraYzer.

[Felix Schmalz]

Hello buendia,

we use openbabel version 2.4.0 .
other packages you need are:
numpy
scipy
lammps (see https://lammps.sandia.gov/doc/Python_shlib.html)
sqlite3
optionally for concentration profiles you'll need:
pygraphviz
matplotlib.pyplot

Best, Felix

[yacare]

Hello! Thank you very much for your hard work.
My question is：
In the reac.rate.tab file, does the "-1" represent "there are one back reaction" ?
If so, my simulation has so many back reactions. Should I increase the value of norec?

[Felix Schmalz]

Hello yacare,
thanks for appreciating!
I guess you mean the "-1" in the reac.reac.tab file. Then you are right. That means the reaction went backwards and the forward direction ("1") was detected first in the simulation. If you get a lot of subsequent back and forth events ("1" and "-1" repeatedly), you can try increasing the recrossing interval, just as you said.
However, this value is system dependent and there is no consent on how to determine it, so you need to try around a little bit.
Best, Felix

[Felix Schmalz]

Hello yacare,
after a second thought it might well be that your system actually shows a lot of back and forth reactions, and that they are not unphysical. For example, if you have equilibrium reactions at low temperatures, you can expect the behaviour you describe.
Best, Felix

[Ning Gu]

Hello Felix Scgmalz,
I found the python script is running is very slow when I deal a 10ns trajectory with 28GB.
I thought it maybe quick if we can run the script in parallel.

[George]

Hi, thank you for your work.
I would like to ask what are the units for the rate constant in the file reac.rate.tab.

[Felix Schmalz]

Hi George,
the rate constants are printed in the unit cm^3/mol/s for bimolecular reactions. Or in general, (cm^3/mol)^(n-1) /s where n is the reaction order.
Best, Felix
EDIT: I need to point out that this applies to version 2.x

[Yi Wang]

Hello,
I have installed it and I can use the processing.py program well. But when I use the analyzing.py program just as the readme file:
python analyzing.py bonds.tatb.reac -main reac -dp 0.1 -gen id
I do get an error like:

1 molecule converted
Traceback (most recent call last):
File "../analyzingi.py", line 1022, in <module>
anly.depictMolecules(Resize=resize)
File "../analyzingi.py", line 673, in depictMolecules
call(['convert','-density','1200','-trim',path+'.svg',path+'.png'])
File "/usr/lib/python2.7/subprocess.py", line 172, in call
return Popen(popenargs, *kwargs).wait()
File "/usr/lib/python2.7/subprocess.py", line 394, in init
errread, errwrite)
File "/usr/lib/python2.7/subprocess.py", line 1047, in _execute_child
raise child_exception
OSError: [Errno 2] No such file or directory</module>

I am using the current version 2.1. Can you please help me?

Many thanks in advance!

[Felix Schmalz]

Dear Yi Wang,
it seems something went wrong while converting the openbabel SVG images to PNG. I'm not sure what exactly. Could you test if you can use the "convert" command outside of ChemTraYzer?
Best, Felix

[Yi Wang]

Dear Felix,

It seems that I have solved this problem by installing all of the packages suggested in the README file in ChemTraYzer1.1.

I have another problem. If I want to analyse Manganese relevant elements, I have to add it to the valence list, but I am not sure the valence of Mn, should I add its highest valence? And what's the difference between self.Val and self.Cor, I can't find a detailed explanation about this anywhere.

Thank you so much for your answer!

Best Reagrds,
Yi Wang

[Felix Schmalz]

Dear Yi Wang,

great. I believe that "convert" is a standard unix program, but I might be wrong. It's only used to have PNG pictures in the end instead of SVGs, so it can be replaced if it causes problems in the future.
Those valence ("Val") and coordination ("Cor") values are used in ChemTraYzer's molecule detection, such that long range interactions between atoms can be filtered out and are not counted as bonded interaction. The valence limits the sum of bond orders an atom can support, and the coordination limits the number of neighbors an atom can have. E.g. for Nitrogen, the maximum valence we use is 5 and the maximum coordination is 3.
Thanks for pointing that out, I'll put a description in the docs.
For Manganese I'd use a high max valency (look at possible oxidation numbers) and try to limit the number of neighbors (".Cor") dependent on your application.

Best, Felix

[Jiuke Chen]

Dear Felix,

Thanks for your effort in this amazing software. Following your instructions, I got a same error message as Sam after I converted the python scripts using '2to3' command. Multiple warnings and an error. Should I use python2 instead of python3?

Python3 processing.py example/bonds.tatb

WARNING
attempt to read from bonds.tatb failed.Please check whether the file is broken or does not exist. File will be ignore...

ERROR
No input file found

Kind regards,
Jiuke

[Yi Wang]

This is a software written in python2, this error is caused by incompatibility between python2 and python3.

[Felix Schmalz]

Dear Jiuke,

thanks!
We plan to release a python3 version of the code in January. Regarding your error, I assume the file exists. The warning you got is then likely printed, because we used the "file" constructor to open the file, which is not there anymore in python 3.

Best, Felix

[Raj maddipati]

Dear Jiuke,
You can install Python2 environment in Anaconda and install all the software necessary for CTY. It worked for me!
Best regards,
Raj

[Stratos]

Dear Felix,

Amazing software; I have already used CTY extensively.

I am having an issue when I run the processing script for a large molecule (number of atoms > 2000). Specifically, I am getting the error "RuntimeError: maximum recursion depth exceeded in cmp" coming from function "setMolIDRecursive".

Is there any solution to this issue? Thank you for your time in advance.

Kind regards,
Stratos

[Felix Schmalz]

Hello Stratos,

first, apologies for the late reply. Second, thanks for using ChemTraYzer and finding issues like this! In Python the recursion limit is 1000 by default, but one can change it to higher numbers:

sys.setrecursionlimit(3000)

One source I found claimed 10^6 is reasonable. For your case, 3e3 is enough I guess. We'll change that in the next version. For now, you need to add the line above to processing.py (e.g. after the imports).

Thanks for the remark and best wishes,
Felix

[Raj maddipati]

[Raj maddipati]

Dear Felix,
I am performing reactive molecular dynamics of H2/O2 system. I am doing post-processing using CTY with bond-order cutoff 0.3. Can I use different bond order cutoffs for different bonds O-H, O-O etc.? I am using bonds file generated using LAMMPS. For LAMMPS simulations I have used bond order cutoff as 0.3 in control file. So, Can I just mention different bond order cutoffs in CTY will work if there is a way to do it? If not, how can I give multiple bond order cutoff values in LAMMPS and post-process in CTY? Please let me know, if there is any solution to it. Thank you!

Best Regards,
Raj

[Felix Schmalz]

Hi Raj,

The short answer is no, unfortunately. But I'll put in on the list for our code refactoring in January.

The medium answer would be to use a ReaxFF parameterization where the bond orders behave better and all atom combinations are definitely unbonded below 0.5 or so.

The long answer is that I looked at the code and seems a quite some of it would need to be changed to implement that. However if you really need it, you could add some code that takes the elements into account, when reading the bonds from the bond file. In processing.py, look for BO[i] > self.static in initSystem() and for BO[i] > Static in extractBonds(). In both contexts, there is access to the atomic number via self.num[TYPE] where TYPE is the LAMMPS element type of one atom. You would probably also split the for-loop into two, and read each timestep twice, so you can look-up the elements of ID (the first atom id) and BP (the bond partner id).

I hope I could help,
Best regards, Felix

[Raj maddipati]

Dear Felix,
Thank you so much for your elaborate answer! I have followed your medium answer by changing the Lammps subroutine, such that it takes different bond orders for different bonds. I will try to add the modifications to the CTY as you suggested. Thank you for adding this aspect in the list.

I have one more query related to 'noreac' and 'skip' flags which is not clear for me with the description given in the code. I have collected the data at every 10 steps with 0.1fs as the time step. In a paper (link below) authors have collected data at each ten steps, but they have mentioned sampling interval as 0.02ps. Can you please elaborate about the use of these two flags and whether sampling interval is related to this? That would be really helpful!

Link for the paper: https://doi.org/10.1002/jcc.25809

Kind regards,
Raj