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.004214295690324429 0.004214295690324429; -0.004938195733055206 0.004938195733055206; -0.0033313265344328147 -0.005612932025066167], [2.6787019672236414 3.2715585786609656; 6.491647899773637 7.186340964532931; 10.080242643089761 9.919757356910239])
 ([0.005968623287510229 -0.005968623287510229; 0.001163055594040892 -0.001163055594040892; -0.004365528313156955 -0.004578730246342027], [-1.4985526082080902 -0.349966322474935; 6.153773031840323 6.377588414469011; 9.989742994123091 10.010257005876909])
 ([0.004654408567604114 -0.004654408567604114; -0.004394465533998505 0.004394465533998505; -0.006435292567863545 -0.0025089659916354375], [4.617593138421673 3.9122867740327067; 6.949754531101049 7.615670396500154; 9.851255907025545 10.148744092974455])
 ([-0.0010722682272641083 0.0010722682272641083; 0.0004904831583424754 -0.0004904831583424754; 0.0015385140976278723 -0.010482772657126854], [1.2766224019674808 1.0727568882613023; 5.4051829812985215 5.498436328266424; 9.428611764912503 10.571388235087497])
 ([0.0010249481276080222 -0.0010249481276080222; -0.005589400732545496 0.005589400732545496; -0.005883240485143305 -0.0030610180743556765], [3.6889547800772164 3.479349731219212; 7.92053001345278 9.063579678186855; 9.855711706190966 10.144288293809034])
 ([0.0016643189787050343 -0.0016643189787050343; 0.0034556998201688074 -0.0034556998201688074; -0.007506124825090841 -0.0014381337344081414], [2.3374577949335036 2.772247627957051; 2.1303986995825492 3.03317224661003; 10.396303362946552 9.603696637053448])
 ([-0.0010674200973184386 0.0010674200973184386; -0.00012453752089082164 0.00012453752089082164; -0.005296728480262817 -0.0036475300792361645], [4.677462366647145 5.778917473524581; 7.419804600125833 7.5483130371766745; 9.574554103480827 10.425445896519173])
 ([0.0061160305172225544 -0.0061160305172225544; -0.000435134807466096 0.000435134807466096; -0.0027601722903283464 -0.0061840862691706355], [6.633657087666473 5.794545918966341; 7.760889457901646 7.820589368784496; 10.11743909927463 9.88256090072537])
 ([0.004746687940464727 -0.004746687940464727; 0.003249083062768142 -0.003249083062768142; -0.007748514956949082 -0.0011957436025498997], [-3.8638409474552353 -3.101953118898964; 7.9782434148224235 8.499751647623803; 10.262945278942233 9.737054721057767])
 ([0.002609733910104062 -0.002609733910104062; 0.0039020192372587084 -0.0039020192372587084; -0.002374413445525468 -0.006569845113973514], [-1.8910147370671861 -1.515614545438642; 6.997823528817108 7.559113997448362; 9.849125820241472 10.150874179758528])
 ⋮
 ([-0.0022980546130067824 0.0022980546130067824; 0.0019086739578502013 -0.0019086739578502013; -0.0015559955657553654 -0.0073882629937436165], [2.8181112414124145 3.1937205512927953; 2.1167179907488523 1.8047516414562417; 10.238316131747434 9.761683868252566])
 ([-0.004463979982608644 0.004463979982608644; 0.0032334919888441136 -0.0032334919888441136; -0.004224669478545723 -0.004719589080953259], [-1.7890033986510558 -2.4116566773683457; 7.30441986143277 7.755439892078904; 9.982741672546634 10.017258327453366])
 ([0.0018934919863588925 -0.0018934919863588925; -0.0016852444428580285 0.0016852444428580285; -0.007396944661900035 -0.0015473138975989472], [5.156294146706142 4.504880352806409; 4.876855555005165 5.456626414390526; 9.496891167676202 10.503108832323798])
 ([0.001446389032504002 -0.001446389032504002; -0.0032793420891516424 0.0032793420891516424; -0.0030113829009334007 -0.005932875658565581], [2.689792175023385 3.1873901442285275; 2.227743928699206 1.0995592375681968; 9.748731352600103 10.251268647399897])
 ([-0.0014224761643337968 0.0014224761643337968; 0.0056566808497787125 -0.0056566808497787125; -0.007351057621202589 -0.0015932009382963932], [-0.11786524319257301 -0.3152710759715423; 1.9065563315993743 2.691568989044072; 10.19976336370518 9.80023663629482])
 ([-0.001825618988283217 0.001825618988283217; 0.004361468907617757 -0.004361468907617757; -0.008905475839527957 -3.8782719971024e-5], [7.369216752320113 7.116512519719344; 6.100556462763236 6.7042758918365; 10.306834407242775 9.693165592757225])
 ([-0.0010459926767460068 0.0010459926767460068; 0.0005813337593278052 -0.0005813337593278052; -0.003397575974496754 -0.005546682585002228], [-0.28212949925346575 0.11181225806648529; 5.892277521401052 5.673335603517415; 10.20234913055133 9.79765086944867])
 ([-0.002792495844415072 0.002792495844415072; 0.003301782549646592 -0.003301782549646592; 0.0005774510935660357 -0.009521709653065018], [1.4785419942212923 1.0373694755123122; 2.6796802424539914 3.2013124327293045; 9.601120965645398 10.398879034354602])
 ([-0.005097746132333547 0.005097746132333547; 0.0022752003030002367 -0.0022752003030002367; -0.005193552877910761 -0.0037507056815882208], [5.335285996581638 4.2301527019167215; 6.81680239180073 7.310039901399168; 10.078198209285878 9.921801790714122])

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.49165  6.15377  6.94975  5.40518  …  6.10056  5.89228  2.67968  6.8168
 7.18634  6.37759  7.61567  5.49844     6.70428  5.67334  3.20131  7.31004

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

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.