8. Rotational Spectra of Methane
The rotational spectrum of methane, a nonlinear molecule, is computed by the phase-space sampling (see also section III B of the paper).
Redundant space |
Potential |
Integrator |
Physical space |
|---|---|---|---|
Quaternion (4 DOFs) with 1 constraint (norm) |
none |
|
Euler angles (3 DOFs) |
Rigid rotor: The rotation is propagated directly by the quaternion, as in the Direct Dynamics tutorial. Since the potential energy is constant for the rigid rotation, no potential is given, and no Hessian evaluation is needed.
Initial states: Coherent states in the Euler angles centered at the reference orientation (\(\psi, \theta, \phi\)) = (0, \(\pi / 2\), 0), with the rotational kinetic energies of 50, 200 and 500 \(\textrm{cm}^{-1}\) and the widths of 0.2 rad.
Sampling: As in the paper, the initial orientation is fixed at the center of the initial states, and only the angular momenta are sampled: the direction is uniform on the unit sphere, and the rotational kinetic energy is uniform.
Coordinate singularities: The trajectories pass through the gimbal lock of the Euler angles (\(\theta = 0, \pi\)) depending on the direction of the angular momentum. The steps where the symplectic condition is not satisfied are excluded.
Reference: Since methane is a spherical top molecule, the rotational levels are given analytically by \(E_J = B J(J+1)\), the same formula as for a linear molecule.
# import necessary libraries
import time
import numpy as np
from pyscivr.basis import get_overlap
from pyscivr.coordinates import QuaternionCoordinates, quaternion_metric
from pyscivr.integrators import ImplicitVerlet
from pyscivr.scivr import get_action, get_lagrangian, get_prefactor, get_ta_spectrum_fft
from pyscivr.units import CmInv, FemtoSec, HBar
from pyscivr.utils import get_geom, get_inertia, get_mass_center
from pyscivr.variations import make_jacobi_func, make_vari_propagator_implicit, propagate_variation_implicit, to_internal
System and coordinates
The system and the coordinates are the same as in the Direct Dynamics tutorial.
The reference structure must be centered at the origin and oriented along the principal axes, so that the moments of inertia can be passed to quaternion_metric.
QuaternionCoordinates defines the Euler angles (yaw, pitch and roll) and the norm condition of the quaternion, where the Euler angles are active. The Jacobian functions are made once and reused for all trajectories.
Sampling of initial conditions
angular_momentum gives the angular momentum in the body frame with the rotational kinetic energy K_rot along the direction L_dir, where \(K_{\text{rot}} = \frac{L^T I^{-1} L}{2}\).
sample_momenta samples the direction uniformly on the unit sphere and K_rot uniformly, and converts the angular momenta into the quaternion momenta by \(p = 2 S(q)^T L_{\text{body}}\).
Calculation of spectra
run computes the spectra for all the initial states by the following steps for each trajectory:
Propagate the trajectory by
ImplicitVerletwithnormalize=True.Propagate the variation and compute the monodromy matrix in the Euler angles.
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
XYZSTR = """
C -0.000000000000 0.000000000000 -0.000000000000
H 0.627615334142 0.627615334142 0.627615334142
H -0.627615334142 -0.627615334142 0.627615334142
H 0.627615334142 -0.627615334142 -0.627615334142
H -0.627615334142 0.627615334142 -0.627615334142
"""
SYMBOLS, POSI_REF, MASS = get_geom(XYZSTR)
assert np.allclose(get_mass_center(POSI_REF, MASS), 0.0) # centered at the origin
INERTIA_TENSOR = get_inertia(POSI_REF, MASS)
assert np.allclose(INERTIA_TENSOR, np.diag(np.diag(INERTIA_TENSOR))) # principal axis frame
INERTIA = np.diag(INERTIA_TENSOR) # (3,) principal moments of inertia
METRIC_FUNC = quaternion_metric(INERTIA)
Q_INIT = np.array([1.0, 0.0, 0.0, 0.0]) # reference orientation (no rotation)
# Coordinates (Euler angles and the norm condition)
COORD_ALL = QuaternionCoordinates()
ACTIVE_IDX = np.array([0, 1, 2]) # yaw, pitch, roll
TOPOLOGIES = [COORD_ALL.topology[i] for i in ACTIVE_IDX] # S1, I-, S1
JACOBI_FUNC, JACOBI_REV_FUNC, INTERNAL_FUNC, MONODROMY_FUNC = make_jacobi_func(COORD_ALL, ACTIVE_IDX)
def get_S(q):
"""omega_body = 2 S(q) dq/dt, p = 2 S(q)^T L_body, L_body = S(q) p / 2 (|q| = 1)."""
q0, q1, q2, q3 = q
return np.array([[-q1, q0, q3, -q2],
[-q2, -q3, q0, q1],
[-q3, q2, -q1, q0]])
def angular_momentum(K_rot, L_dir):
"""Angular momentum in the body frame with the rotational kinetic energy K_rot [cm^-1] along L_dir."""
L_dir = np.atleast_2d(L_dir)
return np.sqrt(2.0 * np.atleast_1d(K_rot)[:, None] * CmInv / np.sum(L_dir**2 / INERTIA, axis=1, keepdims=True)) * L_dir
def sample_momenta(n_sample, K_rot_range, seed=10202):
rng = np.random.default_rng(seed)
L_dir = rng.normal(size=(n_sample, 3))
L_dir /= np.linalg.norm(L_dir, axis=1, keepdims=True) # uniform on the unit sphere
K_rot = rng.uniform(*K_rot_range, size=n_sample)
mome = 2.0 * angular_momentum(K_rot, L_dir) @ get_S(Q_INIT) # (n_sample, 4)
return mome, K_rot
def run(mome_sample_list, K_rot_list, mome_ref_list, gamma_ref_list, fpath_list,
dt, n_step, n_fft, symp_err_tol=1e-2):
# 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 Euler angles
posi_ref_int, mome_ref_int = zip(*[INTERNAL_FUNC(Q_INIT, p) for p in mome_ref_list])
n_sample = len(mome_sample_list)
n_ref = len(mome_ref_list)
spectra = [0.0] * n_ref
vari_propagator = make_vari_propagator_implicit(METRIC_FUNC, dt, normalize=True)
t_start = time.perf_counter()
for i in range(n_sample):
# Get trajectory
traj = ImplicitVerlet(Q_INIT, mome_sample_list[i], dt, METRIC_FUNC, normalize=True).run(n_step)
posi, mome = traj["posi"], traj["mome"]
# Get variation
J_rev_init = JACOBI_REV_FUNC(posi[0], mome[0])
vari = propagate_variation_implicit(vari_propagator, posi, mome, traj["mome_half"], J_rev_init)
posi_int, mome_int, mono_int, symplectic_error = to_internal(MONODROMY_FUNC, posi, mome, vari)
# Points with large symplectic error are treated as singularities
valid = judge_func(symplectic_error)
lag = get_lagrangian(traj["total_energy"], traj["potential_energy"]) # V = 0
action = get_action(lag, dt)
for j in range(n_ref):
gamma = gamma_ref_list[j]
prefactor = get_prefactor(mono_int, gamma, valid=valid)
overlap = get_overlap(TOPOLOGIES, posi_ref_int[j], mome_ref_int[j], posi_int, mome_int, gamma)
ener_grid, 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}, "
f"max symp error = {np.max(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])
Execute calculation
The total propagation time is 2000 fs.
A time step of 1 fs is sufficient for the slow rotation of methane.
if __name__ == "__main__":
# set parameters
n_sample = 2000
dt = 1.0 * FemtoSec
n_step = 2000
n_fft = 2**19
# set reference states
K_rot_ref_list = [50.0, 200.0, 500.0]
L_dir_ref = np.array([1.0, -1.0, -1.0])
mome_ref_list = [2.0 * angular_momentum(k_rot, L_dir_ref)[0] @ get_S(Q_INIT) for k_rot in K_rot_ref_list]
gamma_ref = [1.0 / (2 * 0.2**2)] * 3 # sqrt(1 / 2 gamma) = 0.2 rad for yaw, pitch, roll
gamma_ref_list = np.asarray([gamma_ref] * len(K_rot_ref_list))
fpath_list = [f"methane_pure_rot_{k_rot:.0f}_n_{n_sample}_spectra.npz" for k_rot in K_rot_ref_list]
# sample the initial momenta
mome_sample_list, K_rot_list = sample_momenta(n_sample, K_rot_range=(10.0, 1000.0))
run(mome_sample_list, K_rot_list, mome_ref_list, gamma_ref_list, fpath_list, dt, n_step, n_fft)
1 / 2000, K_rot = 439.3, max symp error = 5.819e-13, (1.1 sec)
11 / 2000, K_rot = 67.9, max symp error = 5.046e-13, (1.9 sec)
21 / 2000, K_rot = 336.4, max symp error = 9.443e-13, (2.6 sec)
31 / 2000, K_rot = 739.7, max symp error = 1.821e-12, (3.3 sec)
41 / 2000, K_rot = 307.7, max symp error = 1.533e-12, (4.0 sec)
51 / 2000, K_rot = 890.7, max symp error = 1.272e-12, (4.7 sec)
61 / 2000, K_rot = 367.5, max symp error = 7.165e-13, (5.4 sec)
71 / 2000, K_rot = 590.6, max symp error = 8.061e-13, (6.1 sec)
81 / 2000, K_rot = 822.9, max symp error = 8.603e-13, (6.8 sec)
91 / 2000, K_rot = 676.4, max symp error = 1.498e-12, (7.5 sec)
101 / 2000, K_rot = 964.9, max symp error = 1.716e-12, (8.2 sec)
111 / 2000, K_rot = 46.4, max symp error = 2.675e-13, (8.9 sec)
121 / 2000, K_rot = 169.4, max symp error = 1.049e-12, (9.6 sec)
131 / 2000, K_rot = 216.6, max symp error = 9.727e-13, (10.3 sec)
141 / 2000, K_rot = 889.7, max symp error = 2.480e-12, (11.0 sec)
151 / 2000, K_rot = 564.2, max symp error = 1.064e-10, (11.7 sec)
161 / 2000, K_rot = 26.2, max symp error = 1.954e-13, (12.4 sec)
171 / 2000, K_rot = 955.0, max symp error = 2.483e-12, (13.2 sec)
181 / 2000, K_rot = 725.0, max symp error = 4.058e-11, (13.9 sec)
191 / 2000, K_rot = 19.3, max symp error = 1.167e-13, (14.6 sec)
201 / 2000, K_rot = 55.8, max symp error = 3.434e-13, (15.3 sec)
211 / 2000, K_rot = 662.7, max symp error = 2.407e-12, (16.1 sec)
221 / 2000, K_rot = 944.3, max symp error = 1.750e-10, (16.8 sec)
231 / 2000, K_rot = 480.9, max symp error = 9.874e-13, (17.5 sec)
241 / 2000, K_rot = 421.0, max symp error = 2.400e-12, (18.2 sec)
251 / 2000, K_rot = 190.7, max symp error = 5.120e-13, (18.9 sec)
261 / 2000, K_rot = 134.6, max symp error = 1.481e-13, (19.6 sec)
271 / 2000, K_rot = 410.8, max symp error = 7.184e-12, (20.3 sec)
281 / 2000, K_rot = 492.0, max symp error = 1.124e-12, (21.0 sec)
291 / 2000, K_rot = 674.8, max symp error = 6.907e-13, (21.7 sec)
301 / 2000, K_rot = 516.9, max symp error = 5.276e-13, (22.5 sec)
311 / 2000, K_rot = 598.9, max symp error = 1.916e-12, (23.2 sec)
321 / 2000, K_rot = 317.9, max symp error = 3.925e-13, (23.9 sec)
331 / 2000, K_rot = 459.9, max symp error = 1.296e-12, (24.6 sec)
341 / 2000, K_rot = 150.7, max symp error = 5.621e-13, (25.3 sec)
351 / 2000, K_rot = 300.8, max symp error = 1.211e-12, (26.0 sec)
361 / 2000, K_rot = 265.6, max symp error = 5.056e-13, (26.7 sec)
371 / 2000, K_rot = 974.8, max symp error = 1.288e-12, (27.4 sec)
381 / 2000, K_rot = 918.6, max symp error = 8.231e-13, (28.2 sec)
391 / 2000, K_rot = 864.9, max symp error = 2.891e-12, (28.9 sec)
401 / 2000, K_rot = 592.6, max symp error = 2.066e-12, (29.6 sec)
411 / 2000, K_rot = 973.3, max symp error = 1.934e-12, (30.3 sec)
421 / 2000, K_rot = 704.8, max symp error = 1.541e-12, (31.0 sec)
431 / 2000, K_rot = 884.7, max symp error = 2.846e-12, (31.8 sec)
441 / 2000, K_rot = 566.6, max symp error = 3.688e-07, (32.5 sec)
451 / 2000, K_rot = 412.0, max symp error = 1.637e-12, (33.2 sec)
461 / 2000, K_rot = 413.9, max symp error = 9.813e-13, (33.9 sec)
471 / 2000, K_rot = 349.1, max symp error = 3.858e-13, (34.6 sec)
481 / 2000, K_rot = 321.3, max symp error = 1.112e-11, (35.3 sec)
491 / 2000, K_rot = 392.0, max symp error = 1.092e-12, (36.1 sec)
501 / 2000, K_rot = 223.7, max symp error = 1.262e-12, (36.8 sec)
511 / 2000, K_rot = 735.3, max symp error = 6.344e-13, (37.4 sec)
521 / 2000, K_rot = 270.9, max symp error = 1.805e-12, (38.1 sec)
531 / 2000, K_rot = 956.4, max symp error = 1.850e-12, (38.8 sec)
541 / 2000, K_rot = 493.2, max symp error = 1.146e-11, (39.4 sec)
551 / 2000, K_rot = 40.8, max symp error = 4.507e-13, (40.1 sec)
561 / 2000, K_rot = 252.0, max symp error = 4.614e-12, (40.8 sec)
571 / 2000, K_rot = 676.3, max symp error = 1.138e-12, (41.4 sec)
581 / 2000, K_rot = 422.3, max symp error = 1.358e-12, (42.1 sec)
591 / 2000, K_rot = 23.1, max symp error = 3.315e-13, (42.8 sec)
601 / 2000, K_rot = 396.8, max symp error = 7.155e-13, (43.6 sec)
611 / 2000, K_rot = 531.2, max symp error = 2.885e-10, (44.3 sec)
621 / 2000, K_rot = 853.8, max symp error = 2.194e-12, (45.1 sec)
631 / 2000, K_rot = 550.3, max symp error = 1.615e-12, (45.9 sec)
641 / 2000, K_rot = 244.1, max symp error = 2.024e-12, (46.6 sec)
651 / 2000, K_rot = 638.0, max symp error = 1.224e-12, (47.3 sec)
661 / 2000, K_rot = 923.9, max symp error = 8.687e-13, (48.0 sec)
671 / 2000, K_rot = 309.8, max symp error = 4.636e-13, (48.7 sec)
681 / 2000, K_rot = 527.0, max symp error = 1.904e-12, (49.4 sec)
691 / 2000, K_rot = 345.0, max symp error = 9.063e-13, (50.1 sec)
701 / 2000, K_rot = 28.2, max symp error = 1.209e-13, (50.8 sec)
711 / 2000, K_rot = 46.5, max symp error = 5.237e-13, (51.5 sec)
721 / 2000, K_rot = 174.9, max symp error = 5.636e-13, (52.2 sec)
731 / 2000, K_rot = 298.7, max symp error = 1.559e-12, (52.9 sec)
741 / 2000, K_rot = 61.9, max symp error = 2.521e-13, (53.6 sec)
751 / 2000, K_rot = 405.6, max symp error = 9.015e-13, (54.3 sec)
761 / 2000, K_rot = 359.3, max symp error = 1.092e-12, (55.0 sec)
771 / 2000, K_rot = 631.8, max symp error = 7.547e-13, (55.7 sec)
781 / 2000, K_rot = 926.7, max symp error = 1.986e-12, (56.4 sec)
791 / 2000, K_rot = 152.1, max symp error = 9.646e-13, (57.2 sec)
801 / 2000, K_rot = 888.9, max symp error = 1.833e-12, (57.9 sec)
811 / 2000, K_rot = 961.4, max symp error = 1.823e-12, (58.6 sec)
821 / 2000, K_rot = 739.9, max symp error = 1.502e-12, (59.4 sec)
831 / 2000, K_rot = 550.3, max symp error = 1.277e-12, (60.1 sec)
841 / 2000, K_rot = 28.5, max symp error = 3.684e-12, (60.8 sec)
851 / 2000, K_rot = 436.9, max symp error = 1.659e-12, (61.5 sec)
861 / 2000, K_rot = 504.7, max symp error = 1.304e-12, (62.2 sec)
871 / 2000, K_rot = 159.4, max symp error = 9.337e-13, (62.9 sec)
881 / 2000, K_rot = 136.8, max symp error = 3.470e-13, (63.6 sec)
891 / 2000, K_rot = 559.4, max symp error = 1.319e-12, (64.3 sec)
901 / 2000, K_rot = 696.9, max symp error = 4.045e-12, (65.1 sec)
911 / 2000, K_rot = 963.7, max symp error = 9.500e-13, (65.8 sec)
921 / 2000, K_rot = 674.8, max symp error = 1.039e-12, (66.5 sec)
931 / 2000, K_rot = 794.3, max symp error = 1.931e-12, (67.3 sec)
941 / 2000, K_rot = 193.7, max symp error = 1.848e-07, (68.0 sec)
951 / 2000, K_rot = 584.7, max symp error = 3.891e-08, (68.7 sec)
961 / 2000, K_rot = 588.5, max symp error = 2.113e-12, (69.4 sec)
971 / 2000, K_rot = 788.1, max symp error = 1.566e-12, (70.2 sec)
981 / 2000, K_rot = 254.7, max symp error = 1.198e-12, (71.0 sec)
991 / 2000, K_rot = 25.3, max symp error = 9.202e-13, (71.7 sec)
1001 / 2000, K_rot = 342.3, max symp error = 8.701e-13, (72.5 sec)
1011 / 2000, K_rot = 892.1, max symp error = 1.037e-12, (73.2 sec)
1021 / 2000, K_rot = 118.1, max symp error = 2.509e-12, (74.0 sec)
1031 / 2000, K_rot = 385.1, max symp error = 1.127e-12, (74.7 sec)
1041 / 2000, K_rot = 652.7, max symp error = 7.690e-13, (75.4 sec)
1051 / 2000, K_rot = 493.2, max symp error = 1.318e-12, (76.2 sec)
1061 / 2000, K_rot = 400.3, max symp error = 8.285e-13, (76.9 sec)
1071 / 2000, K_rot = 688.1, max symp error = 1.435e-12, (77.6 sec)
1081 / 2000, K_rot = 915.9, max symp error = 2.223e-12, (78.3 sec)
1091 / 2000, K_rot = 251.1, max symp error = 1.153e-12, (79.0 sec)
1101 / 2000, K_rot = 866.7, max symp error = 7.419e-13, (79.7 sec)
1111 / 2000, K_rot = 193.5, max symp error = 1.151e-12, (80.4 sec)
1121 / 2000, K_rot = 493.7, max symp error = 3.579e-13, (81.2 sec)
1131 / 2000, K_rot = 602.3, max symp error = 7.773e-13, (81.9 sec)
1141 / 2000, K_rot = 320.5, max symp error = 3.072e-13, (82.6 sec)
1151 / 2000, K_rot = 82.5, max symp error = 6.062e-12, (83.3 sec)
1161 / 2000, K_rot = 468.3, max symp error = 1.050e-12, (84.0 sec)
1171 / 2000, K_rot = 772.7, max symp error = 1.313e-12, (84.7 sec)
1181 / 2000, K_rot = 677.1, max symp error = 1.576e-12, (85.4 sec)
1191 / 2000, K_rot = 406.5, max symp error = 5.472e-13, (86.1 sec)
1201 / 2000, K_rot = 390.3, max symp error = 4.237e-10, (86.9 sec)
1211 / 2000, K_rot = 592.9, max symp error = 1.819e-11, (87.6 sec)
1221 / 2000, K_rot = 826.9, max symp error = 1.264e-12, (88.3 sec)
1231 / 2000, K_rot = 711.1, max symp error = 2.254e-12, (89.0 sec)
1241 / 2000, K_rot = 596.0, max symp error = 3.994e-12, (89.7 sec)
1251 / 2000, K_rot = 178.6, max symp error = 1.025e-12, (90.4 sec)
1261 / 2000, K_rot = 428.3, max symp error = 6.344e-13, (91.1 sec)
1271 / 2000, K_rot = 548.2, max symp error = 2.137e-12, (91.8 sec)
1281 / 2000, K_rot = 630.7, max symp error = 8.693e-13, (92.5 sec)
1291 / 2000, K_rot = 124.0, max symp error = 3.351e-13, (93.2 sec)
1301 / 2000, K_rot = 238.2, max symp error = 1.339e-12, (93.9 sec)
1311 / 2000, K_rot = 809.4, max symp error = 4.970e-13, (94.6 sec)
1321 / 2000, K_rot = 777.6, max symp error = 1.594e-12, (95.3 sec)
1331 / 2000, K_rot = 543.0, max symp error = 1.616e-10, (96.0 sec)
1341 / 2000, K_rot = 935.0, max symp error = 1.266e-12, (96.7 sec)
1351 / 2000, K_rot = 914.0, max symp error = 1.710e-12, (97.4 sec)
1361 / 2000, K_rot = 522.2, max symp error = 4.010e-12, (98.2 sec)
1371 / 2000, K_rot = 479.7, max symp error = 1.568e-12, (98.9 sec)
1381 / 2000, K_rot = 166.4, max symp error = 1.385e-11, (99.6 sec)
1391 / 2000, K_rot = 421.4, max symp error = 9.244e-13, (100.3 sec)
1401 / 2000, K_rot = 389.9, max symp error = 1.806e-12, (101.0 sec)
1411 / 2000, K_rot = 727.5, max symp error = 1.458e-12, (101.8 sec)
1421 / 2000, K_rot = 538.9, max symp error = 1.274e-12, (102.5 sec)
1431 / 2000, K_rot = 756.7, max symp error = 1.519e-12, (103.2 sec)
1441 / 2000, K_rot = 204.9, max symp error = 1.366e-12, (103.9 sec)
1451 / 2000, K_rot = 829.9, max symp error = 5.990e-12, (104.6 sec)
1461 / 2000, K_rot = 758.5, max symp error = 1.374e-12, (105.3 sec)
1471 / 2000, K_rot = 951.8, max symp error = 1.421e-12, (106.0 sec)
1481 / 2000, K_rot = 503.8, max symp error = 1.353e-12, (106.7 sec)
1491 / 2000, K_rot = 508.0, max symp error = 1.102e-12, (107.4 sec)
1501 / 2000, K_rot = 63.6, max symp error = 2.597e-13, (108.2 sec)
1511 / 2000, K_rot = 213.3, max symp error = 9.822e-13, (108.9 sec)
1521 / 2000, K_rot = 316.3, max symp error = 1.023e-12, (109.6 sec)
1531 / 2000, K_rot = 834.3, max symp error = 5.648e-12, (110.3 sec)
1541 / 2000, K_rot = 634.7, max symp error = 1.911e-12, (111.0 sec)
1551 / 2000, K_rot = 931.7, max symp error = 3.078e-12, (111.8 sec)
1561 / 2000, K_rot = 752.3, max symp error = 1.758e-12, (112.5 sec)
1571 / 2000, K_rot = 941.9, max symp error = 1.179e-12, (113.2 sec)
1581 / 2000, K_rot = 425.1, max symp error = 1.168e-11, (113.9 sec)
1591 / 2000, K_rot = 112.6, max symp error = 3.246e-13, (114.6 sec)
1601 / 2000, K_rot = 325.2, max symp error = 1.210e-11, (115.3 sec)
1611 / 2000, K_rot = 768.0, max symp error = 9.784e-13, (116.0 sec)
1621 / 2000, K_rot = 654.5, max symp error = 1.152e-12, (116.7 sec)
1631 / 2000, K_rot = 802.5, max symp error = 1.976e-12, (117.5 sec)
1641 / 2000, K_rot = 927.4, max symp error = 6.437e-13, (118.2 sec)
1651 / 2000, K_rot = 233.0, max symp error = 5.927e-13, (118.9 sec)
1661 / 2000, K_rot = 633.1, max symp error = 1.559e-12, (119.6 sec)
1671 / 2000, K_rot = 143.9, max symp error = 6.444e-13, (120.3 sec)
1681 / 2000, K_rot = 292.2, max symp error = 5.467e-13, (121.0 sec)
1691 / 2000, K_rot = 30.5, max symp error = 4.867e-13, (121.8 sec)
1701 / 2000, K_rot = 55.5, max symp error = 9.172e-13, (122.5 sec)
1711 / 2000, K_rot = 41.2, max symp error = 1.004e-13, (123.2 sec)
1721 / 2000, K_rot = 406.7, max symp error = 9.606e-13, (123.9 sec)
1731 / 2000, K_rot = 586.1, max symp error = 2.692e-12, (124.6 sec)
1741 / 2000, K_rot = 499.9, max symp error = 1.966e-12, (125.3 sec)
1751 / 2000, K_rot = 731.2, max symp error = 3.186e-12, (126.0 sec)
1761 / 2000, K_rot = 10.8, max symp error = 1.192e-11, (126.7 sec)
1771 / 2000, K_rot = 722.8, max symp error = 6.450e-13, (127.5 sec)
1781 / 2000, K_rot = 958.6, max symp error = 2.810e-12, (128.2 sec)
1791 / 2000, K_rot = 626.4, max symp error = 9.325e-13, (128.9 sec)
1801 / 2000, K_rot = 71.5, max symp error = 4.262e-13, (129.6 sec)
1811 / 2000, K_rot = 352.1, max symp error = 5.041e-13, (130.4 sec)
1821 / 2000, K_rot = 380.1, max symp error = 1.133e-12, (131.1 sec)
1831 / 2000, K_rot = 818.7, max symp error = 1.577e-12, (131.8 sec)
1841 / 2000, K_rot = 72.0, max symp error = 1.004e-12, (132.5 sec)
1851 / 2000, K_rot = 112.5, max symp error = 4.005e-13, (133.2 sec)
1861 / 2000, K_rot = 548.9, max symp error = 1.236e-12, (133.9 sec)
1871 / 2000, K_rot = 256.5, max symp error = 1.503e-12, (134.6 sec)
1881 / 2000, K_rot = 166.2, max symp error = 4.501e-13, (135.4 sec)
1891 / 2000, K_rot = 67.2, max symp error = 5.416e-13, (136.1 sec)
1901 / 2000, K_rot = 360.0, max symp error = 6.761e-13, (136.8 sec)
1911 / 2000, K_rot = 745.4, max symp error = 1.371e-12, (137.5 sec)
1921 / 2000, K_rot = 595.9, max symp error = 2.397e-12, (138.2 sec)
1931 / 2000, K_rot = 335.7, max symp error = 2.097e-10, (138.9 sec)
1941 / 2000, K_rot = 400.4, max symp error = 3.903e-13, (139.7 sec)
1951 / 2000, K_rot = 141.0, max symp error = 3.830e-13, (140.5 sec)
1961 / 2000, K_rot = 164.7, max symp error = 2.625e-13, (141.2 sec)
1971 / 2000, K_rot = 169.2, max symp error = 5.982e-13, (142.0 sec)
1981 / 2000, K_rot = 264.3, max symp error = 9.873e-13, (142.7 sec)
1991 / 2000, K_rot = 801.4, max symp error = 7.671e-13, (143.4 sec)
Plot results
Since methane is a spherical top with the moment of inertia \(I\), the rotational levels are given analytically by \(E_J = B J(J+1)\) with \(B = \frac{\hbar^2}{2I}\).
if __name__ == "__main__":
import matplotlib.pyplot as plt
from scipy.signal import find_peaks
B = HBar**2 / (2.0 * INERTIA[0]) / CmInv
analytical_energies = np.array([B * J * (J + 1) for J in range(20)])
analytical_energies = analytical_energies[analytical_energies <= 1000]
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[ener_grid <= 1000])
peaks, _ = find_peaks(spectrum, height=0.1, prominence=0.1)
peaks = peaks[ener_grid[peaks] <= 1000]
print(f"Peaks (cm^-1) for K = {K_rot_ref_list[i]:.0f}: {np.round(ener_grid[peaks], 1)}")
ax[i].plot(ener_grid, spectrum, label=f"K = {K_rot_ref_list[i]:.0f}", color=colors[i])
ax[i].plot(ener_grid[peaks], spectrum[peaks], "x", color="black")
ax[i].vlines(analytical_energies, 0.0, 1.2, color="gray", linestyle="--", label="analytical")
ax[i].set_ylim(0, 1.2)
ax[i].legend()
ax[-1].set_xlim(0, 1000)
ax[-1].set_xlabel("Energy (cm$^{-1}$)")
ax[0].set_ylabel("Intensity (arb. unit)")
plt.tight_layout()
plt.show()
Rotational constant B = 5.307 cm^-1
Analytical rotational energies (cm^-1): [ 0. 10.6 31.8 63.7 106.1 159.2 222.9 297.2 382.1 477.6 583.8 700.5
827.9 965.9]
Peaks (cm^-1) for K = 50: [ 30.7 63.3 104.7 158.5 220.6 295.4]
Peaks (cm^-1) for K = 200: [ 30.5 62.5 104.5 158. 219.8 294.3 382.2 475.6 581.6]
Peaks (cm^-1) for K = 500: [104.6 157.1 219.4 293.5 381.6 475.3 580.9 697.9 826.6 965.7]
The computed peak positions for \(K = 500\) \(\textrm{cm}^{-1}\) are compared with the analytical eigenenergies \(E_J = B J(J+1)\) with \(B = 5.307\) \(\textrm{cm}^{-1}\) (all in \(\textrm{cm}^{-1}\)).
\(J\) |
Computed |
\(B J(J+1)\) |
Deviation |
|---|---|---|---|
4 |
104.6 |
106.1 |
−1.5 |
5 |
157.1 |
159.2 |
−2.1 |
6 |
219.4 |
222.9 |
−3.5 |
7 |
293.5 |
297.2 |
−3.7 |
8 |
381.6 |
382.1 |
−0.5 |
9 |
475.3 |
477.6 |
−2.3 |
10 |
580.9 |
583.8 |
−2.9 |
11 |
697.9 |
700.5 |
−2.6 |
12 |
826.6 |
827.9 |
−1.3 |
13 |
965.7 |
965.9 |
−0.2 |
Quantitative agreement: The computed peaks agree with the analytical eigenenergies with an RMSE of about 2.3 \(\textrm{cm}^{-1}\), well below the intrinsic resolution determined by the propagation time. The same peak positions are obtained for all the initial states, up to \(J = 13\) (about 970 \(\textrm{cm}^{-1}\)) for \(K = 500\) \(\textrm{cm}^{-1}\). The rotational spectrum of the nonlinear molecule is thus reproduced by the quaternion-based propagation, even though the trajectories pass through the gimbal lock of the Euler angles.
Line shape: As in the case of H2, 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 105, 158 and 475 \(\textrm{cm}^{-1}\) for the initial states with \(K\) = 50, 200 and 500 \(\textrm{cm}^{-1}\), respectively.