Managing large reaction networks

go back to Main Page, Group Pages, Núria López and Group, Scripts_for_VASP

This example is centered on the full decomposition/reaction network of ethylene glycol, that was studied by Qiang and Rodrigo during several years.

First part: Define the whole reaction network.

Take the most complex molecule and write down all possible intermediates you can imagine by removing one atom at a time.

Example: For the decomposition of ethylene glycol, you can start breaking a O-H bond, a C-H bond, a C-C bond or a C-O bond. Depending in how are you considering the conformations and the symmetry of your systems, there can be actually two possible O-H bond breakings: In the first one, the H participates in an intramolecular hydrogen bond; in the second one, the O-H is receiving an hydrogen bond, but the leaving H atom is not participating in one. 

Complete your reaction network by drawing all the lateral paths.

Example: If you remove one hydroxyl group from ethylene glycol, you end up with CH2CH2OH, which can be hydrogenated to form ethanol or another derivative (like acetaldehyde). Conversely, if you break all C-C and C-O bonds, you will end up having derivatives of methane and water. Actually, the full reaction network of ethylene glycol decomposition has overlaps with the ones of ethanol, ethane, methanol, methane, and water, considering that some hydrogenations may occur during the decomposition (carbon monoxide and hydrogen are not necessarily the unique decomposition products). If we consider that our reactions may take place even in a sligtly oxidazing atmosphere, the reaction network must be expanded to include O2 adsorption and decomposition, and oxidation of several species. 

Identify all possible products and by-products.

It is a very common mistake in theoretical scientific literature to omit possible products, which normally happens for not writting down the full reaction network. If the reaction network is not complete, all further work will describe Narnia and not the real world. So be sure to triple-check and quadruple-check all your possible products, this is, all the intermediates that may be stable in gas phase. Check experiments: Your reaction must include all products that were observed experimentally; however, some of your reported products may have not been detected if they are not very stable in the reaction conditions. 

Second Part: Define the computational setup Before doing the first calculation, define all the computational details you will use along your study, and make them as cheap as possible, without loosing accuracy. One or two weeks of test calculations may be needed. Don't be affrain of taking many days in testing: each testing day may save you several months of calculations! Some tips:

• Put a vacuum as small as possible: 8 to 10 Å will be enough if you are only interested in the energetics. If you want electronic details for some intermediates, you can further extend the vacuum, put more kpoints and increase the energy cutoff for that specific intermediate.

If you are not convinced, imagine you put 15Å of vacuum instead of 10Å. This will double your calculation time while increasing just sligthly the accuracy (0.05 eV). Because of the excesive computational time, you will be tempted to stop your calculations before full convergence, and some of your transition states may have errors as large as 0.20 to 1.50 eV. Therefore, it is fundamental to work with well-converged results, and to reach them as fast as possible, even if the results may have a slight error for reduced vacuum.

• Make sure your cell is big enough to avoid lateral interactions in "x" and "y" with the periodic molecule. Mind: Your final states will be two molecules instead of one, may the one fragments be closer to the periodic image of the other fragment than from the original one? In this sense, 3x3 cells are simply too small.
• Also, try not to make your cell too big. If you will do your calculations in MareNostrum, think that some cell sizes may paralelize better than others. For instance, 2·sqrt(3)×2·sqrt(3) cells are larger than 3×3, but parallelizes better and run faster in MareNostrum for they full k-points conditions.
• Do intensive parallelization tests, and learn how to touch NBANDS, NPAR, NSIM, and KPAR.

Third Part: Create the structures of all your intermediates.

Fourth Part: Obtain the transition states.

Define the structure of your initial states. They should be the most reactive conformation. They should also be stable (not necessarily the most stable)

If an intermediate can decompose to yield several products, use the same initial state for all these reactions, so the post-processing will be much easier. This is, if you want to monitor the behaviour of certain C-C bond (for obtaining bond orders and similar stuff), you will know that the bond is allways between atom #58 and atom #59. This will greatly simplify further scripts for post-processing

Create the structure of your final states. For these reactions were there is only one leaving atom, you can easily script this process.

Be sure of fully relax all initial and final states before calculate the NEBs, relaxing the surface, and using the full k-points set.

Prepare the NEB. Use between 4 and 6 images. Use up to 8 only if you automate all the process.

go back to Main Page, Group Pages, Núria López and Group, Scripts_for_VASP