Path integral molecular dynamics with surface hopping for thermal equilibrium sampling of nonadiabatic systems

TL;DR

This paper introduces a new ring polymer approach combined with surface hopping dynamics to efficiently sample the thermal equilibrium of nonadiabatic quantum systems, preserving electronic state discreteness.

Contribution

It proposes a novel ring polymer representation that maintains electronic state discreteness and develops a path integral molecular dynamics with surface hopping method for nonadiabatic systems.

Findings

The method accurately samples equilibrium distributions.

The approach is validated through theoretical analysis.

Numerical examples demonstrate effectiveness.

Abstract

In this work, a novel ring polymer representation for multi-level quantum system is proposed for thermal average calculations. The proposed presentation keeps the discreteness of the electronic states: besides position and momentum, each bead in the ring polymer is also characterized by a surface index indicating the electronic energy surface. A path integral molecular dynamics with surface hopping (PIMD-SH) dynamics is also developed to sample the equilibrium distribution of ring polymer configurational space. The PIMD-SH sampling method is validated theoretically and by numerical examples.

Tables

Table 1. Table 1: Errors in numerical empirical averages with 95 % confidence intervals and mean squared errors. The reference value is 0.640172 0.640172 0.640172 .

	${\beta }_N=\frac{1}{4}$	${\beta }_N=\frac{1}{8}$	${\beta }_N=\frac{1}{16}$
Error	5.80e-04	2.89e-04	1.52e-04
95% C.I.	2.23e-04	2.23e-04	2.24e-04
M.S.E.	3.49e-07	9.67e-08	3.61e-08

Table 2. Table 2: Error in numerical empirical averages with 95 % confidence intervals and mean squared errors. The reference value is 0.640172 0.640172 0.640172 .

	$\Delta t=\frac{1}{25}$	$\Delta t=\frac{1}{50}$	$\Delta t=\frac{1}{100}$	$\Delta t=\frac{1}{200}$	$\Delta t=\frac{1}{400}$
Error	6.62e-3	1.98e-3	1.34e-3	5.94e-3	2.08e-3
95% C.I.	7.22e-3	9.54e-3	9.48e-3	5.52e-3	5.07e-3
M.S.E.	5.74e-5	2.76e-5	2.55e-5	4.26e-5	7.77e-6

Table 3. Table 3: Errors in numerical empirical averages with 95 % confidence intervals and mean squared errors. The reference value is − 0.593497 0.593497 -0.593497 .

	$\Delta t=\frac{1}{25}$	$\Delta t=\frac{1}{50}$	$\Delta t=\frac{1}{100}$	$\Delta t=\frac{1}{200}$	$\Delta t=\frac{1}{400}$
Error	1.99e-1	1.87e-1	6.54e-2	1.50e-2	1.85e-3
95% C.I.	2.70e-2	2.09e-2	2.46e-2	2.05e-2	1.14e-2
M.S.E.	3.99e-2	3.51e-2	4.43e-3	3.34e-4	3.74e-5

Table 4. Table 4: Errors in numerical empirical averages with 95 % confidence intervals and mean squared errors. The reference value is − 0.593497 0.593497 -0.593497 .

	$\eta =5$	$\eta =10$	$\eta =20$	$\eta =40$
Error	8.27e-3	1.38e-3	4.97e-3	1.44e-3
95% C.I.	1.09e-2	7.67e-3	3.59e-3	1.87e-3
M.S.E.	9.93e-5	1.32e-5	2.80e-5	2.98e-6

Equations

$$\hat{H}=\hat{T}+\hat{V}=\frac{1}{2M}\left(\begin{array}{cc}{\hat{p}}^2& \\ & {\hat{p}}^2\end{array}\right)+\left(\begin{array}{cc}{V}_{00}(\hat{q})& V_{01}(\hat{q})\\ V_{10}(\hat{q})& V_{11}(\hat{q})\end{array}\right),$$

$$V(q)=\left(\begin{array}{cc}{V}_{00}(q)& V_{01}(q)\\ V_{10}(q)& V_{11}(q)\end{array}\right)$$

$$⟨\hat{A}⟩=\frac{{\mathrm{Tr}}_{ne}[e^{-\beta \hat{H}}\hat{A}]}{{\mathrm{Tr}}_{ne}[e^{-\beta \hat{H}}]},$$

$${\mathrm{Tr}}_{ne}={\mathrm{Tr}}_n{\mathrm{Tr}}_e={\mathrm{Tr}}_{L^2({\mathbb{R}}^d)}{\mathrm{Tr}}_{{\mathbb{C}}^2}.$$

$$\hat{A}=\left(\begin{array}{cc}{A}_{00}(\hat{q})& A_{01}(\hat{q})\\ A_{10}(\hat{q})& A_{11}(\hat{q})\end{array}\right).$$

$$\begin{gathered} {\mathrm{Tr}}_{ne}[e^{-\beta \hat{H}}]\approx {\mathcal{Z}}_N:=\frac{1}{(2\pi {)}^{dN}}{\int }_{{\mathbb{R}}^{2dN}}\mathrm{d}\mathbf{q}\mathrm{d}\mathbf{p}\sum _{\mathbf{\ell }\in \{0,1{\}}^N} \\ \times \exp (-{\beta }_N{H}_N(\mathbf{q},\mathbf{p},\mathbf{\ell })). \end{gathered}$$

$$H_N(\mathbf{q},\mathbf{p},\mathbf{\ell })=\sum _{k=1}^N⟨{\ell }_k|G_k|{\ell }_{k+1}⟩,$$

$$\langle\ell\rvert G_{k}\lvert\ell^{\prime}\rangle=\begin{cases}\displaystyle\frac{p^{2}_{k}}{2M}+\frac{M\left(q_{k}-q_{k+1}\right)^{2}}{2(\beta_{N})^{2}}+V_{\ell\ell}(q_{k})-\frac{1}{\beta_{N}}\ln\Bigl{(}\cosh\bigl{(}\beta_{N}\left\lvert V_{01}(q_{k})\right\rvert\bigr{)}\Bigr{)},&\ell=\ell^{\prime},\\ \displaystyle\frac{p^{2}_{k}}{2M}+\frac{M\left(q_{k}-q_{k+1}\right)^{2}}{2(\beta_{N})^{2}}+\frac{V_{00}(q_{k})+V_{11}(q_{k})}{2}-\frac{1}{\beta_{N}}\ln\Bigl{(}\sinh\bigl{(}\beta_{N}\left\lvert V_{01}(q_{k})\right\rvert\bigr{)}\Bigr{)},&\ell\neq\ell^{\prime}.\end{cases}$$

$$\begin{gathered} {\mathrm{Tr}}_{ne}[e^{-\beta \hat{H}}\hat{A}]\approx \frac{1}{(2\pi {)}^{dN}}{\int }_{{\mathbb{R}}^{2dN}}\mathrm{d}\mathbf{q}\mathrm{d}\mathbf{p}\sum _{\mathbf{l}\in \{0,1{\}}^N} \\ \times \exp (-{\beta }_N{H}_N)W_N[A], \end{gathered}$$

$$W_{N}[A](\boldsymbol{q},\boldsymbol{p},\boldsymbol{\ell})=\frac{1}{N}\sum_{k=1}^{N}\biggl{(}\langle\ell_{k}\rvert A(q_{k})\lvert\ell_{k}\rangle-\\ -e^{\beta_{N}\langle\ell_{k}|G_{k}|\ell_{k+1}\rangle-\beta_{N}\langle\bar{\ell}_{k}\rvert G_{k}\lvert\ell_{k+1}\rangle}\langle\ell_{k}\rvert A(q_{k})\lvert\bar{\ell}_{k}\rangle\frac{V_{\ell_{k}\bar{\ell}_{k}}}{\lvert V_{\ell_{k}\bar{\ell}_{k}}\rvert}\biggr{)},$$

$$⟨\hat{A}⟩\approx \frac{1}{(2\pi {)}^{dN}}{\int }_{{\mathbb{R}}^{2dN}}\mathrm{d}\mathbf{q}\mathrm{d}\mathbf{p}\sum _{\mathbf{\ell }\in \{0,1{\}}^N}\pi (\tilde{\mathbf{z}})W_N[A](\tilde{\mathbf{z}}),$$

$$\pi (\tilde{\mathbf{z}})=\frac{1}{{\mathcal{Z}}_N}\exp (-{\beta }_N{H}_N(\tilde{\mathbf{z}})).$$

$$\frac{1}{(2\pi {)}^N}{\int }_{{\mathbb{R}}^{2dN}}\mathrm{d}\mathbf{q}\mathrm{d}\mathbf{p}\sum _{\mathbf{\ell }\in \{0,1{\}}^N}\pi (\tilde{\mathbf{z}})=1.$$

$$⟨\hat{A}⟩\approx \underset{T\to \infty }{\lim }\frac{1}{T}{\int }_0^T{W}_N[A](\tilde{\mathbf{z}}(t))\mathrm{d}t.$$

$$\left\{\begin{array}[]{l l}\,\mathrm{d}{\boldsymbol{q}}&=\nabla_{\boldsymbol{p}}H_{N}(\boldsymbol{q}(t),\boldsymbol{p}(t),\boldsymbol{\ell}(t))\,\mathrm{d}t,\vspace{0.2cm}\\ \,\mathrm{d}{\boldsymbol{p}}&=-\nabla_{\boldsymbol{q}}H_{N}(\boldsymbol{q}(t),\boldsymbol{p}(t),\boldsymbol{\ell}(t))\,\mathrm{d}t\vspace{0.1cm}\\ &\qquad\qquad-\gamma\boldsymbol{p}\,\mathrm{d}t+\sqrt{2\gamma\beta^{-1}_{N}M}\,\mathrm{d}\boldsymbol{B}.\end{array}\right.$$

$${\nabla }_{\mathbf{p}}{H}_N(\mathbf{q}(t),\mathbf{p}(t),\mathbf{\ell }(t))=\mathbf{p}(t)/M.$$

$$\dot{\mathbf{q}}=\frac{1}{M}\mathbf{p}.$$

$$\mathbb{P}\bigl{(}\boldsymbol{\ell}(t+\delta t)=\boldsymbol{\ell}^{\prime}\mid\boldsymbol{\ell}(t)=\boldsymbol{\ell}\,,\boldsymbol{z}(t)=\boldsymbol{z}\bigr{)}=\\ =\delta_{\boldsymbol{\ell}^{\prime},\boldsymbol{\ell}}+\eta\lambda_{\boldsymbol{\ell}^{\prime},\boldsymbol{\ell}}(\boldsymbol{z})\delta t+o(\delta t).$$

$${\lambda }_{{\mathbf{\ell }}^{\prime },\mathbf{\ell }}(\mathbf{z})=\{\begin{array}{ll}-\sum _{\tilde{\mathbf{\ell }}\in S_{\mathbf{\ell }}}{p}_{\tilde{\mathbf{\ell }},\mathbf{\ell }}(\mathbf{z}),& {\mathbf{\ell }}^{\prime }=\mathbf{\ell },\\ p_{{\mathbf{\ell }}^{\prime },\mathbf{\ell }}(\mathbf{z}),& {\mathbf{\ell }}^{\prime }\in S_{\mathbf{\ell }},\\ 0,& \text{otherwise},\end{array}$$

$$p_{\boldsymbol{\ell}^{\prime},\boldsymbol{\ell}}(\boldsymbol{z})=\exp\left(\frac{\beta_{N}}{2}\bigl{(}H_{N}(\boldsymbol{z},\boldsymbol{\ell})-H_{N}(\boldsymbol{z},\boldsymbol{\ell}^{\prime})\bigr{)}\right),$$

$$p_{\boldsymbol{\ell}^{\prime},\boldsymbol{\ell}}(\boldsymbol{z})e^{-\beta_{N}H_{N}(\boldsymbol{z},\boldsymbol{\ell})}=e^{-\frac{\beta_{N}}{2}\bigl{(}H_{N}(\boldsymbol{z},\boldsymbol{\ell})+H_{N}(\boldsymbol{z},\boldsymbol{\ell}^{\prime})\bigr{)}}\\ =p_{\boldsymbol{\ell},\boldsymbol{\ell}^{\prime}}(\boldsymbol{z})e^{-\beta_{N}H_{N}(\boldsymbol{z},\boldsymbol{\ell}^{\prime})}.$$

$$\frac{\partial}{\partial t}f(t,\boldsymbol{z},\boldsymbol{\ell})=\bigl{\{}H_{N},f(t,\boldsymbol{z},\boldsymbol{\ell})\bigr{\}}_{\boldsymbol{z}}-\gamma\nabla_{\boldsymbol{p}}\cdot\Bigl{(}\boldsymbol{p}f+\frac{M}{\beta_{N}}\nabla_{\boldsymbol{p}}f\Bigr{)}\\ +\sum_{\boldsymbol{\ell}^{\prime}}\lambda_{\boldsymbol{\ell},\boldsymbol{\ell}^{\prime}}(\boldsymbol{z})f(t,\boldsymbol{z},\boldsymbol{\ell}^{\prime}),$$

$$\{h,f{\}}_{\mathbf{z}}={\nabla }_{\mathbf{q}}h\cdot {\nabla }_{\mathbf{p}}f-{\nabla }_{\mathbf{p}}h\cdot {\nabla }_{\mathbf{q}}f.$$

$$\bigl{\{}H_{N}(\boldsymbol{z},\boldsymbol{\ell}),e^{-\beta_{N}H_{N}(\boldsymbol{z},\boldsymbol{\ell})}\bigr{\}}_{\boldsymbol{z}}=0.$$

$$\nabla_{\boldsymbol{p}}\cdot\Bigl{(}\boldsymbol{p}e^{-\beta_{N}H_{N}(\boldsymbol{z},\boldsymbol{\ell})}+\frac{M}{\beta_{N}}\nabla_{\boldsymbol{p}}e^{-\beta_{N}H_{N}(\boldsymbol{z},\boldsymbol{\ell})}\Bigr{)}=0.$$

$$\begin{gathered} \sum _{{\mathbf{\ell }}^{\prime }}{\lambda }_{\mathbf{\ell },{\mathbf{\ell }}^{\prime }}(\mathbf{z})e^{-{\beta }_N{H}_N(\mathbf{z},{\mathbf{\ell }}^{\prime })} \\ =\sum _{{\mathbf{\ell }}^{\prime }\in S_{\mathbf{\ell }}}{p}_{\mathbf{\ell },{\mathbf{\ell }}^{\prime }}{e}^{-{\beta }_N{H}_N(\mathbf{z},{\mathbf{\ell }}^{\prime })}-\sum _{{\mathbf{\ell }}^{\prime }\in S_{\mathbf{\ell }}}{p}_{{\mathbf{\ell }}^{\prime },\mathbf{\ell }}{e}^{-{\beta }_N{H}_N(\mathbf{z},\mathbf{\ell })}=0. \end{gathered}$$

$$\widehat{H}=\frac{1}{2M}\left[\begin{array}[]{c c}\widehat{p}^{2}&\\ &\widehat{p}^{2}\end{array}\right]+\left[\begin{array}[]{c c}V_{00}(\widehat{q})&V_{01}(\widehat{q})\\ V_{10}(\widehat{q})&V_{11}(\widehat{q})\end{array}\right]=:\widehat{T}+\widehat{V}.$$

$$0<{\beta }_N<2{\beta }_N<\cdots <(N-1){\beta }_N<N{\beta }_N=\beta .$$

$${\mathrm{Tr}}_{ne}[e^{-\beta \hat{H}}]$$

$$={\mathrm{Tr}}_{{\mathbb{C}}^2}{\int }_{{\mathbb{R}}^d}\mathrm{d}{q}_1⟨q_1|e^{-{\beta }_N\hat{H}}\cdots e^{-{\beta }_N\hat{H}}|q_1⟩$$

$$={\mathrm{Tr}}_{{\mathbb{C}}^2}{\int }_{{\mathbb{R}}^{dN}}\mathrm{d}\mathbf{q}\prod _{k=1}^N⟨q_k|e^{-{\beta }_N\hat{H}}|q_{k+1}⟩,$$

Full text

I Introduction

The development of efficient methods to simulate complex chemical systems on the quantum level has been a central challenge in theoretical and computational chemistry. Exact quantum simulation of coupled nuclear and electronic system is numerically formidable even for small molecules, and, therefore, approximate methods are needed with a reasonable balance of computational effort and incorporating some quantum mechanical aspects of the dynamics.

Under the classical Born-Oppenheimer approximation, one can separate the degrees of freedom associated to the nuclei and electrons, so that a Hamiltonian operator consists of kinetic term and a potential energy surface can be obtained for the nuclear degrees of freedoms. However, when the nonadiabatic effect can not be neglected (often referred as the regime of beyond Born-Oppenheimer dynamics), we need to explicitly include multi-level electronic states in the Hamiltonian, and thus more than one energy surfaces corresponding to different electronic states have to be incorporated. We refer the readers to the reviews Makri (1999); Stock and Thoss (2005); Kapral (2006) for general discussions on simulation methods in the nonadiabatic regime. In this work, we focus on the thermal averages like $\operatorname{Tr}[e^{-\beta\widehat{H}}\widehat{A}]$, where $\widehat{H}$ is a multi-level Hamiltonian operator, $\widehat{A}$ is an observable, and $\beta$ is the inverse temperature.

For the thermal average calculation, the ring polymer representation, based on the imaginary time path integral, has been a popular approach to map a quantum particle in thermal equilibrium to a fictitious classical ring polymer in copies of the phase space Feynman (1972). The representation is asymptotically exact as the number of beads in the ring polymer goes to infinity. Based on the ring polymer representation, path integral Monte Carlo (PIMC) Chandler and Wolynes (1981); Berne and Thirumalai (1986) and path integral molecular dynamics (PIMD) Markland and Manolopoulos (2008); Ceriotti et al. (2010) sampling techniques have been developed to calculate the quantum statistical average. The ring polymer representation has also been used in the dynamics simulations, such as, the centroid molecular dynamics Cao and Voth (1994a, b); Jang and Voth (1999), the ring polymer molecular dynamics Craig and Manolopoulos (2004); Habershon et al. (2013), Matsubara dynamics Hele et al. (2015a, b), and path integral Liouville dynamics Liu (2014); Liu and Zhang (2016).

The conventional ring polymer representation however only works in the adiabatic regime. To apply methods like path-integral molecular dynamics (or dynamic extensions like ring-polymer molecular dynamics) to multi-level systems when the nonadiabatic effects cannot be neglected, a popular strategy is to use the mapping variable approach Meyer and Miller (1979); Stock and Thoss (1997), see also the review article Stock and Thoss (2005) and more recent developments in Ananth and Miller III (2010); Richardson and Thoss (2013); Ananth (2013); Menzeleev et al. (2014); Cotton and Miller (2015); Hele and Ananth (2016); Liu (2016). The idea is to replace the multi-level system by a single level system with higher dimension by mapping the discrete electronic states to continuous variables using uncoupled harmonic oscillators Stock and Thoss (1997). The ring polymer representation can then be applied to the mapped system Ananth and Miller III (2010); Richardson and Thoss (2013); Ananth (2013); Menzeleev et al. (2014); Hele and Ananth (2016).

In this work, we consider rather an alternative strategy for extending the ring polymer representation to multi-level systems, following the spirit of the pioneering work of Schmidt and Tully Schmidt and Tully (2007). In the proposed representation, each bead in the ring polymer is associated with a surface index indicating which energy surface the bead lies on. The total Hamiltonian, and thus the sampling, is given in the extended phase space of position, momentum and surface index of each bead in a ring polymer. While Schmidt and Tully (2007) uses the adiabatic picture, the idea can be generalized to other basis, and this work uses the diabatic picture (which also recovers the ring polymer representation in Schmidt and Tully (2007) in adiabatic picture).

As another main contribution of our current work, we propose a path integral molecular dynamics with surface hopping method (PIMD-SH) to efficiently sample the ring polymer representation, where the discrete electronic state is sampled via a consistent surface hopping algorithm coupled with Hamiltonian dynamics of the position and momentum with Langevin thermostat. It is shown that the PIMD-SH method ergodically samples the correct equilibrium distribution on the extended phase space and can thus be used to sample multi-level quantum systems. In addition, effective numerical integrators are proposed for PIMD-SH to simulate the ergodic trajectories. The numerical results validate the proposed PIMD-SH method for thermal equilibrium sampling of multi-level quantum systems.

Compared with the approaches based on mapping variables, the PIMD-SH is more direct and treats the discrete electronic variables explicitly in the sampling. We think this is more advantageous than treating the electronic and nuclear degrees of freedom on the same footing, as it allows us to employ numerical strategies that exploit the scale separation between the nuclei and electrons and also offer more flexibility. As another advantage, since we use a surface hopping type dynamics to treat the discrete electronic states, it is more natural to combine the proposed thermal (imaginary time) sampling method with real time surface hopping dynamics Tully (1990); Hammes-Schiffer and Tully (1994); Tully (1998); Kapral (2006); Shenvi et al. (2009); Barbatti (2011); Subotnik and Shenvi (2011); Subotnik et al. (2016), which is one of the motivations of our development following our recent works in surface hopping dynamics Lu and Zhou (2016a, b). Let us remark that there have been recent works trying to combine the path integral formulation and surface hopping dynamics Shushkov et al. (2012), though it is unclear if the trajectory with hopping dynamics can preserve the thermal equilibrium. See also discussions on equilibrium properties of the fewest switches surface hopping dynamics in Schmidt et al. (2008); Sherman and Corcelli (2015).

The paper is outlined as follows. In Section II, we present the ring polymer representations for the thermal averages for multi-level quantum systems. The detailed derivation of the ring polymer representation is given together with the proposed PIMD-SH dynamics to sample the equilibrium distribution on the extended phase space. The numerical integration of the PIMD-SH dynamics is discussed in Section III, where we combine a surface hopping dynamics and the Langevin thermostat to treat the electronic states and phase space variables. The numerical tests are presented in Section IV to validate the performance of the PIMD-SH method. In the conclusion and Appendix, we discuss possible future directions and in particular the generalization to more general Hamiltonians.

II Theory

II.1 The ring polymer representation for canonical distribution

for two-level systems

For simplicity of notation, we will restrict the presentation to two-level systems, while the methodology can be generalized to multi-level systems, which will be deferred to Appendix B. In a diabatic representation, the Hamiltonian operator of a general two-level system can be written as (atomic unit is used)

[TABLE]

where $\widehat{q}$ and $\widehat{p}$ are the nuclear position and momentum operators, and $M$ is the mass of the nuclei (for simplicity we assume all nuclei have the same mass). At any position $q\in{\mathbb{R}}^{d}$, the matrix potential

[TABLE]

is a Hermitian matrix which for example comes from the projection of the electronic Hamiltonian to two low-lying states. Here for simplicity, we will assume that the off-diagonal potential functions $V_{01}=V_{10}$ are real valued ($V_{00}$ and $V_{11}$ are real since $V$ is Hermitian). In addition, we will consider the simpler case that the off-diagonal term $V_{01}(q)$ does not change sign for all $q$, the formulation for the general case is discussed in the Appendix. Let us remark that the same simplifying assumptions are made explicitly or implicitly also in the mapping variable approaches (see e.g., Ananth (2013); Menzeleev et al. (2014)).

The Hilbert space of the system is thus $L^{2}({\mathbb{R}}^{d})\otimes{\mathbb{C}}^{2}$, where $d$ is the spatial dimension of the nuclei position degree of freedom. In this work, we consider the thermal equilibrium average of observables, given by

[TABLE]

for an operator $\widehat{A}$, where $\beta=\frac{1}{k_{B}T}$ with $k_{B}$ the Boltzmann constant and $T$ the absolute temperature, and $\operatorname{Tr}_{ne}$ denotes trace with respect to both the nuclear and electronic degrees of freedom, namely,

[TABLE]

The denominator in (1) is the partition function given by $\mathcal{Z}=\operatorname{Tr}_{ne}[e^{-\beta\widehat{H}}]$. For simplicity, throughout this work, we will assume that the observable $\widehat{A}$ only depends on $q$, but not $p$, in other words, $\widehat{A}$ can be written as

[TABLE]

As we will show in Section II.3, for a sufficiently large $N\in\mathbb{N}$, we may approximate the partition function by a ring polymer representation with $N$ beads

[TABLE]

where $\beta_{N}=\beta/N$. The ring polymer that consists of $N$ beads is prescribed by the configuration $(\boldsymbol{q},\boldsymbol{p},\boldsymbol{\ell})\in{\mathbb{R}}^{2dN}\times\{0,1\}^{N}$, where $\boldsymbol{q}=(q_{1},\cdots,q_{N})$ and $\boldsymbol{p}=(p_{1},\cdots,p_{N})$ are the position and momentum of each bead, and $\boldsymbol{\ell}=(\ell_{1},\cdots,\ell_{N})$ indicates the energy level of the bead (thus each bead in the ring polymer lives on two copies of the classical phase space ${\mathbb{R}}^{2d}$, see Figure 1 for an illustration). For a given ring polymer with configuration $(\boldsymbol{q},\boldsymbol{p},\boldsymbol{\ell})$, the effective Hamiltonian $H_{N}(\boldsymbol{q},\boldsymbol{p},\boldsymbol{\ell})$ is given by

[TABLE]

where we take the convention that $\ell_{N+1}=\ell_{1}$ and matrix elements of $G_{k}$, $k=1,\ldots,N$, are given by

[TABLE]

Compared to the conventional ring polymer for a single potential energy surface, the difference in the two-level case is that now each bead is associated with a level index $\ell_{k}$. In particular, when $\ell_{k}\neq\ell_{k+1}$, two consecutive $k$-th and $(k+1)$-th beads in the ring polymer stay on different energy surfaces. We will call this a kink in the ring polymer. Note that the number of kinks is always even since the beads form a ring. In Figure 1, a schematic plot of a ring polymer of $8$ beads with two kinks is shown, where $5$ beads are on the upper energy surface and $3$ are on the lower energy surface.

Moreover, notice that if the off-diagonal terms of the matrix potential $V_{01}=0$, the diagonal part of $G_{k}$ falls back to the usual term in standard ring polymer representation; the current representation is thus a natural extension to the multi-level case.

For an observable $\widehat{A}$, under the ring polymer representation, we have

[TABLE]

where the weight function associated to the observable is given by (recall that $\widehat{A}$ only depends on position by our assumption)

[TABLE]

where we have introduced the short hand notation $\bar{\ell}_{k}=1-\ell_{k}$, i.e., $\bar{\ell}_{k}$ is the level index of the other potential energy surface than the one corresponds to $\ell_{k}$ in our two-level case. Similar as for the partition function, the ring polymer representation (5) replaces the quantum thermal average by an average over ring polymer configurations on the extended phase space ${\mathbb{R}}^{2dN}\times\{0,1\}^{N}$, which consists of not only the position and momentum of the ring polymer, but also the the level index of each bead. We will make precise the accuracy of the approximation below in §II.3.

Comparing (5) and (2), one observes that the partition function (2) under ring polymer representation can be viewed as the thermal average with respect to the weight function $W_{N}[I]$, where $I$ is the identity operator. Note that, the kinks always show up in pairs, so a direct verification shows that $W_{N}[I]=1$.

Next, we shall construct path integral molecular dynamics with surface hopping (PIMD-SH) method to sample the thermal average based on the above formulations in §II.2. The derivation of the ring polymer representation is given in §II.3.

II.2 PIMD-SH method

We observe that (5) can be viewed as (up to a normalization) an average with respect to the classical Gibbs distribution for ring polymers on the extended phase space with Hamiltonian $H_{N}$:

[TABLE]

with distribution

[TABLE]

To simplify the notation, we have denoted by $\widetilde{\boldsymbol{z}}=(\boldsymbol{z},\boldsymbol{\ell})\in{\mathbb{R}}^{2dN}\times\{0,1\}^{N}$ a state vector on the extended phase space, where $\boldsymbol{z}=(\boldsymbol{q},\boldsymbol{p})$ are the position and momentum variables. Notice that in (8), $\mathcal{Z}_{N}$ introduced in (2) normalizes the distribution in the sense that

[TABLE]

As a result, if we can construct a trajectory $\widetilde{\boldsymbol{z}}(t)$ that is ergodic with respect to the equilibrium distribution $\pi$, we can sample the ensemble average on the right hand side of (7) by a time average to approximate $\langle\widehat{A}\rangle$:

[TABLE]

This is the basis of our path-integral molecular dynamics with surface hopping (PIMD-SH) method.

The dynamics of $\widetilde{\boldsymbol{z}}(t)$ is constructed as follows. The position and momentum part of the trajectory $\boldsymbol{z}(t)=(\boldsymbol{q}(t),\boldsymbol{p}(t))$ evolves according to a Langevin dynamics with Hamiltonian $H_{N}(\boldsymbol{q}(t),\boldsymbol{p}(t),\boldsymbol{\ell}(t))$ given the surface index $\boldsymbol{\ell}(t)$, i.e., a Langevin thermostat is used. More specifically, we have

[TABLE]

Here $\boldsymbol{B}=\boldsymbol{B}(t)$ is a vector of $dN$ independent Brownian motion (thus the derivative of each component is an independent white noise), and $\gamma\in\mathbb{R}_{+}$ denotes the friction constant, as usual in Langevin dynamics. Notice that for $\forall\,\boldsymbol{\ell}(t)\in\{0,1\}^{N}$,

[TABLE]

Thus the evolution of the position just follows as usual

[TABLE]

The force term on the hand, given by $-\nabla_{\boldsymbol{q}}H_{N}(\boldsymbol{q}(t),\boldsymbol{p}(t),\boldsymbol{\ell}(t))$, is in general $\boldsymbol{\ell}(t)$-dependent, as the potential energy landscape depends on the level index $\boldsymbol{\ell}(t)$. Hence, the evolution of the momentum $\boldsymbol{p}(t)$ depends on the level index.

The evolution of $\boldsymbol{\ell}(t)$ follows a surface hopping type dynamics in the spirit of the fewest switches surface hopping Tully (1990) (see also our recent works Lu and Zhou (2016a, b)). In particular, we can take it to be a Markov jump process with infinitesimal transition rate over the time period $(t,t+\delta t)$ for $\delta t\ll 1$ given by

[TABLE]

This means that if the current configuration of the ring polymer is given by $(\boldsymbol{z},\boldsymbol{\ell})$, during the time interval $(t,t+\delta t)$, the level index might change to $\boldsymbol{\ell}^{\prime}\neq\boldsymbol{\ell}$ with probability $\eta\lambda_{\boldsymbol{\ell}^{\prime},\boldsymbol{\ell}}(\boldsymbol{z})\delta t+o(\delta t)$. Here $\eta>0$ is an overall scaling parameter for hopping intensity (the larger $\eta$ is, the more frequent hopping occurs), the coefficients $\lambda_{\boldsymbol{\ell}^{\prime},\boldsymbol{\ell}}$ are specified as

[TABLE]

where $S_{\boldsymbol{\ell}}=\{\boldsymbol{\ell}^{\prime}\mid\lVert\boldsymbol{\ell}^{\prime}-\boldsymbol{\ell}\rVert_{1}=1\;\mbox{or}\;\boldsymbol{\ell}^{\prime}=\boldsymbol{1}-\boldsymbol{\ell}\}$ denotes all allowed configuration $\boldsymbol{\ell}^{\prime}$ after the hopping: $\boldsymbol{1}$ is the vector with all entries $1$, so $\boldsymbol{\ell^{\prime}}=\boldsymbol{1}-\boldsymbol{\ell}$ indicates that the surface index of each bead is flipped; and $\lVert\boldsymbol{\ell^{\prime}}-\boldsymbol{\ell}\rVert_{1}=\sum_{k}\lvert\ell^{\prime}_{k}-\ell_{k}\rvert=1$ indicates that one and only one bead jumps to the opposite energy surface. Here in the rate expression, $p_{\boldsymbol{\ell}^{\prime},\boldsymbol{\ell}}(\boldsymbol{z})$ is defined as

[TABLE]

which is chosen so that the detailed balance relation is satisfied

[TABLE]

This guarantees that the distribution $\pi$ is preserved under the dynamics of the jumping process (as will be further discussed below).

The above choice of $S_{\boldsymbol{\ell}}$ allows only two types of change of level indices: either changing the surface index of one single bead (single hop) or changing the surface index of all beads (total flip). This is chosen for simplicity, as it ensures that any surface index configuration can be reached and at the same time we do not need to consider all possibilities at a single time step, which is combinatorial and inefficient for practical implementation.

To show that $\pi$ as in (8) is indeed the equilibrium distribution corresponds to the dynamics of $\widetilde{\boldsymbol{z}}(t)$, it is more convenient to write down the associated Fokker-Planck equation of the dynamics. Denote the probability distribution on the extended phase space at time $t$ by $f(t,\widetilde{\boldsymbol{z}})$, $f$ satisfies the following Fokker-Planck equation (i.e., forward Kolmogorov equation)

[TABLE]

where on the right hand side, the first term accounts for the Hamiltonian part of the dynamics of $\boldsymbol{z}$, the second term comes from the dissipation and fluctuation due to the Langevin thermostat, and the last term is due to the jumping process of $\boldsymbol{\ell}$. In the above equation, $\{\cdot,\cdot\}_{\boldsymbol{z}}$ stands for the usual Poisson bracket corresponding to the Hamiltonian dynamics

[TABLE]

Since the Boltzmann-Gibbs distribution is proportional to $e^{-\beta_{N}H_{N}}$, which is a function of the Hamiltonian, it follows that

[TABLE]

As usual for Langevin dynamics, the fluctuation-dissipation balance ensures that

[TABLE]

Moreover, using the detailed balance relation (12), we verify using the definition of $\lambda$ as in (11) that

[TABLE]

Therefore, we conclude that the Boltzmann-Gibbs distribution, which is a constant multiple of $e^{-\beta_{N}H_{N}(\boldsymbol{z},\boldsymbol{\ell})}$, is an stationary solution to the Fokker-Planck equation. Note that this remains the case regardless the choice of the hopping intensity parameter $\eta$ and the friction parameter $\gamma$ (the latter is of course a familiar fact for Langevin thermostat). We will study the effect of tuning $\eta$ in our numerical examples in §IV.4.

II.3 Ring polymer representations for the two-level Hamiltonians

Let us now present the derivation of the ring polymer representation for the two-level system, as presented in §II.1. Recall the Hamiltonian in the diabatic picture is given by

[TABLE]

We recall here for convenience that we have assumed $\widehat{V}$ is real symmetric, and the off-diagonal function, which is denoted by $v$ for simplicity, does not change sign for $q\in\mathbb{R}^{d}$. We consider a large fixed $N\in\mathbb{N}^{+}$ and introduce an equispaced partition of $[0,\beta]$ (recall that $\beta_{N}=\beta/N$),

[TABLE]

By inserting resolution of identities with respect to position, we have for the partition function

[TABLE]

where we have used the convention $q_{N+1}=q_{1}$ to simplify the expression. So far, the reformulation is exact. Applying the Strang splitting Strang (1968) to the short imaginary time propagator $e^{-\beta_{N}\widehat{H}}$ and inserting a resolution of identity with respect to momentum, we get

[TABLE]

Direct calculation of the right hand side leads to

[TABLE]

where we suppress in the notation the dependence of $S_{k}$ on $p_{k},q_{k}$ and $q_{k+1}$:

[TABLE]

Note that, if the $V$ matrices are diagonal, we have recovered the familiar terms in the usual ring polymer representation for single energy surface. Let us emphasize though we consider the general case that $V$ contains off-diagonal terms, and in particular, $\exp(-\frac{\beta_{N}}{2}V)$ is a $2\times 2$ matrix and hence $\exp(-\frac{\beta_{N}}{2}V(q_{k}))$ in general does not commute with $\exp(-\frac{\beta_{N}}{2}V(q_{k+1}))$.

Applying the above calculation to every term in the product on the right hand side of (14), we obtain

[TABLE]

We have an extra $N$ factor in the error $\mathcal{O}(N\beta_{N}^{3})$ since the error from each operator splitting adds up.

So far the basis for the discrete electronic states has not been fixed. In below, we will choose to use the diabatic picture, while the adiabatic picture can be also used, as the calculations shown in Appendix C. In the adiabatic picture, we recover the ring polymer representation of Schmidt and Tully in Schmidt and Tully (2007). The sampling strategy based on PIMD-SH can be applied to adiabatic picture as well, which we will leave for future studies.

Let us further simplify the integrand in the above equation in order to arrive at the desired Boltzmann-Gibbs form as in (2). In particular, to deal with the matrix exponential $e^{-\beta_{N}V}$, we split the matrix potential $V$ into diagonal and off-diagonal parts

[TABLE]

Using another Strang splitting, we obtain

[TABLE]

Explicit calculation gives the matrix exponential

[TABLE]

Substitute this into the right hand side of (16) and rewrite the resulting matrix elements as exponentials, we arrive at

[TABLE]

where the prefactor $-V_{01}/\lvert V_{01}\rvert$ when $\ell\neq\ell^{\prime}$ is due to dependence of the sign of $\sinh(-\beta_{N}V_{01})$ on the positivity of $V_{01}$.

Applying the above expression for $e^{-\beta_{N}V}$ to (15), we arrive at the following approximation (recall the definition of matrices $G_{k}$ in (4))

[TABLE]

which justifies the ring polymer representation of the partition function (2). In order to get the above, we realize that the $-V_{01}/\lvert V_{01}\rvert$ factors for all kinks (for $k$ such that $\ell_{k}\neq\ell_{k+1}$) in a ring polymer will cancel since there are even number of kinks and $V_{01}$ does not change sign by assumption (when $V_{01}$ is always negative, the factors are simply always $1$; when $V_{01}$ is always positive, the factors are $-1$, even number of them multiply to $1$).

As a result, each term in the average (2) is positive (as it is an exponential), and thus we can view $\exp(-\beta_{N}H_{N}(\boldsymbol{q},\boldsymbol{p},\boldsymbol{\ell}))$ as a probability density for the ring polymer configuration $(\boldsymbol{q},\boldsymbol{p},\boldsymbol{\ell})$. The PIMD-SH then samples this distribution on the extended phase space. This is no longer true without the assumption that the off-diagonal entry of the matrix potential $V_{01}$ does not change sign. While we can still take the absolute value of the summand as the distribution, we would also need to approximate the partition function by an average of terms which change sign depending on the ring polymer configuration. The sign change in general increases the difficulty of the sampling, this is a manifestation of the familiar “fermionic sign problem” in quantum Monte Carlo simulations, see e.g., Ceperley (2010). Further discussions on the formulation for a general two-level system can be found in the Appendix A and will be explored in future works.

For an observable $\widehat{A}$ that only depends on the position variable, following a similar derivation leads to (15), we have

[TABLE]

By symmetry, we can also move the $A$ matrix to before $e^{-\beta_{N}V(q_{k})}$ and evaluate $A$ at $q_{k}$. Taking an average over all the possibilities, we get

[TABLE]

Again using the expansion (16) of $e^{-\beta_{N}V}$, we arrive at (recall that $\bar{\ell}_{k}=1-\bar{\ell}_{k}$)

[TABLE]

where in the first equality, we have used the short hand notations $J_{k}$:

[TABLE]

The above approximation of $\operatorname{Tr}_{ne}[e^{-\beta\widehat{H}}\widehat{A}]$ is exactly (5) by checking the definition of the weight function in (6).

III Time integrator for PIMD-SH

Recall that given the two level quantum Hamiltonian $\widehat{H}$ and $\beta$, we choose a sufficiently large $N\in\mathbb{N}$, and the effective Hamiltonian for the ring polymer representation $H_{N}$ is given in (3). The PIMD-SH dynamics is ergodic with respect to the corresponding Gibbs distribution of the Hamiltonian. In practical implementations, the PIMD-SH dynamics is discretized in time to sample the equilibrium distribution.

To start with, we specify initial conditions to the sampling trajectory $\widetilde{\boldsymbol{z}}(0)=(\boldsymbol{q}(0),\boldsymbol{p}(0),\boldsymbol{\ell}(0))$. Due to the ergodicity of the dynamics, any initial conditions can be used, while a better initial sampling will accelerate the convergence of the sampling. In our current implementation, for simplicity, we initialize all the beads in the same position, sample their momentum by the Gaussian distribution $\mathcal{N}(0,M\beta_{N}^{-1})$, and take $\boldsymbol{\ell}(0)=\boldsymbol{0}$, where $\boldsymbol{0}$ is a vector of zeros, meaning that initially all beads of the ring polymer stay on the lower energy surface. Possible better initial sampling strategies can also be used.

The overall strategy we take for the time integration is time splitting schemes, by carrying out the jumping step, denoted by J, and the Langevin step denoted by L, in an alternating way. In this work, we apply the Strang splitting, such that the resulting splitting scheme is represented by JLJ. This means that, within the time interval $[t^{n},t^{n}+\Delta t]$ ($\Delta t$ being the time step size), we carry out the following steps in order:

We numerically simulate the jumping process for $\boldsymbol{\ell}$ for $\Delta t/2$ time with fixed position and momentum of the ring polymer; 2. 2.

We propagate numerically the position and momentum of the ring polymer using a discretization of the Langevin dynamics for $\Delta t$ time while fixing the surface index $\boldsymbol{\ell}$ (from the previous sub-step); 3. 3.

The jumping process for $\boldsymbol{\ell}$ is simulated for another $\Delta t/2$ time with fixed position and momentum of the ring polymer; 4. 4.

The weight function $W_{N}[A](\widetilde{\boldsymbol{z}}(t^{n+1}))$ of the observable $\widehat{A}$ is calculated (and stored, if needed) to update the running average of the observable.

The above procedure is repeated for each time step until we reach a prescribed total sampling time $T$ or when the convergence of the sampling is achieved under certain stopping criteria. In this work, we use a standard Monte Carlo scheme for the jumping process and the BAOAB integrator for the Langevin dynamics, the details of both will be further elaborated in subsections.

We remark that from a numerical analysis point of view, the above splitting scheme corresponds to a splitting of the Fokker-Planck equation introduced in (13): The jump process step corresponds to

[TABLE]

while the Langevin step corresponds to

[TABLE]

In particular, this leads to error analysis of the weak order of the proposed scheme. We shall not go into the details of numerical analysis in this work.

In what follows, we present the details of each steps in the splitting scheme.

III.1 Simulation of the jumping process

Within a short time interval $\Delta t$, we consider the possible jumps of the surface index $\boldsymbol{\ell}$, while fixing the position and momentum of the ring polymer. For simplicity, we suppress the appearance of $\boldsymbol{q}$ and $\boldsymbol{p}$ in the notation in this subsection, as they do not change.

We assume that the time interval is so short that we may only consider one or no jump during the interval, i.e.,, two jumps happening at probability $\mathcal{O}(\Delta t)^{2}$ can be neglected. Recall the hopping intensity (11) for the jumping process, in particular, by our choice of $S_{\boldsymbol{\ell}}$, only jumps with $\lVert\boldsymbol{\ell}-\boldsymbol{\ell}^{\prime}\rVert_{1}=1$ or when $\boldsymbol{\ell}=\boldsymbol{1}-\boldsymbol{\ell}^{\prime}$ are allowed. Thus, the probability of a jump occurs from the current level index $\boldsymbol{\ell}$ to $\boldsymbol{\ell^{\prime}}\neq\boldsymbol{\ell}$ during the time interval $\Delta t$ is given by

[TABLE]

where recall that $S_{\boldsymbol{\ell}}=\{\boldsymbol{\ell}^{\prime}\mid\lVert\boldsymbol{\ell}^{\prime}-\boldsymbol{\ell}\rVert_{1}=1\;\mbox{or}\;\boldsymbol{\ell}^{\prime}=\boldsymbol{1}-\boldsymbol{\ell}\}$. We assume that $\Delta t$ is chosen sufficiently small that

[TABLE]

and the probability that the level index is unchanged is then $1-\sum_{\boldsymbol{\ell}^{\prime}\neq\boldsymbol{\ell}}h_{\boldsymbol{\ell}^{\prime},\boldsymbol{\ell}}$. Thus a uniform random number on $[0,1]$ can be drawn to decide which event happens (as in standard kinetic Monte Carlo simulations).

Recall that the choice of the hopping intensity parameter $\eta$ does not change the equilibrium distribution of the PIMD-SH. In practical simulations, one can adjust the parameter $\eta$ to make surface hoppings more often, which in certain cases help accelerate the sampling. This is particularly useful when the observable has off-diagonal elements. We shall numerically verify the effect the this parameter in §IV.4. Of course, a larger $\eta$ requires smaller time step size to ensure the condition (22) and also for accuracy. Therefore, direct simulation using a very large $\eta$ might be difficult. On the other hand, we can explore similar strategies as in the infinite swapping replica exchange molecular dynamics Lu and Vanden-Eijnden (2013); Yu et al. (2016) to simulate the limiting system directly. This will be left for future works.

III.2 The BAOAB scheme for Langevin dynamics

With a fixed level index $\boldsymbol{\ell}$, the Langevin step then corresponds to the propagation of the position and momentum $(\boldsymbol{q},\boldsymbol{p})$ according to

[TABLE]

To numerically integrate the Langevin dynamics, we apply the BAOAB splitting scheme Leimkuhler and Matthews (2013) developed in the context of Langevin thermostat for classical molecular dynamics. The BAOAB scheme has been applied to and analyzed for conventional PIMD simulations (for single energy surface) recently and exhibit advantageous numerical performance over other numerical integrators Liu et al. (2016) 111The work Liu et al. (2016) actually considered two slightly different variants of BAOAB style splitting, the BAOAB splitting we use here, which is consistent with Leimkuhler and Matthews (2013), corresponds to what Liu et al. (2016) referred as BAOAB-num.. In the BAOAB scheme, the Langevin dynamics is divided into three parts, the kinetic part (denoted by “A”),

[TABLE]

the potential part (denoted by “B”)

[TABLE]

and the Langevin thermostat part (denoted by “O”)

[TABLE]

The nice feature of such splitting is that each of these substeps can be integrated explicitly. For example, the O part has the following exact solution in the sense of the generated probability distributions

[TABLE]

where each component of $\boldsymbol{n}$ is an independent standard Gaussian random variable.

The BAOAB scheme stands for a splitting scheme solving B part and A part for a half time step $\Delta t/2$, followed by solving O part for a full time step $\Delta t$, and followed by solving A part and B part for another half time step $\Delta t/2$. Each substep uses the exact time propagations. Compared with other prevailing splitting methods for the Lagevin dynamics, the BAOAB scheme has demonstrated higher accuracy without sacrificing computational efficiency. Also, it is also demonstrated Leimkuhler and Matthews (2013) the BAOAB scheme enables using larger time steps in the simulation while keeping the stability of the integrator.

Combining the BAOAB splitting scheme with the step of jumping process, our overall scheme can be represented as JBAOABJ. While other numerical integrators are possible, in our numerical experiments shown in Section IV, the JBAOABJ scheme seems to perform quite well for the PIMD-SH sampling.

IV Numerical results

IV.1 Test examples

To validate the PIMD-SH method and to understand the choice of parameters in the method, we test with following two potentials. Both potentials are chosen to be one-dimensional and periodic over $[-\pi,\pi]$, so that the reference solutions can be obtained accurately with pseudo-spectral approximations and compared to the PIMD-SH results. The first test potential is given by

[TABLE]

Clearly, $V_{11}\geqslant V_{00}$ and the two energy surfaces only intersect at $x=0$, where the off-diagonal term takes its largest value. The energy surfaces are symmetric with respect to $x=0$. At thermal equilibrium, the density are expected to concentrate around $x=0$, where transition between the two bands is the most noticeable due to the larger off-diagonal coupling terms. The diabatic energy surfaces are plotted in Figure 2(top). Moreover, for $\beta=1$, $M=10$ $N=64$, $\Delta t=1\mathrm{e}-3$ and $T=1\mathrm{e}4$, we plot the (numerically obtained) marginal equilibrium distribution of position variable on each energy surface in Figure 2(bottom).

The other test potential we take is given by

[TABLE]

The potential is designed so that $V_{11}\geqslant V_{00}$, and the two energy surfaces achieve their minima around $x=1$. Moreover, the energy surfaces almost intersect when $x$ is slightly less than [math], and the off-diagonal potential is most noticeable around $x=-0.5$. Thus, in this model, the potential is asymmetric, and the location where the equilibrium distribution is concentrated deviates from the most active hopping area. These make this test model more challenging than the previous one. We plot the diabatic energy surfaces in Figure 3(top); and for $\beta=1$, $M=10$, $N=64$, $\Delta t=0.001$ and $T=10000$, we plot the (numerically obtained) marginal equilibrium distribution on each energy surface in Figure 3(bottom).

Note that in both test cases, the off-diagonal potential $V_{01}$ has the same sign for all $q$, as we have assumed for the proposed method in this work. We will study the cases when the observables are diagonal, or only its off-diagonal elements are nonzero. When, the observables are diagonal, i.e., $\langle\ell|A|\ell^{\prime}\rangle=0$, when $\ell\neq\ell^{\prime}$, the weight function of the observable simplifies to

[TABLE]

From the simplified representation of the weight function, we learn that each beads gives similar contribution to the observations.

When the diagonal elements of the observable are zero, i.e., $\langle\ell|A|\ell\rangle=0$ for $\ell=0,1$, the weight function of the observable simplifies to

[TABLE]

From which we see that, due to the exponential factor in the weight function, the kinks of a ring polymer might give larger contributions to the observables compared to other beads. The uneven contribution to the observable from beads brings in numerical challenges since the sample variance may be large. Numerical study of related issues is presented in §IV.4.

IV.2 Convergence with number of beads

We first test the PIMD-SH method for the Hamiltonian with the test potential (24) with $\beta=1$, $M=10$ and $\eta=1$. We carry out the simulations with the following diagonal observable

[TABLE]

for $\beta_{N}=\frac{1}{4}$, $\frac{1}{8}$ and $\frac{1}{16}$ (and correspondingly different number of beads in the ring polymer), $\Delta t=1\rm e-3$ and $T=1\rm e7$. The time step size is chosen sufficiently small to ensure the accuracy of numerical integration. We plot the running average of each simulation in Figure 4, from which we observe that the running average for each $\beta_{N}$ approaches a steady value and is close to the reference solution. We further compute the mean squared errors for different simulation time for $\beta_{N}=\frac{1}{4}$, $\frac{1}{8}$ and $\frac{1}{16}$, and plot them in Figure 5, from which we observe that all three tests behave similarly, and the mean squared errors decay roughly in proportion to $\mathcal{O}(t^{-1})$ until around $t=1\rm e5$, consistent with the convergence rate of the Monte Carlo sampling error. The asymptotic errors become noticeable when $t=1\rm e5\sim 1\rm e7$, which tell apart the different $\beta_{N}$ with also reduced rates for decrease of the error.

The errors in empirical averages together with the $95\%$ confidence intervals and the mean squared errors are shown in Table 1, from which we see that the result approximates the correct expectation of the observable while both the sampling error and the asymptotic error make contributions to the total numerical error. The mean squared errors of the empirical averages are defined as $\text{M.S.E.}=\text{Bias}^{2}+\text{Var}$, where Bias is calculated using the reference value and Var is estimated using the observed data and the effective sampling size. The result further confirms that increasing number of beads help reduce the error in the PIMD-SH method.

The test for different $\beta_{N}$ is also carried out for the first potential, though the conclusion is pretty similar (and the second potential is more challenging), and hence we omit these results here.

IV.3 Convergence with different time step sizes

In this test, we use the test potential (24), with $\beta=1$, $\beta_{N}=\frac{1}{16}$, $M=10$ and $\eta=1$, and vary the time step sizes for the numerical integration. We choose the same diagonal observable as before in (25). The PIMD-SH method is tested with time step sizes $\Delta t=\frac{1}{25}$, $\frac{1}{50}$, $\frac{1}{100}$, $\frac{1}{200}$ and $\frac{1}{400}$ with total simulation time $T=10000$. The running averages of the sampling for various parameters are plotted in Figure 6. The errors in the empirical averages together with their $95\%$ confidence intervals and mean squared errors are shown in Table 2. We observe that the PIMD-SH captures the correct thermal average of the observable even for a relatively large $\Delta t$. While the MSE decays for smaller time step, the decrease is not too significant, as the sampling error is probably dominant. We also plot the empirical histogram for the spatial distribution on each energy level in Figure 7 for two choices of the time step sizes.

We also carry out the $\Delta t$ test for the test potential (23) with an off-diagonal observable:

[TABLE]

We take $\beta=1$, $M=10$ and $\eta=5$, and test time step sizes $\Delta t=\frac{1}{25}$, $\frac{1}{50}$, $\frac{1}{100}$, $\frac{1}{200}$ and $\frac{1}{400}$ till $T=10000$. The numerical results are plotted in Figure 6. The errors in the empirical averages together with their $95\%$ confidence intervals and mean squared errors are shown in Table 3. In this test, when $\Delta t$ is large, the numerical error seems to be dominated by error in numerical integration of the trajectory, and reducing $\Delta t$ reduces the error of the sampling. Compared with the previous case with diagonal observables, the variance of the observable is larger even when a small time step size is used; sampling the thermal average for observables with off-diagonal entries is more challenging. This is easy to understand as for off-diagonal observables, the ring polymer has to have a kink to make the weight function non-zero, and hence the surface hopping dynamics becomes more important. As we will see in the next subsection, increasing the hopping parameter $\eta$ in the PIMD-SH can improve the results.

IV.4 Effect of the hopping intensity parameter $\eta$

In this test, we implement the PIMD-SH method with the test potential (23), with $\beta=1$, $\beta_{N}=\frac{1}{16}$ and $M=10$. We choose the off-diagonal observable (26) as before.

As we mentioned before in §IV.1 and also have seen in the numerical results in Figures 6 and 8, sampling of the off-diagonal observable is more challenging in PIMD-SH due to the contribution to the weight by ring polymer configurations with kinks and also the sampling variance is larger. In this test, we show that increasing the hopping intensity parameter $\eta$, which makes hopping more frequent, helps sampling off-diagonal observables. For $\eta=5$, $10$, $20$ and $40$, we test the PIMD-SH method with $\Delta t=(100\eta)^{-1}$ respectively till $T=10000$. Recall that a small time step size is required to maintain the accuracy and stability of the integrator.

The results are plotted in Figure 9 and the errors in the empirical averages together with their $95\%$ confidence intervals and mean squared errors are shown in Table 4. Increasing the hopping intensity parameter can effectively reduce the sampling error and variance. We admit though it is computationally more expensive to use a larger $\eta$ since the time step size has to be smaller by directly applying the JBAOABJ scheme. A better numerical scheme in the spirit of Lu and Vanden-Eijnden (2013); Yu et al. (2016) is needed and will be leaved for future works.

V Conclusion

We have proposed in this work the PIMD-SH method for sampling the thermal equilibrium average of multi-level quantum systems. The formulation is justified theoretically and supported by numerical results.

Among the possible future directions based on the current work, the most interesting direction is perhaps to combine this approach with surface hopping method for sampling dynamical correlation functions of the type $\operatorname{Tr}_{ne}[e^{-\beta\widehat{H}}\widehat{A}(t)\widehat{B}]$, for which the ring polymer representation we developed in this work is well suited.

Better numerical integration strategies especially for using a larger hopping intensity parameter $\eta$ is worth exploring. Further investigation is also needed when the off-diagonal components of the potential function changes sign or takes complex values.

Acknowledgements.

This work is partially supported by the National Science Foundation under grant DMS-1454939. J.L. would like to thank Giovanni Ciccotti and John Tully for helpful discussions.

Appendix A Ring polymer representations for a general two-level matrix potential

We have assumed that the off-diagonal terms of the matrix potential in a two-level system $V_{01}=V_{10}$ is real and does not change sign. In this Appendix, we discuss the formulations for a general two-level system. Since the potential matrix is Hermitian, the diagonal potential terms $V_{00}$ and $V_{11}$ are always real, while $V_{01}$ and $V_{10}$ may be complex with $V_{10}=V_{01}^{\ast}$.

If we repeat the derivation in Section II.3 for the general case, everything is parallel till the step when we approximate $e^{-\beta_{N}V}$ with the Strang splitting in (16). In the general case, the off diagonal part of the matrix changes to

[TABLE]

Using again the Strang splitting and the explicit expression of $\exp(-\beta_{N}V_{o})$, we get (cf. (17))

[TABLE]

Different from the case considered in the main text, since $V_{\ell\ell^{\prime}}$ is complex, in general the phase factor does not cancel for all the kinks in the ring polymer. We can still approximate the partition function as an average over the extended phase space, though the quantities to be averages is now complex (in the case where the off-diagonal entry of $V$ is real but changes sign, we will need to average over terms with different signs). If we collect of the phase factors in a single term, we arrive at the following approximation

[TABLE]

where the weight factor becomes

[TABLE]

and $G_{k}$, $J_{k}$ are still defined as before in (4) and (20) respectively. Recall that in the previous case, the weight factor for the partition function, corresponding to the identity operator, is just constant $1$. As a result, the thermal average of a general observable becomes a ratio of ensemble averages

[TABLE]

where the expression of $W_{N}[A]$ is same as before in (5) and the equilibrium distribution is given in (8). This can still be sampled using a PIMD-SH dynamics associated with the Gibbs distribution $\pi$ on ${\mathbb{R}}^{2dN}\times\{0,1\}^{N}$, though some bias will be introduced if we use the same ergodic trajectory to sample both the numerator and denominator.

Another consequence of the general phase of the term $-V_{01}/\lvert V_{01}\rvert$ is that the partition function is determined through averaging of terms that change signs / phases, which is reminiscent of the “sign problems” in quantum Monte Carlo simulations. Further investigations of such issues are needed.

Let us remark that there are other possibilities of choosing the effective Hamiltonian and the corresponding weight functions. For example, after the splitting of $e^{-\beta_{N}V}$, instead rewrite the $\sinh$ and $\cosh$ functions to the exponent, we can put these terms to the weight functions. We choose the present approach since in the case that the off-diagonal matrix potentials are real and do not change sign, it reduces to a nice probabilistic sampling problem; but for general cases, it is worth considering other approaches.

Appendix B Ring polymer representations for $M$-level system

($M\geqslant 3$)

We show the extension of the ring polymer representation to a general $M$-level system ($M\geqslant 2$) in this Appendix. For simplicity of the presentation, let us just focus on the expression for the partition function, cf. the derivation in Section II.3; the expression for the observable follows analogously. For the $M$-level case, The Hamiltonian is given by

[TABLE]

where the kinetic operator is still diagonal, and the potential matrix is a $M\times M$ Hermitian matrix such that $V_{ij}=V_{ji}^{*},\,i,j=1,\cdots,M$. We decompose the matrix potential potential as follows

[TABLE]

where $V_{d}$ is the diagonal part of $V$ and the off-diagonal part has been decomposed pairwisely into a sum of $V^{(i,j)}$ with elements given by

[TABLE]

Note that in reformulating the partition function, one can still apply the Strang splitting to each $e^{-\beta_{N}\widehat{H}}$ and obtain (15). Next, in order to derive an approximate formula for matrix elements of $e^{-\beta_{N}V}$, we applied the Strang splitting multiple times and get

[TABLE]

where the last line defines $M_{v}$, the approximate matrix for $e^{-\beta_{N}V}$. Here, for the diagonal part,

[TABLE]

And for the off-diagonal part, we have

[TABLE]

Next, we observe that

[TABLE]

where $\boldsymbol{n}=(n_{(M,M-1)},\cdots,n_{(2,1)},n_{(2,1)}^{\prime},\cdots,n_{(M,M-1)}^{\prime})$ with each entry takes possible values in $\{1,\ldots,M\}$. We introduce the augmented index $\overline{\boldsymbol{n}}=\{\ell,\boldsymbol{n},\ell^{\prime}\}$, and observe that when two consecutive index in $\overline{\boldsymbol{n}}$ are different, the product in (32) is either [math] or gains a multiplier of order $\mathcal{O}(\beta_{N})$. Therefore, if we omit all the terms which are of order $\mathcal{O}(\beta_{N}^{3})$, we conclude that, when $\ell^{\prime}=\ell$,

[TABLE]

where

[TABLE]

By the Taylor’s expansion, we can further simplify and obtain

[TABLE]

Similarly, for $\ell^{\prime}\neq\ell$, if we omit all the terms which are of order $\mathcal{O}(\beta_{N}^{3})$ in (32), we obtain that

[TABLE]

Note that the two expressions above are natural generalization of the two level counterpart (17).

Appendix C Ring polymer representations for $2$-level system in the adiabatic picture

In this section, we instead use the adiabatic basis to handle discrete electronic states, and we are able to recover the results by Schmidt and Tully in Schmidt and Tully (2007). We start with (15), where the electronic states have not been specified yet. We denote the adiabatic states by $\Phi_{n}(q)$ and the adiabatic surface by $E_{n}(q)$, where $n=0,1$, and they satisfy

[TABLE]

By inserting multiple times the resolution of identity $I=\sum_{n=0,1}\lvert\Phi_{n}(q)\rangle\langle\Phi_{n}(q)\rvert$, we get

[TABLE]

As $\Phi_{n}$ are eigenfunctions of $V$, $e^{-\beta_{N}V(q)}\Phi_{n}(q)=e^{-\beta_{N}E_{n}(q)}\Phi_{n}(q)$, the above can be further simplified as

[TABLE]

If we define the short hand

[TABLE]

we obtain

[TABLE]

with the effective Hamiltonian given by

[TABLE]

This is exactly the ring polymer representation derived in Schmidt and Tully (2007). The PIMD-SH sampling method can be applied to the adiabatic picture as well, which we will consider in future works.