diff --git a/docs/cli.rst b/docs/cli.rst
new file mode 100644
index 0000000..3ae7f8e
--- /dev/null
+++ b/docs/cli.rst
@@ -0,0 +1,270 @@
+======================
+Command Line Interface
+======================
+
+``nonrad`` includes a comprehensive Command Line Interface (CLI) built with `Click <https://click.palletsprojects.com/>`_. This tool allows you to perform CCD calculations, compute capture coefficients, scale coefficients, and extract electron-phonon coupling parameters directly from your terminal.
+
+Usage
+=====
+
+The entry point for all commands is ``nonrad``. You can view the list of commands and general help by running:
+
+.. code-block:: sh
+
+    nonrad --help
+
+To view the help and option details for a specific command, use:
+
+.. code-block:: sh
+
+    nonrad <command> --help
+
+CCD Commands
+============
+
+These commands are used to prepare and process configuration-coordinate diagrams (CCDs).
+
+prep-ccd
+--------
+Prepare displaced structures along the configuration coordinate.
+
+**Usage:**
+
+.. code-block:: sh
+
+    nonrad prep-ccd [OPTIONS] GROUND_PATH EXCITED_PATH CC_DIR
+
+**Arguments:**
+* ``GROUND_PATH``: Path to the directory containing the ground-state calculation (with ``CONTCAR``).
+* ``EXCITED_PATH``: Path to the directory containing the excited-state calculation (with ``CONTCAR``).
+* ``CC_DIR``: Path to the directory where displaced structures should be written.
+
+**Options:**
+* ``-d, --displace <float float float>``: Displacement range specified as ``(min, max, count)``. (Default: ``-0.5 0.5 9``).
+* ``-f, --input-files TEXT``: Comma-separated list of VASP input files to copy into each displacement folder. (Default: ``KPOINTS,POTCAR,INCAR,job_script.sh``).
+
+pes
+---
+Extract potential energy surfaces (PES) from CCD calculations.
+
+**Usage:**
+
+.. code-block:: sh
+
+    nonrad pes [OPTIONS] CC_DIR GROUND_FILES EXCITED_FILES
+
+**Arguments:**
+* ``CC_DIR``: Root directory containing the CCD displacement folders.
+* ``GROUND_FILES``: Path to the ground-state equilibrium calculation directory (containing ``CONTCAR`` and ``vasprun.xml``).
+* ``EXCITED_FILES``: Path to the excited-state equilibrium calculation directory (containing ``CONTCAR`` and ``vasprun.xml``).
+
+**Options:**
+* ``-e, --energy-diff FLOAT``: [Required] Energy difference *dE* between ground and excited states (in eV).
+* ``-p, --plot``: Flag to plot the potential energy surfaces and save the figure.
+* ``-n, --plot-name TEXT``: Filename for the saved plot. (Default: ``pes.png``).
+
+dq
+--
+Compute the configuration-coordinate displacement *dQ* between two structures.
+
+**Usage:**
+
+.. code-block:: sh
+
+    nonrad dq GROUND_PATH EXCITED_PATH
+
+**Arguments:**
+* ``GROUND_PATH``: Path to the ground-state structure file (e.g. ``POSCAR``, ``CONTCAR``).
+* ``EXCITED_PATH``: Path to the excited-state structure file (e.g. ``POSCAR``, ``CONTCAR``).
+
+q-from-struct
+-------------
+Compute the *Q* value of an arbitrary structure relative to ground and excited endpoint structures.
+
+**Usage:**
+
+.. code-block:: sh
+
+    nonrad q-from-struct [OPTIONS] GROUND_PATH EXCITED_PATH STRUCT_PATH
+
+**Arguments:**
+* ``GROUND_PATH``: Path to the ground-state structure file.
+* ``EXCITED_PATH``: Path to the excited-state structure file.
+* ``STRUCT_PATH``: Path to the structure file to evaluate.
+
+**Options:**
+* ``-t, --tol FLOAT``: Distance tolerance for filtering coordinate differences. (Default: ``1e-4``).
+* ``--nround INTEGER``: Number of decimal places to round *Q*. (Default: ``5``).
+
+barrier
+-------
+Compute the classical barrier height using the harmonic approximation.
+
+**Usage:**
+
+.. code-block:: sh
+
+    nonrad barrier [OPTIONS]
+
+**Options (All Required):**
+* ``--dq FLOAT``: Displacement *dQ* (in amu\ :sup:`1/2` Å).
+* ``--de FLOAT``: Energy offset *dE* between states (in eV).
+* ``--wi FLOAT``: Initial-state frequency (in eV).
+* ``--wf FLOAT``: Final-state frequency (in eV).
+
+
+Capture Commands
+================
+
+capture
+-------
+Evaluate the nonradiative capture coefficient *C* (in cm\ :sup:`3` s\ :sup:`-1`).
+
+**Usage:**
+
+.. code-block:: sh
+
+    nonrad capture [OPTIONS]
+
+**Options:**
+* ``--dq FLOAT``: [Required] Configuration displacement *dQ* (in amu\ :sup:`1/2` Å).
+* ``--de FLOAT``: [Required] Energy offset *dE* (in eV).
+* ``--wi FLOAT``: [Required] Initial-state frequency (in eV).
+* ``--wf FLOAT``: [Required] Final-state frequency (in eV).
+* ``--wif FLOAT``: [Required] Electron-phonon coupling matrix element *Wif* (in eV amu\ :sup:`-1/2` Å\ :sup:`-1`).
+* ``--volume FLOAT``: [Required] Supercell volume (in Å\ :sup:`3`).
+* ``--g INTEGER``: Degeneracy factor. (Default: ``1``).
+* ``-T, --temperature FLOAT``: Temperature in Kelvin for a single-point calculation. (Default: ``300.0``).
+* ``--temperature-range <float float float>``: Temperature range as ``(min, max, count)`` for scanning multiple temperatures.
+* ``--sigma TEXT``: Smearing parameter (float) or interpolation method (``pchip`` or ``cubic``). (Default: ``pchip``).
+* ``--occ-tol FLOAT``: Bose-weight occupation tolerance. (Default: ``1e-5``).
+* ``--overlap-method [Analytic|Integral|HermiteGauss]``: Method for evaluating vibrational overlaps. (Default: ``HermiteGauss``).
+
+
+Scaling Commands
+================
+
+sommerfeld
+----------
+Compute the Sommerfeld enhancement factor.
+
+**Usage:**
+
+.. code-block:: sh
+
+    nonrad sommerfeld [OPTIONS]
+
+**Options:**
+* ``-T, --temperature FLOAT``: Temperature in Kelvin. (Default: ``300.0``).
+* ``--temperature-range <float float float>``: Temperature range as ``(min, max, count)``.
+* ``--z INTEGER``: [Required] Charge ratio *Z = Q / q* (negative values represent attractive potential).
+* ``--m-eff FLOAT``: [Required] Carrier effective mass in units of electron mass.
+* ``--eps0 FLOAT``: [Required] Static dielectric constant.
+* ``--dim [1|2|3]``: Dimensionality of the system. (Default: ``3``).
+* ``--method [Integrate|Analytic]``: Evaluation method. (Default: ``Integrate``).
+
+thermal-velocity
+----------------
+Compute the thermal velocity of a carrier.
+
+**Usage:**
+
+.. code-block:: sh
+
+    nonrad thermal-velocity [OPTIONS]
+
+**Options:**
+* ``-T, --temperature FLOAT``: Temperature in Kelvin. (Default: ``300.0``).
+* ``--temperature-range <float float float>``: Temperature range as ``(min, max, count)``.
+* ``--m-eff FLOAT``: [Required] Carrier effective mass in units of electron mass.
+
+charged-supercell
+-----------------
+Estimate the charged-supercell scaling factor for VASP calculations.
+
+**Usage:**
+
+.. code-block:: sh
+
+    nonrad charged-supercell [OPTIONS] WAVECAR_PATH
+
+**Arguments:**
+* ``WAVECAR_PATH``: Path to the VASP WAVECAR file.
+
+**Options:**
+* ``--bulk-index INTEGER``: [Required] 1-based band index of the bulk wavefunction.
+* ``--def-index INTEGER``: 1-based band index of the defect wavefunction. Provide this or ``--def-coord``. (Default: ``-1``).
+* ``--def-coord <float float float>``: Cartesian coordinates (x, y, z) of the defect position.
+* ``--cutoff FLOAT``: Slope cutoff for plateau detection. (Default: ``0.02``).
+* ``--limit FLOAT``: Upper radial limit (in Å) for the fitting window. (Default: ``5.0``).
+* ``--spin INTEGER``: Spin channel (0 = up, 1 = down). (Default: ``0``).
+* ``--kpoint INTEGER``: 1-based k-point index. (Default: ``1``).
+
+
+Electron-Phonon Commands
+========================
+
+These commands are nested under the ``elphon`` subcommand group and compute the electron-phonon matrix element *Wif*.
+
+elphon wavecars
+---------------
+Compute *Wif* using VASP WAVECAR files along the CCD path.
+
+**Usage:**
+
+.. code-block:: sh
+
+    nonrad elphon wavecars [OPTIONS] CC_DIR GROUND_PATH EXCITED_PATH INIT_WAVECAR
+
+**Arguments:**
+* ``CC_DIR``: Directory containing the CCD displacement folders.
+* ``GROUND_PATH``: Path to the ground-state equilibrium directory.
+* ``EXCITED_PATH``: Path to the excited-state equilibrium directory.
+* ``INIT_WAVECAR``: Path to the initial (reference) WAVECAR file.
+
+**Options:**
+* ``--def-index INTEGER``: [Required] 1-based defect wavefunction index.
+* ``-b, --bulk-index INTEGER``: [Required] 1-based bulk wavefunction index (repeatable for multiple bands).
+* ``--spin INTEGER``: Spin channel (0 = up, 1 = down). (Default: ``0``).
+* ``--kpoint INTEGER``: 1-based k-point index. (Default: ``1``).
+
+elphon wswq
+-----------
+Compute *Wif* using VASP WSWQ files along the CCD path.
+
+**Usage:**
+
+.. code-block:: sh
+
+    nonrad elphon wswq [OPTIONS] CC_DIR GROUND_PATH EXCITED_PATH INIT_VASPRUN
+
+**Arguments:**
+* ``CC_DIR``: Directory containing the CCD displacement folders.
+* ``GROUND_PATH``: Path to the ground-state equilibrium directory.
+* ``EXCITED_PATH``: Path to the excited-state equilibrium directory.
+* ``INIT_VASPRUN``: Path to the reference ``vasprun.xml`` (used to extract eigenvalues).
+
+**Options:**
+* ``--def-index INTEGER``: [Required] 1-based defect wavefunction index.
+* ``-b, --bulk-index INTEGER``: [Required] 1-based bulk wavefunction index (repeatable for multiple bands).
+* ``--spin INTEGER``: Spin channel (0 = up, 1 = down). (Default: ``0``).
+* ``--kpoint INTEGER``: 1-based k-point index. (Default: ``1``).
+
+elphon unk
+----------
+Compute *Wif* using UNK files.
+
+**Usage:**
+
+.. code-block:: sh
+
+    nonrad elphon unk [OPTIONS] INIT_UNK
+
+**Arguments:**
+* ``INIT_UNK``: Path to the initial (reference) UNK file.
+
+**Options:**
+* ``--def-index INTEGER``: [Required] 1-based defect wavefunction index.
+* ``-b, --bulk-index INTEGER``: [Required] 1-based bulk wavefunction index (repeatable).
+* ``--eigs TEXT``: [Required] Comma-separated list of eigenvalues (in eV).
+* ``-u, --unk-pair <float text>``: [Required] A ``(Q, unk_path)`` pair. (Repeatable for all displacements).
diff --git a/docs/index.rst b/docs/index.rst
index d9ec2d2..30a43d5 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -2,6 +2,7 @@
    :hidden:
 
    installation
+   cli
    compatibility
    citation
    Tutorial<tutorial>
diff --git a/nonrad/cli.py b/nonrad/cli.py
new file mode 100644
index 0000000..3d6a850
--- /dev/null
+++ b/nonrad/cli.py
@@ -0,0 +1,859 @@
+# Copyright (c) Chris G. Van de Walle
+# Distributed under the terms of the MIT License.
+
+"""Command-line interface for nonrad.
+
+This module provides a Click-based CLI exposing the main functionality of
+the ``nonrad`` package, including configuration-coordinate diagram (CCD)
+utilities, the nonradiative capture coefficient calculation, scaling
+helpers, and electron-phonon coupling tools.
+
+Commands
+--------
+prep-ccd
+    Prepare displaced structures for a CCD calculation.
+pes
+    Extract potential energy surfaces from CCD VASP outputs.
+dq
+    Compute the configuration coordinate displacement *dQ*.
+q-from-struct
+    Compute the *Q* value for an arbitrary structure.
+barrier
+    Compute the classical barrier height in the harmonic approximation.
+capture
+    Evaluate the nonradiative capture coefficient.
+sommerfeld
+    Compute the Sommerfeld enhancement factor.
+thermal-velocity
+    Compute the thermal velocity of a carrier.
+charged-supercell
+    Estimate the charged-supercell scaling factor.
+elphon wavecars
+    Compute Wif from WAVECAR files.
+elphon wswq
+    Compute Wif from WSWQ files.
+elphon unk
+    Compute Wif from UNK files.
+"""
+
+import os
+from glob import glob
+from pathlib import Path
+from shutil import copyfile
+
+import click
+import numpy as np
+
+
+# ---------------------------------------------------------------------------
+# Main CLI group
+# ---------------------------------------------------------------------------
+
+@click.group()
+@click.version_option(package_name='nonrad')
+def nonrad():
+    """Nonradiative recombination coefficient calculator.
+
+    A command-line tool for computing nonradiative capture coefficients
+    and related quantities from first-principles calculations.
+    """
+
+
+# ---------------------------------------------------------------------------
+# CCD commands
+# ---------------------------------------------------------------------------
+
+@nonrad.command(name='prep-ccd')
+@click.argument('ground_path', type=click.Path(exists=True))
+@click.argument('excited_path', type=click.Path(exists=True))
+@click.argument('cc_dir', type=click.Path())
+@click.option(
+    '--displace', '-d',
+    nargs=3, type=float, default=(-0.5, 0.5, 9),
+    show_default=True,
+    help='Displacement range (min, max, count).',
+)
+@click.option(
+    '--input-files', '-f',
+    type=str, default='KPOINTS,POTCAR,INCAR,job_script.sh',
+    show_default=True,
+    help='Comma-separated list of input files to copy into each sub-directory.',
+)
+def prep_ccd(ground_path, excited_path, cc_dir, displace, input_files):
+    """Prepare displaced structures for a configuration-coordinate diagram.
+
+    Reads CONTCAR from GROUND_PATH and EXCITED_PATH, generates displaced
+    structures, writes POSCAR files into CC_DIR/ground/N and
+    CC_DIR/excited/N sub-directories, and copies the requested input files.
+
+    Parameters
+    ----------
+    ground_path : str
+        Directory containing the ground-state CONTCAR.
+    excited_path : str
+        Directory containing the excited-state CONTCAR.
+    cc_dir : str
+        Output directory for the CCD sub-directories.
+    displace : tuple of float
+        ``(min, max, count)`` passed to ``numpy.linspace``.
+    input_files : str
+        Comma-separated filenames to copy from each source directory.
+    """
+    from pymatgen.core import Structure
+
+    from nonrad.ccd import get_cc_structures
+
+    ground_files = Path(ground_path)
+    excited_files = Path(excited_path)
+    ground_struct = Structure.from_file(str(ground_files / 'CONTCAR'))
+    excited_struct = Structure.from_file(str(excited_files / 'CONTCAR'))
+
+    cc_dir = Path(cc_dir)
+    os.makedirs(str(cc_dir), exist_ok=True)
+    os.makedirs(str(cc_dir / 'ground'), exist_ok=True)
+    os.makedirs(str(cc_dir / 'excited'), exist_ok=True)
+
+    displacements = np.linspace(displace[0], displace[1], int(displace[2]))
+    ground, excited = get_cc_structures(ground_struct, excited_struct,
+                                        displacements)
+
+    files_to_copy = [f.strip() for f in input_files.split(',')]
+
+    for i, struct in enumerate(ground):
+        working_dir = cc_dir / 'ground' / str(i)
+        os.makedirs(str(working_dir), exist_ok=True)
+        struct.to(filename=str(working_dir / 'POSCAR'), fmt='poscar')
+        for f in files_to_copy:
+            src = ground_files / f
+            if src.exists():
+                copyfile(str(src), str(working_dir / f))
+            else:
+                click.echo(f'Warning: {src} not found, skipping.')
+
+    for i, struct in enumerate(excited):
+        working_dir = cc_dir / 'excited' / str(i)
+        os.makedirs(str(working_dir), exist_ok=True)
+        struct.to(filename=str(working_dir / 'POSCAR'), fmt='poscar')
+        for f in files_to_copy:
+            src = excited_files / f
+            if src.exists():
+                copyfile(str(src), str(working_dir / f))
+            else:
+                click.echo(f'Warning: {src} not found, skipping.')
+
+    click.echo(f'Prepared {len(ground)} ground and {len(excited)} excited '
+               f'structures in {cc_dir}')
+
+
+@nonrad.command(name='pes')
+@click.argument('cc_dir', type=click.Path(exists=True))
+@click.argument('ground_files', type=click.Path(exists=True))
+@click.argument('excited_files', type=click.Path(exists=True))
+@click.option(
+    '--energy-diff', '-e',
+    type=float, required=True,
+    help='Energy difference dE between ground and excited states (eV).',
+)
+@click.option(
+    '--plot', '-p',
+    is_flag=True, default=False,
+    help='Plot the potential energy surfaces.',
+)
+@click.option(
+    '--plot-name', '-n',
+    type=str, default='pes.png',
+    show_default=True,
+    help='Filename for the PES plot.',
+)
+def pes(cc_dir, ground_files, excited_files, energy_diff, plot, plot_name):
+    """Extract potential energy surfaces from CCD calculations.
+
+    Reads vasprun.xml files from CC_DIR/ground/*/vasprun.xml and
+    CC_DIR/excited/*/vasprun.xml (including ``.gz`` variants), appends
+    the equilibrium vaspruns from GROUND_FILES and EXCITED_FILES, and
+    computes the PES and harmonic phonon frequencies.
+
+    Parameters
+    ----------
+    cc_dir : str
+        Root directory of the CCD calculation tree.
+    ground_files : str
+        Directory with the ground-state equilibrium calculation.
+    excited_files : str
+        Directory with the excited-state equilibrium calculation.
+    energy_diff : float
+        Energy offset *dE* between the two PES (eV).
+    plot : bool
+        Whether to save a PES plot.
+    plot_name : str
+        Output filename for the plot image.
+    """
+    from pymatgen.core import Structure
+
+    from nonrad.ccd import (
+        get_dQ,
+        get_omega_from_PES,
+        get_PES_from_vaspruns,
+    )
+
+    ground_struct = Structure.from_file(str(Path(ground_files) / 'CONTCAR'))
+    excited_struct = Structure.from_file(str(Path(excited_files) / 'CONTCAR'))
+
+    dQ = get_dQ(ground_struct, excited_struct)
+    click.echo(f'dQ = {dQ:.2f} amu^{{1/2}} Angstrom')
+
+    # Collect vasprun.xml (and .gz) from the CCD sub-directories
+    ground_vaspruns = (
+        glob(str(Path(cc_dir) / 'ground' / '*' / 'vasprun.xml')) +
+        glob(str(Path(cc_dir) / 'ground' / '*' / 'vasprun.xml.gz'))
+    )
+    excited_vaspruns = (
+        glob(str(Path(cc_dir) / 'excited' / '*' / 'vasprun.xml')) +
+        glob(str(Path(cc_dir) / 'excited' / '*' / 'vasprun.xml.gz'))
+    )
+
+    # Append equilibrium vaspruns
+    for gf in [str(Path(ground_files) / 'vasprun.xml'),
+               str(Path(ground_files) / 'vasprun.xml.gz')]:
+        if os.path.isfile(gf):
+            ground_vaspruns.append(gf)
+            break
+
+    for ef in [str(Path(excited_files) / 'vasprun.xml'),
+               str(Path(excited_files) / 'vasprun.xml.gz')]:
+        if os.path.isfile(ef):
+            excited_vaspruns.append(ef)
+            break
+
+    Q_ground, E_ground = get_PES_from_vaspruns(
+        ground_struct, excited_struct, ground_vaspruns,
+    )
+    Q_excited, E_excited = get_PES_from_vaspruns(
+        ground_struct, excited_struct, excited_vaspruns,
+    )
+    E_excited = energy_diff + E_excited
+
+    # Always compute omega
+    if plot:
+        import matplotlib
+        matplotlib.use('Agg')
+        import matplotlib.pyplot as plt
+
+        _, ax = plt.subplots(figsize=(5, 5))
+        ax.scatter(Q_ground, E_ground, s=10)
+        ax.scatter(Q_excited, E_excited, s=10)
+        q = np.linspace(-1.0, 3.5, 100)
+        ground_omega = get_omega_from_PES(Q_ground, E_ground, ax=ax, q=q)
+        excited_omega = get_omega_from_PES(Q_excited, E_excited, ax=ax, q=q)
+        ax.set_xlabel(r'$Q$ [amu$^{1/2}$ $\AA$]')
+        ax.set_ylabel(r'$E$ [eV]')
+        plt.savefig(plot_name, dpi=300, format=plot_name.split('.')[-1])
+        click.echo(f'Plot saved to {plot_name}')
+    else:
+        ground_omega = get_omega_from_PES(Q_ground, E_ground)
+        excited_omega = get_omega_from_PES(Q_excited, E_excited)
+
+    click.echo(f'Ground state omega = {ground_omega:.2f} eV')
+    click.echo(f'Excited state omega = {excited_omega:.2f} eV')
+
+
+@nonrad.command(name='dq')
+@click.argument('ground_path', type=click.Path(exists=True))
+@click.argument('excited_path', type=click.Path(exists=True))
+def dq_cmd(ground_path, excited_path):
+    """Compute the configuration-coordinate displacement dQ.
+
+    GROUND_PATH and EXCITED_PATH are structure files (e.g. POSCAR,
+    CONTCAR, CIF) — *not* directories.
+
+    Parameters
+    ----------
+    ground_path : str
+        Path to the ground-state structure file.
+    excited_path : str
+        Path to the excited-state structure file.
+    """
+    from pymatgen.core import Structure
+
+    from nonrad.ccd import get_dQ
+
+    ground_struct = Structure.from_file(ground_path)
+    excited_struct = Structure.from_file(excited_path)
+    dQ = get_dQ(ground_struct, excited_struct)
+    click.echo(f'dQ = {dQ:.6f} amu^{{1/2}} Angstrom')
+
+
+@nonrad.command(name='q-from-struct')
+@click.argument('ground_path', type=click.Path(exists=True))
+@click.argument('excited_path', type=click.Path(exists=True))
+@click.argument('struct_path', type=click.Path(exists=True))
+@click.option(
+    '--tol', '-t',
+    type=float, default=1e-4,
+    show_default=True,
+    help='Distance tolerance for filtering coordinates.',
+)
+@click.option(
+    '--nround',
+    type=int, default=5,
+    show_default=True,
+    help='Number of decimal places for rounding Q.',
+)
+def q_from_struct(ground_path, excited_path, struct_path, tol, nround):
+    """Compute the Q value of an arbitrary structure.
+
+    Given the ground and excited endpoint structures, determine the Q
+    coordinate of STRUCT_PATH assuming linear interpolation.
+
+    Parameters
+    ----------
+    ground_path : str
+        Path to the ground-state structure file.
+    excited_path : str
+        Path to the excited-state structure file.
+    struct_path : str
+        Path to the structure file to evaluate.
+    tol : float
+        Distance cutoff for discarding nearly-stationary sites.
+    nround : int
+        Decimal places used when rounding Q.
+    """
+    from pymatgen.core import Structure
+
+    from nonrad.ccd import get_Q_from_struct
+
+    ground_struct = Structure.from_file(ground_path)
+    excited_struct = Structure.from_file(excited_path)
+    Q = get_Q_from_struct(ground_struct, excited_struct, struct_path,
+                          tol=tol, nround=nround)
+    click.echo(f'Q = {Q}')
+
+
+@nonrad.command(name='barrier')
+@click.option('--dq', type=float, required=True,
+              help='Displacement dQ (amu^{1/2} Angstrom).')
+@click.option('--de', type=float, required=True,
+              help='Energy offset dE (eV).')
+@click.option('--wi', type=float, required=True,
+              help='Initial-state frequency (eV).')
+@click.option('--wf', type=float, required=True,
+              help='Final-state frequency (eV).')
+def barrier_cmd(dq, de, wi, wf):
+    """Compute the classical barrier height in the harmonic approximation.
+
+    Parameters
+    ----------
+    dq : float
+        Configuration-coordinate displacement (amu^{1/2} Angstrom).
+    de : float
+        Energy offset between oscillators (eV).
+    wi : float
+        Initial harmonic-oscillator frequency (eV).
+    wf : float
+        Final harmonic-oscillator frequency (eV).
+    """
+    from nonrad.ccd import get_barrier_harmonic
+
+    result = get_barrier_harmonic(dq, de, wi, wf)
+    if result is not None:
+        click.echo(f'Barrier height = {result:.6f} eV')
+    else:
+        click.echo('No crossing point found.')
+
+
+# ---------------------------------------------------------------------------
+# Core capture coefficient command
+# ---------------------------------------------------------------------------
+
+@nonrad.command(name='capture')
+@click.option('--dq', type=float, required=True,
+              help='Displacement dQ (amu^{1/2} Angstrom).')
+@click.option('--de', type=float, required=True,
+              help='Energy offset dE (eV).')
+@click.option('--wi', type=float, required=True,
+              help='Initial-state frequency (eV).')
+@click.option('--wf', type=float, required=True,
+              help='Final-state frequency (eV).')
+@click.option('--wif', type=float, required=True,
+              help='Electron-phonon coupling Wif (eV amu^{-1/2} Angstrom^{-1}).')
+@click.option('--volume', type=float, required=True,
+              help='Supercell volume (Angstrom^3).')
+@click.option('--g', type=int, default=1, show_default=True,
+              help='Degeneracy factor.')
+@click.option('--temperature', '-T', type=float, default=300.,
+              show_default=True,
+              help='Temperature (K) for a single-point evaluation.')
+@click.option('--temperature-range', nargs=3, type=float, default=None,
+              help='Temperature range as (min, max, count) for scanning.')
+@click.option('--sigma', type=str, default='pchip', show_default=True,
+              help='Smearing parameter (float) or interpolation method '
+                   '("pchip" or "cubic").')
+@click.option('--occ-tol', type=float, default=1e-5, show_default=True,
+              help='Occupation tolerance for Bose weights.')
+@click.option('--overlap-method',
+              type=click.Choice(['Analytic', 'Integral', 'HermiteGauss'],
+                                case_sensitive=False),
+              default='HermiteGauss', show_default=True,
+              help='Method for evaluating vibrational overlaps.')
+def capture_cmd(dq, de, wi, wf, wif, volume, g, temperature,
+                temperature_range, sigma, occ_tol, overlap_method):
+    """Compute the nonradiative capture coefficient.
+
+    Evaluates the capture coefficient following Alkauskas *et al.*,
+    Phys. Rev. B 90, 075202 (2014).  The result is **unscaled** and has
+    units of cm^3 s^{-1}.
+
+    Parameters
+    ----------
+    dq : float
+        Configuration-coordinate displacement (amu^{1/2} Angstrom).
+    de : float
+        Energy offset between oscillators (eV).
+    wi : float
+        Initial-state frequency (eV).
+    wf : float
+        Final-state frequency (eV).
+    wif : float
+        Electron-phonon coupling matrix element
+        (eV amu^{-1/2} Angstrom^{-1}).
+    volume : float
+        Supercell volume (Angstrom^3).
+    g : int
+        Degeneracy factor of the final state.
+    temperature : float
+        Single temperature (K).
+    temperature_range : tuple of float or None
+        ``(min, max, count)`` for a temperature scan.
+    sigma : str
+        Smearing parameter or interpolation keyword.
+    occ_tol : float
+        Bose-weight occupation tolerance.
+    overlap_method : str
+        Overlap evaluation method.
+    """
+    from nonrad.nonrad import get_C
+
+    # Parse sigma: interpret as float when possible
+    try:
+        sigma_val = float(sigma)
+    except ValueError:
+        sigma_val = sigma
+
+    if temperature_range is not None:
+        temps = np.linspace(temperature_range[0], temperature_range[1],
+                            int(temperature_range[2]))
+    else:
+        temps = temperature
+
+    result = get_C(dq, de, wi, wf, wif, volume,
+                   g=g, T=temps, sigma=sigma_val,
+                   occ_tol=occ_tol, overlap_method=overlap_method)
+
+    if isinstance(temps, np.ndarray):
+        click.echo(f'{"T (K)":>12s}  {"C (cm^3/s)":>14s}')
+        click.echo('-' * 28)
+        for t, c in zip(temps, np.atleast_1d(result)):
+            click.echo(f'{t:12.2f}  {c:14.6e}')
+    else:
+        click.echo(f'C = {result:.6e} cm^3/s  (T = {temps} K)')
+
+
+# ---------------------------------------------------------------------------
+# Scaling commands
+# ---------------------------------------------------------------------------
+
+@nonrad.command(name='sommerfeld')
+@click.option('--temperature', '-T', type=float, default=300.,
+              show_default=True, help='Temperature (K).')
+@click.option('--temperature-range', nargs=3, type=float, default=None,
+              help='Temperature range as (min, max, count).')
+@click.option('--z', type=int, required=True,
+              help='Charge ratio Z = Q/q (negative = attractive).')
+@click.option('--m-eff', type=float, required=True,
+              help='Effective mass in units of electron mass.')
+@click.option('--eps0', type=float, required=True,
+              help='Static dielectric constant.')
+@click.option('--dim',
+              type=click.Choice(['1', '2', '3'], case_sensitive=False),
+              default='3', show_default=True,
+              help='Dimensionality of the system.')
+@click.option('--method',
+              type=click.Choice(['Integrate', 'Analytic'],
+                                case_sensitive=False),
+              default='Integrate', show_default=True,
+              help='Evaluation method.')
+def sommerfeld_cmd(temperature, temperature_range, z, m_eff, eps0, dim,
+                   method):
+    """Compute the Sommerfeld enhancement factor.
+
+    Evaluates the temperature-dependent Sommerfeld parameter following
+    Pässler *et al.*, phys. stat. sol. (b) 78, 625 (1976).
+
+    Parameters
+    ----------
+    temperature : float
+        Temperature (K) for a single-point evaluation.
+    temperature_range : tuple of float or None
+        ``(min, max, count)`` for a temperature scan.
+    z : int
+        Charge ratio *Z = Q / q*.
+    m_eff : float
+        Carrier effective mass in units of *m_e*.
+    eps0 : float
+        Static dielectric constant.
+    dim : str
+        System dimensionality (``'1'``, ``'2'``, or ``'3'``).
+    method : str
+        ``'Integrate'`` or ``'Analytic'``.
+    """
+    from nonrad.scaling import sommerfeld_parameter
+
+    dim_int = int(dim)
+
+    if temperature_range is not None:
+        temps = np.linspace(temperature_range[0], temperature_range[1],
+                            int(temperature_range[2]))
+    else:
+        temps = temperature
+
+    result = sommerfeld_parameter(temps, z, m_eff, eps0,
+                                  dim=dim_int, method=method)
+
+    if isinstance(temps, np.ndarray):
+        click.echo(f'{"T (K)":>12s}  {"Sommerfeld":>14s}')
+        click.echo('-' * 28)
+        for t, s in zip(temps, np.atleast_1d(result)):
+            click.echo(f'{t:12.2f}  {s:14.6e}')
+    else:
+        click.echo(f'Sommerfeld factor = {result:.6e}  (T = {temps} K)')
+
+
+@nonrad.command(name='thermal-velocity')
+@click.option('--temperature', '-T', type=float, default=300.,
+              show_default=True, help='Temperature (K).')
+@click.option('--temperature-range', nargs=3, type=float, default=None,
+              help='Temperature range as (min, max, count).')
+@click.option('--m-eff', type=float, required=True,
+              help='Effective mass in units of electron mass.')
+def thermal_velocity_cmd(temperature, temperature_range, m_eff):
+    """Compute the thermal velocity of a carrier.
+
+    Parameters
+    ----------
+    temperature : float
+        Temperature (K) for a single-point evaluation.
+    temperature_range : tuple of float or None
+        ``(min, max, count)`` for a temperature scan.
+    m_eff : float
+        Carrier effective mass in units of *m_e*.
+    """
+    from nonrad.scaling import thermal_velocity
+
+    if temperature_range is not None:
+        temps = np.linspace(temperature_range[0], temperature_range[1],
+                            int(temperature_range[2]))
+    else:
+        temps = temperature
+
+    result = thermal_velocity(temps, m_eff)
+
+    if isinstance(temps, np.ndarray):
+        click.echo(f'{"T (K)":>12s}  {"v_th (cm/s)":>14s}')
+        click.echo('-' * 28)
+        for t, v in zip(temps, np.atleast_1d(result)):
+            click.echo(f'{t:12.2f}  {v:14.6e}')
+    else:
+        click.echo(f'Thermal velocity = {result:.6e} cm/s  (T = {temps} K)')
+
+
+@nonrad.command(name='charged-supercell')
+@click.argument('wavecar_path', type=click.Path(exists=True))
+@click.option('--bulk-index', type=int, required=True,
+              help='Index of the bulk wavefunction (1-based).')
+@click.option('--def-index', type=int, default=-1, show_default=True,
+              help='Index of the defect wavefunction (1-based). '
+                   'Provide this or --def-coord.')
+@click.option('--def-coord', nargs=3, type=float, default=None,
+              help='Cartesian coordinates (x y z) of the defect position.')
+@click.option('--cutoff', type=float, default=0.02, show_default=True,
+              help='Slope cutoff for plateau detection.')
+@click.option('--limit', type=float, default=5., show_default=True,
+              help='Upper radial limit (Angstrom) for the fitting window.')
+@click.option('--spin', type=int, default=0, show_default=True,
+              help='Spin channel (0 = up, 1 = down).')
+@click.option('--kpoint', type=int, default=1, show_default=True,
+              help='k-point index (1-based).')
+def charged_supercell_cmd(wavecar_path, bulk_index, def_index, def_coord,
+                          cutoff, limit, spin, kpoint):
+    """Estimate the charged-supercell scaling factor.
+
+    Compares the radial distribution of the bulk wavefunction around the
+    defect position to a uniform distribution and returns a scaling
+    factor for the capture coefficient.
+
+    Parameters
+    ----------
+    wavecar_path : str
+        Path to the WAVECAR file.
+    bulk_index : int
+        Band index of the bulk wavefunction (1-based).
+    def_index : int
+        Band index of the defect wavefunction (1-based), or -1 if
+        ``def_coord`` is provided instead.
+    def_coord : tuple of float or None
+        Cartesian coordinates of the defect centre.
+    cutoff : float
+        Gradient cutoff for plateau identification.
+    limit : float
+        Upper radius limit for curve fitting (Angstrom).
+    spin : int
+        Spin channel (0 = up, 1 = down).
+    kpoint : int
+        k-point index (1-based).
+    """
+    from nonrad.scaling import charged_supercell_scaling_VASP
+
+    if def_index == -1 and def_coord is None:
+        raise click.UsageError(
+            'Either --def-index (not -1) or --def-coord must be specified.'
+        )
+
+    coord = np.array(def_coord) if def_coord is not None else None
+
+    scaling = charged_supercell_scaling_VASP(
+        wavecar_path, bulk_index,
+        def_index=def_index, def_coord=coord,
+        cutoff=cutoff, limit=limit,
+        spin=spin, kpoint=kpoint,
+    )
+    click.echo(f'Charged-supercell scaling factor = {scaling:.6f}')
+
+
+# ---------------------------------------------------------------------------
+# Electron-phonon sub-group
+# ---------------------------------------------------------------------------
+
+@nonrad.group('elphon')
+def elphon():
+    """Electron-phonon coupling utilities.
+
+    Sub-commands for computing the electron-phonon coupling matrix
+    element Wif using WAVECARs, WSWQ files, or UNK files.
+    """
+
+
+@elphon.command(name='wavecars')
+@click.argument('cc_dir', type=click.Path(exists=True))
+@click.argument('ground_path', type=click.Path(exists=True))
+@click.argument('excited_path', type=click.Path(exists=True))
+@click.argument('init_wavecar', type=click.Path(exists=True))
+@click.option('--def-index', type=int, required=True,
+              help='Defect wavefunction index (1-based).')
+@click.option('--bulk-index', '-b', type=int, required=True, multiple=True,
+              help='Bulk wavefunction index (1-based, repeatable).')
+@click.option('--spin', type=int, default=0, show_default=True,
+              help='Spin channel (0 = up, 1 = down).')
+@click.option('--kpoint', type=int, default=1, show_default=True,
+              help='k-point index (1-based).')
+def elphon_wavecars(cc_dir, ground_path, excited_path, init_wavecar,
+                    def_index, bulk_index, spin, kpoint):
+    """Compute Wif from WAVECAR files along the CCD path.
+
+    Automatically discovers sub-directories inside CC_DIR that contain
+    both a ``vasprun.xml`` (or ``.gz``) and a ``WAVECAR`` file.
+
+    Parameters
+    ----------
+    cc_dir : str
+        Directory tree with CCD calculations.
+    ground_path : str
+        Directory containing the ground-state CONTCAR.
+    excited_path : str
+        Directory containing the excited-state CONTCAR.
+    init_wavecar : str
+        Path to the initial (reference) WAVECAR.
+    def_index : int
+        Band index of the defect wavefunction.
+    bulk_index : tuple of int
+        Band indices of the bulk wavefunctions.
+    spin : int
+        Spin channel.
+    kpoint : int
+        k-point index.
+    """
+    from pymatgen.core import Structure
+    from pymatgen.io.vasp.outputs import Vasprun
+
+    from nonrad.ccd import get_Q_from_struct
+    from nonrad.elphon import get_Wif_from_wavecars
+
+    ground_struct = Structure.from_file(str(Path(ground_path) / 'CONTCAR'))
+    excited_struct = Structure.from_file(str(Path(excited_path) / 'CONTCAR'))
+
+    # Auto-discover sub-directories with vasprun + WAVECAR
+    wavecars = []
+    for subdir in sorted(Path(cc_dir).iterdir()):
+        if not subdir.is_dir():
+            continue
+        wavecar_file = subdir / 'WAVECAR'
+        if not wavecar_file.exists():
+            continue
+        vr_path = None
+        for vr_name in ('vasprun.xml', 'vasprun.xml.gz'):
+            candidate = subdir / vr_name
+            if candidate.exists():
+                vr_path = candidate
+                break
+        if vr_path is None:
+            continue
+        vr = Vasprun(str(vr_path), parse_dos=False, parse_eigen=False)
+        q = get_Q_from_struct(ground_struct, excited_struct,
+                              vr.structures[-1])
+        wavecars.append((q, str(wavecar_file)))
+
+    results = get_Wif_from_wavecars(
+        wavecars, init_wavecar, def_index, list(bulk_index),
+        spin=spin, kpoint=kpoint,
+    )
+
+    click.echo(f'{"bulk_index":>12s}  {"Wif (eV amu^{-1/2} A^{-1})":>30s}')
+    click.echo('-' * 44)
+    for bi, wif_val in results:
+        click.echo(f'{bi:12d}  {wif_val:30.6e}')
+
+
+@elphon.command(name='wswq')
+@click.argument('cc_dir', type=click.Path(exists=True))
+@click.argument('ground_path', type=click.Path(exists=True))
+@click.argument('excited_path', type=click.Path(exists=True))
+@click.argument('init_vasprun', type=click.Path(exists=True))
+@click.option('--def-index', type=int, required=True,
+              help='Defect wavefunction index (1-based).')
+@click.option('--bulk-index', '-b', type=int, required=True, multiple=True,
+              help='Bulk wavefunction index (1-based, repeatable).')
+@click.option('--spin', type=int, default=0, show_default=True,
+              help='Spin channel (0 = up, 1 = down).')
+@click.option('--kpoint', type=int, default=1, show_default=True,
+              help='k-point index (1-based).')
+def elphon_wswq(cc_dir, ground_path, excited_path, init_vasprun,
+                def_index, bulk_index, spin, kpoint):
+    """Compute Wif from WSWQ files along the CCD path.
+
+    Automatically discovers sub-directories inside CC_DIR that contain
+    both a ``vasprun.xml`` (or ``.gz``) and a ``WSWQ`` (or ``.gz``)
+    file.
+
+    Parameters
+    ----------
+    cc_dir : str
+        Directory tree with CCD calculations.
+    ground_path : str
+        Directory containing the ground-state CONTCAR.
+    excited_path : str
+        Directory containing the excited-state CONTCAR.
+    init_vasprun : str
+        Path to the initial vasprun.xml for eigenvalue extraction.
+    def_index : int
+        Band index of the defect wavefunction.
+    bulk_index : tuple of int
+        Band indices of the bulk wavefunctions.
+    spin : int
+        Spin channel.
+    kpoint : int
+        k-point index.
+    """
+    from pymatgen.core import Structure
+    from pymatgen.io.vasp.outputs import Vasprun
+
+    from nonrad.ccd import get_Q_from_struct
+    from nonrad.elphon import get_Wif_from_WSWQ
+
+    ground_struct = Structure.from_file(str(Path(ground_path) / 'CONTCAR'))
+    excited_struct = Structure.from_file(str(Path(excited_path) / 'CONTCAR'))
+
+    # Auto-discover sub-directories with vasprun + WSWQ
+    wswqs = []
+    for subdir in sorted(Path(cc_dir).iterdir()):
+        if not subdir.is_dir():
+            continue
+        wswq_file = None
+        for wswq_name in ('WSWQ', 'WSWQ.gz'):
+            candidate = subdir / wswq_name
+            if candidate.exists():
+                wswq_file = candidate
+                break
+        if wswq_file is None:
+            continue
+        vr_path = None
+        for vr_name in ('vasprun.xml', 'vasprun.xml.gz'):
+            candidate = subdir / vr_name
+            if candidate.exists():
+                vr_path = candidate
+                break
+        if vr_path is None:
+            continue
+        vr = Vasprun(str(vr_path), parse_dos=False, parse_eigen=False)
+        q = get_Q_from_struct(ground_struct, excited_struct,
+                              vr.structures[-1])
+        wswqs.append((q, str(wswq_file)))
+
+    results = get_Wif_from_WSWQ(
+        wswqs, init_vasprun, def_index, list(bulk_index),
+        spin=spin, kpoint=kpoint,
+    )
+
+    click.echo(f'{"bulk_index":>12s}  {"Wif (eV amu^{-1/2} A^{-1})":>30s}')
+    click.echo('-' * 44)
+    for bi, wif_val in results:
+        click.echo(f'{bi:12d}  {wif_val:30.6e}')
+
+
+@elphon.command(name='unk')
+@click.argument('init_unk', type=click.Path(exists=True))
+@click.option('--def-index', type=int, required=True,
+              help='Defect wavefunction index (1-based).')
+@click.option('--bulk-index', '-b', type=int, required=True, multiple=True,
+              help='Bulk wavefunction index (1-based, repeatable).')
+@click.option('--eigs', type=str, required=True,
+              help='Comma-separated eigenvalues (eV).')
+@click.option('--unk-pair', '-u', type=(float, str), required=True,
+              multiple=True,
+              help='(Q, unk_path) pair — repeatable.')
+def elphon_unk(init_unk, def_index, bulk_index, eigs, unk_pair):
+    """Compute Wif from UNK files.
+
+    Parameters
+    ----------
+    init_unk : str
+        Path to the initial (reference) UNK file.
+    def_index : int
+        Band index of the defect wavefunction.
+    bulk_index : tuple of int
+        Band indices of the bulk wavefunctions.
+    eigs : str
+        Comma-separated string of eigenvalues (eV).
+    unk_pair : tuple of (float, str)
+        One or more ``(Q, path)`` pairs specifying displaced UNK files.
+    """
+    from nonrad.elphon import get_Wif_from_UNK
+
+    eigs_array = np.array([float(e) for e in eigs.split(',')])
+    unk_list = [(q, path) for q, path in unk_pair]
+
+    results = get_Wif_from_UNK(
+        unk_list, init_unk, def_index, list(bulk_index),
+        eigs=eigs_array,
+    )
+
+    click.echo(f'{"bulk_index":>12s}  {"Wif (eV amu^{-1/2} A^{-1})":>30s}')
+    click.echo('-' * 44)
+    for bi, wif_val in results:
+        click.echo(f'{bi:12d}  {wif_val:30.6e}')
+
+
+# ---------------------------------------------------------------------------
+# Entry point
+# ---------------------------------------------------------------------------
+
+if __name__ == '__main__':
+    nonrad()
diff --git a/nonrad/tests/test_cli.py b/nonrad/tests/test_cli.py
new file mode 100644
index 0000000..bfd09da
--- /dev/null
+++ b/nonrad/tests/test_cli.py
@@ -0,0 +1,606 @@
+# pylint: disable=C0114,C0115,C0116
+
+import os
+import unittest
+from pathlib import Path
+from unittest.mock import MagicMock, patch
+
+import numpy as np
+from click.testing import CliRunner
+
+from nonrad.cli import nonrad as nonrad_cli
+from nonrad.tests import TEST_FILES
+
+
+# ---------------------------------------------------------------------------
+# Tests for the top-level `nonrad` group
+# ---------------------------------------------------------------------------
+class TestNonradGroup(unittest.TestCase):
+    def setUp(self):
+        self.runner = CliRunner()
+
+    def test_help(self):
+        result = self.runner.invoke(nonrad_cli, ['--help'])
+        self.assertEqual(result.exit_code, 0, result.output)
+        self.assertIn('Nonradiative recombination', result.output)
+
+    def test_no_args(self):
+        result = self.runner.invoke(nonrad_cli, [])
+        self.assertEqual(result.exit_code, 0, result.output)
+
+
+# ---------------------------------------------------------------------------
+# Tests for `prep-ccd`
+# ---------------------------------------------------------------------------
+class TestPrepCcd(unittest.TestCase):
+    def setUp(self):
+        self.runner = CliRunner()
+
+    def test_help(self):
+        result = self.runner.invoke(nonrad_cli, ['prep-ccd', '--help'])
+        self.assertEqual(result.exit_code, 0, result.output)
+        self.assertIn('ground_path', result.output.lower())
+
+    @patch('nonrad.ccd.get_cc_structures')
+    @patch('pymatgen.core.structure.Structure.from_file')
+    def test_basic(self, mock_from_file, mock_get_cc):
+        mock_struct = MagicMock()
+        mock_from_file.return_value = mock_struct
+
+        # get_cc_structures returns (list_of_structs, list_of_structs)
+        struct_a = MagicMock()
+        struct_b = MagicMock()
+        mock_get_cc.return_value = ([struct_a], [struct_b])
+
+        with self.runner.isolated_filesystem():
+            # Create fake ground/excited directories with required files
+            os.makedirs('ground')
+            Path('ground/CONTCAR').touch()
+            for f in ['KPOINTS', 'POTCAR', 'INCAR', 'job_script.sh']:
+                Path(f'ground/{f}').write_text('dummy')
+
+            os.makedirs('excited')
+            Path('excited/CONTCAR').touch()
+            for f in ['KPOINTS', 'POTCAR', 'INCAR', 'job_script.sh']:
+                Path(f'excited/{f}').write_text('dummy')
+
+            result = self.runner.invoke(
+                nonrad_cli,
+                ['prep-ccd', 'ground', 'excited', 'ccd_output']
+            )
+            self.assertEqual(result.exit_code, 0,
+                             f'exit_code={result.exit_code}\n{result.output}')
+            # Verify output directories were created
+            self.assertTrue(os.path.isdir('ccd_output'))
+            self.assertTrue(os.path.isdir('ccd_output/ground'))
+            self.assertTrue(os.path.isdir('ccd_output/excited'))
+            # Verify struct.to() was called for each displacement struct
+            struct_a.to.assert_called()
+            struct_b.to.assert_called()
+
+    @patch('nonrad.ccd.get_cc_structures')
+    @patch('pymatgen.core.structure.Structure.from_file')
+    def test_custom_displacements(self, mock_from_file, mock_get_cc):
+        mock_struct = MagicMock()
+        mock_from_file.return_value = mock_struct
+        mock_get_cc.return_value = ([], [])
+
+        with self.runner.isolated_filesystem():
+            os.makedirs('ground')
+            Path('ground/CONTCAR').touch()
+            for f in ['KPOINTS', 'POTCAR', 'INCAR', 'job_script.sh']:
+                Path(f'ground/{f}').write_text('dummy')
+
+            os.makedirs('excited')
+            Path('excited/CONTCAR').touch()
+            for f in ['KPOINTS', 'POTCAR', 'INCAR', 'job_script.sh']:
+                Path(f'excited/{f}').write_text('dummy')
+
+            result = self.runner.invoke(
+                nonrad_cli,
+                ['prep-ccd', 'ground', 'excited', 'ccd_output',
+                 '-d', '-0.3', '0.3', '5']
+            )
+            self.assertEqual(result.exit_code, 0,
+                             f'exit_code={result.exit_code}\n{result.output}')
+            # Verify get_cc_structures was called
+            mock_get_cc.assert_called_once()
+            # Check displacements were passed correctly (via linspace)
+            call_args = mock_get_cc.call_args
+            displacements = call_args[0][2]
+            self.assertEqual(len(displacements), 5)
+
+    @patch('nonrad.ccd.get_cc_structures')
+    @patch('pymatgen.core.structure.Structure.from_file')
+    def test_nonexistent_ground_path(self, mock_from_file, mock_get_cc):
+        with self.runner.isolated_filesystem():
+            result = self.runner.invoke(
+                nonrad_cli,
+                ['prep-ccd', 'nonexistent', 'also_nonexistent', 'out']
+            )
+            self.assertNotEqual(result.exit_code, 0)
+
+
+# ---------------------------------------------------------------------------
+# Tests for `pes`
+# ---------------------------------------------------------------------------
+class TestPes(unittest.TestCase):
+    def setUp(self):
+        self.runner = CliRunner()
+
+    def test_help(self):
+        result = self.runner.invoke(nonrad_cli, ['pes', '--help'])
+        self.assertEqual(result.exit_code, 0, result.output)
+        self.assertIn('cc_dir', result.output.lower())
+
+    @patch('nonrad.ccd.get_omega_from_PES', return_value=0.05)
+    @patch('nonrad.ccd.get_PES_from_vaspruns',
+           return_value=(np.array([0.0, 1.0]), np.array([0.0, 0.5])))
+    @patch('nonrad.ccd.get_dQ', return_value=1.5)
+    @patch('nonrad.cli.glob', return_value=[])
+    @patch('pymatgen.core.structure.Structure.from_file')
+    def test_basic(self, mock_from_file, mock_glob,
+                   mock_get_dq, mock_get_pes, mock_get_omega):
+        mock_struct = MagicMock()
+        mock_from_file.return_value = mock_struct
+
+        with self.runner.isolated_filesystem():
+            # Create fake dirs / files
+            os.makedirs('cc/ground/0', exist_ok=True)
+            os.makedirs('cc/excited/0', exist_ok=True)
+            os.makedirs('ground_files')
+            Path('ground_files/CONTCAR').touch()
+            Path('ground_files/vasprun.xml').touch()
+            os.makedirs('excited_files')
+            Path('excited_files/CONTCAR').touch()
+            Path('excited_files/vasprun.xml').touch()
+
+            result = self.runner.invoke(
+                nonrad_cli,
+                ['pes', 'cc', 'ground_files', 'excited_files', '--energy-diff', '1.0']
+            )
+            # Verify exit code is 0 and mock gets called
+            self.assertEqual(result.exit_code, 0, result.output)
+            mock_get_dq.assert_called_once()
+
+
+# ---------------------------------------------------------------------------
+# Tests for `dq`
+# ---------------------------------------------------------------------------
+class TestDq(unittest.TestCase):
+    def setUp(self):
+        self.runner = CliRunner()
+
+    def test_help(self):
+        result = self.runner.invoke(nonrad_cli, ['dq', '--help'])
+        self.assertEqual(result.exit_code, 0, result.output)
+
+    def test_basic_real_files(self):
+        ground = str(TEST_FILES / 'POSCAR.C0.gz')
+        excited = str(TEST_FILES / 'POSCAR.C-.gz')
+        result = self.runner.invoke(
+            nonrad_cli, ['dq', ground, excited]
+        )
+        self.assertEqual(result.exit_code, 0,
+                         f'exit_code={result.exit_code}\n{result.output}')
+        # dQ ≈ 1.6859
+        self.assertIn('1.6858', result.output)
+
+    def test_nonexistent_file(self):
+        result = self.runner.invoke(
+            nonrad_cli, ['dq', 'nonexistent_file', 'also_nonexistent']
+        )
+        self.assertNotEqual(result.exit_code, 0)
+
+
+# ---------------------------------------------------------------------------
+# Tests for `q-from-struct`
+# ---------------------------------------------------------------------------
+class TestQFromStruct(unittest.TestCase):
+    def setUp(self):
+        self.runner = CliRunner()
+
+    def test_help(self):
+        result = self.runner.invoke(nonrad_cli, ['q-from-struct', '--help'])
+        self.assertEqual(result.exit_code, 0, result.output)
+
+    def test_basic_real_files(self):
+        ground = str(TEST_FILES / 'POSCAR.C0.gz')
+        excited = str(TEST_FILES / 'POSCAR.C-.gz')
+        struct = str(TEST_FILES / 'POSCAR.C0.gz')
+        result = self.runner.invoke(
+            nonrad_cli, ['q-from-struct', ground, excited, struct]
+        )
+        self.assertEqual(result.exit_code, 0,
+                         f'exit_code={result.exit_code}\n{result.output}')
+        # Q for ground == ground should be ≈ 0
+        self.assertIn('0.0', result.output)
+
+    def test_nonexistent_file(self):
+        result = self.runner.invoke(
+            nonrad_cli,
+            ['q-from-struct', 'nonexistent', 'nonexistent', 'nonexistent']
+        )
+        self.assertNotEqual(result.exit_code, 0)
+
+
+# ---------------------------------------------------------------------------
+# Tests for `barrier`
+# ---------------------------------------------------------------------------
+class TestBarrier(unittest.TestCase):
+    def setUp(self):
+        self.runner = CliRunner()
+
+    def test_help(self):
+        result = self.runner.invoke(nonrad_cli, ['barrier', '--help'])
+        self.assertEqual(result.exit_code, 0, result.output)
+
+    def test_crossing(self):
+        result = self.runner.invoke(
+            nonrad_cli,
+            ['barrier',
+             '--dq', '1.0', '--de', '0.0', '--wi', '0.05', '--wf', '0.05']
+        )
+        self.assertEqual(result.exit_code, 0,
+                         f'exit_code={result.exit_code}\n{result.output}')
+        # Should contain barrier height info
+        output_lower = result.output.lower()
+        self.assertTrue(
+            'barrier' in output_lower or 'height' in output_lower
+            or 'ev' in output_lower,
+            f'Expected barrier output, got: {result.output}'
+        )
+
+    @patch('nonrad.ccd.get_barrier_harmonic', return_value=None)
+    def test_no_crossing(self, mock_barrier):
+        result = self.runner.invoke(
+            nonrad_cli,
+            ['barrier',
+             '--dq', '1.0', '--de', '5.0', '--wi', '0.01', '--wf', '0.01']
+        )
+        # The CLI should handle None return (no crossing) gracefully
+        self.assertEqual(result.exit_code, 0,
+                         f'exit_code={result.exit_code}\n{result.output}')
+        output_lower = result.output.lower()
+        self.assertTrue(
+            'no' in output_lower or 'none' in output_lower
+            or 'crossing' in output_lower,
+            f'Expected no-crossing message, got: {result.output}'
+        )
+
+    def test_missing_options(self):
+        result = self.runner.invoke(
+            nonrad_cli,
+            ['barrier', '--dq', '1.0']  # missing --de, --wi, --wf
+        )
+        self.assertNotEqual(result.exit_code, 0)
+
+
+# ---------------------------------------------------------------------------
+# Tests for `capture`
+# ---------------------------------------------------------------------------
+class TestCapture(unittest.TestCase):
+    def setUp(self):
+        self.runner = CliRunner()
+
+    def test_help(self):
+        result = self.runner.invoke(nonrad_cli, ['capture', '--help'])
+        self.assertEqual(result.exit_code, 0, result.output)
+
+    @patch('nonrad.nonrad.get_C', return_value=1.23e-8)
+    def test_basic_single_temperature(self, mock_get_c):
+        result = self.runner.invoke(
+            nonrad_cli,
+            ['capture',
+             '--dq', '1.0', '--de', '1.0',
+             '--wi', '0.03', '--wf', '0.03',
+             '--wif', '0.01', '--volume', '1000',
+             '-T', '300']
+        )
+        self.assertEqual(result.exit_code, 0,
+                         f'exit_code={result.exit_code}\n{result.output}')
+
+    @patch('nonrad.nonrad.get_C',
+           return_value=np.array([1e-8, 2e-8, 3e-8]))
+    def test_temperature_range(self, mock_get_c):
+        result = self.runner.invoke(
+            nonrad_cli,
+            ['capture',
+             '--dq', '1.0', '--de', '1.0',
+             '--wi', '0.03', '--wf', '0.03',
+             '--wif', '0.01', '--volume', '1000',
+             '--temperature-range', '100', '500', '3']
+        )
+        self.assertEqual(result.exit_code, 0,
+                         f'exit_code={result.exit_code}\n{result.output}')
+
+    def test_missing_required_options(self):
+        result = self.runner.invoke(
+            nonrad_cli,
+            ['capture', '--dq', '1.0']  # missing many required options
+        )
+        self.assertNotEqual(result.exit_code, 0)
+
+
+# ---------------------------------------------------------------------------
+# Tests for `sommerfeld`
+# ---------------------------------------------------------------------------
+class TestSommerfeld(unittest.TestCase):
+    def setUp(self):
+        self.runner = CliRunner()
+
+    def test_help(self):
+        result = self.runner.invoke(nonrad_cli, ['sommerfeld', '--help'])
+        self.assertEqual(result.exit_code, 0, result.output)
+
+    def test_neutral_z0(self):
+        """Z=0 is the neutral special case: sommerfeld_parameter returns 1.0"""
+        result = self.runner.invoke(
+            nonrad_cli,
+            ['sommerfeld',
+             '--z', '0', '--m-eff', '1.0', '--eps0', '10.0',
+             '-T', '300']
+        )
+        self.assertEqual(result.exit_code, 0,
+                         f'exit_code={result.exit_code}\n{result.output}')
+        # Z=0 → sommerfeld = 1.0
+        self.assertIn('1.0', result.output)
+
+    def test_temperature_range(self):
+        result = self.runner.invoke(
+            nonrad_cli,
+            ['sommerfeld',
+             '--z', '0', '--m-eff', '1.0', '--eps0', '10.0',
+             '--temperature-range', '100', '500', '3']
+        )
+        self.assertEqual(result.exit_code, 0,
+                         f'exit_code={result.exit_code}\n{result.output}')
+
+    def test_missing_required_options(self):
+        result = self.runner.invoke(
+            nonrad_cli,
+            ['sommerfeld', '-T', '300']  # missing --z, --m-eff, --eps0
+        )
+        self.assertNotEqual(result.exit_code, 0)
+
+
+# ---------------------------------------------------------------------------
+# Tests for `thermal-velocity`
+# ---------------------------------------------------------------------------
+class TestThermalVelocity(unittest.TestCase):
+    def setUp(self):
+        self.runner = CliRunner()
+
+    def test_help(self):
+        result = self.runner.invoke(
+            nonrad_cli, ['thermal-velocity', '--help']
+        )
+        self.assertEqual(result.exit_code, 0, result.output)
+
+    def test_basic(self):
+        result = self.runner.invoke(
+            nonrad_cli,
+            ['thermal-velocity', '--m-eff', '1.0', '-T', '300']
+        )
+        self.assertEqual(result.exit_code, 0,
+                         f'exit_code={result.exit_code}\n{result.output}')
+        self.assertIn('cm/s', result.output)
+
+    def test_temperature_range(self):
+        result = self.runner.invoke(
+            nonrad_cli,
+            ['thermal-velocity', '--m-eff', '1.0',
+             '--temperature-range', '100', '500', '3']
+        )
+        self.assertEqual(result.exit_code, 0,
+                         f'exit_code={result.exit_code}\n{result.output}')
+
+    def test_missing_meff(self):
+        result = self.runner.invoke(
+            nonrad_cli,
+            ['thermal-velocity', '-T', '300']  # missing --m-eff
+        )
+        self.assertNotEqual(result.exit_code, 0)
+
+
+# ---------------------------------------------------------------------------
+# Tests for `charged-supercell`
+# ---------------------------------------------------------------------------
+class TestChargedSupercell(unittest.TestCase):
+    def setUp(self):
+        self.runner = CliRunner()
+
+    def test_help(self):
+        result = self.runner.invoke(
+            nonrad_cli, ['charged-supercell', '--help']
+        )
+        self.assertEqual(result.exit_code, 0, result.output)
+
+    @patch('nonrad.scaling.charged_supercell_scaling_VASP', return_value=0.95)
+    def test_basic_with_def_index(self, mock_scaling):
+        with self.runner.isolated_filesystem():
+            Path('WAVECAR').write_bytes(b'\x00' * 64)
+            result = self.runner.invoke(
+                nonrad_cli,
+                ['charged-supercell', 'WAVECAR',
+                 '--bulk-index', '189', '--def-index', '192']
+            )
+            self.assertEqual(result.exit_code, 0,
+                             f'exit_code={result.exit_code}\n{result.output}')
+
+    @patch('nonrad.scaling.charged_supercell_scaling_VASP', return_value=0.95)
+    def test_basic_with_def_coord(self, mock_scaling):
+        with self.runner.isolated_filesystem():
+            Path('WAVECAR').write_bytes(b'\x00' * 64)
+            result = self.runner.invoke(
+                nonrad_cli,
+                ['charged-supercell', 'WAVECAR',
+                 '--bulk-index', '189',
+                 '--def-coord', '0.5', '0.5', '0.5']
+            )
+            self.assertEqual(result.exit_code, 0,
+                             f'exit_code={result.exit_code}\n{result.output}')
+
+    def test_missing_def_index_and_def_coord(self):
+        """Both --def-index and --def-coord are omitted → error."""
+        with self.runner.isolated_filesystem():
+            Path('WAVECAR').write_bytes(b'\x00' * 64)
+            result = self.runner.invoke(
+                nonrad_cli,
+                ['charged-supercell', 'WAVECAR', '--bulk-index', '189']
+            )
+            # The CLI or underlying function should flag this as an error
+            self.assertNotEqual(result.exit_code, 0,
+                                f'Expected error, got: {result.output}')
+
+    def test_nonexistent_wavecar(self):
+        result = self.runner.invoke(
+            nonrad_cli,
+            ['charged-supercell', 'nonexistent_WAVECAR',
+             '--bulk-index', '189', '--def-index', '192']
+        )
+        self.assertNotEqual(result.exit_code, 0)
+
+
+# ---------------------------------------------------------------------------
+# Tests for `elphon` group
+# ---------------------------------------------------------------------------
+class TestElphonGroup(unittest.TestCase):
+    def setUp(self):
+        self.runner = CliRunner()
+
+    def test_help(self):
+        result = self.runner.invoke(nonrad_cli, ['elphon', '--help'])
+        self.assertEqual(result.exit_code, 0, result.output)
+        self.assertIn('elphon', result.output.lower())
+
+
+# ---------------------------------------------------------------------------
+# Tests for `elphon wavecars`
+# ---------------------------------------------------------------------------
+class TestElphonWavecars(unittest.TestCase):
+    def setUp(self):
+        self.runner = CliRunner()
+
+    def test_help(self):
+        result = self.runner.invoke(
+            nonrad_cli, ['elphon', 'wavecars', '--help']
+        )
+        self.assertEqual(result.exit_code, 0, result.output)
+
+    @patch('nonrad.elphon.get_Wif_from_wavecars',
+           return_value=[(189, 0.087)])
+    @patch('nonrad.ccd.get_Q_from_struct', return_value=0.5)
+    @patch('pymatgen.io.vasp.outputs.Vasprun')
+    @patch('pymatgen.core.structure.Structure.from_file')
+    def test_basic(self, mock_from_file, mock_vasprun,
+                   mock_get_q, mock_get_wif):
+        mock_struct = MagicMock()
+        mock_from_file.return_value = mock_struct
+
+        mock_vr = MagicMock()
+        mock_vr.structures = [mock_struct]
+        mock_vasprun.return_value = mock_vr
+
+        with self.runner.isolated_filesystem():
+            # Create fake directory structure
+            os.makedirs('cc/ground/0', exist_ok=True)
+            os.makedirs('cc/excited/0', exist_ok=True)
+            Path('cc/ground/0/WAVECAR').touch()
+            Path('cc/ground/0/vasprun.xml').touch()
+            os.makedirs('ground_files')
+            Path('ground_files/CONTCAR').touch()
+            os.makedirs('excited_files')
+            Path('excited_files/CONTCAR').touch()
+            Path('WAVECAR_init').touch()
+
+            result = self.runner.invoke(
+                nonrad_cli,
+                ['elphon', 'wavecars',
+                 'cc', 'ground_files', 'excited_files', 'WAVECAR_init',
+                 '--def-index', '192', '-b', '189']
+            )
+            self.assertEqual(result.exit_code, 0,
+                             f'exit_code={result.exit_code}\n{result.output}')
+
+
+# ---------------------------------------------------------------------------
+# Tests for `elphon wswq`
+# ---------------------------------------------------------------------------
+class TestElphonWswq(unittest.TestCase):
+    def setUp(self):
+        self.runner = CliRunner()
+
+    def test_help(self):
+        result = self.runner.invoke(
+            nonrad_cli, ['elphon', 'wswq', '--help']
+        )
+        self.assertEqual(result.exit_code, 0, result.output)
+
+    @patch('nonrad.elphon.get_Wif_from_WSWQ',
+           return_value=[(189, 0.094)])
+    @patch('nonrad.ccd.get_Q_from_struct', return_value=0.5)
+    @patch('pymatgen.io.vasp.outputs.Vasprun')
+    @patch('pymatgen.core.structure.Structure.from_file')
+    def test_basic(self, mock_from_file, mock_vasprun,
+                   mock_get_q, mock_get_wif):
+        mock_struct = MagicMock()
+        mock_from_file.return_value = mock_struct
+
+        mock_vr = MagicMock()
+        mock_vr.structures = [mock_struct]
+        mock_vasprun.return_value = mock_vr
+
+        with self.runner.isolated_filesystem():
+            os.makedirs('cc/ground/0', exist_ok=True)
+            os.makedirs('cc/excited/0', exist_ok=True)
+            Path('cc/ground/0/WSWQ').touch()
+            Path('cc/ground/0/vasprun.xml').touch()
+            os.makedirs('ground_files')
+            Path('ground_files/CONTCAR').touch()
+            Path('ground_files/vasprun.xml').touch()
+            os.makedirs('excited_files')
+            Path('excited_files/CONTCAR').touch()
+            Path('excited_files/vasprun.xml').touch()
+
+            result = self.runner.invoke(
+                nonrad_cli,
+                ['elphon', 'wswq',
+                 'cc', 'ground_files', 'excited_files',
+                 'ground_files/vasprun.xml',
+                 '--def-index', '192', '-b', '189']
+            )
+            self.assertEqual(result.exit_code, 0,
+                             f'exit_code={result.exit_code}\n{result.output}')
+
+
+# ---------------------------------------------------------------------------
+# Tests for `elphon unk`
+# ---------------------------------------------------------------------------
+class TestElphonUnk(unittest.TestCase):
+    def setUp(self):
+        self.runner = CliRunner()
+
+    def test_help(self):
+        result = self.runner.invoke(
+            nonrad_cli, ['elphon', 'unk', '--help']
+        )
+        self.assertEqual(result.exit_code, 0, result.output)
+
+    @patch('nonrad.elphon.get_Wif_from_UNK',
+           return_value=[(1, 0.5)])
+    def test_basic(self, mock_get_wif):
+        init_unk = str(TEST_FILES / 'UNK.0')
+        result = self.runner.invoke(
+            nonrad_cli,
+            ['elphon', 'unk', init_unk,
+             '--def-index', '2', '-b', '1',
+             '--eigs', '0.0,1.0',
+             '-u', '1.0', str(TEST_FILES / 'UNK.1')]
+        )
+        self.assertEqual(result.exit_code, 0,
+                         f'exit_code={result.exit_code}\n{result.output}')
+
+
+if __name__ == '__main__':
+    unittest.main()
diff --git a/pyproject.toml b/pyproject.toml
index 8926ed0..b896cfb 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -53,6 +53,7 @@ dependencies = [
     "monty>=2023.4.10",
     "numba>=0.57.0",
     "mpmath>=1.3.0",
+    "click>=8.0.1",
 ]
 
 [project.urls]
@@ -83,3 +84,6 @@ extend-exclude = ["docs"]
 [tool.mypy]
 ignore_missing_imports = "True"
 plugins = "numpy.typing.mypy_plugin"
+
+[project.scripts]
+nonrad = "nonrad.cli:nonrad"
\ No newline at end of file