adf2stowf — ADF to CASINO STO Wavefunction Converter

Overview

This module converts a Slater-Type Orbital (STO) wavefunction from an ADF (Amsterdam Density Functional) TAPE21.asc file into the stowfn.data format used by the CASINO quantum Monte Carlo code.

The conversion pipeline covers:

  • Parsing geometry, basis set (valence + core), and MO coefficients from the ADF output via adfread.AdfParser;

  • Remapping ADF’s Cartesian angular basis functions to the real solid-harmonic (spherical) polynomial basis expected by CASINO;

  • Optionally projecting out non-spherical contamination in d/f shells (--cart2harm-projection);

  • Enforcing nuclear cusp conditions on the resulting MO coefficients through one of three strategies: enforce, project, or none;

  • Writing a ready-to-use stowfn.StoWfn file and, optionally, plotting the wavefunction cusp behaviour per atom.

The module is intended to be run as a command-line script. The ADFToStoWF class encapsulates all conversion logic and can also be instantiated programmatically.

ADFToStoWF class

class adf2stowf.ADFToStoWF(plot_cusps, cusp_method, do_dump, cart2harm_projection, only_occupied)

Top-level converter object. Parses TAPE21.asc on construction and exposes the full conversion pipeline through run().

Parameters:

  • plot_cusps (bool) – When True, evaluate MO values and Laplacians along a z-axis line through each nucleus before and after cusp correction, and save a cusp_constraint.svg plot.

  • cusp_method (str) – Strategy for enforcing nuclear cusp conditions. One of 'enforce', 'project', or 'none' (see apply_cusp_correction()).

  • do_dump (bool) – When True, write a human-readable TAPE21.txt dump of the parsed ADF data (via adfread.AdfParser.write_dump()).

  • cart2harm_projection (bool) – When True, project MO coefficients onto the pure spherical-harmonic subspace before the Cartesian-to-spherical transformation, removing s-type contamination in d/f shells.

  • only_occupied (bool) – When True, include only occupied MOs in the output (occupation >= 2/Nspins). When False, virtual orbitals are also written.

Key instance attributes

General / metadata (from the General section of TAPE21)

Nspins: int

Number of spin channels (1 = restricted, 2 = unrestricted).

spin_restricted: bool

True when Nspins == 1.

Nvalence_electrons: int

Total number of valence electrons.

Natoms: int

Number of real (non-dummy) atoms.

Natomtypes: int

Number of distinct atom types.

atyp_idx: numpy.ndarray, shape(Natoms)

Zero-based atom-type index for each atom; maps atom positions to their corresponding atom type in the basis tables.

Valence basis (populated by process_valence_basis())

nbset: int

Total number of valence STO shells across all atom types.

valence_shelltype: numpy.ndarray

Shell-type code for every valence shell (CASINO encoding: 1 = s, 2 = sp, 3 = p, 4 = d, 5 = f).

valence_order_r: numpy.ndarray

Radial power N for every valence shell (nqbas - lqbas - 1).

valence_zeta: numpy.ndarray

Slater exponents for every valence shell.

valence_cartnorm: numpy.ndarray

Cartesian normalization factors (bnorm) for every valence Cartesian basis function.

Core basis (populated by process_core_basis())

ncset: int

Total number of frozen-core STO shells across all atom types.

core_shelltype: numpy.ndarray

Shell-type codes for every core shell.

core_order_r: numpy.ndarray

Radial power N for every core shell.

core_zeta: numpy.ndarray

Slater exponents for every core shell.

Transformation matrices (built in initialize_data())

harm2cart_map: dict[int, numpy.ndarray]

Mapping from shell-type code to the matrix M that converts spherical AO coefficients (rows) to Cartesian AO coefficients (columns). Available for s (1), p (3), d (4), and f (5) shells; d and f matrices follow the CASINO stowfdet.f90 polynomial ordering convention.

cart2harm_map: dict[int, numpy.ndarray]

Inverse of harm2cart_map (computed via numpy.linalg.inv).

cart2harm_matrix: numpy.ndarray, shape(Nharmbasfns, Nvalence_cartbasfn)

Block-diagonal matrix that converts valence MO coefficients from the Cartesian AO basis (as stored in TAPE21) to the real solid-harmonic basis required by CASINO.

cart2harm_constraint: numpy.ndarray

Rows of the harmonic-to-Cartesian matrices corresponding to non-spherical (contaminating) components (e.g. the s-contamination in a d shell). Used to detect and optionally project out violations.

MO coefficients (populated by process_coefficients() and finalize_coefficients())

molorb_cart_coeff: list of numpy.ndarray

molorb_cart_coeff[spin] — valence MO coefficients in the Cartesian AO basis, shape (Nvalence_molorbs[spin], Nvalence_cartbasfn), sorted by eigenvalue.

valence_molorb_harm_coeff: list of numpy.ndarray

Valence MO coefficients after the Cartesian-to-spherical transformation, shape (Nharmbasfns, Nvalence_molorbs[spin]) per spin.

coeff: list of numpy.ndarray

Final full MO coefficient matrices (core + valence combined, spherical basis), shape (Nharmbasfns, Nmolorbs[spin]) per spin.

sto: stowfn.StoWfn

The stowfn.StoWfn object constructed in setup_stowfn() and written to stowfn.data at the end of apply_cusp_correction().

Methods

initialize_data()

Extract the primary data sections from the parsed TAPE21 dictionary and compute derived metadata.

Populates General, Geometry, Basis, Core, and Symmetry section references, scalar counts (Nspins, Natoms, …), the atom-type index atyp_idx, and the harm2cart_map / cart2harm_map transformation matrices for s, p, d, and f shells. Optionally writes TAPE21.txt when DO_DUMP is True.

Raises:

AssertionError – If atom or atom-type counts in the file are internally inconsistent.

process_valence_basis()

Parse the Basis section of TAPE21 and populate all valence-basis attributes.

Reads shell counts (nbset, nbaspt), quantum numbers (nqbas, lqbas), Slater exponents (alfbas), Cartesian function counts (nbptr), and normalization factors (bnorm). Derives per-centre and per-atom-type slices for all of the above.

Raises:

AssertionError – If any basis dimension is inconsistent with the geometry or the total number of AOs (naos).

process_core_basis()

Parse the Core section of TAPE21 and populate all frozen-core-basis attributes.

Reads shell counts (ncset, ncorpt, nrcset), quantum numbers (nqcor, lqcor), Slater exponents (alfcor), and per-shell normalization factors (cornrm). Expands normalization values into per-basis-function arrays for s and p core shells.

Raises:

  • ValueError – If d- or f-type frozen core shells are encountered (not yet implemented).

  • AssertionError – If core-shell counts do not match the nrcset table.

process_shells()

Merge core and valence shell data into per-centre arrays.

Concatenates core shells (from process_core_basis()) before valence shells (from process_valence_basis()) for each atom, producing shelltype_per_centre, order_r_per_centre, and zeta_per_centre.

select_coeff(sp)

Collect and sort MO coefficients for spin channel sp from the symmetry-resolved Eigen-Bas blocks in TAPE21.

Iterates over all symmetry labels (symlab), reads fractional occupations (froc_A / froc_B), eigenvalues (eps_A / eps_B), and coefficient matrices (Eigen-Bas_A / Eigen-Bas_B). Handles partial occupations by accumulating leftover occupation at degenerate eigenvalues. Returns orbitals sorted by eigenvalue.

Prints a Warning: HOMO > LUMO message if the highest occupied eigenvalue exceeds the lowest unoccupied eigenvalue, which can occur for open-shell or fractional-occupation calculations.

Parameters:

sp (int) – Spin channel index (0 = alpha, 1 = beta).

Returns:

MO coefficient matrix in the Cartesian AO basis.

Return type:

numpy.ndarray, shape (Nmolorbs_selected, Nvalence_cartbasfn)

process_coefficients()

Build the Cartesian-to-spherical transformation, detect angular contamination, optionally project it out, and transform MO coefficients to the spherical basis.

Constructs the block-diagonal cart2harm_matrix and the cart2harm_constraint from the per-shell cart2harm_map entries. Checks the constraint violation norm per MO and prints a WARNING for any orbital exceeding 1e-5.

When CART2HARM_PROJECTION is True, computes an SVD-based null-space projector P = Q Qᵀ (where the columns of Q span null_space(cart2harm_constraint)) and applies it to each MO in place, verifying that the residual violation drops below 1e-10.

Finally, applies cart2harm_matrix to produce valence_molorb_harm_coeff.

Raises:

AssertionError – If basis-function counts are inconsistent after the transformation.

process_core_orbitals()

Build the coefficient matrix for frozen-core MOs in the spherical basis.

Reads the core-orbital expansion coefficients (ccor) and the number of core MOs per shell type (nrcorb). Constructs core_molorb_coeff of shape (Nharmbasfns, Ncore_molorbs) by placing each core AO coefficient into the appropriate rows (atom’s harmonic basis functions) and column (core MO index).

Supports s- and p-type core shells only.

finalize_coefficients()

Concatenate core and valence MO coefficient matrices into coeff and build the combined AO normalization array norm_per_harmbasfn.

After this step, coeff[spin] has shape (Nharmbasfns, Nmolorbs[spin]) and is ready to be passed to stowfn.StoWfn.

setup_stowfn()

Construct and populate a stowfn.StoWfn object from all previously computed attributes.

Computes the nuclear repulsion energy from the Atomic Distances matrix (for multi-atom systems), populates all geometry, basis, and orbital attributes on sto, and calls stowfn.StoWfn.check_and_normalize() to validate the result.

apply_cusp_correction()

Enforce nuclear cusp conditions on all MOs and write stowfn.data.

For each spin and each MO the cusp constraint violation A c is evaluated. Any MO with a violation norm above 1e-9 is flagged and corrected according to CUSP_METHOD:

Method

Action

'enforce'

Applies the cusp-enforcing matrix M to pin the coefficient of the highest-exponent s-type AO on each centre to the cusp-satisfying value while leaving all other coefficients unchanged.

'project'

Projects the coefficient vector onto the null space of A, removing all cusp-violating components.

'none'

No correction applied.

When PLOT_CUSPS is True, MO values and Laplacians are evaluated before and after correction for use by plot_cusps(). After processing, stowfn.data is saved.

Raises:

AssertionError – If any residual cusp violation exceeds 1e-8 after correction.

plot_cusps()

Generate and save cusp_constraint.svg.

Does nothing when PLOT_CUSPS is False. Produces a 2 × Natoms grid of subplots:

  • Top row — MO values along the z-axis through each nucleus (dashed = before cusp correction, solid = after).

  • Bottom row — local energy \(E_\mathrm{loc} = -\nabla^2\psi/(2\psi) - Z/|r|\) along the same line, with y-limits set from the post-correction extrema.

Only orbitals with a non-zero cusp violation (stored in self.fixed) are plotted. Sign normalization is applied so all MOs are positive at the midpoint of the z-axis segment.

Raises:

ImportError – If matplotlib is not installed.

run()

Execute the full ADF → CASINO conversion pipeline:

  1. process_valence_basis()

  2. process_core_basis()

  3. process_shells()

  4. process_coefficients()

  5. process_core_orbitals()

  6. finalize_coefficients()

  7. setup_stowfn()

  8. apply_cusp_correction() → writes stowfn.data

  9. plot_cusps() → writes cusp_constraint.svg (optional)

adf2stowf.main()

Command-line entry point.

Parses arguments with argparse and delegates to ADFToStoWF.run().

Command-line interface

Synopsis

python adf2stowf.py [options]

Options

Flag

Default

Description

--plot-cusps

False

Save cusp_constraint.svg with wavefunction values and local energies along the z-axis through each nucleus before and after cusp correction.

--cusp-method {enforce,project,none}

enforce

Nuclear cusp correction strategy (see apply_cusp_correction()).

--cart2harm-projection

False

Project MO coefficients onto the pure spherical-harmonic subspace before the Cartesian-to-spherical transformation, removing angular-momentum contamination.

--all-orbitals

False

Include virtual (unoccupied) MOs in the output.

--dump

False

Write a human-readable TAPE21.txt dump of all parsed ADF data.

Input / output files

File

Description

TAPE21.asc (input)

ASCII export of the ADF binary checkpoint, must be in the working directory.

stowfn.data (output)

CASINO STO wavefunction file, always produced.

TAPE21.txt (optional output)

Human-readable dump of parsed data (requires --dump).

cusp_constraint.svg (optional output)

Cusp diagnostic plot (requires --plot-cusps).

Usage examples

Minimal conversion (defaults)

$ python adf2stowf.py

Debug run — dump raw data and use projection cusp method

$ python adf2stowf.py --dump --cusp-method=project

Inspect cusps for all orbitals including virtuals

$ python adf2stowf.py --plot-cusps --all-orbitals --cusp-method=enforce

Enforce spherical purity, skip cusp correction

$ python adf2stowf.py --cart2harm-projection --cusp-method=none

Programmatic use

from adf2stowf import adf2stowf

conv = adf2stowf.ADFToStoWF(
    plot_cusps=False,
    cusp_method='enforce',
    do_dump=False,
    cart2harm_projection=True,
    only_occupied=True,
)
conv.run()

# Access the resulting StoWfn object directly
sto = conv.sto
print('Number of MOs:', sto.num_molorbs)

Cartesian to spherical-harmonic transformation

ADF stores MO coefficients in a Cartesian STO basis; CASINO requires real solid-harmonic (spherical) basis functions. For shells with angular momentum ℓ ≥ 2, the number of Cartesian functions exceeds the number of spherical ones, introducing redundant (contaminating) components:

Shell

Cartesian fns

Spherical fns

Contaminating component(s)

d (ℓ=2)

6

5

1 s-type: \(x^2 + y^2 + z^2\)

f (ℓ=3)

10

7

3 p-type: \(x(x^2+y^2+z^2)\), \(y(x^2+y^2+z^2)\), \(z(x^2+y^2+z^2)\)

Basis function ordering

The polynomial ordering within each shell type follows stowfdet.f90 from CASINO. Spherical functions are labelled by magnetic quantum number m; Cartesian monomials use exponent triples (i, j, k) with \(x^i y^j z^k\).

d-shell — Cartesian order: \(xy,\ xz,\ yz,\ x^2{-}y^2,\ 2z^2{-}x^2{-}y^2,\ x^2{+}y^2{+}z^2\)

m

Real solid harmonic (unnormalized)

\(-2\)

\(xy\)

\(-1\)

\(yz\)

\(0\)

\(2z^2 - x^2 - y^2\)

\(+1\)

\(xz\)

\(+2\)

\(x^2 - y^2\)

(s-contam)

\(x^2 + y^2 + z^2\)

f-shell — Cartesian order: \(x^3,\ x^2y,\ x^2z,\ xy^2,\ xyz,\ xz^2,\ y^3,\ y^2z,\ yz^2,\ z^3\)

m

Real solid harmonic (unnormalized)

\(0\)

\(2z^3 - 3z(x^2+y^2)\)

\(+1\)

\(4xz^2 - x(x^2+y^2)\)

\(-1\)

\(4yz^2 - y(x^2+y^2)\)

\(+2\)

\((x^2-y^2)z\)

\(-2\)

\(xyz\)

\(+3\)

\(x^3 - 3xy^2\)

\(-3\)

\(3x^2y - y^3\)

(p-contam x)

\(x(x^2+y^2+z^2)\)

(p-contam y)

\(y(x^2+y^2+z^2)\)

(p-contam z)

\(z(x^2+y^2+z^2)\)

Matrix equations

For each shell, the harm2cart_map matrix M (shape \(N_\text{cart} \times N_\text{cart}\)) relates spherical-harmonic coefficients \(\mathbf{c}_\text{sph}\) to Cartesian coefficients \(\mathbf{c}_\text{cart}\):

\[\mathbf{c}_\text{cart} = \mathbf{M}\, \mathbf{c}_\text{sph}\]

The matrix M is partitioned row-wise into a physical block and a constraint block:

\[\begin{split}\mathbf{M} = \begin{pmatrix} \mathbf{M}_\text{phys} \\ \mathbf{A} \end{pmatrix}\end{split}\]

where \(\mathbf{M}_\text{phys}\) has \(N_\text{sph}\) rows (the true spherical harmonics) and \(\mathbf{A}\) has \(N_\text{cart} - N_\text{sph}\) rows (the contaminating components stored in cart2harm_constraint).

The inverse transformation (stored in cart2harm_map) is:

\[\mathbf{c}_\text{sph} = \mathbf{M}^{-1}\, \mathbf{c}_\text{cart}\]

and the physical block of \(\mathbf{M}^{-1}\) gives the global cart2harm_matrix C (shape \(N_\text{harm}^\text{total} \times N_\text{cart}^\text{total}\)), a block-diagonal matrix assembled over all centres and shells:

\[\mathbf{c}_\text{sph}^\text{MO} = \mathbf{C}\, \mathbf{c}_\text{cart}^\text{MO}\]

Contamination check. A Cartesian coefficient vector \(\mathbf{c}_\text{cart}\) is pure spherical if and only if:

\[\mathbf{A}\, \mathbf{c}_\text{cart} = \mathbf{0}\]

Violations are measured as \(\|\mathbf{A}\, \mathbf{c}_\text{cart}\|\) and logged as warnings when they exceed \(10^{-5}\).

Projection (``–cart2harm-projection``). When contamination is present, the coefficient vector is projected onto the null space of A via:

\[\mathbf{c}_\text{proj} = \mathbf{Q}\mathbf{Q}^\top \mathbf{c}_\text{cart}\]

where the columns of Q form an orthonormal basis for \(\ker(\mathbf{A})\), computed by scipy.linalg.null_space. After projection, \(\mathbf{A}\,\mathbf{c}_\text{proj} = \mathbf{0}\) to machine precision.

Nuclear cusp conditions

Each MO must satisfy the cusp condition at every nucleus c with charge Z_c:

\[\left.\frac{\partial \psi}{\partial r}\right|_{r=0} = -Z_c \, \psi(0)\]

This is encoded as a linear constraint A c = 0 on the AO coefficient vector c (see stowfn.StoWfn.cusp_constraint_matrix()). The three --cusp-method strategies differ in how they restore this condition when it is violated:

  • enforce — adjusts only the coefficient of the highest-exponent s-type AO on each centre via the matrix M = I − efix A+ A, leaving all other coefficients unchanged.

  • project — applies the null-space projector Q QT to remove all cusp-violating components simultaneously.

  • none — performs no correction; useful for diagnosis or when the original orbitals already satisfy the cusp.

Dependencies

Package

Usage

numpy

All array arithmetic throughout the pipeline.

scipy.linalg

null_space for the Cartesian-to-spherical projection (imported lazily when --cart2harm-projection is used).

matplotlib

Cusp diagnostic plots (imported lazily when --plot-cusps is used).

argparse

Command-line argument parsing in main().

adf2stowf.adfread

adfread.AdfParser for reading TAPE21.asc.

adf2stowf.stowfn

stowfn.StoWfn for constructing and writing stowfn.data.