6. Rotational Spectra of Hydrogen
The rotational spectrum of hydrogen, a linear molecule, is computed by the phase-space sampling (see also section III A of the paper).
Redundant space |
Potential |
Integrator |
Physical space |
|---|---|---|---|
Cartesian (6 DOFs) with 4 constraints (bond and 3 translations) |
|
|
2D sphere (polar \(\theta\) and azimuthal angle \(\phi\)) |
Rigid rotor: The bond length and the center of mass are fixed by constraints. Since the potential energy is constant for the rigid rotation,
ZeroPotentialis used (see the Constrained Dynamics tutorial), and no Hessian evaluation is needed.Initial states: Coherent states in \((\theta, \phi)\) centered at the reference structure along the x-axis (\(\theta = \pi / 2\), \(\phi = 0\)), with the rotational kinetic energies of 500, 2000 and 5000 \(\textrm{cm}^{-1}\) and the widths of 0.2 rad.
Sampling: The initial positions are fixed at the center of the initial states, and the initial momenta perpendicular to the bond are sampled uniformly in the 2D momentum space (the rotational kinetic energy and the direction of the rotation plane are sampled uniformly).
Coordinate singularities: The trajectories pass through the singular points (\(\theta = 0, \pi\)) depending on the direction of the rotation. The steps where the symplectic condition is not satisfied are excluded.
Reference: The peaks are compared with the analytical rotational levels \(E_J = B J(J+1)\).
# import necessary libraries
import sys
import time
from pathlib import Path
sys.path.insert(0, str(Path.cwd().parent))
import numpy as np
from pyscivr.units import FemtoSec, CmInv, HBar
from pyscivr.potentials import ZeroPotential
from pyscivr.coordinates import InternalCoordinates
from pyscivr.integrators import ShakeRattle
from pyscivr.variations import make_vari_propagator, propagate_variation, make_jacobi_func, to_internal
from pyscivr.basis import get_overlap
from pyscivr.scivr import get_lagrangian, get_action, get_prefactor, get_ta_spectrum_fft
System and coordinates
The system and the coordinates are the same as in the Constrained Dynamics tutorial:
COORD_CONST defines the constraints (bond length and center of mass), and COORD_ALL defines the full set of coordinates, where the polar and azimuthal angles are active.
The Jacobian functions between the redundant and physical spaces are made once and reused for all trajectories.
Sampling of initial conditions
initial_momenta gives the momenta of rigid rotation with the rotational kinetic energy K_rot
in the plane tilted by tilt from the xy-plane, as in the Constrained Dynamics tutorial.
sample_momenta samples K_rot and tilt uniformly, which corresponds to the uniform sampling
of the momenta perpendicular to the bond in the 2D momentum space.
The vibrational kinetic energy K_vib is set to zero for the pure rotation.
Calculation of spectra
run computes the spectra for all the initial states by the following steps for each trajectory:
Propagate the trajectory by
ShakeRattle.Propagate the variation and compute the monodromy matrix in the physical space. The Hessian is zero for
ZeroPotential.Exclude the steps where the symplectic condition is not satisfied.
Compute the action, the prefactor and the overlap with each initial state, and the spectrum by the FFT.
The spectra are averaged over the trajectories and saved for each initial state. Note that each trajectory is used for all the initial states, so that the cost of the dynamics and the monodromy matrix is shared.
# Constants
MASS = np.array([1.008, 1.008]) # mass of hydrogen atoms [amu (atomic mass unit)]
R_E = 0.74160392 # equilibrium bond length [Angs]
POSI_REF = np.array([[-R_E / 2, 0.0, 0.0], [R_E / 2, 0.0, 0.0]])
CALC = ZeroPotential()
ENER_REF = CALC.energy(POSI_REF)
# Coordinate for constraints
COORD_CONST = InternalCoordinates(POSI_REF, MASS)
COORD_CONST.add_bond(0, 1)
COORD_CONST.add_translation()
COORD_CONST.build()
# Coordinate for all DOFs
COORD_ALL = InternalCoordinates(POSI_REF, MASS)
COORD_ALL.add_translation()
COORD_ALL.add_polar(0, 1)
COORD_ALL.add_azimuthal(0, 1)
COORD_ALL.add_bond(0, 1)
COORD_ALL.build()
ACTIVE_IDX = np.array([COORD_ALL.find_coordinate(kind, (0, 1))[0] for kind in ("polar", "azimuthal")])
TOPOLOGIES = [COORD_ALL.topology[i] for i in ACTIVE_IDX]
JACOBI_FUNC, JACOBI_REV_FUNC, INTERNAL_FUNC, MONODROMY_FUNC = make_jacobi_func(COORD_ALL, ACTIVE_IDX)
def initial_momenta(K_rot, K_vib, tilt, sign=1):
p_rot = np.sqrt(K_rot * CmInv * MASS[0])
p_vib = np.sqrt(K_vib * CmInv * MASS[0]) * sign
p_rot_y = p_rot * np.cos(tilt)
p_rot_z = p_rot * np.sin(tilt)
return np.array([[-p_vib, -p_rot_y, -p_rot_z], [p_vib, p_rot_y, p_rot_z]])
def sample_momenta(n_sample, K_rot_range, K_vib_range, seed=10202):
rng = np.random.default_rng(seed)
K_rot = rng.uniform(*K_rot_range, size=n_sample)
K_vib = rng.uniform(*K_vib_range, size=n_sample)
sign = rng.choice([-1, 1], size=n_sample)
tilt = rng.uniform(0.0, 2.0 * np.pi, size=n_sample)
return [initial_momenta(K_rot[i], K_vib[i], tilt[i], sign[i]) for i in range(n_sample)], K_rot, K_vib
def run(posi_sample_list, mome_sample_list, K_rot_list, K_vib_list,
posi_ref_list, mome_ref_list, gamma_ref_list, fpath_list,
dt, n_step, n_fft, symp_err_tol=1e-4):
# Judge function to check if the point is singular or not
judge_func= lambda symp_err : np.isfinite(symp_err) & (symp_err < symp_err_tol)
# Reference positions and momenta in the internal coordinates
posi_ref_int, mome_ref_int = zip(*[INTERNAL_FUNC(p, q) for p, q in zip(posi_ref_list, mome_ref_list)])
n_sample = len(posi_sample_list)
n_ref = len(posi_ref_list)
spectra = [0.0] * n_ref
vari_propagator = make_vari_propagator(COORD_CONST, dt)
t_start = time.perf_counter()
for i in range(n_sample):
# Get trajectory
dyn = ShakeRattle(CALC, posi_sample_list[i], mome_sample_list[i], MASS, dt, COORD_CONST)
traj = dyn.run(n_step)
traj_posi, traj_mome = traj["posi"], traj["mome"]
# get variation
hess = CALC.hessian(traj_posi)
J_rev_init = JACOBI_REV_FUNC(traj_posi[0], traj_mome[0])
vari = propagate_variation(vari_propagator, traj_posi, traj_mome, traj["lam"], traj["mu"], hess, vari_init=J_rev_init)
posi_int, mome_int, mono_int, symplectic_error = to_internal(MONODROMY_FUNC, traj_posi, traj_mome, vari)
# points with large symplectic error are treated as singularities
valid = judge_func(symplectic_error)
ener_grid = 2.0 * np.pi * HBar * np.arange(n_fft) / (dt * n_fft)
lag = get_lagrangian(traj["total_energy"], traj["potential_energy"], ener_ref=ENER_REF)
action = get_action(lag, dt)
for j in range(n_ref):
posi_center, mome_center = posi_ref_int[j], mome_ref_int[j]
gamma = gamma_ref_list[j]
prefactor = get_prefactor(mono_int, gamma, valid=valid)
overlap = get_overlap(TOPOLOGIES, posi_center, mome_center, posi_int, mome_int, gamma)
_, spec = get_ta_spectrum_fft(prefactor, overlap, action, dt, n_fft, valid=valid)
spectra[j] += spec
if i % 10 == 0:
print(f"{i+1:6d} / {n_sample}, K_rot = {K_rot_list[i]:.1f}, K_vib = {K_vib_list[i]:.1f}, "
f"max symp error = {np.max(np.abs(symplectic_error)):.3e}, ({time.perf_counter() - t_start:.1f} sec)", flush=True)
for j in range(n_ref):
spectra[j] /= n_sample
for j in range(n_ref):
np.savez(fpath_list[j], ener_grid=ener_grid, spectrum=spectra[j])
Build complete: 4 coordinates
Build complete: 6 coordinates
Execute calculation
The total propagation time is 500 fs.
if __name__ == "__main__":
# set parameters
n_sample = 2000
dt = 0.1 * FemtoSec
n_step = 5000
n_fft = 2**19
# set reference state
K_rot_ref_list = [500.0, 2000.0, 5000.0]
mome_ref_list = [initial_momenta(k_rot, 0.0, np.pi / 4) for k_rot in K_rot_ref_list]
posi_ref_list = [POSI_REF] * len(K_rot_ref_list)
gamma_ref = [1.0 / (2 * 0.2 ** 2), # sqrt(1 / 2 gamma) = 0.2 rad
1.0 / (2 * 0.2 ** 2)]
gamma_ref_list = [gamma_ref] * len(K_rot_ref_list)
fpath_list = [f"hydrogen_pure_rot_{k_rot:.0f}_n_{n_sample}_spectra.npz" for k_rot in K_rot_ref_list]
posi_ref_list = np.asarray(posi_ref_list)
mome_ref_list = np.asarray(mome_ref_list)
gamma_ref_list = np.asarray(gamma_ref_list)
K_rot_range = (10.0, 10000.0)
K_vib_range = (0.0, 0.0) # No vibrational kinetic energy for pure rotational motion
mome_sample_list, K_rot_list, K_vib_list = sample_momenta(n_sample, K_rot_range, K_vib_range)
posi_sample_list = np.asarray([POSI_REF] * n_sample)
run(posi_sample_list, mome_sample_list, K_rot_list, K_vib_list,
posi_ref_list, mome_ref_list, gamma_ref_list, fpath_list,
dt, n_step, n_fft)
1 / 2000, K_rot = 2788.2, K_vib = 0.0, max symp error = 2.680e-09, (1.2 sec)
11 / 2000, K_rot = 4946.9, K_vib = 0.0, max symp error = 3.705e-10, (5.6 sec)
21 / 2000, K_rot = 9808.5, K_vib = 0.0, max symp error = 2.921e-10, (10.0 sec)
31 / 2000, K_rot = 8930.5, K_vib = 0.0, max symp error = 1.116e-09, (14.7 sec)
41 / 2000, K_rot = 1592.5, K_vib = 0.0, max symp error = 7.402e-10, (19.2 sec)
51 / 2000, K_rot = 7396.6, K_vib = 0.0, max symp error = 1.458e-10, (23.8 sec)
61 / 2000, K_rot = 2192.0, K_vib = 0.0, max symp error = 1.056e-09, (28.5 sec)
71 / 2000, K_rot = 6597.3, K_vib = 0.0, max symp error = 5.435e-10, (33.0 sec)
81 / 2000, K_rot = 2769.8, K_vib = 0.0, max symp error = 1.855e-09, (37.6 sec)
91 / 2000, K_rot = 2161.3, K_vib = 0.0, max symp error = 1.721e-09, (42.1 sec)
101 / 2000, K_rot = 3024.0, K_vib = 0.0, max symp error = 1.338e-09, (46.7 sec)
111 / 2000, K_rot = 4798.8, K_vib = 0.0, max symp error = 8.624e-10, (51.1 sec)
121 / 2000, K_rot = 6286.0, K_vib = 0.0, max symp error = 1.251e-09, (55.6 sec)
131 / 2000, K_rot = 5551.3, K_vib = 0.0, max symp error = 1.501e-09, (60.1 sec)
141 / 2000, K_rot = 4030.5, K_vib = 0.0, max symp error = 1.468e-10, (64.8 sec)
151 / 2000, K_rot = 1705.5, K_vib = 0.0, max symp error = 1.197e-09, (69.5 sec)
161 / 2000, K_rot = 2817.6, K_vib = 0.0, max symp error = 1.689e-09, (74.2 sec)
171 / 2000, K_rot = 3701.8, K_vib = 0.0, max symp error = 2.316e-09, (78.8 sec)
181 / 2000, K_rot = 7288.7, K_vib = 0.0, max symp error = 6.891e-10, (83.5 sec)
191 / 2000, K_rot = 6585.2, K_vib = 0.0, max symp error = 1.260e-09, (88.1 sec)
201 / 2000, K_rot = 7396.7, K_vib = 0.0, max symp error = 1.210e-09, (92.6 sec)
211 / 2000, K_rot = 9867.3, K_vib = 0.0, max symp error = 8.693e-10, (97.1 sec)
221 / 2000, K_rot = 7455.7, K_vib = 0.0, max symp error = 1.696e-09, (101.6 sec)
231 / 2000, K_rot = 3787.1, K_vib = 0.0, max symp error = 9.178e-10, (106.2 sec)
241 / 2000, K_rot = 3470.1, K_vib = 0.0, max symp error = 5.026e-11, (110.9 sec)
251 / 2000, K_rot = 1747.2, K_vib = 0.0, max symp error = 9.645e-11, (115.7 sec)
261 / 2000, K_rot = 9973.6, K_vib = 0.0, max symp error = 5.150e-10, (120.4 sec)
271 / 2000, K_rot = 4377.4, K_vib = 0.0, max symp error = 4.803e-10, (125.2 sec)
281 / 2000, K_rot = 9515.0, K_vib = 0.0, max symp error = 2.860e-10, (130.0 sec)
291 / 2000, K_rot = 6893.6, K_vib = 0.0, max symp error = 9.438e-10, (134.7 sec)
301 / 2000, K_rot = 5376.0, K_vib = 0.0, max symp error = 1.139e-09, (139.3 sec)
311 / 2000, K_rot = 5321.8, K_vib = 0.0, max symp error = 1.367e-09, (144.0 sec)
321 / 2000, K_rot = 9159.7, K_vib = 0.0, max symp error = 9.742e-10, (148.6 sec)
331 / 2000, K_rot = 5392.6, K_vib = 0.0, max symp error = 1.207e-09, (153.4 sec)
341 / 2000, K_rot = 487.6, K_vib = 0.0, max symp error = 3.172e-10, (158.2 sec)
351 / 2000, K_rot = 6850.8, K_vib = 0.0, max symp error = 1.104e-09, (163.0 sec)
361 / 2000, K_rot = 219.6, K_vib = 0.0, max symp error = 1.266e-09, (167.7 sec)
371 / 2000, K_rot = 7582.7, K_vib = 0.0, max symp error = 4.256e-10, (172.5 sec)
381 / 2000, K_rot = 7779.4, K_vib = 0.0, max symp error = 8.252e-10, (177.2 sec)
391 / 2000, K_rot = 4940.4, K_vib = 0.0, max symp error = 6.222e-10, (181.9 sec)
401 / 2000, K_rot = 8660.0, K_vib = 0.0, max symp error = 4.726e-10, (186.6 sec)
411 / 2000, K_rot = 1841.7, K_vib = 0.0, max symp error = 1.018e-09, (191.4 sec)
421 / 2000, K_rot = 2463.3, K_vib = 0.0, max symp error = 5.678e-10, (196.0 sec)
431 / 2000, K_rot = 7594.3, K_vib = 0.0, max symp error = 1.289e-09, (200.7 sec)
441 / 2000, K_rot = 7529.8, K_vib = 0.0, max symp error = 4.533e-10, (205.3 sec)
451 / 2000, K_rot = 3684.4, K_vib = 0.0, max symp error = 8.943e-10, (209.9 sec)
461 / 2000, K_rot = 390.8, K_vib = 0.0, max symp error = 5.925e-10, (214.5 sec)
471 / 2000, K_rot = 1503.1, K_vib = 0.0, max symp error = 1.107e-09, (219.1 sec)
481 / 2000, K_rot = 2764.9, K_vib = 0.0, max symp error = 1.081e-09, (223.7 sec)
491 / 2000, K_rot = 5392.0, K_vib = 0.0, max symp error = 8.092e-11, (228.4 sec)
501 / 2000, K_rot = 2453.7, K_vib = 0.0, max symp error = 2.457e-09, (233.1 sec)
511 / 2000, K_rot = 9722.6, K_vib = 0.0, max symp error = 9.067e-10, (237.8 sec)
521 / 2000, K_rot = 2164.7, K_vib = 0.0, max symp error = 5.723e-10, (242.5 sec)
531 / 2000, K_rot = 1284.9, K_vib = 0.0, max symp error = 7.824e-10, (247.1 sec)
541 / 2000, K_rot = 3925.6, K_vib = 0.0, max symp error = 1.650e-09, (251.9 sec)
551 / 2000, K_rot = 9981.7, K_vib = 0.0, max symp error = 2.264e-10, (256.6 sec)
561 / 2000, K_rot = 5556.3, K_vib = 0.0, max symp error = 1.349e-09, (261.3 sec)
571 / 2000, K_rot = 7794.8, K_vib = 0.0, max symp error = 9.959e-10, (265.9 sec)
581 / 2000, K_rot = 3518.4, K_vib = 0.0, max symp error = 6.393e-10, (270.5 sec)
591 / 2000, K_rot = 3835.0, K_vib = 0.0, max symp error = 5.076e-10, (275.1 sec)
601 / 2000, K_rot = 817.9, K_vib = 0.0, max symp error = 1.835e-09, (279.7 sec)
611 / 2000, K_rot = 2152.6, K_vib = 0.0, max symp error = 6.105e-10, (284.4 sec)
621 / 2000, K_rot = 2828.7, K_vib = 0.0, max symp error = 6.351e-10, (289.2 sec)
631 / 2000, K_rot = 3203.5, K_vib = 0.0, max symp error = 2.871e-10, (293.9 sec)
641 / 2000, K_rot = 4838.7, K_vib = 0.0, max symp error = 7.875e-11, (298.6 sec)
651 / 2000, K_rot = 4176.6, K_vib = 0.0, max symp error = 7.929e-10, (303.2 sec)
661 / 2000, K_rot = 2076.9, K_vib = 0.0, max symp error = 7.046e-10, (307.7 sec)
671 / 2000, K_rot = 7519.9, K_vib = 0.0, max symp error = 1.352e-09, (312.2 sec)
681 / 2000, K_rot = 9171.1, K_vib = 0.0, max symp error = 7.952e-10, (316.9 sec)
691 / 2000, K_rot = 5976.7, K_vib = 0.0, max symp error = 4.052e-10, (321.6 sec)
701 / 2000, K_rot = 5310.1, K_vib = 0.0, max symp error = 1.088e-09, (326.3 sec)
711 / 2000, K_rot = 7895.5, K_vib = 0.0, max symp error = 2.454e-10, (330.9 sec)
721 / 2000, K_rot = 9576.3, K_vib = 0.0, max symp error = 1.343e-09, (335.7 sec)
731 / 2000, K_rot = 9017.8, K_vib = 0.0, max symp error = 1.400e-09, (340.3 sec)
741 / 2000, K_rot = 6847.5, K_vib = 0.0, max symp error = 2.638e-10, (345.0 sec)
751 / 2000, K_rot = 6915.3, K_vib = 0.0, max symp error = 1.021e-09, (349.8 sec)
761 / 2000, K_rot = 2417.0, K_vib = 0.0, max symp error = 6.727e-10, (354.7 sec)
771 / 2000, K_rot = 4178.8, K_vib = 0.0, max symp error = 2.217e-09, (359.4 sec)
781 / 2000, K_rot = 9218.0, K_vib = 0.0, max symp error = 4.434e-11, (364.1 sec)
791 / 2000, K_rot = 2660.2, K_vib = 0.0, max symp error = 9.337e-10, (368.7 sec)
801 / 2000, K_rot = 2507.7, K_vib = 0.0, max symp error = 1.147e-09, (373.6 sec)
811 / 2000, K_rot = 5999.4, K_vib = 0.0, max symp error = 5.702e-10, (378.4 sec)
821 / 2000, K_rot = 282.7, K_vib = 0.0, max symp error = 1.887e-09, (383.3 sec)
831 / 2000, K_rot = 6042.2, K_vib = 0.0, max symp error = 1.155e-09, (388.1 sec)
841 / 2000, K_rot = 3275.4, K_vib = 0.0, max symp error = 5.697e-10, (392.9 sec)
851 / 2000, K_rot = 2054.4, K_vib = 0.0, max symp error = 5.034e-10, (397.6 sec)
861 / 2000, K_rot = 4978.0, K_vib = 0.0, max symp error = 1.589e-09, (402.3 sec)
871 / 2000, K_rot = 4786.3, K_vib = 0.0, max symp error = 1.403e-09, (407.0 sec)
881 / 2000, K_rot = 3213.5, K_vib = 0.0, max symp error = 1.004e-09, (411.8 sec)
891 / 2000, K_rot = 1980.6, K_vib = 0.0, max symp error = 1.929e-10, (416.7 sec)
901 / 2000, K_rot = 9113.5, K_vib = 0.0, max symp error = 1.709e-09, (421.5 sec)
911 / 2000, K_rot = 7299.0, K_vib = 0.0, max symp error = 1.198e-09, (426.2 sec)
921 / 2000, K_rot = 4054.4, K_vib = 0.0, max symp error = 8.651e-10, (430.9 sec)
931 / 2000, K_rot = 9594.0, K_vib = 0.0, max symp error = 5.862e-10, (435.5 sec)
941 / 2000, K_rot = 7537.2, K_vib = 0.0, max symp error = 7.554e-10, (440.2 sec)
951 / 2000, K_rot = 9288.8, K_vib = 0.0, max symp error = 5.175e-10, (444.9 sec)
961 / 2000, K_rot = 2782.4, K_vib = 0.0, max symp error = 1.324e-09, (449.5 sec)
971 / 2000, K_rot = 7003.2, K_vib = 0.0, max symp error = 2.273e-10, (454.2 sec)
981 / 2000, K_rot = 2347.4, K_vib = 0.0, max symp error = 9.900e-10, (458.9 sec)
991 / 2000, K_rot = 3765.6, K_vib = 0.0, max symp error = 1.279e-09, (463.6 sec)
1001 / 2000, K_rot = 9447.7, K_vib = 0.0, max symp error = 1.079e-09, (468.3 sec)
1011 / 2000, K_rot = 4098.7, K_vib = 0.0, max symp error = 3.597e-10, (473.0 sec)
1021 / 2000, K_rot = 2402.9, K_vib = 0.0, max symp error = 7.812e-10, (477.7 sec)
1031 / 2000, K_rot = 7194.7, K_vib = 0.0, max symp error = 1.044e-09, (482.4 sec)
1041 / 2000, K_rot = 2557.9, K_vib = 0.0, max symp error = 2.052e-09, (487.3 sec)
1051 / 2000, K_rot = 7609.2, K_vib = 0.0, max symp error = 3.957e-10, (492.0 sec)
1061 / 2000, K_rot = 4817.1, K_vib = 0.0, max symp error = 1.236e-09, (496.8 sec)
1071 / 2000, K_rot = 5458.6, K_vib = 0.0, max symp error = 3.431e-10, (501.5 sec)
1081 / 2000, K_rot = 4589.1, K_vib = 0.0, max symp error = 1.135e-09, (506.2 sec)
1091 / 2000, K_rot = 5344.1, K_vib = 0.0, max symp error = 1.134e-09, (510.9 sec)
1101 / 2000, K_rot = 9410.6, K_vib = 0.0, max symp error = 5.722e-10, (515.6 sec)
1111 / 2000, K_rot = 3666.5, K_vib = 0.0, max symp error = 8.712e-10, (520.4 sec)
1121 / 2000, K_rot = 240.3, K_vib = 0.0, max symp error = 4.725e-10, (525.2 sec)
1131 / 2000, K_rot = 6294.3, K_vib = 0.0, max symp error = 1.975e-09, (529.9 sec)
1141 / 2000, K_rot = 4357.6, K_vib = 0.0, max symp error = 9.068e-10, (534.8 sec)
1151 / 2000, K_rot = 9056.6, K_vib = 0.0, max symp error = 1.128e-09, (539.6 sec)
1161 / 2000, K_rot = 7166.2, K_vib = 0.0, max symp error = 2.227e-11, (544.4 sec)
1171 / 2000, K_rot = 8916.5, K_vib = 0.0, max symp error = 1.021e-09, (549.1 sec)
1181 / 2000, K_rot = 451.6, K_vib = 0.0, max symp error = 4.935e-10, (553.9 sec)
1191 / 2000, K_rot = 2126.8, K_vib = 0.0, max symp error = 1.139e-09, (558.5 sec)
1201 / 2000, K_rot = 9789.0, K_vib = 0.0, max symp error = 8.513e-10, (563.1 sec)
1211 / 2000, K_rot = 5422.0, K_vib = 0.0, max symp error = 1.167e-10, (567.7 sec)
1221 / 2000, K_rot = 3624.6, K_vib = 0.0, max symp error = 3.619e-10, (572.3 sec)
1231 / 2000, K_rot = 2296.3, K_vib = 0.0, max symp error = 7.671e-10, (576.9 sec)
1241 / 2000, K_rot = 2172.9, K_vib = 0.0, max symp error = 1.472e-09, (581.5 sec)
1251 / 2000, K_rot = 5639.4, K_vib = 0.0, max symp error = 5.967e-10, (586.3 sec)
1261 / 2000, K_rot = 4547.6, K_vib = 0.0, max symp error = 9.726e-10, (591.1 sec)
1271 / 2000, K_rot = 350.1, K_vib = 0.0, max symp error = 2.029e-09, (595.8 sec)
1281 / 2000, K_rot = 3161.5, K_vib = 0.0, max symp error = 8.883e-10, (600.6 sec)
1291 / 2000, K_rot = 2615.2, K_vib = 0.0, max symp error = 9.626e-10, (605.4 sec)
1301 / 2000, K_rot = 3493.6, K_vib = 0.0, max symp error = 1.142e-09, (610.1 sec)
1311 / 2000, K_rot = 2293.5, K_vib = 0.0, max symp error = 6.199e-10, (614.9 sec)
1321 / 2000, K_rot = 1026.4, K_vib = 0.0, max symp error = 2.312e-09, (619.7 sec)
1331 / 2000, K_rot = 8248.1, K_vib = 0.0, max symp error = 3.949e-10, (624.5 sec)
1341 / 2000, K_rot = 9451.1, K_vib = 0.0, max symp error = 1.147e-09, (629.4 sec)
1351 / 2000, K_rot = 6264.2, K_vib = 0.0, max symp error = 6.107e-10, (634.3 sec)
1361 / 2000, K_rot = 4614.2, K_vib = 0.0, max symp error = 3.291e-10, (639.1 sec)
1371 / 2000, K_rot = 8100.1, K_vib = 0.0, max symp error = 8.619e-10, (643.9 sec)
1381 / 2000, K_rot = 6091.0, K_vib = 0.0, max symp error = 3.813e-10, (648.7 sec)
1391 / 2000, K_rot = 3336.7, K_vib = 0.0, max symp error = 2.489e-10, (653.5 sec)
1401 / 2000, K_rot = 9340.1, K_vib = 0.0, max symp error = 2.215e-09, (658.3 sec)
1411 / 2000, K_rot = 3087.2, K_vib = 0.0, max symp error = 1.210e-10, (663.0 sec)
1421 / 2000, K_rot = 6007.2, K_vib = 0.0, max symp error = 2.716e-10, (667.7 sec)
1431 / 2000, K_rot = 8549.9, K_vib = 0.0, max symp error = 1.224e-09, (672.5 sec)
1441 / 2000, K_rot = 3685.4, K_vib = 0.0, max symp error = 1.428e-09, (677.2 sec)
1451 / 2000, K_rot = 1389.7, K_vib = 0.0, max symp error = 1.039e-09, (681.9 sec)
1461 / 2000, K_rot = 5130.7, K_vib = 0.0, max symp error = 7.576e-10, (686.8 sec)
1471 / 2000, K_rot = 866.6, K_vib = 0.0, max symp error = 1.269e-09, (691.4 sec)
1481 / 2000, K_rot = 7205.9, K_vib = 0.0, max symp error = 1.365e-09, (696.1 sec)
1491 / 2000, K_rot = 5494.3, K_vib = 0.0, max symp error = 7.108e-10, (700.9 sec)
1501 / 2000, K_rot = 6408.4, K_vib = 0.0, max symp error = 1.753e-10, (705.6 sec)
1511 / 2000, K_rot = 2938.0, K_vib = 0.0, max symp error = 5.596e-10, (710.5 sec)
1521 / 2000, K_rot = 9885.7, K_vib = 0.0, max symp error = 1.164e-09, (715.1 sec)
1531 / 2000, K_rot = 5770.4, K_vib = 0.0, max symp error = 3.947e-10, (719.9 sec)
1541 / 2000, K_rot = 9285.6, K_vib = 0.0, max symp error = 2.698e-09, (724.5 sec)
1551 / 2000, K_rot = 1669.5, K_vib = 0.0, max symp error = 4.818e-10, (729.3 sec)
1561 / 2000, K_rot = 8796.3, K_vib = 0.0, max symp error = 6.143e-10, (734.1 sec)
1571 / 2000, K_rot = 2098.4, K_vib = 0.0, max symp error = 3.571e-11, (738.8 sec)
1581 / 2000, K_rot = 6521.5, K_vib = 0.0, max symp error = 1.607e-09, (743.4 sec)
1591 / 2000, K_rot = 8306.3, K_vib = 0.0, max symp error = 1.836e-09, (748.2 sec)
1601 / 2000, K_rot = 5014.2, K_vib = 0.0, max symp error = 1.455e-09, (752.9 sec)
1611 / 2000, K_rot = 2860.3, K_vib = 0.0, max symp error = 7.849e-10, (757.6 sec)
1621 / 2000, K_rot = 1863.2, K_vib = 0.0, max symp error = 1.546e-10, (762.4 sec)
1631 / 2000, K_rot = 6431.5, K_vib = 0.0, max symp error = 1.173e-09, (767.0 sec)
1641 / 2000, K_rot = 4589.4, K_vib = 0.0, max symp error = 1.086e-09, (771.8 sec)
1651 / 2000, K_rot = 3933.0, K_vib = 0.0, max symp error = 2.092e-10, (776.6 sec)
1661 / 2000, K_rot = 9232.2, K_vib = 0.0, max symp error = 3.533e-10, (781.4 sec)
1671 / 2000, K_rot = 7627.6, K_vib = 0.0, max symp error = 5.849e-10, (786.2 sec)
1681 / 2000, K_rot = 6536.5, K_vib = 0.0, max symp error = 5.926e-10, (790.8 sec)
1691 / 2000, K_rot = 897.3, K_vib = 0.0, max symp error = 1.001e-09, (795.5 sec)
1701 / 2000, K_rot = 1867.7, K_vib = 0.0, max symp error = 1.550e-09, (800.2 sec)
1711 / 2000, K_rot = 1623.7, K_vib = 0.0, max symp error = 4.588e-10, (805.1 sec)
1721 / 2000, K_rot = 5538.5, K_vib = 0.0, max symp error = 8.589e-10, (809.8 sec)
1731 / 2000, K_rot = 5618.3, K_vib = 0.0, max symp error = 1.106e-09, (814.7 sec)
1741 / 2000, K_rot = 8806.9, K_vib = 0.0, max symp error = 1.843e-09, (819.4 sec)
1751 / 2000, K_rot = 5849.4, K_vib = 0.0, max symp error = 1.731e-09, (824.1 sec)
1761 / 2000, K_rot = 192.3, K_vib = 0.0, max symp error = 1.355e-09, (828.9 sec)
1771 / 2000, K_rot = 4075.3, K_vib = 0.0, max symp error = 2.365e-09, (833.7 sec)
1781 / 2000, K_rot = 5590.7, K_vib = 0.0, max symp error = 1.645e-09, (838.4 sec)
1791 / 2000, K_rot = 2198.4, K_vib = 0.0, max symp error = 1.520e-09, (843.2 sec)
1801 / 2000, K_rot = 6597.9, K_vib = 0.0, max symp error = 1.272e-09, (848.0 sec)
1811 / 2000, K_rot = 4280.0, K_vib = 0.0, max symp error = 1.775e-10, (852.7 sec)
1821 / 2000, K_rot = 4019.1, K_vib = 0.0, max symp error = 1.308e-09, (857.5 sec)
1831 / 2000, K_rot = 3659.6, K_vib = 0.0, max symp error = 6.542e-10, (862.2 sec)
1841 / 2000, K_rot = 4965.5, K_vib = 0.0, max symp error = 2.012e-09, (866.9 sec)
1851 / 2000, K_rot = 2757.4, K_vib = 0.0, max symp error = 1.790e-10, (871.8 sec)
1861 / 2000, K_rot = 9182.1, K_vib = 0.0, max symp error = 1.243e-09, (876.5 sec)
1871 / 2000, K_rot = 8172.8, K_vib = 0.0, max symp error = 6.560e-10, (881.0 sec)
1881 / 2000, K_rot = 3323.9, K_vib = 0.0, max symp error = 5.745e-10, (885.6 sec)
1891 / 2000, K_rot = 5787.3, K_vib = 0.0, max symp error = 5.782e-10, (890.1 sec)
1901 / 2000, K_rot = 128.8, K_vib = 0.0, max symp error = 8.916e-10, (894.6 sec)
1911 / 2000, K_rot = 4775.1, K_vib = 0.0, max symp error = 3.214e-10, (899.2 sec)
1921 / 2000, K_rot = 7999.3, K_vib = 0.0, max symp error = 7.878e-10, (903.9 sec)
1931 / 2000, K_rot = 650.4, K_vib = 0.0, max symp error = 5.829e-10, (908.7 sec)
1941 / 2000, K_rot = 2412.1, K_vib = 0.0, max symp error = 1.661e-09, (913.4 sec)
1951 / 2000, K_rot = 9588.3, K_vib = 0.0, max symp error = 1.619e-10, (918.2 sec)
1961 / 2000, K_rot = 9360.6, K_vib = 0.0, max symp error = 1.845e-09, (922.8 sec)
1971 / 2000, K_rot = 9367.1, K_vib = 0.0, max symp error = 8.983e-10, (927.5 sec)
1981 / 2000, K_rot = 9965.5, K_vib = 0.0, max symp error = 1.212e-09, (932.2 sec)
1991 / 2000, K_rot = 2330.2, K_vib = 0.0, max symp error = 8.930e-10, (936.9 sec)
Plot results
The reference rotational energies can be calculated analytically using the formula \(E_J = B_e J(J+1)\), where \(B_e = \frac{\hbar^2}{2I}\) is the rotational constant and \(J\) is the rotational quantum number.
The moment of inertia \(I\) for H2 can be calculated as \(I = \frac{m_H}{2}r_e^2\).
if __name__ == "__main__":
import matplotlib.pyplot as plt
from scipy.signal import find_peaks
# moment of inertia of H2
I = MASS[0] * R_E**2 / 2.0
# rotational constant
B = HBar**2 / (2.0 * I) / CmInv
analytical_energies = np.array([B * J * (J + 1) for J in range(20) if B * J * (J + 1) <= 10000])
print(f"Rotational constant B = {B:.3f} cm^-1")
print(f"Analytical rotational energies (cm^-1): {np.round(analytical_energies, 1)}")
fig, ax = plt.subplots(3, 1, figsize=(7, 7), sharex=True)
colors = ["red", "green", "blue"]
for i, fpath in enumerate(fpath_list):
data = np.load(fpath)
ener_grid = data["ener_grid"] / CmInv
spectrum = data["spectrum"]
spectrum /= np.max(spectrum)
peaks, _ = find_peaks(spectrum, height=0.1, prominence=0.1)
peaks = peaks[ener_grid[peaks] <= 10000]
print(f"Peaks (cm^-1) for K = {K_rot_ref_list[i]:.0f}: {np.round(ener_grid[peaks], 1)}")
ax[i].plot(ener_grid[peaks], spectrum[peaks], "x", color="black")
ax[i].plot(ener_grid, spectrum, label=f"K = {K_rot_ref_list[i]:.0f}", color=colors[i])
ax[i].set_xlim(0, 10000)
ax[i].set_ylim(0, 1.2)
ax[i].vlines(analytical_energies, 0.0, 1.2, color="gray", linestyle="--", label="analytical")
ax[i].legend()
ax[-1].set_xlabel("Energy (cm^-1)")
ax[1].set_ylabel("intensity (arb. unit)")
ax[-1].set_xlim(0, 10000)
plt.tight_layout()
plt.show()
Rotational constant B = 60.817 cm^-1
Analytical rotational energies (cm^-1): [ 0. 121.6 364.9 729.8 1216.3 1824.5 2554.3 3405.7 4378.8 5473.5
6689.8 8027.8 9487.4]
Peaks (cm^-1) for K = 500: [ 111.3 370.3 730.4 1214.5 1826.6 2553.8 3403.8 4377.9]
Peaks (cm^-1) for K = 2000: [ 369. 730.4 1213.3 1825.3 2553.2 3403.2 4376.6 5475.3 6686.1]
Peaks (cm^-1) for K = 5000: [1212. 1824.7 2552.5 3403.2 4375.9 5474.7 6685.4 8028.5 9494.4]
The computed peak positions for \(K = 2000\) \(\textrm{cm}^{-1}\) are compared with the analytical eigenenergies \(E_J = B_e J(J+1)\) with \(B_e = 60.817\) \(\textrm{cm}^{-1}\) (all in \(\textrm{cm}^{-1}\)).
\(J\) |
Computed |
\(B_e J(J+1)\) |
Deviation |
|---|---|---|---|
2 |
369.0 |
364.9 |
+4.1 |
3 |
730.4 |
729.8 |
+0.6 |
4 |
1213.3 |
1216.3 |
−3.0 |
5 |
1825.3 |
1824.5 |
+0.8 |
6 |
2553.2 |
2554.3 |
−1.1 |
7 |
3403.2 |
3405.7 |
−2.5 |
8 |
4376.6 |
4378.8 |
−2.2 |
9 |
5475.3 |
5473.5 |
+1.8 |
10 |
6686.1 |
6689.8 |
−3.7 |
Quantitative agreement: The computed peaks agree with the analytical eigenenergies with an RMSE of about 2.5 \(\textrm{cm}^{-1}\), without systematic dependence on \(J\). The same peak positions are obtained for all the initial states, up to \(J = 12\) (about 9500 \(\textrm{cm}^{-1}\)) for \(K\) = 5000 \(\textrm{cm}^{-1}\).
Line shape: Each peak has an asymmetric line shape with a sharp edge at the eigenenergy and a tail on the low-energy side, which originates from the time-averaging expression (see Appendix K of the paper).
Dependence on the initial state: The intensity distribution follows the rotational kinetic energy of the initial state: the strongest peaks appear around 730, 2550 and 5470 \(\textrm{cm}^{-1}\) for the initial states with \(K\) = 500, 2000 and 5000 \(\textrm{cm}^{-1}\), respectively.