Modeling the dynamics of photophysical and (photo)chemical reactions in extended molecular systems is a new frontier for quantum chemistry. Many dynamical phenomena, such as intersystem crossing, nonradiative relaxation, and charge and energy transfer, require a nonadiabatic description which incorporate transitions between electronic states. Additionally, these dynamics are often highly sensitive to quantum coherences and interference effects. Several methods exist to simulate nonadiabatic dynamics; however, they are typically either too expensive to be applied to large molecular systems (10's100's of atoms), or they are based on ad hoc schemes which may include severe approximations due to inconsistencies in classical and quantum mechanics. We present, in detail, an algorithm based on Monte Carlo sampling of the semiclassical timedependent wavefunction that involves running simple surface hopping dynamics, followed by a postprocessing step which adds little cost. The method requires only a few quantities from quantum chemistry calculations, can systematically be improved, and provides excellent agreement with exact quantum mechanical results. Here we show excellent agreement with exact solutions for scattering results of standard test problems. Additionally, we find that convergence of the wavefunction is controlled by complex valued phase factors, the size of the nonadiabatic coupling region, and the choice of sampling function. These results help in determining the range of applicability of the method, and provide a starting point for further improvement.
I. INTRODUCTION
Advances in theoretical methods as well as in computational power have positioned quantum chemistry as a powerful tool in studying various properties of molecules and molecular complexes.^{1,2} Accurate, but numerically expensive, wavefunction based approaches enable detailed description of electronic properties in relatively small molecules. In contrast, efficient Density functional theory (DFT) methods in combination with molecular dynamics (MD) allow one to deduce valuable information on ground state properties for molecules with hundreds of atoms in size.^{3} Further advances in timedependent density functional theory (TDDFT) have made it possible to carry out efficient calculations for exited states in these systems and thus predict their susceptibilities, absorption, and emission spectra, etc.^{4}
At present quantum chemistry faces a new frontier: Now it aims not only at computations of the equilibrium (static) molecular properties, but also at modeling dynamics of photophysical processes and (photo)chemical reactions in molecular systems.^{5,6} The latter are characterized by different reaction pathways that lead to different reaction products. For example, upon absorbing a photon, a molecule may undergo radiative relaxation processes,^{7} such as fluorescence,^{8} and phosphorescence; nonradiative relaxation processes, such as photodissociation^{9–11} or photoisomerization;^{12–15} charge;^{16–19} and energy transfer;^{20–22} or intersystem crossing.^{23–25} Such dynamics are obviously accompanied by electronic transitions, involving complex electronphonon interaction^{26,27} and possibly many electronic states, and therefore cannot be described within the standard adiabatic or BornOppenheimer approximation. That is, while in the traditional adiabatic approximation the nuclei move along a given potential energy surface (PES) of a molecule, description of photophysical or photochemical processes requires accounting for the transitions between different PESs. Such transitions typically occur in the relatively small regions of the phase space where the PESs closely approach or cross each other, so that the energy separation between relevant PESs becomes comparable with inverse time scales of the nuclear motion (in units of ℏ), i.e., phonon frequencies. Thus the assumption of separation between electronic and nuclear timescales breaks down in the vicinity of electronic PES crossings and dynamics becomes nonadiabatic.
Several approaches exist in the literature to treat the nonadiabatic dynamics within the framework of quantum chemistry.^{28–36} The most well known of these are mixed quantum classical treatments, the Ehrenfest and surface hopping methods. The Ehrenfest method is easy to implement. It requires one to run a single trajectory of the nuclei, on an average PES, with electronic populations being evaluated on the fly.^{37} The simplicity of the implementation of the method is flawed by its poor accuracy, related to its mean field nature. Usually different electronic configurations are associated (“entangled”) with different paths of the nuclei and so the mean field approach becomes inadequate when these paths are very dissimilar. In order to fix this problem, Tully introduced a surface hopping approach, which accounts for the trajectory branching by running multiple trajectories that can hop between PESs.^{38} Fewest switches surface hopping (FSSH) has become the most popular approach for problems at various scales, from molecules to nanoclusters.^{39–45} However, the method is build on ad hoc assumptions and, as a result, frequently fails to properly describe correlations between the trajectories of nuclei and electronic states. In particularly, it does not take into account decoherence arising due to spacial separation between the components of the nuclear wavefunctions corresponding to different electronic states (PES) and, as a results, does not properly accounts for the interference effects, etc.^{46–48} Many approaches exist in the literature for overcoming this decoherence problem. Some attempt to add decoherence into Tully's surface hopping procedure,^{47,49–61} while others involve rigorous treatments at the density matrix level, e.g., quantum classical Liouville equation (QCLE)^{62–68} and the MeyerMillerStockThoss^{36,69,70} formalism. While these approaches all improve upon the FSSH algorithm, they are often either too costly, lead to complicated algorithms, or may not be accurate in certain scenarios.
In this paper, we describe a method based on Monte Carlo sampling of the semiclassical timedependent wavefunction of a molecule. Unlike common surfaces hopping approaches, the SCMC is not ad hoc by its construction, and is able to account for the quantum coherence effects between different photoinduced pathways. The main idea is based on writing the full molecular wavefunction as a series of terms that can be viewed as contributions of “partial” wavefunctions of the system, each corresponding to a certain number of transitions (hops) between electronic PES. These partial wavefunctions can be sampled by a stochastic process similar to that used in the surface hopping method: “onthefly” propagation of nuclei according to the classical equations of motion without prior knowledge of the PESs involved. The accumulated phase information characterizing each classical trajectory (e.g., the action along the trajectory) is further used to perform “postprocessing” evaluation of multidimensional integrals that correspond to each partial wavefunction using the Monte Carlo technique.
The method was first reported in Ref. 71. It was shown that the results obtained with our method are in very good agreement with numerically exact results for three standard sample problems (the socalled Tully problems).^{38} Additionally, comparison of the results for 1D and 2D test problems,^{71} demonstrated that the convergence of the method does not deteriorate with an increase in the number of nuclear degrees of freedom. This is a consequence of a fact that the stochasticity of the Monte Carlo procedure is related purely to the nondeterminism in the subspace of the electronic states (between the hops nuclei propagate according to the classical equations of motion, i.e., deterministically). In this paper, we present a detailed algorithmic description of the SCMC procedure, which should make its implementation straightforward. Additionally, we show that the Monte Carlo procedure used in our method is reasonably well convergent for the problems involving relatively few level crossings. Finally, we demonstrate and discuss how the convergence of the wavefunction is controlled by complex valued phase factors, the size of the nonadiabatic coupling region, and the choice of sampling function.
The QCLE procedure developed by Kapral et al.^{62–65} is similar in spirit to our method. The largest difference in the approach presented here is that, while Kapral's approach involves mixed classicalquantum dynamics of the density matrix using the WignerLiouville equation, our approach involves direct calculation of the dynamics of the wavefunctions. We believe our approach provides a more “straightforward” path to calculating accurate nonadiabatic molecular dynamics, with minimal deviation from the commonly used FSSH procedure. For example, such as FSSH, the dynamics in our approach always follow a single PES, whereas dynamics used to sample QCLE involve following single and averaged PES.^{62} Additionally, QCLE dynamics may require branched trajectories.^{63}
The paper is organized as follows: In Sec. II, we will formulate the problem and describe the theory behind the numerical method. In Sec. III, we describe, in detail, how to carry out the SCMC simulation, and provide an algorithm for a system with two PES. In Sec. IV, we analyze the results of the SCMC method when applied to Tully's 1D test problems and an additional three electronic levels, 1D, test problem.
II. THEORY
Consider a Hamiltonian with nuclear and electronic degrees of freedom
where
Here
In the absence of the nonadiabatic coupling, the PESs are decoupled and the nuclei propagate on a potential landscape given by a particular E_{i}(R). Formally, in order to describe quantum dynamics for a system of N particles corresponding to Hamiltonian
That is, the gaussian wavepacket with average momentum P_{c}(t) is centered around positions R_{c}(t), is evaluated according to the classical equations of motion, with S_{i}(t) being the classical action of the system produced during a time, t. For simplicity in this paper, we assume that the dispersion of the gaussian wavepacket corresponds to that of a free particle, σ^{2}(t) = σ^{2}(0) + iℏt/M. While this approximation is clearly violated for PESs with sufficient curvature, for dynamics in realistic molecules this is practically never the case. Indeed, a typical spread of the nuclear wavepacket is of the order of nuclear zero point motion associated with the vibrational ground states of a molecule, which is at least an order of magnitude smaller than the size of a typical electronic bond (which sets the length scale of E_{i}(R)) due to large mass, M, of the nuclei. Furthermore, this approximation did not lead to significant deviation from exact solution in the scattering results considered below. Thus, instead of solving a N + 1 dimensional Schrodinger equation, a problem reduces to solving N ordinary second order differential equations.
Unfortunately, for
Equation (4) is an infinite series with second and higher order terms being time ordered integrals and containing second and higher orders of perturbation. Here we neglected the
Suppose that the state Ψ(0)⟩ is a gaussian wavepacket localized on the PES 1. Furthermore, let us assume that the electronic subspace is two dimensional, i.e., there are only two close PES, the remaining PESs are separated by large gaps so that transition probabilities are infinitesimally small due to rapid oscillations of the integrands in Eq. (4). In fact such situation is usually the case; typically only two PES cross in a relevant region of the phase space. Then
Let us consider the second term in the RHS of Eq. (4) which corresponds to a “piece” of the wavefunction that has hopped once at time t_{1} from PES 1 to PES 2. It is clear that the wavepacket right before the hopping event at t_{1} can be approximated as a gaussian. It is reasonable to assume that after the scattering event at t_{1} the wavepacket retains its gaussian shape, but acquires an amplitude c_{12} (proportional to d_{12} at R_{c}(t_{1})). Similar assumptions apply to the higher order integrals. That is, we approximate the wavefunction in Eq. (4) as
where the normalization constant
Equation (5) is incomplete on itself: One needs to specify the c_{12/21}[R_{c}(t)] as well as the value of the momentum right after the hop. (Here we assume that neither the position nor the dispersion of the wavepacket changes immediately after the hop, which is easy to check by a simple calculation similar to the one discussed below.) For that let us divide the time integrals in Eqs. (4) and (5) into integrals over smaller time intervals Δt, where Δt is smaller than the timescale at which the particle moves over the distance at which potential
where ϕ_{2}(R, t) is a normalized wavepacket at the second PES. The t^{2} term in Eq. (6) arises due to two hop contribution in Eq. (5). Note that the wavefunction in Eq. (6) is normalized up to O(t^{2}). Since c_{12/21} and
Equations (7) and (8) provide the matching (boundary) conditions that determine values of c_{12/21} and P_{2}. Note that in the RHS of Eq. (8) we have used the “true” momentum of the system, which is given by its velocity
and
with P_{1} being the momentum right before the hop.
Similarly, the LHS of Eq. (8) gives P_{1} + c_{21}c_{12}t^{2}(P_{2} − P_{1}), where we have assumed that t is small enough and so the momenta P_{1} and P_{2} do not change during the evolution over the time interval t. The RHS yields P_{1} − d_{12}(R)(E_{2} − E_{1})(d_{12}(R) · P_{1})t^{2}/M, so that
Equation (10) gives a condition for the rescaling of momentum of the wavepacket which, at first glance, may seem to be in contradiction with that proposed by Herman^{76} and used by Tully^{38} in his surface hopping algorithm. Indeed, while in Eq. (10) the change in the momentum at the hop is along d_{12}, Eq. (10) does not conserve the total energy at the hop. However, careful inspection reveals that Eq. (10) is simply the zeroth and first order terms in an expansion, in powers of
III. SCMC CALCULATION
The SCMC provides a first principles method for calculating nonadiabatic molecular dynamics, with similar cost to surface hopping methods. The calculation occurs in two main steps: (A) the nonadiabatic propagation of classical trajectories, and (B) the use of trajectory data to sample the integrals in Eq. (5), with c_{ij}[R_{c}(t)] = d_{ij}(t). The classical propagation is similar to standard surface hopping algorithms.^{29,38,76} For simplicity, we assume that we have only two electronic states to consider. Generalization to higher numbers of states is straightforward.^{71} In Subsections III A and III B, we describe, in algorithmic fashion, the surfacehopping dynamics, A, and the Monte Carlo calculation of the timedependent semiclassical wavefunction, B. Figure 1 shows a diagrammatic map of the algorithm.
A. Semiclassical molecular dynamics

All trajectories are initialized with the same position and momentum. Various tracked values are initialized (see Fig. 1). We begin propagation of the nuclear wavepackets.

Begin time step. Supplied with a molecular geometry, R_{c}, quantum chemistry methods (e.g., DFT and TDDFT) are used to calculate the electronic state energies, E_{i}(t), the HellmannFeynman forces (electronic gradients), F_{c}, the first order nonadiabatic coupling vectors (NACV),
${\bf d}_{ij} ({\bf R}_c) = \langle \Psi _i ({\bf R}_c) \vert \frac{\partial }{\partial {\bf R}_c} \Psi _j ({\bf R}_c) \rangle$$dij(Rc)=\u27e8\Psi i(Rc)\u2202\u2202Rc\Psi j(Rc)\u27e9$, and scalars$d_{ij}(t)\break = \langle \Psi _i ({\bf R}_c(t)) \vert \frac{\partial }{\partial t} \Psi _j ({\bf R}_c(t)) \rangle = {\bf d}_{ij}({\bf R}_c) \cdot \dot{{\bf R}}_c(t)$$dij(t)=\u27e8\Psi i(Rc(t))\u2202\u2202t\Psi j(Rc(t))\u27e9=dij(Rc)\xb7R\u0307c(t)$.  The classical particles are propagate for one time step, ∂t, on the current adiabatic surface, i, using Newton/Lagrange equations of motion,by some numerical method, e.g., Verlet, RungaKutta. The classical action along the trajectory,(11)\begin{align} \frac{\partial }{\partial t} {\bf P}_c(t)={\bf F}_c(t)= \frac{\partial }{\partial {\bf R}_c} E_i( {\bf R}_c(t)) , \end{align}$\u2202\u2202tPc(t)=Fc(t)=\u2212\u2202\u2202RcEi(Rc(t)),$is calculated for subsequent use in part B. We assume our propagated gaussian nuclear wavepackets are “free,” i.e., they broaden with time, independent of the PES.(12)\begin{align} S(t)&= \int ^{t}_0 dt^{\prime } \, \sum _{\alpha }^{N} \frac{{P}_\alpha ^2(t^{\prime })}{2M_\alpha }  E(t^{\prime }) \nonumber\\ &= \sum _{\alpha }^{N}\int ^{{R}_c^\alpha (t)}_{{R}^\alpha _{c\,0}} d\bar{{R}}_c^\alpha \,{P}_\alpha (\bar{{R}}_c^\alpha ) , \end{align}$S(t)=\u222b0tdt\u2032\u2211\alpha NP\alpha 2(t\u2032)2M\alpha \u2212E(t\u2032)=\u2211\alpha N\u222bRc0\alpha Rc\alpha (t)dR\xafc\alpha P\alpha (R\xafc\alpha ),$
 At the end of the time step, we determine whether or not to hop to a new state, j. First, we determine whether the hop is allowed or frustrated by, the assumed, conservation of energy, i.e.,If the hop is allowed the probability, γ_{i→j}(t), is calculated. If not γ_{i→j}(t) = 0. Unlike other surface hopping methods,^{29,38,76} the final result is formally independent of the choice of hopping rate, assuming that the hopping rate is nonzero everywhere d_{ij}(t) is nonzero. However, the choice of hopping rate is crucial to achieve rapid convergence of the result.(13)\begin{align} \frac{{\bf P}_c^2(t^{\prime })}{2M} + E_i({\bf R}_c) \ge E_j({\bf R}_c) . \end{align}$Pc2(t\u2032)2M+Ei(Rc)\u2265Ej(Rc).$
 At each time step, the hopping probability to the new state is added to the integralHere, m is the number of hops over the course of the trajectory, and t_{l} is the time at which the lth hop occurs (t_{m + 1} ≡ t, s_{m + 1} ≡ j). If RN < γ_{i→j}, where RN is a generated random number between 0 and 1, then the trajectory hops to state j.(14)\begin{align} \Omega ({\bf t})=\sum _{l=0}^{l=m} \int _{t_l}^{t_{l+1}} dt^{\prime } \, \gamma _{s_l\rightarrow s_{l+1}}[{\bf R}_c(t^{\prime })] . \end{align}$\Omega (t)=\u2211l=0l=m\u222btltl+1dt\u2032\gamma sl\u2192sl+1[Rc(t\u2032)].$
 If a hop occurs, the productsand(15)\begin{align} &D(t_{m1}\ldots t_0)=\prod _{l=0}^{m1} d_{s_{l}\, s_{l+1}}(t_l), \end{align}$D(tm\u22121...t0)=\u220fl=0m\u22121dslsl+1(tl),$are multiplied by d_{ij}(t) and γ_{i→j}(t), respectively. The number of hops, m, is increased by one.(16)\begin{align} &\Gamma (t_{m1}\ldots t_0) =\prod _{l=0}^{m1} \gamma _{s_l\rightarrow s_{l+1}}(t_l) \end{align}$\Gamma (tm\u22121...t0)=\u220fl=0m\u22121\gamma sl\u2192sl+1(tl)$
The procedure is repeated for the next time step on either the original or new state. These dynamics continue for the desired time. Once completed, the final position, R_{c}(t), and momentum, P_{c}(t), of the trajectory are stored. The process is repeated for the specified number of trajectories. In calculations of realistic molecular systems,
B. Postprocessing
Information collected by the surface hopping dynamics is further analyzed in the second step at negligible computational expense when compared to part A. Namely, we calculate the nuclear wavefunctions,
or probabilities, P_{i} = ∫dR Ψ_{i}(R, t)^{2}. Here t ≡ {t_{1}, …t_{m}, t} is the set of time integration variables and
Using the standard Monte Carlo approach, the integral
Here, k is the index of a trajectory from step 1,
The kth value of the probability distribution can be expressed as
This probability can also be expressed as an expectation value, using the same probability distribution,
Ω_{i}(t_{i}) is given by Eq. (14). We seek the value of
and
We first determine
$N_i^m$$Nim$ and calculate the sum in Eq. (21) for every i/m pair.We calculate each
$A_i^m$$Aim$. For a given trajectory${\bf I}_{i}^m ({\bf R}_k,{\bf t}_k)$$Iim(Rk,tk)$ is easily calculated, see Eq. (5), from the surface hopping data, Eqs. (11), (12), and (15),where ϕ is a “free” gaussian wavepacket,(22)\begin{eqnarray} {\bf I}_{i}^m ({\bf R}_k,{\bf t}_k) &=& D_{k}({\bf t}_k) \times e^{i [ S^i ({\bf t}_k) {\bf P}^k (t) {\bf R}_c^k (t) ] } \nonumber\\ & & \times\, \bm{\phi } \big({\bf R},{\bf R}_c^k,{\bf P}_c^k,t\big) , \end{eqnarray}$Iim(Rk,tk)=Dk(tk)\xd7ei[Si(tk)\u2212Pk(t)Rck(t)]\xd7\varphi R,Rck,Pck,t,$with the final position and momentum of the trajectory k, and an initial width, σ. At this point the nuclear wave functions can be calculated by summing over all the trajectories, Eq. (18).(23)\begin{align} \bm{\phi } \big({\bf R},{\bf R}_c^k,&{\bf P}_c^k,t\big)= \root 4 \of {\frac{M^2 \sigma ^2 }{\pi [ M^2\sigma ^4 + t^2]}} e^{i{\bf P}^k(t){\bf R}}e^{\frac{[{\bf R}{\bf R}_c^k (t)]^2}{\sigma ^2 + \frac{it}{M}}} , \end{align}$\varphi (R,Rck,Pck,t)=M2\sigma 2\pi [M2\sigma 4+t2]4eiPk(t)Re[R\u2212Rck(t)]2\sigma 2+itM,$
Once Ψ_{i}(R, t) is known, expectation values (e.g., x, P, density matrix) can be calculated. When the number of trajectories is low, the random error in Ψ_{i}(R, t) can be reduced by normalization. However, the effect is minimal in a well converged calculation.
Some notes on the case of more than 2 electronic states: (1) The rate in Eq. (14) would be replaced with a sum over the rates to all possible new states. (2) Eqs. (15) and (16) remain unchanged. (3) Determination of which state to hop to can be determined as it is in Ref. 38. (4) For every possible path of intermediate electronic excited states, there is an integral which must be sampled independently. Thus each distinct path has a different probability distribution, Eq. (21).
IV. ANALYSIS OF SCMC METHOD
In order to test the SCMC procedure, we consider three twolevel models proposed by Tully.^{38} Figure 2 shows the energy eigenstates and NACVs of these three onedimensional problems. For the exact form of the Hamiltonians see Ref. 38. For these three problems, scattering probabilities calculated using the SCMC, FSSH, and Ehrenfest methods are compared to the exact timedependent Schrödinger equation in Gorshkov et al.^{71} In Figure 3, we represent those scattering results. For Problems 1 and 3: 25 000 and 10 000 trajectories were used for the SCMC and FSSH approaches, respectively. For Problem 2: 75 000 trajectories were used in the SCMC calculation, while only 10 000 were used in the FSSH calculation. Here we seek to provide indepth analysis of these results and insight into the method. All parameters and results in this section are in atomic units.
Problem 3 provides a simple model problem where the standard implementation of the FSSH quantitatively and qualitatively fails to reproduce the exact TDSchrödinger result, see Figs. 3(c) and 3(d). It is well known that this failure is due to the lack of decoherence.^{38,77} The equation of motion for the electronic density matrix, or complex expansion coefficients, does not include the effect of bifurcation of nuclear wavepackets. Corrections can be introduced by hand into the FSSH algorithm to correct for this problem.^{47,49–61} In the SCMC approach, we do not need to propagate the electronic density matrix in order to determine hopping rates. We find the simple hopping rate γ_{1 → 2}(t) = d_{12}(t) leads to accurate results for all of the test problems, see Figs. 3 and 4. Thus, since we do not utilize the electronic density matrix, we do not encounter the “over coherence” problem. The left (right) column of Figure 4 shows the results of the SCMC (timedependent Schrödinger) calculation for Problem 3 with an initial wavefunction
with initial position: x_{0} = −12, momentum: k = 10, and width:
To test the limits of our method, we consider the three level crossing “Model X” problem from Ref. 48 (see Ref. 48 for Hamiltonian). The PES can be seen in Fig. 2(d). We consider a particle entering from negative x on the middle surface. The particle has a welldefined momentum (plane wave for the exact solution,
While the SCMC method does a good job of reproducing the exact result, it does so at the cost of a high number of trajectories. For a problem like “Model X” significant improvement in the convergence could be achieved by using the SCMC only in the vicinity of the strong NAC, then adiabatically propagating wavepackets with the normalization, position, momentum, and phase given by SCMC. This is similar, in essence, to the branching scheme recently proposed by Herman, which significantly improved convergence in a two crossing problem.^{80}
We have shown how the details of the problem, number, and nature of level crossings, can affect the number of trajectories required, i.e., 2.5 × 10^{4} for a simple avoided crossing (∼1%−2% absolute error) with 3 × as many required for the double avoided crossing. In Subsection IV A, we will discuss the factors which lead to this increased cost in trajectories.
A. Role of phases
While the SCMC method does a good job of reproducing the exact result, including when FSSH and Ehrenfest methods fail, it does so at an increased cost in trajectories. In order to determine the areas where SCMC is most applicable, we seek to understand how the convergence rate is affected by different models and parameters. First, we look at the simple single level crossing, Problem 1. Figure 6 (Inset) shows the relative error, δ_{rel}, in the calculated probabilities of transmission on the upper (lower) surfaces, TU (TL), as a function of the total number of trajectories, N. For all momentums and both probabilities, the convergence obeys an approximate power law of the form
This k dependence can be understood by looking at the phase in Eq. (22). Assuming a onedimensional system with two electronic states, we can calculate the phase for a trajectory which begins on the lower surface and hops once,
Since all trajectories have the same final and initial times, we can subtract off any term which does not depend on the time or position of the hop. Additionally, if we assume that the initial momentum, k, is large enough that we do not have any reflection, and t_{f} is large enough that P_{2}(t_{f}) is a constant, we can rewrite Eq. (24) as a relative phase that depends on the position of the hop,
With the same considerations we can write the relative phase for a trajectory with two hops,
In Figure 7, we can see that the number of oscillations, which must be sampled, is reduced in the one (a) and (b) and two (c)–(f) hop integrals when we go from the low momentum, k = 10, to the high momentum k = 30. Indeed for integrals of all numbers of hops the number of oscillations will decrease with an increasing momentum. The “smoother” integrands associated with higher momentum are easier to sample via Monte Carlo methods. For the TL in Figure 6, the effect competes with constant error sources, due to, e.g., the free gaussian approximation, which appears as an increasing relative error with decreasing average probability (probability of TL decreases with increase k). This leads to the nonmonotonic k dependence.
When going from a single (Problem 1) to double (Problem 2) crossing systems the integrated area of d_{12}(x) increases from 1.57 to 2.60. The number of trajectories required to sample an integral with m hops will approximately scale with ∫dx d_{12}(x) as
However, integrals with more hops, generally, will contribute less to the total wavefunction. Additionally, oscillations due to the phases, Eqs. (25) and (26), increase for Problem 2 (see Fig. 8). These two, not independent, effects result in a much slower convergence rate for Problem 2 as compared to Problem 1 or 3. It was shown in Ref. 71 that ∼3 times the number of trajectories is required to obtain similar precision for Problem 2 as for Problems 1 and 3.
B. Choice of hopping rate
In principle, the hopping rate, γ_{i→j}, is arbitrary. It is simply the choice of a probability distribution for sampling the integral, Eq. (5). A better probability distribution leads to less samples being required to accurately calculate the correct result. The closer the probability distribution is to the absolute value of the normalized integrand, the better the sampling. If we knew the normalization factor for each integrand
One choice which may come to the readers mind, is the fewest switches hopping rate from Tully,^{38}
with
Here a is the electronic density matrix. Figure 9 shows the relative error of the SCMC calculation for the different hopping rates,
Additionally, when the integrals are unevenly sampled without consideration of the phases, as they are when using
V. CONCLUSION
The semiclassical Monte Carlo method provides an approach, which is not ad hoc, for calculating nonadiabatic dynamics in extended molecular systems. It is based on an expansion of the full molecular wavefunction into partial wavefunctions corresponding to the number of transitions between electronic states. The procedure involves standard surface hopping dynamics with an additional postprocessing step (calculation of the wave functions from collected data). While in the surface hopping method the hopping probability controls all of the dynamics; in SCMC it is primarily interference between trajectories with different (or similar) phases that determines the final result. For this reason, the SCMC method does not suffer from the same “over coherence” problem as the FSSH approach.
Some notes on the applicability of this method: (1) While the phases lead to difficulty in sampling for small momentum, for higher momentums the sampling becomes easier. (2) Additionally, we believe that convergence rate does not increase with increased number of system coordinates.^{71} (3) Increasing the number of crossing regions increases the number of trajectories required for sampling, see Eq. (27), due to an increase in the area of nonadiabatic coupling and an increase in the number of oscillations due to the phases. (4) While it has been shown that, in some systems, FSSH approach requires complete coupling, i.e., the quantum mechanical electronic dynamics requires calculation of NACV for all combinations of states (
ACKNOWLEDGMENTS
We acknowledge the support of the U.S. Department of Energy through the Los Alamos National Laboratory (LANL) LDRD Program. LANL is operated by Los Alamos National Security, LLC, for the National Nuclear Security Administration of the U.S. Department of Energy under Contract No. DEAC5206NA25396. We acknowledge the support of the Center for Nonlinear Studies (CNLS) and the Center for Integrated Nanotechnology (CINT) at LANL. We also thank J.E. Subotnik for sharing his numerical results reported in Ref. 48.