# 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

Computational results of NH3$$\cdots$$H2O (kcal/mol)

$$\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

Computational results of CH3-NH2 (kcal/mol)

$$\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 Computational results of Ionic bond in NaCl (kcal/mol) $$\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 Computational results of $$\pi$$-$$\pi$$ interaction (kcal/mol) $$\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

\begin{split}\begin{aligned} \Delta E &= E_\textrm{TS} - E_\textrm{R} \\ &= E_\textrm{TS}^\textrm{A} + E_\textrm{TS}^\textrm{B} + \Delta E_\textrm{TS}^\textrm{int} - E_\textrm{R}^\textrm{A} - E_\textrm{R}^\textrm{B} - \Delta E_\textrm{R}^\textrm{int}\\ &= \Delta \Delta E_\textrm{TS-R}^\textrm{int} + \left( E_\textrm{TS}^\textrm{A} - E_\textrm{R}^\textrm{A} \right) + \left( E_\textrm{TS}^\textrm{B} - E_\textrm{R}^\textrm{B} \right)\\ &= \Delta \Delta E_\textrm{TS-R}^\textrm{ele} + \Delta \Delta E_\textrm{TS-R}^\textrm{exrep} + \Delta \Delta E_\textrm{TS-R}^\textrm{pol} + \Delta \Delta E_\textrm{TS-R}^\textrm{corr/disp} + \Delta E_\textrm{TS-R}^\textrm{geo} \end{aligned}\end{split}

where $$\Delta E_\textrm{TS-R}^\textrm{geo}$$ is the geometric relaxation energy.

After all computations are completed, the EDA results are obtained:

Computational results of Menshutkin reaction with XEDA (in kcal/mol)

$$\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

1. 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?

2. 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.

Computational results with PBE and BLYP (kcal/mol)

$$\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

1. What is the difference between interactions in the water dimer and tetramer?

2. How many two-body interactions are there? Please calculate all the two-body interactions.

3. How many three-body interactions are there? Please calculate all the three-body interactions.

4. Subtract all two-body interactions from the results of the four-body interactions and what do you get?

#### 2.2.3.3. Reference

Su, P.; Li, H. Energy Decomposition Analysis of Covalent Bonds and Intermolecular Interactions. The Journal of Chemical Physics 2009, 131 (1), 014102.

### 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.