Reaction pathway analysis¶
In the previous chapters, we have been looking at the behavior near the stable structure. In this chapter, we will analyze the reaction of substances.
The relationship between the activation energy \(E_a\) and the enthalpy of reaction \(\Delta H\) associated with the reaction X→Y. Cited from Wikipedia.
In a reaction process in which product Y is formed from reactant X, the point with the maximum energy in the reaction path is called the Transition State (TS), and the energy difference between reactant/product and the transition state is called the Activation energy \(E_a\).
According to the Arrhenius equation, the reaction rate constant \(k\) is wrriten as
where
\(A\): pre-exponential factor (a constant independent of temperature)
\(E_a\): activation energy
\(k_B\): Boltzmann constant
\(T\): temperature
Thus, the likelihood of a chemical reaction is determined by the activation energy \(E_a\).
The energy of reactants/products does not change, but the activation energy changes as the reaction path changes with the surrounding environment, i.e., the presence of a catalyst. If we can find a catalyst with lower activation energy, we can make the target reaction more likely to be produced.
Reference:
Even if reactant X and product Y are known, finding the transition state is not obvious. In this section, we will try to calculate the activation energy using the NEB (Nudged Elastic Band) method, which searches for reaction pathways when the structures before and after the reaction are known.
NEB¶
Diagram of the NEB method: After generating reaction path candidates (red line) as the initial configuration, structural relaxation is performed to search for reaction paths (blue line) that pass through transition states.
[Note] The following content is based on the NEB Tutorial in the Matlantis.
We will investigate an organic chemical reaction called the Curtius rearrangement.
In the NEB method, we usually start with preparing the structures before and after a particular chemical reaction of interest.
This time, we will write down the structures directly.
[1]:
import numpy as np
from IPython.display import Image
import matplotlib.pyplot as plt
import ase
from ase import Atoms
from ase.visualize import view
from ase.optimize import BFGS
from ase.optimize import FIRE
from ase.neb import NEB
from ase.build.rotate import minimize_rotation_and_translation
from pfp_api_client.pfp.calculators.ase_calculator import ASECalculator
from pfp_api_client.pfp.estimator import Estimator
from pfp_api_client.pfp.estimator import Estimator, EstimatorCalcMode
estimator = Estimator(calc_mode=EstimatorCalcMode.CRYSTAL, model_version="v2.0.0")
calculator = ASECalculator(estimator)
[2]:
react = ase.Atoms(
symbols="C2NON2H3",
positions = [
[ 2.12, -0.48, 0.00],
[ 0.76, 0.16, -0.00],
[-0.28, -0.79, -0.00],
[ 0.57, 1.35, 0.00],
[-1.42, -0.28, 0.00],
[-2.48, 0.07, 0.00],
[ 2.67, -0.14, -0.88],
[ 2.06, -1.57, -0.00],
[ 2.67, -0.14, 0.88],
],
)
prod = ase.Atoms(
symbols="C2NON2H3",
positions=[
[ 2.30, -0.87, 0.00],
[ 0.68, 0.92, -0.00],
[ 0.99, -0.25, 0.00],
[ 0.25, 2.03, 0.00],
[-2.34, 0.11, 0.00],
[-3.13, -0.67, -0.00],
[ 2.87, -0.57, -0.89],
[ 2.19, -1.96, 0.00],
[ 2.87, -0.57, 0.89],
],
)
The first step is to perform structural optimization. The way of structural optimization is the same as before.
To facilitate future work, we will move the structures before and after the reaction so that they are as close as possible. This is done by minimize_rotation_and_translation method.
[3]:
react.calc = calculator
opt = BFGS(react)
opt.run(fmax=0.01)
prod.calc = calculator
opt = BFGS(prod)
opt.run(fmax=0.01)
minimize_rotation_and_translation(react, prod)
Step Time Energy fmax
BFGS: 0 07:11:23 -41.288279 3.3720
BFGS: 1 07:11:23 -41.050347 9.2272
BFGS: 2 07:11:24 -41.352632 0.9379
BFGS: 3 07:11:24 -41.360422 0.8468
BFGS: 4 07:11:24 -41.367067 0.4149
BFGS: 5 07:11:24 -41.370915 0.2225
BFGS: 6 07:11:25 -41.373628 0.1877
BFGS: 7 07:11:25 -41.374702 0.1327
BFGS: 8 07:11:25 -41.375998 0.1432
BFGS: 9 07:11:25 -41.377709 0.1950
BFGS: 10 07:11:26 -41.379433 0.1925
BFGS: 11 07:11:26 -41.380617 0.1081
BFGS: 12 07:11:26 -41.381467 0.0848
BFGS: 13 07:11:27 -41.382043 0.0992
BFGS: 14 07:11:27 -41.382362 0.0638
BFGS: 15 07:11:27 -41.382599 0.0515
BFGS: 16 07:11:27 -41.382840 0.0601
BFGS: 17 07:11:28 -41.383157 0.0738
BFGS: 18 07:11:28 -41.383481 0.0658
BFGS: 19 07:11:28 -41.383667 0.0388
BFGS: 20 07:11:28 -41.383759 0.0254
BFGS: 21 07:11:29 -41.383802 0.0258
BFGS: 22 07:11:29 -41.383865 0.0298
BFGS: 23 07:11:29 -41.383911 0.0271
BFGS: 24 07:11:29 -41.383936 0.0141
BFGS: 25 07:11:30 -41.383955 0.0091
Step Time Energy fmax
BFGS: 0 07:11:30 -42.925378 0.7301
BFGS: 1 07:11:30 -42.911607 1.6637
BFGS: 2 07:11:30 -42.928921 0.2516
BFGS: 3 07:11:31 -42.928744 0.5499
BFGS: 4 07:11:31 -42.931010 0.2054
BFGS: 5 07:11:31 -42.933623 0.2522
BFGS: 6 07:11:32 -42.937579 0.4179
BFGS: 7 07:11:32 -42.942069 0.4193
BFGS: 8 07:11:32 -42.946957 0.2883
BFGS: 9 07:11:32 -42.950890 0.1795
BFGS: 10 07:11:33 -42.954245 0.2710
BFGS: 11 07:11:33 -42.956652 0.2315
BFGS: 12 07:11:33 -42.958064 0.1185
BFGS: 13 07:11:34 -42.959282 0.1460
BFGS: 14 07:11:34 -42.960218 0.1848
BFGS: 15 07:11:34 -42.960760 0.1164
BFGS: 16 07:11:34 -42.961065 0.0687
BFGS: 17 07:11:34 -42.961415 0.1194
BFGS: 18 07:11:35 -42.962055 0.2019
BFGS: 19 07:11:35 -42.963172 0.2617
BFGS: 20 07:11:35 -42.964888 0.2410
BFGS: 21 07:11:36 -42.966962 0.1318
BFGS: 22 07:11:36 -42.969055 0.1111
BFGS: 23 07:11:36 -42.970605 0.1539
BFGS: 24 07:11:36 -42.971380 0.1119
BFGS: 25 07:11:36 -42.971683 0.0615
BFGS: 26 07:11:37 -42.971868 0.0360
BFGS: 27 07:11:37 -42.971954 0.0460
BFGS: 28 07:11:37 -42.972000 0.0245
BFGS: 29 07:11:37 -42.972008 0.0135
BFGS: 30 07:11:38 -42.972006 0.0211
BFGS: 31 07:11:38 -42.972034 0.0406
BFGS: 32 07:11:38 -42.972068 0.0647
BFGS: 33 07:11:38 -42.972157 0.0917
BFGS: 34 07:11:39 -42.972297 0.0956
BFGS: 35 07:11:39 -42.972386 0.0525
BFGS: 36 07:11:39 -42.972556 0.0783
BFGS: 37 07:11:39 -42.972755 0.1364
BFGS: 38 07:11:40 -42.972918 0.1228
BFGS: 39 07:11:40 -42.972963 0.1925
BFGS: 40 07:11:40 -42.973029 0.0464
BFGS: 41 07:11:41 -42.973089 0.0477
BFGS: 42 07:11:41 -42.973177 0.0895
BFGS: 43 07:11:41 -42.973234 0.0430
BFGS: 44 07:11:41 -42.973260 0.0586
BFGS: 45 07:11:42 -42.973297 0.0315
BFGS: 46 07:11:42 -42.973299 0.0083
We check the reactant & product structures by visualization.
[4]:
from ase.io import write
from IPython.display import Image
import matplotlib.pyplot as plt
import matplotlib.image as mpimg
write("output/curtius_react.png", react, rotation="0x,0y,0z")
write("output/curtius_prod.png", prod, rotation="0x,0y,0z")
fig, axes = plt.subplots(1, 2, figsize=(6, 3))
ax0, ax1 = axes
ax0.imshow(mpimg.imread("output/curtius_react.png"))
ax0.set_axis_off()
ax0.set_title("Reactant")
ax1.imshow(mpimg.imread("output/curtius_prod.png"))
ax1.set_axis_off()
ax1.set_title("Product")
fig.show()
[5]:
from pfcc_extras.visualize.view import view_ngl
view_ngl([react, prod], representations=["ball+stick"], replace_structure=True)
We will start to run NEB calculations from here. NEB is a method for finding reaction paths by interpolating many discrete structures between the structures before and after the reaction and then performing structural optimization on the combination of these structures in phase space.
We can use the functions in ASE for both interpolating the intermediate structure and performing NEB calculation. Note that in order to run the NEB calculation in parallel, we create a unique calculator instance for each structure, and specify allow_shared_calculator=False, parallel=True for NEB.
Because NEB contains multiple structures internally, the time cost of one step is greater than that of MD. Nevertheless, such a calculation can be completed in a few minutes with PFP, compared to several days or months with quantum chemical computation, e.g. DFT.
The calculations are performed in the following steps
First, create a list of ASE atoms as
images.The first one is the reactant
react, and the last one is the productprod.The coordinates of intermediate
imagewill be changed later with theneb.interpolate()method. Here, we only need to make a copy here.In this case, we have created 7 intermediate images, for a total of 9 images.
Create and set a
calculatorfor each atoms.Create NEB class
kis the spring constant that represents the strength of the spring connecting eachimages.By setting
climb=True, we use the Climbing Image NEB method, which is a method that reverts the energy gradient to find transition states.
neb.interpolate()Linear interpolation of the
imagesis executed to create a list of coordinates that gradually change from reactant to product.The initial configuration (candidates) of the reaction paths can be confirmed by commenting out the line with
view_ngl.
optimization with
FIREThe obtained reaction pathway candidates are modified to the appropriate reaction pathway by optimizing them using the
FIREmethod.
[6]:
images = [react.copy()]
images += [react.copy() for i in range (7)]
images += [prod.copy()]
for image in images:
estimator = Estimator(calc_mode=EstimatorCalcMode.CRYSTAL, model_version="v2.0.0")
calculator = ASECalculator(estimator)
image.calc = calculator
neb = NEB(images, k=0.1, climb=True, allow_shared_calculator=False, parallel=True)
neb.interpolate()
# Check interpolated images
# view_ngl(images, representations=["ball+stick"], replace_structure=True)
opt = FIRE(neb)
status = opt.run(fmax=0.05, steps=500)
Step Time Energy fmax
*Force-consistent energies used in optimization.
FIRE: 0 07:13:17 -38.199767* 23.4912
FIRE: 1 07:13:17 -39.254614* 6.5052
FIRE: 2 07:13:18 -39.008607* 11.1636
FIRE: 3 07:13:18 -39.405843* 7.0572
FIRE: 4 07:13:18 -39.450036* 6.9048
FIRE: 5 07:13:18 -39.462446* 5.3536
FIRE: 6 07:13:19 -39.478978* 3.7188
FIRE: 7 07:13:19 -39.489785* 3.0161
FIRE: 8 07:13:19 -39.491196* 2.8974
FIRE: 9 07:13:20 -39.488352* 3.9156
FIRE: 10 07:13:20 -39.488611* 4.0391
FIRE: 11 07:13:20 -39.495302* 3.5124
FIRE: 12 07:13:21 -39.507837* 2.8101
FIRE: 13 07:13:21 -39.519245* 1.8177
FIRE: 14 07:13:21 -39.525173* 3.2658
FIRE: 15 07:13:21 -39.539541* 2.8926
FIRE: 16 07:13:22 -39.572465* 1.9634
FIRE: 17 07:13:22 -39.604970* 2.3038
FIRE: 18 07:13:22 -39.606554* 2.0281
FIRE: 19 07:13:23 -39.609300* 1.5137
FIRE: 20 07:13:23 -39.612558* 0.9486
FIRE: 21 07:13:23 -39.615882* 0.9099
FIRE: 22 07:13:24 -39.619306* 1.1980
FIRE: 23 07:13:24 -39.623187* 1.5334
FIRE: 24 07:13:24 -39.627865* 1.5309
FIRE: 25 07:13:24 -39.633946* 1.1657
FIRE: 26 07:13:25 -39.641256* 0.7378
FIRE: 27 07:13:25 -39.649407* 0.8005
FIRE: 28 07:13:25 -39.658478* 1.2002
FIRE: 29 07:13:26 -39.669462* 1.0858
FIRE: 30 07:13:26 -39.682556* 0.6026
FIRE: 31 07:13:26 -39.696465* 0.9078
FIRE: 32 07:13:27 -39.710635* 1.0347
FIRE: 33 07:13:27 -39.725172* 0.6204
FIRE: 34 07:13:27 -39.737158* 1.1204
FIRE: 35 07:13:27 -39.747343* 0.5140
FIRE: 36 07:13:28 -39.752792* 1.0859
FIRE: 37 07:13:28 -39.755747* 0.8419
FIRE: 38 07:13:28 -39.753584* 0.5378
FIRE: 39 07:13:29 -39.748985* 1.0645
FIRE: 40 07:13:29 -39.752538* 0.4342
FIRE: 41 07:13:29 -39.753741* 0.6426
FIRE: 42 07:13:29 -39.754036* 0.5796
FIRE: 43 07:13:30 -39.754571* 0.4683
FIRE: 44 07:13:30 -39.755176* 0.3427
FIRE: 45 07:13:30 -39.755772* 0.3333
FIRE: 46 07:13:31 -39.756359* 0.3786
FIRE: 47 07:13:31 -39.757009* 0.4330
FIRE: 48 07:13:31 -39.757876* 0.4118
FIRE: 49 07:13:32 -39.758968* 0.3195
FIRE: 50 07:13:32 -39.760281* 0.2992
FIRE: 51 07:13:32 -39.761587* 0.3999
FIRE: 52 07:13:33 -39.762970* 0.4432
FIRE: 53 07:13:33 -39.764575* 0.3677
FIRE: 54 07:13:34 -39.766266* 0.2657
FIRE: 55 07:13:34 -39.767656* 0.3065
FIRE: 56 07:13:35 -39.768928* 0.2712
FIRE: 57 07:13:35 -39.769798* 0.3823
FIRE: 58 07:13:35 -39.770114* 0.2827
FIRE: 59 07:13:36 -39.769901* 0.3204
FIRE: 60 07:13:36 -39.769478* 0.2699
FIRE: 61 07:13:37 -39.768488* 0.3123
FIRE: 62 07:13:37 -39.768018* 0.3649
FIRE: 63 07:13:37 -39.768416* 0.2041
FIRE: 64 07:13:38 -39.768642* 0.2207
FIRE: 65 07:13:38 -39.768710* 0.3374
FIRE: 66 07:13:39 -39.769102* 0.2798
FIRE: 67 07:13:39 -39.769614* 0.1618
FIRE: 68 07:13:39 -39.769868* 0.2333
FIRE: 69 07:13:40 -39.770150* 0.2230
FIRE: 70 07:13:40 -39.770635* 0.1537
FIRE: 71 07:13:41 -39.770827* 0.2656
FIRE: 72 07:13:41 -39.771113* 0.1617
FIRE: 73 07:13:42 -39.771178* 0.2339
FIRE: 74 07:13:42 -39.771294* 0.1410
FIRE: 75 07:13:42 -39.770918* 0.3239
FIRE: 76 07:13:43 -39.771114* 0.2307
FIRE: 77 07:13:43 -39.771130* 0.2419
FIRE: 78 07:13:44 -39.770793* 0.3642
FIRE: 79 07:13:44 -39.771259* 0.1196
FIRE: 80 07:13:44 -39.771153* 0.2493
FIRE: 81 07:13:45 -39.771224* 0.2040
FIRE: 82 07:13:45 -39.771311* 0.1449
FIRE: 83 07:13:46 -39.771377* 0.1116
FIRE: 84 07:13:46 -39.771394* 0.1337
FIRE: 85 07:13:46 -39.771365* 0.1640
FIRE: 86 07:13:47 -39.771384* 0.1680
FIRE: 87 07:13:47 -39.771443* 0.1443
FIRE: 88 07:13:48 -39.771532* 0.1119
FIRE: 89 07:13:48 -39.771601* 0.1236
FIRE: 90 07:13:49 -39.771593* 0.1528
FIRE: 91 07:13:49 -39.771657* 0.1381
FIRE: 92 07:13:49 -39.771759* 0.1046
FIRE: 93 07:13:50 -39.771788* 0.1323
FIRE: 94 07:13:50 -39.771831* 0.1146
FIRE: 95 07:13:51 -39.771917* 0.1162
FIRE: 96 07:13:51 -39.771895* 0.1199
FIRE: 97 07:13:51 -39.771969* 0.1137
FIRE: 98 07:13:52 -39.771947* 0.1062
FIRE: 99 07:13:52 -39.771940* 0.1201
FIRE: 100 07:13:53 -39.772026* 0.1290
FIRE: 101 07:13:53 -39.772065* 0.1533
FIRE: 102 07:13:53 -39.772045* 0.2577
FIRE: 103 07:13:54 -39.772143* 0.0641
FIRE: 104 07:13:54 -39.772100* 0.2081
FIRE: 105 07:13:55 -39.772121* 0.1670
FIRE: 106 07:13:55 -39.772161* 0.0945
FIRE: 107 07:13:55 -39.772160* 0.0595
FIRE: 108 07:13:56 -39.772143* 0.1189
FIRE: 109 07:13:56 -39.772150* 0.1140
FIRE: 110 07:13:57 -39.772159* 0.1047
FIRE: 111 07:13:57 -39.772169* 0.0916
FIRE: 112 07:13:57 -39.772171* 0.0762
FIRE: 113 07:13:58 -39.772176* 0.0630
FIRE: 114 07:13:58 -39.772174* 0.0556
FIRE: 115 07:13:59 -39.772165* 0.0558
FIRE: 116 07:13:59 -39.772176* 0.0618
FIRE: 117 07:13:59 -39.772177* 0.0688
FIRE: 118 07:14:00 -39.772181* 0.0723
FIRE: 119 07:14:00 -39.772179* 0.0695
FIRE: 120 07:14:01 -39.772188* 0.0611
FIRE: 121 07:14:01 -39.772200* 0.0539
FIRE: 122 07:14:01 -39.772202* 0.0583
FIRE: 123 07:14:02 -39.772204* 0.0673
FIRE: 124 07:14:02 -39.772207* 0.0658
FIRE: 125 07:14:03 -39.772225* 0.0539
FIRE: 126 07:14:03 -39.772234* 0.0547
FIRE: 127 07:14:04 -39.772249* 0.0589
FIRE: 128 07:14:04 -39.772255* 0.0490
If you want to suppress default logging, please specify the logfile parameter like FIRE(neb, logfile=None).
Let’s look at the results. The Notebook environment is also suitable for the graphical visualization of numerical data. Below is a visualization of an energy profile using matplotlib.
You can see the energy going up once and then going down in the reaction. The energy difference between the point of maximum energy (transition state) and each of the left and right stable points corresponds to the activation energy.
[7]:
energies = np.array([image.get_total_energy() for image in images])
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.plot(energies)
ax.set_xlabel("replica")
ax.set_ylabel("energy")
fig.show()
The activation energy can be calculated as follows.
Here, the activation energy in the forward reaction X→Y is denoted as E_act_forward, and reverse reaction Y→X is denoted as E_act_backward.
[8]:
# Transition state takes maximum energy in the reaction path
ts_index = np.argmax(energies)
E_act_forward = energies[ts_index] - energies[0]
E_act_backward = energies[ts_index] - energies[-1]
print(f"ts_index = {ts_index}")
print(f"E_act_forward = {E_act_forward:.2f} eV")
print(f"E_act_backward = {E_act_backward:.2f} eV")
ts_index = 2
E_act_forward = 1.61 eV
E_act_backward = 3.20 eV
You can visualize the structure in the Notebook as well, or you can save the image as a file. Let’s have a try.
The results are output as a trajectories rather than a single structure since it is NEB calculation. Let’s save it as an animated GIF using ASE.
[9]:
fig = plt.figure(facecolor="white")
ax = fig.add_subplot()
ase.io.write(
"output/curtius_NEB.gif",
images,
format="gif",
ax=ax
)
You can see a static image of the final structure above, but if you double-click on “curtius_NEB.gif” in the file viewer on the left, you can directly see the animated GIF.
If the calculations are successful, you will see the following animated GIFs.
[10]:
Image("output/curtius_NEB.gif", format="gif")
[10]:
<IPython.core.display.Image object>
Finally, let’s compare how the NEB method changed the reaction path before and after optimization.
Call the neb.interpolate() function and compare the energy of the candidate structure interpolated_images of the reaction pathway before FIRE opt, and images after optimization by the NEB.
[11]:
_images = [react.copy()]
_images += [react.copy() for i in range (7)]
_images += [prod.copy()]
for image in _images:
estimator = Estimator()
calculator = ASECalculator(estimator)
image.calc = calculator
neb = NEB(_images, k=0.1, climb=True, allow_shared_calculator=False, parallel=True)
neb.interpolate()
interpolated_images = _images
[12]:
initial_energies = [image.get_total_energy() for image in interpolated_images]
opt_energies = [image.get_total_energy() for image in images]
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.plot(initial_energies, label="neb interpolated", color="red")
ax.plot(energies, label="after neb opt")
ax.set_xlabel("replica")
ax.set_ylabel("energy")
ax.legend()
fig.show()
As can be seen from the above figure, NEB method has led to the discovery of a pathway with lower activation energy than the candidate reaction pathway (red line), i.e., a pathway with lower transition state energy (blue line).
[Column] Comparison of reaction pathway analysis methods and molecular dynamics¶
In this section, the reaction pathway was explored using the NEB method, which specifies the starting and ending states.
Molecular dynamics (MD) methods, described in the next chapter, can follow the time evolution of atoms. So, what about running the MD simulation until the chemical reaction occurs naturally? There are examples where activation energies can be analyzed by MD as well. The diffusion phenomenon of the Li ion is such an example.
However, reaction time is a challenging issue in general. Reactions are rare events. As you can see from the Arrhenius equation, it is a phenomenon whose frequency of occurrence changes exponentially with respect to the activation energy \(E_a\).
Therefore, MD can see the reaction if the activation energy is very low. However, if the activation energy becomes higher, the reaction does not occur in the time scale that MD can handle (~ns) and cannot be detected efficiently. Futhermore, if the activation energy is even higher, it will be a phenomenon that does not occur even in the real world.
For reference, let us look at the relationship between activation energy and the likelihood of a phenomenon. At 300 K, which is about room temperature, the value of \(k_B T\), which is the denominator inside the exponential function of the Arrhenius equation, is roughly 0.026 eV. This means that the occurrence probablilty of a phenomenon with 1 eV higher activation energy is reduced by a factor of \(\exp\left(-1/0.026\right)\), or \(2\times10^{-17}\). Even considering the difference between macroscale and atomic scale, this is still a large differencet and the reaction can be regarded as a rare event. Similarly, \(k_B T\) is 0.086 eV and 0.172 eV at 1000 K and 2000 K, respectively. It means that the likelihood of a phenomenon with 1 eV higher activation energy is about \(9 \times10^{-6}\) and \(0.003\) times lower, respectively. This is significantly larger than the value at 300 K, and can be regarded as a frequently occurring phenomenon on the macroscopic scale. In general, the occurrence probability of the phenomenon is very sensitive to the activation energy and temperature. By using reaction path analysis such as the NEB method, it is possible to efficiently analyze such rare events across exponential scales.
[13]:
from ase.units import kB
from math import exp
for T in [300, 1000, 2000]:
kBT = kB * T
ratio = exp(-1 / (kBT))
print(f"----- T = {T} K -----")
print(f"kB T : {kBT:.3f}")
print(f"ratio: {ratio:.2e}")
----- T = 300 K -----
kB T : 0.026
ratio: 1.59e-17
----- T = 1000 K -----
kB T : 0.086
ratio: 9.12e-06
----- T = 2000 K -----
kB T : 0.172
ratio: 3.02e-03
In this sense, it may be possible to efficiently explore reactions that are difficult to see with MD by using the NEB method.
In addition, several approaches for reaction path searching are described below.
Approach | Characteristic | Methods |
|---|---|---|
Reaction pathway search for a specific reaction pathway | A method to find the minimum energy path (MEP) in a neighborhood starting from a path connecting two stable structures. Although the MEP and the corresponding transition state can be found with relatively low computational cost, the results strongly depend on the reaction path assumed before the search because it is local search method. | NEB, String method etc. |
Reaction path search without assuming reaction paths in advance | A method to find one or more reaction pathways without a priori assumptions about reaction pathways or post-reaction structures. In addition to the standard method of MD, there are other methods such as metadynamics, which aims to efficiently search for energy surfaces, and a method that uses information on local energy surfaces to search (ADDF). | MD, metadynamics, ADDF etc. |
Reference¶
The following references may be useful to those who would like to study further.