TD Coupling Model

A Python package for computing electromagnetic coupling between molecular and plasmonic excitations using transition densities and an FFT-based Coulomb solver.

td-coupling-model builds and diagonalizes the polariton Hamiltonian of coupled molecule–plasmon systems, enabling quantitative analysis of strong-coupling effects such as Rabi splitting and hybrid polariton states. Transition densities must be provided as Gaussian .cube files, as typically produced by TDDFT calculations.

---

Theoretical Background

The coupling between a molecular electronic transition and a plasmonic excitation is described by the following model Hamiltonian:

$$ H = \begin{pmatrix} \omega_\text{mol} & g \\ g & \omega_\text{plas} \end{pmatrix} $$

which generalizes, for $N_m$ molecules and $N_p$ plasmons, to a block matrix:

$$ H = \begin{pmatrix} H_\text{mol} & g \\ g^T & H_\text{plas} \end{pmatrix} $$

where $H_\text{mol}$ and $H_\text{plas}$ contain the bare excitation energies on the diagonal and the molecule–molecule ($V_{ij}$) or plasmon–plasmon ($J_{ij}$) couplings as off-diagonal terms.

The coupling constant $g$ between species A and B is evaluated via the Coulomb integral over their real-space transition densities:

$$ g = \int \rho_A(\mathbf{r}), \Phi_B(\mathbf{r}), dV $$

where $\Phi_B$ is the electrostatic potential generated by $\rho_B$, obtained by solving the Poisson equation in Fourier space:

$$ \tilde{\Phi}_B(\mathbf{k}) = \frac{\tilde{\rho}_B(\mathbf{k})}{\varepsilon_0, \varepsilon_\text{eff}, k^2} $$

Diagonalization of $H$ yields the polariton eigenenergies and Hopfield coefficients $|C_i|^2$, which quantify the molecular vs. plasmonic character of each hybrid state. Strong coupling is identified by the appearance of Rabi splitting in the eigenspectrum, when $g$ exceeds the average linewidth of the coupled transitions.

---

Installation

git clone https://github.com/TUO-USER/TD_coupling_model.git
cd TD_coupling_model
pip install .

Requirements: Python ≥ 3.9, numpy ≥ 1.23, scipy ≥ 1.10, periodictable ≥ 2.1.0

---

Package Structure

src/tdcouplingmodel/
├── io/               # Cube file I/O
├── core/             # Coupling calculation and Hamiltonian
├── preprocessing/    # Grid alignment
├── utils/            # Rotation, summation, conversion, constants
└── cli/              # Command-line interfaces
tests                 # An example of how the code works

---

Module Reference

io — Cube File I/O

Reads and writes Gaussian .cube files, the standard volumetric format for quantum-chemistry data. Unit detection (Bohr vs. Ångström) is automatic from the sign of the grid vectors in the file header.

• parse_cube_file(path) — parses a .cube file and returns the 3D density grid, atomic coordinates, grid vectors, origin, and metadata
• write_cube_file(path, header, density, ...) — writes a 3D density array back to .cube format
• density_3d_from_flat(flat, shape) — reshapes a flat density array into a 3D (NX, NY, NZ) grid

---

core — Coupling & Hamiltonian

The computational core of the package.

• compute_g_fft_coulomb(rho_A, rho_B, steps, units_len, ...) — computes $g$ in eV between two transition densities using the FFT Poisson solver. Supports tunable zero-padding, physical margin extension, and optional bounding-box cropping to reduce memory usage on large grids.
• build_hamiltonian(omega_m, omega_p, g_matrix, V_matrix, J_matrix) — assembles the full $(N_m + N_p) \times (N_m + N_p)$ polariton Hamiltonian from bare energies and all pairwise coupling matrices.
• diagonalize_hamiltonian(H) — diagonalizes the Hamiltonian via numpy.linalg.eigh and returns sorted polariton energies and eigenvectors.
• analyze_polariton_character(states, N_m, N_p) — computes Hopfield coefficients $|C_i|^2$ for each polariton state, returning the percentage molecular and plasmonic character.

---

preprocessing — Grid Alignment

Before computing the coupling integral, all transition density grids must share the same spatial support and voxel spacing.

• align_cube_grids(input_files, output_files, spacing, order, ...) — resamples all input cube files onto a common grid defined by their union bounding box. The target voxel spacing can be set to the minimum, maximum, or mean spacing across the inputs, or specified as a fixed value. Spline interpolation order is configurable (0–5).

---

utils — Utilities

• roto_traslation — applies rigid-body transformations (rotation + translation) to a .cube file, rigidly moving both atomic coordinates and the volumetric density. Rotations can be specified as Euler angles (ZYX convention) or derived automatically from two transition dipole vectors: given the plasmonic dipole $\boldsymbol{\mu}_1$ and the molecular dipole $\boldsymbol{\mu}_2$, the code computes and applies the rotation that aligns $\boldsymbol{\mu}_2$ onto $\boldsymbol{\mu}_1$, maximising the dipole–dipole coupling.
• sum_cubes — sums two .cube files on a common grid. Each input file can be multiplied by a scalar coefficient before summation, producing a weighted linear combination of the form $c_1 \rho_1 + c_2 \rho_2$. Atomic coordinates from both files are merged in the output. Grid compatibility (size, origin, voxel spacing, units) is verified automatically.
• conversion_cube_to_xyz — extracts atomic coordinates from a .cube file and writes them as an .xyz file (always in Ångström) for visualization in VESTA, Avogadro, VMD, or similar tools.
• constants — physical constants used throughout the package: E_CHARGE, EPS0, BOHR_TO_M, ANG_TO_M, BOHR_TO_ANG, ANG_TO_BOHR.

---

Command-Line Interface

All operations are available as standalone terminal commands after installation.

---

tdcoupling-calc — Compute Coupling & Diagonalize

Computes all pairwise coupling constants, builds the full polariton Hamiltonian, diagonalizes it, and saves results to a text file. The output includes polariton energies (eV and nm), Hopfield coefficients $|C_i|^2$, and raw eigenvectors.

tdcoupling-calc \
  --mol-cubes  mol1.cube mol2.cube \
  --mol-energies 2.5 2.6 \
  --plas-cubes plas1.cube plas2.cube \
  --plas-energies 2.4 2.5 \
  --output results

Option	Default	Description
--mol-cubes	(required)	Molecular transition density cube files
--mol-energies	(required)	Molecular excitation energies
--plas-cubes	(required)	Plasmon transition density cube files
--plas-energies	(required)	Plasmon excitation energies
--energy-units	ev	Energy units: ev or nm
--output	coupling_results	Output file prefix
--eps-eff	1.0	Effective dielectric constant $\varepsilon_\text{eff}$
--pad	1.0	FFT zero-padding factor
--margin	0.0	Extra FFT margin in bohr/Å
--use-crop	off	Crop densities to bounding box before FFT
--crop-threshold	1e-6	Density threshold for bounding-box cropping
--crop-margin	2.0	Physical safety margin around bounding box (bohr/Å)
--disable-mol-mol	off	Set all molecule–molecule couplings $V = 0$
--disable-plas-plas	off	Set all plasmon–plasmon couplings $J = 0$
-q, --quiet	off	Suppress verbose output

---

tdcoupling-align — Grid Alignment

Resamples two or more cube files onto a common spatial grid. This step is required before running tdcoupling-calc whenever the input files have different grids.

tdcoupling-align \
  -i mol.cube plas.cube \
  -o mol_aligned.cube plas_aligned.cube \
  --spacing min --verify

Option	Default	Description
-i, --input	(required)	Input cube files (≥ 2)
-o, --output	(required)	Output cube files (same count as input)
--spacing	min	Target grid spacing: min, max, mean, or numeric value
--order	3	Spline interpolation order (0–5)
--verify	off	Check charge conservation after resampling
-q, --quiet	off	Suppress output messages

---

tdcoupling-rotate — Rotation, Translation & Dipole Alignment

Applies a rigid-body transformation to a .cube file, moving both atomic coordinates and volumetric density. Rotations are specified as Euler angles (ZYX convention in degrees). Translations are in Bohr or Ångström depending on --input_unit.

Alternatively, providing two transition dipole vectors via --dipole_1 and --dipole_2 computes the rotation that aligns the molecular dipole $\boldsymbol{\mu}_2$ onto the plasmonic dipole $\boldsymbol{\mu}_1$, maximising the dipole-orientation contribution to the coupling $g$.

# Rotate by Euler angles and translate
tdcoupling-rotate mol.cube mol_rot.cube \
  --angles 0 45 0 --trans 5.0 0.0 0.0 --input_unit ang

# Align molecular dipole onto plasmonic dipole
tdcoupling-rotate mol.cube mol_rot.cube \
  --dipole_1 0 0 1 --dipole_2 1 0 0

Option	Default	Description
cube_in	(required)	Input cube file
cube_out	(required)	Output cube file
--angles α β γ	0 0 0	Euler angles in degrees (ZYX order)
--trans TX TY TZ	0 0 0	Translation vector
--dipole_1 x y z	—	Target dipole vector $\boldsymbol{\mu}_1$ (for example plasmon)
--dipole_2 x y z	—	Source dipole vector $\boldsymbol{\mu}_2$ (for example molecule, to be rotated)
--center CX CY CZ	grid center	Rotation center
--input_unit	ang	Units for translation and center: ang or bohr
-q, --quiet	off	Suppress output messages

---

tdcoupling-sum — Weighted Cube Summation

Sums two .cube files on a common grid, with optional scalar coefficients. The result is $c_1 \cdot \rho_1 + c_2 \cdot \rho_2$, written as a new .cube file with merged atomic coordinates. Grid compatibility is verified automatically.

# Simple sum
tdcoupling-sum mode1.cube mode2.cube combined.cube

# Weighted sum: 0.7 × mode1 + 0.3 × mode2
tdcoupling-sum mode1.cube mode2.cube combined.cube \
  --factor1 0.7 --factor2 0.3

Option	Default	Description
file1	(required)	First input cube file
file2	(required)	Second input cube file
output	(required)	Output cube file
--factor1	1.0	Scalar coefficient for the first cube
--factor2	1.0	Scalar coefficient for the second cube
--verify	on	Re-read output file to verify write integrity
-q, --quiet	off	Suppress output messages

---

tdcoupling-cube2xyz — Cube to XYZ Conversion

Extracts atomic coordinates from a .cube file and writes them as an .xyz file. Units are auto-detected from the sign of the grid vectors (positive → Bohr, negative → Ångström); output is always in Ångström.

tdcoupling-cube2xyz molecule.cube molecule.xyz

---

Quick Start Example

# 1. Align grids
tdcoupling-align \
  -i mol.cube plas.cube \
  -o mol_al.cube plas_al.cube --verify

# 2. Compute coupling and diagonalize
tdcoupling-calc \
  --mol-cubes mol_al.cube --mol-energies 2.5 \
  --plas-cubes plas_al.cube --plas-energies 2.5 \
  --output results

Python API

from tdcouplingmodel import (
    parse_cube_file, align_cube_grids,
    compute_g_fft_coulomb, build_hamiltonian,
    diagonalize_hamiltonian, analyze_polariton_character,
)
import numpy as np

mol  = parse_cube_file("mol_aligned.cube")
plas = parse_cube_file("plas_aligned.cube")

g = compute_g_fft_coulomb(
    mol["density"], plas["density"],
    steps=mol["steps"], units_len="bohr"
)
print(f"g = {g:.4f} eV  (2g = {2*g:.4f} eV)")

H = build_hamiltonian(
    omega_m=np.array([2.5]),
    omega_p=np.array([2.5]),
    g_matrix=np.array([[g]]),
    V_matrix=np.zeros((1, 1)),
    J_matrix=np.zeros((1, 1)),
)
energies, states = diagonalize_hamiltonian(H)
print(f"Rabi splitting: {energies - energies:.4f} eV") [reddit](https://www.reddit.com/r/learnprogramming/comments/vxfku6/how_to_write_a_readme/)

character = analyze_polariton_character(states, N_m=1, N_p=1)
for i, c in enumerate(character):
    print(f"State {i+1}: {c['mol_character']:.1f}% mol  "
          f"{c['plas_character']:.1f}% plas")

---

License

GNU General Public License v3.0 or later (GPLv3+). See LICENSE for details.

Contact

Lucia Cascino – lucia.cascino@unisalento.it
Issues: https://github.com/luciacasc/TD_model/issues