2. Practice in energy decomposition analysis with XEDA
2.1. Examples
2.1.1. H-bond in NH3\(\cdots\)H2O
This is a simple example of calculating the hydrogen bonding interactions between water and ammonia.
2.1.1.1. Input File
$CTR
METHOD=RHF BASIS=aug-cc-pVTZ
NMUL=1 CHARGE=0
DFT=wB97XD
$END
$GEO
N 1.37515100 -0.01976400 -0.00003700
H 1.80210900 -0.41927400 0.82476800
H 1.60816800 0.96438600 -0.01647200
H 1.80722600 -0.44771800 -0.80772500
O -1.54314600 0.10499800 -0.00000100
H -0.57947200 -0.02462600 -0.00024700
H -1.91891900 -0.77439900 -0.00006000
$END
$EDA
NMOL = 2
MATOM = 4 3
MMULT = 1 1
MCHARGE = 0 0
$END
2.1.1.2. Computational Results
\(\Delta\)Eele |
\(\Delta\)Eexrep |
\(\Delta\)Epol |
\(\Delta\)Ecorr/disp |
\(\Delta\)Etotal |
|---|---|---|---|---|
-11.64 |
12.14 |
-4.41 |
-2.88 |
-6.78 |
Generally, in hydrogen bond, electrostatic interaction is dominant rather than polarization.
2.1.1.3. Question
What is the change in results when the basis set is replaced with cc-PVTZ?
2.1.2. Covalent bond in CH3-NH2
2.1.2.1. Input File
$CTR
METHOD=RHF BASIS=6-31G*
NMUL=1 CHARGE=0
DFT=B3LYP
$END
$GEO
C -0.70374300 -0.00000500 0.01799100
H -1.12298600 -0.88527300 -0.48881000
H -1.12282100 0.88596400 -0.48770400
H -1.08818100 -0.00061600 1.06100200
N 0.75325100 -0.00000600 -0.13004000
H 1.14182300 0.81010400 0.35892800
H 1.14186200 -0.81010200 0.35892100
$END
$EDA
NMOL = 2
MATOM = 4 3
MMULT = 2 -2
MCHARGE = 0 0
$END
Tip
The molecule is singlet and separated into two neutral groups CH3 and NH2, so the multiplicities of momoners should be specified as MMULT=2 -2. The number 2 means one extra \(\alpha\) spin electron in the first monomer , and -2 means the extra \(\beta\) spin electron in the second monomer.
2.1.2.2. Computational Results
\(\Delta\)Eele |
\(\Delta\)Eexrep |
\(\Delta\)Epol |
\(\Delta\)Ecorr/disp |
\(\Delta\)Etotal |
|---|---|---|---|---|
-169.68 |
345.89 |
-244.18 |
-28.97 |
-96.93 |
For covalent bonds, there is a large orbital overlap between the two monomers, so polarization and exchange-repulsion interactions are dominant.
Note
When the supramolecule is a closed-shell molecule, the user does not need to change the METHOD to ROHF or UHF, and the program automatically performs the ROHF calculation when calculating the open-shell monomer. See $CTR in manual for details.
Please try recalculating using the ROHF and UHF keywords. Then explain the obtained calculation results according to the descriptions provided in the keywords.
2.1.3. Ionic bond in NaCl
2.1.3.1. Input File
$CTR
METHOD=RHF BASIS=6-31+G*
NMUL=1 CHARGE=0
DFT=B3LYP
$END
$GEO
Na -1.2465865600 0.4204757279 -6.0207462300
Cl -1.2465865600 2.9489348721 -6.0207462300
$END
$EDA
NMOL = 2
MATOM = 1 1
MMULT = 1 1
MCHARGE = 1 -1
$END
Tip
For ionic bond (Na+Cl-), MCHARGE=1 -1 and MMULT=1 1.
2.1.3.2. Computational Results
\(\Delta\)Eele |
\(\Delta\)Eexrep |
\(\Delta\)Epol |
\(\Delta\)Ecorr/disp |
\(\Delta\)Etotal |
|---|---|---|---|---|
-133.92 |
18.57 |
-10.33 |
-3.86 |
-129.53 |
For interactions between ionic compounds, the contribution of electrostatic interactions far exceeds that of other interactions, which is in accordance with our physical intuition.
2.1.3.3. Question
Is there a significant difference in the results obtained by replacing the functinal with PBE0 ?
2.1.4. \(\pi\)-\(\pi\) interaction
2.1.4.1. Input File
$CTR
METHOD=RHF BASIS=6-31+G*
NMUL=1 CHARGE=0
DFT=M06-2X
$END
$GEO
C -2.0787228974 2.7229635393 6.6336036675
C -2.3384441639 3.0120408790 5.2804066468
C -3.4306696328 2.4084131369 4.6291749443
C -4.2610379401 1.5132632818 5.3297470057
C -3.9991720354 1.2207967832 6.6815518176
C -2.9090520542 1.8273866261 7.3341751219
H -1.2288553785 3.1899751017 7.1383108766
H -1.6894330874 3.7015094040 4.7343839033
H -3.6306903605 2.6320648397 3.5777976439
H -5.1071486402 1.0412059601 4.8230040132
H -4.6393490139 0.5196893667 7.2231652503
H -2.7038451909 1.5979210774 8.3831631954
C -1.3788405600 -0.4314556117 4.8424264217
C -2.0431453313 -0.8230204965 3.6654369688
C -1.5413768555 -1.8846113993 2.8896480985
C -0.3710084116 -2.5563209868 3.2910801781
C 0.2950727691 -2.1663503674 4.4683054357
C -0.2096426847 -1.1043713410 5.2424621056
H -2.0585487502 -2.1883410484 1.9752620328
H 0.0200082528 -3.3814097325 2.6894789026
H 1.2041577486 -2.6882912577 4.7795214837
H 0.3061860493 -0.8005718704 6.1577900894
H -1.7716100935 0.3921427230 5.4395136059
H -2.9512037353 -0.2977462047 3.3583626867
$END
$EDA
NMOL = 2
MATOM = 12 12
MMULT = 1 1
MCHARGE = 0 0
$END
2.1.4.2. Computational Results
\(\Delta\)Eele |
\(\Delta\)Eexrep |
\(\Delta\)Epol |
\(\Delta\)Ecorr/disp |
\(\Delta\)Etotal |
|---|---|---|---|---|
-3.57 |
8.73 |
-1.26 |
-5.98 |
-2.08 |
2.1.4.3. Question
Please replace the functional with B3LYP and recalculate this example. Do you think the results are reasonable? Why?
2.2. Exercises
2.2.1. Energy Decomposition Analysis for Menshutkin Reaction
2.2.1.1. Introduction
In this exercise, the Menshutkin reaction NH3 + CH3Cl \(\rightarrow\) [NH3CH3]+ + Cl-will be analyzed with GKS-EDA in XEDA. The reaction is exothermic in solution. In this example, the reaction barrier in gas phase will be calculated with L-VBSCF and L-BOVB. For more information of the reaction and the computations in solution, please refer to:
The computations for reactant and TS are proceeded in the B3LYP-D3(BJ)/6-31G* level.
2.2.1.2. Input Files
Input file for reactant
$CTR
METHOD=RHF BASIS=6-31G*
NMUL=1 CHARGE=0
DFT=B3LYP-D3BJ
$END
$GEO
N 0.000000 0.000000 0.000000
H 0.939678 0.000000 -0.389227
H -0.469839 0.813785 -0.389227
H -0.469839 -0.813785 -0.389227
C 0.000000 0.000000 10.000000
H -1.029985 0.000000 9.647297
H 0.514992 0.891993 9.647297
H 0.514992 -0.891993 9.647297
Cl 0.000000 0.000000 11.778400
$END
$EDA
NMOL = 2
MATOM = 4 5
MMULT = 1 1
MCHARGE = 0 0
$END
Input file for TS
$CTR
METHOD=RHF BASIS=6-31G*
NMUL=1 CHARGE=0
DFT=B3LYP-D3BJ
$END
$GEO
N 0.0000000000 0.0000000000 -2.4403680000
H -0.9550790000 0.0000000000 -2.8031240000
H 0.4775395000 0.8271226766 -2.8031240000
H 0.4775395000 -0.8271226766 -2.8031240000
C 0.0000000000 0.0000000000 -0.6327610000
H 1.0651490000 0.0000000000 -0.4753840000
H -0.5325745000 -0.9224460928 -0.4753840000
H -0.5325745000 0.9224460928 -0.4753840000
Cl 0.0000000000 0.0000000000 1.8067450000
$END
$EDA
NMOL = 2
MATOM = 4 5
MMULT = 1 1
MCHARGE = 0 0
$END
2.2.1.3. Computational Results
After the job is finished, you can see the total interaction energy and different energy components in the output file:
Computational Results for Reactant:
--------------------------------------------------------------------------------
ALL BASIS SET HARTREE KCAL/MOL
--------------------------------------------------------------------------------
ELECTROSTATIC ENERGY ES= -0.000157 -0.10
EXCHANGE ENERGY EX= -0.000000 -0.00
REPULSION ENERGY REP= 0.000000 0.00
POLARIZATION ENERGY POL= -0.000024 -0.02
ELECTRON CORRELATION CORR= 0.000024 0.02
GRIMME DISP CORRECTION DC= -0.000003 -0.00
TOTAL INTERACTION ENERGY E= -0.000160 -0.10
--------------------------------------------------------------------------------
Computational Results for TS:
--------------------------------------------------------------------------------
ALL BASIS SET HARTREE KCAL/MOL
--------------------------------------------------------------------------------
ELECTROSTATIC ENERGY ES= -0.123292 -77.37
EXCHANGE ENERGY EX= -0.242129 -151.94
REPULSION ENERGY REP= 0.473559 297.16
POLARIZATION ENERGY POL= -0.129126 -81.03
ELECTRON CORRELATION CORR= 0.002253 1.41
GRIMME DISP CORRECTION DC= -0.002933 -1.84
TOTAL INTERACTION ENERGY E= -0.021668 -13.60
--------------------------------------------------------------------------------
The user may also remove the $EDA section in input file and use XEDA to calculate the single point energies of NH3 and CHCl3 in the reactants and transition states separately.
Then the user may get the decomposition of reaction energy barrier with expression
where \(\Delta E_\textrm{TS-R}^\textrm{geo}\) is the geometric relaxation energy.
After all computations are completed, the EDA results are obtained:
\(\Delta E^\textrm{geo}\) |
\(\Delta\Delta E^\textrm{ele}\) |
\(\Delta\Delta E^\textrm{exrep}\) |
\(\Delta\Delta E^\textrm{pol}\) |
\(\Delta\Delta E^\textrm{corr/disp}\) |
\(\Delta\Delta E^\textrm{total}\) |
|---|---|---|---|---|---|
43.4 |
-77.3 |
145.2 |
-81.0 |
-0.4 |
29.9 |
From the results, it can be seen that both repulsive and geometric relaxation energies contribute to the energy barrier, while polarization interactions play a significant role in the stability of the transition state.
2.2.1.4. Question
For further investigation, the user may divide the supermolecule into three monomers NH3, CH:sub:3\(\cdot\) and \(\cdot\)Cl and calculate the changes in the three-body interaction during this process. What conclusions would you get?
There is a precursor in the Menshutkin reaction in which NH: sub:3 get closer to CH3Cl than in reactant. The geometry is given below. Try to run XEDA and find the difference between the precursor and reactant in EDA results.
N 0.000000 0.000000 -3.256735 H 0.000000 -0.931171 -3.628844 H 0.806418 0.465586 -3.628844 H -0.806418 0.465586 -3.628844 C 0.000000 0.000000 0.163425 H -0.000000 1.020365 -0.178685 H 0.883662 -0.510183 -0.178685 H -0.883662 -0.510183 -0.178685 CL 0.000000 0.000000 1.955246
2.2.2. Cation-pi interaction
In this exercise, we have prepared the input file, you just need to submit the work, look at the output. The input file as follows
2.2.2.1. Input file
$CTR
METHOD=RHF BASIS=6-31G
NMUL=1 CHARGE=1
DFT=PBE
$END
$GEO
C 0.21878885 -1.05605662 0.98621428
C -0.71192259 1.08998501 2.54244971
C 1.13066745 -0.15188175 1.55317605
C -1.16105485 -0.88495368 1.19366765
C -1.62583232 0.18710911 1.97563076
C 0.66559380 0.91605544 2.33674526
H 2.19866991 -0.30624300 1.42366076
H -1.86558342 -1.60525036 0.78568548
H -2.69005466 0.29701740 2.16732907
H 1.37428355 1.59029782 2.81085777
H -1.06942880 1.89622796 3.17752671
H 0.57966918 -1.90618217 0.41317517
H -0.12747069 1.28426826 -1.57669163
N -0.48181823 1.64109278 -0.68697900
H -1.29492235 2.23728108 -0.85431081
H 0.24909009 2.17709398 -0.21351482
H -0.74412972 0.84957516 -0.06764305
$END
$EDA
NMOL = 2
MATOM = 12 5
MMULT = 1 1
MCHARGE = 0 1
$END
2.2.2.2. Question
Do the calculation again with BLYP and PBE functionals and compare the two results in the table below.
\(\Delta\)Eele |
\(\Delta\)Eexrep |
\(\Delta\)Epol |
\(\Delta\)Ecorr/disp |
\(\Delta\)Etotal |
|
|---|---|---|---|---|---|
PBE |
|||||
BLYP |
2.2.3. Water tetramer
This is four-body interaction, on the basis of hints and coordinate, calculate the four-body interaction. The basis set is cc-pVDZ and DFT functional is B3LYP-D3BJ.
In this exercise, the input file should be prepared by the users and the geometry will be provided.
2.2.3.1. Geometry
O 1.2245363506 1.9369107829 -1.5249005020
H 0.6059851978 2.5732298867 -1.8865861643
H 1.1674934363 2.0403343903 -0.5520881596
O 1.4426170965 2.2936999659 1.1607557179
H 2.4102571835 2.4430589600 1.2006505497
H 1.2658427477 1.5917744265 1.7885628634
O 4.1276499342 2.6315069165 0.8942820200
H 4.5631857354 3.4690433373 1.0588206123
H 4.1895885441 2.4896149845 -0.0733228603
O 3.9157890062 2.2134800316 -1.7828881536
H 2.9430129055 2.1035712740 -1.8276494239
H 4.2809601224 1.4397383938 -2.2142395894
Tip
NMOL = 4
MATOM =3 3 3 3
MCHARGE =0 0 0 0
MMULT =1 1 1 1
2.2.3.2. Question
What is the difference between interactions in the water dimer and tetramer?
How many two-body interactions are there? Please calculate all the two-body interactions.
How many three-body interactions are there? Please calculate all the three-body interactions.
Subtract all two-body interactions from the results of the four-body interactions and what do you get?
2.2.3.3. Reference
2.2.4. Multiple bond
This is a multiple bond, on the basis of hints and coordinate, calculate the multiple bond. The user should prepare the input file themselves with the provided geometry.
The basis set is cc-pVDZ and the DFT functional is B3LYP.
N 0.00000000 0.00000000 0.55222303
N 0.00000000 0.00000000 -0.55222303
Tip
METHOD = RHF
NMOL = 2
MATOM =1 1
MCHARGE =0 0
MMULT =4 -4
2.2.4.1. Question
How to perform EDA calculation of CO molecule?
2.2.5. Using initial guess
Using the prepared vec file, calculate the interaction between benzene and water.
The basis set is 6-31G and DFT functional is B3LYP-D3BJ.
Tip
MATOM=12 3
MCHARGE=0 0
MMULT=1 1
$VEC1
$VEC2
2.2.5.1. Geometry
C -2.22442746 3.10841370 6.79274511
C -2.58046699 3.45113039 5.48894310
C -3.75658441 2.94959259 4.93307257
C -4.57816410 2.10676932 5.68095732
C -4.22316265 1.76416957 6.98474979
C -3.04587722 2.26470447 7.54093933
H -1.31453454 3.50149727 7.22691774
H -1.94533265 4.10915327 4.91047335
H -4.03254080 3.21535087 3.92073798
H -5.49220753 1.71919358 5.25018549
H -4.86301088 1.11526763 7.56816530
H -2.77135658 2.00123024 8.55407429
O -4.88698006 5.16577339 7.71546459
H -5.31795645 5.35227680 6.87650633
H -4.34728813 4.38893747 7.52339935
2.2.6. Base Pair (Guanine\(\cdots\)Cytosine)
This is an interaction between guanine and cytosine, on the basis of hints and coordinate, calculate the interaction of base pairs. The basis set is aug-cc-pVDZ and DFT functional is PBE.
2.2.6.1. Geometry
C -1.60545400 -4.76004200 0.00000000
N -0.43338800 -4.19010400 0.00000000
C -0.70709500 -2.83875700 0.00000000
C 0.15615900 -1.70395500 0.00000000
O 1.39343300 -1.67109900 0.00000000
N -0.57071300 -0.50869100 0.00000000
C -1.93515300 -0.39634200 0.00000000
N -2.43756800 0.85011000 0.00000000
N -2.74962500 -1.43773200 0.00000000
C -2.08429500 -2.60528300 0.00000000
H -1.78573000 -5.82722700 0.00000000
H 0.00000000 0.35440600 0.00000000
H -3.44059800 0.93906600 0.00000000
H -1.85217900 1.69009000 0.00000000
N -2.64385900 -3.85092400 0.00000000
H -3.63383600 -4.04644300 0.00000000
N 1.25536600 4.22216600 0.00000000
C 0.46379200 3.06706400 0.00000000
O -0.76263000 3.19975400 0.00000000
N 1.09991300 1.87493400 0.00000000
C 2.43342300 1.80377400 0.00000000
N 2.98731800 0.59422900 0.00000000
C 3.25406300 2.98666000 0.00000000
C 2.60993100 4.17793500 0.00000000
H 3.12707000 5.13194300 0.00000000
H 4.33581900 2.93065600 0.00000000
H 2.39505600 -0.25873800 0.00000000
H 3.98980100 0.49292400 0.00000000
H 0.75384100 5.09983900 0.00000000
Tip
NMOL = 2
MATOM = 16 13
MCHARGE = 0 0
MMULT =1 1
2.2.6.2. Question
Trying to use the METHOD=FAST in the $EDA section to compare the time consumption with the default case.
2.2.7. Non-covalent interactions in S66
Using the coordinates of S66 provided in the zip file, calculate the interactions with the DFT functional and basis set you prefer.