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.003803695354146845 0.003803695354146845; -5.1752933043987364e-5 5.1752933043987364e-5; -0.004970546651794926 -0.003973711907704056], [5.643989841202298 6.3053129388466385; 8.23659120073855 8.245589137430697; 9.95667175079604 10.04332824920396])
 ([-0.0028378323649324005 0.0028378323649324005; -0.004223599371026677 0.004223599371026677; -0.008817839433323786 -0.00012641912617519617], [2.631011176115023 3.046592644175072; 1.9842586729078677 2.6027764010156647; 9.681799984434795 10.318200015565205])
 ([-0.0008820177926860737 0.0008820177926860737; -0.005439646614311254 0.005439646614311254; -0.002695286191161346 -0.006248972368337636], [-0.6489830801336319 -0.5277777347858713; 3.2006289891728605 3.948135731014784; 10.122085337714815 9.877914662285185])
 ([-0.003717083913330497 0.003717083913330497; 0.002900911535159273 -0.002900911535159273; -0.003813225563634484 -0.005131032995864497], [-0.43273343866135444 -0.9903769947806287; 6.371373304245582 6.806573224982428; 9.950575044852114 10.049424955147886])
 ([0.0031344382278381882 -0.0031344382278381882; 0.0013189304402121682 -0.0013189304402121682; 0.0012938088953786953 -0.010238067454877677], [-0.09981448840947252 0.37627815709729906; 7.197943053858677 7.398276586337672; 9.562103228699044 10.437896771300956])
 ([0.00040824617809901006 -0.00040824617809901006; 0.0019525749716809353 -0.0019525749716809353; -0.0073494067669173085 -0.001594851792581674], [-2.3625752766084918 -2.5195567872700972; 7.398640621114229 6.647823625933013; 9.44680514998148 10.55319485001852])
 ([0.0015269374556054305 -0.0015269374556054305; -0.0001650203751682248 0.0001650203751682248; -0.007618423923068603 -0.0013258346364303781], [2.733670376804001 3.3208184294395235; 2.3692439422734464 2.305789220163176; 10.604916974844238 9.395083025155762])
 ([0.004166277451446088 -0.004166277451446088; 0.0026184084350370828 -0.0026184084350370828; -0.008922823603878118 -2.1434955620863554e-5], [-1.425733329005312 -0.7624177905952276; 6.749577725541579 7.166456107047662; 10.354298859871106 9.645701140128894])
 ([0.0018674210558226934 -0.0018674210558226934; -0.00021452300699438472 0.00021452300699438472; -0.002432041345597661 -0.006512217213901321], [6.391778737881378 7.313832205486662; 1.57521961894781 1.4692972390627106; 9.496345468516823 10.503654531483177])
 ([0.0010029526203614012 -0.0010029526203614012; -0.00196419160192784 0.00196419160192784; -0.006154769375420007 -0.0027894891840789756], [3.817421233410004 3.32220568544103; 5.0638000497612685 6.033634716641399; 9.584591774279534 10.415408225720466])
 ⋮
 ([-0.00031159024571322043 0.00031159024571322043; -0.001782228953076496 0.001782228953076496; -0.0010823518471349209 -0.007861906712364062], [6.31086739492197 6.202793957254822; 8.330924500680572 7.712767763366737; 9.412136764373631 10.587863235626369])
 ([-0.0012688211102988521 0.0012688211102988521; -0.002830901699204945 0.002830901699204945; -0.006739309935226539 -0.002204948624272443], [6.951703095335747 7.391787023061944; 2.368617108429848 3.3505004714040374; 9.60682015785023 10.39317984214977])
 ([0.0015989984301171063 -0.0015989984301171063; -0.006168967243607155 0.006168967243607155; -0.002419655968482741 -0.00652460259101624], [4.542833994966838 4.3004241703347965; 0.595482928639753 1.530704777215361; 10.155578544067012 9.844421455932988])
 ([-0.006085675885408606 0.006085675885408606; -0.00011121085993219125 0.00011121085993219125; -0.0020043594768794615 -0.00693989908261952], [4.100603253187802 3.248315942888208; 8.214208087346869 8.198633218347524; 9.827196770964296 10.172803229035704])
 ([0.006255127793536046 -0.006255127793536046; 0.0013253903592505202 -0.0013253903592505202; -0.0029706555637935687 -0.005973602995705414], [2.842934180944996 3.7189529080539114; 5.4465510354470155 5.632169429468141; 9.894860733716868 10.105139266283132])
 ([-0.005056835800998405 0.005056835800998405; -0.001192125782774056 0.001192125782774056; -0.008122740871438928 -0.000821517688060055], [3.4098803801066966 2.7166736299232808; 3.5053302624719893 3.341909959300685; 10.25021858502478 9.74978141497522])
 ([-0.0033919133406433025 0.0033919133406433025; -0.0003518177677732181 0.0003518177677732181; -0.004244823786784387 -0.004699434772714595], [-3.2097718800291015 -2.5640066467753457; 6.418255022520928 6.485235413837705; 10.021637637808375 9.978362362191625])
 ([0.004439807056430802 -0.004439807056430802; 0.0023517496961292634 -0.0023517496961292634; -0.0037581532042903027 -0.005186105355208679], [6.199550926158379 5.556657217069198; 8.659595521403975 8.3190570552551; 10.051692643567478 9.948307356432522])
 ([0.0014951104816556457 -0.0014951104816556457; 0.004372854971362113 -0.004372854971362113; -0.001432822374438058 -0.007511436185060924], [1.8592186715719667 1.5269678693283377; 9.466705457880316 8.494948126127; 10.337704861929698 9.662295138070302])

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}:
 8.23659  1.98426  3.20063  6.37137  …  3.50533  6.41826  8.6596   9.46671
 8.24559  2.60278  3.94814  6.80657     3.34191  6.48524  8.31906  8.49495

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.