3.26 Extended Transition State - Natural Orbitals for Chemical Valence (ETS-NOCV) analysis¶
The Extended Transition State - Natural Orbitals for Chemical Valence (ETS-NOCV) method was proposed in J. Chem. Theory Comput., 5, 962 (2009) by Ziegler et al., it is a combination of the NOCV method proposed in J. Phys. Chem. A, 112, 1933 (2008) and the idea of ETS proposed in Theor. Chim. Acta, 46, 1 (1977). ETS-NOCV has been widely employed in studying chemical bonds and weak interactions between fragments. Briefly speaking, ETS-NOCV enables one to gain a deep insight into interfragment interaction by visual decomposition analysis of the variation of density matrix caused by the interaction.
I will first introduce the theory of ETS-NOCV as clearly and detailedly as possible in Section 3.26.1, and then some implementation details of this analysis in Multiwfn will be mentioned in Section 3.26.2. In Section 3.26.3 I will describe the functions and usage of the ETS-NOCV module in Multiwfn. Application examples of ETS-NOCV analysis are provided in Section 4.23.
3.26.1 Theory¶
The ETS-NOCV analysis in Multiwfn supports arbitrary number of fragments. However, for simplicity, only two fragments are involved when I introduce the theory of ETS-NOCV below.
Physical components of interfragment interaction energy¶
I first review the physical components of interfragment interaction energy, this is important preliminary knowledge for understanding ETS-NOCV.
Total energy variation due to interfragment interaction between fragment A and \(B\) may be expressed as (fragments mentioned below are in complex structure, so deformation energy due to distortion of fragment structure during combination is not taken into account in this context):
where \(E\)AB is the complex electronic energy, \(E\)A and \(E\)B are electronic energies of \(A\) and \(B\) at complex geometry. The four physical components are defined as follows
-
\(∆{E}_{\mathrm{els}}\): Interfragment electrostatic interaction energy. It is evaluated as classical Coulomb interaction energy between original wavefunctions of A and B (\({\Psi}_{A}\) and \({\Psi}_{B}\)).
-
\(∆{E}_{\mathrm{XC}}^{0}\): Change in exchange-correlation (XC) energy during combination from \({\Psi}_{A}\) and \({\Psi}_{B}\) to \({\Psi}_{A}{\Psi}_{B}\). It also accounts for interfragment dispersion interaction since dispersion effect is essentially interfragment Coulomb correlation.
-
\(∆{E}_{\mathrm{Pauli}}\): Energy increase due to Pauli repulsion between electrons of the two fragments. It is also known as exchange-repulsion term.
-
\(∆{E}_{\mathrm{orb}}\): Orbital interaction energy due to mix of fragment orbitals, which accounts for polarization effect (intrafragment mix of occupied and unoccupied orbitals) and charge-transfer effect (interfragment mix of occupied and unoccupied orbitals).
Specifically, the Δ\(E_{\mathrm{Pauli}}\) is expressed as
\[∆{E}_{\mathrm{Pauli}}=E\left[ {\Psi}_{AB}^{0} \right]-E\left[ {\Psi}_{A}{\Psi}_{B} \right]\]where \(E\)[ ] denotes the energy functional of the complex. \({\Psi}_{A}\) and \({\Psi}_{B}\) are original wavefunctions of \(A\) and \(B\). \({\Psi}_{A}{\Psi}_{B}\) may be referred to as promolecular wavefunction, it is the Hartree product of \({\Psi}_{A}\) and \({\Psi}_{B}\), and it simply corresponds to the Slater determinant directly constructed by all occupied orbitals of \(A\) and \(B\). \({\Psi}_{AB}^{0}\) is denoted as frozen state wavefunction, it is defined as antisymmetric product of \({\Psi}_{A}\) and \({\Psi}_{B}\); specifically, it corresponds to the Slater determinant constructed by all occupied orbitals of \(A\) and \(B\), and Löwdin orthonormalization has been employed among these orbitals to make them to be orthonormal with each other. The frozen state can be viewed as an artificial intermediate state comes from the fact that when two fragments are combined together to form an identical particle system, all electron orbitals must constitute an orthonormal set, or occupied orbitals of a fragment must be orthogonal to those of another fragment to comply Pauli exclusion principle. The interfragment Löwdin orthonormalization brings additional nodal plane(s) to original fragment orbitals, this phenomenon evidently makes kinetic energy of the orbitals increase, this is mainly why Δ\(E_{\mathrm{Pauli}}\) must be a positive value and plays a destabilization effect.
The Δ\(E_{\mathrm{orb}}\) is expressed as
\[∆{E}_{\mathrm{orb}}=E\left[ {\Psi}_{AB} \right]-E\left[ {\Psi}_{AB}^{0} \right]\]where \({\Psi}_{AB}\) is the actual complex wavefunction yielded after self-consistent field (SCF) convergence. The change from \({\Psi}_{AB}^{0}\) to \({\Psi}_{AB}\) arises from intrafragment and interfragment orbital mixing (also known as orbital relaxation). Δ\(E_{\mathrm{orb}}\) is always a negative term and thus stabilizes the complex. The ETS-NOCV method focuses on gaining deep chemical insights into the Δ\(E_{\mathrm{orb}}\) term.
In summary, aforementioned terms can be organized as the following relationship
where Δ\(E_{\mathrm{prep}}\) is known as preparation energy, it includes distortion energy of fragment A and \(B\) from their isolated geometries to complex geometry, and it also includes the energy change of their electronic states from most stable states to reference state ΨA and ΨB (for example, to use ETS-NOCV to study the double bond of H2Ge=GeH2, the reasonable fragment reference state should be triplet, however most stable state of GeH2 in its isolated status is singlet. This difference should be incorporated into the Δ\(E_{\mathrm{prep}}\)). Clearly the choice of reference states of the fragments affects result of ETS-NOCV analysis while it is somewhat arbitrary in certain cases.
Frankly speaking, in my viewpoint, the above commonly accepted partition of interaction energy is not completely rigorous. Because during transformation of complex wavefunction from promolecular (reference) state ΨAΨB to actual state ΨAB, the electrostatic interaction energy and exchange-correlation energy must also markedly change, therefore the Δ\(E_{\mathrm{orb}}\) term should not be regarded as solely reflecting the contribution to interaction energy due to orbital mix effect.
The \(∆{E}_{\mathrm{els}}+∆{E}_{\mathrm{XC}}^{0}+∆{E}_{\mathrm{Pauli}}\) is sometimes referred to as steric term Δ\(E_{\mathrm{steric}}\) in literature for convenience in discussion. Multiwfn is unable to directly evaluate it or its component, but you can calculate Δ\(E_{\mathrm{steric}}\) by Gaussian in combination with Multiwfn, see Section 3.100.8. The Δ\(E_{\mathrm{prep}}\) can be directly manually calculated via any quantum chemistry code.
NOCV theory¶
First, we look at the natural orbitals for chemical valence (NOCV) theory. The orbital interaction between the two fragments leads to a density matrix difference
where P and P0 are density matrices of actual complex state and frozen state, respectively; they can be easily generated based on coefficient matrix of occupied orbitals in corresponding state. The NOCV method diagonalizes the \(∆{\mathbf{P}}^{\mathrm{orb}}\) to solve its eigenvalues and eigenvectors, namely one has the following relationship (the matrices are expressed in Löwdin orthogonalized basis functions)
where CNOCV is the coefficient matrix of NOCV orbitals, each of its columns corresponds to expansion coefficients of a NOCV orbital with respect to Löwdin orthogonalized basis functions. v is a diagonal matrix, \(v\)i,i corresponds to eigenvalue of the \(i\)th NOCV orbital. One can also say NOCV orbitals are eigenfunctions of density matrix difference operator, namely
Note that the number of NOCV orbitals (\(N\)) is equal to the number of basis functions. Hence, usually the \(N\) is large, but only very few NOCV orbitals have notable magnitude of eigenvalues and which are what one should focus on during analysis.
Orbital interaction results in variation of electron density, which can be represented as “orbital deformation density”
The \(∆{\rho}^{\mathrm{orb}}\) can be decomposed as contributions of various NOCV densities {\({v}_{i}{\varphi}_{i}^{2}\)}
Because \(∆{\mathbf{P}}^{\mathrm{orb}}\) is a traceless matrix represented under a set of orthonormal basis, a noteworthy feature of NOCV orbitals is that they occur in pairs, that is if the eigenvalues are sorted from most positive to most negative, then \(v\)N+1-i = −\(v\)i. For simplicity, \(N+1\)−\(i\) will be abbreviated to −\(i\) later. So, the \(∆{\rho}^{\mathrm{orb}}\) can also be decomposed as NOCV pair contributions for easier analysis and discussion
Obviously, the NOCV analysis allows us to look into the details of the orbital deformation density in terms of NOCV orbitals or densities to better understand the nature of the orbital interaction.
ETS and ETS-NOCV theory¶
The extended transition state (ETS) theory shows that
where µ and ν are indices of Löwdin orthogonalized basis functions. FTS is the so-called extended transition state Fock matrix (or Kohn-Sham matrix in KS-DFT case), which is the Fock matrix constructed using average of \({\Psi}_{AB}\) and \({\Psi}_{AB}^{0}\). Note that the “extended transition state” in this context is very different to the transition state in common sense, in the present context it refers to the artificial electronic structure at midpoint between \({\Psi}_{AB}\) and \({\Psi}_{AB}^{0}\), where orbital interaction takes place only half.
ETS-NOCV theory shows that
where \(i\) is index of NOCV, and \({\tilde{F}}_{i,i}^{\mathrm{TS}}\) is the \(i\)th diagonal term of Fock matrix in the basis of NOCV orbitals; in other words, it corresponds to energy of the \(i\)th NOCV orbital estimated using \({\hat{F}}^{\mathrm{TS}}\) operator, namely \({\tilde{F}}_{i,i}^{\mathrm{TS}}=\left\langle {\varphi}_{i}\middle|{\hat{F}}^{\mathrm{TS}}\middle|{\varphi}_{i} \right\rangle\).
Again, because NOCV orbitals are paired, the \(∆{E}_{\mathrm{orb}}\) can be decomposed as contributions of NOCV pairs
From the energy contributions, we can determine which NOCV pairs play a major role in orbital interaction and then analyze their characters emphatically. The NOCV pairs with very small eigenvalues or energies can be ignored during discussion.
Deformation density¶
Finally, three kinds of deformation densities are summarized below, they are involved in ETS-NOCV analysis in Multiwfn, note that ρ, ρ0, ρA and ρB correspond to electron density of \({\Psi}_{AB}\), \({\Psi}_{AB}^{0}\), \({\Psi}_{A}\) and \({\Psi}_{B}\), respectively.
Pauli deformation density:
Orbital deformation density:
Total deformation density:
3.26.2 Implementation details¶
Features¶
The ETS-NOCV module in Multiwfn has the following capabilities:
-
Calculate NOCV orbital wavefunctions, eigenvalues and energies
-
Visualize NOCV orbitals and export corresponding cube files
-
Calculate energy contributions to Δ\(E_{\mathrm{orb}}\) of NOCV pairs
-
Visualize NOCV pair densities and export corresponding cube files
-
Calculate composition of NOCV pair and orbitals using SCPA method
-
Visualize promolecular, frozen state and actual complex orbitals
-
Visualize ΔρPauli, Δρorb, and Δρ isosurfaces
The ETS-NOCV module supports defining arbitrary number of fragments. If you defined \(M\) fragments, then the ETS-NOCV will analyze the total interaction between all the \(M\) fragments.
Only Hartree-Fock and Kohn-Sham DFT wavefunctions are supported. Multiconfiguration wavefunctions such as coupled-cluster, CASSCF, and double-hybrid functional wavefunctions are not supported.
Only restricted closed-shell or unrestricted open-shell wavefunctions are acceptable. Restricted open-shell wavefunction is not supported.
If complex is open-shell, or any fragment is open-shell, the ETS-NOCV analysis will be automatically conducted in open-shell form. In this case, the alpha and beta NOCV orbitals are solved independently, and their energies are estimated using alpha and beta Fock matrices respectively.
The energies of NOCV orbitals in Multiwfn are not calculated in the strict way of the standard ETS-NOCV method as described above! This is because the FTS in ETS-NOCV analysis currently is not available in Multiwfn. In the post-processing menu of ETS-NOCV module, you can choose to load a file containing the actual Fock matrix of the complex outputted by a quantum chemistry code, or you can also choose to let Multiwfn directly generate the actual Fock matrix of the complex based on the orbital energies and coefficients recorded in the input file via F=SCEC-1 relationship, then for example the \(i\)th NOCV orbital energy will be evaluated as \(\left\langle {\varphi}_{i}\middle|\hat{F}\middle|{\varphi}_{i} \right\rangle\), where \(\hat{F}\) is the Fock operator corresponding to the loaded or generated Fock matrix. According to my comparison with some published ETS-NOCV data and result of ORCA code, the NOCV energies evaluated in this approximated way using \(\hat{F}\) is close to the NOCV energies in strict sense derived based on \({\hat{F}}^{\mathrm{TS}}\) (especially for weak interaction case), at least this discrepancy does not qualitatively affect your identification of dominant NOCV orbitals/pairs. Notice that due to this difference, the sum of energies of all NOCV orbitals or pairs given by Multiwfn is not exactly equal to Δ\(E_{\mathrm{orb}}\).
Input files¶
You need to provide wavefunction files of the complex and that of each fragment. The file should contain basis function information, for example you can use .fch, .mwfn, .molden, and so on; however, .wfn and .wfx cannot be used since they do not contain basis function information. If you are confused, see Section 2.5.
To prepare the wavefunction files needed by ETS-NOCV analysis for a system, commonly you should do following steps
-
(1) Optimize the geometry of the complex first by your favourite quantum chemistry code, and meantime obtain wavefunction file of the complex
-
(2) Extract coordinate of each fragment from the optimized complex, and then save them as input files of single point task. If you are using Gaussian, do not forget to add nosymm keyword to avoid automatic reorientation during the calculation.
-
(3) Run input file of every fragment to obtain their wavefunction files.
Evidently, the calculation level of the fragments must be identical to that of complex.
If there is no special reason, using diffuse functions is not only completely unnecessary but also deprecated. Usually, using a 3-zeta basis set such as def2-TZVP is absolutely adequate, while using a 2-zeta basis set such as 6-31G* and def2-SVP is also sufficient for qualitative study.
After arranging the atoms in each fragment according to the loading order of the fragments in ETS-NOCV module, the order of the atoms must be consistent with that in the complex. Obviously, this requires that the order of atoms in fragment files must be the same as that in the complex, and the atoms of each fragment must occur consecutively in the complex.
Usually, using an implicit solvation model should be avoided, it makes situation much more complicated.
3.26.3 Usage¶
The general process of using the ETS-NOCV module is:
-
(1) Boot up Multiwfn, and load wavefunction file of the complex
-
(2) Enter main function 23
-
(3) Input the number of fragments
-
(4) Input file path of wavefunction file of each fragment. The loading order should be in line with the occurrence order of the fragments in the complex
Note
Note that if a fragment is found to be unrestricted open-shell, you will also be asked to choose if flipping spin. If you choose y, then the wavefunctions and occupancies of alpha and beta orbitals will be exchanged. You should properly set spin flipping so that the sum of number of alpha (beta) electrons of all fragments is identical to that of the complex.
(5) Now the NOCV orbital wavefunctions and eigenvalues are calculated, and NOCVs are automatically paired. From the NOCV information on screen, you can directly find NOCV pair indices and corresponding NOCV indices as well as their eigenvalues. Currently NOCV energies have not been evaluated.
After that a post-processing menu occurs, and you can choose corresponding option according to your requirement, they are briefly mentioned below.
If you want to view NOCV information again, choose “0 Print NOCV information”. By default, the printing threshold of NOCV eigenvalues is 0.001, the NOCVs with absolute eigenvalue less than this threshold are not shown to avoid too lengthy output. You can manually change the threshold via “-3 Set printing threshold of NOCV eigenvalues”.
To obtain NOCV energies, you need to select either “-1 Load Fock/KS matrix and evaluate NOCV orbital energies” to load a file containing Fock matrix of the current complex calculated at the same level (see Appendix 7 of this manual for detail), or select “-2 Generate Fock/KS matrix and evaluate NOCV orbital energies” to generate Fock matrix of current complex based on the energies and coefficients of molecular orbitals in the complex wavefunction file (i.e. the file loaded after booting up Multiwfn). Then Multiwfn will calculate NOCV energies, and then NOCV information with energies will be shown on screen.
“1 Show isosurface of NOCV orbitals” is used to visualize NOCV orbitals. After selecting this option, NOCV information will be shown on screen for user’s convenience, and a GUI window appears, you can select corresponding NOCV orbital from the list at the bottom right corner of the window, you can also directly input the index of the NOCV orbital that you want to view in the text box at the bottom right corner of the window.
“2 Show isosurface of NOCV pair density” is used to visually examine density of NOCV pair(s). After choosing it, NOCV information will be shown on screen for user’s convenience, and then you can input the NOCV pair index of interest. For example, if then you input 6, then grid data of \({v}_{6}\left[ {\varphi}_{6}^{2}\left( \mathbf{r} \right)-{\varphi}_{-6}^{2}\left( \mathbf{r} \right) \right]\) will be calculated and its isosurface map will be shown in a GUI window. It is important to note that you can input a range of pair indices to obtain their total density; for example, if you input 2,4-6,9, then the grid data of \(\sum_{i=2,4,5,6,9} {v}_{i}\left[ {\varphi}_{i}^{2}\left( \mathbf{r} \right)-{\varphi}_{-i}^{2}\left( \mathbf{r} \right) \right]\) will be calculated and plotted as isosurface. If there are e.g. 50 NOCVs in total and you input 1-50, then the grid data will simply correspond to the orbital deformation density Δρorb. In addition, it is worth to note that for open-shell case, pair index of alpha and beta spins is different (as shown on screen), therefore you can input proper indices to view sum of specific alpha and beta NOCV pair(s).
If you want to visualize the ΔρPauli, Δρorb, and Δρ mentioned above, you can respectively select “3 Show isosurface of Pauli deformation density”, “4 Show isosurface of orbital deformation density” and “5 Show isosurface of total deformation density”. The ΔρPauli allows you to vividly understand how electrons are repulsed in the region between fragments due to Pauli repulsion. Δρorb enables you to graphically examine electron concentration between fragments due to formation of covalent interaction or electron transfer between fragments due to mix of occupied orbitals of fragment(s) and unoccupied orbitals of other fragments. The Δρ is simply the sum of ΔρPauli and Δρorb, in fact it can also be calculated by making use of the custom operation feature of main function 5 (Section 3.7.1).
The grid data of a NOCV orbital, NOCV pair(s), ΔρPauli, Δρorb, and Δρ can also be exported as .cub file so that you can visualize them in third-part visualization softwares, you just need to select corresponding one of options 6~10.
When selecting options 2 to 8 to visualize or export various kinds of densities, if you have not defined grid setting by “-5 Set grid for calculation of various densities”, you will be automatically asked to set grid first. Next time you will not need to set the grid again, but you can still change grid setting anytime via option -5.
If you are interested in the occupied orbitals of \({\Psi}_{A}{\Psi}_{B}\), \({\Psi}_{AB}^{0}\) and \({\Psi}_{AB}\) mentioned in Section 3.26.1, you can respectively select “11 Visualize promolecular orbitals”, “12 Visualize frozen state orbitals” and “13 Visualize actual complex orbitals” to visualize the corresponding orbital isosurfaces. The orbitals of \({\Psi}_{A}{\Psi}_{B}\) are simply the orbitals in the fragment wavefunction files, the orbitals of \({\Psi}_{AB}\) are simply those in the complex wavefunction file. By comparing the orbitals in \({\Psi}_{AB}^{0}\) and those in \({\Psi}_{A}{\Psi}_{B}\), you can examine how Pauli repulsion between electrons of the fragments (realized by Löwdin orthogonalization as mentioned earlier) deform the occupied fragment molecular orbitals. If isovalue has been set to a small enough value, you will be able to observe the additional nodal plane(s) on the fragment orbitals of \({\Psi}_{AB}^{0}\) caused by the Löwdin orthogonalization. Note that the orthogonalization is not applied to unoccupied orbitals during ETS-NOCV analysis in Multiwfn, hence the unoccupied orbitals in \({\Psi}_{AB}^{0}\) and \({\Psi}_{A}{\Psi}_{B}\) are exactly the same.
Studying composition of NOCV pairs and that of corresponding NOCV orbitals is often useful when you want to better understand their nature. If you choose option “14 Calculate composition of NOCV orbitals and pairs”, you will be asked to choose a NOCV pair, then contributions from every basis function, shell, angular moment, and atom to the NOCV pair and corresponding NOCV orbitals will be printed on screen. The compositions of NOCV orbitals are evaluated using the SCPA method mentioned in Section 3.10.3. The composition of NOCV pair is calculated as \({v}_{i}{\Theta}_{i}+{v}_{j}{\Theta}_{j}\), where \(i\) and \(j\) are indices of the two paired NOCV orbitals, \({v}_{j}=-{v}_{i}\), and Θ is orbital composition evaluated by SCPA method. This option is useful in studying electron transfer. For example, if using this option you find contribution of atom A and \(B\) to NOCV pair 1 is 34.12% and -20.53%, respectively, that means atom A gained 0.3412 electron and atom B lost 0.2053 electron due to the interaction characterized by NOCV pair 1. Note that diffuse functions should never be used since SCPA method is incompatible with them. If you really have to use diffuse functions, after performing ETS-NOCV analysis, you should return to main menu (in this case the orbitals in memory are NOCV orbitals) and then use main function 8 to compute orbital composition via e.g. Hirshfeld or Becke method, which works normally when diffuse functions are presented.
Option “-6 Manually define correspondence between NOCV pairs and orbitals” is sometimes useful. Multiwfn automatically constructs NOCV pairs according to NOCV orbitals sorted by eigenvalues. For symmetric systems (e.g. linear system Ne...BeO), some NOCV orbitals may be degenerated, namely having exactly the same eigenvalues; in this case their order (index) is arbitrary, and the automatically determined correspondence between NOCV pairs and orbitals may be unexpected. In this case you can use option -6 to choose a NOCV pair, and then manually input indices of the two orbitals that the pair should correspond to. This option can be used multiple times to redefine multiple NOCV pairs.
Abundant examples of ETS-NOCV analysis are given in Section 4.23.