GitHub

Semiclassical EBK quantisation

In surface science, it is often of interest to investigate how collisions with surfaces can perturb the quantum states of molecules. In particular, for diatomic molecules, the rotational and vibrational quantum numbers can undergo significant changes when the molecule impacts the surface.

Einstein-Brillouinn-Keller (EBK) quantisation allows for a semiclassical investigation into these phenomena by providing a link between the quantum numbers and classical positions and velocities. The quantisation procedure allows the user to generate a classical distribution with a given set of quantum numbers, then perform semiclassical dynamics and extract the quantum numbers at the end by reversing the procedure.

These three steps can be applied to give insight into the processes taking place during surface scattering and allow us to attempt to predict the experimentally observed change in the quantum numbers.

A detailed yet approachable description of the theory is given by [4] so we shall not delve into the theory here. Briefly, the procedure for a diatomic molecule involves an optimisation process to find the bounds of an integral, then computing the integral to obtain the vibrational quantum number. The rotational quantum number comes directly from the classical angular momentum of the molecule.

Configurations can be generated by randomly selecting bond lengths from the appropriate probability distribution and selecting a matching radial velocity.

Example

In this example we will create a quantised distribution suitable for use as initial conditions for hydrogen scattering simulations.

The simulation can be set up in the usual way, by specifying the atoms along with the model and the simulation cell.

using NQCDynamics
using Unitful, UnitfulAtomic

atoms = Atoms([:H, :H])
model = DiatomicHarmonic()
cell = PeriodicCell(austrip.([5.883 -2.942 0; 0 5.095 0; 0 0 20] .* u"Å"))
sim = Simulation(atoms, model; cell=cell)
Simulation{Classical}:
  Atoms{Float64}([:H, :H], [1, 1], [1837.4715941070515, 1837.4715941070515])
  DiatomicHarmonic
  r₀: Float64 1.0

The distribution is generated using the QuantisedDiatomic.generate_configurations function. We have to provide the desired vibrational ν and rotational J quantum numbers, along with the number of samples and some other options as keyword arguments. In addition to the rotational and vibrational energy we have applied a translational impulse of 1 eV and positioned the molecule at a height of 10 bohr.

using NQCDynamics.InitialConditions: QuantisedDiatomic

ν, J = 2, 0
nsamples = 150

configurations = QuantisedDiatomic.generate_configurations(sim, ν, J;
    samples=nsamples, translational_energy=1u"eV", height=10)
150-element Vector{Tuple{Matrix{Float64}, Matrix{Float64}}}:
 ([-0.000468154383111016 0.000468154383111016; -0.0057044541955887365 0.0057044541955887365; -0.007856726622633696 -0.0010875319368652868], [2.3966809227234855 2.3227791328788325; 6.335870543651092 5.435378333497253; 10.267142432674056 9.732857567325944])
 ([-0.0026656809485983294 0.0026656809485983294; 0.004419927672108907 -0.004419927672108907; -0.00037933024289343344 -0.008564928316605549], [7.238777327384728 7.6734332904909195; 6.051923670077572 5.331226798912108; 10.33367824984835 9.66632175015165])
 ([0.004384765362670697 -0.004384765362670697; -0.0037698812088642363 0.0037698812088642363; -0.007333388769276697 -0.0016108697902222852], [2.1466249004684035 1.3925908906216544; 7.298007833839069 7.946302191563701; 9.753979199837255 10.246020800162745])
 ([-0.0014352253965986877 0.0014352253965986877; 0.0016158790419848648 -0.0016158790419848648; 0.0007253888815863582 -0.00966964744108534], [2.219626292892629 2.4184375476116333; 4.366347626319716 4.142511746204655; 10.359987047863301 9.640012952136699])
 ([-0.001747852894395606 0.001747852894395606; -0.005915681624703305 0.005915681624703305; -0.006844970168551906 -0.0020992883909470767], [9.405847059417546 9.652700169718063; 0.4091258941917879 1.244610627082068; 9.832439230004061 10.167560769995939])
 ([0.0017416290279954887 -0.0017416290279954887; 0.005314573056509905 -0.005314573056509905; -0.005163443000892728 -0.003780815558606254], [8.17610631472349 7.799154010757816; 4.183292241498323 3.0330240006485303; 9.925187196653482 10.074812803346518])
 ([-0.0025295327460284358 0.0025295327460284358; 0.003683115142782889 -0.003683115142782889; -0.009426036859631602 0.00048177830013262056], [-3.4391105961530726 -3.0457374485871456; 9.15424196300173 8.581472703999637; 9.614803520431076 10.385196479568924])
 ([-0.0004515222805755135 0.0004515222805755135; -0.0012319741685254857 0.0012319741685254857; -0.010600118031019823 0.0016558594715208406], [2.3755039354067424 2.2931072604703218; 4.859366721162387 4.634548191798397; 10.559137299391828 9.440862700608172])
 ([0.0027351441185419745 -0.0027351441185419745; 0.0019232288917905376 -0.0019232288917905376; -0.002369826035279446 -0.006574432524219536], [5.931275321341408 6.850973770454511; 2.803315488621059 3.4500055906447; 9.646547869196883 10.353452130803117])
 ([-0.0007176053817506439 0.0007176053817506439; -0.0005215194764194578 0.0005215194764194578; -0.011072017175984467 0.0021277586164854846], [-2.5483294230908147 -2.6608967550634315; 6.635122661480696 6.553314389796504; 10.51764646049972 9.48235353950028])
 ⋮
 ([0.003790046179709312 -0.003790046179709312; -2.5824925546042113e-5 2.5824925546042113e-5; -0.005404330540678482 -0.0035399280188205], [0.5359362468572302 1.1864253142726224; 2.5355595905619515 2.531127235139885; 10.079997274454273 9.920002725545727])
 ([-0.0038271164037234347 0.0038271164037234347; 0.000651399587160879 -0.000651399587160879; -0.004068427939419969 -0.004875830620079013], [2.1656504690438427 2.8225019326687915; 9.200566578586514 9.088766275931746; 10.034643813813034 9.965356186186966])
 ([-0.0020293251034335037 0.0020293251034335037; 0.004507392693807035 -0.004507392693807035; -0.00014969074134335483 -0.008794567818155626], [4.5374998883113715 4.821629592011588; 4.882529351791418 4.2514406472831725; 10.302596459064592 9.697403540935408])
 ([-0.00023543556780646616 0.00023543556780646616; -0.0003040150577992492 0.0003040150577992492; 0.0009612424486640046 -0.009905501008162986], [8.349325262624664 8.295880494673987; 0.9069590351258119 0.8379464625184789; 9.38330199230235 10.61669800769765])
 ([0.00046259560973824324 -0.00046259560973824324; 0.006060182754158149 -0.006060182754158149; -0.004357698154154981 -0.004586560405344001], [1.0711112718208853 1.0080066702435235; 1.6148171499626467 0.7881222783802246; 10.007805014181583 9.992194985818417])
 ([0.0011359332337231855 -0.0011359332337231855; -0.0009914573531539642 0.0009914573531539642; 0.0019007023842411847 -0.010844960943740167], [2.4847430778896973 2.643328102891478; 0.9270628877025695 0.7886478079511166; 9.555151817570477 10.444848182429523])
 ([-0.0015578138811444967 0.0015578138811444967; 0.0007044204117113229 -0.0007044204117113229; -0.009847306398033502 0.0009030478385345193], [4.2560264624511674 4.59274027111004; 7.756794492690342 7.604537485229275; 9.419090951229766 10.580909048770234])
 ([0.0011576545928499201 -0.0011576545928499201; -0.006139119309815578 0.006139119309815578; -0.004918549322504346 -0.004025709236994636], [9.314633129509287 9.103146357791742; 0.5954804918908222 1.7170089098150911; 9.959222753377668 10.040777246622332])
 ([-0.0064472126224431785 0.0064472126224431785; 0.0005884836260338818 -0.0005884836260338818; -0.003598430638308026 -0.005345827921190956], [6.544185547937426 7.4421462747554745; 2.839559526449789 2.757596170267865; 10.060843895884375 9.939156104115625])

The output contains both the positions and velocities, these can be passed directly to the DynamicalDistribution for use with dynamics. Here however, let's focus on the positions and visualise the distribution.

This collects the x and y coordinate for each atom:

r = last.(configurations)
x = hcat([i[1,:] for i in r]...)
y = hcat([i[2,:] for i in r]...)
2×150 Matrix{Float64}:
 6.33587  6.05192  7.29801  4.36635  …  0.927063  7.75679  0.59548  2.83956
 5.43538  5.33123  7.9463   4.14251     0.788648  7.60454  1.71701  2.7576

Now we can plot the distribution:

using Plots

plt = plot(
    xlabel="x coordinate / bohr",
    ylabel="y coordinate / bohr",
    legend=false,
)

for i=1:nsamples
    plot!([x[1,i], x[2,i]], [y[1,i], y[2,i]], linewidth=5, color=:black)
    scatter!([x[1,i], x[2,i]], [y[1,i], y[2,i]])
end
plt
Example block output

Here we can see that the molecule is randomly distributed within the unit cell. Since we have used a harmonic potential, this could have been produced without using the EBK procedure, but this technique can use any arbitrary potential. In the hydrogen scattering example we build on this example and use the sample procedure to perform scattering simulations starting from this distribution.