Ctrl+K

Tutorial: Analysis#

Open in Bohrium

In this tutorial, we analyze a simple trajectroy as an example. The trajectory is taken from the real case the thermal decomposition of CL-20.

Download this trajectory using wget:

%%bash
wget https://raw.githubusercontent.com/tongzhugroup/TRAJREAX/f10a5c2cab77d3f3b659d9dd08256ae7b27c2820/cl20/cl20.lammpstrj -O cl20.lammpstrj

--2022-11-05 08:51:40--  https://raw.githubusercontent.com/tongzhugroup/TRAJREAX/f10a5c2cab77d3f3b659d9dd08256ae7b27c2820/cl20/cl20.lammpstrj
Resolving raw.githubusercontent.com (raw.githubusercontent.com)... 185.199.109.133, 185.199.110.133, 185.199.111.133, ...
Connecting to raw.githubusercontent.com (raw.githubusercontent.com)|185.199.109.133|:443... connected.
HTTP request sent, awaiting response... 200 OK
Length: 920127 (899K) [text/plain]
Saving to: ‘cl20.lammpstrj’

     0K .......... .......... .......... .......... ..........  5% 4.24M 0s
    50K .......... .......... .......... .......... .......... 11% 5.15M 0s
   100K .......... .......... .......... .......... .......... 16% 24.7M 0s
   150K .......... .......... .......... .......... .......... 22% 20.2M 0s
   200K .......... .......... .......... .......... .......... 27% 7.99M 0s
   250K .......... .......... .......... .......... .......... 33% 42.7M 0s
   300K .......... .......... .......... .......... .......... 38% 38.7M 0s
   350K .......... .......... .......... .......... .......... 44% 40.9M 0s
   400K .......... .......... .......... .......... .......... 50% 54.6M 0s
   450K .......... .......... .......... .......... .......... 55% 58.8M 0s
   500K .......... .......... .......... .......... .......... 61% 9.66M 0s
   550K .......... .......... .......... .......... .......... 66% 45.7M 0s
   600K .......... .......... .......... .......... .......... 72% 38.5M 0s
   650K .......... .......... .......... .......... .......... 77% 80.2M 0s
   700K .......... .......... .......... .......... .......... 83% 70.4M 0s
   750K .......... .......... .......... .......... .......... 89% 68.7M 0s
   800K .......... .......... .......... .......... .......... 94% 67.0M 0s
   850K .......... .......... .......... .......... ........  100% 79.9M=0.05s

2022-11-05 08:51:41 (18.5 MB/s) - ‘cl20.lammpstrj’ saved [920127/920127]

Then use reacnetgenerator to analyze it. Note the order of elements (H C N O) MUST map to those in the trajectory.

%%bash
reacnetgenerator -i cl20.lammpstrj -a H C N O --type dump --nohmm

2022-11-05 08:51:46,154 - ReacNetGenerator 1.6.6.dev24+g068dac4 - INFO: ReacNetGenerator: an automatic reaction network generator for reactive
molecular dynamics simulation.

Please cite: ReacNetGenerator: an automatic reaction network generator
for reactive molecular dynamic simulations, Phys. Chem. Chem. Phys.,
2020, 22 (2): 683-691, doi: 10.1039/C9CP05091D

Jinzhe Zeng (jinzhe.zeng@rutgers.edu), Tong Zhu (tzhu@lps.ecnu.edu.cn)

==================
Features
==================
* Processing of MD trajectory containing atomic coordinates or bond orders
* Hidden Markov Model (HMM) based noise filtering
* Isomers identifying accoarding to SMILES
* Generation of reaction network for visualization using force-directed
  algorithm
* Parallel computing

==================
Simple example
==================
ReacNetGenerator can process any kind of trajectory files containing 
atomic coordinates, e.g. a LAMMPS dump file prepared by running “dump 1
all custom 100 dump.reaxc id type x y z” in LAMMPS:
$ reacnetgenerator --type dump -i dump.reaxc -a C H O --nohmm
where C, H, and O are atomic names in the input file. Analysis report
will be generated automatically.

Also, ReacNetGenerator can process files containing bond information, 
e.g. LAMMPS bond file:
$ reacnetgenerator --type bond -i bonds.reaxc -a C H O --nohmm

You can running the following script for help:
$ reacnetgenerator -h

2022-11-05 08:51:46,154 - ReacNetGenerator 1.6.6.dev24+g068dac4 - INFO: Version: 1.6.6.dev24+g068dac4  Creation date: 2018-03-11

Read bond information and Detect molecules: 0timestep [00:00, ?timestep/s]
Read bond information and Detect molecules: 2timestep [00:01,  1.40timestep/s]
Read bond information and Detect molecules: 101timestep [00:01, 70.84timestep/s]

Save molecules:   0%|          | 0/11 [00:00<?, ?molecule/s]
Save molecules: 100%|██████████| 11/11 [00:00<00:00, 3376.31molecule/s]
2022-11-05 08:51:47,881 - ReacNetGenerator 1.6.6.dev24+g068dac4 - INFO: Step 1: Done! Time consumed (s): 1.726 (Read bond information and detect molecules)
2022-11-05 08:51:47,941 - ReacNetGenerator 1.6.6.dev24+g068dac4 - INFO: Step 2: Done! Time consumed (s): 0.060 (Merge isomers)

HMM filter: 0molecule [00:00, ?molecule/s]
HMM filter: 11molecule [00:00, 2628.61molecule/s]
2022-11-05 08:51:48,139 - ReacNetGenerator 1.6.6.dev24+g068dac4 - INFO: Step 3: Done! Time consumed (s): 0.198 (HMM filter)

Indentify isomers:   0%|          | 0/11 [00:00<?, ?molecule/s]
Indentify isomers: 100%|██████████| 11/11 [00:00<00:00, 585.62molecule/s]

Analyze atoms:   0%|          | 0/11 [00:00<?, ?molecule/s]
Analyze atoms: 100%|██████████| 11/11 [00:00<00:00, 7775.08molecule/s]

Collect reaction paths:   0%|          | 0/288 [00:00<?, ?atom/s]
Collect reaction paths: 100%|██████████| 288/288 [00:00<00:00, 11114.42atom/s]

Analyze reactions (A+B->C+D):   0%|          | 0/100 [00:00<?, ?timestep/s]
Analyze reactions (A+B->C+D): 100%|██████████| 100/100 [00:00<00:00, 12790.24timestep/s]
2022-11-05 08:51:48,583 - ReacNetGenerator 1.6.6.dev24+g068dac4 - INFO: Step 4: Done! Time consumed (s): 0.444 (Indentify isomers and collect reaction paths)
2022-11-05 08:51:48,651 - ReacNetGenerator 1.6.6.dev24+g068dac4 - INFO: Step 5: Done! Time consumed (s): 0.068 (Reaction matrix generation)
2022-11-05 08:51:48,667 - ReacNetGenerator 1.6.6.dev24+g068dac4 - INFO: Species are:
2022-11-05 08:51:48,667 - ReacNetGenerator 1.6.6.dev24+g068dac4 - INFO: 1 [H]C12N(N(=O)=O)C3([H])N(N(=O)=O)C1([H])N(N(=O)=O)C1([H])N(N(=O)=O)C3([H])N(N(=O)=O)C1([H])N2N(=O)=O
2022-11-05 08:51:48,667 - ReacNetGenerator 1.6.6.dev24+g068dac4 - INFO: 2 [H]C12NC3([H])N(N(=O)=O)C1([H])N(N(=O)=O)C1([H])N(N(=O)=O)C3([H])N(N(=O)=O)C1([H])N2N(=O)=O
2022-11-05 08:51:48,667 - ReacNetGenerator 1.6.6.dev24+g068dac4 - INFO: 3 O=NO
2022-11-05 08:51:48,672 - ReacNetGenerator 1.6.6.dev24+g068dac4 - INFO: The position of the species in the network is:
2022-11-05 08:51:48,672 - ReacNetGenerator 1.6.6.dev24+g068dac4 - INFO: {1: array([-5.95629110e-05, -1.00129074e-04]), 2: array([-0.59480018, -0.99989987]), 3: array([0.59485974, 1.        ])}
2022-11-05 08:51:48,855 - ReacNetGenerator 1.6.6.dev24+g068dac4 - INFO: Step 6: Done! Time consumed (s): 0.204 (Draw reaction network)
2022-11-05 08:51:48,878 - ReacNetGenerator 1.6.6.dev24+g068dac4 - INFO: Report is generated. Please see cl20.lammpstrj.html for more details.
2022-11-05 08:51:48,942 - ReacNetGenerator 1.6.6.dev24+g068dac4 - INFO: Step 7: Done! Time consumed (s): 0.087 (Generate analysis report)
2022-11-05 08:51:48,942 - ReacNetGenerator 1.6.6.dev24+g068dac4 - INFO: ====== Summary ======
2022-11-05 08:51:48,942 - ReacNetGenerator 1.6.6.dev24+g068dac4 - INFO: Total time(s): 2.787 s

The results are shown in cl20.lammpstrj.html. In addition, we can view an overview of reactions in the trajectory in the command line:

%%bash
cat cl20.lammpstrj.reactionabcd

1 [H]C12N(N(=O)=O)C3([H])N(N(=O)=O)C1([H])N(N(=O)=O)C1([H])N(N(=O)=O)C3([H])N(N(=O)=O)C1([H])N2N(=O)=O->O=NO+[H]C12NC3([H])N(N(=O)=O)C1([H])N(N(=O)=O)C1([H])N(N(=O)=O)C3([H])N(N(=O)=O)C1([H])N2N(=O)=O

Here we can clearly see that the CL-20 molecule is decomposited into a nitrogen dioxide molecule.

Edit on GitHub
Show Source