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 14:43:15 -41.288286 3.372026
BFGS: 1 14:43:15 -41.050344 9.227204
BFGS: 2 14:43:15 -41.352642 0.937929
BFGS: 3 14:43:15 -41.360438 0.846813
BFGS: 4 14:43:15 -41.367057 0.414904
BFGS: 5 14:43:16 -41.370927 0.222537
BFGS: 6 14:43:16 -41.373629 0.187670
BFGS: 7 14:43:16 -41.374702 0.132663
BFGS: 8 14:43:16 -41.376004 0.143219
BFGS: 9 14:43:16 -41.377716 0.194960
BFGS: 10 14:43:16 -41.379438 0.192537
BFGS: 11 14:43:16 -41.380632 0.108121
BFGS: 12 14:43:16 -41.381465 0.084829
BFGS: 13 14:43:16 -41.382037 0.099253
BFGS: 14 14:43:16 -41.382358 0.063755
BFGS: 15 14:43:16 -41.382590 0.051461
BFGS: 16 14:43:16 -41.382839 0.060071
BFGS: 17 14:43:16 -41.383148 0.073773
BFGS: 18 14:43:16 -41.383478 0.065854
BFGS: 19 14:43:16 -41.383675 0.038820
BFGS: 20 14:43:16 -41.383762 0.025375
BFGS: 21 14:43:16 -41.383797 0.025825
BFGS: 22 14:43:16 -41.383857 0.029769
BFGS: 23 14:43:17 -41.383902 0.027067
BFGS: 24 14:43:17 -41.383944 0.014089
BFGS: 25 14:43:17 -41.383969 0.009115
Step Time Energy fmax
BFGS: 0 14:43:17 -42.925393 0.730117
BFGS: 1 14:43:17 -42.911635 1.663660
BFGS: 2 14:43:17 -42.928937 0.251647
BFGS: 3 14:43:17 -42.928764 0.549965
BFGS: 4 14:43:17 -42.931020 0.205362
BFGS: 5 14:43:17 -42.933635 0.252208
BFGS: 6 14:43:17 -42.937590 0.417864
BFGS: 7 14:43:17 -42.942083 0.419333
BFGS: 8 14:43:17 -42.946968 0.288300
BFGS: 9 14:43:17 -42.950906 0.179478
BFGS: 10 14:43:17 -42.954262 0.270988
BFGS: 11 14:43:17 -42.956673 0.231530
BFGS: 12 14:43:17 -42.958089 0.118511
BFGS: 13 14:43:17 -42.959289 0.146021
BFGS: 14 14:43:17 -42.960235 0.184809
BFGS: 15 14:43:17 -42.960773 0.116383
BFGS: 16 14:43:18 -42.961079 0.068664
BFGS: 17 14:43:18 -42.961431 0.119403
BFGS: 18 14:43:18 -42.962067 0.201879
BFGS: 19 14:43:18 -42.963184 0.261755
BFGS: 20 14:43:18 -42.964901 0.241058
BFGS: 21 14:43:18 -42.966974 0.131766
BFGS: 22 14:43:18 -42.969068 0.111116
BFGS: 23 14:43:18 -42.970628 0.153895
BFGS: 24 14:43:18 -42.971380 0.111868
BFGS: 25 14:43:18 -42.971700 0.061562
BFGS: 26 14:43:18 -42.971887 0.035948
BFGS: 27 14:43:18 -42.971988 0.045980
BFGS: 28 14:43:18 -42.972005 0.024488
BFGS: 29 14:43:18 -42.972040 0.013546
BFGS: 30 14:43:18 -42.972020 0.021100
BFGS: 31 14:43:18 -42.972060 0.040577
BFGS: 32 14:43:18 -42.972098 0.064734
BFGS: 33 14:43:18 -42.972179 0.091686
BFGS: 34 14:43:19 -42.972329 0.095636
BFGS: 35 14:43:19 -42.972412 0.052399
BFGS: 36 14:43:19 -42.972583 0.078331
BFGS: 37 14:43:19 -42.972768 0.136407
BFGS: 38 14:43:19 -42.972924 0.122819
BFGS: 39 14:43:19 -42.972981 0.193897
BFGS: 40 14:43:19 -42.973057 0.047662
BFGS: 41 14:43:19 -42.973107 0.047561
BFGS: 42 14:43:19 -42.973201 0.088515
BFGS: 43 14:43:19 -42.973237 0.042605
BFGS: 44 14:43:19 -42.973302 0.058457
BFGS: 45 14:43:19 -42.973299 0.032923
BFGS: 46 14:43:19 -42.973306 0.007837
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
FIRE: 0 14:43:27 -38.199809 23.490956
FIRE: 1 14:43:27 -39.254621 6.505189
/tmp/ipykernel_25517/2722927025.py:8: FutureWarning: Please import NEB from ase.mep, not ase.neb.
neb = NEB(images, k=0.1, climb=True, allow_shared_calculator=False, parallel=True)
FIRE: 2 14:43:27 -39.008617 11.163567
FIRE: 3 14:43:27 -39.405844 7.057152
FIRE: 4 14:43:27 -39.450020 6.904893
FIRE: 5 14:43:27 -39.462435 5.353684
FIRE: 6 14:43:27 -39.478987 3.718861
FIRE: 7 14:43:27 -39.489786 3.016063
FIRE: 8 14:43:27 -39.491198 2.897458
FIRE: 9 14:43:27 -39.488352 3.915679
FIRE: 10 14:43:27 -39.488606 4.039161
FIRE: 11 14:43:28 -39.495308 3.512413
FIRE: 12 14:43:28 -39.507829 2.810146
FIRE: 13 14:43:28 -39.519259 1.817654
FIRE: 14 14:43:28 -39.525174 3.265853
FIRE: 15 14:43:28 -39.539539 2.892644
FIRE: 16 14:43:28 -39.572475 1.963444
FIRE: 17 14:43:28 -39.604964 2.303778
FIRE: 18 14:43:28 -39.606562 2.028061
FIRE: 19 14:43:28 -39.609304 1.513654
FIRE: 20 14:43:28 -39.612562 0.948627
FIRE: 21 14:43:28 -39.615873 0.909884
FIRE: 22 14:43:28 -39.619310 1.198012
FIRE: 23 14:43:29 -39.623196 1.533371
FIRE: 24 14:43:29 -39.627874 1.530900
FIRE: 25 14:43:29 -39.633955 1.165640
FIRE: 26 14:43:29 -39.641255 0.737797
FIRE: 27 14:43:29 -39.649404 0.800478
FIRE: 28 14:43:29 -39.658483 1.200173
FIRE: 29 14:43:29 -39.669464 1.085851
FIRE: 30 14:43:29 -39.682563 0.602644
FIRE: 31 14:43:29 -39.696479 0.907790
FIRE: 32 14:43:29 -39.710647 1.034686
FIRE: 33 14:43:29 -39.725190 0.620371
FIRE: 34 14:43:29 -39.737168 1.120402
FIRE: 35 14:43:29 -39.747340 0.514008
FIRE: 36 14:43:30 -39.752792 1.085785
FIRE: 37 14:43:30 -39.755755 0.841935
FIRE: 38 14:43:30 -39.753593 0.537654
FIRE: 39 14:43:30 -39.748983 1.064867
FIRE: 40 14:43:30 -39.752546 0.434176
FIRE: 41 14:43:30 -39.753732 0.642694
FIRE: 42 14:43:30 -39.754038 0.579657
FIRE: 43 14:43:30 -39.754581 0.468322
FIRE: 44 14:43:30 -39.755184 0.342757
FIRE: 45 14:43:30 -39.755784 0.333293
FIRE: 46 14:43:30 -39.756348 0.378653
FIRE: 47 14:43:30 -39.757035 0.433039
FIRE: 48 14:43:31 -39.757874 0.411865
FIRE: 49 14:43:31 -39.758987 0.319442
FIRE: 50 14:43:31 -39.760280 0.299208
FIRE: 51 14:43:31 -39.761602 0.399914
FIRE: 52 14:43:31 -39.762984 0.443207
FIRE: 53 14:43:31 -39.764582 0.367706
FIRE: 54 14:43:31 -39.766270 0.265688
FIRE: 55 14:43:31 -39.767651 0.306568
FIRE: 56 14:43:31 -39.768935 0.271192
FIRE: 57 14:43:31 -39.769791 0.382343
FIRE: 58 14:43:31 -39.770105 0.282708
FIRE: 59 14:43:31 -39.769909 0.320444
FIRE: 60 14:43:32 -39.769471 0.269958
FIRE: 61 14:43:32 -39.768488 0.312202
FIRE: 62 14:43:32 -39.768029 0.365170
FIRE: 63 14:43:32 -39.768429 0.204112
FIRE: 64 14:43:32 -39.768648 0.220764
FIRE: 65 14:43:32 -39.768724 0.337588
FIRE: 66 14:43:32 -39.769099 0.279886
FIRE: 67 14:43:32 -39.769610 0.161808
FIRE: 68 14:43:32 -39.769876 0.233502
FIRE: 69 14:43:32 -39.770145 0.223181
FIRE: 70 14:43:32 -39.770639 0.153745
FIRE: 71 14:43:33 -39.770829 0.265734
FIRE: 72 14:43:33 -39.771117 0.161783
FIRE: 73 14:43:33 -39.771184 0.234144
FIRE: 74 14:43:33 -39.771296 0.141023
FIRE: 75 14:43:33 -39.770925 0.323859
FIRE: 76 14:43:33 -39.771124 0.230935
FIRE: 77 14:43:33 -39.771124 0.242028
FIRE: 78 14:43:33 -39.770795 0.364439
FIRE: 79 14:43:33 -39.771262 0.119569
FIRE: 80 14:43:33 -39.771155 0.249417
FIRE: 81 14:43:33 -39.771227 0.204118
FIRE: 82 14:43:33 -39.771324 0.144965
FIRE: 83 14:43:33 -39.771389 0.111566
FIRE: 84 14:43:34 -39.771388 0.133739
FIRE: 85 14:43:34 -39.771367 0.164100
FIRE: 86 14:43:34 -39.771395 0.168005
FIRE: 87 14:43:34 -39.771444 0.144330
FIRE: 88 14:43:34 -39.771531 0.111891
FIRE: 89 14:43:34 -39.771593 0.123567
FIRE: 90 14:43:34 -39.771598 0.152820
FIRE: 91 14:43:34 -39.771657 0.138157
FIRE: 92 14:43:34 -39.771756 0.104620
FIRE: 93 14:43:34 -39.771798 0.132303
FIRE: 94 14:43:34 -39.771839 0.114626
FIRE: 95 14:43:34 -39.771912 0.116256
FIRE: 96 14:43:35 -39.771909 0.119960
FIRE: 97 14:43:35 -39.771974 0.113715
FIRE: 98 14:43:35 -39.771948 0.106299
FIRE: 99 14:43:35 -39.771956 0.120135
FIRE: 100 14:43:35 -39.772026 0.129044
FIRE: 101 14:43:35 -39.772069 0.153474
FIRE: 102 14:43:35 -39.772048 0.258043
FIRE: 103 14:43:35 -39.772141 0.064037
FIRE: 104 14:43:35 -39.772104 0.207688
FIRE: 105 14:43:35 -39.772121 0.166595
FIRE: 106 14:43:35 -39.772154 0.094351
FIRE: 107 14:43:35 -39.772153 0.059529
FIRE: 108 14:43:36 -39.772165 0.118646
FIRE: 109 14:43:36 -39.772160 0.113806
FIRE: 110 14:43:36 -39.772162 0.104494
FIRE: 111 14:43:36 -39.772173 0.091514
FIRE: 112 14:43:36 -39.772172 0.076141
FIRE: 113 14:43:36 -39.772176 0.063003
FIRE: 114 14:43:36 -39.772170 0.055623
FIRE: 115 14:43:36 -39.772164 0.055830
FIRE: 116 14:43:36 -39.772183 0.061785
FIRE: 117 14:43:36 -39.772179 0.068822
FIRE: 118 14:43:36 -39.772177 0.072309
FIRE: 119 14:43:36 -39.772180 0.069552
FIRE: 120 14:43:37 -39.772191 0.061102
FIRE: 121 14:43:37 -39.772186 0.053883
FIRE: 122 14:43:37 -39.772210 0.058257
FIRE: 123 14:43:37 -39.772206 0.067317
FIRE: 124 14:43:37 -39.772221 0.065830
FIRE: 125 14:43:37 -39.772235 0.053900
FIRE: 126 14:43:37 -39.772238 0.054710
FIRE: 127 14:43:37 -39.772259 0.058955
FIRE: 128 14:43:37 -39.772262 0.049006
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
/tmp/ipykernel_25517/4187733795.py:8: FutureWarning: Please import NEB from ase.mep, not ase.neb.
neb = NEB(_images, k=0.1, climb=True, allow_shared_calculator=False, parallel=True)
[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.