chemistrykit.quantum

chemistrykit.quantum hub · History · Examples · API reference

chemistrykit.quantum: quantum chemistry – pedagogical but numerically real.

Exactly solvable models (particle in a 1D/3D box, applied to the Kuhn free-electron model of conjugated-dye UV-Vis absorption; the quantum harmonic oscillator compared against the exact Morse-potential vibrational levels; the rigid rotor; hydrogen-like radial wavefunctions and orbital shapes) and genuinely variational ones (Huckel molecular-orbital theory for conjugated pi systems, including Huckel’s 4n+2 aromaticity rule checked against the computed spectrum; a minimal Gaussian-basis variational treatment of H2+; the effective-nuclear-charge variational helium atom; Rayleigh-Schrodinger perturbation theory for the anharmonic oscillator, checked against exact numerical diagonalization in a truncated basis).

Every variational model here reduces to the same secular equation \(HC=SCE\) (chemistrykit.quantum.utils.secular_equation), which chemistrykit.quantum.systems.huckel, chemistrykit.quantum.systems.hartree_fock, and chemistrykit.quantum.systems.perturbation all build a matrix Hamiltonian for and hand off to it – mirroring the tight-binding Hamiltonian-diagonalization pattern physicskit’s physicskit.condensed.tight_binding uses for solid-state band structure.

class chemistrykit.quantum.EigenstateResult(energies, coefficients, basis_labels=<factory>, overlap=None, extra=<factory>)

Bases: object

Container for the output of a VariationalSolver.solve() call.

Mirrors chemistrykit.kinetics.core.base_system.KineticsResult and chemistrykit.md.core.base_system.MDResult (a stable dataclass return type consumed by visualizers/tests), specialized to a diagonalized-Hamiltonian spectrum: ordered eigenvalues plus the matrix of eigenvectors (molecular-orbital coefficients), each column one eigenstate.

Parameters:

• energies (ndarray)

• coefficients (ndarray)

• basis_labels (Sequence[str])

• overlap (ndarray | None)

• extra (dict)

basis_labels: Sequence[str]

Human-readable label for each basis function (e.g. atom labels for Huckel pi orbitals, center labels for an LCAO basis), same order as coefficients’ rows.

Type:

tuple of str

coefficients: ndarray

Eigenvector matrix; column i is the basis-function expansion coefficients of the state with energy energies[i].

Type:

ndarray, shape (n_basis, n_states)

energies: ndarray

Eigenvalues (orbital/state energies), ascending.

Type:

ndarray, shape (n_states,)

extra: dict

Free-form slot for additional diagnostics a solver chooses to attach (e.g. total electronic energy, nuclear repulsion).

Type:

dict

is_degenerate(i, j, tol=1e-09)

Return whether states i and j are degenerate to within tol.

Parameters:

• i (int)

• j (int)

• tol (float) – Absolute energy tolerance, in the same units as energies.

Return type:

bool

Returns:

bool

Examples

>>> import numpy as np
>>> result = EigenstateResult(energies=np.array([-2.0, -1.0, -1.0, 1.0]), coefficients=np.eye(4))
>>> result.is_degenerate(1, 2)
True
>>> result.is_degenerate(0, 1)
False

orbital_coefficients(i)

Return the basis-expansion coefficients of the i-th state.

Parameters:

i (int)

Return type:

ndarray

Returns:

ndarray, shape (n_basis,)

orbital_energy(i)

Return the energy of the i-th state (0-indexed, ascending).

Parameters:

i (int)

Return type:

float

Returns:

float

overlap: ndarray | None = None

The overlap matrix S used to obtain this result, or None when the basis was assumed orthonormal (S = I), as in standard Huckel theory.

Type:

ndarray, shape (n_basis, n_basis), or None

class chemistrykit.quantum.ExponentOptimizationResult(optimized_alpha, optimized_energy, naive_alpha, naive_energy, converged)

Bases: object

Result of H2PlusVariational.optimize_exponent().

Parameters:

• optimized_alpha (float)

• optimized_energy (float)

• naive_alpha (float)

• naive_energy (float)

• converged (bool)

converged: bool

Whether the underlying 1D optimizer reported success.

Type:

bool

property improvement: float

How much lower the optimized energy is than the naive guess’s, in J (>= 0 by the variational theorem).

Type:

float

naive_alpha: float

The starting-guess exponent, in m^-2.

Type:

float

naive_energy: float

Total energy at naive_alpha, in J.

Type:

float

optimized_alpha: float

The variationally optimized Gaussian exponent, in m^-2.

Type:

float

optimized_energy: float

Total energy at optimized_alpha, in J.

Type:

float

class chemistrykit.quantum.H2PlusVariational(bond_length)

Bases: VariationalSolver

H2+ (one electron, two protons separated by bond_length) in a minimal 2-Gaussian LCAO basis.

Parameters:

bond_length (float) – Proton-proton separation R, in m.

Examples

At a plausible bond length, the bonding MO is lower in energy than the antibonding one – the basic LCAO-MO picture of a covalent bond:

>>> h2plus = H2PlusVariational(bond_length=106e-12)
>>> result = h2plus.solve(alpha=1.0 / (5.29e-11 ** 2))
>>> bool(result.energies[0] < result.energies[1])
True

hamiltonian(alpha=1.0)

Build the 2x2 core Hamiltonian (kinetic + both nuclear attractions) at exponent alpha.

Parameters:

alpha (float) – Shared Gaussian orbital exponent, in m^-2.

Return type:

ndarray

Returns:

ndarray, shape (2, 2) – Energy, in J.

property nuclear_repulsion: float

The proton-proton repulsion energy, \(e^2/(4\pi\varepsilon_0R)\), in J.

Type:

float

optimize_exponent(alpha_guess=3.5734577849564575e+20)

Variationally optimize the shared Gaussian exponent to minimize the total energy.

By the variational theorem, the computed ground-state energy at any value of alpha is an upper bound on the true H2+ ground-state energy, so minimizing over alpha with scipy.optimize.minimize_scalar() produces a strictly better (lower or equal) energy than any single fixed guess – a genuine, if minimal, variational calculation (Szabo & Ostlund, Modern Quantum Chemistry, 1st ed. rev., Ch. 1.3).

Parameters:

alpha_guess (float) – Initial bracket center for the 1D minimization, in m^-2.

Return type:

ExponentOptimizationResult

Returns:

ExponentOptimizationResult

Examples

Optimizing the exponent never makes the energy worse than the naive starting guess – the variational principle in action:

>>> h2plus = H2PlusVariational(bond_length=106e-12)
>>> naive_energy = h2plus.total_energy(alpha_guess := 1.0 / (5.29e-11 ** 2))
>>> result = h2plus.optimize_exponent(alpha_guess)
>>> bool(result.optimized_energy <= naive_energy)
True

overlap_matrix(alpha=1.0)

Build the 2x2 overlap matrix at exponent alpha.

Parameters:

alpha (float)

Return type:

ndarray

Returns:

ndarray, shape (2, 2)

solve(alpha=1.0)

Solve the secular equation \(HC=SCE\) at a given (fixed) orbital exponent.

Parameters:

alpha (float) – Shared Gaussian orbital exponent, in m^-2.

Return type:

EigenstateResult

Returns:

EigenstateResult – Two states (bonding, antibonding); extra["nuclear_repulsion"] and extra["total_energy"] (electronic ground state + nuclear repulsion) are also attached.

total_energy(alpha)

Total energy (bonding-orbital electronic energy + nuclear repulsion) at exponent alpha.

Parameters:

alpha (float) – Shared Gaussian orbital exponent, in m^-2 (must be positive).

Return type:

float

Returns:

float – Energy, in J.

class chemistrykit.quantum.HeliumVariationalResult(Z, effective_charge, energy, first_order_energy, independent_electron_energy)

Bases: object

Result of optimize_helium_like_effective_charge().

Parameters:

• Z (float)

• effective_charge (float)

• energy (float)

• first_order_energy (float)

• independent_electron_energy (float)

Z: float

True nuclear charge.

Type:

float

effective_charge: float

Optimal variational charge \(\zeta^*=Z-5/16\).

Type:

float

energy: float

Optimized variational energy, in J.

Type:

float

first_order_energy: float

First-order perturbation-theory energy (\(\zeta=Z\)), in J.

Type:

float

independent_electron_energy: float

Energy with electron-electron repulsion ignored, \(-Z^2E_h\), in J.

Type:

float

property ionization_energy: float

Energy to remove one electron, \(-Z^2E_h/2-E^*\), in J (negative means predicted unbound).

Type:

float

property screening_constant: float

\(Z-\zeta^*\) (always \(5/16\)).

Type:

float

class chemistrykit.quantum.HuckelSystem(n_atoms, bonds, alpha=0.0, beta=-1.0, labels=None)

Bases: VariationalSolver

A conjugated pi system’s Huckel Hamiltonian, built from an explicit bond list.

Parameters:

• n_atoms (int) – Number of conjugated (sp2, pi-contributing) atoms.

• bonds (sequence of tuple(int, int)) – Pi bonds, as 0-indexed atom pairs (i, j).

• alpha (float) – The Coulomb integral \(\alpha\), in energy units (every atom is assumed identical – carbon 2p_z – so a single scalar suffices; heteroatom Huckel parameters are not modeled here).

• beta (float) – The (negative, by convention) resonance integral \(\beta\) for a pi bond.

• labels (sequence of str, optional) – Human-readable atom labels (defaults to C1, C2, …).

Examples

Ethene (a single pi bond, 2 atoms) has the trivial 2x2 Huckel result \(\alpha\pm\beta\):

>>> ethene = HuckelSystem(n_atoms=2, bonds=[(0, 1)])
>>> result = ethene.solve()
>>> np.round(result.energies, 6)
array([-1.,  1.])

classmethod cyclic_polyene(n_atoms, alpha=0.0, beta=-1.0)

Build a cyclic conjugated ring of n_atoms sp2 carbons (e.g. benzene, n_atoms=6).

Parameters:

• n_atoms (int) – Must be >= 3.

• alpha (float)

• beta (float)

Return type:

HuckelSystem

Returns:

HuckelSystem

delocalization_energy(n_pi_electrons)

Huckel delocalization (resonance) energy relative to isolated (localized) pi bonds.

\(E_{deloc}=E_\pi(\text{Huckel})-n_{bonds}\times2(\alpha+\beta)\), comparing the delocalized Huckel pi-energy to the energy of the same number of pi electrons confined to isolated, non-interacting double bonds (each an ethene-like 2-orbital system contributing \(2(\alpha+\beta)\) when filled with 2 electrons) – the standard measure of aromatic/conjugative stabilization (Atkins & de Paula, Physical Chemistry, 11th ed., Ch. 9.5(a)).

Parameters:

n_pi_electrons (int) – Must be even (whole pi bonds) and equal to the number of pi bonds in bonds times 2.

Return type:

float

Returns:

float – More negative means more delocalization-stabilized.

frontier_electron_density(n_pi_electrons, orbital='homo', tol=1e-06)

Fukui’s frontier electron density on each atom, \(f_r=2c_{r,\mathrm{F}}^2\).

Fukui, Yonezawa, and Shingu’s frontier-orbital reactivity index (K. Fukui, T. Yonezawa, H. Shingu, J. Chem. Phys. 20, 722 (1952)): an electrophile attacks the atom with the largest HOMO density and a nucleophile the atom with the largest LUMO density. The factor 2 makes the densities sum to 2 over all atoms (two frontier electrons).

Parameters:

• n_pi_electrons (int) – Total pi electrons (must be even, so the HOMO is doubly occupied).

• orbital (str) – Which frontier orbital to use.

• tol (float) – A frontier orbital degenerate with its neighbor within tol is ambiguous and raises ValueError.

Return type:

ndarray

Returns:

ndarray, shape (n_atoms,)

Examples

Butadiene’s HOMO is concentrated on the terminal carbons, the sites where electrophiles add:

>>> density = HuckelSystem.linear_polyene(4).frontier_electron_density(4)
>>> np.round(density, 3)
array([0.724, 0.276, 0.276, 0.724])

hamiltonian()

Build the Huckel Hamiltonian matrix: \(\alpha\) on the diagonal, \(\beta\) for bonded pairs.

Return type:

ndarray

Returns:

ndarray, shape (n_atoms, n_atoms)

classmethod linear_polyene(n_atoms, alpha=0.0, beta=-1.0)

Build an open (acyclic) conjugated chain of n_atoms sp2 carbons (e.g. butadiene, n_atoms=4).

Parameters:

• n_atoms (int)

• alpha (float)

• beta (float)

Return type:

HuckelSystem

Returns:

HuckelSystem

pi_electron_energy(n_pi_electrons)

Total pi-electron energy, filling MOs lowest-energy-first, two electrons per orbital.

\(E_\pi=\sum_kn_k\varepsilon_k\), with occupation numbers \(n_k\in\{0,1,2\}\) (Aufbau + Pauli, the same filling rule used throughout this subpackage, e.g. chemistrykit.quantum.systems.particle_in_box.conjugated_dye_absorption_wavelength()).

Parameters:

n_pi_electrons (int) – Total number of pi electrons to place, <= 2 * n_atoms.

Return type:

float

Returns:

float – Total pi-electron energy, in the same units as alpha/beta.

solve()

Diagonalize the Huckel Hamiltonian (orthonormal basis: S=I).

Return type:

EigenstateResult

Returns:

EigenstateResult – energies ascending (most bonding, i.e. most negative for beta<0, first); coefficients columns are the molecular-orbital pi-electron-density expansion coefficients.

class chemistrykit.quantum.HydrogenLikeAtom(Z=1, reduced_mass=9.1093837139e-31)

Bases: QuantumSystem

A one-electron atom/ion of nuclear charge Z (H, He+, Li2+, …).

\[E_n = -\frac{Z^2\mu e^4}{2(4\pi\varepsilon_0)^2\hbar^2n^2}, \qquad n=1,2,3,\dots\]

(Atkins & de Paula, Physical Chemistry, 11th ed., eq. 8.2.) The radial wavefunctions are

\[R_{n,l}(r) = \sqrt{\left(\frac{2Z}{na_0}\right)^3\frac{(n-l-1)!}{2n\,(n+l)!}} e^{-\rho/2}\rho^lL_{n-l-1}^{2l+1}(\rho), \qquad \rho=\frac{2Zr}{na_0}\]

with \(a_0=4\pi\varepsilon_0\hbar^2/(\mu e^2)\) the (reduced-mass) Bohr radius and \(L_{n-l-1}^{2l+1}\) an associated Laguerre polynomial (Atkins & de Paula, Physical Chemistry, 11th ed., Table 8.1; Levine, Quantum Chemistry, 7th ed., eq. 6.94), and the full wavefunction is \(\psi_{n,l,m}(r,\theta,\phi)=R_{n,l}(r)Y_l^m(\theta,\phi)\).

By default reduced_mass is the bare electron mass – the infinite-nuclear-mass approximation, giving the textbook \(-13.6\,\text{eV}\) hydrogen ground state to 3 significant figures. Passing the true two-body reduced mass \(\mu=m_em_{nuc}/(m_e+m_{nuc})\) instead removes this approximation and reproduces hydrogen’s ground-state energy to 5-6 figures (\(-13.598\,\text{eV}\)) – the small residual difference is exactly the nuclear recoil the infinite-mass approximation neglects.

Parameters:

• Z (int) – Nuclear charge (in units of the elementary charge).

• reduced_mass (float) – Reduced mass of the electron-nucleus system, in kg.

Examples

Hydrogen’s ground-state ionization energy is the textbook 13.6 eV:

>>> h_atom = HydrogenLikeAtom(Z=1)
>>> round(float(-h_atom.energy(1) / 1.602176634e-19), 1)
13.6

He+ (Z=2) is bound 4x more tightly than H at the same n (\(E\propto Z^2\)):

>>> he_plus = HydrogenLikeAtom(Z=2)
>>> round(float(he_plus.energy(1) / h_atom.energy(1)), 6)
4.0

angular_wavefunction(theta, phi, l, m)

Return the (generally complex) spherical harmonic \(Y_l^m(\theta,\phi)\).

Parameters:

• theta (float or array-like of float) – Polar angle, in radians, in [0, pi].

• phi (float or array-like of float) – Azimuthal angle, in radians, in [0, 2*pi).

• l (int) – Orbital angular momentum quantum number, >= 0.

• m (int) – Magnetic quantum number, -l <= m <= l.

Returns:

complex or ndarray of complex

property bohr_radius: float

The (reduced-mass) Bohr radius \(a_0=4\pi\varepsilon_0\hbar^2/(\mu e^2)\), in m.

Type:

float

check_radial_normalization(n, l, r_max_bohr_radii=60.0)

Numerically verify \(\int_0^\infty R_{n,l}(r)^2r^2\,dr=1\).

A direct, formula-agnostic sanity check of radial_wavefunction() (independent of any particular normalization-constant convention): integrates the radial distribution function (radial_distribution_function()) out to r_max_bohr_radii Bohr radii, where the (exponentially decaying) integrand is already negligible.

Parameters:

• n (int)

• l (int)

• r_max_bohr_radii (float) – Upper integration limit, in units of bohr_radius.

Return type:

float

Returns:

float – Should equal 1.0 to within numerical-quadrature error.

Examples

>>> h_atom = HydrogenLikeAtom(Z=1)
>>> round(h_atom.check_radial_normalization(n=1, l=0), 6)
1.0
>>> round(h_atom.check_radial_normalization(n=2, l=1), 6)
1.0

energy(n)

Return \(E_n=-Z^2\mu e^4/(2(4\pi\varepsilon_0)^2\hbar^2n^2)\).

Parameters:

n (int or array-like of int) – Principal quantum number(s), each >= 1.

Returns:

float or ndarray – Energy, in J (negative – bound).

radial_distribution_function(r, n, l)

Return the radial distribution function \(P(r)=r^2R_{n,l}(r)^2\).

The probability of finding the electron in a thin spherical shell between r and r+dr, obtained by integrating \(|\psi|^2\) over all angles (Atkins & de Paula, Physical Chemistry, 11th ed., eq. 8.4).

Parameters:

• r (float or array-like of float)

• n (int)

• l (int)

Returns:

float or ndarray – Probability density, in m^-1.

radial_wavefunction(r, n, l)

Return the radial wavefunction \(R_{n,l}(r)\).

Parameters:

• r (float or array-like of float) – Radial distance(s) from the nucleus, in m.

• n (int) – Principal quantum number, >= 1.

• l (int) – Orbital angular momentum quantum number, 0 <= l <= n - 1.

Returns:

float or ndarray – Amplitude, in m^-3/2.

Examples

The 1s radial wavefunction is a simple decaying exponential, \(R_{1,0}(r)=2a_0^{-3/2}e^{-r/a_0}\) (Atkins & de Paula, Physical Chemistry, 11th ed., Table 8.1):

>>> h_atom = HydrogenLikeAtom(Z=1)
>>> a0 = h_atom.bohr_radius
>>> R10 = h_atom.radial_wavefunction(a0, n=1, l=0)
>>> expected = 2.0 * a0 ** -1.5 * np.exp(-1.0)
>>> round(float(R10 / expected), 6)
1.0

wavefunction(r, theta, phi, n, l, m)

Return the full (generally complex) wavefunction \(\psi_{n,l,m}(r,\theta,\phi)=R_{n,l}(r)Y_l^m(\theta,\phi)\).

Parameters:

• r (float or array-like of float) – Radial distance, in m.

• theta (float or array-like of float) – Polar angle, in radians.

• phi (float or array-like of float) – Azimuthal angle, in radians.

• n (int) – Quantum numbers, 0 <= l <= n - 1, -l <= m <= l.

• l (int) – Quantum numbers, 0 <= l <= n - 1, -l <= m <= l.

• m (int) – Quantum numbers, 0 <= l <= n - 1, -l <= m <= l.

Returns:

complex or ndarray of complex – Amplitude, in m^-3/2.

class chemistrykit.quantum.MorseOscillator(mass, force_constant, dissociation_energy)

Bases: QuantumSystem

A diatomic-like vibration in the Morse potential, the standard anharmonic model of a real chemical bond.

\[V(x) = D_e\left(1-e^{-ax}\right)^2, \qquad E_v = \hbar\omega\left(v+\frac12\right) - \frac{(\hbar\omega)^2}{4D_e}\left(v+\frac12\right)^2\]

(P. M. Morse, Phys. Rev. 34, 57 (1929); Atkins & de Paula, Physical Chemistry, 11th ed., eq. 8.9, giving the exact vibrational eigenvalues of the Morse potential – an important contrast with QuantumHarmonicOscillator, whose energies are exact only for the harmonic (quadratic) approximation to the true potential, not for any real bond.) The bond dissociates once the levels stop being bound – v_max gives the last bound vibrational level, beyond which this closed form no longer applies (Atkins & de Paula, Physical Chemistry, 11th ed., Ch. 8.1).

omega is fixed by the curvature at the bottom of the well, exactly as for the harmonic oscillator, so the two are directly comparable at a shared force constant (see compare_harmonic_vs_morse()).

Parameters:

• mass (float) – Reduced mass, in kg.

• force_constant (float) – Force constant at the potential minimum, k, in N/m (fixes \(\omega=\sqrt{k/m}\), identical to the harmonic oscillator’s).

• dissociation_energy (float) – Well depth \(D_e\) (measured from the potential minimum, not from the zero-point level), in J.

Examples

Level spacing decreases monotonically with v (anharmonicity), unlike the harmonic oscillator’s perfectly even spacing:

>>> morse = MorseOscillator(mass=1.6e-27, force_constant=500.0, dissociation_energy=7.0e-19)
>>> spacing = [morse.energy(v + 1) - morse.energy(v) for v in range(4)]
>>> bool(all(spacing[i] > spacing[i + 1] for i in range(3)))
True

property angular_frequency: float

\(\omega=\sqrt{k/m}\), in rad/s – same formula as the harmonic oscillator.

Type:

float

property anharmonicity_constant: float

The dimensionless anharmonicity constant \(x_e=\hbar\omega/(4D_e)\).

Type:

float

energy(v)

Return the exact Morse vibrational energy \(E_v\).

Parameters:

v (int or array-like of int) – Vibrational quantum number(s), each >= 0 (and, physically, <= v_max, though this is not enforced).

Returns:

float or ndarray – Energy, in J, measured from the bottom of the well.

property v_max: int

The highest bound vibrational quantum number.

Found from \(dE_v/dv=0\) at \(v=v_{max}+\frac12\), i.e. the last integer v before the Morse levels turn over and start decreasing (an artifact of the quadratic-in-v formula beyond the true dissociation limit; Atkins & de Paula, Physical Chemistry, 11th ed., Ch. 8.1).

Type:

int

class chemistrykit.quantum.ParticleInBox1D(length, mass=9.1093837139e-31)

Bases: QuantumSystem

A particle of mass m confined to an infinite 1D box of length L, \(0\le x\le L\).

\[E_n = \frac{n^2h^2}{8mL^2}, \qquad \psi_n(x) = \sqrt{\frac2L}\sin\!\left(\frac{n\pi x}{L}\right), \qquad n=1,2,3,\dots\]

(Atkins & de Paula, Physical Chemistry, 11th ed., eq. 7.6a-b.)

Parameters:

• length (float) – Box length L, in m.

• mass (float) – Particle mass, in kg.

Examples

Doubling the box length cuts every energy level by a factor of 4 (\(E\propto1/L^2\)):

>>> box1 = ParticleInBox1D(length=1.0e-9)
>>> box2 = ParticleInBox1D(length=2.0e-9)
>>> round(float(box1.energy(1) / box2.energy(1)), 6)
4.0

energy(n)

Return \(E_n=n^2h^2/(8mL^2)\).

Parameters:

n (int or array-like of int) – Quantum number(s), each >= 1.

Returns:

float or ndarray – Energy, in J.

Examples

>>> box = ParticleInBox1D(length=1.0e-9)  # a 1 nm box
>>> round(float(box.energy(1) / 1.602176634e-19), 4)  # in eV
0.376

probability_density(x, n)

Return \(|\psi_n(x)|^2\).

Parameters:

• x (float or array-like of float)

• n (int)

Returns:

float or ndarray – Probability density, in m^-1.

wavefunction(x, n)

Return \(\psi_n(x)=\sqrt{2/L}\sin(n\pi x/L)\).

Parameters:

• x (float or array-like of float) – Position(s), in m; should lie in [0, length] (zero outside, per the infinite-wall boundary condition, but this is not enforced – callers restrict the domain themselves).

• n (int) – Quantum number, >= 1.

Returns:

float or ndarray – Amplitude, in m^-1/2.

Examples

The wavefunction vanishes at both walls (the boundary condition that quantizes n in the first place):

>>> box = ParticleInBox1D(length=1.0e-9)
>>> bool(abs(box.wavefunction(0.0, n=1)) < 1e-9)
True
>>> bool(abs(box.wavefunction(1.0e-9, n=1)) < 1e-9)
True

class chemistrykit.quantum.ParticleInBox3D(Lx, Ly, Lz, mass=9.1093837139e-31)

Bases: QuantumSystem

A particle confined to an infinite rectangular box of dimensions (Lx, Ly, Lz).

\[E_{n_x,n_y,n_z} = \frac{h^2}{8m}\left(\frac{n_x^2}{L_x^2}+\frac{n_y^2}{L_y^2}+\frac{n_z^2}{L_z^2}\right)\]

(Atkins & de Paula, Physical Chemistry, 11th ed., eq. 7.10.) A cubic box (\(L_x=L_y=L_z\)) shows genuine degeneracy – e.g. the \((2,1,1)\), \((1,2,1)\), \((1,1,2)\) states are degenerate – a textbook illustration of how symmetry produces degenerate levels (Atkins & de Paula, Physical Chemistry, 11th ed., Ch. 7.3(c)).

Parameters:

• Lx (float) – Box dimensions, in m.

• Ly (float) – Box dimensions, in m.

• Lz (float) – Box dimensions, in m.

• mass (float) – Particle mass, in kg.

Examples

>>> box = ParticleInBox3D(Lx=1.0e-9, Ly=1.0e-9, Lz=1.0e-9)
>>> e211 = box.energy(2, 1, 1)
>>> e121 = box.energy(1, 2, 1)
>>> e112 = box.energy(1, 1, 2)
>>> round(e211, 30) == round(e121, 30) == round(e112, 30)
True

degeneracy(n_max)

Enumerate energies (rounded) and their degeneracies up to nx, ny, nz <= n_max.

Parameters:

n_max (int) – Largest quantum number to include along each axis.

Return type:

dict

Returns:

dict – Maps a (relative-tolerance-rounded) energy value, in J, to the number of (nx, ny, nz) triples sharing it.

Examples

A cubic box has a 3-fold-degenerate first excited level:

>>> box = ParticleInBox3D(Lx=1.0e-9, Ly=1.0e-9, Lz=1.0e-9)
>>> degeneracies = box.degeneracy(n_max=2)
>>> sorted(degeneracies.values())[-2]
3

energy(nx, ny, nz)

Return \(E_{n_x,n_y,n_z}\).

Parameters:

• nx (int) – Quantum numbers, each >= 1.

• ny (int) – Quantum numbers, each >= 1.

• nz (int) – Quantum numbers, each >= 1.

Returns:

float – Energy, in J.

class chemistrykit.quantum.QuantumHarmonicOscillator(mass, force_constant)

Bases: QuantumSystem

A 1D quantum harmonic oscillator, potential \(V(x)=\frac12kx^2\).

\[\omega = \sqrt{k/m}, \qquad E_v = \hbar\omega\left(v+\frac12\right), \qquad v=0,1,2,\dots\]

(Atkins & de Paula, Physical Chemistry, 11th ed., eq. 7.16-7.17.) Energy levels are exactly evenly spaced by \(\hbar\omega\) – the defining feature that MorseOscillator below (a real bond’s anharmonic potential) violates.

Parameters:

• mass (float) – Oscillator mass (the reduced mass, for a two-body vibration), in kg.

• force_constant (float) – Force constant k, in N/m.

Examples

Consecutive levels are always separated by exactly \(\hbar\omega\):

>>> ho = QuantumHarmonicOscillator(mass=1.6e-27, force_constant=500.0)
>>> spacing1 = ho.energy(1) - ho.energy(0)
>>> spacing2 = ho.energy(5) - ho.energy(4)
>>> bool(round(float(spacing1), 30) == round(float(spacing2), 30))
True
>>> round(float(spacing1 / (HBAR * ho.angular_frequency)), 9)
1.0

property angular_frequency: float

The classical angular frequency \(\omega=\sqrt{k/m}\), in rad/s.

Type:

float

energy(v)

Return \(E_v=\hbar\omega(v+1/2)\).

Parameters:

v (int or array-like of int) – Vibrational quantum number(s), each >= 0.

Returns:

float or ndarray – Energy, in J.

wavefunction(x, v)

Return the harmonic-oscillator eigenfunction \(\psi_v(x)\).

\[\psi_v(x) = \left(\frac{m\omega}{\pi\hbar}\right)^{1/4} \frac{1}{\sqrt{2^vv!}}H_v(\xi)e^{-\xi^2/2}, \qquad \xi=x\sqrt{m\omega/\hbar}\]

where \(H_v\) is the physicists’ Hermite polynomial (Atkins & de Paula, Physical Chemistry, 11th ed., eq. 7.18; here evaluated via scipy.special.hermite()).

Parameters:

• x (float or array-like of float) – Displacement from equilibrium, in m.

• v (int) – Vibrational quantum number, >= 0.

Returns:

float or ndarray – Amplitude, in m^-1/2.

Examples

The ground state is a normalized Gaussian, maximal at x=0:

>>> ho = QuantumHarmonicOscillator(mass=1.6e-27, force_constant=500.0)
>>> psi0_at_0 = ho.wavefunction(0.0, v=0)
>>> psi0_away = ho.wavefunction(1.0e-11, v=0)
>>> bool(psi0_at_0 > psi0_away > 0)
True

property zero_point_energy: float

The zero-point energy \(E_0=\frac12\hbar\omega\), in J.

Type:

float

class chemistrykit.quantum.QuantumSystem

Bases: ABC

Common interface for an exactly solvable bound-state model with a discrete spectrum.

Concrete subclasses (ParticleInBox1D, QuantumHarmonicOscillator, RigidRotor, HydrogenLikeAtom, …) implement energy() with whatever quantum number(s) are natural for that model (a single principal quantum number n, or (n, l), etc.) – mirroring chemistrykit.kinetics.core.base_system.RateLaw’s closed-form concentration(t), there is nothing to numerically integrate or diagonalize here.

abstractmethod energy(*quantum_numbers)

Return the energy eigenvalue for the given quantum number(s).

Parameters:

*quantum_numbers – One or more quantum numbers, model-specific (e.g. a single n for a 1D box, (n, l) for a hydrogen-like atom).

Returns:

float or ndarray – Energy, in J.

class chemistrykit.quantum.RigidRotor(moment_of_inertia)

Bases: QuantumSystem

A linear rigid rotor of moment of inertia I.

\[E_J = J(J+1)\frac{\hbar^2}{2I}, \qquad g_J = 2J+1, \qquad J=0,1,2,\dots\]

(Atkins & de Paula, Physical Chemistry, 11th ed., eq. 7.20-7.21.) Each level is \((2J+1)\)-fold degenerate (one state per value of the projection quantum number \(M_J=-J,\dots,J\)), and allowed microwave (pure rotational) transitions follow the selection rule \(\Delta J=\pm1\) (Atkins & de Paula, Physical Chemistry, 11th ed., Ch. 12.2), giving evenly spaced absorption lines at \(2B(J+1)\) where \(B=\hbar^2/(2I)\) is the rotational constant.

Parameters:

moment_of_inertia (float) – Moment of inertia I, in kg m^2.

Examples

>>> rotor = RigidRotor(moment_of_inertia=1.45e-46)  # ~HCl-like
>>> rotor.degeneracy(0)
1
>>> rotor.degeneracy(2)
5

degeneracy(J)

Return the degeneracy \(g_J=2J+1\) of level J.

Parameters:

J (int)

Return type:

int

Returns:

int

energy(J)

Return \(E_J=J(J+1)\hbar^2/(2I)\).

Parameters:

J (int or array-like of int) – Rotational quantum number(s), each >= 0.

Returns:

float or ndarray – Energy, in J.

classmethod from_diatomic(mass1, mass2, bond_length)

Build a RigidRotor from a diatomic’s atomic masses and bond length.

Uses the reduced mass \(\mu=m_1m_2/(m_1+m_2)\) and \(I=\mu r^2\) (Atkins & de Paula, Physical Chemistry, 11th ed., eq. 7.19).

Parameters:

• mass1 (float) – Atomic masses, in kg.

• mass2 (float) – Atomic masses, in kg.

• bond_length (float) – Bond length r, in m.

Return type:

RigidRotor

Returns:

RigidRotor

Examples

>>> import scipy.constants as sc
>>> rotor = RigidRotor.from_diatomic(mass1=1.008 * sc.atomic_mass, mass2=34.97 * sc.atomic_mass, bond_length=127.5e-12)  # HCl
>>> bool(1.0e-47 < rotor.moment_of_inertia < 5.0e-47)
True

property rotational_constant: float

The rotational constant \(B=\hbar^2/(2I)\), in J.

Type:

float

transition_energy(J)

Return the \(J\to J+1\) absorption transition energy, \(2B(J+1)\).

Parameters:

J (int) – Initial (lower) rotational quantum number, >= 0.

Return type:

float

Returns:

float – Energy, in J.

Examples

Successive rotational-spectrum lines are evenly spaced by \(2B\) – the textbook rigid-rotor selection-rule result:

>>> rotor = RigidRotor(moment_of_inertia=1.45e-46)
>>> spacing1 = rotor.transition_energy(1) - rotor.transition_energy(0)
>>> spacing2 = rotor.transition_energy(2) - rotor.transition_energy(1)
>>> bool(abs(spacing1 - spacing2) < 1e-30)
True

class chemistrykit.quantum.VariationalSolver

Bases: ABC

Common interface for a model whose spectrum comes from diagonalizing a matrix Hamiltonian.

Concrete subclasses (HuckelSystem, H2PlusVariational) build a Hamiltonian matrix H (and, for a non-orthogonal basis, an overlap matrix S) from their physical parameters and delegate the actual eigenvalue problem to chemistrykit.quantum.utils.secular_equation.solve_secular_equation().

abstractmethod hamiltonian()

Build and return this model’s Hamiltonian matrix.

Return type:

ndarray

Returns:

ndarray, shape (n_basis, n_basis)

abstractmethod solve()

Diagonalize the Hamiltonian and return the resulting spectrum.

Return type:

EigenstateResult

Returns:

EigenstateResult

chemistrykit.quantum.anharmonic_energy_levels(n_basis, mass, omega, b, c=0.0, n_levels=3)

Exact (within the truncated basis) anharmonic energy levels, by direct numerical diagonalization.

Diagonalizes anharmonic_hamiltonian_matrix() via chemistrykit.quantum.utils.secular_equation.solve_secular_equation() (an ordinary eigenvalue problem – the harmonic-oscillator basis is orthonormal) and returns the lowest n_levels eigenvalues, the numerically “exact” benchmark against which quartic_perturbation_first_order_correction()’s perturbative estimate is checked.

Parameters:

• n_basis (int) – Truncated basis size; should be well above n_levels for the returned levels to be trustworthy (see anharmonic_hamiltonian_matrix()’s truncation caveat).

• mass (float) – Oscillator mass, in kg.

• omega (float) – Angular frequency, in rad/s.

• b (float) – Quartic perturbation strength, in J/m^4.

• c (float) – Cubic perturbation strength, in J/m^3.

• n_levels (int) – Number of lowest levels to return.

Return type:

ndarray

Returns:

ndarray, shape (n_levels,) – Ascending.

Examples

For a weak quartic perturbation, exact diagonalization matches the first-order perturbative estimate to within a small relative error:

>>> mass, omega, b = 1.6e-27, 1.0e14, 1.0e18  # a weak perturbation
>>> exact = anharmonic_energy_levels(n_basis=40, mass=mass, omega=omega, b=b, n_levels=1)
>>> E0_harmonic = 0.5 * HBAR * omega
>>> E0_perturbative = E0_harmonic + quartic_perturbation_first_order_correction(0, mass, omega, b)
>>> relative_error = abs(exact[0] - E0_perturbative) / abs(E0_perturbative)
>>> bool(relative_error < 1.0e-3)
True

chemistrykit.quantum.anharmonic_hamiltonian_matrix(n_basis, mass, omega, b, c=0.0)

Build the matrix of \(\hat H=\hat H_0+c\hat x^3+b\hat x^4\) in the truncated harmonic-oscillator basis.

\(\hat H_0\) is diagonal (\(\hbar\omega(n+1/2)\)); the perturbation terms are built from powers of position_operator_matrix(). This is the unperturbed-basis matrix-mechanics approach to the anharmonic oscillator: truncating at finite n_basis is itself an approximation (states near the top of the truncated basis are contaminated by the missing higher states), so only the lowest handful of eigenvalues (well below n_basis) should be trusted – anharmonic_energy_levels() exists specifically to hand back only those.

Parameters:

• n_basis (int) – Truncated basis size.

• mass (float) – Oscillator mass, in kg.

• omega (float) – Angular frequency, in rad/s.

• b (float) – Quartic perturbation strength, in J/m^4.

• c (float) – Cubic perturbation strength, in J/m^3.

Return type:

ndarray

Returns:

ndarray, shape (n_basis, n_basis) – Symmetric, in J.

chemistrykit.quantum.compare_harmonic_vs_morse(mass, force_constant, dissociation_energy, v_max)

Tabulate harmonic-oscillator vs. Morse-potential vibrational energies at a shared force constant.

Both models share the same curvature at the well minimum (same mass, force_constant, hence the same \(\omega\)), so they agree closely for the lowest levels (\(v\approx0\), where the Morse potential is well approximated by its harmonic term) and diverge increasingly as v grows: the Morse level spacing shrinks toward zero as \(v\to v_{max}\) (approaching dissociation, where the vibrational levels become a near-continuum), while the harmonic spacing stays exactly \(\hbar\omega\) at every v – the harmonic approximation’s defining failure mode for a real bond (Atkins & de Paula, Physical Chemistry, 11th ed., Ch. 8.1).

Parameters:

• mass (float) – Reduced mass, in kg.

• force_constant (float) – Shared force constant k, in N/m.

• dissociation_energy (float) – Morse well depth \(D_e\), in J.

• v_max (int) – Highest vibrational quantum number to tabulate (may safely exceed the Morse oscillator’s own bound-level cutoff for comparison purposes – the harmonic levels stay well-defined at any v).

Returns:

• v (ndarray, shape (v_max + 1,)) – Vibrational quantum numbers 0, 1, ..., v_max.

• harmonic_energies (ndarray, shape (v_max + 1,)) – Harmonic-oscillator energies, in J.

• morse_energies (ndarray, shape (v_max + 1,)) – Morse-potential energies, in J.

Examples

The two models agree closely at v=0 (small relative difference) and disagree much more by v=5 (anharmonicity accumulates):

>>> v, E_harmonic, E_morse = compare_harmonic_vs_morse(mass=1.6e-27, force_constant=500.0, dissociation_energy=7.0e-19, v_max=5)
>>> relative_diff_0 = abs(E_harmonic[0] - E_morse[0]) / E_harmonic[0]
>>> relative_diff_5 = abs(E_harmonic[5] - E_morse[5]) / E_harmonic[5]
>>> bool(relative_diff_5 > relative_diff_0)
True

chemistrykit.quantum.conjugated_dye_absorption_wavelength(box_length, n_pi_electrons, mass=9.1093837139e-31)

Predict a linear conjugated dye’s UV-Vis absorption wavelength via Kuhn’s free-electron model.

Models the delocalized pi electrons of a linear conjugated chain (e.g. a cyanine dye) as free particles in a 1D box spanning the conjugated system (H. Kuhn, J. Chem. Phys. 17, 1198 (1949); Levine, Quantum Chemistry, 7th ed., Ch. 2.2c). With n_pi_electrons electrons filling the box levels two at a time (Pauli exclusion), the HOMO is level \(n_{HOMO}=n_{\pi}/2\) (only defined here for even n_pi_electrons, the closed-shell case) and the LUMO is \(n_{LUMO}=n_{HOMO}+1\). The HOMO->LUMO transition energy is then

\[\Delta E = E_{n_{HOMO}+1}-E_{n_{HOMO}} = \frac{h^2}{8mL^2}\left(2n_{HOMO}+1\right) = \frac{h^2}{8mL^2}(n_\pi+1)\]

and the predicted absorption wavelength follows from \(\Delta E=hc/\lambda\). This is a crude model (it ignores the actual sigma-bond framework’s effect on the box walls and treats the pi electrons as fully independent/non-interacting), but it correctly predicts the qualitative trend that longer conjugated chains (more pi electrons, longer effective box) absorb at longer wavelength – the basis of dye chemistry’s color tuning by conjugation length.

Parameters:

• box_length (float) – Effective 1D box length spanned by the conjugated pi system, in m (conventionally taken as the number of conjugated bonds times a typical C-C bond length, plus one bond length of “overhang” at each end).

• n_pi_electrons (int) – Number of pi electrons, must be a positive even integer (closed-shell filling).

• mass (float)

Return type:

float

Returns:

float – Predicted absorption wavelength, in m.

Raises:

ValueError – If n_pi_electrons is not a positive even integer.

Examples

A cyanine-dye-like chain with 8 pi electrons over a ~1.2 nm box absorbs in the visible, as observed for real cyanine dyes:

>>> wavelength = conjugated_dye_absorption_wavelength(box_length=1.2e-9, n_pi_electrons=8)
>>> bool(3.8e-7 < wavelength < 8.0e-7)  # visible range
True

Adding more conjugation (more pi electrons over a longer box) red-shifts the absorption:

>>> short = conjugated_dye_absorption_wavelength(box_length=0.8e-9, n_pi_electrons=6)
>>> long = conjugated_dye_absorption_wavelength(box_length=1.4e-9, n_pi_electrons=10)
>>> bool(long > short)
True

chemistrykit.quantum.cubic_perturbation_first_order_correction(n)

First-order energy correction from a cubic perturbation \(\hat H'=c\hat x^3\): always exactly zero.

\(\langle n|\hat x^3|n\rangle=0\) for every harmonic-oscillator eigenstate n, by parity: \(\hat x^3\) is an odd function of \(\hat x\), while \(|n\rangle\) has definite parity \((-1)^n\), so the integrand \(\psi_n\hat x^3\psi_n\) is always odd and integrates to zero over the symmetric domain \((-\infty,\infty)\) (Griffiths, Introduction to Quantum Mechanics, 2nd ed., Ch. 6.3). A real cubic term (e.g. the leading anharmonic correction to a Morse-like bond potential) therefore only shifts the energy levels at second order in perturbation theory – not implemented here, but this is why quartic_perturbation_first_order_correction() is the one that matters at first order.

Parameters:

n (int) – Vibrational quantum number, >= 0 (unused – the result is zero for every n; kept as a parameter for a uniform calling convention with quartic_perturbation_first_order_correction()).

Return type:

float

Returns:

float – Exactly 0.0.

Examples

>>> cubic_perturbation_first_order_correction(5)
0.0

chemistrykit.quantum.cyclic_polyene_eigenvalues(n_atoms, alpha=0.0, beta=-1.0)

Closed-form Huckel eigenvalues of a cyclic conjugated ring.

\[E_k = \alpha+2\beta\cos\!\left(\frac{2\pi k}{n}\right), \qquad k=0,1,\dots,n-1\]

(the “Frost circle” mnemonic; A. A. Frost & B. Musulin, J. Chem. Phys. 21, 572 (1953); Levine, Quantum Chemistry, 7th ed., eq. 16.62.) Provided as an independent closed-form cross-check of HuckelSystem.solve()’s numerical diagonalization for HuckelSystem.cyclic_polyene().

Parameters:

• n_atoms (int)

• alpha (float)

• beta (float)

Return type:

ndarray

Returns:

ndarray, shape (n_atoms,) – Ascending.

Examples

Benzene’s famous Huckel pattern – \(\alpha+2\beta\), \(\alpha+\beta\) (doubly degenerate), \(\alpha-\beta\) (doubly degenerate), \(\alpha-2\beta\):

>>> np.round(cyclic_polyene_eigenvalues(6), 4)
array([-2., -1., -1.,  1.,  1.,  2.])

chemistrykit.quantum.helium_like_variational_energy(zeta, Z=2.0)

Energy of the screened-1s-product trial wavefunction for a two-electron atom or ion.

\(E(\zeta)=(\zeta^2-2Z\zeta+5\zeta/8)\,E_h\).

Parameters:

• zeta (float or array-like of float) – Variational effective nuclear charge (dimensionless, positive).

• Z (float) – True nuclear charge (2 for helium, 1 for H-, 3 for Li+, …).

Returns:

float or ndarray – Energy, in J.

Examples

Setting \(\zeta=Z\) recovers first-order perturbation theory, \(E=(-Z^2+5Z/8)E_h=-2.75\,E_h\) for helium:

>>> round(float(helium_like_variational_energy(2.0, Z=2.0) / HARTREE_ENERGY), 6)
-2.75

chemistrykit.quantum.is_aromatic_by_huckel_rule(n_pi_electrons, energies, tol=1e-06)

Check Huckel’s \(4n+2\) rule against an actual computed Huckel spectrum.

A cyclic, fully conjugated, planar system is predicted aromatic when (Atkins & de Paula, Physical Chemistry, 11th ed., Ch. 9.5(b); E. Huckel, Z. Phys. 70, 204 (1931)):

1. it has \(4n+2\) pi electrons for some non-negative integer n, and

2. that electron count exactly fills a set of Huckel MOs with no unpaired electron in a partially-filled degenerate level (the closed-shell condition) – checked here directly against the degeneracies of the computed energies, not assumed.

Condition 2 is what makes this a genuine check against the diagonalized spectrum rather than a bare electron-counting rule: a system could satisfy \(N=4n+2\) yet still land mid-degenerate-level for an unusual (non-uniform-ring) Huckel Hamiltonian, which this function would correctly flag as not closed-shell.

Parameters:

• n_pi_electrons (int) – Number of pi electrons.

• energies (array-like of float) – The full computed Huckel spectrum (e.g. from HuckelSystem.solve()’s .energies), any order.

• tol (float) – Absolute energy tolerance for treating two levels as degenerate.

Return type:

bool

Returns:

bool

Examples

Benzene: 6 pi electrons, \(4(1)+2\), and the computed spectrum’s filling is closed-shell (both the doubly degenerate HOMO orbitals are fully occupied):

>>> benzene_energies = cyclic_polyene_eigenvalues(6)
>>> is_aromatic_by_huckel_rule(6, benzene_energies)
True

Cyclobutadiene: 4 pi electrons is not \(4n+2\) for any integer n (it is \(4n\), the antiaromatic count), so this correctly returns False regardless of the degeneracy structure:

>>> cbd_energies = cyclic_polyene_eigenvalues(4)
>>> is_aromatic_by_huckel_rule(4, cbd_energies)
False

chemistrykit.quantum.linear_polyene_eigenvalues(n_atoms, alpha=0.0, beta=-1.0)

Closed-form Huckel eigenvalues of a linear (acyclic) conjugated chain.

\[E_k = \alpha+2\beta\cos\!\left(\frac{k\pi}{n+1}\right), \qquad k=1,\dots,n\]

(Coulson’s formula; C. A. Coulson, Proc. R. Soc. Lond. A 169, 413 (1939); Levine, Quantum Chemistry, 7th ed., eq. 16.36.) Provided as an independent closed-form cross-check of HuckelSystem.solve()’s numerical diagonalization for HuckelSystem.linear_polyene().

Parameters:

• n_atoms (int)

• alpha (float)

• beta (float)

Return type:

ndarray

Returns:

ndarray, shape (n_atoms,) – Ascending.

Examples

Butadiene’s Huckel eigenvalues (\(n=4\)):

>>> np.round(linear_polyene_eigenvalues(4), 4)
array([-1.618, -0.618,  0.618,  1.618])

chemistrykit.quantum.optimize_helium_like_effective_charge(Z=2.0)

Minimize helium_like_variational_energy() over the effective charge (closed form).

\(dE/d\zeta=0\) gives \(\zeta^*=Z-5/16\) and \(E^*=-(Z-5/16)^2E_h\).

Parameters:

Z (float) – True nuclear charge; must exceed 5/16.

Return type:

HeliumVariationalResult

Returns:

HeliumVariationalResult

Examples

Helium: \(\zeta^*=27/16\) and \(E^*=-(27/16)^2E_h\approx-2.8477\,E_h\):

>>> result = optimize_helium_like_effective_charge(2)
>>> result.effective_charge
1.6875
>>> round(result.energy / HARTREE_ENERGY, 5)
-2.84766

chemistrykit.quantum.position_operator_matrix(n_basis, mass, omega)

Matrix of the position operator \(\hat x\) in the truncated harmonic-oscillator number basis.

\[\hat x = \sqrt{\frac{\hbar}{2m\omega}}\left(\hat a+\hat a^\dagger\right)\]

with \(\langle n-1|\hat a|n\rangle=\langle n|\hat a^\dagger|n-1\rangle=\sqrt n\) the standard ladder-operator matrix elements (Levine, Quantum Chemistry, 7th ed., eq. 2.44), giving a tridiagonal matrix in the \(\{|0\rangle,\dots,|n_{basis}-1\rangle\}\) basis. Powers of this matrix (via ordinary matrix multiplication) give the matrix elements of \(x^2\), \(x^3\), \(x^4\), etc. – exact within the truncated basis, and increasingly accurate for the lowest states as n_basis grows.

Parameters:

• n_basis (int) – Truncated basis size (>= 2).

• mass (float) – Oscillator mass, in kg.

• omega (float) – Angular frequency, in rad/s.

Return type:

ndarray

Returns:

ndarray, shape (n_basis, n_basis) – Symmetric tridiagonal, in m.

Examples

>>> X = position_operator_matrix(3, mass=1.6e-27, omega=1.0e14)
>>> bool(np.allclose(X, X.T))
True
>>> bool(np.allclose(np.diag(X), 0.0))  # <n|x|n> = 0 for every n (parity)
True

chemistrykit.quantum.quartic_perturbation_first_order_correction(n, mass, omega, b)

First-order energy correction from a quartic perturbation \(\hat H'=b\hat x^4\).

\[E_n^{(1)} = b\langle n|\hat x^4|n\rangle = \frac{3b\hbar^2}{4m^2\omega^2}\left(2n^2+2n+1\right)\]

using \(\langle n|\hat x^4|n\rangle=3\left(\frac{\hbar}{2m\omega}\right)^2(2n^2+2n+1)\) (Griffiths, Introduction to Quantum Mechanics, 2nd ed., Ch. 6 problem on the quartic perturbation; at \(n=0\) this reduces to the Gaussian ground state’s fourth moment \(\langle x^4\rangle=3\langle x^2\rangle^2\)). Unlike the cubic term, this survives at first order because \(\hat x^4\) is even, matching \(|n\rangle\)’s squared (always-even) probability density.

Parameters:

• n (int) – Vibrational quantum number, >= 0.

• mass (float) – Oscillator mass, in kg.

• omega (float) – Angular frequency, in rad/s.

• b (float) – Quartic perturbation strength, in J/m^4.

Return type:

float

Returns:

float – Energy correction, in J.

Examples

The correction grows with n (higher states are more delocalized, probing the quartic term’s steep walls more):

>>> correction_0 = quartic_perturbation_first_order_correction(0, mass=1.6e-27, omega=1.0e14, b=1.0e18)
>>> correction_5 = quartic_perturbation_first_order_correction(5, mass=1.6e-27, omega=1.0e14, b=1.0e18)
>>> bool(correction_5 > correction_0 > 0)
True

chemistrykit.quantum.solve_secular_equation(H, S=None)

Solve the secular equation \(HC=SCE\) for a Hermitian H and (optional) overlap S.

Parameters:

• H (array-like, shape (n_basis, n_basis)) – Hamiltonian matrix in the chosen basis. Must be Hermitian (only the values are checked for approximate symmetry, not enforced).

• S (array-like, shape (n_basis, n_basis), optional) – Overlap matrix. If omitted (or None), the basis is assumed orthonormal (\(S=I\)) and the ordinary eigenvalue problem is solved via numpy.linalg.eigh() – the approximation standard Huckel theory makes (neglect of differential overlap between atomic pi orbitals). If given, the true generalized problem is solved via scipy.linalg.eigh(), which requires S to be symmetric positive-definite; a non-positive-definite S means the basis is linearly dependent (e.g. two basis functions placed on top of each other) and raises ValueError.

Returns:

• energies (ndarray, shape (n_basis,)) – Eigenvalues, ascending.

• coefficients (ndarray, shape (n_basis, n_basis)) – Eigenvectors as columns, S-orthonormal (\(C^TSC=I\)) when S is given, orthonormal (\(C^TC=I\)) otherwise.

Raises:

ValueError – If S is given but is not symmetric positive-definite (a linearly dependent basis).

Examples

A trivial 2x2 orthonormal-basis case (S omitted) recovers the plain eigenvalues of a symmetric matrix:

>>> import numpy as np
>>> H = np.array([[0.0, -1.0], [-1.0, 0.0]])
>>> energies, C = solve_secular_equation(H)
>>> np.round(energies, 6)
array([-1.,  1.])

A non-orthogonal basis (S != I) genuinely changes the eigenvalues relative to the naive (ordinary) eigenvalue problem:

>>> S = np.array([[1.0, 0.5], [0.5, 1.0]])
>>> energies_S, _ = solve_secular_equation(H, S)
>>> bool(np.any(np.abs(np.sort(energies_S) - np.sort(energies)) > 1e-9))
True