Non-adbatic dynamics using trajectories

class: center, middle

## **Non-adiabatic Dynamics Using Trajectories**: Ehrenfest Dynamics and the Exact Factorization Method

##### *Shiling Liang* | *Section of Physics, EPFL*

- Born-Huang expansion and coupled-channel equation
- Born-Oppenheimer limit
- Ehrenfest dynamics
- Exact fractorization

---
## Background knowledge

### Classical: Netwon's equation

$$
M\_\gamma \frac{d^2\mathbf{R}\_{\gamma}}{dt^2}=-\nabla\_\gamma E(\mathbf{R})
$$

where $\mathbf{R}_{\gamma}$ is the coordinates of nucleus $\gamma$, $E(\mathbf{R})$ is the empirical potential energy surface.
### Quantum: Schrödinger equation for molecules
The Schrödinger equation for nucleus with electrons reads
$$
i\hbar\partial\_t\Psi(\mathbf{r},\mathbf{R},t)=\hat{H}\Psi(\mathbf{r},\mathbf{R},t)
$$

where $\Psi(\mathbf{r},\mathbf{R},t)$ is the wave function of the electron-nuclies, $\mathbf{r}$ is the electron positions, $\mathbf{R}$ is the nuclei positions. $\hat{H}$ is the Hamiltonian of the system, which consists the following parts

$$
\hat{H}=\hat{T}\_{nu}+ \underbrace{\hat{T}\_{el} + \hat{V}\_{ee}+\hat{V}\_{en}+\hat{V}\_{nn}}\_{\hat{H}\_{el}(\mathbf{R},\mathbf{r})}
$$

---
## The derivation of coupled-channel equation using Born-Huang expansion
The Born-Huang expand the total wave function as
$$
\Psi(\mathbf{r},\mathbf{R},t)=\sum\_i^{\infty}\Omega\_i(\mathbf{R},t)\Phi\_i(\mathbf{r};\mathbf{R}).
$$
where $\Phi\_i(\mathbf{r};\mathbf{R})$ is the eigenwavefunction of electronic Hamiltonian $\hat{H}\_{e}(\mathbf{r};\mathbf{R})$ with eigenenergy $E^{el}\_i$.

Inserting this expansion back to the TDSE gives the time evolution equation of the coefficient $\Omega\_i(\mathbf{R},t)$

$$
i\hbar \partial\_t\Omega\_j(\mathbf{R},t)=\left[-\hat{T}\_{nu}+E\_j^{el}\right]\Omega\_j(\mathbf{R},t)+\sum\_i^{\infty}\mathcal{F}\_{ij}(\mathbf{R})\Omega\_i (\mathbf{R},t)
$$

where $\mathbf{F}\_{ji}(\mathbf{R})$ is the coupling matrix between channel $i$ and $j$

$$
\mathcal{F}\_{ji}(\mathbf{R})
=\int d\mathbf{r}\Phi\_j^\* (\mathbf{r},\mathbf{R})\hat{T}\_{nu}\Phi\_i(\mathbf{r};\mathbf{R})+\sum\_{\gamma}\frac{1}{M\_\gamma}\left[\int d\mathbf{r}\Phi\_j^\*   (\mathbf{r};\mathbf{R})\hat{p}\_{\gamma}\Phi\_i(\mathbf{r};\mathbf{R}) \right]\hat{p}\_\gamma
$$

---
## Consider the Born-Oppenheimer limit of coupled-channel equation

If we neglect the coupling terms, $\mathcal{F}\_{ji}(\mathbf{R})$, the Born-Huang expansion is reduced to the Born-Oppenheimer approximation, i.e. adiabatic approximation. In the adiabatic approximation, the total wavefunction can be simply written as the product of electronic wavefunction and the nuclear wave function

$$
\Psi(\mathbf{r},\mathbf{R},t)=\Omega\_j(\mathbf{R},t)\Phi\_j(\mathbf{r};\mathbf{R})
$$

To get the equation of motion of the nucleus, we need to find the ground state of electrons to give the potential energy surface (PES), $E_j^{el}(\mathbf{R})$.

$$
M\_\gamma\ddot{\mathbf{R}}^\gamma\_j(t)=-\nabla\_\gamma E\_j^{el}(\mathbf{R}(t))
$$
---
## Ehrenfest dynamics (mean-field)

### Key assumptions
The molecular wavefunction in Ehrenfest dynamics reads

$$
\Psi(\mathbf{r},\mathbf{R},t)=\Phi(\mathbf{r},t)\Omega(\mathbf{R},t)\exp\left[\frac{i}{\hbar}\int\_{t\_0}^{t}dt'E\_{el}(t')\right]
$$

where

$$
E\_{el}(t)=\iint d\mathbf{r} d\mathbf{R} \Phi^\*   (\mathbf{r},t)\Omega^\*  (\mathbf{R},t)\hat{H}\_{el}(\mathbf{r},\mathbf{R}) \Phi (\mathbf{r},t)\Omega(\mathbf{R},t)
$$

subsituting the anstaz wavefunction backto the time-dependent Schrödinger equation one can finally obtains

One need to use these two relations

$$
\begin{aligned}
\int d\mathbf{r}\Phi^\*  (\mathbf{r},t)\partial\_t\Phi(\mathbf{r},t)=E\_{el}(t)\\\
\int d\mathbf{R}\Omega^\*  (\mathbf{R},t)\partial\_t\Omega(\mathbf{R},t)= E
\end{aligned}
$$
---
## Ehrenfest dynamics (mean-field)
### Results
the equation of mothion for nuclear and electron wavefunctions reads

$$
\begin{aligned}
i\hbar\partial\_t\Phi &= \hat{T}\_{el}\Phi+\underbrace{\left[\int d\mathbf{R}\Omega^\*  (\hat{V}\_{ee}+\hat{V}\_{en}+\hat{V}\_{nn})\Omega\right]}\_{\langle \hat{V}\rangle\_{nu}}\Phi\\\
i\hbar\partial\_t\Omega &= \hat{T}\_{nu}\Omega + \underbrace{\left[\int d\mathbf{R}\Phi^\*  \hat{H}\_{el}\Phi\right]}\_{\langle \hat{H}\_{el}\rangle\_{el}}\Omega\\\
\end{aligned}
$$

and then taking **classical limit**, one can finnal obtain

$$
\begin{aligned}
i\hbar\partial\_t\Phi &= \hat{H}\_{el}\Phi\\\
M\_\gamma \ddot{R}\_{\gamma}&=-\nabla\_{\gamma}\left\langle \hat{H}\_{el}(\mathbf{r},\mathbf{R})\right\rangle\_t
\end{aligned}
$$
---
## Exact factorization method
Total wavefunction = electronic wavefucntion $\times$ nuclear wavefunction
$$
\Phi(\mathbf{r},\mathbf{R},t)=\Phi(\mathbf{r};\mathbf{R},t)\Omega(\mathbf{R},t)
$$

The advantage is that both two terms have physical meaning: the **marginal nuclear probability distribution** and the **conditional electronic distribution**

$$
\begin{aligned}
 \mathrm{P}\_{el}(\mathbf{r}|\mathbf{R},t) = |\Phi(\mathbf{r};\mathbf{R},t)|^2, \quad
\mathrm{P}\_{nu}(\mathbf{R}) = |\Omega(\mathbf{R},t)|^2 
\end{aligned}
$$

The key assumption to make these interpretation ture is the normalization condition of $\mathrm{P}\_{el}$
$$
\int d\mathbf{r} \mathrm{P}\_{el}(\mathbf{r}|\mathbf{R},t)=1
$$

---
## Comparing the Ehrenfest dynamics and the exact factorization method

</center>
- Exact factorization gives a rigrous defination of the time-dependent potential energy surface, while the Ehrenfest dyanmics gives mean-field potential energy surface.

- More accurate.

- Exact factorization is a rigorous starting point for making approximations.

- Exact factorization can capture tunneling (tunneling is the leading mechanism for the dissociation)
---
## (pseudo)Code implementation of Ehrenfest dynamics

```python
initialize phi, R, v_R, t, dt 
# calculate the potential energy surface using mean-field appraoch
E_el(t) = <phi(t)|H_el|phi(t)>
for k steps
  # update the position of nucleus
  R(t+dt) = R(t-dt/2) + v(t)*dt  # Euler integral
  # update the electronic wave function using
  phi(t+dt) = U(E_el*dt)*phi(r,R,t) #M-SOFT
  # update potential energy
  E_el(t+dt) = <phi(t+dt)|H_el|phi(t+dt)|phi>
  # update velocity
  v_R(t+dt) = v_R(t) + (E_el(t+dt)-E_el(t))/((R(t+dt)-R(t))*M)
  #update time
  t = t + dt 
end for
```

## References
[1] De Carvalho, F.F., Bouduban, M.E., Curchod, B.F. and Tavernelli, I., 2014. Nonadiabatic molecular dynamics based on trajectories. Entropy, 16(1), pp.62-85.

[2] Tully, J.C., 1998. Nonadiabatic dynamics. In Modern methods for multidimensional dynamics computations in chemistry (pp. 34-72).

[3] Li, X., Tully, J.C., Schlegel, H.B. and Frisch, M.J., 2005. Ab initio Ehrenfest dynamics. The Journal of chemical physics, 123(8), p.084106.

[4] Abedi, A., Maitra, N.T. and Gross, E.K., 2010. Exact factorization of the time-dependent electron-nuclear wave function. Physical review letters, 105(12), p.123002.

---
# Thanks!