7. Ro-vibrational Spectra of Hydrogen

The ro-vibrational spectrum of hydrogen, a linear molecule, is computed by the phase-space sampling (see also section III D of the paper).

Redundant space	Potential	Integrator	Physical space
Cartesian (6 DOFs) with 3 constraints (3 translations)	H2SplinePotential	ShakeRattle	bond length \(r\), polar \(\theta\) and azimuthal angle \(\phi\)

• Ro-vibration: Only the center of mass is fixed by constraints, so that the molecule vibrates while rotating. The analytical Hessian of H2SplinePotential is used at every step.

• Initial states: Coherent states centered at the equilibrium structure, in the vibrational ground state (\(v = 0\)) with the rotational kinetic energies of 500, 2000 and 5000 \(\textrm{cm}^{-1}\).

• Sampling: As in the pure rotation, the initial positions are fixed at the equilibrium structure, and only the momenta are sampled, including the vibrational momenta (see below).

• Reference: The peaks are compared with the Dunham expansion for \(v = 0\) using the experimental parameters.

# 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 H2SplinePotential
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 Hessian and Variation tutorial: COORD_CONST defines the constraints (center of mass), and COORD_ALL defines the full set of coordinates, where the bond length and the polar and azimuthal angles are active.

Sampling of initial conditions

sample_momenta samples K_rot, tilt, the vibrational kinetic energy K_vib and the sign of the vibrational momentum uniformly.

The vibrational momenta are sampled even for the vibrational ground state. In the ground state, the molecule still vibrates with the zero-point motion, and the rotational levels are determined by the rotational constant averaged over this motion, \(B_0\), which is smaller than that at the equilibrium structure, \(B_e\). Sampling the vibrational momenta takes this ro-vibrational coupling into account.

Calculation of spectra

run computes the spectra for all the initial states by the following steps for each trajectory:

1. Propagate the trajectory by ShakeRattle.

2. Propagate the variation and compute the monodromy matrix in the physical space. The Hessian is evaluated at every step using H2SplinePotential.

3. Exclude the steps where the symplectic condition is not satisfied.

4. 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)]
OMEGA_E = 4401.226  # vibrational frequency of bond stretching [cm^-1]
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 = H2SplinePotential()
ENER_REF = CALC.energy(POSI_REF)

# Coordinate for constraints
COORD_CONST = InternalCoordinates(POSI_REF, MASS)
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 ("bond", "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: 3 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]

    # no vibrational kinetic energy for the reference state to target the vibrational ground state (v=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)

    omega_e = OMEGA_E * CmInv / HBar  # convert to internal unit
    gamma_ref = [MASS[0] * omega_e / (2.0 * HBar), # harmonic width (m_H * omega / 2 hbar)
                 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_rovib_{k_rot:.0f}_v0_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, OMEGA_E)
    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 = 2483.7, max symp error = 1.547e-09, (1.2 sec)
    11 / 2000, K_rot = 4946.9, K_vib = 1226.6, max symp error = 3.852e-10, (5.9 sec)
    21 / 2000, K_rot = 9808.5, K_vib = 3522.1, max symp error = 4.342e-10, (10.7 sec)
    31 / 2000, K_rot = 8930.5, K_vib = 2197.7, max symp error = 4.411e-10, (15.3 sec)
    41 / 2000, K_rot = 1592.5, K_vib = 1338.6, max symp error = 1.241e-09, (19.7 sec)
    51 / 2000, K_rot = 7396.6, K_vib = 900.2, max symp error = 1.002e-09, (24.1 sec)
    61 / 2000, K_rot = 2192.0, K_vib = 3063.7, max symp error = 7.122e-10, (28.6 sec)
    71 / 2000, K_rot = 6597.3, K_vib = 3239.4, max symp error = 2.918e-10, (33.2 sec)
    81 / 2000, K_rot = 2769.8, K_vib = 3814.3, max symp error = 1.893e-09, (38.0 sec)
    91 / 2000, K_rot = 2161.3, K_vib = 2645.7, max symp error = 1.017e-09, (42.5 sec)
   101 / 2000, K_rot = 3024.0, K_vib = 1407.6, max symp error = 1.212e-09, (47.0 sec)
   111 / 2000, K_rot = 4798.8, K_vib = 2901.8, max symp error = 1.571e-09, (52.1 sec)
   121 / 2000, K_rot = 6286.0, K_vib = 521.7, max symp error = 3.860e-10, (60.1 sec)
   131 / 2000, K_rot = 5551.3, K_vib = 304.8, max symp error = 5.869e-10, (65.1 sec)
   141 / 2000, K_rot = 4030.5, K_vib = 101.7, max symp error = 1.311e-09, (70.2 sec)
   151 / 2000, K_rot = 1705.5, K_vib = 3088.0, max symp error = 1.029e-09, (75.2 sec)
   161 / 2000, K_rot = 2817.6, K_vib = 3.2, max symp error = 7.516e-11, (80.0 sec)
   171 / 2000, K_rot = 3701.8, K_vib = 1368.5, max symp error = 5.450e-10, (85.0 sec)
   181 / 2000, K_rot = 7288.7, K_vib = 1800.7, max symp error = 3.011e-10, (89.8 sec)
   191 / 2000, K_rot = 6585.2, K_vib = 134.9, max symp error = 5.543e-10, (94.8 sec)
   201 / 2000, K_rot = 7396.7, K_vib = 1439.0, max symp error = 6.702e-10, (99.6 sec)
   211 / 2000, K_rot = 9867.3, K_vib = 2770.9, max symp error = 3.976e-10, (104.8 sec)
   221 / 2000, K_rot = 7455.7, K_vib = 3599.9, max symp error = 5.752e-10, (109.6 sec)
   231 / 2000, K_rot = 3787.1, K_vib = 3608.2, max symp error = 5.457e-10, (114.5 sec)
   241 / 2000, K_rot = 3470.1, K_vib = 4191.9, max symp error = 1.220e-09, (119.3 sec)
   251 / 2000, K_rot = 1747.2, K_vib = 3873.8, max symp error = 3.657e-10, (124.4 sec)
   261 / 2000, K_rot = 9973.6, K_vib = 3839.2, max symp error = 6.849e-10, (129.3 sec)
   271 / 2000, K_rot = 4377.4, K_vib = 2260.9, max symp error = 1.189e-09, (134.5 sec)
   281 / 2000, K_rot = 9515.0, K_vib = 4348.2, max symp error = 4.160e-10, (139.3 sec)
   291 / 2000, K_rot = 6893.6, K_vib = 1516.3, max symp error = 1.276e-09, (143.9 sec)
   301 / 2000, K_rot = 5376.0, K_vib = 4257.9, max symp error = 5.235e-10, (148.3 sec)
   311 / 2000, K_rot = 5321.8, K_vib = 1053.0, max symp error = 5.579e-10, (152.8 sec)
   321 / 2000, K_rot = 9159.7, K_vib = 2362.1, max symp error = 5.578e-10, (157.2 sec)
   331 / 2000, K_rot = 5392.6, K_vib = 3831.9, max symp error = 3.662e-10, (161.7 sec)
   341 / 2000, K_rot = 487.6, K_vib = 2345.0, max symp error = 1.201e-09, (166.2 sec)
   351 / 2000, K_rot = 6850.8, K_vib = 3830.6, max symp error = 7.851e-10, (170.7 sec)
   361 / 2000, K_rot = 219.6, K_vib = 3923.8, max symp error = 5.729e-10, (175.2 sec)
   371 / 2000, K_rot = 7582.7, K_vib = 3769.7, max symp error = 8.757e-10, (179.7 sec)
   381 / 2000, K_rot = 7779.4, K_vib = 1281.0, max symp error = 3.308e-10, (184.3 sec)
   391 / 2000, K_rot = 4940.4, K_vib = 493.9, max symp error = 2.802e-10, (188.7 sec)
   401 / 2000, K_rot = 8660.0, K_vib = 458.3, max symp error = 7.543e-10, (193.1 sec)
   411 / 2000, K_rot = 1841.7, K_vib = 2777.0, max symp error = 1.155e-10, (197.7 sec)
   421 / 2000, K_rot = 2463.3, K_vib = 3043.4, max symp error = 2.053e-10, (202.2 sec)
   431 / 2000, K_rot = 7594.3, K_vib = 676.2, max symp error = 3.350e-10, (206.9 sec)
   441 / 2000, K_rot = 7529.8, K_vib = 1375.6, max symp error = 6.600e-10, (211.7 sec)
   451 / 2000, K_rot = 3684.4, K_vib = 2426.1, max symp error = 1.455e-09, (216.2 sec)
   461 / 2000, K_rot = 390.8, K_vib = 1589.5, max symp error = 2.390e-09, (220.6 sec)
   471 / 2000, K_rot = 1503.1, K_vib = 2211.0, max symp error = 5.834e-10, (225.1 sec)
   481 / 2000, K_rot = 2764.9, K_vib = 717.1, max symp error = 1.334e-09, (229.8 sec)
   491 / 2000, K_rot = 5392.0, K_vib = 2641.8, max symp error = 7.976e-10, (234.5 sec)
   501 / 2000, K_rot = 2453.7, K_vib = 3233.2, max symp error = 4.626e-10, (239.1 sec)
   511 / 2000, K_rot = 9722.6, K_vib = 4054.7, max symp error = 5.191e-10, (243.6 sec)
   521 / 2000, K_rot = 2164.7, K_vib = 2386.7, max symp error = 8.077e-10, (248.5 sec)
   531 / 2000, K_rot = 1284.9, K_vib = 100.8, max symp error = 1.945e-10, (253.1 sec)
   541 / 2000, K_rot = 3925.6, K_vib = 1523.9, max symp error = 7.747e-10, (257.6 sec)
   551 / 2000, K_rot = 9981.7, K_vib = 240.5, max symp error = 4.452e-10, (262.2 sec)
   561 / 2000, K_rot = 5556.3, K_vib = 328.1, max symp error = 1.488e-09, (266.8 sec)
   571 / 2000, K_rot = 7794.8, K_vib = 4344.3, max symp error = 2.749e-10, (271.4 sec)
   581 / 2000, K_rot = 3518.4, K_vib = 4241.0, max symp error = 1.184e-09, (275.9 sec)
   591 / 2000, K_rot = 3835.0, K_vib = 775.0, max symp error = 1.092e-09, (280.5 sec)
   601 / 2000, K_rot = 817.9, K_vib = 1561.8, max symp error = 1.009e-09, (285.1 sec)
   611 / 2000, K_rot = 2152.6, K_vib = 3333.5, max symp error = 1.169e-09, (289.7 sec)
   621 / 2000, K_rot = 2828.7, K_vib = 2157.5, max symp error = 1.876e-10, (294.3 sec)
   631 / 2000, K_rot = 3203.5, K_vib = 994.2, max symp error = 1.074e-09, (298.9 sec)
   641 / 2000, K_rot = 4838.7, K_vib = 4078.5, max symp error = 1.174e-09, (303.4 sec)
   651 / 2000, K_rot = 4176.6, K_vib = 3858.7, max symp error = 2.300e-10, (307.9 sec)
   661 / 2000, K_rot = 2076.9, K_vib = 578.9, max symp error = 3.775e-10, (312.7 sec)
   671 / 2000, K_rot = 7519.9, K_vib = 3621.0, max symp error = 6.933e-10, (317.5 sec)
   681 / 2000, K_rot = 9171.1, K_vib = 2870.4, max symp error = 6.783e-10, (322.2 sec)
   691 / 2000, K_rot = 5976.7, K_vib = 432.7, max symp error = 1.181e-09, (326.7 sec)
   701 / 2000, K_rot = 5310.1, K_vib = 3735.4, max symp error = 6.075e-10, (331.2 sec)
   711 / 2000, K_rot = 7895.5, K_vib = 1044.6, max symp error = 4.220e-10, (335.7 sec)
   721 / 2000, K_rot = 9576.3, K_vib = 786.6, max symp error = 4.420e-10, (340.1 sec)
   731 / 2000, K_rot = 9017.8, K_vib = 1438.3, max symp error = 5.052e-10, (344.5 sec)
   741 / 2000, K_rot = 6847.5, K_vib = 1233.1, max symp error = 4.541e-10, (349.0 sec)
   751 / 2000, K_rot = 6915.3, K_vib = 2984.2, max symp error = 8.296e-10, (353.4 sec)
   761 / 2000, K_rot = 2417.0, K_vib = 614.9, max symp error = 7.797e-10, (357.9 sec)
   771 / 2000, K_rot = 4178.8, K_vib = 1514.2, max symp error = 6.346e-10, (362.4 sec)
   781 / 2000, K_rot = 9218.0, K_vib = 1418.5, max symp error = 6.936e-10, (366.8 sec)
   791 / 2000, K_rot = 2660.2, K_vib = 734.1, max symp error = 3.763e-10, (371.3 sec)
   801 / 2000, K_rot = 2507.7, K_vib = 2022.8, max symp error = 1.487e-09, (375.9 sec)
   811 / 2000, K_rot = 5999.4, K_vib = 2258.1, max symp error = 3.067e-10, (380.6 sec)
   821 / 2000, K_rot = 282.7, K_vib = 1278.3, max symp error = 1.292e-09, (385.5 sec)
   831 / 2000, K_rot = 6042.2, K_vib = 3129.6, max symp error = 9.717e-10, (390.2 sec)
   841 / 2000, K_rot = 3275.4, K_vib = 1746.8, max symp error = 1.292e-09, (394.8 sec)
   851 / 2000, K_rot = 2054.4, K_vib = 997.1, max symp error = 5.412e-10, (399.4 sec)
   861 / 2000, K_rot = 4978.0, K_vib = 3446.3, max symp error = 9.314e-10, (404.2 sec)
   871 / 2000, K_rot = 4786.3, K_vib = 1854.4, max symp error = 4.216e-10, (408.9 sec)
   881 / 2000, K_rot = 3213.5, K_vib = 148.8, max symp error = 1.112e-09, (413.5 sec)
   891 / 2000, K_rot = 1980.6, K_vib = 4212.3, max symp error = 1.628e-09, (418.2 sec)
   901 / 2000, K_rot = 9113.5, K_vib = 1305.4, max symp error = 7.734e-10, (422.9 sec)
   911 / 2000, K_rot = 7299.0, K_vib = 1744.8, max symp error = 8.285e-10, (427.5 sec)
   921 / 2000, K_rot = 4054.4, K_vib = 3002.2, max symp error = 5.457e-10, (432.3 sec)
   931 / 2000, K_rot = 9594.0, K_vib = 652.2, max symp error = 3.325e-10, (437.1 sec)
   941 / 2000, K_rot = 7537.2, K_vib = 2377.7, max symp error = 8.212e-10, (441.9 sec)
   951 / 2000, K_rot = 9288.8, K_vib = 278.3, max symp error = 2.818e-10, (446.6 sec)
   961 / 2000, K_rot = 2782.4, K_vib = 3178.6, max symp error = 1.426e-09, (451.1 sec)
   971 / 2000, K_rot = 7003.2, K_vib = 407.0, max symp error = 2.431e-10, (455.7 sec)
   981 / 2000, K_rot = 2347.4, K_vib = 2065.2, max symp error = 1.072e-09, (460.3 sec)
   991 / 2000, K_rot = 3765.6, K_vib = 1960.8, max symp error = 4.861e-10, (464.8 sec)
  1001 / 2000, K_rot = 9447.7, K_vib = 2569.5, max symp error = 4.796e-10, (469.4 sec)
  1011 / 2000, K_rot = 4098.7, K_vib = 2410.4, max symp error = 5.334e-10, (474.0 sec)
  1021 / 2000, K_rot = 2402.9, K_vib = 3137.9, max symp error = 1.211e-09, (478.7 sec)
  1031 / 2000, K_rot = 7194.7, K_vib = 3809.3, max symp error = 5.450e-10, (483.3 sec)
  1041 / 2000, K_rot = 2557.9, K_vib = 4334.3, max symp error = 4.228e-10, (487.9 sec)
  1051 / 2000, K_rot = 7609.2, K_vib = 1718.9, max symp error = 3.706e-10, (492.5 sec)
  1061 / 2000, K_rot = 4817.1, K_vib = 2678.6, max symp error = 5.726e-10, (497.2 sec)
  1071 / 2000, K_rot = 5458.6, K_vib = 3853.6, max symp error = 1.011e-09, (501.7 sec)
  1081 / 2000, K_rot = 4589.1, K_vib = 2305.1, max symp error = 4.308e-10, (506.2 sec)
  1091 / 2000, K_rot = 5344.1, K_vib = 1818.7, max symp error = 3.584e-10, (510.9 sec)
  1101 / 2000, K_rot = 9410.6, K_vib = 2516.0, max symp error = 7.857e-10, (515.5 sec)
  1111 / 2000, K_rot = 3666.5, K_vib = 3507.4, max symp error = 8.588e-10, (520.0 sec)
  1121 / 2000, K_rot = 240.3, K_vib = 3517.0, max symp error = 6.272e-10, (524.6 sec)
  1131 / 2000, K_rot = 6294.3, K_vib = 3297.4, max symp error = 6.860e-10, (529.2 sec)
  1141 / 2000, K_rot = 4357.6, K_vib = 1741.9, max symp error = 6.422e-10, (533.8 sec)
  1151 / 2000, K_rot = 9056.6, K_vib = 2342.0, max symp error = 5.668e-10, (538.4 sec)
  1161 / 2000, K_rot = 7166.2, K_vib = 60.6, max symp error = 1.908e-10, (542.9 sec)
  1171 / 2000, K_rot = 8916.5, K_vib = 325.7, max symp error = 4.829e-10, (547.6 sec)
  1181 / 2000, K_rot = 451.6, K_vib = 2794.4, max symp error = 1.549e-09, (552.3 sec)
  1191 / 2000, K_rot = 2126.8, K_vib = 777.3, max symp error = 5.189e-10, (557.2 sec)
  1201 / 2000, K_rot = 9789.0, K_vib = 3297.3, max symp error = 6.385e-10, (562.1 sec)
  1211 / 2000, K_rot = 5422.0, K_vib = 3283.8, max symp error = 5.213e-10, (566.7 sec)
  1221 / 2000, K_rot = 3624.6, K_vib = 2727.0, max symp error = 2.683e-10, (571.3 sec)
  1231 / 2000, K_rot = 2296.3, K_vib = 1739.8, max symp error = 5.492e-10, (575.9 sec)
  1241 / 2000, K_rot = 2172.9, K_vib = 880.1, max symp error = 8.172e-10, (580.4 sec)
  1251 / 2000, K_rot = 5639.4, K_vib = 1741.8, max symp error = 5.743e-10, (585.0 sec)
  1261 / 2000, K_rot = 4547.6, K_vib = 1369.5, max symp error = 2.408e-10, (589.5 sec)
  1271 / 2000, K_rot = 350.1, K_vib = 1926.1, max symp error = 8.248e-10, (594.0 sec)
  1281 / 2000, K_rot = 3161.5, K_vib = 743.2, max symp error = 9.574e-10, (598.5 sec)
  1291 / 2000, K_rot = 2615.2, K_vib = 460.4, max symp error = 8.505e-10, (603.1 sec)
  1301 / 2000, K_rot = 3493.6, K_vib = 1263.6, max symp error = 9.196e-10, (607.9 sec)
  1311 / 2000, K_rot = 2293.5, K_vib = 1831.8, max symp error = 7.507e-10, (612.7 sec)
  1321 / 2000, K_rot = 1026.4, K_vib = 635.0, max symp error = 1.644e-09, (617.3 sec)
  1331 / 2000, K_rot = 8248.1, K_vib = 4326.5, max symp error = 3.777e-10, (621.9 sec)
  1341 / 2000, K_rot = 9451.1, K_vib = 4342.0, max symp error = 1.887e-10, (626.5 sec)
  1351 / 2000, K_rot = 6264.2, K_vib = 2723.8, max symp error = 7.169e-10, (631.1 sec)
  1361 / 2000, K_rot = 4614.2, K_vib = 1380.0, max symp error = 4.753e-10, (635.7 sec)
  1371 / 2000, K_rot = 8100.1, K_vib = 343.3, max symp error = 5.250e-10, (640.3 sec)
  1381 / 2000, K_rot = 6091.0, K_vib = 1225.5, max symp error = 1.007e-09, (644.8 sec)
  1391 / 2000, K_rot = 3336.7, K_vib = 3789.3, max symp error = 1.017e-09, (649.4 sec)
  1401 / 2000, K_rot = 9340.1, K_vib = 1438.4, max symp error = 2.340e-10, (654.0 sec)
  1411 / 2000, K_rot = 3087.2, K_vib = 4160.6, max symp error = 6.153e-10, (658.6 sec)
  1421 / 2000, K_rot = 6007.2, K_vib = 3360.5, max symp error = 1.029e-09, (663.5 sec)
  1431 / 2000, K_rot = 8549.9, K_vib = 975.1, max symp error = 4.291e-10, (668.3 sec)
  1441 / 2000, K_rot = 3685.4, K_vib = 2469.4, max symp error = 9.055e-10, (672.9 sec)
  1451 / 2000, K_rot = 1389.7, K_vib = 509.4, max symp error = 3.514e-10, (677.4 sec)
  1461 / 2000, K_rot = 5130.7, K_vib = 3560.5, max symp error = 7.755e-10, (681.9 sec)
  1471 / 2000, K_rot = 866.6, K_vib = 854.5, max symp error = 1.545e-09, (686.3 sec)
  1481 / 2000, K_rot = 7205.9, K_vib = 2784.5, max symp error = 7.075e-10, (690.8 sec)
  1491 / 2000, K_rot = 5494.3, K_vib = 2530.7, max symp error = 5.348e-10, (695.4 sec)
  1501 / 2000, K_rot = 6408.4, K_vib = 111.1, max symp error = 5.867e-10, (700.2 sec)
  1511 / 2000, K_rot = 2938.0, K_vib = 4235.7, max symp error = 2.310e-10, (705.0 sec)
  1521 / 2000, K_rot = 9885.7, K_vib = 1728.6, max symp error = 4.125e-10, (709.5 sec)
  1531 / 2000, K_rot = 5770.4, K_vib = 795.9, max symp error = 4.784e-10, (714.1 sec)
  1541 / 2000, K_rot = 9285.6, K_vib = 1281.8, max symp error = 4.266e-10, (718.9 sec)
  1551 / 2000, K_rot = 1669.5, K_vib = 1074.7, max symp error = 1.146e-09, (723.6 sec)
  1561 / 2000, K_rot = 8796.3, K_vib = 2569.7, max symp error = 6.655e-10, (728.1 sec)
  1571 / 2000, K_rot = 2098.4, K_vib = 776.9, max symp error = 2.883e-09, (732.6 sec)
  1581 / 2000, K_rot = 6521.5, K_vib = 1206.9, max symp error = 2.943e-10, (737.4 sec)
  1591 / 2000, K_rot = 8306.3, K_vib = 3139.3, max symp error = 8.134e-10, (742.0 sec)
  1601 / 2000, K_rot = 5014.2, K_vib = 4385.6, max symp error = 5.634e-10, (746.8 sec)
  1611 / 2000, K_rot = 2860.3, K_vib = 3898.4, max symp error = 1.571e-09, (751.5 sec)
  1621 / 2000, K_rot = 1863.2, K_vib = 3050.5, max symp error = 5.345e-10, (756.5 sec)
  1631 / 2000, K_rot = 6431.5, K_vib = 1493.9, max symp error = 3.303e-10, (761.4 sec)
  1641 / 2000, K_rot = 4589.4, K_vib = 1564.1, max symp error = 1.162e-09, (766.1 sec)
  1651 / 2000, K_rot = 3933.0, K_vib = 1785.4, max symp error = 5.184e-10, (770.8 sec)
  1661 / 2000, K_rot = 9232.2, K_vib = 474.1, max symp error = 3.365e-10, (775.7 sec)
  1671 / 2000, K_rot = 7627.6, K_vib = 2090.2, max symp error = 7.010e-10, (780.6 sec)
  1681 / 2000, K_rot = 6536.5, K_vib = 1957.2, max symp error = 4.526e-10, (785.5 sec)
  1691 / 2000, K_rot = 897.3, K_vib = 1412.1, max symp error = 1.139e-09, (790.3 sec)
  1701 / 2000, K_rot = 1867.7, K_vib = 3648.6, max symp error = 2.817e-10, (795.2 sec)
  1711 / 2000, K_rot = 1623.7, K_vib = 2682.8, max symp error = 1.224e-09, (799.9 sec)
  1721 / 2000, K_rot = 5538.5, K_vib = 1783.7, max symp error = 6.541e-10, (804.5 sec)
  1731 / 2000, K_rot = 5618.3, K_vib = 3934.1, max symp error = 1.223e-09, (809.2 sec)
  1741 / 2000, K_rot = 8806.9, K_vib = 2356.1, max symp error = 4.143e-10, (813.8 sec)
  1751 / 2000, K_rot = 5849.4, K_vib = 3275.8, max symp error = 6.418e-10, (818.3 sec)
  1761 / 2000, K_rot = 192.3, K_vib = 1805.9, max symp error = 6.478e-11, (822.9 sec)
  1771 / 2000, K_rot = 4075.3, K_vib = 1004.7, max symp error = 8.884e-10, (827.5 sec)
  1781 / 2000, K_rot = 5590.7, K_vib = 3830.3, max symp error = 4.418e-10, (832.1 sec)
  1791 / 2000, K_rot = 2198.4, K_vib = 3373.9, max symp error = 1.232e-09, (836.7 sec)
  1801 / 2000, K_rot = 6597.9, K_vib = 3443.8, max symp error = 1.006e-09, (841.2 sec)
  1811 / 2000, K_rot = 4280.0, K_vib = 976.2, max symp error = 9.927e-10, (845.8 sec)
  1821 / 2000, K_rot = 4019.1, K_vib = 807.5, max symp error = 2.038e-09, (850.5 sec)
  1831 / 2000, K_rot = 3659.6, K_vib = 1709.2, max symp error = 4.123e-10, (855.1 sec)
  1841 / 2000, K_rot = 4965.5, K_vib = 2212.8, max symp error = 6.215e-10, (859.6 sec)
  1851 / 2000, K_rot = 2757.4, K_vib = 3108.4, max symp error = 6.709e-10, (864.5 sec)
  1861 / 2000, K_rot = 9182.1, K_vib = 3572.3, max symp error = 2.572e-10, (869.3 sec)
  1871 / 2000, K_rot = 8172.8, K_vib = 4193.5, max symp error = 5.430e-10, (874.1 sec)
  1881 / 2000, K_rot = 3323.9, K_vib = 2253.9, max symp error = 5.853e-10, (879.0 sec)
  1891 / 2000, K_rot = 5787.3, K_vib = 3324.0, max symp error = 1.406e-09, (883.8 sec)
  1901 / 2000, K_rot = 128.8, K_vib = 3748.5, max symp error = 6.703e-10, (888.7 sec)
  1911 / 2000, K_rot = 4775.1, K_vib = 1557.6, max symp error = 8.290e-10, (893.5 sec)
  1921 / 2000, K_rot = 7999.3, K_vib = 1498.7, max symp error = 2.526e-10, (898.2 sec)
  1931 / 2000, K_rot = 650.4, K_vib = 1933.2, max symp error = 1.772e-10, (902.8 sec)
  1941 / 2000, K_rot = 2412.1, K_vib = 1297.5, max symp error = 4.961e-10, (907.6 sec)
  1951 / 2000, K_rot = 9588.3, K_vib = 3472.6, max symp error = 4.916e-10, (912.4 sec)
  1961 / 2000, K_rot = 9360.6, K_vib = 4297.3, max symp error = 1.074e-09, (917.1 sec)
  1971 / 2000, K_rot = 9367.1, K_vib = 3811.7, max symp error = 2.225e-10, (921.7 sec)
  1981 / 2000, K_rot = 9965.5, K_vib = 4154.6, max symp error = 5.324e-10, (926.2 sec)
  1991 / 2000, K_rot = 2330.2, K_vib = 1885.9, max symp error = 2.760e-09, (930.8 sec)

Plot results

The reference levels are the experimental ro-vibrational levels of H2 in the vibrational ground state, given by the Dunham expansion

\[ E_J = E_{\text{ZPE}} + B J(J+1) - D [J(J+1)]^2 + H [J(J+1)]^3 + L [J(J+1)]^4 \]

with the coefficients from D. E. Jennings and J. W. Brault, J. Mol. Spectrosc. 102, 265-272 (1983).

The zero-point energy is set to \(E_{\text{ZPE}}\) = 2183 \(\textrm{cm}^{-1}\) here.

def f(x, B, D, H, L):
    return B*x*(x+1)-D*x*x*(x+1)*(x+1)+H*(x**3)*((x+1)**3)+L*(x**4)*((x+1)**4)

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

    # Dunham expansion coefficients for H2 vibrational ground state
    # D. E. Jennings and J. W. Brault, "The ground state of molecular hydrogen", Journal of Molecular Spectroscopy 102, 265–272 (1983).
    B = 59.334547
    D = 0.0456821
    H = 4.8542e-5
    L = -5.411e-8
    ZPE = 2183.0  # zero-point energy [cm^-1]

    analytical_energies = np.array([f(n, B, D, H, L) + ZPE for n in range(20)])
    analytical_energies = analytical_energies[analytical_energies <= 10000]

    print(f"Analytical 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)
        mask = (ener_grid[peaks] <= 10000)
        peaks = peaks[mask]

        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="reference")
        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()

Analytical energies (cm^-1): [2183.  2301.5 2537.4 2888.5 3351.8 3923.2 4597.9 5370.5 6234.9 7184.9
 8213.7 9314.4]
Peaks (cm^-1) for K = 500: [2289.1 2537.9 2883.4 3341.4 3915.3 4594.2 5362.1]
Peaks (cm^-1) for K = 2000: [2541.7 2882.1 3340.2 3913.4 4592.3 5362.1 6226.1 7174.7 8204.1]
Peaks (cm^-1) for K = 5000: [3338.3 3912.1 4591.  5361.5 6225.4 7174.1 8203.5 9308.6]

Since the peak corresponding to the vibrational zero-point energy is not clearly resolved, we evaluate the successive intervals \(\Delta_J \equiv E_J - E_{J-1}\), which are independent of the choice of the zero-point energy, as in the paper. The computed intervals from the spectrum for \(K\) = 2000 \(\textrm{cm}^{-1}\) are compared with the reference values \(\Delta_J^{(\text{ref})}\) from the Dunham expansion and the rigid-rotor intervals \(\Delta_J^{(\text{rigid})} = 2 B_e J\) (all in \(\textrm{cm}^{-1}\)).

\(J\)	\(\Delta_J\)	\(\Delta_J^{(\text{ref})}\)	\(\Delta_J - \Delta_J^{(\text{ref})}\)	\(\Delta_J^{(\text{rigid})}\)	\(\Delta_J - \Delta_J^{(\text{rigid})}\)
3	340.4	351.1	−10.7	364.9	−24.5
4	458.1	463.3	−5.2	486.6	−28.5
5	573.2	571.4	+1.8	608.2	−35.0
6	678.9	674.7	+4.2	729.8	−50.9
7	769.8	772.6	−2.8	851.5	−81.7
8	864.0	864.5	−0.5	973.1	−109.1
9	948.6	950.0	−1.4	1094.8	−146.2
10	1029.4	1028.8	+0.6	1216.4	−187.0

• The computed intervals reproduce the reference values with an RMSE of about 5 \(\textrm{cm}^{-1}\), without systematic dependence on \(J\).

• In contrast, the computed intervals deviate from the rigid-rotor values increasingly with \(J\), reaching \(-187\) \(\textrm{cm}^{-1}\) at \(J = 10\). The computed spectrum thus follows the Dunham expansion rather than the rigid rotor, showing that the ro-vibrational coupling is properly accounted for.

• Unlike the single-trajectory spectrum (Spectrum Calculation tutorial), which agrees with the reference levels only near the energy of the trajectory, the spectrum with the phase-space sampling reproduces the levels over the whole range.