NAC Calculation for 1D Model Hamiltonian
Post
Cancel

# NAC Calculation for 1D Model Hamiltonian

## Introduction

Let’s start with a parametric 1D Hamiltonian $H(r)$

$$$H(r) = \begin{pmatrix} \varepsilon_1^d(r) & \lambda(r) \\[6pt] \lambda(r) & \varepsilon_2^d(r) \end{pmatrix}$$$

where the dependence of the the matrix elements on the parameter $r$ are assumed to be

\begin{align} \varepsilon_1^d(r) &= \phantom{-}\alpha\cdot r \\[6pt] \varepsilon_2^d(r) &= -\alpha\cdot r \\[6pt] \lambda(r) &= \gamma \cdot e^{-r^2 / \beta^2} \end{align}

Obviously, the derivative of the Hamiltonian with respect to $r$ can be written as

$$$\nabla_r H(r) = \begin{pmatrix} \alpha & -\frac{2\gamma}{\beta^2}\cdot r \cdot e^{-r^2/\beta^2} \\[6pt] -\frac{2\gamma}{\beta^2}\cdot r\cdot e^{-r^2/\beta^2} & -\alpha \end{pmatrix}$$$

The eigenvalues of the Hamiltonian can be solved analytically,

\begin{align} \varepsilon_i^a &= \frac{(\varepsilon_1^d + \varepsilon_2^d)}{2} \pm \frac{ \sqrt{(\varepsilon_1^d - \varepsilon_2^d)^2 + 4\lambda^2} }{2} \\[6pt] &= \pm \sqrt{\alpha^2\cdot r^2 + \gamma^2\cdot e^{-2r^2/\beta^2}} \qquad (i=1,2) \label{eq:eig_val} \end{align}

Apparently, the energy gap at $r=0$ is $\Delta \varepsilon = 2\lambda$.

The corresponding eigenvectors can also be obtained

\begin{align} \psi_i^a &= \begin{pmatrix} \dfrac{-2\lambda}{ (\varepsilon_1^d - \varepsilon_2^d) \mp \sqrt{(\varepsilon_1^d - \varepsilon_2^d)^2 + 4\lambda^2} } \\[12pt] 1 \end{pmatrix} \\[12pt] &= \begin{pmatrix} \dfrac{-\gamma\cdot e^{-r^2/\beta^2}}{ \alpha r \mp \sqrt{\alpha^2\cdot r^2 + \gamma^2\cdot e^{-2r^2/\beta^2}} } \\[12pt] 1 \end{pmatrix} \label{eq:eig_vec} \end{align}

Note that the above eigenvectors is NOT normalized and should therefore be properly normalized when implemented in codes.

 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 #!/usr/bin/env python import sympy as sp from sympy.matrices import Matrix E1, E2, L = sp.symbols("E1 E2 L") Hamil = Matrix([ [E1, L], [L, E2], ]) sp.pprint(Hamil.eigenvects()) # U, T = Hamil.diagonalize(normalize=True) # sp.pprint(U)

The nonadiabatic coupling vector (NAC), which is a scalar in 1D, can now be written as

\begin{align} d_{ij} &= \frac{ \langle \psi_i |\nabla_r H(r)|\psi_j \rangle }{ \varepsilon_j^a -\varepsilon_i^a } \\[9pt] & = \dfrac{ \begin{pmatrix} c_{i1}^* & c_{i2}^* \end{pmatrix} \cdot \begin{pmatrix} \alpha & -\frac{2\gamma}{\beta^2}\cdot r \cdot e^{-r^2/\beta^2} \\[6pt] -\frac{2\gamma}{\beta^2}\cdot r\cdot e^{-r^2/\beta^2} & -\alpha \end{pmatrix} \cdot \begin{pmatrix} c_{j1} \\[6pt] c_{j2} \end{pmatrix} }{ (\varepsilon_j^a - \varepsilon_i^a) } \end{align}

In practice, the derivative of the Hamiltonian with respect to $r$ is hard to get and we often evaluate another quantity $D_{ij} = d_{ij}\cdot\dot{r}$, which is also referred to as nonadiabatic coupling.

\begin{align} D_{ij} & = \frac{ \langle \psi_i |\nabla_r H(r)|\psi_j \rangle }{ \varepsilon_j^a -\varepsilon_i^a }\cdot\dot r = \langle \psi_i | \frac{\partial}{\partial t} | \psi_j \rangle \\[9pt] &\approx \dfrac{ \langle \psi_i(t) |\psi_j(t+\Delta t)\rangle - \langle \psi_i(t+\Delta t) |\psi_j(t)\rangle }{ 2\Delta t } \end{align}

where the finite difference approximation in the second line was proposed by Hammes-Schiffer and Tully (HST).1

Below, we show an example results of the model Hamiltonian:

 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 #!/usr/bin/env python import numpy as np rmax = 0.75 nr = 501 rr = np.linspace(-rmax, rmax, nr, endpoint=True) alpha = 2.0; beta = 0.20; gamma = 0.20 E1 = alpha * rr E2 = -alpha * rr W0 = gamma * np.exp(-rr**2 / beta**2) # The Hamiltonian Ham = np.zeros((nr, 2, 2)) Ham[:,0,0] = E1 Ham[:,1,1] = E2 Ham[:,0,1] = W0 Ham[:,1,0] = W0 # The derivative of Hamiltonian to R Hderv = np.zeros((nr, 2, 2)) Hderv[:,0,0] = alpha Hderv[:,1,1] = -alpha Hderv[:,0,1] = -2 * rr * W0 / beta**2 Hderv[:,1,0] = -2 * rr * W0 / beta**2 E3 = -np.sqrt(alpha**2 * rr**2 + W0**2) E4 = -E3 V3 = np.array([ -W0 / (alpha*rr + np.sqrt(alpha**2 * rr**2 + W0**2)), np.ones(nr) ]).T # normalization V3 /= np.linalg.norm(V3, axis=1)[:,None] V4 = np.array([ -W0 / (alpha*rr - np.sqrt(alpha**2 * rr**2 + W0**2)), np.ones(nr) ]).T # normalization V4 /= np.linalg.norm(V4, axis=1)[:,None] nac = np.sum( V3 * np.matmul(Hderv, V4[...,None])[:,:,0], axis=1 ) / (E4 - E3) # |c_11|^2 rho_1 = V3[:,0]**2 # |c_21|^2 rho_2 = V4[:,0]**2
• Let us now focus on the NAC at $r = 0$, where

$$$\psi_1 = \begin{pmatrix} \frac{1}{\sqrt{2}} \\[6pt] \frac{1}{\sqrt{2}} \end{pmatrix}; \qquad \psi_2 = \begin{pmatrix} \frac{-1}{\sqrt{2}} \\[6pt] \frac{1}{\sqrt{2}} \end{pmatrix}; \qquad \nabla_r H(r) = \begin{pmatrix} \alpha & 0 \\[6pt] 0 & -\alpha \end{pmatrix}$$$

One can show that at $r = 0$,

$$$d_{12}(r=0) = \frac{\alpha}{(\varepsilon_1^a - \varepsilon_2^a)} = \frac{\alpha}{2\lambda}$$$

i.e. the NAC is inversely propotional to the energy gap at $r = 0$, which is consistent with the adiabatic theorem.2

One should not confuse the external perturbation $\lambda$ with NAC.

• Comparison of $D_{ij}$ from two methods. Suppose that the parameter $r$ follows the motion of a harmonic oscillator, i.e.

\begin{align} r &= -r_0 \cdot\cos(\omega t) \\[6pt] \dot r &= \omega\,r_0\cdot\sin(\omega t) \end{align}

## Band Ordering and Phase Correction

When the eigenvalues and eigenvectors are obtained numerically, there can be problems with the ordering of the band (eigenvalues may not be properly ordered) and the phase of the wavefunction (arbitrary phase can be assigned to the eigenvectors), both of which may affect the calculation of the NAC.

• No band reordering and NO phase correction.

 1 2 3 4 vals, vecs = np.linalg.eig(Ham) E3, E4 = vals.T V3 = vecs[...,0] V4 = vecs[...,1]

Eigenvalues returned by numpy.linalg.eig are not necessarily ordered, use numpy.linalg.eigh instead, which arrange the eigenvalues in ascending order. Here, we use numpy.linalg.eig only to show the effects of band ordering.

What if we introduce more artificial band crossings?

 1 2 3 4 5 6 7 8 9 10 11 12 # use np.linalg.eigh to get the correct band order and then reshuffle it vals, vecs = np.linalg.eigh(Ham) # randomize the band order oo = np.random.randint(2, size=(nr)) xx = oo - 1 E3 = vals[range(nr),oo] E4 = vals[range(nr),xx] V3 = vecs[range(nr),:,oo] V4 = vecs[range(nr),:,xx]
• With more artificial band crossings, as shown in the upper left panel, the wavefunction coefficients $|c_{i1}|^2$ in the upper right panel shows many abrupt changes. However, the overall envelope shape is correct.

• The NACs in the lower panel also become unsmooth with many abrupt changes. Moreover, there are envelopes enclosing the NAC curve. The lower and upper envelope are the correct NAC and that differ by a sign factor, respectively.

• With band reordering and NO phase correction.

 1 2 3 4 5 6 7 8 9 10 11 12 13 14 vals, vecs = np.linalg.eig(Ham) E3, E4 = vals.T V3 = vecs[...,0] V4 = vecs[...,1] oo = np.asarray(E3 <= E4, dtype=int) xx = oo - 1 # band reorder --- energy E3 = vals[range(nr),xx] E4 = vals[range(nr),oo] # band reorder --- wavefunction V3 = vecs[range(nr),:,xx] V4 = vecs[range(nr),:,oo]
• With band reordering and With phase correction.

 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 vals, vecs = np.linalg.eig(Ham) E3, E4 = vals.T V3 = vecs[...,0] V4 = vecs[...,1] oo = np.asarray(E3 <= E4, dtype=int) xx = oo - 1 # band reorder --- energy E3 = vals[range(nr),xx] E4 = vals[range(nr),oo] # band reorder --- wavefunction V3 = vecs[range(nr),:,xx] V4 = vecs[range(nr),:,oo] # phase correction V3[V3[:,1] < 0] *= -1 V4[V4[:,1] < 0] *= -1