.. tip:: All input files can be downloaded: :download:`Files `. .. tip:: For more information of this section, please refer to these keywords: - :doc:`../keywords/scf` - :doc:`../keywords/scfguess` TSO-DFT (2): Diabatic States ============================================================== .. contents:: :local: This tutorial will lead you step by step to do target state optimization DFT (TSO-DFT) using Qbics. TSO-DFT is an originally developed method in Qbics. It is a powerful and accurate method for **excited** and **diabatic** states. In many physical processes like charge transfer, **TSO-DFT** can give **much more reasonable** results than TDDFT. Therefore, it is worth learning this method. .. hint:: If you use Qbics to do TSO-DFT in you paper, please cite the following references: - `J. Chem. Theory Comput. 2023, 19, 1777 `_ In this tutorial, we will introduce how to use TSO-DFT to study **dibatic states.** Example: Definition of Diabatic States of NaCl -------------------------------------------------- What is a diabatic state? Its exact definition can be found in textbooks but are not relevant to us here. Here we just use an example to show what a diabatic state is. We will take NaCl **molecule** as an example. Adiabatic State ^^^^^^^^^^^^^^^^^^ We will calculate the ground state of NaCl at B3LYP/def2-TZVP level of theory: .. code-block:: bash :linenos: :caption: nacl-1.inp basis def2-tzvp end scf charge 0 spin2p1 1 type U end mol Na 0.0 0.0 -0.68381052 Cl 0.0 0.0 1.68381052 end task energy b3lyp end Let us check the Mulliken charges in ``nacl-1.out``: .. code-block:: bash :linenos: :caption: nacl-1.out Mulliken Populations ==================== # Symbol Charge Spin ---------------------------------------------- 1 Na 0.64773687 -0.00000889 2 Cl -0.64773687 0.00000889 ---------------------------------------------- Sum -0.00000000 0.00000000 ---------------------------------------------- ...omitted... Final total energy: -622.60503718 Hartree Therefore, the Mulliken charges are ``0.6477`` and ``-0.6477``, respectively, so the charges are delocalzied over the whole molecule. This state is also called **adiabatic state.** Diabatic State: [Na\ :sup:`+`][Cl\ :sup:`-`] ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Intuitively, one may think that NaCl is composed of a sodium cation Na\ :sup:`+` and a chloride antion Cl\ :sup:`-`. However, as shown above, in the adiabatic state, charges are delocalized over the molecule. Is it possible to calculate a state that charges are forced to localize on some atoms? Such state is called **diabatic state**. In the following, we will calculate a diabatic state of [Na\ :sup:`+`][Cl\ :sup:`-`]: .. code-block:: bash :linenos: :caption: nacl-2.inp basis def2-tzvp end scf charge 0 spin2p1 1 type U # This is needed for TSO-DFT. no_scf tso end scfguess type fragden frag +1 1 1 frag -1 1 2 end mol Na 0.0 0.0 -0.68381052 Cl 0.0 0.0 1.68381052 end task energy b3lyp end You can see that, in ``scfguess...end``, we are setting a superposition of fragment density (``fragden``) as SCF initial guess. The only difference is that we set ``no_scf tso`` in ``scf...end``: - ``frag +1 1 1`` We are seeting the fragment of atom ``1`` with charge ``+1`` and spin multiplicity ``1``. See :doc:`../keywords/scfguess` for details. - ``frag -1 1 2`` We are seeting the fragment of atom ``2`` with charge ``-1`` and spin multiplicity ``1``. After calculation, in ``nacl-2.out``, we can find this: .. code-block:: bash :linenos: :caption: nacl-2.out Mulliken Populations ==================== # Symbol Charge Spin ---------------------------------------------- 1 Na 1.00000000 0.00000000 2 Cl -1.00000000 -0.00000000 ---------------------------------------------- Sum -0.00000000 -0.00000000 ---------------------------------------------- ...omitted... Final total energy: -622.58940107 Hartree Indeed, at Na and Cl atom, the Mulliken charge is ``1.00`` and ``-1.00``, respectively, in line with our assumed diabatic state [Na\ :sup:`+`][Cl\ :sup:`-`]. Its is higher than the adiabatic state by -622.58940107 - (-622.60503718) = 0.0156 Hartree = 9.81 kcal/mol Actually, this is the so-called **charge-transfer interaction energy.** By the way, if you check the SCF process of this diabatic state: .. code-block:: bash :linenos: :caption: nacl-2.out #it SCF energy 2e energy |deltaE| |deltaP| Orbital Time/second --------------------------------------------------------------------------------------------- 1 -622.57645781 318.90071152 1.950E-01 1.100E-05 guess 0.046 2 -622.58412894 319.42032296 7.671E-03 3.461E-05 diis 0.054 3 -622.58935124 319.21103043 5.222E-03 3.376E-05 diis 0.052 4 -622.58939707 319.19096764 4.584E-05 5.708E-06 diis 0.052 5 -622.58939958 319.18902943 2.506E-06 9.746E-07 diis 0.051 6 -622.58940082 319.19159833 1.235E-06 3.969E-07 diis 0.053 7 -622.58940107 319.19461843 2.492E-07 9.841E-08 diis 0.052 --------------------------------------------------------------------------------------------- Great (^ _ ^) SCF has converged in 7 iterations. (0.367 seconds) The energy difference of the first and last cycle -622.58940107 - (-622.57645781) = 0.0129 Hartree = 8.12 kcal/mol is the so-called **polarization interaction enerrgy**. Energy Decomposition Analysis of [Na\ :sup:`+`][Cl\ :sup:`-`] ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ To check this, if you do an EDA (see :doc:`eda`) for [Na\ :sup:`+`][Cl\ :sup:`-`] with the following input: .. code-block:: bash :linenos: :caption: nacl-eda.inp basis def2-tzvp end scf charge 0 spin2p1 1 type U # This is needed for TSO-DFT. end eda type tso frag +1 1 1 frag -1 1 2 end mol Na 0.0 0.0 -0.68381052 Cl 0.0 0.0 1.68381052 end task eda b3lyp end In the output ``nacl-eda.out``: .. code-block:: bash :linenos: :caption: nacl-eda.out WITHOUT BSSE correction: Electrostatic interaction energy: -146.30 kcal/mol Exchange-correlation interaction energy: 23.91 kcal/mol Polarization interaction energy: -8.12 kcal/mol Charge transfer interaction energy: -9.81 kcal/mol Grimme's dispersion interaction: 0.00 kcal/mol ---------------------------------------------------------------- Total interaction energy: -140.33 kcal/mol Here, ``Charge transfer interaction energy: -9.81 kcal/mol`` and ``Polarization interaction energy: -8.12 kcal/mol`` confirm our previous calculations. Diabatic State: [Na][Cl] ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Besides [Na\ :sup:`+`][Cl\ :sup:`-`], we can also assume another diabatic state: NaCl is composed of 2 neutral radical: .. code-block:: bash :linenos: :caption: nacl-3.inp basis def2-tzvp end scf charge 0 spin2p1 1 type U # This is needed for TSO-DFT. no_scf tso end scfguess type fragden frag 0 2 1 frag 0 -2 2 end mol Na 0.0 0.0 -0.68381052 Cl 0.0 0.0 1.68381052 end task energy b3lyp end Here, ``frag 0 2 1`` means that Na has a single alpha electron, while ``frag 0 2 1`` means that Cl has a single **beta** electron, like this: .. figure:: figs/a30.jpg After calculation, we can find this in ``nacl-3.out``: .. code-block:: bash :linenos: :caption: nacl-3.out --------------------------------------------------------------------------------------------- #it SCF energy 2e energy |deltaE| |deltaP| Orbital Time/second --------------------------------------------------------------------------------------------- 1 -622.42824748 316.12630594 3.002E-02 3.077E-04 guess 0.051 2 -622.43470909 315.62361860 6.462E-03 2.354E-04 diis 0.063 3 -622.43589868 315.69293084 1.190E-03 2.159E-05 diis 0.057 4 -622.43598379 315.69186515 8.511E-05 2.740E-05 diis 0.059 5 -622.43599795 315.67167823 1.417E-05 5.266E-07 diis 0.057 6 -622.43600008 315.68068350 2.120E-06 9.311E-06 diis 0.071 7 -622.43600052 315.68009817 4.407E-07 5.151E-07 diis 0.057 --------------------------------------------------------------------------------------------- ...omitted... Mulliken Populations ==================== # Symbol Charge Spin ---------------------------------------------- 1 Na -0.00000000 0.50000000 2 Cl 0.00000000 -0.50000000 ---------------------------------------------- Sum -0.00000000 -0.00000000 ---------------------------------------------- ...omitted.. Final total energy: -622.43600052 Hartree So, the polarization and charge transfer interaction energy is: -622.43600052 - (-622.42824748) Hartree = -4.88 kcal/mol -622.60503718 - (-622.43600052) Hartree = -106.07 kcal/mol Compared with [Na\ :sup:`+`][Cl\ :sup:`-`], the charge transfer energy of [Na][Cl] is extremely high. You can also try to compute a diabatic state of [Na\ :sup:`-`][Cl\ :sup:`+`] if you wish. In practise, to study diabatic states of molecular dimer or clusters, you can simply use :doc:`../keywords/eda` since all needed calculations will be done automatically. You just need to set the fragment charges and spin multiplicities. Example: Reactant and Product Diabatic States of F\ :sup:`-`\ +CH\ :sub:`3`\ Cl ----------------------------------------------------------------------------------------- Reactant Complex ^^^^^^^^^^^^^^^^^^^ Now we consider the reaction F\ :sup:`-`\ +CH\ :sub:`3`\ Cl. First, we calculate the reactant complex. Now we calculate the **reactant** and **product** diabatic state: .. tabs:: .. tab:: fch3cl-reac-reac.inp .. code-block:: bash :linenos: :caption: fch3cl-reac-reac.inp basis def2-svp end scf charge 0 spin2p1 -1 type U # This is needed for TSO-DFT. no_scf tso end scfguess type fragden frag 0 1 1-5 frag -1 1 6 end mol C -0.43654823 1.13197968 0.00000000 H -0.07987539 1.63637787 0.87365150 H -0.07987539 1.63637787 -0.87365150 H -1.50654823 1.13199286 0.00000000 Cl 0.15009830 -0.52737135 0.00000000 F -1.63124586 2.75454539 -0.44024554 end task energy b3lyp end .. tab:: fch3cl-reac-prod.inp .. code-block:: bash :linenos: :caption: fch3cl-reac-prod.inp basis def2-svp end scf charge 0 spin2p1 -1 type U # This is needed for TSO-DFT. no_scf tso end scfguess type fragden frag 0 1 1-4 6 frag -1 1 5 end mol C -0.43654823 1.13197968 0.00000000 H -0.07987539 1.63637787 0.87365150 H -0.07987539 1.63637787 -0.87365150 H -1.50654823 1.13199286 0.00000000 Cl 0.15009830 -0.52737135 0.00000000 F -1.63124586 2.75454539 -0.44024554 end task energy b3lyp end The two dibatic states can be set using ``frag``: .. figure:: figs/a31.jpg Transition State ^^^^^^^^^^^^^^^^^^^ In the same way, we calculate the reactant and product diabatic states, but at the transition state: .. tabs:: .. tab:: fch3cl-ts-reac.inp .. code-block:: bash :linenos: :caption: fch3cl-ts-reac.inp basis def2-svp end scf charge 0 spin2p1 -1 type U # This is needed for TSO-DFT. no_scf tso end scfguess type fragden frag 0 1 1-5 frag -1 1 6 end mol C -0.59663490 1.29366556 -0.07171644 H -0.18159761 1.68969864 0.84430018 H -0.12252870 1.52677601 -1.01427388 H -1.59298434 0.87532545 -0.06612308 Cl 0.43269073 -0.72181312 0.13298792 F -1.52293997 3.10024978 -0.26542025 end task energy b3lyp end .. tab:: fch3cl-ts-prod.inp .. code-block:: bash :linenos: :caption: fch3cl-ts-prod.inp basis def2-svp end scf charge 0 spin2p1 -1 type U # This is needed for TSO-DFT. no_scf tso end scfguess type fragden frag 0 1 1-4 6 frag -1 1 5 end mol C -0.59663490 1.29366556 -0.07171644 H -0.18159761 1.68969864 0.84430018 H -0.12252870 1.52677601 -1.01427388 H -1.59298434 0.87532545 -0.06612308 Cl 0.43269073 -0.72181312 0.13298792 F -1.52293997 3.10024978 -0.26542025 end task energy b3lyp end .. figure:: figs/a32.jpg After 4 calculations, we can get the following energies in ``fch3cl-reac-reac/prod.out`` and ``fch3cl-ts-reac/prod.out``: .. list-table:: :stub-columns: 1 * - Structure - Reactant Complex - Transition State * - Reactant Dibatic State - -599.62530856 - -599.62955746 * - Product Dibatic State - -599.49007047 - -599.62959766 Interestingly, we observe that the reactant and product diabatic energy are the same for transition state. Acutally, this is the basis for dibatic MECP method for transition state search. See :doc:`ts` and :doc:`../keywords/mecp` for more information.