1. Manual
1.1. Keywords for XEDA
All keywords are case-insensitive.
1.1.1. $CTR section
This group specifies the type of wavefunction, basis set and density functional theory.
1.1.1.1. METHOD=option
Computational method used in the calculation. Available methods are:
RHF: Restricted Hartree-Fock
UHF: Unrestricted Hartree-Fock
ROHF: Restricted Open-shell Hartree-Fock
Note
Basically, the same method will be used for both supermolecules and each monomer. For example, if METHOD=UHF is specified, UHF computations will be proceeded for both supermolecule and each monomer.
Specially, if the supermolecule is singlet while some monomers are open-shell molecules/fragments, METHOD=RHF will make the program run RHF computation for supermolecule but ROHF for monomers.
This keyword is essential even DFT is specified. See the DFT keyword for details.
1.1.1.2. BASIS=option
Specify the basis set for the computation. Currently supported basis set are listed below.
Pople’s basis sets (with * for polarization functions and + for diffusion functions):
3-21G
6-31G
6-311G
Dunning-type correction consistent basis sets (with prefix aug- for diffusion functions):
cc-pVDZ
cc-pVTZ
cc-pVQZ
def2 basis sets
def2-SVP
def2-SVPP
def2-TZVP
def2-TZVPP
def2-QZVP
def2-QZVPP
def2-QZVPD
1.1.1.3. CHARGE=n
Total charge of the supermolecule. Default is 0, meaning the neutral molecule or radical.
1.1.1.4. NMUL=n
Spin multiplicity of the supermolecule, 1 for singlet, 2 for doublet and so on. Default is 1.
1.1.1.5. DFT
The name of density functional (M06-2X, B3LYP(-D3/-D3BJ), cam-B3LYP, PBE0, PBE, BLYP, \(\omega\)B97X-D).
Keyword |
Functional |
|---|---|
BLYP |
BLYP |
B3LYP |
B3LYP |
B3LYP-D3 |
B3LYP with D3 dispersion correction |
B3LYP-D3BJ |
B3LYP with D3BJ dispersion correction |
CAM-B3LYP |
CAM-B3LYP |
PBE |
PBE |
PBE0 |
Hybrid PBE with 25% HF exchange |
M06-2X |
M06-2X |
wB97XD |
\(\omega\)B97X-D |
Note
If this keyword is absent, HF method will be used by default.
The DFT keyword is only used to specify the functional. Keyword METHOD=option should also be specified.
1.1.1.6. MAX_ITER=n
Maximum number of iterations. Default value is 50.
1.1.2. $GEO section
This group shows molecular cartesian coordinates. Only Cartesian format is supported so far. An example for water molecule is given below:
$GEO
O 0.452988940 -0.000000000 -0.000009816
H -0.126498394 0.000000000 0.748296646
H -0.126490546 0.000000000 -0.748286838
$END
1.1.3. $EDA section
This group defines the information of each monomer during EDA calculations. Currently up to 10 monomers can be supported.
Note
In this section, each line should contain only one keyword.
1.1.3.1. NMOL=n
Number of monomers.
1.1.3.2. MATOM=list
A list giving the number of atoms in each monomer. The ordering of atoms in $GEO section should be the same as the ordering of monomers. The sum of the MATOM list must be equal to the total number in the supermolecule. Numbers should be seperated by space.
1.1.3.3. MMULT=list
A list giving the multiplicity of each monomer. A positive integer means alpha spin, a negative integer means beta spin. Numbers should be seperated by space.
Tip
The summation of MMULT should be equal to the total multiplicity specified in NMUL=n.
For example, if an ethane molecule is separated into two neutral CH3 groups, MMULT should be specified as MMULT=2 -2 or MMULT=-2 2.
1.1.3.4. MCHARGE=list
A list giving the charge of each monomer. The sum of the charges in the monomers must be equal to the total charge of the supermolecule as specified in CHARGE=n. Numbers should be seperated by space.
1.1.3.5. METHOD=algorithm
Algorithms used for EDA evaluation. Available options are:
EXACT: Coulumb and exchange terms are computed directly. (default).
FAST: Coulumb and exchange terms are computed with RIJ-COSK algorithm.
Note
METHOD=FAST is faster for large system with huge basis sets.
For pure functional (such as PBE and BLYP), METHOD=FAST is highly recommended.
1.1.4. $VEC section
For the system that is difficult to converge, it is helpful to provide appropriate initial guess. The initial guess of supermolecule is given by $VEC0, and $VEC1, $VEC2, … for monomers with no BSSE correction.
Note
We can use the orbitals from one functional as the initial guess to other functional.
1.1.4.1. Bibliography
LMOEDA
GKSEDA
2. Tutorial
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.
3. Citation
You need to cite the references of XEDA as following formats:
A) J. Chem. Phys. format: The ab initio Valence Bond calculations are performed with the XEDA program.
P. Su, Z. Tang, W. Wu, WIREs Comput Mol Sci. 2020, 10:e1460;
Z. Tang, Y. Song, S. Zhang, W. Wang, Y. Xu, D.Wu, W. Wu, P. Su, J. Comput. Chem. 2021, 42, 2341.
B) American Chemical Society format: The ab initio Valence Bond calculations are performed with the XEDA program.