physikmdb

Photoemission momentum maps and molecular orbitals from PhysikMDB, in Python.

The same C physics kernels the website runs, plus a client for the database behind it. Look up a calculation, ask for an orbital by name, get a numpy array.

pip install physikmdb

numpy is the only dependency, and no compiler is needed — a prebuilt kernel library ships for Linux, macOS and Windows, and is built from source only on a platform without one.

Five lines

import physikmdb

db = physikmdb.Database()
calc = db.calculation(21)

image = calc.momentum_map("HOMO")   # (200, 200); every argument but the
                                    # orbital has a default (hnu=30 eV, k_max=3 Å⁻¹, ...)

basis.bin and the HOMO's coefficients are downloaded on first use and cached, so the obvious loop costs one small request per orbital:

for name in ("HOMO", "HOMO-1", "HOMO-2"):
    image = calc.momentum_map(name)

Finding a calculation

db.systems()                                   # every molecule
db.systems(short_name="2A")
db.systems(xc_functional="B3LYP", basis_set="cc-pVTZ", charge=0)

db.calculations()                              # every calculation
db.calculations(system_id=7)                   # one molecule's
db.calculations(code="Orca", xc_functional="B3LYP")

db.filters()                                   # what those filters can be set to

calc = db.calculation(21)                      # by database id — the number in
                                               # an entry page's URL

Both listings return plain records and do no downloading.

One calculation

calc.formula, calc.code, calc.xc_functional, calc.basis_set
calc.charge, calc.spin, calc.spin_restricted
calc.total_energy, calc.homo_energy, calc.lumo_energy, calc.gap    # eV

calc.orbitals                        # every orbital
calc.orbital("HOMO-2")               # by name, case-insensitive
calc.orbital(21)                     # by row index
calc.homo, calc.lumo

calc.basis                           # the basis set
calc.coefficients                    # the full (Nmo, Nbasis) matrix

An orbital is a record — index, name, energy (eV) and energy_hartree, occupation, symmetry, spin. For an unrestricted calculation a bare "HOMO" is the up channel, the one the website shows first; pass calc.orbital("HOMO", spin="down") for the other.

Computing

Every compute method takes an orbital name, an index, an Orbital, or a list of them. Every other argument has a default:

calc.momentum_map("HOMO", hnu=30.0,           # photon energy [eV]
                  k_max=3.0, points=200,
                  angles=(0, 30, 0),          # orientation (phi, theta, psi) [deg]
                  substrate="fcc110",         # average over its domains
                  polarisation=(45, 0),       # (polar, azimuth) of A [deg]
                  polarisation_type="linear", # or "circular", "toroid"
                  s_share=0.0,                # s-polarised share [%]
                  handedness="left",          # or "right", "cd"
                  gamma=None,                 # IMFP damping [Å⁻¹]; None computes
                                              # it from the kinetic energy, as the
                                              # website itself does
                  normalise=False)

calc.momentum_map(["HOMO", "HOMO-1"], weights=[1.0, 0.5])   # incoherent sum

calc.momentum_map("HOMO", kinetic_energy=15.5)   # E_kin directly, instead of hν

calc.wavefunction("HOMO", extent=8.0, points=64)       # signed ψ(r)
calc.density(["HOMO", "HOMO-1"], extent=8.0)           # Σ|ψ(r)|²
calc.momentum_density("HOMO", k_max=3.0)               # |ψ̃(k)|²

kinetic_energy, if given, is used directly instead of deriving it from hnu. E_kin = hnu + orbital.energy, with the binding energy negative, as on the entry page. Several orbitals are summed incoherently, each at its own kinetic energy; one the photon cannot emit contributes nothing.

gamma (the inelastic-mean-free-path damping) is computed automatically from each orbital's kinetic energy unless you override it — the same "universal curve" the website's own JS uses. handedness="cd" depends on it, so it is zero everywhere only if you explicitly pass gamma=0.

Energy spectrum

A broadened density of states, like the entry page's energy plot:

energies, intensity = calc.energy_spectrum(fwhm=0.15, shape="gaussian")   # both eV

shape is "gaussian" or "lorentzian"; energy_range=(low, high) restricts the window, otherwise it is sized around the orbital energies automatically.

Excited states

A TD-DFT (casida) calculation carries its excited states, and each is a coherent sum of one-electron transitions. calc.calculation_type is "groundstate" or "casida"; find casida calculations directly rather than checking .excitations for a non-empty list as a side effect:

casida_calcs = db.calculations(calculation_type="casida")
state = calc.excitation(5)             # the 5th root, as the code numbered it
state.energy, state.oscillator_strength, state.tda
state.holes()                          # every occupied orbital it empties

Photoemission from an exciton is one map per photohole, each at its own kinetic energy hν + ε_j + Ω — an entangled state genuinely looks different at each of them, which is the point of measuring it:

for hole in state.holes():
    image = calc.exciton_momentum_map(state, hole, hnu=35.0)

Underneath, the coherent sum over conduction orbitals is a linear combination of coefficient rows, so it is a coefficient row — one field evaluation per photohole, not one per transition:

row = calc.dyson_row(state, hole=32)   # same shape as any coefficients.bin row

cutoff= is the amplitude floor |X_vc| a transition has to clear (0.05 by default, matching the website). Omit hole= to sum every photohole of the state, which is what an analyser with no energy resolution would see.

Working offline

calc.save("naphthalene/")              # basis.bin + coefficients.bin + calculation.json
calc = physikmdb.load("naphthalene/")  # same object, no network at all

load() also opens a folder holding just basis.bin and coefficients.bin — the website's download button, or your own writer. Without calculation.json there are no orbital names or energies, so address orbitals by index and pass kinetic_energy=.

Units

eV and Ångström, the same units the website's own controls are labelled in. Energies in eV, k_max and gamma in Å⁻¹, extent in Å.

Orbitals carry both: orbital.energy is eV, orbital.energy_hartree is Hartree. physikmdb.units holds the two constants and the four conversions, and is the only place in the package where a number changes meaning.

Plotting

Optional, and deliberately small — enough to see whether a map looks right:

pip install physikmdb[plot]

from physikmdb import plot

ax = plot.momentum_map(image, k_max=3.0, title="HOMO")
ax.figure.savefig("homo.png")   # it's a plain matplotlib Axes - use it as usual

Each function takes and returns an ordinary Axes (ax= to draw into an existing one), so once matplotlib is installed you drive it directly - import matplotlib.pyplot as plt for multi-panel figures, ax.figure for anything else. Nothing outside physikmdb.plot imports matplotlib.

The kernels, unconverted

For your own basis and coefficients, or when you want nothing at all between you and the C:

from physikmdb import kernels        # Hartree, Bohr, Bohr⁻¹ throughout

basis = kernels.read_basis("basis.bin")
rows  = kernels.read_coefficients("coefficients.bin", basis)
image = kernels.momentum_map(basis, rows[21], E_kin=0.779, k_max=1.59)

kernels.Basis is nine plain numpy arrays, so a basis you built yourself works the same way. physikmdb.binary reads and writes the basis.bin / coefficients.bin format the website serves.

Examples

examples/ in the source distribution, simplest first:

Licence

EUPL-1.2. See LICENSE.

The scientific data served by a PhysikMDB instance is licensed separately — see physikmdb.uni-graz.at/license.