Elastic scattering¶
This example illustrates how the QDYN library can be used for scattering calculations using elastic scattering of Ne with the HD molecules.
In order to run the example, QDYN need to be installed. If everything is set up correctly, the following commands should succeed and display the installed QDYN version.
[1]:
! qdyn_version
QDYN 23.04+dev revision 266eb1a578d191e9e85a582473afcff1423ed041 (scattering)
features: no-check-cheby, no-check-newton, no-parallel-ham, no-parallel-oct, backtraces, debug, no-no-ipo, no-shared-lib
compiled with gfortran on Fri May 26 16:53:27 2023 on host romema
Additionally, the following python packages will be used in this example:
[2]:
from pathlib import Path
import matplotlib.pyplot as plt
import numpy as np
In order to clean up the folder of this example, run
[3]:
! make clean
Physical background¶
For the collision with the HD molecule, only the isotropic part of the interaction potential is taken into account. Hence, the HD molecule effectively behaves like an atom and the scattering problem can be formulated in terms of a one-dimensional potential energy surface with the distance between Ne and HD being the reaction coordinate \(r\). As a consequence, the collision has cylindrical symmetry.
We will solve the scattering problem for partial waves with quantum number \(l\) up to \(l_{\mathrm{max}}\). Due to the cylindrical symmetry, we only need to consider partial waves with \(m_l=0\). Without external fields, spin-orbit or hyperfine interactions, the individual partial waves are uncoupled and the collision is purely elastic. However, the procedure followed in this example in principle also apply for inelastic collisions.
Setting up the simulation¶
We will solve the scattering problem with the renormalized Numerov method. The directory of this example already contains a runfolder with a config file set up for \(l_{\mathrm{max}} = 8\):
[4]:
! cat ./rf_numerov/config
grid: dim = 1, coord_type = cartesian
* r_min = 3.5, r_max = 20000.00, nr = 1000, &
maptype = doubling, E_max = 2_K, npt = 20
ham: mass = 4791.3134, type= op
* op_type = pot, op_form = file, op_surf = 1, filename = NeHD_isotropic
scattering: system = atom-atom, method = numerov, l_max = 8, &
E_min = 2.0d-3_K, E_max = 2_K, nE = 101, E_grid_type = log
The grid section defines the spatial grid for the reaction coordinate. The grid needs to be large enough, such that all potentials are decayed to zero in the asymptotic limit. Since this can lead to very large grids, a mapped grid is used here with maptype = doubling. This means, that the step size of the grid will be adjusted according to the local de Broglie wavelength, with npt points per period for a particle with energy E_max. The parameter E_max should be set to the
highest scattering energy of interest. If not given, E_max from the scattering section is used.
The ham section contains the isotropic potential as a single potential energy surface. Since, system is set to atom-atom, the partial waves will be initialized automatically as individual channels based on this potential.
We use the renormalized Numerov method to solve the scattering problem, selected with method = numerov, in the scattering section. Lastly, the set of collision energies is specified as a grid from E_min to E_max with nE steps equally distributed in log-space. Note that the collision energy is defined with respect to the asymptotic limit of the entrance channel potential.
Since we only have a single potential energy surface, we don’t need to specify the index of the entrance channel. In a more complex scattering situation, this can be done by setting entrance_surf or entrance_channel in the scattering section. See the documentation of the scattering config file section for details.
The calculation can be performed with the qdyn_scattering utility program:
[5]:
! qdyn_scattering rf_numerov
QDYN 23.04+dev revision 266eb1a578d191e9e85a582473afcff1423ed041 (scattering)
features: no-check-cheby, no-check-newton, no-parallel-ham, no-parallel-oct, backtraces, debug, no-no-ipo, no-shared-lib
compiled with gfortran on Fri May 26 16:53:27 2023 on host romema
***** START OF PROGRAM qdyn_scattering ******
Fri May 26 16:57:37 +0200 2023
*** Read config file rf_numerov/config ***
*** Done reading config file ***
*** Initializing system ***
*** Initializing grid ***
Initializing grid for label `` as 1D cartesian grid with mapping
Original grid parameters were changed by mapping
original after mapping
r_min 3.50000000000000000E+00 3.50000000000000000E+00
r_max 2.00000000000000000E+04 7.71852720075120033E+02
nr 1000 1000
*** Initializing pulses ***
*** Initializing dynamical generator ***
*** Done with initialization ***
*** Initializing scattering calculation ***
*** Initializing channels ***
Initializing coupled channels for atom-atom scattering
Projection of total angular momentum is conserved: M = 0.0
Channels coupled with entrance channel by present couplings will be initialized.
List of channels:
# 1: l= 0, ml= 0, i1= 0.0, mi1= 0.0, s1= 0.0, ms1= 0.0, i2= 0.0, mi2= 0.0, s2= 0.0, ms2= 0.0
# 2: l= 1, ml= 0, i1= 0.0, mi1= 0.0, s1= 0.0, ms1= 0.0, i2= 0.0, mi2= 0.0, s2= 0.0, ms2= 0.0
# 3: l= 2, ml= 0, i1= 0.0, mi1= 0.0, s1= 0.0, ms1= 0.0, i2= 0.0, mi2= 0.0, s2= 0.0, ms2= 0.0
# 4: l= 3, ml= 0, i1= 0.0, mi1= 0.0, s1= 0.0, ms1= 0.0, i2= 0.0, mi2= 0.0, s2= 0.0, ms2= 0.0
# 5: l= 4, ml= 0, i1= 0.0, mi1= 0.0, s1= 0.0, ms1= 0.0, i2= 0.0, mi2= 0.0, s2= 0.0, ms2= 0.0
# 6: l= 5, ml= 0, i1= 0.0, mi1= 0.0, s1= 0.0, ms1= 0.0, i2= 0.0, mi2= 0.0, s2= 0.0, ms2= 0.0
# 7: l= 6, ml= 0, i1= 0.0, mi1= 0.0, s1= 0.0, ms1= 0.0, i2= 0.0, mi2= 0.0, s2= 0.0, ms2= 0.0
# 8: l= 7, ml= 0, i1= 0.0, mi1= 0.0, s1= 0.0, ms1= 0.0, i2= 0.0, mi2= 0.0, s2= 0.0, ms2= 0.0
# 9: l= 8, ml= 0, i1= 0.0, mi1= 0.0, s1= 0.0, ms1= 0.0, i2= 0.0, mi2= 0.0, s2= 0.0, ms2= 0.0
Initialising wigner symbols on 16 threads:
*** Using renormalized Numerov method ***
Collision energy: 6.33362431804951774E-09 [iu]
Collision energy: 6.78660115805814544E-09 [iu]
Collision energy: 7.27197458695294708E-09 [iu]
Collision energy: 7.79206158626306716E-09 [iu]
Collision energy: 8.34934484642695936E-09 [iu]
Collision energy: 8.94648461817835619E-09 [iu]
Collision energy: 9.58633141153512004E-09 [iu]
Collision energy: 1.02719396030112880E-08 [iu]
Collision energy: 1.10065820160070447E-08 [iu]
Collision energy: 1.17937655439783122E-08 [iu]
Collision energy: 1.26372478909645975E-08 [iu]
Collision energy: 1.35410555093875286E-08 [iu]
Collision energy: 1.45095028207485246E-08 [iu]
Collision energy: 1.55472128109772577E-08 [iu]
Collision energy: 1.66591390987453028E-08 [iu]
Collision energy: 1.78505895820899486E-08 [iu]
Collision energy: 1.91272517762284147E-08 [iu]
Collision energy: 2.04952199635147439E-08 [iu]
Collision energy: 2.19610242851429793E-08 [iu]
Collision energy: 2.35316619134683483E-08 [iu]
Collision energy: 2.52146304537515420E-08 [iu]
Collision energy: 2.70179637347725859E-08 [iu]
Collision energy: 2.89502701591641127E-08 [iu]
Collision energy: 3.10207737965342201E-08 [iu]
Collision energy: 3.32393584155410420E-08 [iu]
Collision energy: 3.56166146651100728E-08 [iu]
Collision energy: 3.81638906300200372E-08 [iu]
Collision energy: 4.08933460021887319E-08 [iu]
Collision energy: 4.38180101262520795E-08 [iu]
Collision energy: 4.69518441965217808E-08 [iu]
Collision energy: 5.03098079022271227E-08 [iu]
Collision energy: 5.39079308391774971E-08 [iu]
Collision energy: 5.77633890287381677E-08 [iu]
Collision energy: 6.18945869093888722E-08 [iu]
Collision energy: 6.63212451922626394E-08 [iu]
Collision energy: 7.10644950000519760E-08 [iu]
Collision energy: 7.61469787386628239E-08 [iu]
Collision energy: 8.15929581831403830E-08 [iu]
Collision energy: 8.74284302938254378E-08 [iu]
Collision energy: 9.36812513156001095E-08 [iu]
Collision energy: 1.00381269752627221E-07 [iu]
Collision energy: 1.07560468853349664E-07 [iu]
Collision energy: 1.15253119285917998E-07 [iu]
Collision energy: 1.23495942732859614E-07 [iu]
Collision energy: 1.32328287185922754E-07 [iu]
Collision energy: 1.41792314777887931E-07 [iu]
Collision energy: 1.51933203047980297E-07 [iu]
Collision energy: 1.62799360601641848E-07 [iu]
Collision energy: 1.74442658194143312E-07 [iu]
Collision energy: 1.86918676341135502E-07 [iu]
Collision energy: 2.00286970638132384E-07 [iu]
Collision energy: 2.14611356055462350E-07 [iu]
Collision energy: 2.29960211565795457E-07 [iu]
Collision energy: 2.46406806558419495E-07 [iu]
Collision energy: 2.64029650598440811E-07 [iu]
Collision energy: 2.82912868200515156E-07 [iu]
Collision energy: 3.03146600406139817E-07 [iu]
Collision energy: 3.24827435081475137E-07 [iu]
Collision energy: 3.48058867989761891E-07 [iu]
Collision energy: 3.72951796839322583E-07 [iu]
Collision energy: 3.99625050665533827E-07 [iu]
Collision energy: 4.28205957073822037E-07 [iu]
Collision energy: 4.58830950051491177E-07 [iu]
Collision energy: 4.91646221249830715E-07 [iu]
Collision energy: 5.26808417845459179E-07 [iu]
Collision energy: 5.64485390312235085E-07 [iu]
Collision energy: 6.04856993673291325E-07 [iu]
Collision energy: 6.48115946058057247E-07 [iu]
Collision energy: 6.94468748662689576E-07 [iu]
Collision energy: 7.44136671505421210E-07 [iu]
Collision energy: 7.97356809682458475E-07 [iu]
Collision energy: 8.54383215166568414E-07 [iu]
Collision energy: 9.15488109551100093E-07 [iu]
Collision energy: 9.80963183528628785E-07 [iu]
Collision energy: 1.05112098930741324E-06 [iu]
Collision energy: 1.12629643261249132E-06 [iu]
Collision energy: 1.20684837139367514E-06 [iu]
Collision energy: 1.29316132887203680E-06 [iu]
Collision energy: 1.38564732910224849E-06 [iu]
Collision energy: 1.48474786381313282E-06 [iu]
Collision energy: 1.59093599991521499E-06 [iu]
Collision energy: 1.70471863773579096E-06 [iu]
Collision energy: 1.82663893076140044E-06 [iu]
Collision energy: 1.95727887843856710E-06 [iu]
Collision energy: 2.09726210440987753E-06 [iu]
Collision energy: 2.24725683344757948E-06 [iu]
Collision energy: 2.40797908129536070E-06 [iu]
Collision energy: 2.58019607264543659E-06 [iu]
Collision energy: 2.76472990356698620E-06 [iu]
Collision energy: 2.96246146586887908E-06 [iu]
Collision energy: 3.17433465213014339E-06 [iu]
Collision energy: 3.40136086147130476E-06 [iu]
Collision energy: 3.64462382757525940E-06 [iu]
Collision energy: 3.90528479200501571E-06 [iu]
Collision energy: 4.18458804751335287E-06 [iu]
Collision energy: 4.48386687780623855E-06 [iu]
Collision energy: 4.80454992211400465E-06 [iu]
Collision energy: 5.14816799495214365E-06 [iu]
Collision energy: 5.51636139362673943E-06 [iu]
Collision energy: 5.91088772836762752E-06 [iu]
Collision energy: 6.33363031246712307E-06 [iu]
*** Done Numerov propagation ***
*** Calculating scattering observables ***
**** WARNING **** in thread 001 module qdyn_scattering_mod, routine init_cross_sections:
Total cross sections are calculated for the entrance channel only and with certain assumptions about the channel quantum numbers. Check, if these really are the quantities that you would consider to be total cross sections.
*** Writing scattering observables ***
Fri May 26 16:57:38 +0200 2023
***** END OF PROGRAM qdyn_scattering ******
After the calculation has finished, the scattering observables are stored inside the runfolder. Below, a plot of the elastic cross section is shown, displaying a shape resonance in the partial wave \(l=1\).
[6]:
runfolder = Path('./rf_numerov')
E_iu, sigma = np.loadtxt(runfolder/'total_elastic_cross_section.dat', unpack=True)
sigma_ch = np.loadtxt(runfolder/f'elastic_cross_sections.dat', unpack=True)[1:]
E = E_iu * 3.157746662555312504e5 # conversion from atomic units to K
[7]:
fig, ax = plt.subplots()
ax.plot(E, sigma, 'k-', label='total')
for l in range(len(sigma_ch)):
ax.plot(E, sigma_ch[l], label=f'l={l}')
ax.set_xscale('log')
ax.set_xlabel('Collision energy [K]')
ax.set_ylabel('Cross section (arb. units)')
ax.legend()
plt.show()
