Asymmetric top molecule: Excitation of rotational states with microwave fields¶
This example illustrates the use of the QDYN library to simulate the rotational dynamics of symmetric or asymmetric top molecules.
In order to run the example, QDYN and the QDYN-pylib need to be installed. If everything is set up correctly, the following commands should succeed and display the installed QDYN version.
[1]:
import qdyn
[2]:
! qdyn_version
QDYN 23.04+dev revision 1b80bbb91a2b8be8e979590d43fc56ff1184a868 (package1)
features: no-check-cheby, no-check-newton, no-parallel-ham, use-mkl=sequential, parallel-oct, backtraces, no-debug, no-no-ipo, no-shared-lib
compiled with ifort on Fri Jan 19 11:25:40 2024 on host balmora
Additionally, the following python packages will be used in this example:
[3]:
from pathlib import Path
from subprocess import run
from multiprocessing import Pool
import matplotlib.pyplot as plt
import numpy as np
In order to clean up the folder of this example, run
[4]:
! make clean
Physical background¶
The eigenfunctions of an asymmetic top are determined by \(H_{rot} |j,\tau,m\rangle = E_{j,\tau} |j,\tau,m \rangle\) with \(H_{rot} = A J^2_a + B J^2_b + C J^2_c\), where \(J_a\), \(J_b\) and \(J_c\) are the angular momentum operators with respect to the principle molecular axis and \(A > B > C\) are the rotational constants. The energy eigenvalues \(E_{j,\tau}\) depend on the quantum numbers \(j=0,1,2,...\) and \(\tau=0,1,..,(2j+1)\). The eigenvalues are \((2 j+1)\)-fold degenerate. The eigenvalues and eigenstates of an asymmetric top can be written as a superposition of symmetric top eigenstates
\begin{equation} |j,\tau,m \rangle = \sum_k c_k^j(\tau) |j,k,m\rangle. \end{equation}
Note that only states with same \(j\) and \(m\) mix. The eigenvalues are thus calculated by diagonalizing the rotational Hamiltonian in the symmetric top basis for every value of \(j\) and \(m\).
The interaction of a rigid asymmetric top with an electromagnetic field in dipole approximation is described by \(H_{int} = \vec{\mu} \cdot \vec{E}(t)\) with the electric field
\begin{equation} \vec{E}(t) = \vec{e} E(t) \cos(\omega_L t + \phi), \end{equation}
where \(\vec{e}\) denotes the polarization direction, \(E(t)\) the temporal shape, \(\omega_L\) the frequency and \(\phi\) the phase of the field. Note, that \(\vec{\mu}^T = (\mu_x,\mu_y,\mu_z)\) is the molecular dipole moment in space fixed coordinates, The transformation of the interaction Hamiltonian into molecule fixed coordinates can be expressed in terms of Wigner-D-matrix-elements (see Ref. [1]). In the symmetric top basis, the interaction Hamiltonian, i.e. the transition matrx elements can then be expressed by Wigner-3j-symbols (see Ref. [1]). The eigenvalues of the asymmetric top consitute the transformation matrix \(U\) between the symmetric and asymmetric top basis. Here, the propagation is carried out in the asymmetric top basis.
All energy eigenvalues and transition matrix elements up to \(j=j_{max}\) are calculated. Interaction with microwave radiation typically affects only a few rotational levels. The relevant rotational levels can be selected from the complete Hamiltonian and the propagation can be performed within this subsystem.
Setting up the simulation¶
Identifying relevant rotational levels¶
As a first step, we will determine which part of the rotational spectrum is relevant for our simulation.
In this example, we consider the asymmetric top molecule carvone. The runfolder rf_spectrum contains a config file that defines the field free rotational Hamiltonian, \(H_{rot}\), up to a given maximum rotational quantum number \(j_{max}\).
[5]:
! cat ./rf_spectrum/config
ham: A = 3.4089, B = 1, C = 0.8832, jmax = 5
* type = matrix, op_type = pot, op_form = top, sparsity_model = indexed
In qdyn, the field free Hamiltonian of a rigid rotor is stored as type = matrix with op_form = top and op_type = pot. As molecular parameters, we have to specify the rotational constants A, B, C and the maximum rotational quantum number jmax. Qdyn will use this value to determine the dimension of the Hamiltonian matrix automatically. Addtionally, we need to specify how the matrix shall be stored by choosing a sparsity_model.
In this example, we use carvone with the rotational constants \(A=2237.21\) MHz, \(B=656.28\) MHz and \(C=579.64\) MHz [2]. Here, we use scaled variables, such that the rotational constants \(A,\) \(B\) and \(C\) are given in units of \(B\). The interaction Hamiltonian \(H_{int}\), which we will use later on for the propagation, depends in addition on the components \(\mu_a\), \(\mu_b\) and \(\mu_c\) of the molecular dipole moments. In the following we will extract the molecular parameters from the config file. The values of the dipole moment components (see Ref. [2]) are given in Debye.
[6]:
config_spectrum = qdyn.config.read_config_file('./rf_spectrum/config')
A = qdyn.config.get_config_value(config_spectrum, ('ham', 0, 'A'))
B = qdyn.config.get_config_value(config_spectrum, ('ham', 0, 'B'))
C = qdyn.config.get_config_value(config_spectrum, ('ham', 0, 'C'))
jmax = qdyn.config.get_config_value(config_spectrum, ('ham', 0, 'jmax'))
mu_a = 2.0
mu_b = 3.0
mu_c = 0.5
We can construct the Hamiltonian matrix using the `write_ham utility <../../doc/utils/write_ham.rst>`__:
[7]:
! write_ham rf_spectrum
*** Read config file rf_spectrum/config ***
*** Done reading config file ***
*** Initializing system ***
*** Initializing grid ***
No explicit grid found in config file! Initializing grid as one dimensional 1x1 fake grid
Initializing grid as 1D cartesian grid without mapping
Number of grid points was set to nr = 1 for dimension dim = 1! Initializing dimension 1 as a 1x1 fake-grid...
*** Initializing pulses ***
*** Initializing dynamical generator ***
*** Done with initialization ***
Written op_matrix 1 to rf_spectrum/op_matrix1.dat
Now, \(H_{rot}\) is saved in op_matrix1.dat inside the runfolder. Since \(H_{rot}\) is diagonal and stored in a sparse matrix format, the file only contains the rotational energies. By default, rotational levels are sorted by quantum number \(j\) (in increasing order) and by \(m = -j,-j+1,...,j\).
We can load the energies into a NumPy array and visualize it using the matplotlib package.
[8]:
row, col, Erot = np.loadtxt('./rf_spectrum/op_matrix1.dat',
unpack=True)
fig, ax = plt.subplots()
ax.plot(col, Erot, '.')
ax.set_xlabel('Level index')
ax.set_ylabel('Rotational energy (in units of B)')
plt.show()
For each rotational quantum number \(j\), there are \((2 j +1)\) levels with different energy, denoted by \(\tau=1,...,(2j+1)\). Each of these states is \((2j+1)\)-fold degenerate.
Note that the labeling \((j,\tau)\) can be related to the spectroscopic notation \(j_{K_a,K_c}\) with \(K_a=|k_a|\) and \(K_b=|k_b|\) as follows:
\begin{align*} j_{j,0} &\Leftrightarrow (j,2j+1) \\ &... \\ j_{1,j-1} &\Leftrightarrow (j,3) \\ j_{1,j} &\Leftrightarrow (j,2) \\ j_{0,j} &\Leftrightarrow (j,1) \\ \end{align*}
For the interaction with resonant microwaves, only those states which are in resonance (or in near resonance) with the frequency of the field are driven by the interaction. In this example, transition between the rotational levels \(4_{0,4}=(4,1)\) and \(5_{1,5}=(5,2)\) are driven by a single microwave field with frequency \(\omega_L=E_{5,2}-E_{4,1}\). We therefore only select the rotational states of these two states and the corresponding transition matrix elements.
[9]:
levels = [
[4, 1], # j, tau
[5, 2],
]
n_level = np.sum(2*np.transpose(levels)[0]+1) # = sum_j 2j+1
print('Number of levels in subspace:', n_level)
Number of levels in subspace: 20
Below, we define some utility functions that help us to find the selected rotational levels in the larger rotational space.
[10]:
def idx(j, tau, m=0):
"""Index of the given level in the total rotational space."""
nj = j * (2*j-1) * (2*j+1) // 3
return nj + (m+j) * (2*j+1) + tau - 1
def get_subspace_energies(Erot, levels):
"""
Extract energies of levels in the subspace from the list of
all rotational energies.
"""
Erot_sub = np.zeros(np.sum(2*np.transpose(levels)[0]+1))
i = 0
for j, tau in levels:
for m in range(-j, j+1):
Erot_sub[i] = Erot[idx(j, tau, m)]
i += 1
return Erot_sub
Setting up the propagation¶
Now, we need to prepare the runfolder for the propagation, which we call ./rf_prop. We can use the model module of the QDYN-python package, to define the parameters for the propagation and to create the config file.
[11]:
rf_prop = Path('./rf_prop')
rf_prop.mkdir(exist_ok=True)
The following parameters define the initial and final time of the propagation, as well as the number of time steps.
Note that the time is given in units of \(1/(2\pi B)=\) for \(B\) given as a frequency. In the following, this unit is denoted as internal unit and abbreviated with iu. For carvone, the conversion factor is given as iu \(=0.24251\) ns.
[12]:
time_to_ns = 0.24251
t_start = 0
t_stop = 1100
nt = 11000 #50000
As a first step, we will define the initial state that we want to propagate: We start with propagating a single initial state (for propagation of an incoherent mixture of initial states see below). As an example we choose the state with \(m=0\) of the level \(4_{0,4}=(4,1)\).
[13]:
psi = np.zeros(n_level, dtype=np.complex128)
psi[4] = 1
Now we need to define the pulses used for the propagation. Here, we use a Gaussian shaped pulse, linerarly polarised along the \(z\)-direction.
We choose the frequency w_L such that it is resonant to the transition between the two levels, i.e. \(\omega=E_{5,2}-E_{4,1}\).
The maximum amplitude of the electric field is \(E_0\), given in units of \(Bh/\mu_0\), with \(\mu_0=1D\). In the following, this unit is also denoted as internal unit and abbreviated with iu. For carvone, the electric field in units of V/m is \(E_0 (\mbox{V/m}) = 1.3037 \cdot 10^5 E_0 (\mbox{iu})\).
The next cell defines all parameters for our Gaussian pulse.
[14]:
w_L = Erot[idx(5,2)] - Erot[idx(4,1)]
t_fwhm = 400
t_0 = 500
E_0 = 0.02
pulse_pol = 'z'
For technical reasons, QDYN stores the pulses on the intervals of the time grid. That means, the time grid for the pulse needs to be shifted by \(dt/2\) with respect to the time grid used for the propagation. We can let pulse_tgrid take care of the details.
[15]:
pgrid = qdyn.pulse.pulse_tgrid(T=t_stop, nt=nt, t0=t_start)
Next, we calculate the pulse amplitude. You can try different pulses, by adjusting the code below. Instead of a Gaussian, there are other pulse shapes predefined, for example flattop, which has a more or less rectangular shape. See the pulse module for available options. Of course, you can also specify arbitray pulse shapes by yourself.
[16]:
def gaussian(tgrid, t_0, t_fwhm, E_0, w_L, phi=0.0):
sigma = t_fwhm/(2*np.sqrt(2*np.log(2)))
shape = E_0 * qdyn.pulse.gaussian(tgrid, t_0, sigma)
#shape = E_0 * qdyn.pulse.flattop(tgrid, 100, 1000, 100)
return shape * np.cos(w_L * (tgrid-t_0) + phi)
my_amplitude = gaussian(pgrid, t_0, t_fwhm, E_0, w_L)
After caluclating the pulse amplitude, we store it as a Pulse object, which comes shipped with some handy features.
[17]:
my_pulse = qdyn.pulse.Pulse(pgrid,
my_amplitude,
time_unit='iu', ampl_unit='iu', freq_unit='iu'
)
[18]:
my_pulse.show()
Now, we define our Hamiltonian. It consists of one operator, which represents \(H_{rot}\), and additional operators for the dipole interaction. QDYN internally stores the dipole interaction in the form \begin{equation}
H_{int}(t) = \sum_i H_i \, \epsilon_i(t) \;,
\end{equation} where \(H_i\) are the time-independent parts of \(H_{int}(t)\), which are determined by the dipole moment components, and \(\epsilon_i(t)\) are the cooresponding pulse amplitudes. The index \(i\) is called the pulse_id and keeps track about which operator couples to which pulses. In our case, we only have a single laser pulse and hence only one operator that couples to a pulse.
If you want to include more pulses, you also need to include more operators with the corresponding pulses below.
[19]:
model = qdyn.model.LevelModel()
model.add_ham('matrix', op_type='pot', op_form='top', sparsity_model='indexed',
A=A, B=B, C=C, jmax=jmax, mu_a=mu_a, mu_b=mu_b, mu_c=mu_c,
n_surf=n_level, levelfile='./levels')
model.add_ham('matrix', op_type='dip', op_form='top', sparsity_model='indexed',
A=A, B=B, C=C, jmax=jmax, mu_a=mu_a, mu_b=mu_b, mu_c=mu_c,
pulse_pol=pulse_pol, pulse=my_pulse,
n_surf=n_level, levelfile='./levels')
Next, we define the observables that we want to calculate during the propagation. In our case, we are interested in the population of all states in our selected subspace. The following cell specifies that the populations shall be written as individual columns to the file pop.dat.
More details about defining observables are given in the how-to about expectation values and in the observable config file section.
[20]:
for i in range(n_level):
model.add_observable('pop', outfile='pop.dat', time_unit='iu', exp_unit='iu',
exp_surf=i+1, col_label=f'{i+1}')
As a final step, we need to set the parameters for the propagation …
[21]:
model.set_propagation(T=t_stop, nt=nt, t0=t_start, prop_method='cheby', time_unit='iu',
initial_state=psi)
… and save the config file together with the file that specifies the rotational subspace to the runfolder.
[22]:
np.savetxt(rf_prop/'levels', levels, fmt='%5i',
header=f"{'j':>3} {'tau':>5}")
model.write_to_runfolder(rf_prop, config='config')
If we take a look into the runfolder, we see that write_to_runfolder has also saved the initial state and the pulse amplitude:
[23]:
! ls rf_prop
config levels psi.dat pulse1.dat
Hence, we have everything we need to start the propagation.
Propgation and results¶
For the propagation, we use the qdyn_prop_traj utility that comes as part of QDYN.
[24]:
! qdyn_prop_traj rf_prop
***** START OF PROGRAM qdyn_prop_traj ******
Mon Jul 19 11:20:22 +0200 2021
*** Read config file rf_prop/config ***
*** Done reading config file ***
*** Initializing system ***
*** Initializing grid ***
No explicit grid found in config file! Initializing grid as one dimensional 1x1 fake grid
Initializing grid as 1D cartesian grid without mapping
Number of grid points was set to nr = 1 for dimension dim = 1! Initializing dimension 1 as a 1x1 fake-grid...
*** Initializing pulses ***
1 pulses in config file
Initializing pulse 1
*** Initializing dynamical generator ***
Calculating rigid rotor dipole interaction for z-polarized field
Initialising wigner symbols on 16 threads:
*** Done with initialization ***
Completed traj 1/1 in proc 1/1
Writing avg in proc 1/1
Mon Jul 19 11:20:27 +0200 2021
***** END OF PROGRAM qdyn_prop_traj ******
After the propagation, the population of the individual levels are stored in pop.dat inside the runfolder. We can load and visualize the results together with the pulse:
[25]:
pop_data = np.loadtxt(rf_prop/'pop.dat', unpack=True, encoding='latin1',
converters={i: qdyn.io.fix_fortran_exponent for i in range(1, n_level+1)})
tgrid = pop_data[0]
pop = pop_data[1:]
fig, ax = plt.subplots()
ax.plot(pgrid, np.abs(my_amplitude/(my_amplitude.max())), color='lightgray')
for i in range(len(pop)):
ax.plot(tgrid, pop[i], label=f'{i+1}')
ax.autoscale(enable=True, axis='x', tight=True)
ax.set_xlabel('time (iu)')
ax.set_ylabel('population')
plt.show()
Ensemble propagation¶
Up to now, we have propagated a single pure state. However, in many situation we have to simulate an incoherent/probabilistic mixture of rotational states. For example, consider an ensemble where initially each state of the lower level is populated with the same probability, or a thermic ensemble where each level is populated according to the Boltzmann statistic.
We have basically two options: 1. Propagate the density matrix of the complete ensemble. At the moment, this would require some Fortran programming and the propagation may take longer if you have a larger number of states. 2. Simulate each rotational level individually and average the results with the corresponding initial probability distribution. This approach has the advantage that individual propagations can run in parallel, but requires more data management.
In the following, we use option 2.
We define two different initial probability distributions:
All (degenerate) states of level \(4_{0,4}=(4,1)\) are equally populated and the states of level \(5_{1,5}=(5,2)\) are empty, i.e. \(p_n=1/9\) for \(n=1,...,9\) and \(p_n=0\) for \(n>9\).
[26]:
p_equal = np.zeros(n_level)
p_equal[0:9] = 1/9
All levels are occupied according to the Boltzmann distribution at temperature T. The normalized population of each level is then \begin{equation} p_n = \exp \left [ - \frac{E_{rot}}{k_B T} \right ] / Q \end{equation} with \(Q = \sum_{n=1}^{n_{level}} p_n\). Here, \(k_B\) is the Boltzmann constant, \(T\) is the temperature in Kelvin, and \(E_{rot}(n) = \mbox{Erot_sub(n)} B h\) is the rotational energy of state \(n\) in Joule.
[27]:
T = 0.5 # in K
Erot_sub = get_subspace_energies(Erot, levels)
p_boltzmann = np.exp(-0.031497*Erot_sub/T)
p_boltzmann = p_boltzmann / p_boltzmann.sum()
You can select the different contributions by uncommenting the corresponding line in the next cell:
[28]:
p = p_equal
#p = p_boltzmann
Next, we set up the runfolders for the propagation of each rotational level. We can re-use the LevelModel from the last section and only need to update the initial state of the propagation.
[29]:
def prepare_ensemble_runfolders():
runfolders = []
for i in range(n_level):
psi = np.zeros(n_level, dtype=np.complex128)
psi[i] = 1
runfolder = Path(f'r{i:03}')
model.set_propagation(T=t_stop, nt=nt, t0=t_start,
prop_method='cheby', time_unit='iu',
initial_state=psi)
model.write_to_runfolder(runfolder, config='config')
np.savetxt(runfolder/'levels', levels, fmt='%5i',
header=f"{'j':>3} {'tau':>5}")
runfolders.append(runfolder)
return runfolders
[30]:
runfolders = prepare_ensemble_runfolders()
Now, we use the multiprocessing package to propagate all runfolders in parallel. With nprocs, you can control the number of processors that shall be used.
[31]:
def prop(runfolder):
"""Propagate a single runfolder."""
return run(['qdyn_prop_traj', runfolder],
capture_output=True, text=True)
def prop_ensemble(runfolders, nprocs):
"""Propagate a list of runfolders in parallel.
`nprocs` specifies the number of processors that shall be used.
"""
with Pool(nprocs) as pool:
procs = pool.map(prop, runfolders)
for proc in procs:
if proc.returncode != 0: print(proc.stdout)
[32]:
prop_ensemble(runfolders, nprocs=4)
After the propagation has finished (which may take a few minutes), we load all populations from the given runfolders and perform the ensemble average. The plot below again shows the population of each level.
[33]:
ensemble_pop = np.zeros((len(runfolders), n_level, nt))
for i, rf in enumerate(runfolders):
ensemble_pop[i] = np.loadtxt(rf/'pop.dat', usecols=range(1,n_level+1), unpack=True, encoding='latin1',
converters={i: qdyn.io.fix_fortran_exponent for i in range(1, n_level+1)})
# Average populations with initial distribution
average_pop = np.sum(p[:,None,None]*ensemble_pop, axis=0)
# Plot averaged populations
fig, ax = plt.subplots()
for i in range(len(average_pop)):
ax.plot(tgrid, average_pop[i], label=f'{i+1}')
ax.autoscale(enable=True, axis='x', tight=True)
ax.set_xlabel('Time (iu)')
ax.set_ylabel('Population')
plt.show()
