Surface hopping is a mixed quantum-classical technique that incorporates quantum mechanical effects into molecular dynamics simulations.[1] Traditional molecular dynamics assume the Born-Oppenheimer approximation, where the lighter electrons adjust instantaneously to the motion of the nuclei. Though the Born-Oppenheimer approximation is applicable to a wide range of problems, there are several applications, such as photoexcited dynamics, electron transfer, and surface chemistry where this approximation falls apart. Surface hopping partially incorporates the non-adiabatic effects by including excited adiabatic surfaces in the calculations, and allowing for 'hops' between these surfaces, subject to certain criteria.
Molecular dynamics simulations numerically solve the classical equations of motion. These simulations, though, assume that the forces on the electrons are derived solely by the ground adiabatic surface. Solving the time-dependent Schrödinger equation numerically incorporates all these effects, but is computationally unfeasible when the system has many degrees of freedom. To tackle this issue, one approach is the mean field or Ehrenfest method, where the molecular dynamics is run on the average potential energy surface given by a linear combination of the adiabatic states. This was applied successfully for some applications, but has some important limitations. When the difference between the adiabatic states is large, then the dynamics must be primarily driven by only one surface, and not an average potential. In addition, this method also violates the principle of microscopic reversibility.
Surface hopping accounts for these limitations by propagating an ensemble of trajectories, each one of them on a single adiabatic surface at any given time. The trajectories are allowed to 'hop' between various adiabatic states at certain times such that the quantum amplitudes for the adiabatic states follow the time dependent Schrödinger equation. The probability of these hops are dependent on the coupling between the states, and is generally significant only in the regions where the difference between adiabatic energies is small.
The formulation described here is in the adiabatic representation for simplicity. It can easily be generalized to a different representation.The coordinates of the system are divided into two categories: quantum (
q
R
mn
H=\sumn-
\hbar2 | |
2mn |
2 | |
\nabla | |
qn |
+V(q,R)
V
H
R
\phin(q;R)
q
\begin{align} FR&=-\nablaR\langle\phin|H|\phin\rangle\\ &=-\langle\phin|\nablaRH|\phin\rangle, \end{align}
n
R
\psi(q;R,t)=\sumncn(t)\phin(q;R)
cn(t)
i\hbarcj |
=\sumncn\left(Vjn-i\hbar
R |
.djn\right)
Vjn
djn
\begin{align} Vjn&=\langle\phij|H|\phin\rangle=\langle\phij|H|\phij\rangle\deltajn\\ djn&=\langle\phij|\nablaR\phin\rangle \end{align}
2 | |
|c | |
j(t)| |
2 | |
|c | |
j(t)| |
|
=\sumn
2 | |
\hbar |
Im(anjVjn)-2Re(anj
R |
.djn)
anj=cnc
* | |
j |
2 | |
|c | |
j(t)| |
| |||||||
|
≈
dt | |
ajj |
\sumn
2 | |
\hbar |
Im(anjVjn)-2Re(anj
R |
.djn)
j
Pj\to=
dt | \left( | |
ajj |
2 | |
\hbar |
Im(anjVjn)-2Re(anj
R |
.djn)\right)
Whenever a hop takes place, the velocity is adjusted to maintain conservation of energy. To compute the direction of the change in velocity, the nuclear forces in the transition is
\begin{align} \langle\phij|\nablaRH|\phin\rangle&=\nablaR\langle\phij|H|\phin\rangle-\langle\nablaR\phij|H|\phin\rangle-\langle\phij|H|\nablaR\phin\rangle\\ &=\nablaREj\deltajn+(Ej-En)djn, \end{align}
Ej=\langle\phij|H|\phij\rangle
djn=-dnj
djn
djn
\deltat=\hbar/2\DeltaE
\DeltaE
Surface hopping can develop nonphysical coherences between the quantum coefficients over large time which can degrade the quality of the calculations, at times leading the incorrect scaling for Marcus theory. To eliminate these errors, the quantum coefficients for the inactive state can be damped or set to zero after a predefined time has elapsed after the trajectory crosses the region where hopping has high probabilities.
The state of the system at any time
t
Step 1. Initialize the state of the system. The classical positions and velocities are chosen based on the ensemble required.
Step 2. Compute forces using Hellmann-Feynman theorem, and integrate the equations of motion by time step
\Deltat
t+\Deltat
Step 3. Integrate the Schrödinger equation to evolve quantum amplitudes from time
t
t+\Deltat
\deltat
\deltat
\Deltat
Step 4. Compute probability of hopping from current state to all other states. Generate a random number, and determine whether a switch should take place. If a switch does occur, change velocities to conserve energy. Go back to step 2, till trajectories have been evolved for the desired time.
The method has been applied successfully to understand dynamics of systems that include tunneling, conical intersections and electronic excitation.
In practice, surface hopping is computationally feasible only for a limited number of quantum degrees of freedom. In addition, the trajectories must have enough energy to be able to reach the regions where probability of hopping is large.
Most of the formal critique of the surface hopping method comes from the unnatural separation of classical and quantum degrees of freedom. Recent work has shown, however, that the surface hopping algorithm can be partially justifiedby comparison with the Quantum Classical Liouville Equation. It has further been demonstrated that spectroscopic observables can be calculated in close agreement with the formally exact hierarchical equations of motion.