Scattering of a mono-energetic particle from a surface

### Loading required packages

using NQCDynamics, NQCCalculators
using Unitful, UnitfulAtomic
using Statistics

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)

### Defining termination criterion

(NQCDynamics.DynamicsOutputs.OutputVelocity, NQCDynamics.DynamicsOutputs.OutputPosition, NQCDynamics.DynamicsOutputs.OutputKineticEnergy)

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 ingeredients 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. (not so satisfied with the explanation here) 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}(Main.termination_condition, SciMLBase.terminate!, SciMLBase.INITIALIZE_DEFAULT, SciMLBase.FINALIZE_DEFAULT, Bool[1, 1], nothing)

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’s purpose is only to accelerate the actual simulations (blame Diffeq.jl for this lol). The third line launches our IESH simulations. :)

### Run single Ehrenfest calculation

sim = Simulation{EhrenfestNA}(atoms, model)

kick = run_dynamics(sim,
                    (0.0, dt),
                    dist;
                    trajectories=1,
                    callback=terminate,
                    output=output,
                    dt
                   )

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

┌ Warning: No friction evaluation algorithm has been specified, the friction tensor will be empty

└ @ NQCCalculators ~/.julia/packages/NQCCalculators/t5bm5/src/Caches.jl:425

[ Info: Pre-compiled dynamics in 6.671770086 seconds.

[ Info: Performing 1 trajectory.

[ Info: Finished after 0.00196775 seconds.

[ Info: Pre-compiled dynamics in 1.096891014 seconds.

[ Info: Performing 1 trajectory.

[ Info: Finished after 1.655033598 seconds.

4-element Dictionaries.Dictionary{Symbol, Any}:
                :Time │ [0.0, 4553.338859154286, 4553.338859154286]
      :OutputVelocity │ [[-0.0033941074835368125;;], [0.002989208858542591;;], …
      :OutputPosition │ [[9.448630623128851;;], [10.394314983127463;;], [10.394…
 :OutputKineticEnergy │ [0.07349864435103717, 0.057008638273709694, 0.057008638…

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)

kick = run_dynamics(sim,
                    (0.0, dt),
                    dist;
                    trajectories=1,
                    callback=terminate,
                    output=output,
                    dt
                   )

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

┌ Warning: No friction evaluation algorithm has been specified, the friction tensor will be empty

└ @ NQCCalculators ~/.julia/packages/NQCCalculators/t5bm5/src/Caches.jl:425

[ Info: Pre-compiled dynamics in 8.984125601 seconds.

[ Info: Performing 1 trajectory.

[ Info: Finished after 0.00293141 seconds.

[ Info: Sampling randomly from provided distribution.

[ Info: Pre-compiled dynamics in 4.484304186 seconds.

[ Info: Performing 10 trajectories.

[ Info: Sampling randomly from provided distribution.

[ Info: Finished after 30.387933773 seconds.

10-element Vector{Dictionaries.Dictionary{Symbol, Any}}:
 {:Time = [0.0, 4755.498174763943, 4755.498174763943], :OutputVelocity = [[-0.0033941074835368125;;], [0.003390554049472842;;], [0.003390554049472842;;]], :OutputPosition = [[9.448630623128851;;], [10.394067981186835;;], [10.394067981186835;;]], :OutputKineticEnergy = [0.07349864435103717, 0.07334482725981725, 0.07334482725981725]}
 {:Time = [0.0, 6908.143484333439, 6908.143484333439], :OutputVelocity = [[-0.0033941074835368125;;], [0.0019137599364812902;;], [0.0019137599364812902;;]], :OutputPosition = [[9.448630623128851;;], [10.393977928595763;;], [10.393977928595763;;]], :OutputKineticEnergy = [0.07349864435103717, 0.023367005642973365, 0.023367005642973365]}
 {:Time = [0.0, 5149.894876382783, 5149.894876382783], :OutputVelocity = [[-0.0033941074835368125;;], [0.003391372268450932;;], [0.003391372268450932;;]], :OutputPosition = [[9.448630623128851;;], [10.39460863098261;;], [10.39460863098261;;]], :OutputKineticEnergy = [0.07349864435103717, 0.0733802311318035, 0.0733802311318035]}
 {:Time = [0.0, 4779.476171298422, 4779.476171298422], :OutputVelocity = [[-0.0033941074835368125;;], [0.003389471655311061;;], [0.003389471655311061;;]], :OutputPosition = [[9.448630623128851;;], [10.39457013031362;;], [10.39457013031362;;]], :OutputKineticEnergy = [0.07349864435103717, 0.07329800580242475, 0.07329800580242475]}
 {:Time = [0.0, 4755.498174763943, 4755.498174763943], :OutputVelocity = [[-0.0033941074835368125;;], [0.003390554049472842;;], [0.003390554049472842;;]], :OutputPosition = [[9.448630623128851;;], [10.394067981186835;;], [10.394067981186835;;]], :OutputKineticEnergy = [0.07349864435103717, 0.07334482725981725, 0.07334482725981725]}
 {:Time = [0.0, 4755.498174763943, 4755.498174763943], :OutputVelocity = [[-0.0033941074835368125;;], [0.003390554049472842;;], [0.003390554049472842;;]], :OutputPosition = [[9.448630623128851;;], [10.394067981186835;;], [10.394067981186835;;]], :OutputKineticEnergy = [0.07349864435103717, 0.07334482725981725, 0.07334482725981725]}
 {:Time = [0.0, 6867.215524731484, 6867.215524731484], :OutputVelocity = [[-0.0033941074835368125;;], [0.0019159104212118413;;], [0.0019159104212118413;;]], :OutputPosition = [[9.448630623128851;;], [10.394006780484904;;], [10.394006780484904;;]], :OutputKineticEnergy = [0.07349864435103717, 0.02341954997830035, 0.02341954997830035]}
 {:Time = [0.0, 6920.1324826006785, 6920.1324826006785], :OutputVelocity = [[-0.0033941074835368125;;], [0.0019131900032905525;;], [0.0019131900032905525;;]], :OutputPosition = [[9.448630623128851;;], [10.393625057469878;;], [10.393625057469878;;]], :OutputKineticEnergy = [0.07349864435103717, 0.02335308994876223, 0.02335308994876223]}
 {:Time = [0.0, 6932.121480867918, 6932.121480867918], :OutputVelocity = [[-0.0033941074835368125;;], [0.001912343383677912;;], [0.001912343383677912;;]], :OutputPosition = [[9.448630623128851;;], [10.39359885110002;;], [10.39359885110002;;]], :OutputKineticEnergy = [0.07349864435103717, 0.02333242623069217, 0.02333242623069217]}
 {:Time = [0.0, 5371.898051193388, 5371.898051193388], :OutputVelocity = [[-0.0033941074835368125;;], [0.0033903822208739354;;], [0.0033903822208739354;;]], :OutputPosition = [[9.448630623128851;;], [10.393584023900164;;], [10.393584023900164;;]], :OutputKineticEnergy = [0.07349864435103717, 0.07333739341898975, 0.07333739341898975]}

### 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 7.318874167 seconds.

[ Info: Performing 10 trajectories.

[ Info: Sampling randomly from provided distribution.

[ Info: Finished after 11.450519627 seconds.

10-element Vector{Dictionaries.Dictionary{Symbol, Any}}:
 {:Time = [0.0, 4034.918037529455, 4034.918037529455], :OutputVelocity = [[-0.0033941074835368125;;], [0.003193553989325685;;], [0.003193553989325685;;]], :OutputPosition = [[9.448630623128851;;], [10.394478255783389;;], [10.394478255783389;;]], :OutputKineticEnergy = [0.07349864435103717, 0.06506938041276734, 0.06506938041276734]}
 {:Time = [0.0, 4160.182398735506, 4160.182398735506], :OutputVelocity = [[-0.0033941074835368125;;], [0.0031109426375222366;;], [0.0031109426375222366;;]], :OutputPosition = [[9.448630623128851;;], [10.39424342624923;;], [10.39424342624923;;]], :OutputKineticEnergy = [0.07349864435103717, 0.06174647260912539, 0.06174647260912539]}
 {:Time = [0.0, 8351.500948888295], :OutputVelocity = [[-0.0033941074835368125;;], [-0.0008992790747204714;;]], :OutputPosition = [[9.448630623128851;;], [3.630990717748359;;]], :OutputKineticEnergy = [0.07349864435103717, 0.005159612926114991]}
 {:Time = [0.0, 4202.35059953752, 4202.35059953752], :OutputVelocity = [[-0.0033941074835368125;;], [0.003067233235301346;;], [0.003067233235301346;;]], :OutputPosition = [[9.448630623128851;;], [10.394359440789202;;], [10.394359440789202;;]], :OutputKineticEnergy = [0.07349864435103717, 0.060023559877279885, 0.060023559877279885]}
 {:Time = [0.0, 8351.500948888295], :OutputVelocity = [[-0.0033941074835368125;;], [0.0031150647333923914;;]], :OutputPosition = [[9.448630623128851;;], [8.07614955243684;;]], :OutputKineticEnergy = [0.07349864435103717, 0.061910213015883495]}
 {:Time = [0.0, 8351.500948888295], :OutputVelocity = [[-0.0033941074835368125;;], [-0.0010694066620026246;;]], :OutputPosition = [[9.448630623128851;;], [3.607434621185741;;]], :OutputKineticEnergy = [0.07349864435103717, 0.007296488742024349]}
 {:Time = [0.0, 8351.500948888295], :OutputVelocity = [[-0.0033941074835368125;;], [0.0009207799154834855;;]], :OutputPosition = [[9.448630623128851;;], [3.5398661907586115;;]], :OutputKineticEnergy = [0.07349864435103717, 0.005409284473658406]}
 {:Time = [0.0, 4326.788133276797, 4326.788133276797], :OutputVelocity = [[-0.0033941074835368125;;], [0.003122846551089521;;], [0.003122846551089521;;]], :OutputPosition = [[9.448630623128851;;], [10.393550854060857;;], [10.393550854060857;;]], :OutputKineticEnergy = [0.07349864435103717, 0.06221991814122039, 0.06221991814122039]}
 {:Time = [0.0, 5793.166645480178, 5793.166645480178], :OutputVelocity = [[-0.0033941074835368125;;], [0.0030642385556808034;;], [0.0030642385556808034;;]], :OutputPosition = [[9.448630623128851;;], [10.394637039995343;;], [10.394637039995343;;]], :OutputKineticEnergy = [0.07349864435103717, 0.059906409619686254, 0.059906409619686254]}
 {:Time = [0.0, 4223.84811367188, 4223.84811367188], :OutputVelocity = [[-0.0033941074835368125;;], [0.003054213014402486;;], [0.003054213014402486;;]], :OutputPosition = [[9.448630623128851;;], [10.393969614940438;;], [10.393969614940438;;]], :OutputKineticEnergy = [0.07349864435103717, 0.05951504865810761, 0.05951504865810761]}

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
                   )

[ Info: Sampling randomly from provided distribution.

[ Info: Pre-compiled dynamics in 23.242684995 seconds.

[ Info: Performing 10 trajectories.

[ Info: Sampling randomly from provided distribution.

[ Info: Finished after 2.196294334 seconds.

10-element Vector{Dictionaries.Dictionary{Symbol, Any}}:
 {:Time = [0.0, 4484.063935063259, 4484.063935063259], :OutputVelocity = [[-0.0033941074835368125;;], [0.0035949399176046074;;], [0.0035949399176046074;;]], :OutputPosition = [[9.448630623128851;;], [10.590909554623925;;], [10.590909554623925;;]], :OutputKineticEnergy = [0.07349864435103717, 0.08245394115227975, 0.08245394115227975]}
 {:Time = [0.0, 4812.37182697698, 4812.37182697698], :OutputVelocity = [[-0.0033941074835368125;;], [0.003294608187157551;;], [0.003294608187157551;;]], :OutputPosition = [[9.448630623128851;;], [10.996908751639259;;], [10.996908751639259;;]], :OutputKineticEnergy = [0.07349864435103717, 0.069252537773449, 0.069252537773449]}
 {:Time = [0.0, 5844.0894748216915, 5844.0894748216915], :OutputVelocity = [[-0.0033941074835368125;;], [0.00245585220525486;;], [0.00245585220525486;;]], :OutputPosition = [[9.448630623128851;;], [11.212913643901299;;], [11.212913643901299;;]], :OutputKineticEnergy = [0.07349864435103717, 0.03847978177923337, 0.03847978177923337]}
 {:Time = [0.0, 5281.757657152533, 5281.757657152533], :OutputVelocity = [[-0.0033941074835368125;;], [0.0031512592107842234;;], [0.0031512592107842234;;]], :OutputPosition = [[9.448630623128851;;], [11.199215174953348;;], [11.199215174953348;;]], :OutputKineticEnergy = [0.07349864435103717, 0.06335726222062493, 0.06335726222062493]}
 {:Time = [0.0, 5056.621156805089, 5056.621156805089], :OutputVelocity = [[-0.0033941074835368125;;], [0.0033187663627387903;;], [0.0033187663627387903;;]], :OutputPosition = [[9.448630623128851;;], [10.97440389512251;;], [10.97440389512251;;]], :OutputKineticEnergy = [0.07349864435103717, 0.07027186916569567, 0.07027186916569567]}
 {:Time = [0.0, 5411.944547755553, 5411.944547755553], :OutputVelocity = [[-0.0033941074835368125;;], [0.0030283382572178245;;], [0.0030283382572178245;;]], :OutputPosition = [[9.448630623128851;;], [11.19285857013145;;], [11.19285857013145;;]], :OutputKineticEnergy = [0.07349864435103717, 0.058510918045305955, 0.058510918045305955]}
 {:Time = [0.0, 5546.172057013184, 5546.172057013184], :OutputVelocity = [[-0.0033941074835368125;;], [0.002856339750440248;;], [0.002856339750440248;;]], :OutputPosition = [[9.448630623128851;;], [10.921674310694197;;], [10.921674310694197;;]], :OutputKineticEnergy = [0.07349864435103717, 0.05205325281345526, 0.05205325281345526]}
 {:Time = [0.0, 5883.752328813726, 5883.752328813726], :OutputVelocity = [[-0.0033941074835368125;;], [0.0030119621079222466;;], [0.0030119621079222466;;]], :OutputPosition = [[9.448630623128851;;], [11.265588538394013;;], [11.265588538394013;;]], :OutputKineticEnergy = [0.07349864435103717, 0.05787981762351772, 0.05787981762351772]}
 {:Time = [0.0, 5555.143439497792, 5555.143439497792], :OutputVelocity = [[-0.0033941074835368125;;], [0.0030619046335191315;;], [0.0030619046335191315;;]], :OutputPosition = [[9.448630623128851;;], [10.831537877460322;;], [10.831537877460322;;]], :OutputKineticEnergy = [0.07349864435103717, 0.05981518718533298, 0.05981518718533298]}
 {:Time = [0.0, 5495.630579476603, 5495.630579476603], :OutputVelocity = [[-0.0033941074835368125;;], [0.0029345289845891256;;], [0.0029345289845891256;;]], :OutputPosition = [[9.448630623128851;;], [11.097839346861583;;], [11.097839346861583;;]], :OutputKineticEnergy = [0.07349864435103717, 0.054942061798054644, 0.054942061798054644]}