physikmdb

Photoemission momentum maps and molecular orbitals in the style of PhysikMDB, in Python.

The same C physics kernels the website runs, plus a API wrapper for the database behind it. Look up a calculations, download data and calculate what the website does locally in your own script — or point it at the output of your own quantum chemistry calculation and get the full PhysikMDB interface in your local browser without anything leaving your machine.

pip install physikmdb

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. The dependencies are numpy and, for the browser viewer, jinja2. Optional dependecies are matplotlib and hdf5.

Four lines

import physikmdb

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

image = calc.momentum_map("HOMO") # numpy array of momenutm map with default settings

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.

calc.basis/calc.coefficients cover basis.bin/coefficients.bin. Any other stored file — the archival HDF5, the .xyz geometry, or the code's own input/output — downloads (and caches) with calc.download():

calc.download(".hdf5")   # geometry + full basis + MO coefficients
calc.download(".xyz")    # geometry
calc.download(".out")    # ORCA's output; NWChem uses ".nwo", ".molden"
calc.download(".inp")    # ORCA's input; NWChem uses ".nwi"

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

coverage= is the target fraction of the state's weight its kept pairs must cover (0.95 by default, matching the website; the stored file itself covers 0.99, so that is the ceiling). 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=.

Your own calculations

Point the package at a finished ORCA or NWChem run and look at it in your own browser, with the same pages the website uses:

python3 -m physikmdb view ./my-orca-run

It parses the folder, serves it on 127.0.0.1, and opens your browser. Nothing leaves your machine, and no dev server, Node or extra browser is involved. (physikmdb view … does the same, if the script is on your PATH.) It also opens a folder you already have — one calc.save() wrote, or one the website's download button gave you.

To write that folder without opening anything:

python3 -m physikmdb parse ./my-orca-run -o ./naphthalene

ORCA is licensed separately, and using data it produced stays subject to ORCA's own EULA. This package ships none of ORCA's code and sends none of your data anywhere.

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

GPL-3.0-only. See LICENSE, and THIRD-PARTY-LICENSES.md for the JavaScript this package bundles for view.

The website itself is EUPL-1.2; this package is GPL-3.0 because its ORCA .cis reader may derive from TheoDORE (GPL-3.0, © Felix Plasser). EUPL-1.2 Article 5 permits the combination on those terms. See THIRD-PARTY-LICENSES.md.

The scientific data served by a PhysikMDB instance is licensed separately — see physikmdb.uni-graz.at/license. Data produced by ORCA remains subject to ORCA's own EULA, whoever holds it.