stowfn — Slater-Type Orbital Wavefunction I/O and Evaluation

Overview

This module provides the StoWfn class for working with Slater-Type Orbital (STO) wavefunction files in the format used by CASINO. It handles:

Parsing and serialising the structured plain-text .data file format;
Evaluating atomic orbitals (AOs), molecular orbitals (MOs), their gradients and Laplacians at arbitrary positions in space;
Computing normalization factors for the STO basis;
Enforcing nuclear cusp conditions through projection and linear transformation matrices.

The module supports restricted and spin-unrestricted (alpha/beta) calculations and basis sets up to g-type shells (angular momentum quantum number ℓ = 4).

Module-level constants

stowfn.F2P_bool

dict — Mapping from Fortran logical literals to Python booleans.

{'.false.': False, '.true.': True}

stowfn.P2F_bool

dict — Inverse mapping from Python booleans to Fortran logical literals.

{False: '.false.', True: '.true.'}

stowfn.num_orbs_per_shelltype

numpy.ndarray of shape (7,) — Number of angular basis functions (polynomials) for each shell-type code.

Shell code	Type	Number of functions
0	—	0 (unused placeholder)
1	s	1
2	sp	4 (one s + three p combined)
3	p	3
4	d	5
5	f	7
6	g	9

StoWfn class

class stowfn.StoWfn(fname=None)

Container and evaluator for a Slater-Type Orbital wavefunction.

Parameters:: fname (str or None) – Path to the wavefunction data file. When supplied the file is parsed immediately. When None an empty object is created via initempty().

Instance attributes (populated by readfile())

General / metadata

title: str: Title string from the first line of the file.

code: str: Name of the quantum-chemistry code that generated the file.

periodicity: int: Dimensionality / periodicity of the system (0 = molecule, 1/2/3 = periodic).

spin_unrestricted: bool: True if the calculation is spin-unrestricted (separate alpha and beta orbitals).

nuclear_repulsion_energy: float: Nuclear repulsion energy in atomic units per atom.

num_elec: int: Total number of electrons.

Geometry

num_atom: int: Number of atoms in the system.

atompos: numpy.ndarray: Atomic positions in bohr, shape (num_atom, 3).

atomnum: numpy.ndarray: Atomic numbers (nuclear charges), shape (num_atom,).

atomcharge: numpy.ndarray: Valence charges for each atom, shape (num_atom,).

Basis set

num_centres: int: Number of STO expansion centres.

centrepos: numpy.ndarray: Cartesian coordinates of each centre in bohr, shape (num_centres, 3).

num_shells: int: Total number of STO shells across all centres.

idx_first_shell_on_centre: numpy.ndarray: Zero-based index of the first shell on each centre, plus a sentinel equal to num_shells appended at the end. Shape (num_centres + 1,).

shelltype: numpy.ndarray: Shell-type code for every shell, shape (num_shells,). See num_orbs_per_shelltype for the encoding.

order_r_in_shell: numpy.ndarray: Power N of the radial prefactor r^N in each shell, shape (num_shells,).

zeta: numpy.ndarray: Slater exponent ζ for each shell, shape (num_shells,).

num_atorbs: int: Total number of atomic-orbital basis functions (AOs).

num_molorbs: numpy.ndarray: Number of molecular orbitals per spin channel. Shape (1,) for restricted or (2,) for unrestricted calculations.

Derived attributes (computed during file read / validation)

num_shells_on_centre: numpy.ndarray: Number of shells on each centre, shape (num_centres,).

max_order_r_on_centre: numpy.ndarray: Maximum radial order N present on each centre (at least 2), shape (num_centres,).

max_order_r: int: Global maximum radial order N across all centres.

max_shell_type_on_centre: numpy.ndarray: Highest shell-type code present on each centre, shape (num_centres,). Used to gate computation of higher angular- momentum polynomials.

num_spins: int: Number of spin channels (1 for restricted, 2 for unrestricted).

Orbital coefficients

coeff: list of numpy.ndarray: Expansion coefficients in the un-normalized AO basis. coeff[spin] has shape (num_molorbs[spin], num_atorbs).

coeff_norm: list of numpy.ndarray: Expansion coefficients in the normalized AO basis, obtained by multiplying each column of coeff by the corresponding AO normalization factor from get_norm(). coeff_norm[spin] has shape (num_molorbs[spin], num_atorbs).

footer: list of str: Raw lines after the orbital-coefficient block, preserved verbatim on write-back.

Methods

initempty(): Create an empty StoWfn object without reading any file. Currently a no-op placeholder; subclasses may override it.

readfile(fname)

Parse a wavefunction data file and populate all instance attributes.

Parameters:: fname (str) – Path to the file to read.
Raises:: AssertionError – If the file structure is inconsistent (e.g. shell counts do not match basis function counts).

The file is expected to follow a strict section-based format with labelled headers (BASIC INFO, GEOMETRY, BASIS SET, MULTIDETERMINANT INFORMATION, ORBITAL COEFFICIENTS). Floating- point data blocks use fixed-width 20-character fields.

writefile(fname)

Write all wavefunction data back to a file in the canonical format understood by readfile().

Parameters:: fname (str) – Destination file path.

check_and_normalize()

Validate all attributes for type and shape consistency, recompute the derived attributes (num_shells_on_centre, max_order_r_on_centre, etc.) and refresh coeff_norm.

Raises:: AssertionError – If any attribute has an unexpected type, dtype, or shape.

Call this method after modifying attributes programmatically to ensure the object remains self-consistent before evaluation or file output.

read_molorbmods(fname='correlation.data')

Read molecular orbital modifications from a CASINO-style correlation.data file. The START MOLORBMODS / END MOLORBMODS block is located and its START MOLECULAR ORBITAL COEFFICIENTS and START STO EXPONENT ZETAS sub-blocks are consumed (currently parsed but not applied; values are discarded).

Parameters:: fname (str) – Path to the correlation data file.

get_norm()

Compute the normalization factor for every atomic orbital in the basis.

The normalization constant for orbital i with angular polynomial p, radial order N, and exponent ζ is:

\[\mathcal{N}_i = Y_p \cdot (2\zeta)^{n} \sqrt{\frac{2\zeta}{(2n)!}}, \quad n = \ell_p + N + 1\]

where \(Y_p\) is the angular normalization factor for polynomial p (tabulated internally for s through g shells) and \(\ell_p\) is the polynomial power (0 for s, 1 for p, 2 for d, …).

Returns:: Normalization factors, one per AO.
Return type:: numpy.ndarray, shape (num_atorbs,)

iter_atorbs()

Iterate over all atomic orbitals in canonical order.

Yields:

(atorb, centre, nshell, N, pl) tuples where

atorb (int) — global AO index;
centre (int) — centre index;
nshell (int) — shell index within shelltype;
N (int) — radial order order_r_in_shell[nshell];
pl (int) — polynomial index within the shell.

eval_atorbs(pos)

Evaluate all atomic basis functions at a set of points.

The STO basis function for shell s at centre c evaluated at position r is:

\[\phi_{s,p}(\mathbf{r}) = r^N \, P_p(\hat{\mathbf{r}}) \, e^{-\zeta |\mathbf{r} - \mathbf{R}_c|}\]

where \(P_p\) is the real solid-harmonic polynomial for component p, N is the radial order, and ζ is the Slater exponent.

Parameters:: pos (numpy.ndarray, shape (3, num_points)) – Evaluation points in Cartesian coordinates (bohr).
Returns:: AO values at each point.
Return type:: numpy.ndarray, shape (num_points, num_atorbs)

Note

Shells whose contribution is below the numerical cutoff exp(-ζr) < exp(-746) are skipped to avoid floating-point underflow.

eval_molorbs(pos, spin=0)

Evaluate molecular orbitals at a set of points.

MOs are linear combinations of normalized AOs:

\[\psi_j(\mathbf{r}) = \sum_i C_{ji}^{\text{norm}} \, \phi_i(\mathbf{r})\]

Parameters:

pos (numpy.ndarray, shape (3, num_points)) – Evaluation points in Cartesian coordinates (bohr).
spin (int) – Spin channel index. 0 = alpha (default); 1 = beta (only valid for unrestricted calculations).

Returns:

MO values at each point.

Return type:

numpy.ndarray, shape (num_points, num_molorbs[spin])

eval_molorb_derivs(pos, spin=0)

Evaluate molecular orbitals together with their first derivatives (gradient) and second derivatives (Laplacian) at a set of points.

Parameters:

pos (numpy.ndarray, shape (3, num_points)) – Evaluation points in Cartesian coordinates (bohr).
spin (int) – Spin channel index (see eval_molorbs()).

Returns:

A three-element tuple (val, grad, lap):

val — orbital values, shape (num_points, num_molorbs[spin]);
grad — orbital gradients ∂ψ/∂x, ∂ψ/∂y, ∂ψ/∂z, shape (3, num_points, num_molorbs[spin]);
lap — orbital Laplacians ∇²ψ, shape (num_points, num_molorbs[spin]).

Return type:

tuple(numpy.ndarray, numpy.ndarray, numpy.ndarray)

Raises:

NotImplementedError – If the basis contains g-type shells (ℓ = 4); gradients and Laplacians are not yet implemented for those angular functions.

cusp_constraint_matrix()

Build the linear constraint matrix A that encodes the nuclear cusp conditions.

For each nucleus c with charge Z_c and Slater exponent ζ, the cusp condition requires

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

The resulting matrix has one row per nucleus and one column per AO. For an s-type STO centred on nucleus c:

If N = 0: coefficient is \(\mathcal{N}(Z_c - \zeta)\);
If N = 1: coefficient is \(\mathcal{N}\).

Contributions from AOs on other centres are evaluated at the nuclear position.

Returns:: Cusp constraint matrix.
Return type:: numpy.ndarray, shape (num_centres, num_atorbs)

cusp_projection_matrix()

Construct the orthogonal projector Q onto the null space of the cusp constraint matrix A.

Any orbital-coefficient vector in the image of Q satisfies A x = 0, i.e., it imposes the cusp conditions exactly.

\[\mathbf{Q} = \mathbf{V}_{\text{null}} \mathbf{V}_{\text{null}}^{\top}\]

where the columns of \(\mathbf{V}_{\text{null}}\) form an orthonormal basis for null_space(A).

Returns:: Cusp projection matrix.
Return type:: numpy.ndarray, shape (num_atorbs, num_atorbs)

cusp_fixed_atorbs()

Identify the AO index that is “fixed” (pinned) by the cusp constraint on each centre.

For each centre the s-type shell with the largest Slater exponent is selected as the representative AO; its coefficient is determined by the cusp condition rather than optimized freely.

Returns:: Array of fixed AO indices, one per centre.
Return type:: numpy.ndarray, shape (num_centres,), dtype int

cusp_enforcing_matrix()

Build the linear transformation M that maps any AO coefficient vector to one satisfying the cusp conditions.

The matrix is constructed so that the “fixed” AO coefficients (see cusp_fixed_atorbs()) are set by a least-squares inversion of the cusp constraint, leaving all other AOs unchanged:

\[\mathbf{M} = \mathbf{I} - \mathbf{e}_{\text{fix}} \, \mathbf{A}^{+} \mathbf{A}\]

where \(\mathbf{A}^{+}\) is the pseudo-inverse restricted to the fixed AO columns.

Returns:: Cusp-enforcing transformation matrix.
Return type:: numpy.ndarray, shape (num_atorbs, num_atorbs)

Usage examples

Reading and writing a wavefunction file

from stowfn import StoWfn

# Read wavefunction
wfn = StoWfn('stowfn.data')

# Inspect basis
print(f'Number of atoms:   {wfn.num_atom}')
print(f'Number of AOs:     {wfn.num_atorbs}')
print(f'Number of MOs:     {wfn.num_molorbs}')
print(f'Spin unrestricted: {wfn.spin_unrestricted}')

# Round-trip to a new file
wfn.writefile('stowfn_copy.data')

Evaluating MOs and their derivatives on a grid

import numpy as np
from stowfn import StoWfn

wfn = StoWfn('stowfn.data')

# 100 random points in a 10 bohr cube
rng = np.random.default_rng(42)
pos = rng.uniform(-5, 5, size=(3, 100))

# Orbital values (alpha spin)
mo_vals = wfn.eval_molorbs(pos, spin=0)          # (100, n_mo)

# Values + gradient + Laplacian
val, grad, lap = wfn.eval_molorb_derivs(pos, spin=0)

Applying cusp conditions

from stowfn import StoWfn

wfn = StoWfn('stowfn.data')

# Projection matrix — project coefficients into the cusp-satisfying subspace
Q = wfn.cusp_projection_matrix()          # (num_atorbs, num_atorbs)
cusp_coeff = wfn.coeff[0] @ Q.T

# Or use the enforcing matrix to pin fixed AOs
M = wfn.cusp_enforcing_matrix()           # (num_atorbs, num_atorbs)
enforced_coeff = wfn.coeff[0] @ M.T

File format

The wavefunction file is an ASCII text file divided into labelled sections. Each section begins with a header line and a separator, followed by labelled fields, one per line, with the value on the next line. Floating-point arrays use 20-character fixed-width columns (four per line); integer arrays use 10- character columns (eight per line).

<title>

BASIC INFO
----------
Generated by:
 <code_name>
Periodicity:
 <int>
Spin unrestricted:
 <.true.|.false.>
Nuclear repulsion energy (au/atom):
 <float>
Number of electrons
 <int>

GEOMETRY
--------
Number of atoms
 <int>
Atomic positions (au)
 <3·num_atom floats, 3 per line>
...

BASIS SET
---------
Number of STO centres
 <int>
...

MULTIDETERMINANT INFORMATION
----------------------------
GS

ORBITAL COEFFICIENTS (normalized AO)
-------------------------------------
 <num_molorbs[0] · num_atorbs floats>
[<num_molorbs[1] · num_atorbs floats>  (unrestricted only)]

<footer lines>

Dependencies

Package	Usage
`numpy`	Array arithmetic; used throughout for all numerical data.
`scipy.linalg`	`scipy.linalg.null_space()` for `StoWfn.cusp_projection_matrix()`.
`math`	`math.factorial()` and `math.sqrt()` for normalization.