Part 2: Static approaches ¶
Note: The machine-learned force fields (MLFF) used in this tutorial can be downloaded from here. The mlff directory should be placed in the transition_states directory.
By the end of this tutorial, you will be able to:
- perform a vibrational analysis on a trial transition state (TS)
- find a trial unstable direction for TS relaxation
5.1 Task¶
Perform vibrational analysis on a trial transition state for the carbonylation of methylated acidic chabazite during ethanol synthesis from syngas
In part 1, the nudged elastic band (NEB), improved dimer method (IDM), and intrinsic reaction coordinate (IRC) have been covered for separate reactions. A real system, such as a zeolite-catalyzed reaction [Nat. Rev. Mater. 6, 1156 (2021)], is studied by using a combination of these methods. A model reaction is that of the conversion of syngas (CO and H2) to ethanol (CH3CH2OH), catalyzed by chabazite. Chabazite is a type of zeolite Z:
The conversion of syngas to ethanol CH3CH2OH proceeds via six elementary steps:
- CO + 2H2 → CH3OH
- CH3OH + H-Z → H2O + CH3-Z
- CO + CH3-Z → CH3CO+ + Z-
- CH3CO+ + Z- → CH3CO-Z
- CH3OH + CH3CO-Z → CH3COOCH3 + H-Z
- CH3COOCH3 + 2H2 → CH3CH2OH + CH3OH
The rate-determining step (RDS) for this reaction is step 3, the carbonylation of the methylated chabazite CH3-Z [PCCP 13, 2603 (2011)]. During the RDS, a free acetyl cation CH3-C=O+ is created as a product. Depending on the thermodynamic conditions, the acetyl cation attaches to the zeolite framework and forms CH3-C=O-Z, or transforms into a free CH2-C=O+ cation via protonation of the zeolite [J. Catalysis 396, 166 (2021)]. We shall focus on calculating the free energy of activation $\Delta A_{R\rightarrow P}^{\ddagger}$ for the forward RDS, meaning that the state of the product after the carbonylation is not relevant. Note that $\ddagger$ denotes the TS.
In the first exercise, an initial guess for the transition state of the carbonylation step has been made based on the geometry of the reactant and the expected product. The planar methyl CH3 group is assumed to lie halfway between the O on the zeolite and the C in the CO molecule. Both distances have been arbitrarily chosen and set at 2.1 Å. If the initial guess is close to a saddle point, then it will have only one imaginary frequency. If there are multiple imaginary frequencies, then the strongest imaginary frequency can be chosen as a trial direction for IDM.
5.2 Input¶
The input files to run this example are prepared at $TUTORIALS/transition_states/e05_Static_Vib_analysis.
Click here to see the POSCAR!
Trial TS structure
1.0
7.948754845851 -0.000049967724 4.901433149086
-3.974564067268 6.883894310696 4.901254785119
-3.974438768620 -6.883777142027 4.901245706606
Al C H O Si
1 2 3 25 11
Direct
0.62848612 0.39020196 0.33941230
0.36288164 0.94589396 0.59231366
0.42758700 0.08606810 0.44090475
0.38521995 0.16746240 0.50939566
0.35387699 0.01396078 0.36712590
0.53394788 0.05555254 0.46943474
0.23383670 0.75907096 0.99926471
0.17582313 0.47557306 0.95376242
0.25143428 0.62492408 0.73841260
0.51763691 0.75477069 0.73236766
0.02106794 0.34175363 0.14047271
0.35804237 0.65514737 0.49107130
0.24853205 0.51747545 0.24123557
0.40287652 0.76921155 0.25065278
0.46403989 0.47054569 0.66761133
0.75864353 0.75765792 0.91311194
0.32386156 0.85942886 0.68905833
0.58499859 0.24499254 0.74979789
0.51845405 0.83441880 0.01369959
0.27664375 0.26778763 0.08842135
0.74293017 0.50117935 0.76502281
0.53321230 0.53566626 0.30950861
0.79243426 0.25040725 0.96486076
0.66218192 0.35287152 0.51459701
0.48986783 0.22401662 0.28665802
0.65940134 0.00054233 0.84406385
0.99725337 0.65154960 0.83499460
0.33973023 0.00395889 0.13856936
0.51860538 0.17877751 0.00498267
0.76337083 0.36476411 0.23182904
0.85543558 0.53458809 0.03798978
0.16449753 0.62721478 0.88098785
0.39854648 0.62563842 0.65745894
0.61309740 0.83529930 0.87512616
0.37427231 0.83711848 0.09872232
0.38977842 0.61630578 0.32395651
0.18006949 0.40397070 0.10694487
0.61556880 0.39365097 0.67253106
0.83798928 0.61036800 0.88777481
0.64096697 0.16866303 0.89047296
0.40196983 0.16620429 0.11945109
0.85497734 0.37373749 0.09656525
! Vibrational analysis
IBRION = 5 ! Finite differences
NFREE = 2 ! Number of displacements
POTIM = 0.01 ! Width of displacement
NSW=1 ! Number of ionic steps
! Machine learning
ML_LMLFF = T ! Enables the use of MLFF
ML_MODE = run ! Prediction-only mode
K-Points
0
Gamma
1 1 1
0 0 0
Pseudopotentials of Al, C, H, O, and Si.
Check the tags set in the INCAR file!
The first four tags are for calculating the vibrational frequencies. IBRION = 5 sets up a finite differences calculation, with a displacement of POTIM and NFREE number of displacements. NSW usually sets the number of ionic steps; for IBRION = 5, it just has to be positive which is then reset by VASP to three times the number of ions in the cell.
The remaining two tags turn on the machine learning force field (MLFF). Notice that there are no additional tags for the electronic structure, as we are using MLFF. ML_LMLFF enables the use of MLFF, while ML_MODE = run sets the ML to run only, as opposed to additional training.
The KPOINTS are set to the Γ-point. POTCAR defines the pseudopotentials used. The POSCAR file defines the transition state inside the chabazite.
The settings used in this tutorial are sufficient to run quickly within the allotted time. To obtain higher-quality results, calculations may be run without an MLFF. However, this will increase the time significantly.
5.3 Calculation¶
Open a terminal, navigate to this example's directory, and create a symbolic link to the machine-learned force field (MLFF). Set up the directory to calculate the vibrational (phonon) frequencies of the trial transition state structure and run the calculation with the following:
cd $TUTORIALS/transition_states/e05_*
ln -s $TUTORIALS/transition_states/mlff/MLFF_e05-12/ML_FFN ML_FF
mpirun -np 4 vasp_gam
The calculation runs quickly, within a few seconds. If this were done ab initio, it would take several hours.
Let's take a look at the vibrational modes. vibFreq.py extracts the vibrational frequencies and prepares files to visualize them. Run it using the following:
cd $TUTORIALS/transition_states/e05_*
python3 ../scripts/vibFreq/vibFreq.py 0.0 0.2 0.0 > freq.dat
There are 126 different frequency modes calculated. Run the following to look at the first eight:
grep cm-1 $TUTORIALS/transition_states/e05_*/freq.dat | head -n 8
%%bash
grep cm-1 e05_*/freq.dat | head -n 8
Alternatively, you can obtain similar results with py4vasp
from py4vasp import Calculation
calculation = Calculation.from_path("e05_Static_Vib_analysis")
calculation.phonon.mode.print()
The f/i at the start of the first five modes indicates that these are imaginary frequencies, i.e. our structure is not a minimum. The structure has been chosen to approximate the expected transition state but still contains some imaginary frequency modes.
from py4vasp import Calculation
calculation = Calculation.from_path("e05_Static_Vib_analysis")
with open("e05_Static_Vib_analysis/transition_state.molden", "w") as file:
file.write(calculation.force_constant.to_molden())
In the next exercise, you will relax the transition state using the results from the frequency calculation, reducing the number of imaginary frequencies from five to one.
5.4 Questions¶
- Why are there imaginary frequencies?
- How many imaginary frequencies should there be in the transition state?
By the end of this tutorial, you will be able to:
- optimize a guess structure of the transition state
- compare the new transition state's geometry to the initial guess
6.1 Task¶
Relax the transition state structure along the hardest imaginary frequency for carbonylation of methylated chabazite using the improved dimer method
The transition state is optimized by following the direction along the hardest imaginary frequency mode using the improved dimer method (IDM) [J. Chem. Phys. 111 (1999) 7010]. IDM is a local saddle-point search algorithm using forces, as opposed to the more computationally demanding Hessian. It has been described in detail in Example 2, so will only be briefly reviewed here. From the initial trial structure, two points, one point forward and the other backward, are defined along a dimer axis N. The dimer is then rotated on the potential energy surface (PES) about the midpoint to find the direction with the lowest curvature λ on the PES. After rotation, the forces f along the dimer axis are calculated, and the midpoint of the dimer is translated in the reverse direction toward the saddle point. Each IDM step takes approximately four ionic steps and is repeated iteratively until convergence is reached and the transition state is obtained.
You have already encountered this approach in an earlier exercise. In this exericse, you will optimize the transition state for the rate-determining step, the carbonylation of a methyl group in chabazite, using IDM.
6.2 Input¶
The input files to run this example are prepared at $TUTORIALS/transition_states/e06_Static_TS_Opt.
POSCAR.init
Click here to reveal the POSCAR file!
Trial TS structure
1.0
7.948754845851 -0.000049967724 4.901433149086
-3.974564067268 6.883894310696 4.901254785119
-3.974438768620 -6.883777142027 4.901245706606
Al C H O Si
1 2 3 25 11
Direct
0.62848612 0.39020196 0.33941230
0.36288164 0.94589396 0.59231366
0.42758700 0.08606810 0.44090475
0.38521995 0.16746240 0.50939566
0.35387699 0.01396078 0.36712590
0.53394788 0.05555254 0.46943474
0.23383670 0.75907096 0.99926471
0.17582313 0.47557306 0.95376242
0.25143428 0.62492408 0.73841260
0.51763691 0.75477069 0.73236766
0.02106794 0.34175363 0.14047271
0.35804237 0.65514737 0.49107130
0.24853205 0.51747545 0.24123557
0.40287652 0.76921155 0.25065278
0.46403989 0.47054569 0.66761133
0.75864353 0.75765792 0.91311194
0.32386156 0.85942886 0.68905833
0.58499859 0.24499254 0.74979789
0.51845405 0.83441880 0.01369959
0.27664375 0.26778763 0.08842135
0.74293017 0.50117935 0.76502281
0.53321230 0.53566626 0.30950861
0.79243426 0.25040725 0.96486076
0.66218192 0.35287152 0.51459701
0.48986783 0.22401662 0.28665802
0.65940134 0.00054233 0.84406385
0.99725337 0.65154960 0.83499460
0.33973023 0.00395889 0.13856936
0.51860538 0.17877751 0.00498267
0.76337083 0.36476411 0.23182904
0.85543558 0.53458809 0.03798978
0.16449753 0.62721478 0.88098785
0.39854648 0.62563842 0.65745894
0.61309740 0.83529930 0.87512616
0.37427231 0.83711848 0.09872232
0.38977842 0.61630578 0.32395651
0.18006949 0.40397070 0.10694487
0.61556880 0.39365097 0.67253106
0.83798928 0.61036800 0.88777481
0.64096697 0.16866303 0.89047296
0.40196983 0.16620429 0.11945109
0.85497734 0.37373749 0.09656525
1 f = -11.656880025272354 eV/A^2
-0.003479 -0.000154 -0.071385
0.198313 0.189240 0.024508
-0.230294 -0.821940 -0.116528
-0.057661 -0.023041 0.015731
0.007894 -0.001845 -0.087772
0.033910 -0.040434 0.024005
-0.001020 0.004827 0.005196
-0.002863 0.005201 0.001563
-0.014190 0.027443 0.015277
-0.038084 0.044912 0.005268
0.002526 0.005067 0.001345
-0.033048 0.043260 0.032444
-0.000214 -0.001082 0.000813
-0.015209 0.002622 0.007661
-0.011088 0.019041 0.009966
-0.003106 0.003992 -0.007329
0.104476 0.047465 -0.010943
0.009181 0.009746 -0.022807
-0.000116 -0.001618 0.002201
-0.001211 0.013858 -0.041936
-0.001139 -0.001120 0.001266
-0.008478 0.015654 0.038834
0.002254 0.004035 -0.000964
-0.001452 0.009314 0.005405
0.072678 0.329333 0.048268
0.001941 0.014120 -0.013819
-0.000743 0.001393 0.001505
-0.011640 -0.010003 0.040838
0.017790 0.024851 -0.061398
0.010982 0.006181 0.034219
0.002596 0.003796 0.001566
-0.000052 0.003052 0.004344
-0.019100 0.033341 0.014489
-0.007872 0.006197 -0.001482
-0.010095 0.011718 0.013034
0.003427 -0.001429 0.003368
0.014367 -0.005886 -0.002331
-0.000096 0.012640 0.000489
-0.000033 -0.000137 0.000440
-0.006826 0.004068 -0.000523
-0.000259 0.014785 0.081331
-0.002964 -0.002463 0.003843
! Vibrational analysis
NSW=5000 ! Number of ionic steps
EDIFFG = -0.005 ! Break condition for structure optimization
IBRION = 44 ! Improved dimer method
! Machine learning
ML_LMLFF = T ! Enables the use of MLFF
ML_MODE = run ! Prediction-only mode
K-Points
0
Gamma
1 1 1
0 0 0
Pseudopotentials of Al, C, H, O, and Si.
The KPOINTS and POTCAR files remain as in Example 5. The POSCAR.init file is identical to the POSCAR in Example 5, except for the addition of the forces from the strongest imaginary frequency, which is followed toward the transition state. The forces are taken from the hesseModes.dat file that was produced by vibFreq.py in the previous exercise.
Most of the INCAR tags are the same as in Example 5. The IDM calculation is defined using NSW where it now specifies the number of ionic steps in the structure relaxation and IBRION is set to 44, which means that the improved dimer method (IDM) is used. There is one new tag: EDIFFG, which defines the break condition for structure optimization; when negative, it means that the forces are used to achieve convergence, set in terms of eV Å$^{-1}$.
6.3 Calculation¶
Open a terminal and navigate to this example's directory. Prepare your own POSCAR file, using the vibrational mode for the hardest negative frequency that you calculated in Example 5. The initial forces associated with the hardest negative frequency are taken from hesseModes.dat generated by vibFreq.py in the previous exercise. These forces are then added to the end of the POSCAR file using the following script:
cd $TUTORIALS/transition_states/e06_*
ln -s $TUTORIALS/transition_states/mlff/MLFF_e05-12/ML_FFN ML_FF
cp ../e05_*/POSCAR .
head -43 ../e05_*/hesseModes.dat >> POSCAR
from py4vasp import Calculation
calculation = Calculation.from_path("e05_Static_Vib_analysis")
with open("e06_Static_TS_Opt/POSCAR", "w") as file:
initial_forces = "\n".join(
" ".join(f"{x:12.8f}" for x in displ)
for displ in calculation.force_constant.eigenvectors()[0]
)
file.write(f"""\
{calculation.structure.to_POSCAR()})
initial forces
{initial_forces}""")
Run the IDM calculation using the following:
mpirun -np 4 vasp_gam
During the transition state relaxation, the hardest imaginary mode is followed to the transition state. Let us take a look at how the maximum gradient of the PES changes during the relaxation. First, take the gradient from the OUTCAR file using the following:
cd $TUTORIALS/transition_states/e06_*
grep grad OUTCAR | awk '{print $4}' > grad
In most structure relaxations, it is the forces, i.e. the gradient, that is minimized to achieve convergence. It is easier to see how the gradient changes by plotting it:
import numpy as np
import py4vasp
grad = np.loadtxt("./e06_Static_TS_Opt/grad")
step = np.arange(len(grad))
py4vasp.plot(step, grad, xlabel="Ionic step", ylabel="Max. gradient in eV/Å", label="Gradient")
The gradient rapidly decreases in the first 100 steps. It oscillates as it converges, with little change after 200 ionic steps. As the calcualtion goes on for 800 steps, this raises the question as to why. Using the following, you can find how the curvature λ of the PES changes. The curvature is only calculated every four steps, after which the gradient changes significantly. At the TS, the curvature should be minimised. First, find the curvature in the OUTCAR file using the following:
cd $TUTORIALS/transition_states/e06_*
grep -A 10 curvature OUTCAR > out
grep curvature out | awk '{print $6}' > curv
grep Ionic out | awk '{print $4}' > steps
Then plot the curvature with the following script:
import numpy as np
import py4vasp
grad = np.loadtxt("./e06_Static_TS_Opt/grad")
step = np.arange(len(grad))
curv = np.loadtxt("./e06_Static_TS_Opt/curv")
step2 = np.loadtxt("./e06_Static_TS_Opt/steps")
py4vasp.plot(step, grad, y2=True, y2label="Max. gradient in eV Å<sup>-1</sup>", label="Gradient")\
+ py4vasp.plot(step2, curv, xlabel="Ionic step", ylabel="Curvature of PES in eV Å<sup>-2</sup>", label="Curvature")
In contrast to the gradient, the curvature increases in the first few updates. The oscillation of the curvature coincides with the oscillation of the gradient. However, while the gradient rapidly decreases, the curvature takes many more steps. The large number of steps required to converge the curvature explains why the calculation continues to ~600 steps, even though the forces rapidly converge within 200 steps. The curvature is the limiting factor. Eventually, the curvature plateaus and the calculation converges.
Congratulations! You have converged the transition state. Take a look at the converged transition state structure using py4vasp, and then compare to the initial guess:
import py4vasp
mycalc = py4vasp.Calculation.from_path("./e05_Static_Vib_analysis/")
mycalc.structure.plot(supercell=(2,2,1))
import py4vasp
mycalc = py4vasp.Calculation.from_path("./e06_Static_TS_Opt/")
mycalc.structure.plot(supercell=(2,2,1))
By converging the transition state, you have confirmed that the initial guess was a reasonable one. Note that, using the Hessian, instead of performing IDM, would have required a frequency calculation at every step and vastly increased the calculation time. IDM, on the other hand, only requires a few single point calculations per step, making it far more efficient.
Visually, there is only a small difference between the initial and converged structures. Calculate the C-C and C-O distances to make the difference clearer:
from ase.io.vasp import read_vasp
from ase import geometry
def bond_length(POSCAR):
cell = POSCAR.get_cell()
c1 = (POSCAR.get_positions()[1])
c2 = (POSCAR.get_positions()[2])
o = (POSCAR.get_positions()[24])
return(geometry.get_distances(c1, c2, cell=cell, pbc=True), geometry.get_distances(c2, o, cell=cell, pbc=True))
#P = read_vasp("./e05_Static_Vib_analysis/CONTCAR") # Read the structure to POSCAR
calc = py4vasp.Calculation.from_path("./e05_Static_Vib_analysis")
P = calc.structure.to_ase()
cc_P, co_P = bond_length(P)
print("The C-C bond length in the initial structure is: " + str("{:.2f}".format(cc_P[1][0][0])) + " Å")
print("The C-O bond length in the initial structure is: " + str("{:.2f}".format(co_P[1][0][0])) + " Å")
#C = read_vasp(loc + "/e06_Static_TS_Opt/CONTCAR") # Read the structure to POSCAR
calc = py4vasp.Calculation.from_path("./e06_Static_TS_Opt")
C = calc.structure.to_ase()
cc_C, co_C = bond_length(C)
print("The C-C bond length in the converged structure is: " + str("{:.2f}".format(cc_C[1][0][0])) + " Å")
print("The C-O bond length in the converged structure is: " + str("{:.2f}".format(co_C[1][0][0])) + " Å")
The bond lengths both contract from ~2.1 Å to ~2.0 Å. The contraction brings the reactants closer together in the transition state.
In the next exercise, you will check that your converged structure is a transition state, which should contain only one imaginary vibrational mode, by running a frequency calculation, as in Example 5.
6.4 Questions¶
- How many ionic steps does it take before the curvature is updated?
- What advantage does IDM have over recalculating the Hessian each time?
By the end of this tutorial, you will be able to:
- check that the relaxed structure is a transition state
- calculate an unstable trial direction for calculating an intrinsic reaction coordinate (IRC)
- calculate the free energy of the transition state
7.1 Task¶
Perform vibrational analysis on a relaxed transition state for the carbonylation of methylated acidic chabazite to check for a first-order saddle point and calculate the free energy of the transition state at 440 K
A first-order saddle point should contain only one imaginary mode, excluding vanishing translational modes. The imaginary mode may then be used as an unstable trial direction for finding the intrinsic reaction coordinate (IRC). To this end, vibrational analysis must be performed. Additionally, the total energy, including the vibrational contribution, is then calculated, which is important for comparison to experiment. One important thermodynamic property will be calculated, the zero-point vibrational energy $E_{ZPV}$:
$$\tag{7.1} E_{ZPV} = \sum_{i}\frac{\hbar}{2}\nu_i $$
where $\nu_i$ is the frequency of vibrational mode $i$.
Other important vibrational properties are also calculated, including the vibrational enthalpy $H_{vib}$ at temperature $T$, and the vibrational free energy $A_{vib}$. Note that, here, the free energy is Helmholtz, not Gibbs, as a canonical (NVT) ensemble is used so the temperature is constant. The canonical ensemble also fixes the number of particles $N$ and the volume $V$, alongside the temperature $T$. These thermodynamic properties will be calculated and added to the total energy $E_{tot}$ to obtain the total Helmholtz free energy of the transition state $A_{tot}$, which is required later for calculating the rate constant.
In this task, the optimized structure from IDM in Example 6 will be tested to see if it is truly a transition state. This is done by performing vibrational analysis, as in Example 5. There should be only one imaginary mode, which will be used to calculate the intrinsic reaction coordinate (IRC) in the next exercise. Several thermodynamic properties, including the Helmholtz free energy will also be calculated.
7.2 Input¶
The input files to run this example are prepared at $TUTORIALS/transition_states/e07_Static_Saddle_point.
POSCAR.init
Click here to see the POSCAR file!
Trial TS structure
1.00000000000000
7.9487548458510000 -0.0000499677240000 4.9014331490860004
-3.9745640672680000 6.8838943106960002 4.9012547851189998
-3.9744387686199998 -6.8837771420270002 4.9012457066059998
Al C H O Si
1 2 3 25 11
Direct
0.6207250673754521 0.3699890231127426 0.3368050128480131
0.3518672888673914 0.0092222881033168 0.6127479085897781
0.4287803573375186 0.1041848337269776 0.4472913878224521
0.3728475189445214 0.1963560296615889 0.4842470059654446
0.3742489080687871 0.0172358644669198 0.3715904391464536
0.5405873763029050 0.1000472131776645 0.4890137942784961
0.2393597916469515 0.7434406779337375 0.9997950787522755
0.1771256833804921 0.4606038234744174 0.9514145657813884
0.2554573625922323 0.6124761433147156 0.7381789652423850
0.5186279850644584 0.7491539604912466 0.7299961558202187
0.0187253486587985 0.3357869103147398 0.1404148732436977
0.3555496783712448 0.6531039250594723 0.4900712904111116
0.2454726650308809 0.5143872754641907 0.2388682553068515
0.4078382092232326 0.7625878423429384 0.2511057532749806
0.4674618342053629 0.4632347070226297 0.6576694642513465
0.7651963406602207 0.7555455130890456 0.9072880305638334
0.2818667830153339 0.9842650615926847 0.7032369301338822
0.5839152254814785 0.2388916071603572 0.7463975085851604
0.5238099343981981 0.8220041275022747 0.0127693448356100
0.2739768288373348 0.2566915598560014 0.0966914094606334
0.7425425164389633 0.4961346350681660 0.7673974125174176
0.5289314853972662 0.5218894548799863 0.3146372581307010
0.7907409470183825 0.2474051672494902 0.9621552251444515
0.6739299377530010 0.3480019031201989 0.5147975285743343
0.4933925128976741 0.2108488740121813 0.2854916803184149
0.6561358360120118 0.9958149671152241 0.8423207084968691
0.0015282256971498 0.6413780847725588 0.8343363134240118
0.3407926216169423 0.9934575600094236 0.1284113709199748
0.5138647351344992 0.1752869986704866 0.0007349453151985
0.7596335668804375 0.3598422121023206 0.2295632966832555
0.8564433162545910 0.5312341698323186 0.0385841279427200
0.1680717732208178 0.6135308418499742 0.8808418438404609
0.4007141979740237 0.6181451829996698 0.6544318079249416
0.6161521377862831 0.8288877685613528 0.8732754990082748
0.3786252380728058 0.8275249542157562 0.0966907180338860
0.3893340994836090 0.6081040192621107 0.3231393064146926
0.1797122314228319 0.3940255391973485 0.1081099497353717
0.6206017430523201 0.3880399602702920 0.6692452820262623
0.8413953207637753 0.6051396022864584 0.8873054337114012
0.6387159188343304 0.1647709620108561 0.8906870003945679
0.4063202519763102 0.1602970653458495 0.1262618622144166
0.8530290947887168 0.3694224373307244 0.0965837819449348
! Vibrational analysis
IBRION = 5 ! Finite differences
NFREE = 2 ! Number of displacements
POTIM = 0.01 ! Width of displacement
NSW=1 ! Number of ionic steps
! Machine learning
ML_LMLFF = T ! Enables the use of MLFF
ML_MODE = run ! Prediction-only mode
K-Points
0
Gamma
1 1 1
0 0 0
Pseudopotentials of Al, C, H, O, and Si.
7.3 Calculation¶
Open a terminal and navigate to this example's directory. Copy the transition state structure that you relaxed in Example 6.
cd $TUTORIALS/transition_states/e07_*
ln -s $TUTORIALS/transition_states/mlff/MLFF_e05-12/ML_FFN ML_FF
cp ../e06_*/CONTCAR POSCAR
The vibrational analysis is performed identically to Example 5. Run it using the following:
mpirun -np 4 vasp_gam
The following script extracts the frequencies from the output and then determines the Hessian and dynamical matrices. The Hessian and dynamical matrices' corresponding eigenvalues and eigenfunctions are written in several different formats. Process the frequencies using the following script:
cd $TUTORIALS/transition_states/e07_*
python3 ../scripts/vibFreq/vibFreq.py 0.0 0.2 0.0 > freq.dat
%%bash
echo "Frequencies before IDM:"
grep cm-1 ./e05*/freq.dat | head -n 8
echo " "
echo "Frequencies after IDM:"
grep cm-1 ./e07*/freq.dat | head -n 8
If you look at the imaginary modes, there are now fewer. There were three imaginary modes above 100 cm-1 before the IDM and two additional smaller ones. After the IDM, there is only one imaginary mode above 100 cm-1 and a few below 10 cm-1. These final remaining modes are vanishing frequencies, due to translations. The translational modes should be removed manually before calculating vibrational thermodynamic properties. However, first let us take a look at the remaining vibrational modes. vibFreq.py also extracted the vibrational frequencies and prepared files to visualize them. One of these files is the hesseModes_shifted.molden. It contains the information to visualize the vibrational modes. Alternatively, look at the vibrational frequencies and prepare molden files using py4vasp:
from py4vasp import Calculation
calculation = Calculation.from_path("e07_Static_Saddle_point")
calculation.phonon.mode.print()
from py4vasp import Calculation
calculation = Calculation.from_path("e07_Static_Saddle_point")
with open("e07_Static_Saddle_point/transition_state.molden", "w") as file:
file.write(calculation.force_constant.to_molden())
Now let us consider the thermodynamics properties. The script vibpartition_eV.py removes the translational modes and calculates the harmonic vibrational contributions to the enthalpy, entropy, and free energy at 440 K (close to experiment):
cd $TUTORIALS/transition_states/e07_*
python3 ../scripts/thermo/vibpartition_eV.py 440. freq.dat > thermo.dat
Take a look at thermo.dat to see the calculated theromdynamic properties:
cat $TUTORIALS/transition_states/e07_*/thermo.dat
%%bash
cat ./e07_*/thermo.dat
The file format of thermo.dat is as follows:
| tag | meaning |
|---|---|
nDOF |
Degrees of freedom. This is three times the number of atoms |
nVIB |
The number of vibrational modes, once the imaginary frequencies have been removed. |
line beneath nVIB |
The vibrational partition function Q_vib is shown in scientific notation without being explicitly labeled. |
T |
The temperature in K |
ZPV |
The zero-point vibrational energy |
Hvib - ZPV |
The thermal contribution to the vibrational energy, i.e. how the vibrational energy changes due to vibrational modes above the ground state being populated at temperature. |
Hvib |
The vibrational enthalpy |
Svib |
The vibrational entropy |
Gvib |
The vibrational Helmholtz free energy |
These thermodynamic contributions can also be calculated using the thermochemistry module in ASE:
from ase.thermochemistry import HarmonicThermo
from ase.units import invcm
def parse_frequencies_cm1(filepath):
frequencies = []
a = 0
with open(filepath, 'r') as f:
for line in f:
a += 1
if 'f =' in line:
try:
parts = line.split()
cm1 = float(parts[3]) # 4th column: cm⁻¹
if cm1 > 10**(-3): # skip imaginary and vanishing modes
frequencies.append(cm1)
except (ValueError, IndexError):
continue
return frequencies, a
# Step 1: Read frequencies in cm⁻¹
freq_cm1, n = parse_frequencies_cm1("./e07_Static_Saddle_point/freq.dat")
# Step 2: Convert cm⁻¹ to eV
freq_eV = [f * invcm for f in freq_cm1]
# Step 3: Use ASE to compute vibrational free energy
thermo = HarmonicThermo(freq_eV)
T=440.0
F_vib = thermo.get_helmholtz_energy(temperature=T)
ZPE = thermo.get_ZPE_correction()
U = thermo.get_internal_energy(temperature=T)
S = thermo.get_entropy(temperature=T)
# Result
print("\nResults in 'thermo.dat' style:")
print("nDOF: " +str(n))
print("nDOF: " +str(len(freq_cm1)))
print("N/A\n")
print("T [K] ZPE [eV] Hvib-ZPE [eV] Hvib [eV] Svib [eV/K] Gvib [eV]")
print(str(T) + " " + str(ZPE) + " " + str(U-ZPE) + " " + str(U) + " " + str(S) + " " + str(F_vib))
Gvib is the vibrational free energy. Combining the vibrational free energy with the total energy from the OUTCAR gives the Helmholtz free energy of the transition state $\Delta A_{R\rightarrow P}^{\ddagger}$, calculated using the following script:
import py4vasp
def thermo(path, state, data):
i = 0
A_vib = 0.0
with open(path + data) as f:
for line in f:
pass
A_vib = float(line.split()[5])
calc = py4vasp.Calculation.from_path(path)
E_tot = float(calc.energy.read()["free energy TOTEN"])
A_tot = E_tot + A_vib
print("Total energy of " + str(state) +": " + str("{:.3f}".format(E_tot)) + " eV")
print("Vibrational free energy of " + str(state) +": " + str("{:.3f}".format(A_vib)) + " eV")
print("Helmholtz free energy of " + str(state) +": " + str("{:.3f}".format(A_tot)) + " eV")
return(A_vib, E_tot, A_tot)
A_tot_TS = thermo("./e07_Static_Saddle_point/", "TS", "thermo.dat") # TS
Congratulations! You have successfully calculated the free energy for the transition state. Free energy in isolation is useful but is not measured in experiment. Instead, the free energy of the transition state is needed to calculate an energy difference, the free energy of activation, which can be measured experimentally.
7.4 Questions¶
- What type of mode are the remaining vanishing vibrational modes?
- How is the zero-point vibrational energy calculated?
By the end of this tutorial, you will be able to:
- follow an unstable trial direction along the intrinsic reaction coordinate (IRC)
- find and visualize the reactant and product states
8.1 Task¶
Follow the imaginary vibrational mode of the transition state during carbonylation step in chabazite as a trial vector toward the reactants and products to produce the intrinsic reaction coordinate
The transition state sits in-between the reactants and products. The saddle point contains only one imaginary frequency, a single unstable vibrational mode. The TS structure can be propagated along this unstable mode to generate the intrinsic reaction coordinate (IRC). Following the unstable vibrational mode forward leads to the product and backward to the reactant. This effectively checks if the transition state connects our expected reactant and product. A damped-velocity Verlet algorithm is used to propagate along the unstable mode [J. Chem. Phys. A 106, 165 (2002)]. More detail about IRC is given in Example 4.
In this exercise, you will produce the IRC for the carbonylation reaction in chabazite from Example 7. By checking that the reactants and products are what we expect, the transition state is found to link the acetyl cation & Z-, and CO & CH3-Z.
8.2 Input¶
The input files to run this example are prepared at $TUTORIALS/transition_states/e08_Static_Reaction_pathway.
POSCAR.init
Click here to reveal the POSCAR file!
Trial TS structure
1.00000000000000
7.9487548458510000 -0.0000499677240000 4.9014331490860004
-3.9745640672680000 6.8838943106960002 4.9012547851189998
-3.9744387686199998 -6.8837771420270002 4.9012457066059998
Al C H O Si
1 2 3 25 11
Direct
0.6207250673754521 0.3699890231127426 0.3368050128480131
0.3518672888673914 0.0092222881033168 0.6127479085897781
0.4287803573375186 0.1041848337269776 0.4472913878224521
0.3728475189445214 0.1963560296615889 0.4842470059654446
0.3742489080687871 0.0172358644669198 0.3715904391464536
0.5405873763029050 0.1000472131776645 0.4890137942784961
0.2393597916469515 0.7434406779337375 0.9997950787522755
0.1771256833804921 0.4606038234744174 0.9514145657813884
0.2554573625922323 0.6124761433147156 0.7381789652423850
0.5186279850644584 0.7491539604912466 0.7299961558202187
0.0187253486587985 0.3357869103147398 0.1404148732436977
0.3555496783712448 0.6531039250594723 0.4900712904111116
0.2454726650308809 0.5143872754641907 0.2388682553068515
0.4078382092232326 0.7625878423429384 0.2511057532749806
0.4674618342053629 0.4632347070226297 0.6576694642513465
0.7651963406602207 0.7555455130890456 0.9072880305638334
0.2818667830153339 0.9842650615926847 0.7032369301338822
0.5839152254814785 0.2388916071603572 0.7463975085851604
0.5238099343981981 0.8220041275022747 0.0127693448356100
0.2739768288373348 0.2566915598560014 0.0966914094606334
0.7425425164389633 0.4961346350681660 0.7673974125174176
0.5289314853972662 0.5218894548799863 0.3146372581307010
0.7907409470183825 0.2474051672494902 0.9621552251444515
0.6739299377530010 0.3480019031201989 0.5147975285743343
0.4933925128976741 0.2108488740121813 0.2854916803184149
0.6561358360120118 0.9958149671152241 0.8423207084968691
0.0015282256971498 0.6413780847725588 0.8343363134240118
0.3407926216169423 0.9934575600094236 0.1284113709199748
0.5138647351344992 0.1752869986704866 0.0007349453151985
0.7596335668804375 0.3598422121023206 0.2295632966832555
0.8564433162545910 0.5312341698323186 0.0385841279427200
0.1680717732208178 0.6135308418499742 0.8808418438404609
0.4007141979740237 0.6181451829996698 0.6544318079249416
0.6161521377862831 0.8288877685613528 0.8732754990082748
0.3786252380728058 0.8275249542157562 0.0966907180338860
0.3893340994836090 0.6081040192621107 0.3231393064146926
0.1797122314228319 0.3940255391973485 0.1081099497353717
0.6206017430523201 0.3880399602702920 0.6692452820262623
0.8413953207637753 0.6051396022864584 0.8873054337114012
0.6387159188343304 0.1647709620108561 0.8906870003945679
0.4063202519763102 0.1602970653458495 0.1262618622144166
0.8530290947887168 0.3694224373307244 0.0965837819449348
-0.39544783E-03 -0.31588388E-02 -0.62676293E-01
0.69700337E-01 0.31046523E+00 0.51024466E-01
-0.32033798E+00 -0.79704463E+00 -0.30556313E-01
-0.24372420E-02 -0.37441346E-02 -0.79073743E-02
-0.14900716E-02 -0.16413596E-01 -0.28841066E-01
-0.16790196E-01 -0.46140567E-01 0.13350218E-02
0.98888602E-03 0.55323617E-02 0.29823640E-02
-0.21022934E-02 0.37726453E-02 -0.71679080E-03
-0.59507096E-02 0.13858925E-01 0.74783275E-02
-0.76190413E-02 0.12542317E-01 0.99451166E-03
0.20372207E-02 0.43807918E-02 0.96260335E-03
-0.11590822E-01 0.12467761E-01 0.10392734E-01
-0.83560198E-03 -0.75644598E-03 -0.54237396E-03
-0.95769089E-02 -0.10495125E-02 0.12904105E-02
-0.24540580E-02 -0.25671757E-02 0.10883893E-03
0.14890613E-02 0.10163308E-02 -0.41785947E-03
0.19613558E+00 0.16751863E+00 0.11594418E-01
0.57738947E-02 0.40463377E-02 -0.21835018E-01
0.11751591E-02 -0.16283564E-02 -0.11283298E-02
-0.99130142E-04 0.15902318E-01 -0.42495500E-01
-0.13635877E-02 -0.11271898E-02 0.67467050E-03
-0.25995086E-02 0.13308049E-01 0.35478105E-01
0.33858715E-02 0.49288244E-02 -0.18371861E-02
0.43203675E-04 0.28468862E-03 -0.36608873E-02
0.64540177E-01 0.24578122E+00 0.20915201E-01
0.23777336E-02 0.90000794E-02 -0.12908098E-01
0.76893028E-03 0.28983856E-02 0.15231467E-02
-0.77832548E-02 -0.61445446E-02 0.22289668E-01
0.14253429E-01 0.30174731E-01 -0.68771139E-01
0.91471680E-02 0.55764505E-02 0.33381680E-01
0.37034454E-02 0.32814988E-02 0.97077706E-03
-0.52749231E-03 -0.14345845E-02 0.29907403E-02
0.22117960E-02 0.74846929E-03 -0.62356316E-04
-0.30623346E-03 -0.33004676E-02 0.28001697E-02
-0.83263851E-02 0.78906662E-02 0.10928800E-01
0.10125808E-01 -0.28894395E-02 -0.19362441E-02
0.11978396E-01 -0.53758879E-02 -0.21288116E-02
0.66990984E-02 0.83533895E-02 0.67978420E-03
0.16413788E-02 -0.47433499E-03 0.10692095E-02
-0.46741496E-02 -0.12486586E-02 0.27702860E-03
0.21378568E-02 0.14270053E-01 0.63286845E-01
-0.30542756E-02 -0.35017912E-02 0.29921201E-02
! IRC calculation
NSW=5000 ! Number of ionic steps
IBRION = 40 ! IRC
IRC_VNORM0 = 0.0005 ! constant velocity vector in IRC
IRC_DIRECTION = -1 ! negative initial direction of movement
! Machine learning
ML_LMLFF = T ! Enables the use of MLFF
ML_MODE = run ! Prediction-only mode
K-Points
0
Gamma
1 1 1
0 0 0
Pseudopotentials of Al, C, H, O, and Si.
The KPOINTS and POTCAR files remain as in Example 5. The POSCAR.init file is taken from an IDM relaxation, equivalent to Example 6. The INCAR file is similar to that used in Example 6, with the change to IBRION = 40, which sets the calculation to follow the energy profile along the IRC from the transition state using a damped-velocity Verlet algorithm.
There is one new tag: IRC_VNORM0, which defines the constant value that the velocity vector is rescaled to, thereby providing the damping. IRC_DIRECTION defines whether the unstable mode is followed forward to the product (IRC_DIRECTION = 1) or backward to the reactant (IRC_DIRECTION = -1).
8.3 Calculation¶
Starting from the transition state, the intrinsic reaction coordinate (IRC) models the reaction pathway. It is calculated by following a trial unstable mode, taken from Example 7. The unstable mode is followed in two directions, forward to the product (plus, in the p directory) and backward to the reactant (minus, in the m directory). A separate calculation is required for each, the reactant is calculated in m (minus) and the product in p (plus). First, run the calculation for the reactant with the following in your terminal:
mkdir m
cp POSCAR POTCAR KPOINTS INCAR m
cd ./m
ln -s $TUTORIALS/transition_states/mlff/MLFF_e05-12/ML_FFN ML_FF
mpirun -np 4 vasp_gam
Second, calculate the product in the same manner, using the following:
cd $TUTORIALS/transition_states/e08_*
mkdir p
cp POSCAR POTCAR KPOINTS INCAR p
cd ./p
sed -i 's/IRC_DIRECTION = -1/IRC_DIRECTION = 1/g' INCAR
ln -s $TUTORIALS/transition_states/mlff/ML_FFN ML_FF
mpirun -np 4 vasp_gam
You have succeeded in calculating the IRC! Let us extract this potential profile using the ircShift.sh script. This script extracts the IRC coordinate in Å IRC (A): from the OUTCAR files in the m and p directories and merges them to get the entire reaction profile. Extract the IRC using the following:
cd $TUTORIALS/transition_states/e08_*
../scripts/bashScripts/ircShift.sh
ircShift.sh outputs ircShift.dat with all of the IRC points. Extract the IRC data and plot it with:
import py4vasp
import numpy as np
irc = np.genfromtxt("./e08_Static_Reaction_pathway/ircShift.dat",unpack=True,usecols=range(2))
py4vasp.plot(irc[0], irc[1], xlabel="Intrinsic reaction coordinate in Å", ylabel="Energy in eV", label="IRC")
Can you label the different parts of the profile?
At the beginning, ~-3 Å, is the reactant with the initial CO and CH${3}$-Z. As the CO comes closer and the CH${3}$ inverts, the energy increases to a maximum, at 0 Å, the transition state. The free acetyl cation forms before the energy rapidly decreases to the product, ~4 Å.
The reaction profile is made clearer by visualizing the structures of the reactant, transition state, and product. You can use py4vasp to look at each of the structures:
import py4vasp
R = py4vasp.Calculation.from_path("./e08_Static_Reaction_pathway/m/")
TS = py4vasp.Calculation.from_path("./e06_Static_TS_Opt/")
P = py4vasp.Calculation.from_path("./e08_Static_Reaction_pathway/p/")
# Take a look at the structure that your interested in by unhashing the relevant line
#R.structure.plot(supercell=(2,2,1)) # Reactant
#TS.structure.plot(supercell=(2,2,1)) # Transition state
R.structure.plot(supercell=(2,2,1)) # Product
import py4vasp
R = py4vasp.Calculation.from_path("./e08_Static_Reaction_pathway/m/")
TS = py4vasp.Calculation.from_path("./e06_Static_TS_Opt/")
P = py4vasp.Calculation.from_path("./e08_Static_Reaction_pathway/p/")
TS.structure.plot(supercell=(2,2,1)) # Product
import py4vasp
R = py4vasp.Calculation.from_path("./e08_Static_Reaction_pathway/m/")
TS = py4vasp.Calculation.from_path("./e06_Static_TS_Opt/")
P = py4vasp.Calculation.from_path("./e08_Static_Reaction_pathway/p/")
P.structure.plot(supercell=(2,2,1)) # Product
The scripts pos2xyzAnim_rev.py and pos2xyzAnim_for.py extract the IRC coordinate every 10 steps for the backward and forward reaction, respectively. These scripts create a video of the reaction. Run them one after the other and stick them together using the following:
cd $TUTORIALS/transition_states/e08_*
python3 ../scripts/pos2xyzNVT/pos2xyzAnim_rev.py m/XDATCAR 10 0 0.2 0 > anim.xyz
python3 ../scripts/pos2xyzNVT/pos2xyzAnim_for.py p/XDATCAR 10 0 0.2 0 >> anim.xyz
You can view the animated the reaction using the anim.xyz file in JMol. A recording of anim.xyz is given below:
Click here to watch a video of the reaction from reactant to transition state to product!
In the case of both the reactant and product, the IRC calculation terminates once the energy monotonously increases in a sequence of IRC_STOP number of steps. IRC_STOP is set to 20 by default, so the 20th structure from the end of the calculation is the closest to the potential energy minimum. We can extract these from XDATCAR, where they are stored, using the following pickMin_irc.sh:
cd $TUTORIALS/transition_states/e08_*/m
../../scripts/bashScripts/pickMin_irc.sh 20
cd ../p
../../scripts/bashScripts/pickMin_irc.sh 20
This creates the POSCAR.pick files in the m and p directories. Due to the small IRC_VNORM0 value used in our calculation, 0.0005 eV Å-1, this lies very close to the minimum and so no further relaxation is needed. In the next exercise, you will calculate the free energy of the reactant state and the activation free energy of the reaction.
8.4 Questions¶
- What structures are given at the endpoints of
IRC_DIRECTION = 1andIRC_DIRECTION = -1?
By the end of this tutorial, you will be able to:
- do vibrational analysis on the reactant structure from an IRC calculation
- compare the harmonic to semi-classical approaches
- calculate the free energy of the reaction
- calculate the rate constant
9.1 Task¶
Perform vibrational analysis on the reactant (methylated chabazite and CO) and calculate the corresponding activation free energy and rate constant for the carbonylation reaction
Total energies are calculated ab initio. The free energy, or another thermodynamic property, is measured by experiment. The comparison of the two requires the inclusion of a few different thermodynamic terms that can be calculated from the vibrational frequencies. The thermodynamic terms including the zero-point vibrational energy, the thermal contributions to the vibrational energy, and the entropic terms must be accounted for. The free energy difference between the reactant and transition states, i.e. the free energy of reaction $\Delta A_{R\rightarrow P}^{\ddagger}$, is a commonly measured property. If $\Delta A_{R\rightarrow P}^{\ddagger}$ is inserted into the Eyring equation, then a common descriptor of a reaction is found, the rate constant $k_{R\rightarrow P}$:
$$\tag{9.1} \large k_{R\rightarrow P} = \frac{k_{B}T}{h}e^{-\left(\Delta A_{R\rightarrow P}^{\ddagger} / k_B T \right)} $$
where $k_{B}$ is the Boltzmann constant, $h$ is the Planck constant, and $T$ the temperature
In this exercise, you will perform a vibrational analysis on the final structure from the backward propagation in Example 8 to check for imaginary modes. If there are no imaginary modes, then the final structure is a minimum and therefore the reactant. Once the reactant has been confirmed, the thermodynamic properties of the reaction will then be calculated, including the free energy of activation and the rate constant.
9.2 Input¶
The input files to run this example are prepared at $TUTORIALS/transition_states/e09_Static_Thermodynamics.
POSCAR.init
Click here to see the POSCAR!
Reactant
1
7.948755 -0.000050 4.901433
-3.974564 6.883894 4.901255
-3.974439 -6.883777 4.901246
Al C H O Si
1 2 3 25 11
Direct configuration= 1886
0.62952771 0.37745230 0.33739735
0.25191620 0.96168432 0.64117336
0.39172183 0.18563249 0.40461206
0.33338056 0.28035823 0.43170585
0.31736040 0.08893041 0.36882835
0.46548653 0.16358203 0.49640954
0.23530815 0.74715842 0.00005085
0.17906592 0.46352972 0.95104251
0.25856085 0.61712179 0.73953578
0.52467333 0.74540501 0.72744140
0.02313926 0.33701354 0.14234199
0.35773985 0.64394317 0.48891507
0.25821808 0.50115223 0.23828089
0.40429114 0.75927752 0.24917835
0.46806931 0.46075972 0.66268681
0.76729001 0.75405826 0.90906636
0.29012599 0.96936993 0.76193801
0.58983685 0.23821733 0.74827852
0.52208925 0.81114731 0.00984135
0.27114317 0.24921266 0.08195343
0.74486799 0.49531144 0.76650684
0.54082448 0.52953054 0.31392885
0.80009719 0.24362299 0.95949031
0.67299378 0.34901524 0.51407640
0.48225178 0.21962111 0.28622622
0.65496241 0.99357097 0.84772717
0.00200725 0.64029850 0.83182829
0.35054171 0.99285230 0.12607265
0.52908697 0.18468256 0.01096371
0.76107650 0.35158722 0.22556486
0.85797485 0.52866751 0.03797078
0.16890258 0.61623378 0.88027912
0.40347826 0.61670822 0.65521549
0.61715759 0.82524264 0.87300247
0.37898630 0.82422021 0.09494387
0.39467014 0.60632644 0.32261103
0.18296570 0.39105114 0.10384741
0.62162725 0.38679851 0.67049796
0.84289471 0.60395996 0.88703743
0.64523535 0.16388889 0.89087957
0.40586901 0.15890126 0.11699808
0.85715870 0.36630157 0.09419744
! Vibrational analysis
IBRION = 5 ! Finite differences
NFREE = 2 ! Number of displacements
POTIM = 0.01 ! Width of displacement
NSW = 1 ! Number of ionic steps
! Machine learning
ML_LMLFF = T ! Enables the use of MLFF
ML_MODE = run ! Prediction-only mode
K-Points
0
Gamma
1 1 1
0 0 0
Pseudopotentials of Al, C, H, O, and Si.
9.3 Calculation¶
Open a terminal and navigate to this example's directory. Copy the POSCAR.pick file that that you calculated in Example 8:
cd $TUTORIALS/transition_states/e09_*
ln -s $TUTORIALS/transition_states/mlff/MLFF_e05-12/ML_FFN ML_FF
cp ../e08_*/m/POSCAR.pick POSCAR
With the calculation prepared, you can now check that m/POSCAR.pick is an energetic minimum by running the following:
mpirun -np 4 vasp_gam
Once the calculation has finished, vibFreq.py extracts the frequencies from the output and then determines the Hessian and dynamical matrices. The matrices' corresponding eigenvalues and eigenfunctions are written out in several different formats. Process and view the frequencies using the following script:
cd $TUTORIALS/transition_states/e09_*
python3 ../scripts/vibFreq/vibFreq.py 0. 0.2 0. > freq.dat
echo "Frequencies of reactant:"
grep cm-1 freq.dat | head -n 5
%%bash
echo "Frequencies of reactant:"
grep cm-1 ./e09_*/freq.dat | head -n 5
Except for vanishing imaginary frequencies corresponding to translations, there are only positive frequencies, so this is, in fact, an energetic minimum! Congratulations on finding the reactant!
Using the same procedure to extract the energies as in Example 7, vibpartition_eV.py calculates the harmonic vibrational contributions to the enthalpy, entropy, and free energy. These properties are calculated at 440 K, which lies close to the experiment using the following script:
cd $TUTORIALS/transition_states/e09_*
python3 ../scripts/thermo/vibpartition_eV.py 440. freq.dat > thermo.dat
cat thermo.dat
%%bash
cat e09_*/thermo.dat
These thermodynamic contributions can also be calculated using the thermochemistry module in ASE:
from ase.thermochemistry import HarmonicThermo
from ase.units import invcm
def parse_frequencies_cm1(filepath):
frequencies = []
a = 0
with open(filepath, 'r') as f:
for line in f:
a += 1
if 'f =' in line:
try:
parts = line.split()
cm1 = float(parts[3]) # 4th column: cm⁻¹
if cm1 > 10**(-3): # skip imaginary and vanishing modes
frequencies.append(cm1)
except (ValueError, IndexError):
continue
return frequencies, a
# Step 1: Read frequencies in cm⁻¹
freq_cm1, n = parse_frequencies_cm1("./e09_Static_Thermodynamics/freq.dat")
# Step 2: Convert cm⁻¹ to eV
freq_eV = [f * invcm for f in freq_cm1]
# Step 3: Use ASE to compute vibrational free energy
thermo = HarmonicThermo(freq_eV)
T=440.0
F_vib = thermo.get_helmholtz_energy(temperature=T)
ZPE = thermo.get_ZPE_correction()
U = thermo.get_internal_energy(temperature=T)
S = thermo.get_entropy(temperature=T)
# Result
print("\nResults in 'thermo.dat' style:")
print("nDOF: " +str(n))
print("nDOF: " +str(len(freq_cm1)))
print("N/A\n")
print("T [K] ZPE [eV] Hvib-ZPE [eV] Hvib [eV] Svib [eV/K] Gvib [eV]")
print(str(T) + " " + str(ZPE) + " " + str(U-ZPE) + " " + str(U) + " " + str(S) + " " + str(F_vib))
The file format of thermo.dat is as follows:
| tag | meaning |
|---|---|
nDOF |
Degrees of freedom. This is three times the number of atoms |
nVIB |
The number of vibrational modes, once the imaginary frequencies have been removed. |
line beneath nVIB |
The vibrational partition function Q_vib is shown without an introduction. |
T |
The temperature in K |
ZPV |
The zero-point vibrational energy |
Hvib - ZPV |
The thermal contribution to the vibrational energy, i.e. how the vibrational energy changes due to vibrational modes above the ground state being populated at temperatire. |
Hvib |
The vibrational enthalpy |
Svib |
The vibrational entropy |
Gvib |
The vibrational Helmholtz free energy |
Gvib is the vibrational free energy. Recalling from thermodynamics that the Helmholtz free energy $A$ is defined from the partition function $Q$:
$$\tag{9.2} A = - k_B T \ln{Q} $$
The quantum-mechanical vibrational partition functional is defined as:
$$\tag{9.3} Q_{\mu,vib} = \prod_{i=a}^{N_{vib}} \Large{\frac{e^{(-h \nu_i)/(2 k_B T)}}{1 - e^{(-h \nu_i)/(k_B T)}}} $$
Combining the vibrational free energy with the total energy from the OUTCAR gives the Helmholtz free energy of the transition state $\Delta A_{R\rightarrow P}^{\ddagger}$. This is calculated using the following script:
import py4vasp
def thermo(path, state, data):
i = 0
A_vib = 0.0
with open(path + data) as f:
for line in f:
pass
A_vib = float(line.split()[5])
calc = py4vasp.Calculation.from_path(path)
E_tot = float(calc.energy.read()["free energy TOTEN"])
A_tot = E_tot + A_vib
print("Total energy of " + str(state) +": " + str("{:.3f}".format(E_tot)) + " eV")
print("Vibrational free energy of " + str(state) +": " + str("{:.3f}".format(A_vib)) + " eV")
print("Helmholtz free energy of " + str(state) +": " + str("{:.3f}".format(A_tot)) + " eV")
return(A_vib, E_tot, A_tot)
A_tot_TS_QM = thermo("./e07_Static_Saddle_point/", "TS", "thermo.dat") # TS
print(" ")
A_tot_R_QM = thermo("./e09_Static_Thermodynamics/", "reactant", "thermo.dat") # reactant
The free energy of activation for the reaction from reactant to transition state $\Delta A_{R\rightarrow P}^{\ddagger}$ is the difference between the free energy of the TS and the reactant $A_{TS} - A_{R}$:
delta_A_QM = A_tot_TS_QM[2] - A_tot_R_QM[2]
print("ΔA^‡ is: " + str(delta_A_QM) + " eV")
Congratulations! You have calculated the free energy of activation for the carbonylation of methylated chabazite!
$\Delta A_{R\rightarrow P}^{\ddagger}$ was calculated assuming the quantum behavior of harmonic oscillators. The semi-classical vibrational partition function is defined as:
$$\tag{9.4} Q_{\mu,vib}^{SC} = \prod_{i=a}^{N_{vib}} \frac{h \nu_i}{2 k_B T} $$
By inputting $Q_{\mu,vib}^{SC}$ into $A = - k_B T \ln{Q}$, the semi-classical vibrational free energy for each structure can be calculated. The free energy difference for the reaction can then be found and the difference between the semi-classical and the quantum mechanical contributions calculated. The semi-classical vibrational free energy is calculated using the above equation implemented using vibpartitionClassic_eV.py:
cd $TUTORIALS/transition_states
echo "Reactant: "
cd ./e09_*
python3 ../scripts/thermo/vibpartitionClassic_eV.py 440. freq.dat > thermo2.dat
cat thermo2.dat
cd ..
echo " "
echo "TS: "
cd ./e07_*
python3 ../scripts/thermo/vibpartitionClassic_eV.py 440. freq.dat > thermo2.dat
cat thermo2.dat
%%bash
echo "Reactant: "
cd ./e09_*
python3 ../scripts/thermo/vibpartitionClassic_eV.py 440. freq.dat > thermo2.dat
cat thermo2.dat
cd ..
echo " "
echo "TS: "
cd ./e07_*
python3 ../scripts/thermo/vibpartitionClassic_eV.py 440. freq.dat > thermo2.dat
cat thermo2.dat
cd ..
Bringing these thermodynamic terms together, $\Delta A_{R\rightarrow P}^{\ddagger}$ is calculated for the quantum-mechanical and semi-classical approximations, alongside their difference with the following script:
print("Semi-classical approximation:")
A_tot_TS_SC = thermo("./e07_Static_Saddle_point/", "TS", "thermo2.dat") # TS
print("The difference in A_tot for TS is: " + str("{:.3f}".format(A_tot_TS_QM[2] - A_tot_TS_SC[2])) + " eV")
print(" ")
A_tot_R_SC = thermo("./e09_Static_Thermodynamics/", "reactant", "thermo2.dat") # reactant
print("The difference in A_tot for reactant is: " + str("{:.3f}".format(A_tot_R_QM[2] - A_tot_R_SC[2])) + " eV")
print(" ")
delta_A_SC = A_tot_TS_SC[2] - A_tot_R_SC[2]
print("The ΔA^‡ quantum mechanically is: " + str("{:.3f}".format(delta_A_QM)) + " eV")
print("The ΔA^‡ semi-classically is: " + str("{:.3f}".format(delta_A_SC)) + " eV")
print("The difference in ΔA^‡ is: " + str("{:.3f}".format(delta_A_SC - delta_A_QM)) + " eV")
The difference between the semi-classical and quantum mechanical Helmholtz free energies is ~1 eV. Considering this, the difference between the semi-classical and quantum mechanical $\Delta A_{R\rightarrow P}^{\ddagger}$ at ~0.02 eV is extremely small.
What does this suggest about the importance of quantum effects here?
They are small in magnitude, so for the conditions tested, a semi-classical approach is valid. This is an important validation for using molecular dynamics (MD), which uses a semi-classical approach for calculating the forces, failing to include zero-point vibrational contributions to the energy.
Now that the $\Delta A_{R\rightarrow P}^{\ddagger}$ is calculated, the rate constant can be calculated using the Eyring equation:
import scipy
import math
from decimal import Decimal
def eyring(T, delta_A):
k_b = 8.617333262E-5
h = 4.1356692E-15
A = k_b * T / h
B = math.exp(-delta_A/(k_b * T))
return A*B
k_SC = eyring(440., delta_A_SC)
k_QM = eyring(440., delta_A_QM)
print("Semi-classical rate constant: " + str("{:.3f}".format(k_SC)) + " s^-1")
print("Static QM rate constant: " + str("{:.3f}".format(k_QM)) + " s^-1")
Congratulations on calculating the rate constant! You can see that the difference between the vibrational energies makes a ~40% reduction to the rate constant. Small energy differences, when exponentiated, result in big changes in measurable properties.
You could also calculate the rate constant from the products to the transition state by calculating the vibrational frequencies for the products and following the procedure through again.
9.4 Questions¶
- How important are vibrational quantum effects for calculating $\Delta A_{R\rightarrow P}^{\ddagger}$?
- How big a difference does this make in calculating the rate constant?