qscat.dvr

The FEM-DVR grid with an exterior-complex-scaled tail: grid construction, kinetic-energy assembly (dense and sparse), and the N-dimensional tensor layer. See docs/physics/femdvr-ecs.md and docs/physics/nd-tensor-hamiltonian.md.

ElementSpec	One finite element of the radial grid.
GridSpec	Full grid specification: quadrature order (shared by all elements), the ordered list of elements, and the inner boundary x_min.
FemDvrEcsGrid	FEM-DVR-ECS grid geometry built from a validated GridSpec.
kinetic	Assemble the (n, n) complex FEM-DVR kinetic-energy matrix (dense).
kinetic_sparse	Sparse (CSR) FEM-DVR kinetic-energy matrix -- the sparse sibling of kinetic.
dvr_first_derivative_at_node	Row d (length grid.n) s.t.
dvr_interpolation_matrix	Sparse operator P s.t.
hamiltonian	Assemble H = T + diag(V) on the FEM-DVR basis (diagonal-potential approx).
eigen	Eigendecompose H, sorted by ascending Re(E).
gll_nodes_weights	Gauss-Lobatto-Legendre nodes and weights on (-1, 1), n points incl.
diff_matrix	Collocation differentiation matrix D with (D @ f)[j] = f'(x_j) for f sampled at x.
TensorGrid	Tensor product of D FEM-DVR-ECS grids (C order: last axis fastest).
kinetic_nd	sum_d I x … x T_d x … x I, with T_d built at mass masses[d].
potential_nd	Evaluate V at the D-dimensional COMPLEX points, flattened.
hamiltonian_nd	H = kinetic_nd(tgrid, masses) + diag(potential_nd(tgrid, V)) as CSR.

FEM-DVR-ECS: finite-element DVR radial grid with exterior complex scaling.

Public API:

• ElementSpec, GridSpec – grid specification dataclasses.

• FemDvrEcsGrid – grid geometry (points, bridge-summed weights, the local-to-global element_maps) built from a validated GridSpec.

• kinetic, kinetic_sparse – assemble the FEM-DVR kinetic-energy matrix on a grid, dense or sparse (CSR). The dense one is the sparse one’s differential oracle.

• dvr_first_derivative_at_node – the FEM-DVR first-derivative operator’s row at a real grid node (d @ psi_coeffs ~= d/dx psi(x_node)), built from the same element-local Lagrange-derivative matrix kinetic uses; the flow (flux) time-dependent energy extractor’s numerical primitive.

• hamiltonian, eigen – diagonal-potential Hamiltonian assembly and the complex-symmetric (non-Hermitian) eigensolver.

• gll_nodes_weights, diff_matrix – the underlying Gauss-Lobatto- Legendre quadrature/differentiation building blocks (reusable outside FEM-DVR-ECS).

• TensorGrid – tensor product of D FEM-DVR-ECS grids (C order, last axis fastest), with the ECS real-region mask and separable-state construction.

• kinetic_nd, potential_nd, hamiltonian_nd – the N-dimensional Kronecker-sum Hamiltonian assembled on a TensorGrid, sparse (CSR).

See docs/physics/femdvr-ecs.md for the method and its validation benchmarks. The method is Rescigno & McCurdy, Phys. Rev. A 62, 032706 (2000); this implementation was ported from eMoScat.

class qscat.dvr.ElementSpec(length, angle_deg=0.0)

One finite element of the radial grid.

length is the (real, positive) physical length of the element along the unscaled radial coordinate. angle_deg is the exterior-complex-scaling rotation angle (degrees); 0.0 means a real (unscaled) element.

Parameters:

• length (float)

• angle_deg (float)

class qscat.dvr.FemDvrEcsGrid(spec)

FEM-DVR-ECS grid geometry built from a validated GridSpec.

See the module docstring for the element_maps local->global index convention that the kinetic-operator assembly relies on.

Parameters:

spec (GridSpec)

real_index_near(r_value)

Index of the real-region DVR node nearest r_value (bohr).

The ECS complex-tail points (unscaled coordinate > R0) are masked out first, so the returned index is always a genuine real-region node – the point a fixed-position/fixed-surface time-dependent extractor (Dirac/Flux) or a wavepacket launch reads at. Shared primitive for what was copy-pasted as _index_near/_real_index_near/ _nuclear_index_near across the runner and the TD validation drivers.

Parameters:

r_value (float)

Return type:

int

class qscat.dvr.GridSpec(quadrature, elements, x_min=0.0)

Full grid specification: quadrature order (shared by all elements), the ordered list of elements, and the inner boundary x_min.

Validates that complex (ECS) elements form a contiguous tail at the end of the element list – eMoScat’s complex_negative (ECS at the inner boundary) is deliberately not supported here. Computes the ECS pivot R0 = x_min + sum(real element lengths), which by construction sits exactly on an element boundary.

Caveat: using multiple different nonzero angle_deg values across tail elements (a bent/graded ECS contour) is REJECTED because it has never been validated here – the validated, actually-used case is a single ECS tail angle shared by all complex elements.

Parameters:

• quadrature (int)

• elements (list[ElementSpec])

• x_min (float)

class qscat.dvr.TensorGrid(grids)

Tensor product of D FEM-DVR-ECS grids (C order: last axis fastest).

Parameters:

grids (Sequence[FemDvrEcsGrid])

property grids: tuple[FemDvrEcsGrid, ...]

The per-dimension FemDvrEcsGrid objects, in tensor order.

property ndim: int

Number of tensor dimensions D.

property shape: tuple[int, ...]

Per-dimension point counts (n_1, …, n_D).

property size: int

Total number of tensor-product points, prod(shape).

points()

The D complex (ECS-scaled) coordinate arrays, shaped to broadcast.

For D = 2 these are (n_0, 1) and (1, n_1), so a potential can be written naturally as V(r, R) without materializing a full meshgrid. For D = 1 the single array is plain 1-D of shape (n_0,).

Return type:

tuple[NDArray[complex128], …]

weights()

The D bridge-summed, COMPLEX (ECS-scaled) quadrature weight arrays, shaped to broadcast exactly like points() (same _broadcast_shape convention: axis d’s weight varies along axis d, size 1 elsewhere).

Converting a function f to FEM-DVR basis coefficients (not evaluating an already-DVR-diagonal potential, which is what potential_nd does) requires c_j = f(x_j) * sqrt(w_j), using the GLOBAL, bridge-summed, COMPLEX weight FemDvrEcsGrid.weights – qscat.dvr.kinetic’s module docstring calls getting this wrong (using a real or per-element weight instead) “the classic assembly trap”. TensorGrid already spares a caller from remembering real_mask(); the same principle applies here, since a driving term or an asymptotic channel function must go through this exact factor at every coordinate to land correctly on the pre-normalized DVR basis. See also sqrt_weights(), which supplies that factor directly.

Return type:

tuple[NDArray[complex128], …]

sqrt_weights()

sqrt(weights()) per axis – the exact factor a basis-coefficient conversion c_j = f(x_j) * sqrt(w_j) needs, so a caller never has to write np.sqrt on the bridge-summed complex weight itself (a detail easy to get subtly wrong under ECS, since the weight is complex and sqrt of a complex number requires choosing a branch – NumPy’s principal branch, Re(sqrt(z)) >= 0, is the one used here and is the correct one for this weight, which never crosses the negative real axis for a valid ECS tail angle).

Return type:

tuple[NDArray[complex128], …]

real_mask()

Flat boolean mask, True where EVERY coordinate is in the unscaled region.

Under exterior complex scaling a driving term or channel projection is only meaningful on the unscaled region, so anything of that kind must be masked with this before use. Making it a property of the grid is deliberate: the physics layer should not have to remember.

Return type:

NDArray[bool]

outer(vectors)

Separable state ⊗_d vectors[d], flattened to length size.

Parameters:

vectors (Sequence[ArrayLike])

Return type:

NDArray[complex128]

qscat.dvr.diff_matrix(x)

Collocation differentiation matrix D with (D @ f)[j] = f’(x_j) for f sampled at x.

D[j, i] = L_i’(x_j), via barycentric weights (robust for any node set).

Parameters:

x (NDArray[float64])

Return type:

NDArray[float64]

qscat.dvr.dvr_first_derivative_at_node(grid, node_index)

Row d (length grid.n) s.t. d @ psi_coeffs ~= d/dx psi(x_node).

node_index must be a real (unscaled) grid node (grid.real_points[node_index] <= grid.R0) – see module docstring.

Parameters:

• grid (FemDvrEcsGrid)

• node_index (int)

Return type:

NDArray[complex128]

qscat.dvr.dvr_interpolation_matrix(grid, x)

Sparse operator P s.t. (P @ state)[k] is the field value at x[k].

Parameters:

• grid (FemDvrEcsGrid) – The FEM-DVR-ECS radial grid the state is represented on.

• x (array_like) – Real sample coordinates (unscaled), each within [x_min, x_max].

Returns:

Complex matrix of shape (len(x), grid.n). P @ state gives the state’s values at x; P has at most nq nonzeros per row (the nodes of the containing element) and is built once, reused per frame.

Return type:

scipy.sparse.csr_matrix

Raises:

GridError – If any sample point lies outside [x_min, x_max].

qscat.dvr.eigen(H)

Eigendecompose H, sorted by ascending Re(E).

Eigenvectors carry numpy’s v†v=1 (Hermitian) normalization; for ECS c-product observables (resonance residues and widths), re-normalize to vᵀv=1.

Parameters:

H (NDArray[complex128])

Return type:

tuple[NDArray[complex128], NDArray[complex128]]

qscat.dvr.gll_nodes_weights(n)

Gauss-Lobatto-Legendre nodes and weights on (-1, 1), n points incl. endpoints.

Interior nodes are the roots of P’_{n-1}; endpoints are +/-1. Weights: w_i = 2 / (n(n-1) [P_{n-1}(x_i)]^2).

Parameters:

n (int)

Return type:

tuple[NDArray[float64], NDArray[float64]]

qscat.dvr.hamiltonian(grid, V, mass)

Assemble H = T + diag(V) on the FEM-DVR basis (diagonal-potential approx).

Parameters:

• grid (FemDvrEcsGrid)

• V (Callable[[NDArray[complex128]], ArrayLike] | ArrayLike)

• mass (float)

Return type:

NDArray[complex128]

qscat.dvr.hamiltonian_nd(tgrid, masses, V)

H = kinetic_nd(tgrid, masses) + diag(potential_nd(tgrid, V)) as CSR.

Complex symmetric (H = H^T), NOT Hermitian, whenever any grid has an ECS tail. Use general algorithms only.

Parameters:

• tgrid (TensorGrid)

• masses (Sequence[float])

• V (Callable[[...], ArrayLike])

Return type:

csr_matrix

qscat.dvr.kinetic(grid, mass)

Assemble the (n, n) complex FEM-DVR kinetic-energy matrix (dense).

Parameters:

• grid (FemDvrEcsGrid)

• mass (float)

Return type:

NDArray[complex128]

qscat.dvr.kinetic_nd(tgrid, masses)

sum_d I x … x T_d x … x I, with T_d built at mass masses[d].

Parameters:

• tgrid (TensorGrid)

• masses (Sequence[float])

Return type:

csr_matrix

qscat.dvr.kinetic_sparse(grid, mass)

Sparse (CSR) FEM-DVR kinetic-energy matrix – the sparse sibling of kinetic.

Identical mathematics to the dense kinetic(), which is retained as this function’s differential oracle. The only structural difference is that per-element blocks are emitted as COO triplets instead of scatter-added into a dense array: coo_matrix SUMS duplicate (row, col) entries on conversion, which reproduces the dense version’s += bridge accumulation exactly. No bridge special-casing is needed or wanted.

Nonzero count is nq**2 * tnel - 4*nq + 3 - tnel (eMoScat KineticEnergy.cpp:95) – the union of the per-element nq x nq blocks, overlapping by one index at each bridge, less the two dropped Dirichlet endpoints.

Parameters:

• grid (FemDvrEcsGrid)

• mass (float)

Return type:

csr_matrix

qscat.dvr.potential_nd(tgrid, V)

Evaluate V at the D-dimensional COMPLEX points, flattened.

V is called as V(x_0, …, x_{D-1}) with the broadcastable arrays from TensorGrid.points(). It MUST NOT coerce its arguments to a real dtype: the points are complex on the ECS tail, and discarding the imaginary part silently destroys the analytic continuation the method depends on.

TWO BROADCASTING TRAPS to know about, because np.broadcast_to will happily paper over both rather than raising:

• A RANK-DEFICIENT return value. If V returns an array with fewer axes than tgrid.ndim (e.g. a plain (n,) array at D=2, perhaps because V evaluated only one coordinate, or returned a precomputed 1-D array by mistake), NumPy broadcasting aligns it against the LAST axis and tiles it silently along every other axis – a wrong potential with a plausible shape, not a crash. This function guards against exactly that: any non-scalar V result whose ndim is not tgrid.ndim raises ValueError. A genuine scalar (ndim == 0, a spatially constant V) is fine and is deliberately exempt, since a true scalar broadcasts unambiguously to every axis.

• AXIS TRANSPOSITION when two axes have the SAME size. This one is NOT caught by any shape check, because the shapes are literally equal: np.meshgrid defaults to indexing=”xy” (the first two returned arrays are transposed relative to the “ij”/broadcastable convention TensorGrid.points() uses), so a V that builds its own meshgrid internally – or receives an externally precomputed surface of shape (n1, n0) instead of (n0, n1) – silently transposes the potential whenever n0 == n1. There is no shape-based guard for this: call V(*tgrid.points()) (or pass indexing=”ij” explicitly to any internal np.meshgrid call) rather than reconstructing the grid by hand.

Parameters:

• tgrid (TensorGrid)

• V (Callable[[...], ArrayLike])

Return type:

NDArray[complex128]