Scattering probabilities for TullyModelTwo

In this section we aim to reproduce the results of Fig. 5 from [1]. This figure presents the scattering outcomes when a particle interacts with Tully's model 2 with an increasing magnitude of incident kinetic energy.

To reproduce this figure, first, let's set up our system parameters:

using NQCDynamics

sim = Simulation{FSSH}(Atoms(2000), TullyModelTwo())

Simulation{FSSH{Float64}}:
  Atoms{Float64}([:X], [0], [2000.0])
  TullyModelTwo{Float64, Float64, Float64, Float64, Float64}
  a: Float64 0.1
  b: Float64 0.28
  c: Float64 0.015
  d: Float64 0.06
  e: Float64 0.05

!!! note Atomic units

Recall that all of are units are atomic by default, this mass of 2000 is similar to
that of a hydrogen atom.

Each data point in the figure is obtained from an ensemble average of trajectories. We can use our Ensembles setup to run a set of trajectories for every single momentum value. Firstly, we can prepare the parts that will be the same for every ensemble:

using ComponentArrays: ComponentVector

output = OutputStateResolvedScattering1D(sim, :adiabatic)

OutputStateResolvedScattering1D{Simulation{FSSH{Float64}, NQCDynamics.Calculators.DiabaticCalculator{Float64, TullyModelTwo{Float64, Float64, Float64, Float64, Float64}, 2, 4}, Float64, Unitful.Quantity{Int64, 𝚯, Unitful.FreeUnits{(K,), 𝚯, nothing}}, InfiniteCell}}(Simulation{FSSH{Float64}}:
  Atoms{Float64}([:X], [0], [2000.0])
  TullyModelTwo{Float64, Float64, Float64, Float64, Float64}
  a: Float64 0.1
  b: Float64 0.28
  c: Float64 0.015
  d: Float64 0.06
  e: Float64 0.05
, :adiabatic)

Here, we are using the OutputStateResolvedScattering1D along with the MeanReduction which will give us the average scattering outcome from the entire ensemble. Each trajectory outputs the scattering outcome along with its final adiabatic state, and the reduction computes the average over all trajectories.

Next, we can choose how many trajectories we want to perform for each ensemble, and choose the range of momentum values:

trajectories = 500
momenta = 9:2:50

9:2:49

Since each ensemble requires different initial conditions, we will specify the trajectory timespan and the distribution inside the loop. Before the loop begins, we will create an empty list to store the results, and append to this list after every iteration. The distribution we create produces initial conditions where each trajectory has momentum k and starts at a position of -5.

result = []
for k in momenta # Iterate through each momentum value
    v = k / sim.atoms.masses[1] # Starting velocity
    r = -5 # Starting position
    tspan = (0, 2abs(r)/v)
    distribution = DynamicalDistribution(v, -5, size(sim)) * PureState(1, Adiabatic())

    out = run_dynamics(sim, tspan, distribution;
        saveat=tspan[end], trajectories, output, reduction=MeanReduction()
    )

    push!(result, out[:OutputStateResolvedScattering1D])
end

result

21-element Vector{Any}:
 ComponentVector{Float64}(reflection = [0.008, 0.0], transmission = [0.992, 0.0])
 ComponentVector{Float64}(reflection = [0.02, 0.0], transmission = [0.976, 0.004])
 ComponentVector{Float64}(reflection = [0.0, 0.0], transmission = [0.878, 0.122])
 ComponentVector{Float64}(reflection = [0.0, 0.0], transmission = [0.878, 0.122])
 ComponentVector{Float64}(reflection = [0.0, 0.0], transmission = [0.924, 0.076])
 ComponentVector{Float64}(reflection = [0.0, 0.0], transmission = [0.96, 0.04])
 ComponentVector{Float64}(reflection = [0.0, 0.0], transmission = [0.932, 0.068])
 ComponentVector{Float64}(reflection = [0.0, 0.0], transmission = [0.83, 0.17])
 ComponentVector{Float64}(reflection = [0.0, 0.0], transmission = [0.656, 0.344])
 ComponentVector{Float64}(reflection = [0.0, 0.0], transmission = [0.486, 0.514])
 ⋮
 ComponentVector{Float64}(reflection = [0.0, 0.0], transmission = [0.324, 0.676])
 ComponentVector{Float64}(reflection = [0.0, 0.0], transmission = [0.41, 0.59])
 ComponentVector{Float64}(reflection = [0.0, 0.0], transmission = [0.51, 0.49])
 ComponentVector{Float64}(reflection = [0.0, 0.0], transmission = [0.634, 0.366])
 ComponentVector{Float64}(reflection = [0.0, 0.0], transmission = [0.75, 0.25])
 ComponentVector{Float64}(reflection = [0.0, 0.0], transmission = [0.876, 0.124])
 ComponentVector{Float64}(reflection = [0.0, 0.0], transmission = [0.946, 0.054])
 ComponentVector{Float64}(reflection = [0.0, 0.0], transmission = [0.978, 0.022])
 ComponentVector{Float64}(reflection = [0.0, 0.0], transmission = [1.0, 0.0])

Now we can plot our simulation results. We format this plot to match Fig. 3 from [18] which also reproduces Fig. 5 from [1]. We manage to reproduce the FSSH results quite accurately by visual comparison, though a larger number of trajectories would lead to better convergence. Since all of the examples run during the documentation build, we use a minimal number of trajectories to optimise the build time.

using Plots

plot(
    xlabel="Incidence momentum / a.u.",
    ylabel="Scattering probability"
)

r1 = [r.reflection[1] for r in result]
t1 = [r.transmission[1] for r in result]
t2 = [r.transmission[2] for r in result]

scatter!(momenta, r1; label="R1", color=:red)
scatter!(momenta, t1; label="T1", color=:green)
scatter!(momenta, t2; label="T2", color=:blue)

As in [18], R1, T1, T2 refer to reflection on state 1, transmission on state 1 and transmission on state 2 respectively. For this model, surface hopping is successful in closely approximating the exact quantum result, especially at higher momentum values. Refer to [1] and [18] for a detailed discussion of the results.