#!/usr/bin/env python3
# SPDX-License-Identifier: MIT
"""Verified normal-mode and anharmonic vibrational-spectrum benchmarks.

The program has two deliberately small parts:

1. A calibrated, stretch-only, mass-weighted Hessian for linear CO2.
2. A sinc-DVR solution of a Morse model for H35Cl, including transition
   moments from a polynomial dipole surface.

The calculations are pedagogical numerical benchmarks, not ab initio
predictions. NumPy is the only dependency. Atomic units are used internally.
"""

from __future__ import annotations

import argparse
import csv
import json
import math
import platform
from dataclasses import dataclass
from pathlib import Path
from typing import Any

import numpy as np


HARTREE_TO_WAVENUMBER = 219_474.631_363_20
ATOMIC_MASS_UNIT_TO_ELECTRON_MASS = 1_822.888_486_209
BOHR_TO_ANGSTROM = 0.529_177_210_903

MASS_H1_U = 1.007_825_032_23
MASS_C12_U = 12.0
MASS_O16_U = 15.994_914_619_57
MASS_CL35_U = 34.968_852_682

CO2_SYMMETRIC_TARGET_CM = 1_354.0
CO2_ANTISYMMETRIC_TARGET_CM = 2_396.0

HCL_WE_CM = 2_990.946_3
HCL_WEXE_CM = 52.818_6
HCL_WEYE_CM = 0.224_37
HCL_WEZE_CM = -0.012_18

HCL_DIPOLE_REFERENCE_D = 0.0
HCL_DIPOLE_DERIVATIVES = {
    1: 0.925,
    2: 0.16,
    3: -3.83,
    4: -9.3,
}

COMPILED_HCL_BAND_INTENSITIES = {
    1: 130.0,
    2: 2.9,
    3: 0.023,
}


@dataclass(frozen=True)
class CO2Benchmark:
    hessian: np.ndarray
    mass_weighted_hessian: np.ndarray
    frequencies_cm: np.ndarray
    mass_weighted_modes: np.ndarray
    cartesian_modes: np.ndarray
    force_constant_bond: float
    force_constant_outer: float
    rows: list[dict[str, Any]]
    validation: dict[str, Any]


@dataclass(frozen=True)
class MorseParameters:
    reduced_mass: float
    dissociation_energy: float
    range_parameter: float


@dataclass(frozen=True)
class DVRResult:
    grid: np.ndarray
    spacing: float
    potential: np.ndarray
    hamiltonian: np.ndarray
    energies: np.ndarray
    eigenvectors: np.ndarray


def wavenumber_to_angular_frequency(value_cm: float) -> float:
    """Convert a spectroscopic wavenumber to angular frequency in a.u."""
    return value_cm / HARTREE_TO_WAVENUMBER


def co2_benchmark() -> CO2Benchmark:
    """Build and diagonalize a calibrated one-dimensional CO2 Hessian."""
    masses = (
        np.array(
            [MASS_O16_U, MASS_C12_U, MASS_O16_U],
            dtype=float,
        )
        * ATOMIC_MASS_UNIT_TO_ELECTRON_MASS
    )
    mass_o, mass_c, _ = masses
    omega_s = wavenumber_to_angular_frequency(
        CO2_SYMMETRIC_TARGET_CM
    )
    omega_a = wavenumber_to_angular_frequency(
        CO2_ANTISYMMETRIC_TARGET_CM
    )

    force_constant_bond = omega_a**2 / (
        1.0 / mass_o + 2.0 / mass_c
    )
    force_constant_outer = 0.5 * (
        mass_o * omega_s**2 - force_constant_bond
    )

    bond_gradient = np.array(
        [
            [-1.0, 1.0, 0.0],
            [0.0, -1.0, 1.0],
        ]
    )
    outer_gradient = np.array([-1.0, 0.0, 1.0])
    hessian = (
        force_constant_bond * bond_gradient.T @ bond_gradient
        + force_constant_outer
        * np.outer(outer_gradient, outer_gradient)
    )

    inverse_sqrt_mass = np.diag(masses ** -0.5)
    mass_weighted = inverse_sqrt_mass @ hessian @ inverse_sqrt_mass
    eigenvalues, mass_weighted_modes = np.linalg.eigh(mass_weighted)
    eigenvalues = np.maximum(eigenvalues, 0.0)
    eigenvalues[
        eigenvalues < np.max(eigenvalues) * 1.0e-14
    ] = 0.0
    frequencies_cm = (
        np.sqrt(eigenvalues) * HARTREE_TO_WAVENUMBER
    )
    cartesian_modes = inverse_sqrt_mass @ mass_weighted_modes

    for index in range(3):
        vector = cartesian_modes[:, index]
        if index == 0 and float(np.sum(vector)) < 0.0:
            mass_weighted_modes[:, index] *= -1.0
            cartesian_modes[:, index] *= -1.0
        elif index == 1 and vector[2] < 0.0:
            mass_weighted_modes[:, index] *= -1.0
            cartesian_modes[:, index] *= -1.0
        elif index == 2 and float(vector[0] + vector[2]) < 0.0:
            mass_weighted_modes[:, index] *= -1.0
            cartesian_modes[:, index] *= -1.0

    labels = [
        "translation",
        "symmetric stretch",
        "antisymmetric stretch",
    ]
    irreducible_labels = ["Sigma_u+", "Sigma_g+", "Sigma_u+"]
    inversion_labels = ["ungerade", "gerade", "ungerade"]
    ir_labels = ["not a vibration", "forbidden", "allowed"]
    raman_labels = ["not a vibration", "allowed", "forbidden"]

    rows: list[dict[str, Any]] = []
    for index, label in enumerate(labels):
        vector = cartesian_modes[:, index]
        scaled = vector / np.max(np.abs(vector))
        inverted = -vector[::-1]
        parity_expectation = float(
            vector @ inverted / (vector @ vector)
        )
        rows.append(
            {
                "mode_index": index,
                "mode": label,
                "frequency_cm-1": frequencies_cm[index],
                "O_left_displacement": scaled[0],
                "C_displacement": scaled[1],
                "O_right_displacement": scaled[2],
                "inversion_expectation": parity_expectation,
                "inversion_label": inversion_labels[index],
                "linear_molecule_label": irreducible_labels[index],
                "IR_activity_in_stretch_model": ir_labels[index],
                "Raman_activity_in_stretch_model": raman_labels[index],
            }
        )

    analytic_frequencies = np.array(
        [
            0.0,
            math.sqrt(
                (
                    force_constant_bond
                    + 2.0 * force_constant_outer
                )
                / mass_o
            )
            * HARTREE_TO_WAVENUMBER,
            math.sqrt(
                force_constant_bond
                * (1.0 / mass_o + 2.0 / mass_c)
            )
            * HARTREE_TO_WAVENUMBER,
        ]
    )
    centre_of_mass_residuals = [
        abs(float(masses @ cartesian_modes[:, index]))
        for index in (1, 2)
    ]
    validation = {
        "hessian_row_sum_max_Eh_per_bohr2": float(
            np.max(np.abs(np.sum(hessian, axis=1)))
        ),
        "mass_weighted_symmetry_error": float(
            np.max(np.abs(mass_weighted - mass_weighted.T))
        ),
        "orthonormality_error": float(
            np.max(
                np.abs(
                    mass_weighted_modes.T @ mass_weighted_modes
                    - np.eye(3)
                )
            )
        ),
        "analytic_frequency_max_error_cm-1": float(
            np.max(np.abs(frequencies_cm - analytic_frequencies))
        ),
        "vibrational_centre_of_mass_residual_max": max(
            centre_of_mass_residuals
        ),
        "inversion_expectation_max_error": float(
            max(
                abs(abs(row["inversion_expectation"]) - 1.0)
                for row in rows
            )
        ),
    }
    checks = {
        "translation_invariance": (
            validation["hessian_row_sum_max_Eh_per_bohr2"] < 1.0e-14
        ),
        "symmetric_matrix": (
            validation["mass_weighted_symmetry_error"] < 1.0e-14
        ),
        "orthonormal_modes": (
            validation["orthonormality_error"] < 1.0e-12
        ),
        "analytic_frequencies": (
            validation["analytic_frequency_max_error_cm-1"] < 1.0e-7
        ),
        "centre_of_mass_fixed": (
            validation[
                "vibrational_centre_of_mass_residual_max"
            ]
            < 1.0e-10
        ),
        "definite_inversion_parity": (
            validation["inversion_expectation_max_error"] < 1.0e-12
        ),
    }
    validation["checks"] = checks
    if not all(checks.values()):
        raise RuntimeError(f"CO2 validation failed: {validation}")

    return CO2Benchmark(
        hessian=hessian,
        mass_weighted_hessian=mass_weighted,
        frequencies_cm=frequencies_cm,
        mass_weighted_modes=mass_weighted_modes,
        cartesian_modes=cartesian_modes,
        force_constant_bond=force_constant_bond,
        force_constant_outer=force_constant_outer,
        rows=rows,
        validation=validation,
    )


def hcl_morse_parameters() -> MorseParameters:
    """Infer Morse parameters from the H35Cl spectroscopic constants."""
    reduced_mass_u = (
        MASS_H1_U
        * MASS_CL35_U
        / (MASS_H1_U + MASS_CL35_U)
    )
    reduced_mass = (
        reduced_mass_u * ATOMIC_MASS_UNIT_TO_ELECTRON_MASS
    )
    dissociation_energy_cm = (
        HCL_WE_CM**2 / (4.0 * HCL_WEXE_CM)
    )
    dissociation_energy = (
        dissociation_energy_cm / HARTREE_TO_WAVENUMBER
    )
    omega = wavenumber_to_angular_frequency(HCL_WE_CM)
    range_parameter = omega * math.sqrt(
        reduced_mass / (2.0 * dissociation_energy)
    )
    return MorseParameters(
        reduced_mass=reduced_mass,
        dissociation_energy=dissociation_energy,
        range_parameter=range_parameter,
    )


def morse_potential(
    coordinate: np.ndarray,
    parameters: MorseParameters,
) -> np.ndarray:
    """Evaluate V(q) = De [1 - exp(-a q)]^2."""
    return parameters.dissociation_energy * (
        1.0 - np.exp(-parameters.range_parameter * coordinate)
    ) ** 2


def sinc_dvr_kinetic(
    points: int,
    spacing: float,
    reduced_mass: float,
) -> np.ndarray:
    """Return the Colbert-Miller sinc-DVR kinetic-energy matrix."""
    indices = np.arange(points)
    difference = indices[:, None] - indices[None, :]
    kinetic = np.empty((points, points), dtype=float)
    diagonal = difference == 0
    kinetic[diagonal] = (
        math.pi**2 / (6.0 * reduced_mass * spacing**2)
    )
    off_diagonal = ~diagonal
    separation = difference[off_diagonal]
    alternating_sign = np.where(
        separation % 2 == 0,
        1.0,
        -1.0,
    )
    kinetic[off_diagonal] = alternating_sign / (
        reduced_mass * spacing**2 * separation**2
    )
    return kinetic


def solve_morse_dvr(
    points: int,
    coordinate_min: float,
    coordinate_max: float,
    parameters: MorseParameters,
    *,
    eigenvectors: bool = True,
) -> DVRResult:
    """Diagonalize the Morse Hamiltonian on a uniform sinc-DVR grid."""
    if points < 3:
        raise ValueError("points must be at least 3")
    if coordinate_max <= coordinate_min:
        raise ValueError("coordinate-max must exceed coordinate-min")
    grid = np.linspace(coordinate_min, coordinate_max, points)
    spacing = float(grid[1] - grid[0])
    potential = morse_potential(grid, parameters)
    kinetic = sinc_dvr_kinetic(
        points,
        spacing,
        parameters.reduced_mass,
    )
    hamiltonian = kinetic + np.diag(potential)
    if eigenvectors:
        energies, vectors = np.linalg.eigh(hamiltonian)
    else:
        energies = np.linalg.eigvalsh(hamiltonian)
        vectors = np.empty((points, 0))
    return DVRResult(
        grid=grid,
        spacing=spacing,
        potential=potential,
        hamiltonian=hamiltonian,
        energies=energies,
        eigenvectors=vectors,
    )


def morse_exact_term_value(vibrational_quantum_number: int) -> float:
    """Return the exact Morse term value in inverse centimetres."""
    s = vibrational_quantum_number + 0.5
    return HCL_WE_CM * s - HCL_WEXE_CM * s**2


def dunham_term_value(vibrational_quantum_number: int) -> float:
    """Return the H35Cl term value from the selected NIST constants."""
    s = vibrational_quantum_number + 0.5
    return (
        HCL_WE_CM * s
        - HCL_WEXE_CM * s**2
        + HCL_WEYE_CM * s**3
        + HCL_WEZE_CM * s**4
    )


def dipole_surface(coordinate_bohr: np.ndarray) -> np.ndarray:
    """Evaluate the selected HCl dipole polynomial in debye."""
    coordinate_angstrom = coordinate_bohr * BOHR_TO_ANGSTROM
    result = np.full_like(
        coordinate_angstrom,
        HCL_DIPOLE_REFERENCE_D,
        dtype=float,
    )
    for order, derivative in HCL_DIPOLE_DERIVATIVES.items():
        result += (
            derivative
            * coordinate_angstrom**order
            / math.factorial(order)
        )
    return result


def convergence_rows(
    parameters: MorseParameters,
    coordinate_min: float,
    coordinate_max: float,
    final_points: int,
) -> list[dict[str, Any]]:
    """Measure convergence against the analytic Morse energies."""
    point_counts = sorted(
        {
            81,
            101,
            121,
            151,
            final_points,
        }
    )
    rows: list[dict[str, Any]] = []
    exact = np.array(
        [morse_exact_term_value(v) for v in range(8)]
    )
    for points in point_counts:
        result = solve_morse_dvr(
            points,
            coordinate_min,
            coordinate_max,
            parameters,
            eigenvectors=False,
        )
        numerical = (
            result.energies[:8] * HARTREE_TO_WAVENUMBER
        )
        errors = np.abs(numerical - exact)
        rows.append(
            {
                "grid_points": points,
                "q_min_bohr": coordinate_min,
                "q_max_bohr": coordinate_max,
                "spacing_bohr": result.spacing,
                "max_abs_error_v0_to_v7_cm-1": float(
                    np.max(errors)
                ),
                "fundamental_error_cm-1": float(
                    (
                        numerical[1] - numerical[0]
                    )
                    - (
                        exact[1] - exact[0]
                    )
                ),
            }
        )
    return rows


def level_rows(
    result: DVRResult,
    levels: int,
) -> list[dict[str, Any]]:
    """Assemble harmonic, Morse, DVR, and Dunham level comparisons."""
    numerical = result.energies * HARTREE_TO_WAVENUMBER
    rows: list[dict[str, Any]] = []
    for v in range(levels):
        harmonic = HCL_WE_CM * (v + 0.5)
        exact = morse_exact_term_value(v)
        dunham = dunham_term_value(v)
        rows.append(
            {
                "v": v,
                "harmonic_term_cm-1": harmonic,
                "morse_exact_term_cm-1": exact,
                "dvr_term_cm-1": numerical[v],
                "nist_constants_term_cm-1": dunham,
                "dvr_minus_morse_cm-1": numerical[v] - exact,
                "nist_constants_minus_morse_cm-1": dunham - exact,
                "dvr_origin_from_v0_cm-1": (
                    numerical[v] - numerical[0]
                ),
                "nist_origin_from_v0_cm-1": (
                    dunham - dunham_term_value(0)
                ),
            }
        )
    return rows


def transition_rows(
    result: DVRResult,
    upper_levels: int,
) -> list[dict[str, Any]]:
    """Compute v=0 band origins and dipole-intensity proxies."""
    dipole = dipole_surface(result.grid)
    ground_vector = result.eigenvectors[:, 0]
    numerical_cm = result.energies * HARTREE_TO_WAVENUMBER
    raw_rows: list[dict[str, Any]] = []
    for upper in range(1, upper_levels + 1):
        transition_moment = float(
            result.eigenvectors[:, upper]
            @ (dipole * ground_vector)
        )
        dvr_origin = float(
            numerical_cm[upper] - numerical_cm[0]
        )
        intensity_proxy = (
            dvr_origin * abs(transition_moment) ** 2
        )
        compiled_intensity = COMPILED_HCL_BAND_INTENSITIES.get(
            upper
        )
        raw_rows.append(
            {
                "lower_v": 0,
                "upper_v": upper,
                "harmonic_origin_cm-1": upper * HCL_WE_CM,
                "morse_exact_origin_cm-1": (
                    morse_exact_term_value(upper)
                    - morse_exact_term_value(0)
                ),
                "dvr_origin_cm-1": dvr_origin,
                "nist_constants_origin_cm-1": (
                    dunham_term_value(upper)
                    - dunham_term_value(0)
                ),
                "transition_moment_D": abs(transition_moment),
                "intensity_proxy_cm-1_D2": intensity_proxy,
                "compiled_band_intensity_cm-2_atm-1": (
                    compiled_intensity
                ),
            }
        )

    proxy_reference = raw_rows[0]["intensity_proxy_cm-1_D2"]
    compiled_reference = COMPILED_HCL_BAND_INTENSITIES[1]
    for row in raw_rows:
        row["relative_intensity_proxy"] = (
            row["intensity_proxy_cm-1_D2"] / proxy_reference
        )
        compiled = row["compiled_band_intensity_cm-2_atm-1"]
        row["compiled_relative_intensity"] = (
            None
            if compiled is None
            else compiled / compiled_reference
        )
    return raw_rows


def potential_rows(
    result: DVRResult,
    parameters: MorseParameters,
    points: int = 351,
) -> list[dict[str, Any]]:
    """Create plot-ready potential and low-state density data."""
    coordinate = np.linspace(-0.75, 5.0, points)
    potential = morse_potential(coordinate, parameters)
    harmonic = (
        0.5
        * parameters.reduced_mass
        * wavenumber_to_angular_frequency(HCL_WE_CM) ** 2
        * coordinate**2
    )

    rows: list[dict[str, Any]] = []
    for q, morse_value, harmonic_value in zip(
        coordinate,
        potential,
        harmonic,
        strict=True,
    ):
        rows.append(
            {
                "q_bohr": q,
                "q_angstrom": q * BOHR_TO_ANGSTROM,
                "morse_potential_cm-1": (
                    morse_value * HARTREE_TO_WAVENUMBER
                ),
                "harmonic_potential_cm-1": (
                    harmonic_value * HARTREE_TO_WAVENUMBER
                ),
            }
        )
    return rows


def hcl_validation(
    result: DVRResult,
    parameters: MorseParameters,
    convergence: list[dict[str, Any]],
) -> dict[str, Any]:
    """Return validation diagnostics and fail on violated checks."""
    exact = np.array(
        [morse_exact_term_value(v) for v in range(8)]
    )
    numerical = (
        result.energies[:8] * HARTREE_TO_WAVENUMBER
    )
    expected_bound_states = math.floor(
        HCL_WE_CM / (2.0 * HCL_WEXE_CM) - 0.5
    ) + 1
    observed_bound_states = int(
        np.count_nonzero(
            result.energies < parameters.dissociation_energy
        )
    )
    orthonormality_error = float(
        np.max(
            np.abs(
                result.eigenvectors.T @ result.eigenvectors
                - np.eye(result.eigenvectors.shape[1])
            )
        )
    )
    validation: dict[str, Any] = {
        "hamiltonian_symmetry_error_Eh": float(
            np.max(
                np.abs(
                    result.hamiltonian
                    - result.hamiltonian.T
                )
            )
        ),
        "eigenvector_orthonormality_error": orthonormality_error,
        "max_abs_dvr_minus_exact_v0_to_v7_cm-1": float(
            np.max(np.abs(numerical - exact))
        ),
        "fundamental_dvr_minus_exact_cm-1": float(
            (numerical[1] - numerical[0])
            - (exact[1] - exact[0])
        ),
        "expected_bound_states": expected_bound_states,
        "observed_eigenvalues_below_De": observed_bound_states,
        "last_convergence_row": convergence[-1],
    }
    checks = {
        "symmetric_hamiltonian": (
            validation["hamiltonian_symmetry_error_Eh"] < 1.0e-14
        ),
        "orthonormal_eigenvectors": (
            validation["eigenvector_orthonormality_error"] < 1.0e-11
        ),
        "analytic_morse_levels": (
            validation[
                "max_abs_dvr_minus_exact_v0_to_v7_cm-1"
            ]
            < 1.0e-5
        ),
        "analytic_morse_fundamental": (
            abs(
                validation["fundamental_dvr_minus_exact_cm-1"]
            )
            < 1.0e-6
        ),
        "bound_state_count": (
            observed_bound_states == expected_bound_states
        ),
    }
    validation["checks"] = checks
    if not all(checks.values()):
        raise RuntimeError(f"HCl validation failed: {validation}")
    return validation


def format_csv_value(value: Any) -> Any:
    if value is None:
        return ""
    if isinstance(value, (float, np.floating)):
        return f"{float(value):.16g}"
    if isinstance(value, (int, np.integer)):
        return int(value)
    return value


def write_csv(path: Path, rows: list[dict[str, Any]]) -> None:
    if not rows:
        raise ValueError(f"cannot write empty CSV: {path}")
    with path.open("w", newline="", encoding="utf-8") as handle:
        writer = csv.DictWriter(handle, fieldnames=list(rows[0]))
        writer.writeheader()
        for row in rows:
            writer.writerow(
                {
                    key: format_csv_value(value)
                    for key, value in row.items()
                }
            )


def parse_args() -> argparse.Namespace:
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument("--grid-points", type=int, default=201)
    parser.add_argument("--q-min", type=float, default=-1.5)
    parser.add_argument("--q-max", type=float, default=12.0)
    parser.add_argument("--levels", type=int, default=10)
    parser.add_argument("--transitions", type=int, default=5)
    parser.add_argument(
        "--output-dir",
        type=Path,
        default=Path("vibrational-spectra-output"),
    )
    return parser.parse_args()


def main() -> None:
    args = parse_args()
    if args.levels < 8:
        raise ValueError("levels must be at least 8 for validation")
    if args.transitions < 3:
        raise ValueError("transitions must be at least 3")
    if args.transitions >= args.grid_points:
        raise ValueError("transitions must be smaller than grid-points")

    co2 = co2_benchmark()
    parameters = hcl_morse_parameters()
    dvr = solve_morse_dvr(
        args.grid_points,
        args.q_min,
        args.q_max,
        parameters,
    )
    convergence = convergence_rows(
        parameters,
        args.q_min,
        args.q_max,
        args.grid_points,
    )
    levels = level_rows(dvr, args.levels)
    transitions = transition_rows(dvr, args.transitions)
    potential = potential_rows(dvr, parameters)
    hcl_checks = hcl_validation(dvr, parameters, convergence)

    args.output_dir.mkdir(parents=True, exist_ok=True)
    paths = {
        "co2_modes": (
            args.output_dir / "vibrational-co2-normal-modes.csv"
        ),
        "hcl_potential": (
            args.output_dir / "vibrational-h35cl-potential.csv"
        ),
        "hcl_levels": (
            args.output_dir / "vibrational-h35cl-levels.csv"
        ),
        "hcl_transitions": (
            args.output_dir / "vibrational-h35cl-transitions.csv"
        ),
        "hcl_convergence": (
            args.output_dir / "vibrational-h35cl-convergence.csv"
        ),
        "metadata": (
            args.output_dir / "vibrational-spectra-metadata.json"
        ),
    }
    write_csv(paths["co2_modes"], co2.rows)
    write_csv(paths["hcl_potential"], potential)
    write_csv(paths["hcl_levels"], levels)
    write_csv(paths["hcl_transitions"], transitions)
    write_csv(paths["hcl_convergence"], convergence)

    metadata = {
        "scope": {
            "co2": (
                "calibrated one-dimensional stretch-only harmonic model"
            ),
            "h35cl": (
                "spectroscopic Morse model solved by sinc DVR"
            ),
            "claim": (
                "verification and model-error benchmark, not an ab initio "
                "prediction"
            ),
        },
        "constants": {
            "hartree_to_wavenumber_cm-1": HARTREE_TO_WAVENUMBER,
            "atomic_mass_unit_to_electron_mass": (
                ATOMIC_MASS_UNIT_TO_ELECTRON_MASS
            ),
            "bohr_to_angstrom": BOHR_TO_ANGSTROM,
            "isotopic_masses_u": {
                "H1": MASS_H1_U,
                "C12": MASS_C12_U,
                "O16": MASS_O16_U,
                "Cl35": MASS_CL35_U,
            },
        },
        "co2": {
            "target_frequencies_cm-1": {
                "symmetric_stretch": CO2_SYMMETRIC_TARGET_CM,
                "antisymmetric_stretch": (
                    CO2_ANTISYMMETRIC_TARGET_CM
                ),
            },
            "force_constants_Eh_per_bohr2": {
                "bond": co2.force_constant_bond,
                "outer_atom_coupling": co2.force_constant_outer,
            },
            "hessian_Eh_per_bohr2": co2.hessian.tolist(),
            "mass_weighted_hessian": (
                co2.mass_weighted_hessian.tolist()
            ),
            "validation": co2.validation,
            "limitations": [
                "bending modes and rotations are absent",
                "the force field is calibrated rather than predicted",
                "anharmonicity and Fermi resonance are absent",
                "activity labels are symmetry allowances, not intensities",
            ],
        },
        "h35cl": {
            "spectroscopic_constants_cm-1": {
                "omega_e": HCL_WE_CM,
                "omega_e_x_e": HCL_WEXE_CM,
                "omega_e_y_e": HCL_WEYE_CM,
                "omega_e_z_e": HCL_WEZE_CM,
            },
            "morse_parameters": {
                "reduced_mass_electron_masses": (
                    parameters.reduced_mass
                ),
                "dissociation_energy_Eh": (
                    parameters.dissociation_energy
                ),
                "dissociation_energy_cm-1": (
                    parameters.dissociation_energy
                    * HARTREE_TO_WAVENUMBER
                ),
                "range_parameter_per_bohr": (
                    parameters.range_parameter
                ),
            },
            "dipole_surface": {
                "constant_offset_D": (
                    HCL_DIPOLE_REFERENCE_D
                ),
                "constant_offset_note": (
                    "Set to zero because a constant shift does not change "
                    "off-diagonal transition moments between orthogonal "
                    "vibrational eigenstates."
                ),
                "derivatives_D_per_angstrom_power": {
                    str(order): value
                    for order, value
                    in HCL_DIPOLE_DERIVATIVES.items()
                },
            },
            "grid": {
                "points": args.grid_points,
                "q_min_bohr": args.q_min,
                "q_max_bohr": args.q_max,
                "spacing_bohr": dvr.spacing,
            },
            "validation": hcl_checks,
            "limitations": [
                "Morse parameters are fitted from spectroscopic constants",
                "one coordinate and one electronic surface are retained",
                "rotation, temperature, line shape, and isotope mixture are absent",
                "the dipole polynomial is local and intensities are relative proxies",
            ],
        },
        "provenance": {
            "h35cl_constants": {
                "source": "NIST Chemistry WebBook, hydrogen chloride",
                "url": (
                    "https://webbook.nist.gov/cgi/"
                    "cbook.cgi?ID=C7647010&Mask=1000"
                ),
            },
            "sinc_dvr": {
                "authors": "Colbert and Miller",
                "doi": "10.1063/1.462100",
            },
            "hcl_dipole_derivatives": {
                "authors": "Kaiser",
                "journal": (
                    "Journal of Chemical Physics 53, 1686 (1970)"
                ),
                "comment_doi": "10.1063/1.3124083",
            },
        },
        "runtime": {
            "python": platform.python_version(),
            "numpy": np.__version__,
            "platform": platform.platform(),
            "random_seed": None,
        },
        "license": "MIT",
        "outputs": [path.name for path in paths.values()],
    }
    paths["metadata"].write_text(
        json.dumps(metadata, indent=2, sort_keys=True) + "\n",
        encoding="utf-8",
    )

    print(
        "CO2 symmetric stretch     = "
        f"{co2.frequencies_cm[1]:.9f} cm^-1"
    )
    print(
        "CO2 antisymmetric stretch = "
        f"{co2.frequencies_cm[2]:.9f} cm^-1"
    )
    print(
        "H35Cl Morse De            = "
        f"{parameters.dissociation_energy * HARTREE_TO_WAVENUMBER:.9f} "
        "cm^-1"
    )
    print(
        "H35Cl DVR fundamental     = "
        f"{transitions[0]['dvr_origin_cm-1']:.9f} cm^-1"
    )
    print(
        "H35Cl NIST-constant value = "
        f"{transitions[0]['nist_constants_origin_cm-1']:.9f} cm^-1"
    )
    print(
        "max DVR analytic error    = "
        f"{hcl_checks['max_abs_dvr_minus_exact_v0_to_v7_cm-1']:.3e} "
        "cm^-1"
    )
    print(f"outputs                   = {args.output_dir.resolve()}")
    print("validation                = all checks passed")


if __name__ == "__main__":
    main()