.. tip:: All input files can be downloaded: :download:`Files `. .. tip:: For more information of this section, please refer to these keywords: - :doc:`../keywords/eda` Energy Decomposition Analysis ========================================= .. contents:: :local: This tutorial will lead you step by step to do energy decomposition analysis (EDA) using Qbics. Currently, EDA can only be used with DFT. .. hint:: If you use Qbics to do EDA in you paper, please cite the following references: - `J. Chem. Theory Comput. 2023, 19, 1777 `_ - `Phys. Chem. Chem. Phys. 2024, 26, 17549 `_ Example: EDA for (H\ :sub:`2`\ SO\ :sub:`4`)(NH\ :sub:`3`) ------------------------------------------------------------- Here, we will do an EDA for the molecular dimer (H\ :sub:`2`\SO\ :sub:`4`)(NH\ :sub:`3`). The input is given below: .. code-block:: :caption: 1s1n.inp :linenos: basis def2-svp end scf charge 0 spin2p1 1 type U # This is needed for EDA. end grimmedisp type bj # Options: none, zero, bj. end eda type tso frag 0 1 1-7 frag 0 1 8-11 end mol S -0.62161000 -0.11181000 0.09167200 O 0.03138700 -0.14316200 1.35576500 O -1.82120800 -0.82692400 -0.15912300 O 0.38678100 -0.46814500 -1.02324300 O -0.91128600 1.43385100 -0.17168300 H 1.34667100 -0.25497400 -0.71202200 H -1.58729700 1.50230600 -0.85583700 N 2.78190500 0.03915200 -0.05091700 H 3.20716100 0.92239700 -0.30232700 H 3.46881800 -0.69082400 -0.19060100 H 2.55167500 0.07103700 0.93671700 end task eda b3lyp end - Line 8: For EDA, you must explicitly set the ``type`` to ``U`` in ``scf...end``. - Line 16: Do TSO-type EDA. You can also use ``gks`` to do generalized Kohn-Sham (GKS) type EDA. - Line 17-18: The fragments in EDA. Each line defines the charge, spin multiplicity, and atom indices. See below: .. figure:: figs/a25.jpg After running this calculation, we have the following output: .. code-block:: :caption: 1s1n.out :linenos: WITHOUT BSSE correction: Electrostatic interaction energy: -30.95 kcal/mol Exchange-correlation interaction energy: 27.55 kcal/mol Polarization interaction energy: -5.44 kcal/mol Charge transfer interaction energy: -15.25 kcal/mol Grimme's dispersion interaction: -2.22 kcal/mol ---------------------------------------------------------------- Total interaction energy: -26.31 kcal/mol WITH BSSE correction: Electrostatic interaction energy: -30.95 kcal/mol Exchange-correlation interaction energy: 27.55 kcal/mol Polarization interaction energy: -5.44 kcal/mol Charge transfer interaction energy: -10.44 kcal/mol Grimme's dispersion interaction: -2.22 kcal/mol ---------------------------------------------------------------- Total interaction energy: -21.50 kcal/mol In this output, the interaction energy of the following process H\ :sub:`2`\ SO\ :sub:`4` + NH\ :sub:`3` = (H\ :sub:`2`\ SO\ :sub:`4`)(NH\ :sub:`3`) is decomposed into several components. The interaction energy WITH and WITHOUT BSSE is ``-21.50 kcal/mol`` and ``-26.31 kcal/mol``. The BSSE energy is included in ``Charge transfer interaction energy``. The ``Electrostatic interaction energy``, ``Polarization interaction energy`` and ``Charge transfer interaction energy`` are all stablizing the dimer, the first contributing most. The ``Exchange-correlation interaction energy`` can also be recognized as Pauli exclusion energy, which destablizing the dimer. For different kinds of molecular dimer, the weights of different components are quiet different, indicating the nature of chemical interactions. Example: EDA for (Ca\ :sup:`2+`)(SO\ :sub:`4`\ :sup:`2-`)(H\ :sub:`2`\ O)\ :sub:`5` ---------------------------------------------------------------------------------------- The basic MB-EDA formular is: .. math:: \Delta E^{\text{int}} = \frac{1}{2!} \sum_{I_1 \neq I_2}^{N} \Delta E_{I_1 I_2}^{(2)} +\frac{1}{3!} \sum_{I_1 \neq I_2 \neq I_3}^{N} \Delta E_{I_1 I_2 I_3}^{(3)} + \cdots + \frac{1}{n!} \sum_{I_1 \neq \cdots \neq I_n}^{N} \Delta E_{I_1 \cdots I_n}^{(n)} + \cdots + \Delta E_{I_1 \cdots I_N}^{(N)} \equiv \sum_{n=2}^{N} \Delta E^{(n)} The terms higher than :math:`\Delta E^{(2)}` is the many-body interaction term. Usually the most important one is the three-body effect :math:`\Delta E^{(3)}`, the effects of which can be decomposed into three ones: - :math:`\Delta E^{(3)} < 0`: Indicate a **cooperative effect** of the monomers in a cluster. The many-body interaction is stabilizing the cluster. This is often seen in hydrogen bonding clusters, like water clusters. - :math:`\Delta E^{(3)} > 0`: Indicate an **anti-cooperative effect** of the monomers in a cluster. The many-body interaction is destabilizing the cluster. This is often seen in a cluster of charged species, like ionic liquid clusters. Also see below. - :math:`\Delta E^{(3)} \approx 0`: Indicate a **non-cooperative effect** of the monomers in a cluster. There is little many-body interaction in the cluster. This is often seen in a cluster of molecules without charges or hydrogen bonds. Each order can be decomposed into electrostatic, exchange, polarization, charge transfer, and dispersion terms: .. math:: \Delta E_X^{(n)} = \Delta E_X^{(n)\text{-el}} + \Delta E_X^{(n)\text{-ex/xc}} + \Delta E_X^{(n)\text{-pl}} + \Delta E_X^{(n)\text{-ct}} + \Delta E_X^{(n)\text{-disp}} Usually, electrostatic and exchange terms are highly additive, while polarization and charge transfer terms are non-additive. The dispersion term is always additive. Here we will do an EDA for the microsolvated cluster (Ca\ :sup:`2+`)(SO\ :sub:`4`\ :sup:`2-`)(H\ :sub:`2`\ O)\ :sub:`5`. Since there are more than 2 species, we will use many-body EDA (MB-EDA) to decompose the interaction energy components into contributions from different many-body level. Usually, contributions beyond three-body are very small, so we will truncate the computation at four-body. The input is given below: .. code-block:: :caption: caso4h2o5.inp :linenos: basis def2-svp end scf charge 0 spin2p1 1 type U # This is needed for EDA. end grimmedisp type bj end eda type mb_tso mb_level 4 # The many-body truncation level. frag +2 1 1 frag -2 1 2-6 frag 0 1 7-9 frag 0 1 10-12 frag 0 1 13-15 frag 0 1 16-18 frag 0 1 19-21 end mol Ca 1.98604937 1.24582501 1.89905383 S 0.69103489 3.45226578 0.98237460 O 2.08708072 3.44779165 1.50308184 O 0.24826049 2.05102931 0.73622996 O -0.21032250 4.08940001 1.98323165 O 0.63912085 4.22084216 -0.29304506 O 3.02175028 0.83377391 -0.05175855 H 3.13194669 0.09601115 -0.65158745 H 2.90415260 1.59241194 -0.62348488 O 2.51307423 2.92493798 -1.72944079 H 1.80535522 3.44397361 -1.34738364 H 2.91500306 3.50585612 -2.37536585 O 0.60914213 -0.51387571 1.58896733 H 0.07645358 -1.26917200 1.83796991 H 0.01571971 0.04244367 1.08439778 O 0.73702072 1.93958625 3.70546127 H 0.73393089 2.03509947 4.65787899 H 0.20818330 2.67203308 3.38909995 O 3.87639195 1.93942131 2.91131853 H 4.78820918 1.95286631 3.20224881 H 3.69756418 2.83995160 2.64058309 end task eda b3lyp end - Line 16: To do MB-EDA, you should use ``mb_tso``. - Line 17: ``mb_level 4`` means you truncate the MB-EDA at four-body. - Line 18-24: Set 7 chemical fragements. Note that some species like Ca\ :sup:`2+` and SO\ :sub:`4`\ :sup:`2-` are charged. You should give the correct charge. The atomic indices are shown below: .. figure:: figs/a26.jpg After running this calculation, we have the following output: .. code-block:: :caption: caso4h2o5.out :linenos: Table 5. Summary (kcal/mol). --------------------------------------------------------------------------------------------------------------------------------- Interactions delE_el delE_xc delE_pl delE_ct delE_bsse delE_disp delE_tot --------------------------------------------------------------------------------------------------------------------------------- SUM of 2-body -6.73148230E+02 1.47696370E+02 -1.11395126E+02 -2.04547116E+02 1.03684957E+02 -2.06082308E+01 -7.58317375E+02 SUM of 3-body 1.42501290E-08 -2.67700179E+00 6.09855572E+01 1.30716087E+02 -4.86421459E+01 1.07484123E+01 1.51130908E+02 SUM of 4-body -2.82504918E-08 5.22296932E-01 -4.66678754E+00 -4.78862109E+01 1.47634652E+01 -1.52639538E+01 -5.25311901E+01 Remain 1.63009297E-08 -9.27025377E-02 2.21610464E-01 1.02616142E+01 -2.89029982E+00 7.47122033E+00 1.49714426E+01 --------------------------------------------------------------------------------------------------------------------------------- SUM -6.73148230E+02 1.45448963E+02 -5.48547462E+01 -1.11455626E+02 6.69159769E+01 -1.76525519E+01 -6.44746214E+02 --------------------------------------------------------------------------------------------------------------------------------- We can see that the total interaction energy is -644.75 kcal/mol, which is decomposed into 2-body, 3-body, 4-body, and remaining terms. The 2-body term is the most important one (-758.32 kcal/mol), while the 3-body term is also significant but **anti-cooperative** (destablizing the cluster) (+151.13 kcal/mol). The 4-body term is small (-52.53 kcal/mol, slightly cooperative). The remaining term (sum of 5-, 6-, and 7-body) is very small (+14.97 kcal/mol) compared to the low-order terms. We can also see that the electrostatic and exchange energy are highly **additive** (beyond 3-body is almost zero), while the polarization and charge-transfer energy are **non-additive**: .. list-table:: :stub-columns: 1 * - Type - 2-body - 3-body - 4-body - Total * - Polarization - -111.40 kcal/mol - +60.99 kcal/mol - -4.67 kcal/mol - -54.85 kcal/mol * - Charge transfer - -204.55 kcal/mol - +130.72 kcal/mol - -47.89 kcal/mol - -111.45 kcal/mol For different kinds of clusters, the 3-body effects (many-body interactions) can be quite different, see `Phys. Chem. Chem. Phys. 2024, 26, 17549 `_ for more information.