Programs you build with one cc line

pwdft built by viv export; not released yet

Plane-wave density-functional theory in the local-density approximation: the total energy, the Fermi level and the band structure of silicon, germanium, diamond, gallium arsenide, aluminum and copper.

Checked against: An independent numpy implementation to 1e-7 Ha and 1e-5 eV, and the LDA band energies of Marini, Onida and Del Sole, PRL 88, 016403 (2002).

Viv source: dft/

cc -O2 -ffp-contract=off -shared -fPIC pwdft.c -o libpwdft.so -lm

import pwdft
r = pwdft.run(material="si", ecut=12.0, nkmesh=3)
print(r["etot"], r["evb"].shape)       # hartree; path points x bands, eV

vivscf built by viv export; not released yet

Restricted and unrestricted Hartree–Fock for molecules in Gaussian basis sets: the total energy, the orbital energies and the density matrix.

Checked against: PySCF 2.5.0 on eleven cases (water and its cation, by UHF, in STO-3G to cc-pVTZ; one to eight waters in 6-31G*; benzene in cc-pVDZ): total energies within 1.6e-11 Ha (tolerance 1e-10), orbital energies within 3.0e-9 Ha (tolerance 1e-8).

Viv source: study/ports/pyscf/vivscf.viv over lib/hf

cc -O2 -ffp-contract=off -shared -fPIC vivscf.c -o libvivscf.so -lm

import vivscf
r = vivscf.run(basis="cc-pvdz")       # water by default; atoms= and coords= for another molecule
print(r["e_tot"], r["mo_energy"][:5])

vlasov coming soon

The 1D1V Vlasov–Poisson equations by the Cheng–Knorr split-step method: a Langmuir wave in a Maxwellian plasma and its Landau damping rate.

Checked against: The root of the Langmuir dispersion relation, ω = 1.41566 − 0.15336i at kλ_D = 0.5, from scipy.special.wofz and mpmath. Not yet within a stated tolerance: the damping rate is 12% off at the default grid.

Viv source: lang/viv/examples/vlasov_landau.viv

cc -O2 -ffp-contract=off -shared -fPIC vlasov.c -o libvlasov.so -lm

import vlasov
r = vlasov.run(Nx=128, Nv=512, dt=0.0125)
print(r["gamma_meas"], r["omega_meas"])

skyrmehf built by viv export; not released yet

Skyrme Hartree–Fock of oxygen-16 with the SLy5 force on a radial mesh: the binding energy per nucleon, the proton, neutron and charge radii, and the six occupied orbital energies.

Checked against: A numpy twin with the same discretization to 1e-7 on two meshes, and skyrme_rpa (Colò, Cao, Van Giai and Capelli, CPC 184, 142 (2013)): 8.0250 against 8.0261 MeV per nucleon after extrapolating the two meshes.

Viv source: lib/skyrme/skyrmehf.viv over lib/skyrme/skyrme.viv

cc -O2 -ffp-contract=off -shared -fPIC skyrmehf.c -o libskyrmehf.so -lm

import skyrmehf
r = skyrmehf.run(spacing=0.05)        # SLy5; t0= ... w0= for another force
print(r["e_per_a"], r["r_ch"])

boys(n_max, T) built by viv export; not released yet

The Boys functions F₀(T) … Fₙ(T), which every Gaussian integral with a 1/r in it needs, for orders up to 16.

Checked against: mpmath at 50 digits, through the incomplete gamma function: worst relative error 8.0e-16, tolerance 5e-15.

Viv source: lib/hf/boys.viv

import k1
F = k1.boys(8, 2.5)          # F_0 .. F_8
print(F[0], F[8])

plasma_z(zeta) built by viv export; not released yet

The plasma dispersion function Z(ζ) of Fried and Conte and its derivative Z′(ζ), anywhere in the complex plane.

Checked against: i√π · scipy.special.wofz: worst relative error 7.4e-15 for Z (tolerance 2e-14) and 5.8e-13 for Z′ (tolerance 1e-12).

Viv source: lib/zfunc.viv

import k1
Z, Zp = k1.plasma_z(2.0 - 0.2j)
print(Z, Zp)

faddeeva_w(z) built by viv export; not released yet

The Faddeeva function w(z) = e^(−z²) erfc(−iz), by Weideman’s rational expansion with 40 terms.

Checked against: scipy.special.wofz on −6 ≤ Re z ≤ 6, −1 ≤ Im z ≤ 4: worst relative error 7.4e-15, tolerance 2e-14.

Viv source: lib/zfunc.viv

import k1
w = k1.faddeeva_w(1.0 + 0.5j)
print(w)

enorm(x) built by viv export; not released yet

The Euclidean norm without overflow or underflow, as MINPACK computes it: entries from subnormal to 1e300 in one vector.

Checked against: mpmath at 50 digits, and scipy.linalg.norm where its sum of squares stays finite: worst relative error 1.7e-16, tolerance 1e-14.

Viv source: lib/minpack/enorm.viv

import k1, numpy as np
x = np.array([3e300, 4e300])
print(k1.enorm(x))           # 5e300; the plain sum of squares overflows

schur(A) built by viv export; not released yet

The real Schur decomposition A = Z T Zᵀ by Hessenberg reduction and Francis’s double-shift QR.

Checked against: Residuals, since T is not unique: max|A − Z T Zᵀ| / max|A| at most 3.2e-15 and max|ZᵀZ − I| at most 2.4e-15 (tolerance 2e-13); the eigenvalues match scipy.linalg.eigvals.

Viv source: lib/schur.viv

import k1, numpy as np
T, Z = k1.schur(np.random.default_rng(0).standard_normal((5, 5)))
print(np.diag(T))

solve_care(A, B, Q, R) built by viv export; not released yet

The stabilizing solution of the continuous algebraic Riccati equation, by Laub’s Schur method as in SLICOT’s SB02MD.

Checked against: scipy.linalg.solve_continuous_are, by the relative residual: at most 6.3e-15, tolerance 1e-13.

Viv source: lib/are.viv

import k1, numpy as np
A = np.array([[0.0, 1.0], [0.0, 0.0]]); B = np.array([[0.0], [1.0]])
X = k1.solve_care(A, B, np.eye(2), np.eye(1))

solve_lyapunov(A, Q) built by viv export; not released yet

The solution of A X + X Aᵀ = Q, by Bartels and Stewart as in SLICOT’s SB03MD.

Checked against: scipy.linalg.solve_continuous_lyapunov, by the relative residual: at most 8.9e-16, tolerance 1e-13.

Viv source: lib/lyap.viv

import k1, numpy as np
A = np.array([[-1.0, 2.0], [0.0, -3.0]])
X = k1.solve_lyapunov(A, np.eye(2))

lqr(A, B, Q, R) built by viv export; not released yet

The linear-quadratic regulator: the gain K with u = −K x, and the Riccati solution P.

Checked against: K from scipy.linalg.solve_continuous_are and R⁻¹BᵀP: worst relative difference 2.8e-13, tolerance 1e-9.

Viv source: lib/are.viv (lib/lqr.viv is the single-input front end)

import k1, numpy as np
A = np.array([[0.0, 1.0], [0.0, 0.0]]); B = np.array([[0.0], [1.0]])
K, P = k1.lqr(A, B, np.eye(2), np.eye(1))     # K = [1, sqrt(3)]