Visualizing Complex-Valued Molecular Orbitals

TL;DR

This paper introduces a visualization tool for complex-valued molecular orbitals that encodes phase information on mesh vertices, revealing rich chemical and physical insights often overlooked in traditional methods.

Contribution

The work presents a novel program for visualizing complex-valued molecular orbitals, highlighting phase information and its significance in chemical and physical analysis.

Findings

Visualizations reveal phase information in molecular orbitals.

Application to systems with spin-orbit coupling and magnetic fields.

Enhanced understanding of chemical and physical properties.

Abstract

We report an implementation of a program for visualizing complex-valued molecular orbitals. The orbital phase information is encoded on each of the vertices of triangle meshes using the standard color wheel. Using this program, we visualized the molecular orbitals for systems with spin-orbit couplings, external magnetic fields, and complex absorbing potentials. Our work has not only created visually attractive pictures, but also clearly demonstrated that the phases of the complex-valued molecular orbitals carry rich chemical and physical information of the system, which has often been unnoticed or overlooked.

Figures

Equations

$$|{\psi }_i(\mathbf{r})|=\left|\sum _{\mu }{C}_{\mu i}{\varphi }_{\mu }(\mathbf{r})\right|$$

$${\partial }_w|{\psi }_i(\mathbf{r})|=\Re [{\partial }_w{\psi }_i(\mathbf{r})]\cos \theta (\mathbf{r})+\text{imaginary}[{\partial }_w{\psi }_i(\mathbf{r})]\sin \theta (\mathbf{r}),$$

$${\partial }_w{\psi }_i(\mathbf{r})=\sum _{\mu }{C}_{\mu i}{\partial }_w{\varphi }_{\mu }(\mathbf{r}),$$

$$\theta (\mathbf{r})=\arg [{\psi }_i(\mathbf{r})],$$

$$\displaystyle\begin{array}[]{lllll}-3\leq p<-2&&R=1&G=p+3&B=0\\ -2\leq p<-1&&R=-1-p&G=1&B=0\\ -1\leq p<0&&R=0&G=1&B=p+1\\ 0\leq p<1&&R=0&G=1-p&B=1\\ 1\leq p<2&&R=p-1&G=0&B=1\\ 2\leq p<3&&R=1&G=0&B=3-p\\ \end{array}$$

$$\hat{\mathbf{p}}\to \hat{\mathbf{\pi }}=-i\nabla +\mathbf{A}$$

$${\varphi }_i(\mathbf{r})\to {\varphi }_i^{\prime }(\mathbf{r})=\exp \left[-i\mathbf{A}(\mathbf{R})\cdot \mathbf{r}\right]{\varphi }_i(\mathbf{r})$$

$$(H_{\mathrm{SO},1}^{\lambda }{)}_{pq}=-\frac{i}{4c^2}\sum _A{ϵ}_{\mu \nu }^{\lambda }\text{matrixelement}\frac{\partial {\varphi }_p}{\partial \mu }\frac{Z_A}{r_A}\frac{\partial {\varphi }_q}{\partial \nu }$$

$$\hat{H}=\sum _i\hat{h}(i)+\sum _{i<j}\hat{g}(i,j),$$

$$\hat{h}(i)=c^2(\beta -I_4)+c(\mathbf{\alpha }\cdot {\hat{\mathbf{p}}}_i)-\sum _A^{\mathrm{atoms}}\frac{Z_A}{r_{iA}},$$

$${\psi }_{L+}^{+}(\mathbf{r})=[({\psi }_{L-}^{-}(\mathbf{r}){]}^{\ast },$$

$${\psi }_{L-}^{+}(\mathbf{r})=-[{\psi }_{L+}^{-}(\mathbf{r}){]}^{\ast }.$$

Full text

I Introduction

The phase of molecular orbitals is fundamentally important in predicting and understanding chemical reactions. One of the earliest examples in the literature is the demonstration that the selectivity of the Diels–Alder reactions can be explained by the phases of the highest-occupied and lowest-unoccupied molecular orbitals.Woodward and Hoffmann (1965) Molecular orbitals have been ever since considered as a key descriptor for chemical reaction mechanisms. Many of the research articles in computational chemistry thus report the graphic pictures of molecular orbitals, which are often obtained by the standard visualization softwares such as Molden,Schaftenaar and Noordik (2000) Avogadro,Hanwell et al. (2012) IQMol,Gilbert IboViewKnizia (2013); Knizia and Klein (2015) to name a few.

To the best of our knowledge, visualization of molecular orbitals (and associated properties) has been thus far limited to that of real-valued orbitals due to the lack of implementation. This is partly because the software for real-valued orbitals does suffice for traditional electronic structure simulations in which the non-relativistic time-independent Schrödinger equation is solved for bound states with open boundary conditions, owing to the fact that the Hamiltonian for such systems is real and the phase of the orbitals can be chosen to be either $+1$ or $-1$. There are, however, various situations in which the molecular orbitals become complex, including simulation of systems with relativisticReiher and Wolf (2009)/magneticTellgren et al. (2008); Reynolds and Shiozaki (2015) contributions, that with periodicManby (2010)/absorbingRiss and Meyer (1993); Jagau et al. (2017) boundary conditions, and that with explicit time dependence.Goings et al. (2018) Another example is complex generalized Hartree–Fock wave functions that can describe non-coplaner spin polarization.Fukutome (1981); Henderson et al. (2018) Since computational tools for such systems are becoming widely available in recent years, it is of great importance to be able to visualize such orbitals.

In this work, we have developed a new computer program that allows for visualizing complex-valued molecular orbitals (and other properties). We apply this program to systems under an external magnetic field, systems with spin–orbit coupling, and systems with the absorbing boundary condition. Numerical examples will be shown to demonstrate surprisingly rich information that is pertained in the phase of those orbitals.

II Technical details

To compute the isosurface of a complex-valued molecular orbital, we first take the norm of the orbitals and reuse the components in the IboView programsKnizia (2013); Knizia and Klein (2015) that have been written for real-valued molecular orbitals. To do this, we simply evaluate the complex orbitals in real space and take the norm at each point:

[TABLE]

where $\psi_{i}$ is an $i$-th molecular orbital, $\phi_{\mu}$ are basis functions that are real valued, and $C_{\mu i}$ are the molecular-orbital coefficients that are complex. We use the so-called Marching Cubes algorithmLorensen and Cline (1987) for computing the triangular (polygon) meshes of the orbital isosurface. Following the default setting with IboView, the isosurfaces are constructed such that 80% of the molecular orbitals are encapsulated within the surface.

The normals of the isosurfaces are computed exactly from the derivative of the wave functions at each of the vertices. This has also been implemented in IboView and reused in this work.Knizia (2013); Knizia and Klein (2015) The derivative of the norm can be computed as

[TABLE]

where $\partial_{w}=\partial/\partial w$ and $w=x$, $y$, and $z$.

Subsequently, we calculate the phase angle $\theta$ of the molecular orbital at each of the vertices of the triangular meshes. The phase angle is then converted to a color code based on the RGB color wheel that maps the hue of the colors to [0, $2\pi$], as shown in Fig. 1, or, more quantitatively,

[TABLE]

with $p=3\theta/\pi$. Here the RGB scale takes the values between 0 and 1. The resulting color codes are added to the vertex information before the surfaces are rendered. Note that the color wheel is routinely used for visualizing complex functions in standard software programs, e.g., Mathematica.Wolfram Research, Inc.

The above algorithms are implemented using the locally-modified IboView program. The computation of Eqs. (1) and (2) are performed in the modified version of the BAGEL programs.bag ; Shiozaki (2018)

III Examples

III.1 Molecules in external magnetic fields

When molecules are placed under an external magnetic field, the wave function becomes complex, because the momentum operator in the Hamiltonian is modified as

[TABLE]

where $\mathbf{A}$ is a vector potential generated by the external magnetic field. Furthermore, one often uses the so-called gauge-including atomic orbitals (GIAO) to retain gauge-origin invariance,London (1937); Wolinski et al. (1990)

[TABLE]

where $\mathbf{R}$ is the location of the atom that the basis function $i$ belongs to. Note that the use of GIAO-based programs for studying systems under an external magnetic field is relatively new, pioneered by Tellgren and co-workers,Tellgren et al. (2008) who have reported interesting chemical applications in astrochemistry,Lange et al. (2012); Stopkowicz (2018) followed by one of the authors.Reynolds and Shiozaki (2015); Reynolds et al. (2018) In this section, we used the GIAO-based Hartree–Fock program in the BAGEL packagebag as described in Ref. Reynolds and Shiozaki, 2015.

In Fig. 2, we compare the $\pi$ orbitals of a benzene molecule with and without an external magnetic field. The benzene molecule is at its equilibrium geometry on the $xy$-plane, and a magnetic field of 5 T was applied in the $z$ direction. The def2-SVP basis set and corresponding fitting basis sets were used.Schäfer et al. (1992) In the absence of magnetic fields, there are two sets of doubly-degenerate $\pi$ orbitals (one bonding and one anti-bonding orbitals). The degenerate sets become split when the magnetic field is applied; this is because the magnetic field favors counter-clockwise electronic current in the direction of the field. Therefore, the degenerate occupied orbitals take linear combinations such that the phase changes $\pm 2n\pi$ when rotating around the principal axis (where $n$ is the number of nodes perpendicular to the $xy$-plane for the orbitals in the absence of magnetic fields). The hybridization phenomena is clearly depicted in the figures, facilitating intuitive understanding of the molecular orbitals in the presence of the external magnetic field. The amount of phase changes $\pm 2n\pi$ and $\pm 4n\pi$ cannot be read from the phaseless counterparts.

III.2 Spin–orbit coupled wave functions

One of the approaches to taking into account spin–orbit couplings is to use the so-called state interaction method.Malmqvist et al. (2002); Sayfutyarova and Chan (2016) In this approach, one first obtains a number of pure spin states in the absence of spin–orbit couplings, which are then used as a basis to construct an effective Hamiltonian that includes spin–orbit couplings. The spin–orbit coupling elements are complex: for instance, the one-body part of the spin–orbit coupling matrix elements from the Breit–Pauli Hamiltonian areReiher and Wolf (2009)

[TABLE]

where the left-hand side are the matrix elements associated with the Pauli operator $\hat{\sigma}^{\lambda}_{pq}$, $\lambda$, $\mu$, and $\nu$ label Cartesian components, $c$ is the speed of light, $\epsilon_{\mu\nu}^{\lambda}$ is the Levi–Civita symbol, $p$ and $q$ label basis functions, and $A$ labels nuclei.

As a simple example, we considered PbH2 that exhibits spin–orbit coupling of about $2500$ cm*-1*,Matsunaga et al. (1996) which is a heavy-element analogue of the carbene biradical. The geometry was set to $r=1.880$ Å and $\theta=91.5^{\circ}$ with the $C_{2v}$ symmetry. One singlet state and one triplet state were included in the effective Hamiltonian. We first performed a state-averaged CASSCF calculation with the spin-free DKH2 Hamiltonian,Hess (1986) averaging over both the singlet and triplet states. The full-valence active space (6 electrons in 6 orbitals) and the ANO-RCC basis setRoos et al. (2004) were used. Subsequently, the spin-nonconserving density matrices associated with the Pauli operator $\hat{\sigma}^{\lambda}_{pq}$ were computed. The effective Hamiltonian was then formed and diagonalized to obtain the spin–orbit coupled states (total of 4 spin–orbit states). In Fig. 3, we present the natural orbitals for one of the triplet states that is coupled to the singlet state (i.e., the state that is predominantly $|1,1\rangle+|1,-1\rangle$ using the $|L,L_{z}\rangle$ notation). From this figure, one can visually determine the orbitals that play an important role in spin–orbit interaction and their phase relationshiops. The third and fourth orbitals in each column are the hybridized singly occupied orbitals that mix $\alpha$ and $\beta$ components owing to the spin–orbit coupling.

Another approach to accouting for spin–orbit interaction is the four-component relativistic approaches that seek to directly solve the Dirac equations that include spin–orbit interaction a priori.Reiher and Wolf (2009) The Dirac equation reads

[TABLE]

in which $\bm{\alpha}$ and $\beta$ are Dirac’s 4$\times$4 matrices; $I_{4}$ is a 4$\times$4 identity matrix; $c$ denotes the speed of light; $\hat{\mathbf{p}}=-i\nabla$ is the momentum operator; and, $\hat{g}(i,j)$ is the electron–electron interaction (in this work, we used the standard Coulomb operator) with $i$ and $j$ labeling electrons. Such four-component approaches have become applicable to large molecules of chemical interest in the past decade.Kelley and Shiozaki (2013); Storchi et al. (2013); Repiský (2014) One of the consequences of using the four-component wave function approaches is that the orbitals are to be represented by complex-valued 4-spinors; therefore, our visualization code for complex-valued molecular orbitals are particularly useful for the analyses of four-component wave functions.

As a simple example, the molecular spinors are presented for the triplet O2 molecule computed by the Dirac CASSCF method with density fitting as described in Ref. Reynolds et al., 2018 (note that other implementations of Dirac CASSCF includes those reported in Ref. Aa. Jensen et al., 1996; Thyssen et al., 2008). The def2-SVP basis set and corresponding fitting basis sets were used.Schäfer et al. (1992) CAS(2$e$,2$o$) was used and three states were averaged. Figure 4 shows the valence canonical orbitals of this system, in which only the large components (L$+$ and L$-$) of the 4-spinors are included. Both Kramers $+$/$-$ orbitals are shown. Note that the corresponding Kramers $+$/$-$ orbitals are degenerate due to the time-reversal symmetry; therefore, rotations between them are arbitrary. It is shown that both bonding and anti-bonding $\pi$ orbitals are hybridized in the presence of spin–orbit coupling that breaks the degeneracy. In addition, one can clearly observe the time-reversal relation between the Kramers $+$/$-$ pairs,

[TABLE]

Visually confirming this time-reversal symmetry is only possible using visualization code for complex-valued orbitals with phases. These relationships are especially non-trivial for hybridized orbitals; when checking these relationships, notice in Fig. 1 that the phase angles 0 and $\pi$ are cyan and red, respectively.

III.3 Molecules with metastable states

The renewed interest in metastable states in the past decade has inspired the development of various advanced computational tools.Riss and Meyer (1993); Jagau et al. (2017) When modeling metastable states, in which the bound-state wave functions are mixed with continuum wave functions, one often uses an effective Hamiltonian that truncates the continuum degrees of freedom in one way or another; and this truncation leads to Hamiltonians whose eigenvalues are complex. The imaginary part of the eigenvalues is related to the lifetime of the state. Since the associated wave functions are complex, analysis of excited-state wave functions warrants visualization of complex-valued orbitals. Usually, the real and imaginary parts of the orbitals are presented separately.Jagau et al. (2017)

Here we show an example that employs the so-called complex absorbing potential (CAP) using the second-order multireference perturbation theory, which has been reported by Bravaya and co-workersKunitsa et al. (2017) using an uncontracted multireference perturbation theory called the XMCQDPT2 methodGranovsky (2011) (in this work, we use instead the fully contracted variant, XMS-CASPT2, as implemented in BAGELMacLeod and Shiozaki (2015); Vlaisavljevich and Shiozaki (2016)). In short, this approach includes both dynamical electron correlation and complex absorbing potentials as perturbations in the post-CASSCF process. For simplicity, the natural orbitals are computed from the CASSCF part of the wave functions.

We computed the XMS-CASPT2 effective Hamiltonian $H_{\mathrm{eff}}$ and the CAP matrix $W_{\mathrm{cap}}$ within the model space for N2 anion ($r=2.074$ bohr). N2 anion has been theoretically studied previously (Ref. Berman et al., 1983; see also Ref. Zuev et al., 2014 and references therein). The aug-cc-pVQZ basis set was used together with three even-tempered diffuse functions in both $s$ and $p$ shells. The active space included 16 orbitals and 5 electrons, and 14 states were averaged in the calculation without imposing spacial symmetry. The spherical CAP with 4.88 bohr radius was used. The so-called $\eta$ trajectory was obtained by diagonalizing $H_{\mathrm{eff}}-\eta W_{\mathrm{cap}}$ with various $\eta$; the result is compiled in Fig. 5. The stationary point was obtained by minimizing $\eta|dE/d\eta|$ to be $\eta=0.000206$, $E=-109.30623-0.00971i$ (shown as a green dot in Fig. 5). Note that the electronic resonance width $\Gamma$ computed from this result (0.528 eV) is in good agreement with the experimentally derived value (0.41 eV).Berman et al. (1983)

The natural orbitals in the active space with non-zero occupation numbers are shown in Fig. 6. We here use the convention in which the molecule is aligned with the $z$ axis, and the $x$ axis is in-plane. There are two (almost) doubly occupied orbitals ($\pi_{x}$ and $\pi_{y}$) and a singly occupied orbital ($\pi_{x}^{\ast}$). The $\pi_{y}$ and $\pi_{y}^{\ast}$ orbitals resemble very closely the neutral N2 orbitals, and their occupation numbers sum up to 2. It is also noted that the imaginary part of the occupation numbers for $\pi_{y}$ and $\pi_{y}^{\ast}$ is nearly zero. In contrast, the singly occupied orbital $\pi_{x}^{\ast}$ is partially hybridized with the diffuse orbitals that describe continuum states. Furthermore, there is another hybridized orbital $\pi_{x}^{\ast\prime}$ that is a mixture of $\pi_{x}^{\ast}$ and the continuum. The occupation number for this $\pi_{x}^{\ast\prime}$ orbital is $(-0.005,0.006)$, which suggests that it plays a role in the auto-detachment process. Quantitative interpretation of these occupation numbers requires further research.

IV Conclusions

In this work, we developed a software program that allows for visualizing complex-valued molecular orbitals, in which the phases of the complex-valued orbitals are represented by the RGB color wheel. The implementation has been on the basis of the IboView programKnizia (2013); Knizia and Klein (2015) and the BAGEL package.bag Examples were presented for systems under a magnetic field, those with spin–orbit interactions, and those with metastable states. These examples have shown that the phases of the complex-valued orbitals do carry chemical and physical information (hybridization, time-reversal symmetry, and decay channels in these examples, respectively).

Acknowledgements.

Professor Ksenia Bravaya is thanked for help with the calculations for metastable states. We appreciate useful discussions with Professor Sandeep Sharma. R.A.-S. and T.S. have been supported by Air Force Office of Scientific Research (AFOSR FA9550-18-1-0252) and by the National Science Foundation CAREER Award (CHE-1351598), respectively. T.S. is a Sloan Research Fellow.