Scattering of a mono-energetic particle from a surface

### Loading required packages

using NQCDynamics, NQCCalculators
using Unitful, UnitfulAtomic
using Statistics
using GLMakie

This tutorial covers how to run atomic scattering simulations from a metal surface with NQCDynamics where electronic excitations in various flavours are included. To demonstrate how this is done, we will use a simple one-dimensional model, a Newns-Anderson Hamiltonian constructed from a simple one-dimensional, diabatic two-state model, for a neutral and a charged particle interacting with a substrate.

With intention to convey the concepts in the most simplifying manner, we will focus onto a monoenergetic ensemble of impinging Li atoms. At the very beginning, we will define the initial conditions and the other propagation settings which is done in the cell below.

This toy model is supposed to mimic the scattering of a single, mono-energetic Li atom from a metal surface. #7up!

### Define initial conditions and output

atoms = Atoms(7u"u") # Our projectile is a Li atom; #7up!
r0 = austrip(5u"Å") # Projectile is 5A away from surface
m = atoms.masses[1] # Mass of projectile
ke = austrip(2u"eV") # Initial kinetic energy
v0 = -sqrt(2*ke/m) # Velocity direction towards the surface with magnitude defined by kinetic energy
dt = 0.01u"fs"

0.01 fs

Termination criteria & Output

We are interested in scattering simulations. Typcially for this type of simulations, we need to ensure that the run time is sufficiently long which we ensure by setting the maximum distance the particle is allowed to travel to three times of the initial distance from the surface. This is very sufficient for a simple 1D model (in full dimensions, where subsurface scattering can occur, we might need to be more tolerant to longer trajectories). With velocity calculated from the initial kinetic energy, this gives us our maximum run time parameter, tcut.

The parameter ‘tcut’ defines the maximum runtime. Yet, this single termination criterion is not enough because the scattered particle will travel through the vacuum until the run time meets our termination criterion (in periodic surface cells this is even worse because the particle might impact with the periodic images of the surface in perpendicular direction). Lest to artificially increase the required time for our batch of trajectories, we insert a second criterion namely that the trajectory gets terminated once the particle’s position is higher than it’s initial position and the sign of the velocity has changed.

Next, we also set the number of trajectories we want to run: the amount of required trajectories are typically determined by the quantity of interest and sensitivity of the quantity of interest against the number of trajectories should always be checked! Here, we will use only 10 trajectories to ensure very quick runtimes so that one gets a feeling for the simulation and can easily play with other parameters (:

Finally, we also set the required output. For scattering simulations, all quantities of interest can typically extracted from the initial & final positions along with the inital & final velocities of the system. For convenience, we also set the kinetic energy directly albeit not being necessary from a technical point of view.

dcut = 3*r0
tcut = abs(dcut/v0)

function termination_condition(u, t, integrator)::Bool

    return ((t > tcut) || ((mean(DynamicsUtils.get_positions(u)) > austrip(5.5u"Å")) && (mean(DynamicsUtils.get_velocities(u)) > 0)))
end

ntrajs = 10
output= (OutputVelocity, OutputPosition, OutputKineticEnergy, OutputDynamicsVariables)

### Defining termination criterion

(OutputVelocity, OutputPosition, OutputKineticEnergy, OutputDynamicsVariables)

Definition of the model

Now it’s getting more technical. The first two lines of the cell below define the parameterisation of the diabatic two-state model. Γ is a measure for the coupling between the neutral state and the charged state of the incoming atom (the larger Γ, the stronger the avoided crossing of the adiabats and thus the less non-adiabatic the simulations are going to be

The next lines describe the discretisation of the electronic band from the spectral properties of the bath to the method of the discretisation. Here, we’re using the ‘ShenviGaussLegendre’ method, the method of choice in NQCDynamics for metal substrates.

Next, we construct our Newns-Anderson Hamiltonian from the two diabates, the finite number of bathstates and the couplings between them which is done with the keyword ‘AndersonHolstein’. The following two keywords define the temperature of the electronic subsystem of the solid, i.e., the population of the substrate’s states with the electrons.

#### Using model and setting its parameters

Γ = austrip(0.2u"eV") # Defines coupling strength between diabatic ground state and diabatic excited state
thossmodel = ErpenbeckThoss(;Γ, m=atoms.masses[1]) # This model defines the position dependence of the ingredients for the Newns-Anderson Hamiltonian U0, U1, and Vk. Here, we're using the Erpenbeck-Thoss model.

fermi_level = 0*u"eV" # We reference the spectrum of electronic states of the surface to the Fermi level
nstates = 20 # Number of states into which the electronic continuum of the substrate is discretised into
bandwidth = austrip(50*u"eV") #Spectral width of the electronic bands
bandmin = -bandwidth / 2
bandmax = bandwidth / 2
bath = ShenviGaussLegendre(nstates, bandmin, bandmax) #Discretisation method

model = AndersonHolstein(thossmodel, bath)
temperature = 300u"K"
β = 1/austrip(temperature)

1052.5834160174613

Combining the intial conditions for both the electronic and nuclear subspaces

Since we’re evolving both the electrons and the nuclei in mixed-quantum classical dynamics, we need to bring together both the electronic space, and the nuclear space. The first keyword in the cell below, defines the dimensionality of the nuclear phasespace, and the second line defines the dimensionality of the electronic state with the third line folding both subspaces into one big space. The final line sets the termination criterion.

### Setting up distributions

nuclear_distribution = DynamicalDistribution(v0, r0, (1,1))
electronic_distribution = FermiDiracState(fermi_level, temperature)
dist = electronic_distribution * nuclear_distribution

terminate = DynamicsUtils.TerminatingCallback(termination_condition)

SciMLBase.DiscreteCallback{typeof(termination_condition), typeof(SciMLBase.terminate!), typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(SciMLBase.FINALIZE_DEFAULT), Nothing, Tuple{}}(termination_condition, SciMLBase.terminate!, SciMLBase.INITIALIZE_DEFAULT, SciMLBase.FINALIZE_DEFAULT, Bool[1, 1], nothing, (), true)

Ehrenfest

We have now set everything up to successfully start mixed-quantum classical scattering trajectories. We first start with Ehrenfest simulations. Ehrenfest dynamics is governed by a potential energy surfaces which is averaged from the entire spectrum of potential energy surfaces which are obtained from the Newns-Anderson Hamiltonian. Since Ehrenfest dynamics is entirely determinestic, we only need a single trajectory given our deterministic initial conditions (if we had an ensemble of different positions & velocities, things would be different).

The first line defines the simulation environment. The second line launches our IESH simulations. :)

### Run single Ehrenfest calculation

sim = Simulation{EhrenfestNA}(atoms, model)

result = run_dynamics(sim,
                    (0.0, tcut),
                    dist;
                    trajectories=1,
                    callback=terminate,
                    output=output,
                    saveat=[0,tcut],
                    dt
                   )

Independent electron surface hopping

Next, we will run independent electron surface hopping (IESH) simulations. Since IESH is stochastic, we need to run a batch of trajectories. Here, we do only 10 trajectories for briefness’s sake but in practice one should do several hundreds for a 1D model to ensure numerical convergence!

### Run IESH simulations
sim = Simulation{AdiabaticIESH}(atoms, model)

result = run_dynamics(
    sim,
    (0.0, tcut),
    dist;
    trajectories=ntrajs,
    callback=terminate,
    output=(output..., OutputDiscreteState),
    # saveat=[0,tcut],
    dt = dt,
    reduction = SortByOutputReduction(),
)

[ Info: Sampling randomly from provided distribution.
[ Info: Pre-compiled dynamics in 6.153458666 seconds.
[ Info: Performing 10 trajectories.
[ Info: Sampling randomly from provided distribution.
[ Info: Finished after 19.464154958 seconds.

6-element Dictionaries.Dictionary{Symbol, Vector}:
                    :Time │ [[0.0, 0.4134137333533518, 0.8268274667067036, 1.24…
          :OutputVelocity │ [[[-0.0033941074835368125;;], [-0.00339411452744389…
          :OutputPosition │ [[[9.448630623128851;;], [9.447227451027294;;], [9.…
     :OutputKineticEnergy │ [[0.07349864435103717, 0.07349894941984327, 0.07349…
 :OutputDynamicsVariables │ Vector{NQCDynamics.DynamicsMethods.SurfaceHoppingMe…
     :OutputDiscreteState │ [[1, 1, 1, 1, 1, 1, 1, 1, 1, 1  …  1, 1, 1, 1, 1, 1…

function plot_trajectory(idx)
    colors = [colorant"#201e50", colorant"#ed254e"]
    GLMakie.activate!(; px_per_unit=4.0,)
    fig = Figure(
        size = 28.1 .* (18,12),
        fontsize = 20,
    )
    potential_ax = Axis(
        fig[1,1],
        xgridvisible = false,
        xminorgridvisible = false,
        ygridvisible = false,
        yminorgridvisible = false,
        xlabel = "Distance to surface / a.u.",
        ylabel = "Energy / Hartree",
        xautolimitmargin = (0.0,0.0),
    )
    r_values = 0:0.02:austrip(5.0u"Å") |> collect
    model_potential = [NQCDynamics.NQCModels.potential(thossmodel, hcat(R)) for R in r_values]
    lines!(
        potential_ax,
        r_values,
        getindex.(model_potential, 1, 1),
        label = "ground state",
        color = colors[1],
        linewidth = 2.0,
    )
    lines!(
        potential_ax,
        r_values,
        getindex.(model_potential, 2, 2),
        label = "excited state",
        color = colors[2],
        linewidth = 2.0,
    )
    particle_pos = scatter!(
        potential_ax,
        result[:OutputPosition][idx][1500] |> first,
        NQCModels.potential(thossmodel, result[:OutputPosition][idx][1500])[
            result[:OutputDiscreteState][idx][1500] |> last,
            result[:OutputDiscreteState][idx][1500] |> last,
        ],
        markersize = 24,
        strokewidth = 0.75,
        strokecolor = colorant"black",
        # marker = :circle,
    )
    xlims!(0.5, 8)
    ylims!(-0.2,0.75)
    # Electrons in IESH states
    electrons_ax = Axis(
        fig[1,2],
        xgridvisible = false,
        xminorgridvisible = false,
        ygridvisible = false,
        yminorgridvisible = false,
        width = 50,
        xticklabelsvisible = false,
        yticklabelsvisible = false,
        xticksvisible = false,
        yticksvisible = false,
    )
    colgap!(fig.layout, 1,4)
    hlines!(bath.bathstates;color = colorant"black", xmin = 0.25, xmax = 0.75)
    el_sc = scatter!(
        electrons_ax,
        [1 for i in 1:nstates/2],
        bath.bathstates[result[:OutputDynamicsVariables][idx][1500].state .|> Int],
        color = bath.bathstates[result[:OutputDynamicsVariables][idx][1500].state .|> Int],
        # strokewidth = 0.75,
        markersize = 20,
        marker = '↑',
        colormap = :turbo
    )
    ylims!(electrons_ax, bath.bathstates[7] - 0.025, bath.bathstates[13] + 0.025)
    axislegend(potential_ax, position = :rt)
    record(fig, "../resources/Tutorial_MQCD_methods/trj2.mp4", 1:15:length(result[:OutputPosition][idx])) do i
        particle_pos.positions[] = [
            Point2f(
                result[:OutputPosition][idx][i] |> first,
                NQCModels.potential(thossmodel, result[:OutputPosition][idx][i])[
                    result[:OutputDiscreteState][idx][i] |> last,
                    result[:OutputDiscreteState][idx][i] |> last,
                ],
            )]
        particle_pos.color = colors[result[:OutputDiscreteState][idx][i] |> first]
        el_sc.positions[] = [
            Point2f(1, bath.bathstates[Int(el)]) for el in result[:OutputDynamicsVariables][idx][i].state
        ]
    end
end
plot_trajectory(2)
nothing

Example IESH trajectory animation. The electronic ground and excited states for the particle scattering from the surface are shown in the left panel. The right panel shows impurity and bath electronic states and their current occupation with independent electrons.

Molecular Dynamics with Electronic Friction

### Run MDEF calculations
ρ = (nstates+1)/(bandwidth)
dist = nuclear_distribution
bath = TrapezoidalRule(nstates, bandmin, bandmax)
model = AndersonHolstein(thossmodel, bath)

sim = Simulation{DiabaticMDEF}(atoms, model; temperature,
    friction_method=DynamicsMethods.ClassicalMethods.WideBandExact(ρ, β)
)


result = run_dynamics(sim,
                    (0.0, tcut),
                    dist;
                    trajectories=ntrajs,
                    callback=terminate,
                    output=output,
                    saveat=[0,tcut],
                    dt
                   )

[ Info: Sampling randomly from provided distribution.
[ Info: Pre-compiled dynamics in 2.92208725 seconds.
[ Info: Performing 10 trajectories.
[ Info: Sampling randomly from provided distribution.
[ Info: Finished after 5.722769666 seconds.

10-element Vector{Dictionaries.Dictionary{Symbol, Any}}:
 {:Time = [0.0, 8351.500948888295], :OutputVelocity = [[-0.0033941074835368125;;], [-3.67352697973401e-5;;]], :OutputPosition = [[9.448630623128851;;], [3.733364305552911;;]], :OutputKineticEnergy = [0.07349864435103717, 8.609830740728939e-6], :OutputDynamicsVariables = RecursiveArrayTools.ArrayPartition{Float64, Tuple{Matrix{Float64}, Matrix{Float64}}}[RecursiveArrayTools.ArrayPartition{Float64, Tuple{Matrix{Float64}, Matrix{Float64}}}(([-0.0033941074835368125;;], [9.448630623128851;;])), RecursiveArrayTools.ArrayPartition{Float64, Tuple{Matrix{Float64}, Matrix{Float64}}}(([-3.67352697973401e-5;;], [3.733364305552911;;]))]}
 {:Time = [0.0, 7730.836813706072, 7730.836813706072], :OutputVelocity = [[-0.0033941074835368125;;], [0.003042455053074378;;], [0.003042455053074378;;]], :OutputPosition = [[9.448630623128851;;], [10.394725376115256;;], [10.394725376115256;;]], :OutputKineticEnergy = [0.07349864435103717, 0.05905769440279816, 0.05905769440279816], :OutputDynamicsVariables = RecursiveArrayTools.ArrayPartition{Float64, Tuple{Matrix{Float64}, Matrix{Float64}}}[RecursiveArrayTools.ArrayPartition{Float64, Tuple{Matrix{Float64}, Matrix{Float64}}}(([-0.0033941074835368125;;], [9.448630623128851;;])), RecursiveArrayTools.ArrayPartition{Float64, Tuple{Matrix{Float64}, Matrix{Float64}}}(([0.003042455053074378;;], [10.394725376115256;;])), RecursiveArrayTools.ArrayPartition{Float64, Tuple{Matrix{Float64}, Matrix{Float64}}}(([0.003042455053074378;;], [10.394725376115256;;]))]}
 {:Time = [0.0, 7178.5160659463545, 7178.5160659463545], :OutputVelocity = [[-0.0033941074835368125;;], [0.003197021146561491;;], [0.003197021146561491;;]], :OutputPosition = [[9.448630623128851;;], [10.394247291570297;;], [10.394247291570297;;]], :OutputKineticEnergy = [0.07349864435103717, 0.06521074532533636, 0.06521074532533636], :OutputDynamicsVariables = RecursiveArrayTools.ArrayPartition{Float64, Tuple{Matrix{Float64}, Matrix{Float64}}}[RecursiveArrayTools.ArrayPartition{Float64, Tuple{Matrix{Float64}, Matrix{Float64}}}(([-0.0033941074835368125;;], [9.448630623128851;;])), RecursiveArrayTools.ArrayPartition{Float64, Tuple{Matrix{Float64}, Matrix{Float64}}}(([0.003197021146561491;;], [10.394247291570297;;])), RecursiveArrayTools.ArrayPartition{Float64, Tuple{Matrix{Float64}, Matrix{Float64}}}(([0.003197021146561491;;], [10.394247291570297;;]))]}
 {:Time = [0.0, 4338.777131544037, 4338.777131544037], :OutputVelocity = [[-0.0033941074835368125;;], [0.0030350946811241748;;], [0.0030350946811241748;;]], :OutputPosition = [[9.448630623128851;;], [10.393850011843464;;], [10.393850011843464;;]], :OutputKineticEnergy = [0.07349864435103717, 0.05877229278586498, 0.05877229278586498], :OutputDynamicsVariables = RecursiveArrayTools.ArrayPartition{Float64, Tuple{Matrix{Float64}, Matrix{Float64}}}[RecursiveArrayTools.ArrayPartition{Float64, Tuple{Matrix{Float64}, Matrix{Float64}}}(([-0.0033941074835368125;;], [9.448630623128851;;])), RecursiveArrayTools.ArrayPartition{Float64, Tuple{Matrix{Float64}, Matrix{Float64}}}(([0.0030350946811241748;;], [10.393850011843464;;])), RecursiveArrayTools.ArrayPartition{Float64, Tuple{Matrix{Float64}, Matrix{Float64}}}(([0.0030350946811241748;;], [10.393850011843464;;]))]}
 {:Time = [0.0, 4140.338539534558, 4140.338539534558], :OutputVelocity = [[-0.0033941074835368125;;], [0.003101770562619906;;], [0.003101770562619906;;]], :OutputPosition = [[9.448630623128851;;], [10.393833158673782;;], [10.393833158673782;;]], :OutputKineticEnergy = [0.07349864435103717, 0.06138291181411692, 0.06138291181411692], :OutputDynamicsVariables = RecursiveArrayTools.ArrayPartition{Float64, Tuple{Matrix{Float64}, Matrix{Float64}}}[RecursiveArrayTools.ArrayPartition{Float64, Tuple{Matrix{Float64}, Matrix{Float64}}}(([-0.0033941074835368125;;], [9.448630623128851;;])), RecursiveArrayTools.ArrayPartition{Float64, Tuple{Matrix{Float64}, Matrix{Float64}}}(([0.003101770562619906;;], [10.393833158673782;;])), RecursiveArrayTools.ArrayPartition{Float64, Tuple{Matrix{Float64}, Matrix{Float64}}}(([0.003101770562619906;;], [10.393833158673782;;]))]}
 {:Time = [0.0, 8351.500948888295], :OutputVelocity = [[-0.0033941074835368125;;], [-0.0005183610570836693;;]], :OutputPosition = [[9.448630623128851;;], [3.7122680138595845;;]], :OutputKineticEnergy = [0.07349864435103717, 0.0017143239001597547], :OutputDynamicsVariables = RecursiveArrayTools.ArrayPartition{Float64, Tuple{Matrix{Float64}, Matrix{Float64}}}[RecursiveArrayTools.ArrayPartition{Float64, Tuple{Matrix{Float64}, Matrix{Float64}}}(([-0.0033941074835368125;;], [9.448630623128851;;])), RecursiveArrayTools.ArrayPartition{Float64, Tuple{Matrix{Float64}, Matrix{Float64}}}(([-0.0005183610570836693;;], [3.7122680138595845;;]))]}
 {:Time = [0.0, 8351.500948888295], :OutputVelocity = [[-0.0033941074835368125;;], [-0.000937834821437652;;]], :OutputPosition = [[9.448630623128851;;], [3.149035876931468;;]], :OutputKineticEnergy = [0.07349864435103717, 0.005611524378125343], :OutputDynamicsVariables = RecursiveArrayTools.ArrayPartition{Float64, Tuple{Matrix{Float64}, Matrix{Float64}}}[RecursiveArrayTools.ArrayPartition{Float64, Tuple{Matrix{Float64}, Matrix{Float64}}}(([-0.0033941074835368125;;], [9.448630623128851;;])), RecursiveArrayTools.ArrayPartition{Float64, Tuple{Matrix{Float64}, Matrix{Float64}}}(([-0.000937834821437652;;], [3.149035876931468;;]))]}
 {:Time = [0.0, 4206.898150604404, 4206.898150604404], :OutputVelocity = [[-0.0033941074835368125;;], [0.0030241025340672447;;], [0.0030241025340672447;;]], :OutputPosition = [[9.448630623128851;;], [10.393815978889709;;], [10.393815978889709;;]], :OutputKineticEnergy = [0.07349864435103717, 0.058347354595318486, 0.058347354595318486], :OutputDynamicsVariables = RecursiveArrayTools.ArrayPartition{Float64, Tuple{Matrix{Float64}, Matrix{Float64}}}[RecursiveArrayTools.ArrayPartition{Float64, Tuple{Matrix{Float64}, Matrix{Float64}}}(([-0.0033941074835368125;;], [9.448630623128851;;])), RecursiveArrayTools.ArrayPartition{Float64, Tuple{Matrix{Float64}, Matrix{Float64}}}(([0.0030241025340672447;;], [10.393815978889709;;])), RecursiveArrayTools.ArrayPartition{Float64, Tuple{Matrix{Float64}, Matrix{Float64}}}(([0.0030241025340672447;;], [10.393815978889709;;]))]}
 {:Time = [0.0, 8351.500948888295], :OutputVelocity = [[-0.0033941074835368125;;], [-0.00020084611683819292;;]], :OutputPosition = [[9.448630623128851;;], [3.759249858792797;;]], :OutputKineticEnergy = [0.07349864435103717, 0.0002573682829776626], :OutputDynamicsVariables = RecursiveArrayTools.ArrayPartition{Float64, Tuple{Matrix{Float64}, Matrix{Float64}}}[RecursiveArrayTools.ArrayPartition{Float64, Tuple{Matrix{Float64}, Matrix{Float64}}}(([-0.0033941074835368125;;], [9.448630623128851;;])), RecursiveArrayTools.ArrayPartition{Float64, Tuple{Matrix{Float64}, Matrix{Float64}}}(([-0.00020084611683819292;;], [3.759249858792797;;]))]}
 {:Time = [0.0, 4235.423698205766, 4235.423698205766], :OutputVelocity = [[-0.0033941074835368125;;], [0.0030712458819135073;;], [0.0030712458819135073;;]], :OutputPosition = [[9.448630623128851;;], [10.394250624970008;;], [10.394250624970008;;]], :OutputKineticEnergy = [0.07349864435103717, 0.06018071185202984, 0.06018071185202984], :OutputDynamicsVariables = RecursiveArrayTools.ArrayPartition{Float64, Tuple{Matrix{Float64}, Matrix{Float64}}}[RecursiveArrayTools.ArrayPartition{Float64, Tuple{Matrix{Float64}, Matrix{Float64}}}(([-0.0033941074835368125;;], [9.448630623128851;;])), RecursiveArrayTools.ArrayPartition{Float64, Tuple{Matrix{Float64}, Matrix{Float64}}}(([0.0030712458819135073;;], [10.394250624970008;;])), RecursiveArrayTools.ArrayPartition{Float64, Tuple{Matrix{Float64}, Matrix{Float64}}}(([0.0030712458819135073;;], [10.394250624970008;;]))]}

println(typeof(atoms))

Atoms{Float64}

Broadened Classical Master Equation

Finally, we can also use the Broadened Classical Master Equation (BCME) which is governed by the two diabatic energy curves and the hybridisation function Γ themselves instead of a Newns-Anderson Hamiltonian. We therefore have to adapt the electronic subspace and thence the total space which is done in the first and second line, respectively. The rest goes as ever. Albeit one might consider BCME calculations to be gimmicky, an advantageous attribute of them is that they are computationally extremely efficient and can therefore serve as a nice test system for analysis tools which require numerical statics on the one hand, and electronically non-adiabatic effects on the other hand. I definitely learned to love them ;P

### Calculation of BCME
electronic_distribution = PureState(1, Diabatic())
dist = electronic_distribution * nuclear_distribution

sim = Simulation{BCME}(atoms, thossmodel;
    temperature, bandwidth=bandwidth)


result = run_dynamics(sim,
                    (0.0, tcut),
                    dist;
                    trajectories=ntrajs,
                    callback=terminate,
                    output=output,
                    saveat=[0,tcut],
                    dt,
                    abstol = 1e-6,
                    reltol = 1e-6,
                   )

[ Info: Sampling randomly from provided distribution.
[ Info: Pre-compiled dynamics in 10.725150084 seconds.
[ Info: Performing 10 trajectories.
[ Info: Sampling randomly from provided distribution.
[ Info: Finished after 0.031583292 seconds.

10-element Vector{Dictionaries.Dictionary{Symbol, Any}}:
 {:Time = [0.0, 8351.500948888295], :OutputVelocity = [[-0.0033941074835368125;;], [0.00011269537185680292;;]], :OutputPosition = [[9.448630623128851;;], [3.6183214707408227;;]], :OutputKineticEnergy = [0.07349864435103717, 8.102896806551502e-5], :OutputDynamicsVariables = NQCDynamics.DynamicsMethods.SurfaceHoppingMethods.SurfaceHoppingVariables{Float64, RecursiveArrayTools.ArrayPartition{Float64, Tuple{Matrix{Float64}, Matrix{Float64}, Vector{Float64}}}, @NamedTuple{r::Int64, v::Int64, state::Int64}}[[9.448630623128851, -0.0033941074835368125, 1.0], [3.6183214707408227, 0.00011269537185680292, 1.0]]}
 {:Time = [0.0, 5710.510668421319, 5710.510668421319], :OutputVelocity = [[-0.0033941074835368125;;], [0.0028766890190744542;;], [0.0028766890190744542;;]], :OutputPosition = [[9.448630623128851;;], [10.509327924996953;;], [10.509327924996953;;]], :OutputKineticEnergy = [0.07349864435103717, 0.0527975751848534, 0.0527975751848534], :OutputDynamicsVariables = NQCDynamics.DynamicsMethods.SurfaceHoppingMethods.SurfaceHoppingVariables{Float64, RecursiveArrayTools.ArrayPartition{Float64, Tuple{Matrix{Float64}, Matrix{Float64}, Vector{Float64}}}, @NamedTuple{r::Int64, v::Int64, state::Int64}}[[9.448630623128851, -0.0033941074835368125, 1.0], [10.509327924996953, 0.0028766890190744542, 2.0], [10.509327924996953, 0.0028766890190744542, 2.0]]}
 {:Time = [0.0, 6110.643955889959, 6110.643955889959], :OutputVelocity = [[-0.0033941074835368125;;], [0.0027265091980518664;;], [0.0027265091980518664;;]], :OutputPosition = [[9.448630623128851;;], [11.448835660121953;;], [11.448835660121953;;]], :OutputKineticEnergy = [0.07349864435103717, 0.04742879386374899, 0.04742879386374899], :OutputDynamicsVariables = NQCDynamics.DynamicsMethods.SurfaceHoppingMethods.SurfaceHoppingVariables{Float64, RecursiveArrayTools.ArrayPartition{Float64, Tuple{Matrix{Float64}, Matrix{Float64}, Vector{Float64}}}, @NamedTuple{r::Int64, v::Int64, state::Int64}}[[9.448630623128851, -0.0033941074835368125, 1.0], [11.448835660121953, 0.0027265091980518664, 2.0], [11.448835660121953, 0.0027265091980518664, 2.0]]}
 {:Time = [0.0, 5441.992254139246, 5441.992254139246], :OutputVelocity = [[-0.0033941074835368125;;], [0.002826213639873677;;], [0.002826213639873677;;]], :OutputPosition = [[9.448630623128851;;], [10.585515200826077;;], [10.585515200826077;;]], :OutputKineticEnergy = [0.07349864435103717, 0.05096102121477403, 0.05096102121477403], :OutputDynamicsVariables = NQCDynamics.DynamicsMethods.SurfaceHoppingMethods.SurfaceHoppingVariables{Float64, RecursiveArrayTools.ArrayPartition{Float64, Tuple{Matrix{Float64}, Matrix{Float64}, Vector{Float64}}}, @NamedTuple{r::Int64, v::Int64, state::Int64}}[[9.448630623128851, -0.0033941074835368125, 1.0], [10.585515200826077, 0.002826213639873677, 2.0], [10.585515200826077, 0.002826213639873677, 2.0]]}
 {:Time = [0.0, 4776.181735345739, 4776.181735345739], :OutputVelocity = [[-0.0033941074835368125;;], [0.0034105201618294455;;], [0.0034105201618294455;;]], :OutputPosition = [[9.448630623128851;;], [11.208942998061884;;], [11.208942998061884;;]], :OutputKineticEnergy = [0.07349864435103717, 0.0742111888119235, 0.0742111888119235], :OutputDynamicsVariables = NQCDynamics.DynamicsMethods.SurfaceHoppingMethods.SurfaceHoppingVariables{Float64, RecursiveArrayTools.ArrayPartition{Float64, Tuple{Matrix{Float64}, Matrix{Float64}, Vector{Float64}}}, @NamedTuple{r::Int64, v::Int64, state::Int64}}[[9.448630623128851, -0.0033941074835368125, 1.0], [11.208942998061884, 0.0034105201618294455, 2.0], [11.208942998061884, 0.0034105201618294455, 2.0]]}
 {:Time = [0.0, 4036.239702391223, 4036.239702391223], :OutputVelocity = [[-0.0033941074835368125;;], [0.003528624476762427;;], [0.003528624476762427;;]], :OutputPosition = [[9.448630623128851;;], [10.67617270886745;;], [10.67617270886745;;]], :OutputKineticEnergy = [0.07349864435103717, 0.07943996257082518, 0.07943996257082518], :OutputDynamicsVariables = NQCDynamics.DynamicsMethods.SurfaceHoppingMethods.SurfaceHoppingVariables{Float64, RecursiveArrayTools.ArrayPartition{Float64, Tuple{Matrix{Float64}, Matrix{Float64}, Vector{Float64}}}, @NamedTuple{r::Int64, v::Int64, state::Int64}}[[9.448630623128851, -0.0033941074835368125, 1.0], [10.67617270886745, 0.003528624476762427, 2.0], [10.67617270886745, 0.003528624476762427, 2.0]]}
 {:Time = [0.0, 4776.181735345739, 4776.181735345739], :OutputVelocity = [[-0.0033941074835368125;;], [0.0034105201618294455;;], [0.0034105201618294455;;]], :OutputPosition = [[9.448630623128851;;], [11.208942998061884;;], [11.208942998061884;;]], :OutputKineticEnergy = [0.07349864435103717, 0.0742111888119235, 0.0742111888119235], :OutputDynamicsVariables = NQCDynamics.DynamicsMethods.SurfaceHoppingMethods.SurfaceHoppingVariables{Float64, RecursiveArrayTools.ArrayPartition{Float64, Tuple{Matrix{Float64}, Matrix{Float64}, Vector{Float64}}}, @NamedTuple{r::Int64, v::Int64, state::Int64}}[[9.448630623128851, -0.0033941074835368125, 1.0], [11.208942998061884, 0.0034105201618294455, 2.0], [11.208942998061884, 0.0034105201618294455, 2.0]]}
 {:Time = [0.0, 6190.869315291103, 6190.869315291103], :OutputVelocity = [[-0.0033941074835368125;;], [0.002716746803962239;;], [0.002716746803962239;;]], :OutputPosition = [[9.448630623128851;;], [11.427785760321408;;], [11.427785760321408;;]], :OutputKineticEnergy = [0.07349864435103717, 0.04708975987306522, 0.04708975987306522], :OutputDynamicsVariables = NQCDynamics.DynamicsMethods.SurfaceHoppingMethods.SurfaceHoppingVariables{Float64, RecursiveArrayTools.ArrayPartition{Float64, Tuple{Matrix{Float64}, Matrix{Float64}, Vector{Float64}}}, @NamedTuple{r::Int64, v::Int64, state::Int64}}[[9.448630623128851, -0.0033941074835368125, 1.0], [11.427785760321408, 0.002716746803962239, 2.0], [11.427785760321408, 0.002716746803962239, 2.0]]}
 {:Time = [0.0, 8351.500948888295], :OutputVelocity = [[-0.0033941074835368125;;], [0.0005291239536066064;;]], :OutputPosition = [[9.448630623128851;;], [3.1313588097958585;;]], :OutputKineticEnergy = [0.07349864435103717, 0.0017862530832669066], :OutputDynamicsVariables = NQCDynamics.DynamicsMethods.SurfaceHoppingMethods.SurfaceHoppingVariables{Float64, RecursiveArrayTools.ArrayPartition{Float64, Tuple{Matrix{Float64}, Matrix{Float64}, Vector{Float64}}}, @NamedTuple{r::Int64, v::Int64, state::Int64}}[[9.448630623128851, -0.0033941074835368125, 1.0], [3.1313588097958585, 0.0005291239536066064, 1.0]]}
 {:Time = [0.0, 4275.104951975047, 4275.104951975047], :OutputVelocity = [[-0.0033941074835368125;;], [0.003910988269645446;;], [0.003910988269645446;;]], :OutputPosition = [[9.448630623128851;;], [10.715278214687448;;], [10.715278214687448;;]], :OutputKineticEnergy = [0.07349864435103717, 0.09758906856443929, 0.09758906856443929], :OutputDynamicsVariables = NQCDynamics.DynamicsMethods.SurfaceHoppingMethods.SurfaceHoppingVariables{Float64, RecursiveArrayTools.ArrayPartition{Float64, Tuple{Matrix{Float64}, Matrix{Float64}, Vector{Float64}}}, @NamedTuple{r::Int64, v::Int64, state::Int64}}[[9.448630623128851, -0.0033941074835368125, 1.0], [10.715278214687448, 0.003910988269645446, 2.0], [10.715278214687448, 0.003910988269645446, 2.0]]}