Part 3: Dynamic approaches ¶
By the end of this tutorial, you will be able to:
- test a trial reaction coordinate through a molecular dynamics (MD) simulation
- determine the corresponding probability density to find a test reactant structure
10.1 Task¶
Trial a reaction coordinate for the carbonylation of methylated chabazite using molecular dynamics (MD) for 50 ps with a CSVR thermostat at 440 K
Static approaches can be successfully used to model transition states. However, they come with some limitations. The potential energy surface (PES) is covered with many local minima and, for most systems, no individual stationary state will describe the complete system; this limitation is true for reactant, transition state, and product. Instead, a collection of structures, an ensemble, is required. Using an ensemble can make a large difference, e.g. up to ~20 kJ mol$^{-1}$ for the alkene isomerization barrier over acid chabazite [J. Catal. 373, 361 (2019)]. Occasionally, the mechanism itself changes when the temperature increases, e.g. if entropy plays a significant role [J. Chem. Phys. 131, 214508 (2009)], which is not accounted for when using a static approach. Omitting the ensemble by using static approaches would lead to an incorrect mechanism.
Dynamic methods, i.e. molecular dynamics (MD), describe the system as an ensemble. The goal of this exercise is to calculate the free energy of activation from reactant $R$ to product $P$, i.e. the forward reaction, $\Delta A_{R\rightarrow P}^{\ddagger}$ of carbonylation of methylated chabazite using an average of the ensemble. For modeling the reaction, a reaction coordinate $\xi$ must be defined. $\xi$ should be chosen such that it is a continuous function of positions q and that an infinitesimally slow change in $\xi$ should be reversible [Annu. Rev. Phys. Chem. 67, 669 (2016)]. For our test reaction, $\xi$ is the difference between the C-C distance $r_1$ (the CO and the methyl's C atom) and the C-O distance $r_2$ (C of methyl group $r_1$ and the zeolite O atom that it is initially bonded to):
$ \xi = r_{2} - r_{1} $
$\xi$ is defined in the ICONST file as the coordinate of interest to be monitored. The reaction coordinate of the target transitions state (TS) is denoted $\xi^*$. Many different structures will be described by the ensemble produced by MD. For this system and length of MD simulation, the most probable structure in a distribution of the probability density $P(\xi)$ is taken as a candidate reactant structure $\xi_{ref,R}$ from the ensemble with which to initiate constrained MD simulations in the next exercise.
In the first exercise, you will test the trial $\xi$ throughout an MD simulation and find the corresponding probability density $P(\xi)$.
10.2 Input¶
The input files to run this example are prepared at $TUTORIALS/transition_states/e10_Dynamic_Coordinates.
POSCAR.init
Click here to reveal the POSCAR!
Trial reactant state
1.00000000000000
7.9487550000000002 -0.0000500000000000 4.9014329999999999
-3.9745640000000000 6.8838939999999997 4.9012549999999999
-3.9744389999999998 -6.8837770000000003 4.9012460000000004
Al C H O Si
1 2 3 25 11
Direct
0.6309346970627223 0.3723713210786789 0.3463999640746260
1.0848109095396736 0.9260278535224171 0.2926399489611942
0.4099804020649828 0.1635203676555287 0.4242317412149454
0.3217639325302051 0.2348576083179665 0.4566396286325596
0.3571207839297778 0.0569135138957275 0.3880912662684614
0.5042637473104756 0.1522796345324226 0.5127718153497565
0.2785301849746668 0.7602849027699632 0.0046895740075484
0.1772302022874687 0.4693952000609656 0.9852581067086047
0.2701895253319148 0.5861736361576427 0.7750715916929146
0.4999550776871463 0.7590524526462237 0.7787200954501512
0.0140860260835989 0.3546741314241051 0.1696029735737132
0.3557593935467726 0.6611532190177529 0.5154041033865417
0.2544540206424660 0.5021210680612611 0.2757996298531374
0.4336105837783332 0.7628783191716760 0.2759815074597886
0.5047228982323143 0.4861184669416140 0.6596589197893293
0.7601968957625573 0.7366368704305295 0.8992844184732824
0.9712107641816842 0.9043484375197114 0.2439125282507523
0.6171437387491777 0.2497465581602643 0.7534029161450518
0.5739296819827301 0.8278279551700365 0.0494861809032877
0.2710042303195807 0.2546226406335215 0.0984134113610155
0.7888426250824750 0.4845915011078290 0.7660663861216227
0.5327368875966848 0.5192481298621753 0.3290590563052356
0.8066267552126245 0.2238772764998835 0.9705464699735484
0.7108436841145985 0.3597097154262990 0.5210001414153450
0.4786096839937157 0.2170430498565461 0.2980267575983134
0.6542812898781195 -0.0052757178458664 0.8657157594602567
0.0137154582286103 0.6558600497849584 0.8553106868941003
0.3729927855127887 -0.0060397381058633 0.1501283281707725
0.5247804151255646 0.1933846308614665 1.0081949067191767
0.7485685642234065 0.3548953524999690 0.2178644511158051
0.8619588244859279 0.5180438246315547 0.0500553772661750
0.1769139224733230 0.6246561021678009 0.9086556027213539
0.4093006730411545 0.6265414077660237 0.6819688324030556
0.6250863063940451 0.8291988869736581 0.8926121172812201
0.4135027480658851 0.8262439027347905 0.1166284452618589
0.3971418346322438 0.6124033984750867 0.3488076895944575
0.1812613691210096 0.3944416726227042 0.1373879781765745
0.6550981440839965 0.3961083928036578 0.6695112657124653
0.8539075591819076 0.6054829339157256 0.9054596290126368
0.6417172649716902 0.1673432828528636 0.8991724568616257
0.4086811511995426 0.1627581573867222 0.1272466890763508
0.8515158967503381 0.3511106682195257 0.0926967776942586
INCAR.init
! Molecular dynamics
IBRION = 0 ! Molecular dynamics
ISYM = 0 ! Symmetry setting
NPAR = 2 ! Number of bands treated in parallel
NSW = 10000 ! Number of steps in MD run
POTIM = 1.0 ! Time step
MDALGO = 5 ! CSVR thermostat
TEBEG = 440. ! Temperature
CSVR_PERIOD = 20 ! Time scale of CSVR in terms of the number of MD steps
! Machine learning
ML_LMLFF = T ! Enables the use of MLFF
ML_MODE = run ! Prediction-only mode
ML_OUTPUT_MODE = 0 ! Reduced file output mode
ML_OUTBLOCK = 10 ! Output distance in number of steps
K-Points
0
Gamma
1 1 1
0 0 0
Pseudopotentials of Al, C, H, O, and Si.
R 3 25 0
R 3 2 0
S -1 1 7
Check the tags set in the INCAR file!
The first group of tags defines the molecular dynamics (MD) simulation. The first tag, IBRION = 0, sets up the MD calculation. The next two are not specific to MD: ISYM = 0 switches off symmetry except for the reduced k-point mesh due to time-reversal symmetry, and NPAR determines the number of bands that are treated in parallel. NSW sets the number of steps in an MD run, while POTIM sets the corresponding time step in femtoseconds. MDALGO = 5 sets the thermostat for MD, specifically to the canonical sampling through velocity-rescaling (CSVR) thermostat. CSVR_PERIOD determines the time step for the CSVR thermostat in femtoseconds. TEBEG sets the initial temperature for the MD simulation and the target temperature for the thermostat.
The second set of tags defines the machine-learned force fields (MLFF) used for MD. ML_LMLFF = T turns on MLFF and ML_MODE = run means that a previously trained MLFF is used. ML_OUTPUT_MODE reduces specific file output for MD to increase the performance when using MLFF. ML_OUTBLOCK sets the output for MLFF MD to be printed every n steps. Finally, RANDOM_SEED will be added in by the run scripts; it is a random number generator defining the initial velocities, ensuring that the first few MD steps are reproducible.
The POTCAR file specifies the pseudopotentials used, KPOINTS uses the Γ-point, and ICONST monitors the bonds during an MD simulation, in this case the reaction coordinate $\xi$. The R's are the interatomic distances between atoms C(3) and O(25), and C(2) and C(3); the S is the difference between these two R's. 0 means the coordinate is constrained, while 7 indicates that it is monitored. The POSCAR is the CONTCAR taken from Example 9, i.e. the reactant.
10.3 Calculation¶
Open a terminal, navigate to this example's directory, and prepare the directories using:
cd $TUTORIALS/transition_states/e10_*
ln -s ../mlff/MLFF_e05-12/ML_FFN ML_FF
cp $TUTORIALS/transition_states/e09_*/CONTCAR POSCAR
cp INCAR.init INCAR
Note. If you have not done part2.ipynb, you will not be able to copy the structure from cp $TUTORIALS/transition_states/e09_*/CONTCAR. In this case, use the structure from POSCAR.init used here to POSCAR.
Click here if you have not done part2!
cd e10_*
cp POSCAR.init POSCAR
The directory is now prepared to run a molecular dynamics (MD) simulation of the reaction. The following run script is prepared to run 50000 steps across five separate runs.
Click here to reveal the run script!
#!/usr/bin/bash
runvasp="mpirun -np 4 vasp_gam"
#c initialize RNG
aaa=$(echo "RANDOM_SEED = 817605623 0 0" )
step=1
while [ $step -le 5 ]
do
cp POSCAR POSCAR.$step
cp INCAR.init INCAR
echo $aaa >> INCAR
$runvasp
#c extract seed of RNG from the last step
aaa=$(grep RANDOM_SEED REPORT|tail -1 )
#c save important output files
cp REPORT report.$step
cp vasprun.xml vasprun.xml.$step
cp CONTCAR POSCAR
let step+=1
done
This calculation will take around 8-10 minutes and can be run with the following:
bash run
If you are short on time, you may perform a single MD run instead and then use the reference files provided (make sure to copy POSCAR.init to POSCAR) for the remaining MD runs:
Click here for a script that will perform only a single MD run!
#!/usr/bin/bash
cp INCAR.init INCAR
#c initialize RNG
aaa=$(echo "RANDOM_SEED = 817605623 0 0" )
echo $aaa >> INCAR
mpirun -np 4 vasp_gam
cp CONTCAR POSCAR.1
cp REPORT report.1
cp vasprun.xml vasprun.xml.1
cd refdata
cp POSCAR.2 POSCAR.3 POSCAR.4 POSCAR.5 ../
cp report.2 report.3 report.4 report.5 ../
cp vasprun.xml.2 vasprun.xml.3 vasprun.xml.4 vasprun.xml.5 ../
cd ../
Note, that if the calculation is taking a long time, you can analyze the reference data found in ./e10_*/refdata while it is running before repeating the analysis with your own results. Once your MD simulation has finished, take a look at the directory ./e10_*. There are five files report.1 to report.5 that have been copied from the REPORT output from VASP. REPORT files begin with a description of the initial MD conditions:
MDALGO = 5
CSVR_PERIOD = 20.000
SCALING = 0
CNEXP = 9.000 14.000
RANDOM_SEED = 817605623 0 0
original number of atomic DOF: 123
number of independent constraints: 0
number of active DOF: 123
And then the MD output every 10 steps, which concerns us more.
========================================
MD step No. 10
========================================
>Monit_coord
mc> S 1.94477
>Energies
E_tot E_pot E_kin EPS ES
e_b> -0.31931873E+03 -0.31977016E+03 0.45143057E+00 0.00000000E+00 0.00000000E+00
>Temperature
T_sim T_inst
tmprt> 440.000 85.181
RANDOM_SEED = 817605623 2460 0
The meaning of each term of the output is summarised in the following table:
| tag | meaning |
|---|---|
MDALGO |
Sets the thermostat for MD, here the thermostat is CSVR, hence MDALGO = 5. |
CSVR_PERIOD |
Defines the time scale of the CSVR thermostat in terms of the number of MD steps |
SCALING |
Relates to the scaling the velocities. It is not used in this calculation. |
CNEXP |
Are settings for the ICONST flag = D. It is not used in this calculation |
RANDOM_SEED |
A random number generator that defines the velocities. |
DOF |
Degrees of freedom of the ions. |
mc> |
This is the monitored reaction coordinate, defined in the introduction for this exercise. |
e_b> |
MD energies: E_tot is the total energy, E_pot is the potential energy, E_kin is the kinetic energy, EPS is the potential energy of the Nosé-mass (not used here), and ES is the kinetic energy of the Nosé-mass (not used here). |
tmprt> |
T_sim is the initial and target simulation temperature, while T_inst is the temperature of the current MD step. |
We are interested in the monitored reaction coordinate mc> in the report files. The script monitorVal.sh extracts mc> for us. Tracking mc> will allow us to follow the reaction coordinate $\xi$ throughout the MD simulation.
cd $TUTORIALS/transition_states/e10_*
bash ../scripts/bashScripts/monitorVal.sh
The script monitorVal.sh outputs an mc.dat file, which contains all $\xi$ from the MD simulation. The probability density distribution of these $\xi$ can be plotted in a histogram. The script normalizedHistogram.py reads mc.dat and outputs the probability density distribution for 100 intervals to the mc.dat.hist file.
cd $TUTORIALS/transition_states/e10_*
python3 ../scripts/normalizedHistogram/normalizedHistogram.py mc.dat 100
Plot the histogram stored in the mc.dat.hist file using the following script:
import plotly.graph_objects as go
import pandas as pd
# Load data from mc.dat.hist
#df1 = pd.read_csv("./e10_Dynamic_Coordinates/refdata/mc.dat.hist", sep="\t", header=None) # reference data
df2 = pd.read_csv("./e10_Dynamic_Coordinates/mc.dat.hist", sep="\t", header=None) # uncomment to plot your own data
# Find maximum of the probability density
calc_max = df2[0][df2[1] == df2[1].max()].values[0] # Get the x-coordinate of the max
static = 1.95 # Static reaction coordinate for reactant
# Create the histogram
histogram = go.Bar(
x=df2[0],
y=df2[1],
marker=dict(color='#4c265f', line=dict(color='#4c265f', width=1)),
name='calc.'
)
# Add vertical lines for the maxima
dynamic_line = go.Scatter(
x=[calc_max, calc_max],
y=[0, df2[1].max()],
mode='lines',
line=dict(color='#2fb5ab', dash='dot'),
name="\u03BE<sub>ref,R</sub><sup>dynamic</sup>",
showlegend=False,
)
static_line = go.Scatter(
x=[static, static],
y=[0, df2[1].max()],
mode='lines',
line=dict(color='#2fb5ab', dash='dot'),
name="\u03BE<sub>ref,R</sub><sup>static</sup>",
showlegend=False,
)
# Add annotations
annotations = [
dict(
x=calc_max - 0.02,
y=df2[1].max() + 0.03,
text="\u03BE<sub>ref,R</sub><sup>dynamic</sup>",
showarrow=False,
),
dict(
x=static - 0.07,
y=df2[1].max(),
text="\u03BE<sub>ref,R</sub><sup>static</sup>",
showarrow=False,
),
]
# Build the layout
layout = go.Layout(
#title="Histogram of Reaction Coordinate",
xaxis=dict(title="\u03BE"+" in Å"),
yaxis=dict(title=r"Probability density of "+"\u03BE"+" in Å<sup>-1</sup>", range=[0.0, 1.2], dtick = 0.2),
annotations=annotations,
showlegend=True,
font_color="black",
template="ggplot2",
)
# Combine the plots
fig = go.Figure(data=[histogram, dynamic_line, static_line], layout=layout)
# Show the plot
fig.show()
The area of highest probability density is where we will take our trial reactant structure from. The histogram produced by normalizedHistogram.py is somewhat jagged but provides a clear image of the distribution. To more precisely locate the reactant state, the script normalizedKDE.py is used. It smooths out the histogram using the kernel density estimation (KDE) with the following script, in this case smoothed to 0.05 Å. The smoothing and number of intervals (100 here) should be selected for the individual problem:
cd $TUTORIALS/transition_states/e10_*
python3 ../scripts/normalizedKDE/normalizedKDE.py mc.dat 0.05 100
import plotly.graph_objects as go
import pandas as pd
# Load data from mc.dat.hist
#df1 = pd.read_csv("./e10_Dynamic_Coordinates/refdata/mc.dat.kde", sep="\t", header=None) # reference data
df2 = pd.read_csv("./e10_Dynamic_Coordinates/mc.dat.kde", sep="\t", header=None) # uncomment to plot your own data
# Find maximum of the probability density
calc_max = df2[0][df2[1] == df2[1].max()].values[0] # Get the x-coordinate of the max
static = 1.95 # Static reaction coordinate for reactant
# Create the histogram
histogram = go.Bar(
x=df2[0],
y=df2[1],
marker=dict(color='#4c265f', line=dict(color='#4c265f', width=1)),
name='calc.'
)
# Add vertical lines for the maxima
dynamic_line = go.Scatter(
x=[calc_max, calc_max],
y=[0, df2[1].max()],
mode='lines',
line=dict(color='#2fb5ab', dash='dot'),
name="\u03BE<sub>ref,R</sub><sup>dynamic</sup>",
showlegend=False,
)
static_line = go.Scatter(
x=[static, static],
y=[0, df2[1].max()],
mode='lines',
line=dict(color='#2fb5ab', dash='dot'),
name="\u03BE<sub>ref,R</sub><sup>static</sup>",
showlegend=False,
)
# Add annotations
annotations = [
dict(
x=calc_max,
y=df2[1].max() + 0.07,
text="\u03BE<sub>ref,R</sub><sup>dynamic</sup>",
showarrow=False,
),
dict(
x=static - 0.07,
y=df2[1].max(),
text="\u03BE<sub>ref,R</sub><sup>static</sup>",
showarrow=False,
),
]
# Build the layout
layout = go.Layout(
#title="Histogram of Reaction Coordinate",
xaxis=dict(title="\u03BE"+" in Å"),
yaxis=dict(title=r"Probability density of "+"\u03BE"+" in Å<sup>-1</sup>", range=[0.0, 1.2], dtick = 0.2),
annotations=annotations,
showlegend=True,
font_color="black",
template="ggplot2",
)
# Combine the plots
fig = go.Figure(data=[histogram, dynamic_line, static_line], layout=layout)
# Show the plot
fig.show()
The area around the maximum of the probability density $P(\xi)$ ≈ 1.05 Å $^{-1}$ (reference data) is the best-sampled region, $\xi$ ≈ 2.4 Å, that is the reactant state. We select the reference reactant state from the best-sampled region to keep the relative error in $P(\xi)$ values as low as possible.
Note, the reactant state initial POSCAR.1 was taken from the static approach, so $\xi_{ref,R}$ can be taken from the value of mc> S given at the start of the first MD run, i.e. 1.95 Å (taken from the CONTCAR in Example 9). The change in $\xi_{ref,R}$ going from static to dynamic approaches is 1.95 Å compared to 2.40 Å. This is a significant difference of 0.45 Å, illustrating the difference between the two approaches. In the static approach, a single structure is used to describe the reactant state, whereas in MD, the reactant ensemble is described by a distribution of many different structures.
Now that we have chosen a $\xi_{ref,R}$ value from the probability distribution, we can check the $\xi$ for each of the POSCAR.* files to see which lies the closest to our desired value. The structure of the stored POSCAR files can be found in $TUTORIALS/transition_states/e10_*/POSCAR.*. Run the following script to find $\xi$ for the five steps:
cd $TUTORIALS/transition_states
for i in {1..5}; do echo report.$i "POSCAR.$((i+1))" `grep mc e10_*/report.$i | tail -1 | awk '{print $3}'`; done
echo "Note: POSCAR.6 is the final CONTCAR"
%%bash
for i in {1..5}; do echo report.$i "POSCAR.$((i+1))" `grep mc e10_*/report.$i | tail -1 | awk '{print $3}'`; done
echo "Note: POSCAR.6 is the final CONTCAR"
Note that POSCAR.6 is the final CONTCAR.
With the reference data (e10_*/refdata), we find that the structure that lies closest to this maximum is the CONTCAR of the first run, i.e. POSCAR.2, with $\xi$ = ~2.32 Å and $P(\xi)$ = ~0.87 Å $^{-1}$), which will be used in a later exercise. Take care when selecting a POSCAR to continue investigating the reaction with, it may not be POSCAR.2 for you and you will have to edit the following scripts to account for this discrepancy.
Congratulations, you have now found a good approximation of the reactant state using MD! You will use this structure in the next exercise to test the reversibility of the reaction coordinate.
10.4 Questions¶
- What structure does the probability peak coincide with?
- How do the static and dynamic reaction coordinates differ?
By the end of this tutorial, you will be able to:
- transform reactant (R) to transition state (TS) using constrained MD
- test for the continuity of the free energy gradient
- (optional) test the reversibility of the process from TS to R
11.1 Task¶
Scan the free energy gradient along the reaction coordinate from reactant to transition state during the carbonylation of methylated chabazite using a Blue moon ensemble in molecular dynamics (MD) for 24 ps with a CSVR thermostat at 440 K
The suitability of our reaction coordinate $\xi$ can be tested by following the change in free energy gradient $({\partial A} / {\partial \xi_k})_{\xi^*}$ from reactant R to transition state TS, where $\xi_{ref,R}$ and $\xi^*$ are R and TS, respectively. Bear in mind that R and TS are ensembles in MD, not individual structures. The free energy gradient determines the free energy of activation. Tracking it using constrained molecular dynamics generates biased statistical averages.
One suitable approach is using a Blue moon ensemble, which yields the average for a quantity, in our case the free energy gradient. The Blue moon ensemble is relatively simple to use and allows for easy statistical error estimation and to diagnose problems relating to the choice of $\xi$. Using the slow-growth approach, the free energy gradient will be scanned along $\xi$ from reactant to product with a velocity of transformation $\dot{\xi}$, i.e. the change of $\xi$ per step; TS will lie at the maximum of this free energy plot.
A sufficiently small value of $\dot{\xi}$ must be chosen such that the change in free energy gradient is continuous from R to TS. $\dot{\xi}$ is calculated based on the distance to be covered from R to TS, $\xi_{ref,R} \rightarrow \xi^{*}$, dividing this distance by the number of MD steps to yield $\dot{\xi}$. If the free energy is continuous from R to TS, then our choice of $\xi$ is affirmed. A final check of the quality of $\xi$ is that this process is reversible, which requires the procedure to be followed in reverse from TS to R.
In this exercise, you will use constrained MD to scan the free energy gradient along the reaction coordinate from reactant to transition state. If you have sufficient time, the reversibility of the process may be checked by repeating the process in reverse.
11.2 Input¶
The input files to run this example are prepared at $TUTORIALS/transition_states/e11_Dynamic_Validation.
POSCAR.init
Click here to see the initial POSCAR!
Initial carbonylation step ch3-chabazite system
1.00000000000000
7.9487550000000002 -0.0000500000000000 4.9014329999999999
-3.9745640000000000 6.8838939999999997 4.9012549999999999
-3.9744389999999998 -6.8837770000000003 4.9012460000000004
Al C H O Si
1 2 3 25 11
Direct
0.6387710208557332 0.3907175316431041 0.3555272802644143
-0.2437288530469274 0.7221304421154419 0.3394681113161917
0.6196484383307015 0.0888195592882838 0.3273090915859819
0.6667162005136688 0.0480052933182199 0.2239622317982071
0.7103459923079960 0.1127555039761568 0.4079258429185892
0.5403796471339859 0.0017509353146397 0.3553646828991581
0.2772993598473935 0.7637916742247260 -0.0092738633486380
0.1953472143574169 0.4941125387847758 0.9878722142753384
0.2621072217372572 0.5955468900036330 0.7332383630127673
0.5072025566008931 0.7676278687110434 0.7630114506520603
0.0310288745319219 0.3317903704331803 0.1424800858866574
0.3847280079993154 0.6642928115334695 0.5052103910843849
0.2197101874021943 0.5276326368298239 0.2829052903144220
0.4012034667970906 0.7590710552794698 0.2551900470908230
0.5176850915224630 0.4884784250859319 0.6871450003805950
0.7647057523945355 0.7509277495235895 0.8972445522475090
-0.1946088376521494 0.7194819687405966 0.4580200998367944
0.6053902989350126 0.2557840185206084 0.7520655206925642
0.5557278262274221 0.8431519511563210 0.0544658968911583
0.2947920558499533 0.2805256011050637 0.1534962875026886
0.8015134935631735 0.4881764570963363 0.7688312566132961
0.5054372664631622 0.5001105777211905 0.3163255653649657
0.7957093163806136 0.2617859831939608 0.9637430623123721
0.7047569860995276 0.3673222911461803 0.5204833457937326
0.5337095128642886 0.2052350526477613 0.2986282763465009
0.6719746223655358 1.0013390719871922 0.8601720365932347
0.0190721421869335 0.6633304364822351 0.8626289308777665
0.3589445244927258 1.0190927540477548 0.1350708626174435
0.5116627137462921 0.1854456507795967 -0.0021364865223991
0.7715275318430667 0.3880157255258864 0.2406494588495530
0.8569116535063408 0.5402990545881321 0.0466398885327390
0.1803524848515072 0.6255251907351898 0.8866067657003128
0.4245244908746708 0.6299234254334347 0.6746745495256028
0.6254684977756348 0.8399816209395478 0.8942971595143633
0.3927740308938757 0.8484870965385048 0.1150995539269583
0.3820828018169702 0.6135923749031230 0.3390490873837098
0.1872076357444247 0.4098487369655657 0.1368874697513225
0.6651400968443657 0.4014383452741130 0.6809322276452606
0.8548762829843090 0.6101109814018703 0.8971139196993789
0.6470841373481644 0.1704202098270244 0.8857821998710125
0.4220184547303953 0.1801607968001600 0.1329544255406982
0.8672256524014095 0.3801671038798720 0.0992538627633380
INCAR.init
! Molecular dynamics
IBRION = 0 ! Molecular dynamics
ISYM = 0 ! Symmetry setting
NPAR = 2 ! Number of bands treated in parallel
NSW = 4000 ! Number of steps in MD run
POTIM = 1.0 ! Time step
MDALGO = 5 ! CSVR thermostat
TEBEG = 440. ! Temperature
CSVR_PERIOD = 20 ! Time scale of CSVR in terms of the number of MD steps
! Machine learning
ML_LMLFF = T ! Enables the use of MLFF
ML_MODE = run ! Prediction-only mode
ML_OUTPUT_MODE = 0 ! Reduced file output mode
ML_OUTBLOCK = 10 ! Output distance in number of steps
INCREM = -1.1e-4 ! Slow-growth simulation protocol
LBLUEOUT = T ! Blue moon ensemble
K-Points
0
Gamma
1 1 1
0 0 0
Pseudopotentials of Al, C, H, O, and Si.
R 3 25 0
R 3 2 0
S -1 1 7
The POSCAR, POTCAR, KPOINTS, and ICONST files are identical to in Example 10. The INCAR file is largely the same, with the exception of two new tags: INCREM and LBLUEOUT. INCREM controls the transformation velocity using a slow-growth simulation. A negative value is chosen so that the reaction coordinate $\xi$ is followed from reactant to transition state. LBLUEOUT = T performs a constrained MD calculation, specifically a Blue moon ensemble, which writes the free energy gradient to REPORT.
11.3 Calculation¶
Select an initial POSCAR file for the MD simulation using the following script:
cd $TUTORIALS/transition_states
for i in {1..5}; do echo report.$i "POSCAR.$((i+1))" `grep mc e10_*/report.$i | tail -1 | awk '{print $3}'`; done
echo "Note: POSCAR.6 is the final CONTCAR"
If you are having difficulty, copy the POSCAR.init structure to POSCAR and start from there:
cd $TUTORIALS/transition_states/e11_*
cp POSCAR.init POSCAR
ln -s $TUTORIALS/transition_states/mlff/MLFF_e05-12/ML_FFN ML_FF
Then skip the next three lines (cd, ln, and cp).
If things are working as expected, open a terminal and navigate to this example's directory. Note: replace POSCAR.2 with the POSCAR.* that is closest to the reactant state's reaction coordinate $\xi$.
cd $TUTORIALS/transition_states/e11_*
ln -s ../mlff/MLFF_e05-12/ML_FFN ML_FF
cp $TUTORIALS/transition_states/e10_*/POSCAR.2 POSCAR
The Blue moon simulation models the free energy gradient as the reaction coordinate $\xi$ changes from reactant to transition state. We are checking the continuity of the free energy gradient here. The MD simulation is split into 6 runs of 4000 steps each (cf. INCAR and run). Running this MD simulation will provide initial structures for a regular grid of 7 reaction coordinates $\xi$, which are the POSCAR files of each run: POSCAR.1 to POSCAR.6 and the final CONTCAR file. Once the simulation is finished, the structural changes from the change in $\xi$ should be checked to ensure that it is indeed the expected mechanism. Run the MD simulation with the following script; the calculation will take approximately five minutes:
Click here to reveal the run script!
#!/usr/bin/bash
runvasp="mpirun -np 4 vasp_gam"
#c initialize RNG
aaa=$(echo "RANDOM_SEED = 817605623 0 0" )
step=1
while [ $step -le 6 ]
do
cp POSCAR POSCAR.$step
cp INCAR.init INCAR
echo $aaa >> INCAR
$runvasp
#c extract seed of RNG from the last step
aaa=$(grep RANDOM_SEED REPORT|tail -1 )
#c save important output files
cp REPORT report.$step
cp vasprun.xml vasprun.xml.$step
cp CONTCAR POSCAR
let step+=1
done
bash run
Note! If you see Error: SHAKE algorithm did not converge! err= in your stdout, there is a problem. The POSCAR.* file you have chosen is too far away from the reactant state - it is a poor starting choice. Try reselecting the initial POSCAR so that $\xi$ is closer to the maximum of the probability density distribution in the previous exercise.
Congratulations on completing a Blue moon ensemble calculation!
You have produced seven structures (POSCAR.1 to POSCAR.6 and the final CONTCAR file) tracing $\xi$ evenly throughout the reaction. The script pos2xyzStatic.py takes these seven structure files and puts them together to create an animation of the reaction, similar to in Example 8, using the following:
python3 ../scripts/pos2xyzNVT/pos2xyzStatic.py POSCAR.1 0 0.2 0 > anim.xyz
python3 ../scripts/pos2xyzNVT/pos2xyzStatic.py POSCAR.2 0 0.2 0 >> anim.xyz
python3 ../scripts/pos2xyzNVT/pos2xyzStatic.py POSCAR.3 0 0.2 0 >> anim.xyz
python3 ../scripts/pos2xyzNVT/pos2xyzStatic.py POSCAR.4 0 0.2 0 >> anim.xyz
python3 ../scripts/pos2xyzNVT/pos2xyzStatic.py POSCAR.5 0 0.2 0 >> anim.xyz
python3 ../scripts/pos2xyzNVT/pos2xyzStatic.py POSCAR.6 0 0.2 0 >> anim.xyz
python3 ../scripts/pos2xyzNVT/pos2xyzStatic.py CONTCAR 0 0.2 0 >> anim.xyz
pos2xyzStatic.py converts the input structure to xyz coordinates, which is written to the anim.xyz file. You can visualize it in JMol once downloaded. A recording of anim.xyz is provided below:
Click here to watch a video of the reaction from reactant to transition state to product!
The animation of anim.xyz is brief, consisting of only 7 images. In the first image, the methyl CH3 group is still bonded to a zeolite-O atom (ZO) and the CO molecule is elsewhere in the cell. From the second step onwards, the CO and CH3 come closer together, with the ZO-CH3 bond breaking and the CH3 inverting until the CH3 and CO form a bond to create an acetyl cation H3C-CO+ in the final image (CONTCAR). The progression from reactant to product via the transition state is as we expected it to be. Note that the 7 images are individual structures taken from their own ensemble; they are representatives of each ensemble which contains many structures.
Having confirmed that the transformation is as expected, we can confirm the suitability of our guess of $\xi$ using the estimated free energy gradients. The $\xi$ and free energy gradient for each of the MD calculations are found in the report.* files, cc> and b_m>. report.* files are copied from the REPORT file from each calculation. They contain information about the MD ensemble every 10 steps.
Click here to take a look at an MD step from the report.1 file!
========================================
MD step No. 10
========================================
>Energies
E_tot E_pot E_kin EPS ES
e_b> -0.31593178E+03 -0.31851585E+03 0.25840637E+01 0.00000000E+00 0.00000000E+00
>Temperature
T_sim T_inst
tmprt> 440.000 491.584
SHAKE converged in 8 steps
>Const_coord
cc> S 2.60744 2.60744 0.423625E-05
>Blue_moon
lambda |z|^(-1/2) GkT |z|^(-1/2)*(lambda+GkT)
b_m> -0.191337E+00 0.185300E+01 -0.125565E-02 -0.356874E+00
RANDOM_SEED = 817605623 2440 0
Note, the SHAKE algorithm is used to calculate the blue moon ensemble.
The script fgrad.sh extracts $\xi$ and the free energy gradient from the report.* files and outputs grad.dat using the following:
../scripts/bashScripts/fgrad.sh
The resulting grad.dat file can be visualized using the following script:
import pandas as pd
import plotly.graph_objects as go
from scipy.interpolate import Rbf
# Load the free energy gradient data
df = pd.read_csv("./e11_Dynamic_Validation/grad.dat", sep="\t", header=None)
# Extract x and y values
x, y = df[0], df[1]
# Compute the moving average using a radial basis function (RBF)
rbf = Rbf(x, y, function='multiquadric', smooth=1000) # RBF interpolation
y_rbf = rbf(x) # Smoothed y values
# Create the plotly traces
raw_data_trace = go.Scatter(
x=x,
y=y,
mode='lines',
line=dict(color='#2c68fc'),
name='forward',
)
smoothed_trace = go.Scatter(
x=x,
y=y_rbf,
mode='lines',
line=dict(color='#4c265f', dash='solid'),
name='forward avg.',
)
# Create the layout
layout = go.Layout(
xaxis=dict(title="ξ in Å"),
yaxis=dict(title="Free energy gradient in eV Å⁻¹"),
legend=dict(x=0.05, y=0.95), # Position the legend
font_color="black",
template="ggplot2",
)
# Combine traces into a figure
fig = go.Figure(data=[raw_data_trace, smoothed_trace], layout=layout)
# Show the plot
fig.show()
It is clear that there are no obvious discontinuities in the free energy gradient. The continuity of the free energy gradient during the transformation validates our choice of $\xi$.
Alongside the expected reaction being followed, the continuity is the second validation of the $\xi$ that we have chosen.
(OPTIONAL) A final test for the suitability of $\xi$ is to repeat the process in reverse. This calculation will take approximately five minutes. The reverse transformation may be done by starting from the final TS structure generated already in this exercise and changing the sign of INCREM from -1.1e-4 to 1.1e-4 in the INCAR file. A requirement of $\xi$ is that it drives the process reversibly. Reversibility means that the free energy gradient profile should be qualitatively similar in reverse as forward. The reverse slow growth transformation is prepared in ./optional:
cd $TUTORIALS/transition_states/e11_*
cp POTCAR KPOINTS ICONST run ./optional
cp CONTCAR ./optional/POSCAR
cp INCAR.init ./optional/
cd ./optional
ln -s ../../mlff/MLFF_e05-12/ML_FFN ML_FF
sed -i 's/INCREM = -1.1e-4/INCREM = 1.1e-4/g' INCAR.init
And then run in the same directory:
bash run
Once this calculation has finished, check the similarity of free energy gradient with that obtained the first part using fgrad.sh.
../../scripts/bashScripts/fgrad.sh
import pandas as pd
import plotly.graph_objects as go
from scipy.interpolate import Rbf
# Load the free energy gradient data
df1 = pd.read_csv("./e11_Dynamic_Validation/grad.dat", sep="\t", header=None)
df2 = pd.read_csv("./e11_Dynamic_Validation/optional/grad.dat", sep="\t", header=None)
# Extract x and y values
x, y = df[0], df[1]
# Generates a moving average
x1, y1 = df1[0], df1[1]
rbf1 = Rbf(x1, y1, function = 'multiquadric', smooth = 1000) # df1 - calculates a moving average of the free energy gradient using a radial basis function (RBF)
y_rbf1 = rbf1(x1) # df1 - y-coordinates for RBF
x2, y2 = df2[0], df2[1]
rbf2 = Rbf(x2, y2, function = 'multiquadric', smooth = 1000) # df2 - calculates a moving average of the free energy gradient using a radial basis function (RBF)
y_rbf2 = rbf2(x2) # df1 - y-coordinates for RBF
# Create the plotly traces
smoothed_forward = go.Scatter(
x=x1,
y=y_rbf1,
mode='lines',
line=dict(color='#4c265f', dash='solid'),
name='forward avg.',
)
smoothed_reverse = go.Scatter(
x=x2,
y=y_rbf2,
mode='lines',
line=dict(color='#a82c35', dash='solid'),
name='forward avg.',
)
raw_data_forward = go.Scatter(
x=x1,
y=y1,
mode='lines',
line=dict(color='#2c68fc'),
name='forward',
opacity=0.5,
)
raw_data_reverse = go.Scatter(
x=x2,
y=y2,
mode='lines',
line=dict(color='#2fb5ab'),
name='forward',
opacity=0.5,
)
# Create the layout
layout = go.Layout(
xaxis=dict(title="ξ in Å"),
yaxis=dict(title="Free energy gradient in eV Å⁻¹"),
legend=dict(x=0.05, y=0.95), # Position the legend
font_color="black",
template="ggplot2",
)
# Combine traces into a figure
fig = go.Figure(data=[raw_data_forward, smoothed_forward, raw_data_reverse, smoothed_reverse], layout=layout)
# Show the plot
fig.show()
The average free energy gradient changes similarly for both the forward and backward directions, meaning that the process is reversible.
You have now validated your choice of reaction coordinate $\xi$ in two (or three with the optional exercise) times. In the next exercise, you will determine the free energy barrier from reactant to transition state and the corresponding free energy of activation.
11.4 Questions¶
- What tag sets the step size in the slow-growth approximation?
- How should the free energy gradient change from the reactant to the transition state?
- How similar should the free energy gradient be in the forward and reverse directions?
By the end of this tutorial, you will be able to:
- find the free energy gradient along the reaction coordinate from R to TS
- derive the free energy of the reaction
- calculate the rate constant within the dynamic approximation
- compare results using static and dynamic approaches
12.1 Task¶
Calculate the free energy gradient using a Blue moon simulation for seven structures along the the reaction coordinate of the carbonylation of methylated chabazite using constrained molecular dynamics (MD) for 140 ps with a CSVR thermostat at 440 K
The transformation from reactant R to transition state TS produces work in the form of the free energy of activation $\Delta A_{R \rightarrow TS}^{\ddagger}$. The transformation is propagated in a sequence of constrained molecular dynamics (MD) runs, each starting from a different reaction coordinate $\xi$. As done previously, a Blue moon ensemble is used. The free energy of activation is then calculated:
$\normalsize{\Delta A_{R \rightarrow TS}^{\ddagger}} = \normalsize{{\Delta A_{\xi_{ref,R} \rightarrow \xi^{*}}}} - k_BT \ln[\Large{\frac{h}{k_bT} \frac{\langle | \dot\xi^* | \rangle}{2}} \normalsize{P(\xi_{ref,R})}] $
where $\Delta A_{\xi_{ref,R} \rightarrow \xi^{*}}$ is the reversible work done to shift the reaction coordinate from $\xi_{ref,R}$ (the value of $\xi$ for a reference point among the R configurations) to $\xi^*$ (the value of $\xi$ for TS configurations), $\langle ... \rangle$ denotes an ensemble average over all reactant configurations, $\dot\xi^*$ is the velocity associated with the reaction coordinate for configurations at the TS, and $P$ is the probability density within the ensemble of reactant configurations.
The reversible work done is calculated with the following:
$ \normalsize {\Delta A_{\xi_{ref,R} \rightarrow \xi^{*}}} = \Large {\int_{\xi_{ref,R}}^{\xi^{*}} {(\frac{\partial A}{\partial \xi})}_{\xi^*} \,} \normalsize {d\xi} $
where $({\partial A} / {\partial \xi})_{\xi^*}$ is the free energy gradient, defined using the following ensemble average:
$ \Large{(\frac{\partial A}{\partial \xi})}_{\xi^*} = \frac{1}{\langle|Z|^{-1/2}\rangle_{\xi^*}} \langle \large{|Z|^{-1/2} [\lambda + \textit{G} k_{B}T ]} \rangle _{\xi^*} $
where $\langle ... \rangle_{\xi^*}$ denotes a constrained ensemble average where $\xi = \xi^* $ , $Z$ is a mass metric tensor, $\lambda$ are the Lagrange multipliers associated with $\xi$ used in the SHAKE equation, and $\textit{G}$ is derived from $|Z|$, $Z^{-1}$, $\nabla\xi$, and $\nabla|Z|$; further details are available elsewhere and references therein.
Using the free energy gradient, the Eyring equation describes a rate constant $k$:
$ k_{R\rightarrow P} = \Large{\frac{k_{B}T}{h}e^{-\Delta A_{R\rightarrow P}^{\ddagger} / k_B T}} $
In this exercise, you will calculate the free energy gradient along $\xi$ between R and TS using constrained MD. You will calculate the free energy of reaction and the rate constant for the carbonylation of methylated-chabazite within a dynamic approximation.
12.2 Input¶
The input files to run this example are prepared at $TUTORIALS/transition_states/e12_Dynamic_Thermodynamics.
Click here to reveal the seven initial input POSCAR files!
POSCAR.1_init
Reaction pathway test
1.00000000000000
7.9487550000000002 -0.0000500000000000 4.9014329999999999
-3.9745640000000000 6.8838939999999997 4.9012549999999999
-3.9744389999999998 -6.8837770000000003 4.9012460000000004
Al C H O Si
1 2 3 25 11
Direct
0.6387710208557332 0.3907175316431041 0.3555272802644143
-0.2437288530469274 0.7221304421154419 0.3394681113161917
0.6196484383307015 0.0888195592882838 0.3273090915859819
0.6667162005136688 0.0480052933182199 0.2239622317982071
0.7103459923079960 0.1127555039761568 0.4079258429185892
0.5403796471339859 0.0017509353146397 0.3553646828991581
0.2772993598473935 0.7637916742247260 -0.0092738633486380
0.1953472143574169 0.4941125387847758 0.9878722142753384
0.2621072217372572 0.5955468900036330 0.7332383630127673
0.5072025566008931 0.7676278687110434 0.7630114506520603
0.0310288745319219 0.3317903704331803 0.1424800858866574
0.3847280079993154 0.6642928115334695 0.5052103910843849
0.2197101874021943 0.5276326368298239 0.2829052903144220
0.4012034667970906 0.7590710552794698 0.2551900470908230
0.5176850915224630 0.4884784250859319 0.6871450003805950
0.7647057523945355 0.7509277495235895 0.8972445522475090
-0.1946088376521494 0.7194819687405966 0.4580200998367944
0.6053902989350126 0.2557840185206084 0.7520655206925642
0.5557278262274221 0.8431519511563210 0.0544658968911583
0.2947920558499533 0.2805256011050637 0.1534962875026886
0.8015134935631735 0.4881764570963363 0.7688312566132961
0.5054372664631622 0.5001105777211905 0.3163255653649657
0.7957093163806136 0.2617859831939608 0.9637430623123721
0.7047569860995276 0.3673222911461803 0.5204833457937326
0.5337095128642886 0.2052350526477613 0.2986282763465009
0.6719746223655358 1.0013390719871922 0.8601720365932347
0.0190721421869335 0.6633304364822351 0.8626289308777665
0.3589445244927258 1.0190927540477548 0.1350708626174435
0.5116627137462921 0.1854456507795967 -0.0021364865223991
0.7715275318430667 0.3880157255258864 0.2406494588495530
0.8569116535063408 0.5402990545881321 0.0466398885327390
0.1803524848515072 0.6255251907351898 0.8866067657003128
0.4245244908746708 0.6299234254334347 0.6746745495256028
0.6254684977756348 0.8399816209395478 0.8942971595143633
0.3927740308938757 0.8484870965385048 0.1150995539269583
0.3820828018169702 0.6135923749031230 0.3390490873837098
0.1872076357444247 0.4098487369655657 0.1368874697513225
0.6651400968443657 0.4014383452741130 0.6809322276452606
0.8548762829843090 0.6101109814018703 0.8971139196993789
0.6470841373481644 0.1704202098270244 0.8857821998710125
0.4220184547303953 0.1801607968001600 0.1329544255406982
0.8672256524014095 0.3801671038798720 0.0992538627633380
POSCAR.2_init
Reaction pathway test
1.00000000000000
7.9487550000000002 -0.0000500000000000 4.9014329999999999
-3.9745640000000000 6.8838939999999997 4.9012549999999999
-3.9744389999999998 -6.8837770000000003 4.9012460000000004
Al C H O Si
1 2 3 25 11
Direct
0.6384233230045857 0.3635179895068407 0.3503296690868836
0.7664241800249388 0.9941963253978896 0.2738126364973068
0.4730213225941887 0.1153322424537904 0.4262223808573459
0.3796645950189537 0.1536992785826949 0.4676840642614477
0.4499009943123430 -0.0044574742144416 0.3945323390042894
0.5676498319127665 0.1341265853278163 0.4979090952714264
0.2702507271254301 0.7681438592712677 0.9763816536165626
0.2230289256586666 0.4844959869796193 0.9965200746499459
0.2725631580203081 0.5941595244639706 0.7433020008485752
0.5180789534142199 0.7619741281633250 0.7584733015646742
0.0362898170130195 0.3663527092369872 0.1720999140078053
0.3785768604756903 0.6660253451248084 0.4997725702914820
0.2711245156141834 0.4756878486628307 0.2909610002859832
0.3890134793764573 0.7496709993088647 0.2513388594613383
0.4896002668583118 0.4594818647034135 0.6536952942971421
0.7876708901316891 0.7594389235716750 0.9131194474365465
0.7578224836550844 1.0675266847288403 0.3754269666820085
0.6126242122541289 0.2371295961869578 0.7613985005702495
0.5512007813947946 0.8093332473575384 0.0412050483127940
0.2854536512556322 0.2478648912536432 0.1134266339905133
0.7623522889013794 0.4985597540866310 0.7608591244186536
0.5513316618760031 0.5238445391108352 0.3437299548056363
0.8286981317339803 0.2463080627870407 0.9729130331295007
0.6942041661092634 0.3239371629299048 0.5233344168591992
0.5185310168666739 0.1940941188832715 0.2963652305201390
0.6804057840962242 -0.0009585910650768 0.8865902916488048
0.0270551393905135 0.6340015634659456 0.8464937975291433
0.3767030537361343 -0.0125456377016687 0.1387127209798370
0.5359580092262458 0.1875349140681897 1.0122344893043556
0.7675993000082955 0.3536314361269847 0.2327594253164944
0.8710301055681686 0.5369349657269955 0.0457718778299470
0.1937153636589095 0.6143580764735017 0.8915386099387406
0.4165569803083761 0.6174369396604461 0.6612633607054685
0.6363013555727809 0.8276841327428833 0.8967614552127811
0.3907700160229702 0.8210455075376615 0.1023154819762266
0.4028141731668105 0.5996622533344264 0.3384335113293067
0.2016220881127892 0.3894394754553366 0.1410873996556337
0.6366196687757727 0.3819379384428966 0.6722106163754520
0.8628462002387601 0.6065189901076877 0.8856551864897692
0.6638254622184195 0.1671302177460708 0.9056847788294697
0.4275706653554721 0.1564001893374446 0.1331146116929661
0.8713447072231915 0.3735442027299686 0.1038227844078963
POSCAR.3_init
Reaction pathway test
1.00000000000000
7.9487550000000002 -0.0000500000000000 4.9014329999999999
-3.9745640000000000 6.8838939999999997 4.9012549999999999
-3.9744389999999998 -6.8837770000000003 4.9012460000000004
Al C H O Si
1 2 3 25 11
Direct
0.6593679998853460 0.4106543551089780 0.3497985170603597
0.8205228583083806 0.8713796434489578 0.3102997793591317
0.6369366971443284 0.1041174599612637 0.3049838023181419
0.7348069459739432 0.1548837740686913 0.3651080930869013
0.5777920899992666 1.0315507012238316 0.3693785011715264
0.6683465717919475 0.0360591243671101 0.2052175704554934
0.2461635354286980 0.7646904079192439 1.0088582212635258
0.1929460786015471 0.4695788181732047 0.9713205045347785
0.2544767474283657 0.6042464935788584 0.7524738134375027
0.4999761863989385 0.7515317994519944 0.7665001210207051
0.0309180624066225 0.3163339995046774 0.1347567766248343
0.3990104937141098 0.6006789121661388 0.5058773087184222
0.2576858555569930 0.5235309594581270 0.2416907439307473
0.4144971971941396 0.7686322625981743 0.2810869260959064
0.4943123046718796 0.4766100916293045 0.7292847824974398
0.7671116528768356 0.7505861438867950 0.8995511938265869
0.7714594379354743 -0.1546114721232754 0.4118297917934841
0.6283128392236679 0.2311667717042626 0.7460434956289651
0.5275680358540571 0.8061579569132280 0.0429691384248918
0.2968865970536227 0.2624288551673952 0.1321171012590240
0.7754769616700143 0.4911785016026120 0.7769835929218475
0.5408379239759858 0.5186213832884250 0.2838091624448219
0.7565175078652816 0.2740359040090061 1.0031983605845767
0.6481715631514586 0.3809817005427188 0.5214958996079175
0.5602367809462726 0.2211393597638454 0.2570644664216546
0.6608275132469794 0.0015187153856674 0.9010628713444417
0.0198916040265591 0.6492298332472982 0.8468598238097034
0.3816524924537844 1.0051725332567569 0.1486110640393346
0.4863846869445549 0.1914222540095063 -0.0266452673157057
0.8225005190139213 0.4172594014525411 0.2808679139549766
0.8856962798265164 0.5351541251082332 0.0475187903703036
0.1806197042381862 0.6225689820073963 0.8984511784040251
0.4076220786951815 0.6068328529013632 0.6775009732270258
0.6153884877400786 0.8232232122687801 0.8953881250309655
0.3815606975119766 0.8369572897692317 0.1216199438412733
0.4034546096799853 0.6040095855114854 0.3278637646322450
0.1915059761588622 0.3922017886197338 0.1139889571015823
0.6364140125469393 0.3935889238678646 0.6951049951672011
0.8608648291991334 0.6122376681441427 0.8985422307313435
0.6350467485694524 0.1738765626303231 0.9073598305913078
0.4257035566719548 0.1703988869323599 0.1266349459854781
0.8685124222531690 0.3807300481379052 0.1178832273821222
POSCAR.4_init
Reaction pathway test
1.00000000000000
7.9487550000000002 -0.0000500000000000 4.9014329999999999
-3.9745640000000000 6.8838939999999997 4.9012549999999999
-3.9744389999999998 -6.8837770000000003 4.9012460000000004
Al C H O Si
1 2 3 25 11
Direct
0.6387756944561276 0.4012120876851310 0.3532734792269315
0.7392343593244277 0.9037506852494360 0.4727532507676848
0.5853525082706864 0.0832877936108952 0.3587998921755486
0.6117276047767487 0.0222282646495838 0.2611754023467227
0.6858441531998257 0.1402468495791714 0.3981392285899252
0.5157470600221065 0.0279714507590067 0.4258491472348496
0.2685775062272701 0.7373695209230580 -0.0073427159564949
0.1846392768271981 0.4637592773959613 0.9668262173438991
0.2622581576714780 0.5883192719928664 0.7407423947828007
0.5259844617132652 0.7532627959750469 0.7385926402177541
0.0358740885425032 0.3280248573530538 0.1479747583747389
0.3504199597113813 0.6610231358646341 0.5041405600600197
0.2474096197762262 0.5327501790952712 0.2463219783564613
0.4114357070599692 0.7798458969713520 0.2578836353184217
0.5032375050041068 0.4784564125315423 0.6761623616365177
0.7733271621699099 0.7460983286866631 0.9138668421783351
0.8431121546289226 0.8606309883257968 0.5142858996674572
0.5961783925857385 0.2283420817121579 0.7475530371163474
0.5542462030547020 0.8265939395664245 0.0346023981383468
0.3124385339968538 0.2839056807596789 0.1172521374686578
0.7822986922774250 0.4768020594008841 0.7736197844608504
0.5291257994243568 0.5419010688464461 0.3255884750864870
0.7993192179824421 0.2444965761521856 -0.0060179676566297
0.6937276002119118 0.3727081908605321 0.5210271280549834
0.5120416972421128 0.2100918920042797 0.2991193158633122
0.6735094749840188 -0.0062699305063081 0.8518818284055982
0.0146571466006417 0.6597218345636776 0.8287313924791077
0.3459886347546634 -0.0029496686151194 0.1278319257316917
0.5280582808742778 0.1851368263021938 1.0059916428324733
0.7672133162716521 0.4014369192874723 0.2365201250777685
0.8896990363069027 0.5396344651599320 0.0492726405130849
0.1711342839035838 0.6116758034219585 0.8873145563665112
0.4101847477077599 0.6174595613888860 0.6612715163575200
0.6324928371741062 0.8282507332618047 0.8757528064059821
0.3956581716046995 0.8357686611053713 0.0958405331776555
0.3947252240856924 0.6300470052109007 0.3338405891717605
0.1928641844839093 0.4040670786636100 0.1209468498553622
0.6430222881148094 0.3782655919260429 0.6857704451475711
0.8604595151127654 0.6010769276650000 0.8905420685780717
0.6574064095757256 0.1623457484888207 0.9013356984790826
0.4217782217579634 0.1591943128603318 0.1354733720559321
0.8729090052964774 0.3767673883767818 0.1109387177939303
POSCAR.5_init
Reaction pathway test
1.00000000000000
7.9487550000000002 -0.0000500000000000 4.9014329999999999
-3.9745640000000000 6.8838939999999997 4.9012549999999999
-3.9744389999999998 -6.8837770000000003 4.9012460000000004
Al C H O Si
1 2 3 25 11
Direct
0.6394726655281282 0.3965033118891959 0.3435765943359740
0.7050550646684141 0.8649386767332974 0.4415312700225982
0.6402564564936992 0.0805437513428449 0.3635250977157916
0.7433944952795252 0.1231866947929046 0.3226673040389457
0.6419818023195263 0.0955066011057857 0.4818246248112532
0.5722357206608255 -0.0068165977971404 0.2982245821491078
0.2487985096946962 0.7470141012659125 1.0117629304159572
0.2092797619141816 0.4712293306642353 0.9710853802608856
0.2742063778943514 0.6069089249010319 0.7592903966196974
0.5192150775420471 0.7643040915402454 0.7523873280493134
0.0396986711827433 0.3319430980419775 0.1269514780623325
0.3583370291335146 0.6717332531842007 0.5050650615151714
0.2506182710437612 0.5323237574814585 0.2640786935945655
0.4168516154774761 0.7736441815584135 0.2587009365626186
0.4772481222073754 0.4614693002545006 0.6715463771134714
0.7878619599326696 0.7595012202244232 0.9188860574357237
0.7656877373515107 0.7836583795333585 0.3787935282607061
0.6233893237686640 0.2442331639096545 0.7548305799586166
0.5467104541105027 0.8266837890482943 0.0287646298568407
0.2960383694972200 0.2739597354924117 0.1282323649220623
0.7549089904412198 0.5040114972024541 0.7836032965788985
0.5297467574307748 0.5300134105129843 0.3223988185582514
0.7898664319080168 0.2566073345225813 0.9930291923409488
0.6786788316104344 0.3733389026578841 0.5171309625275771
0.5344624097302459 0.2064537673692249 0.2902190484148430
0.6807403134420077 0.9986683119404198 0.8615395351406613
0.0167568854347997 0.6366352211167789 0.8593563227324797
0.3495104915500716 1.0048970608035224 0.1115493952110750
0.5107463007727520 0.1986585298992333 -0.0157353911917621
0.7961073734050081 0.3819439472668675 0.2610677094599224
0.8787772332641948 0.5286898804370725 0.0573320797964962
0.1822874555767301 0.6177169262757385 0.8948172903451975
0.4147252172873378 0.6150836035186108 0.6658054486035624
0.6380542369014253 0.8263548552528495 0.8874824750059908
0.3870405664498269 0.8361356768526701 0.1001853242394696
0.3988202123980464 0.6230479428166796 0.3433213120422873
0.2028990980552815 0.4084709458888915 0.1267012524668510
0.6349164092278642 0.3992807125436610 0.6756560667537945
0.8568267020800896 0.6072747709011483 0.9082292677149003
0.6518477507660498 0.1649247800617258 0.8947984529637162
0.4269847175234723 0.1679864582623027 0.1232702755036939
0.8734698709186332 0.3715059939032876 0.1206404718203736
POSCAR.6_init
Reaction pathway test
1.00000000000000
7.9487550000000002 -0.0000500000000000 4.9014329999999999
-3.9745640000000000 6.8838939999999997 4.9012549999999999
-3.9744389999999998 -6.8837770000000003 4.9012460000000004
Al C H O Si
1 2 3 25 11
Direct
0.6423672312925792 0.3914850453159373 0.3400425466748950
0.7113955037438986 0.9150346283499408 0.5361637305991392
0.6256044516886795 0.0551662134801337 0.3943352247655982
0.5162749964367519 0.0035701314701795 0.3935168700605537
0.7140728048283090 0.0121380363125870 0.3365274866481253
0.6582941106261271 1.1363464635128184 0.4835136982702529
0.2903500707381216 0.7731800188676460 -0.0423139976648281
0.2359661582119096 0.4900679214577685 0.9749210933728444
0.2473561763220278 0.6041119970489047 0.7186092659427901
0.5070876796866493 0.7667293667833623 0.7534764536379063
0.0291740012153105 0.3385269947399128 0.1052594622204845
0.3849917104141334 0.6588292726067998 0.4806452863580031
0.2303919561968580 0.5199981562883159 0.2487505339741046
0.3940817805007507 0.7533982823857752 0.2401389309157462
0.4823263318487423 0.4814926580815018 0.6668456994274861
0.7857445213159805 0.7620239941269132 0.9032104297935665
0.7792843611069514 0.8276324475300932 0.5724090899382076
0.6215124569151401 0.2452469315075300 0.7645437599201244
0.5595808980484352 0.8414657934259560 0.0322634754262991
0.2977831592119001 0.2692436135739456 0.1512228646664003
0.7768084519776864 0.5145058889740337 0.7473047926463240
0.5139981545212536 0.4990967394432917 0.2862923083749087
0.7845312466551515 0.2221409139531066 1.0178061593251491
0.6825052105781074 0.3720488125926638 0.5145622471629652
0.5624743513578218 0.2076204801002911 0.2759986731429730
0.6697362517015619 1.0025198840712872 0.8448228609250398
0.0199224900817955 0.6435612991509839 0.8721030983421870
0.3702033645187299 -0.0036364203820513 0.1611345496412056
0.5043183328393032 0.1639627043729344 0.9949858005571807
0.8051450166676021 0.4101843159508998 0.2692457920482184
0.8371688201593205 0.5138718190062973 0.0295552914911029
0.1983553210518001 0.6258601726506766 0.8813551306128060
0.4058381274673998 0.6344847956179060 0.6507502529876412
0.6338185145645514 0.8350277098113338 0.8853588055181932
0.3991311588377134 0.8356495551712334 0.0953513064550183
0.3911171834396501 0.6002193105864239 0.3085370048234105
0.1989202342917561 0.3999680217879398 0.1133380067565643
0.6347848313609805 0.4018949490736654 0.6729534395936464
0.8546608072772587 0.6062024594534280 0.8976671529835299
0.6470230906492567 0.1603855806068904 0.9047280238924675
0.4268292237329298 0.1625310599657356 0.1422557548225977
0.8637725024576033 0.3709009062296315 0.1197255101943291
POSCAR.7_init
Reaction pathway test
1.00000000000000
7.9487550000000002 -0.0000500000000000 4.9014329999999999
-3.9745640000000000 6.8838939999999997 4.9012549999999999
-3.9744389999999998 -6.8837770000000003 4.9012460000000004
Al C H O Si
1 2 3 25 11
Direct
0.6423876815227452 0.3704470785418887 0.3410101904614561
0.6522087642498382 0.9032233041437587 0.4749145040174509
0.6038660421620968 0.0523365887012187 0.3926370121623233
0.6649259180985949 0.0323027163400948 0.2973669308750807
0.6424514935176628 0.1380611208214780 0.4838000933086251
0.4854003459885129 0.0242014462205816 0.3876469120288183
0.2600258526685770 0.7720354340238853 1.0050687127778255
0.2064495120980419 0.4841845675455626 0.9715195901435231
0.2653176222444578 0.6210747453017822 0.7421568863347959
0.5340503981241044 0.7288462106864937 0.7655010771927696
0.0417590756614686 0.3231019962381220 0.1332421850681333
0.3793467781084565 0.6690254955028867 0.5009781960612938
0.2522781247008773 0.5446264698136632 0.2426699847831768
0.4327794465028360 0.7617624561515465 0.2425167383181072
0.4777982565293319 0.4594405640026710 0.6530986437421354
0.7747538928864915 0.7538452578298634 0.9393772916972319
0.7006972396637633 0.8110700972135022 0.5306441956410269
0.6149741694183757 0.2495594862007427 0.7460284235994441
0.5516139849001016 0.8552598802531722 0.0374985027423479
0.3048469148471732 0.2797770190641420 0.1054894222085869
0.7721432749067342 0.5012892654614191 0.7661279613492491
0.5333804887863809 0.5170801731596528 0.3136735627923074
0.8006306662174166 0.2419356240906301 -0.0057921741220461
0.6893926835316863 0.3589731828300982 0.5191018814904984
0.5315803373657549 0.2082946959068813 0.2726074900725021
0.6861175111781188 0.0020157529126551 0.8276471213134596
0.0155302997862507 0.6613528127563832 0.8594672516249577
0.3497152519844348 1.0100842137656876 0.1412315593873151
0.5145725865046604 0.1654380193824488 0.9781904112371749
0.7983886850162647 0.3732223742807453 0.2603023472089125
0.8813282059198417 0.5216719788780378 0.0464179379486477
0.1899354107544340 0.6297383564568408 0.8897287324111213
0.4135127748372226 0.6155399173271637 0.6659078616820278
0.6357890357178202 0.8382337478422746 0.8823281515593231
0.3976363157607455 0.8476602181020508 0.1063331272318617
0.4051811707472145 0.6102291598847323 0.3278317029804668
0.1949103484030432 0.4019071466862751 0.1167170076797476
0.6394027760091172 0.4009429307522819 0.6740700745668246
0.8613136225000289 0.6077349236638317 0.9059736570728176
0.6651523555964665 0.1695995116661321 0.8892214373157702
0.4210449314788592 0.1658346899295627 0.1217754555106620
0.8768946068406328 0.3692625260144742 0.1166532496020828
INCAR.init
! Molecular dynamics
IBRION = 0 ! Molecular dynamics
ISYM = 0 ! Symmetry setting
NPAR = 2 ! Number of bands treated in parallel
NSW = 10000 ! Number of steps in MD run
POTIM = 1.0 ! Time step
MDALGO = 5 ! CSVR thermostat
TEBEG = 440. ! Temperature
CSVR_PERIOD = 20 ! Time scale of CSVR in terms of the number of MD steps
! Machine learning
ML_LMLFF = T ! Enables the use of MLFF
ML_MODE = run ! Prediction-only mode
ML_OUTPUT_MODE = 0 ! Reduced file output mode
ML_OUTBLOCK = 10 ! Output distance in number of steps
LBLUEOUT = T ! Blue moon ensemble
K-Points
0
Gamma
1 1 1
0 0 0
Pseudopotentials of Al, C, H, O, and Si.
R 3 25 0
R 3 2 0
S -1 1 7
The KPOINTS, POTCAR, and ICONST files remain the same as in Example 10. The POSCAR files are taken from the seven structures generated in Example 11. The INCAR file is similar to that used in Example 11, with the increase of NSW from 4000 to 10000.
12.3 Calculation¶
Open a terminal, navigate to this example's directory.
cd $TUTORIALS/transition_states/e12_*
Seven constrained MD runs are necessary to describe the seven different configurations of $\xi$ along the reaction profile. Note, that the full calculation would take approximately 25 minutes. Instead, perform a single MD run using the reference POSCAR.1_init and later take the remaining six calculations from the reference data provided. Create the directory 1, navigate to it, and then run the first MD:
mkdir 1
cp POSCAR.1_init ./1/POSCAR
cp INCAR.init ./1/
cp ICONST KPOINTS POTCAR ./1
cd ./1
ln -s ../../mlff/MLFF_e05-12/ML_FFN ML_FF
cp POSCAR POSCAR.1
cp ../single_run .
bash single_run
Click here to reveal the single_run script for the first configuration!
#!/usr/bin/bash
ulimit -s unlimited
runvasp="mpirun -np 4 vasp_gam"
drct=$(pwd)
stepmax=2
#c initialize RNG
aaa=$(echo "RANDOM_SEED = 817605623 0 0" )
step=1
while [ $step -le $stepmax ]
do
cp POSCAR POSCAR.$step
cp INCAR.init INCAR
echo $aaa >> INCAR
$runvasp
#c extract seed of RNG from the last step
aaa=$(grep RANDOM_SEED REPORT|tail -1 )
#c save important output files
cp REPORT report.$step
cp vasprun.xml vasprun.xml.$step
cp CONTCAR POSCAR
let step+=1
done
Click here to reveal the run script for all seven configurations!
Copy your POSCAR files from Exercise 11 using the following script:
for a in 1 2 3 4 5 6
do
cp ../../e11_*/POSCAR.$a POSCAR.${a}_init
done
cp ../../e11_*/CONTCAR POSCAR.7_init
Then run the calculations:
bash run
#!/usr/bin/bash
ulimit -s unlimited
export FILES_PATH="../../mlff/MLFF_e05-12/"
#export FILES_PATH="/srv/data/transition_states/mlff/MLFF_e05-12/ML_FFN ML_FF"
runvasp="mpirun -np 4 vasp_gam"
drct=$(pwd)
stepmax=2
#c initialize RNG~
aaa=$(echo "RANDOM_SEED = 817605623 0 0" )
#c loop over integration points
for p in 1 2 3 4 5 6 7
do
cd $drct
if ! test -d $p
then
mkdir $p
fi
cp POSCAR.${p}_init ${p}/POSCAR
cp ICONST INCAR.init $p
cd $p
ln -s $FILES_PATH/ML_FFN ML_FF
cp $FILES_PATH/KPOINTS .
cp $FILES_PATH/POTCAR .
step=1
while [ $step -le $stepmax ]
do
cp POSCAR POSCAR.$step
cp INCAR.init INCAR
echo $aaa >> INCAR
$runvasp
#c extract seed of RNG from the last step
aaa=$(grep RANDOM_SEED REPORT|tail -1 )
#c save important output files
cp REPORT report.$step
cp vasprun.xml vasprun.xml.$step
cp CONTCAR POSCAR
let step+=1
done
done
Congratulations on performing a constrained MD run to model the reaction!
Once this calculation has finished, if you have performed only the first step, then you must copy the remaining results from ./refdata to the $TUTORIALS/transition_states/e12_* directory for analysis:
cd $TUTORIALS/transition_states/e12_*/refdata
cp -r 2 3 4 5 6 7 ../
With all seven calculation results, you can analyze the MD simulation for the entire reaction. For each constrained MD step, the REPORT files are saved as report.* and contain information about the simulation:
>Blue_moon
lambda |z|^(-1/2) GkT |z|^(-1/2)*(lambda+GkT)
b_m> 0.568493E+00 0.145516E+01 -0.387243E-04 0.827192E+00
To obtain the free energy gradient for a given $\xi$ value, the ensemble average of the fourth column |z|^(-1/2)\*(lambda+GkT) must be divided by the ensemble average of the third column |z|^(-1/2). The relevant equations in the introduction are repeated here for ease of use:
$ \Large{(\frac{\partial A}{\partial \xi})}_{\xi^*} = \frac{1}{\langle|Z|^{-1/2}\rangle_{\xi^*}} \langle \large{|Z|^{-1/2} [\lambda + \textit{G} k_{B}T ]} \rangle _{\xi^*} $
The script fgradBM.sh calculates this equation for all $\xi$ points, using:
cd $TUTORIALS/transition_states/e12_*/
../scripts/bashScripts/fgradBM.sh > grad.dat
fgradBM.sg creates a file, grad.dat that contains the free energy gradient vs. $\xi$. Plot grad.dat using the following script. The $\xi$ axis has been inverted such that the $R \rightarrow TS$ direction goes from left to right, as is conventional. This convention is followed throughout:
import pandas as pd
import plotly.graph_objects as go
# Load the data
df1 = pd.read_csv("./e12_Dynamic_Thermodynamics/grad.dat", sep=" ", header=None)
#df2 = pd.read_csv("./e12_Dynamic_Thermodynamics/refdata/grad.dat", sep=" ", header=None)
# Create the trace for the calculated data
calc_trace = go.Scatter(
x=df1[0],
y=df1[1],
mode='lines',
line=dict(color='#4c265f', width=2),
name='calc.',
)
# Add vertical lines for POSCAR.1 and POSCAR.7
poscar1_line = go.Scatter(
x=[df1[0][0], df1[0][0]],
y=[min(df1[1]), max(df1[1])],
mode='lines',
line=dict(color='black', dash='dot'),
name='POSCAR.1',
showlegend=False,
)
poscar7_line = go.Scatter(
x=[df1[0][6], df1[0][6]],
y=[min(df1[1]), max(df1[1])],
mode='lines',
line=dict(color='black', dash='dot'),
name='POSCAR.7',
showlegend=False,
)
# Add annotations for POSCAR.1 and POSCAR.7
annotations = [
dict(
x=df1[0][0] + 0.02,
y=-1,
text="POSCAR.1",
showarrow=False,
font=dict(size=12),
textangle=90,
),
dict(
x=df1[0][6] - 0.02,
y=-1,
text="POSCAR.7",
showarrow=False,
font=dict(size=12),
textangle=90,
),
]
# Invert the x-axis by reversing its range
layout = go.Layout(
xaxis=dict(title="ξ in Å", autorange="reversed"), # Reverse x-axis
yaxis=dict(title="Free energy gradient in eV Å⁻¹"),
annotations=annotations,
legend=dict(x=0.05, y=0.95),
font_color="black",
template="ggplot2",
)
# Combine all traces
fig = go.Figure(data=[calc_trace, poscar1_line, poscar7_line], layout=layout)
# Show the plot
fig.show()
The free energy gradient is largely constant between $\xi$ ~2.5 and 1.5 Å. This $\xi$ range corresponds to POSCAR.1, POSCAR.2, and POSCAR.3 files for the reference data and the reactant ensemble. The free energy gradient then decreases to 0.5 Å before increasing rapidly to POSCAR.7. The TS lies around $\xi \approx$ 0.0 Å.
We want to calculate the free energy of activation. For the present step, we will repeat the equation from the introduction to this exercise:
$ \normalsize {\Delta A_{\xi_{ref,R} \rightarrow \xi^{*}}} = \Large {\int_{\xi_{ref,R}}^{\xi^{*}} {(\frac{\partial A}{\partial \xi})}_{\xi^*} \,} \normalsize {d\xi} $
The change in free energy from $R \rightarrow TS$ requires the integration of each point over the free energy gradient with respect to $\xi$. The Simpson method can be used to integrate the above equation. Integrate grad.dat along $\xi$ using the Simpson method with the script simpsonI_non_cum.py:
cd $TUTORIALS/transition_states/e12_*/
python3 ../scripts/simpsonIntegration_nonuniform/simpsonI_non_cum.py grad.dat > dA.dat
The resulting points dA.dat may be interpolated using cubic splines with the splinesCubic.py script:
cd $TUTORIALS/transition_states/e12_*/
python3 ../scripts/splinesCubic/splinesCubic.py dA.dat 1000 > dA_splines.dat
Plot dA.dat and dA_splines using the following script:
import pandas as pd
import plotly.graph_objects as go
# Load the data
df1 = pd.read_csv("./e12_Dynamic_Thermodynamics/dA.dat", sep=" ", header=None)
df2 = pd.read_csv("./e12_Dynamic_Thermodynamics/dA_splines.dat", sep=" ", header=None)
# Find coordinates of the maximum and minimum free energy
xi_max = df2[df2[1] == df2[1].max()][0].values[0]
A_max = df2[1].max()
xi_min = df2[df2[1] == df2[1].min()][0].values[0]
A_min = df2[1].min()
print(f"Max. free energy = {A_max:.3f} eV")
print(f"Corresponding reaction coordinate = {xi_max:.3f} Å")
# Create traces for the plot
calc_trace = go.Scatter(
x=df1[0],
y=df1[1],
mode='lines',
line=dict(color='#4c265f', width=2),
name='calc.',
)
spline_fit_trace = go.Scatter(
x=df2[0],
y=df2[1],
mode='lines',
line=dict(color='#2fb5ab', width=2),
name='spline fit',
)
# Add annotations for maximum and minimum points
annotations = [
dict(
x=xi_min,
y=A_min + 0.1,
text="\u03BE<sub>ref,R</sub>",
showarrow=False,
font=dict(size=12),
),
dict(
x=xi_max + 0.08,
y=A_max + 0.08,
text="\u03BE<sup>*</sup>",
showarrow=False,
font=dict(size=12),
),
]
# Set layout with inverted x-axis
layout = go.Layout(
xaxis=dict(title="ξ in Å", autorange="reversed"), # Reverse x-axis
yaxis=dict(title="Free energy in eV", range=[0.0, 1.2], dtick = 0.2),
annotations=annotations,
legend=dict(x=0.05, y=0.95),
font_color="black",
template="ggplot2",
)
# Combine traces into a figure
fig = go.Figure(data=[calc_trace, spline_fit_trace], layout=layout)
# Show the plot
fig.show()
The integrated free energy gradient provides the free energy $A$ with respect to the reaction coordinate $\xi$. For the reference data, $A$ increases from ~2.4 Å (the reactant state $\xi_{R,ref}$) to a maximum at $\xi$ = 0.043 Å. The spline fit smooths out the jagged plot. The maximum of the free energy profile $A$ corresponds to the transition state, i.e. $\xi^*$ (cf. dA_splines.dat). The value of $\xi^*$ found for the TS using the static approach was (0.01 Å), meaning that the value from the dynamic approach is similar and the two approaches are equivalent for modeling the reaction. Remember that, in Example 10, the reaction coordinate for the reactant was very similar
The difference between the reactant and transition state is the key part of modeling the reaction. The value of $\Delta A_{\xi_{ref,R} \rightarrow \xi^{*}}$ is obtained as a difference in $A$ between TS and the reference reactant R (i.e. the initial state), which gives a value of 1.11 eV.
The final component required to calculate the free energy of activation is the velocity associated with the reaction coordinate $\dot{\xi}$, which is defined as:
$\langle | \dot\xi^* | \rangle = \sqrt{\Large{\frac{2k_BT}{\pi}}} \Large{\frac{1}{\langle|Z|^{-1/2}\rangle_{\xi^*}}} $
This velocity is fully determined by an ensemble average of |z|^(-1/2) from constrained MD for TS. Since $\xi$ = 0.043 Å corresponding to the TS was not used as an integration point in the Blue moon simulation, we cannot directly access |z|^(-1/2) for the TS.
Instead, either an extra constrained MD calculation would need to be run for that point, or the average |z|^(-1/2) can be estimated from the values computed for the nearest integration points. The script zetBM.sh extracts |z|^(-1/2) from each of the seven MD runs using the following:
cd $TUTORIALS/transition_states/e12_*/
echo "Calc."
echo "ξ <|z|^(-1/2)>"
cd ./e12_*
../scripts/bashScripts/zetBM.sh > z_m_hlf.dat
cat z_m_hlf.dat
cd ..
echo ""
echo "Ref."
echo "ξ <|z|^(-1/2)>"
cd ./e12_*/refdata
../../scripts/bashScripts/zetBM.sh > z_m_hlf.dat
cat z_m_hlf.dat
%%bash
cd ./e12_*
echo "Calc."
echo "ξ <|z|^(-1/2)>"
../scripts/bashScripts/zetBM.sh > z_m_hlf.dat
cat z_m_hlf.dat
It does not change much after the first few points, as was also seen for $\langle|Z|^{-1/2}\rangle$. However, this is clearer when plotted:
import pandas as pd
import plotly.graph_objects as go
# Load the data
df1 = pd.read_csv("./e12_Dynamic_Thermodynamics/z_m_hlf.dat", sep=" ", header=None)
#df2 = pd.read_csv("./e12_Dynamic_Thermodynamics/refdata/z_m_hlf.dat", sep=" ", header=None)
# Create traces for the calculated and reference data
calc_trace = go.Scatter(
x=df1[0],
y=df1[1],
mode='lines',
line=dict(color='#4c265f', width=2),
name='calc.',
)
# Add annotations for POSCAR points
annotations = [
dict(x=df1[0][i], y=1.42, text="ξ<sub>"+str(i+1)+"</sub>", showarrow=False, font=dict(size=12, color='black'), opacity=0.7)
for i in range(7)
]
# Set layout with inverted x-axis
layout = go.Layout(
xaxis=dict(title="ξ in Å", autorange="reversed"), # Reverse x-axis
yaxis=dict(title="⟨|Z|<sup>-1/2</sup>⟩ in amu<sup>-1/2</sup>", range=[1.4, 1.75]),
yaxis_range=[1.4,1.75],
annotations=annotations,
legend=dict(x=0.05, y=0.95),
font_color="black",
template="ggplot2",
)
# Combine traces into a figure
fig = go.Figure(data=[calc_trace], layout=layout)
# Show the plot
fig.show()
In our case, |z|^(-1/2) is nearly constant in the vicinity of the transition state ($\xi \approx$ 0.0 Å). It is clear from the output that the average |z|^(-1/2) barely changes for each point (based on POSCAR.n structures and denoted $\xi_n$ in the plot) from 3 to 7, so point 6 (cf. POSCAR.6) is a very good approximation for our TS. We now have a reference structure for the reactant state (POSCAR.1) and for the transition state (POSCAR.6). All that remains is to calculate the average ${\dot\xi}$. The script aux_vel.sh extracts $\dot{\xi}$. Note: Change the POSCAR.6 file in the scripts below to your POSCAR.* that lies closest to your TS reaction coordinate, i.e. $\xi \approx$ 0.0 Å. Enter the following script into your terminal to extract $\dot{\xi}$.
cd $TUTORIALS/transition_states/e12_*/6
../../scripts/bashScripts/aux_vel.sh 440 > xi_dot.dat # at 440 K
cd ../..
cd ./e12_*/refdata/6
../../../scripts/bashScripts/aux_vel.sh 440 > xi_dot.dat # at 440 K
cd ../../..
# Probability for calculating delta A
grep mc e10_*/refdata/report.1 | tail -1 | awk '{print $3}' > e10_Dynamic_Coordinates/refdata/xi.dat
grep mc e10_*/report.1 | tail -1 | awk '{print $3}' > e10_Dynamic_Coordinates/xi.dat
The value of ${\dot\xi}$ obtained in this way is 1.05285e+13 Å s$^{-1}$, which is the last remaining term required to calculate $\Delta A_{R\rightarrow TS}^{\ddagger}$. Remind yourself of the equation to calculate $\Delta A_{R\rightarrow TS}^{\ddagger}$:
$\normalsize{\Delta A_{R \rightarrow TS}^{\ddagger}} = \normalsize{{\Delta A_{\xi_{ref,R} \rightarrow \xi^{*}}}} - k_BT \ln[\Large{\frac{h}{k_BT} \frac{\langle | \dot\xi^* | \rangle}{2}} \normalsize{P(\xi_{ref,R})}] $
$\Delta A_{\xi_{ref,R} \rightarrow \xi^{*}}$ is the difference between the free of the transition state and that of the reactant state ($A_{TS} - A_R$). We have just obtained $\dot\xi^*$. The final value is the probability density of the reactant state $P(\xi_{ref,R})$, which we have found in Example 10. The following script extracts the free energies and calculates the difference of ($A_{TS} - A_R$) from ./e12_*/dA_splines.dat, and reads $\dot\xi^*$ from ./e12_*/6/xi_dot.dat. It also finds the value of $\xi_{ref,R}$, which is the first structure in ./e10_*/xi.dat and the corresponding $P(\xi_{ref,R})$. Finally, the script brings these terms together and outputs the energy from the above equation for $\Delta A_{\xi_{ref,R} \rightarrow \xi^{*}}$.
Input your $\Delta A_{R \rightarrow TS}^{\ddagger}$ from part2 of the exercises for the static and semi-classical values into the following script where it says delta_A_stat and delta_A_SC, respectively.
import math
import pandas as pd
# Function to calculate the free energy of activation according to the above equation
# Takes the free energy difference between R and TS delta_A_MD, the temperature, the probability density of the reactant state P, and the velocity associated with the reaction coordinate xi_dot as inputs
def calc_delta_A(delta_A_MD, T, P, xi_dot):
# Settings
k_b = 8.617333262E-5
h = 4.1356692E-15
T = 440. # T in Kelvin
# Calculate correction
A = h / (k_b * T)
corr = k_b * T * math.log(A * (xi_dot/2) * P)
# Free energy of activation
delta_A = delta_A_MD - corr
return(delta_A, corr)
P, xi_dot = 0.0, 0.0
# Get xi_dot from file
df1 = pd.read_csv("./e12_Dynamic_Thermodynamics/6/xi_dot.dat", sep=" ", header=None)
#df2 = pd.read_csv("./e12_Dynamic_Thermodynamics/refdata/6/xi_dot.dat", sep=" ", header=None) # reference data
xi_dot = df1[0][0]
#xi_dot_ref = df2[0][0]
# Unhash for using the reference data - note, you must have already for xi.dat and mc.dat.hist from the reference data in Exercise 10
# Get P(xi_ref,R) from file - assume from POSCAR.2 in Exercise 10:
#df3 = pd.read_csv("./e10_Dynamic_Coordinates/refdata/xi.dat", sep=" ", header=None) # reference data
#df4 = pd.read_csv("./e10_Dynamic_Coordinates/refdata/mc.dat.hist", sep="\t", header=None) # reference data
#xi_temp1 = df3[0][0] # the reaction coordinate for the reactant state \xi_{ref,R}
#P1 = float(df4[df4[0].between(xi_temp1 - 0.01, xi_temp1 + 0.01)][1].iloc[0]) # Probability density of the reactant state P(\xi_{ref,R})
df5 = pd.read_csv("./e10_Dynamic_Coordinates/xi.dat", sep=" ", header=None) # calculated data from Exercise 10 -
df6 = pd.read_csv("./e10_Dynamic_Coordinates/mc.dat.hist", sep="\t", header=None) # calculated data from Exercise 10
xi_temp2 = df5[0][0] # the reaction coordinate for the reactant state \xi_{ref,R}
diff = 0.03 # THIS CHANGES HOW BROAD THE FILTER FOR FINDING THE REACTION COORDINATE FROM xi.dat IN xi.dat.hist - IF YOU CANNOT FIND IT, INCREASE THE DEGREE OF UNCERTAINTY
P2 = float(df6[df6[0].between(xi_temp2 - diff, xi_temp2 + diff)][1].iloc[0]) # Probability density of the reactant state P(\xi_{ref,R})
# Free energy difference between R and TS - dynamic thermodynamics
df7 = pd.read_csv("./e12_Dynamic_Thermodynamics/dA_splines.dat", sep=" ", header=None) #
# Finds coordinate and free energy of the maximum and minimum
xi_max = float(df7[df7[1]==df7[1].max()][0].to_string()[5:])
A_max = float(max(df7[1]))
xi_min = float(df7[df7[1]==df7[1].min()][0].to_string()[5:])
A_min = float(min(df7[1]))
A_diff = (A_max - A_min) # free energy difference between R and TS
delta_A, corr = calc_delta_A(A_diff, 440., P2, xi_dot)
# Free energies of activation for the semi-classical 'delta_A_SC', static 'delta_A_stat', or dynamic 'delta_A_dyn' approaches.
delta_A_SC = 1.0666106712102987
delta_A_stat = 1.0484333798001444
delta_A_dyn = delta_A
print("ΔA^‡ - semi-classical: " + str("{:.3f}".format(delta_A_SC)) + " eV")
print("ΔA^‡ - static QM: " + str("{:.3f}".format(delta_A_stat)) + " eV")
print("ΔA^‡ - dynamic QM: " + str("{:.3f}".format(delta_A_dyn)) + " eV")
print("Difference between static QM and semi-classical approaches: " + str("{:.3f}".format(delta_A_stat - delta_A_SC)) + " eV")
print("Difference between static and dynamic QM approaches: " + str("{:.3f}".format(delta_A_dyn - delta_A_stat)) + " eV")
print("Difference between dynamics QM and semi-classical approaches: " + str("{:.3f}".format(delta_A_dyn - delta_A_SC)) + " eV")
We obtain a value of 1.14 eV, which is ~0.1 eV higher than the static result. Recalling the difference between static and semi-classical approaches was 0.02 eV, this difference is much more significant.
Can you (qualitatively) explain what causes this increase?
The increase is due to the inclusion of the ensemble average, rather than relying on only one state.
import scipy
import math
# define the Eyring equation, taking in the temperature T and the free energy of activation delta A as input to calculate the rate constant k
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_stat = eyring(440., delta_A_stat)
k_QM_dyn = eyring(440., delta_A_dyn)
print("Semi-classical rate constant: " + str("{:.3f}".format(k_SC)) + " s^-1")
print("Static QM rate constant: " + str("{:.3f}".format(k_QM_stat)) + " s^-1")
print("Dynamic QM rate constant: " + str("{:.3f}".format(k_QM_dyn)) + " s^-1")
You can see that the small difference of ~0.1 eV in the free energy of the reaction makes a much larger difference in the rate constant, around an order of magnitude.
Congratulations on completing the exercises. You have succeeded in calculating the free energy of activation using both static and dynamic approaches!
12.4 Questions¶
- How big a difference is there between the free energy of activation using static and dynamics approaches?
- How big an impact on the rate constant does a small difference in the free energy of reaction make?