physikmdb

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

The same C physics kernels the website runs, plus a client for the database behind it. Look up a calculation, download what it stores, and compute everything the website computes — locally, in your own script — or point it at the output of your own quantum-chemistry run and get the full PhysikMDB interface in your browser with nothing 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. physikmdb[plot] adds matplotlib, physikmdb[archive] adds h5py.

Four lines

import physikmdb

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

image = calc.momentum_map("HOMO")   # defaults: 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)

Worked examples

Ten scripts ship with the package, simplest first. They are the fastest way in:

physikmdb examples             # list them
physikmdb examples --copy .    # write them here, to run and edit

Finding a calculation

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

db.system("2A")                                # one, by short name
db.system(7)                                   # or by database id

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

A System carries what the landing page's table shows, including the ranking that decides which of a molecule's calculations it links to:

system.short_name, system.full_name, system.formula
system.point_group, system.molecular_mass                # e.g. "D2h", 128.17 u
system.fundamental_gap, system.ionisation_potential      # eV, of its best
system.electron_affinity, system.best_calculation_id     # calculation
system.excitation_count                                  # 0 = no exciton page
system.papers                                            # see below
system.aliases                                           # other names, e.g. ("Naphthalin",)
system.pubchem_cid                                       # e.g. 931, or None

Those five are None for a molecule with no calculation carrying both a HOMO and a LUMO — which is a fact about the data, not about how you looked it up.

One calculation

calc.short_name, calc.full_name, calc.formula, calc.point_group
calc.code, calc.code_version, calc.xc_functional, calc.basis_set
calc.charge, calc.spin, calc.spin_restricted, calc.author, calc.orientation
calc.total_energy, calc.homo_energy, calc.lumo_energy, calc.gap    # eV
calc.calculation_type                # "groundstate" or "casida"

calc.system                          # the molecule, as above
calc.url, calc.exciton_url           # back to the pages this came from

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

calc.structure                       # geometry, in Å
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.

Geometry

calc.structure reads the atoms out of basis.bin, so it costs nothing beyond the basis download, and writes the four formats the site's Export tab offers:

calc.structure.symbols, calc.structure.positions   # Å
calc.structure.formula()                           # Hill notation

open("mol.xyz", "w").write(calc.structure.to_xyz())
open("POSCAR", "w").write(calc.structure.to_poscar())
open("mol.cif", "w").write(calc.structure.to_cif())
open("mol.pdb", "w").write(calc.structure.to_pdb())

Stored files

calc.basis / calc.coefficients cover basis.bin / coefficients.bin. What else a calculation keeps depends on the code that produced it, so read the list rather than guessing at an extension:

for file in calc.files:
    print(file.label, file.filename, file.group)

calc.download(".hdf5")     # geometry + full basis + MO coefficients
calc.download("out")       # ORCA's output; NWChem writes ".nwo", ".molden"
calc.download("2A.inp")    # by name, if you prefer

Asking for something this calculation does not keep says so, and names what it does keep.

References

The molecule's literature, as the site's Further Reading popup lists it:

for paper in calc.papers:
    print(paper.title, paper.authors)
    print(paper.reference)         # "Science 326(5954), 702 (2009)"
    print(paper.link)              # the DOI resolver, or a stated URL

open("references.bib", "w").write(calc.bibtex())

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
                  detector_resolution=0.0,    # analyser k-resolution, FWHM [Å⁻¹]
                  horizon=True,               # absent past |k| = √(2 E_kin)
                  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.

The photoemission horizon

A momentum map is NaN outside |k| = √(2 E_kin). There is no intensity to report out there, and the kernel's 0 is a value: it colours as the bottom of a scale and averages into any mean or integral taken over the map. Use np.nanmax, np.nansum and friends, or pass horizon=False for the raw kernel output.

A map summed over several orbitals reaches as far as its widest hemisphere — a pixel is absent only when no orbital could have emitted into it.

Cutting a map

from physikmdb import KPath

path = KPath.straight((-2.5, 0.0), (2.5, 0.0))   # or KPath(k1, k2, k3) to bend it

profile = calc.line_profile("HOMO", path, hnu=30.0)
profile.along, profile.kpar, profile.intensity

The axis walks the path, so the curve is continuous even where the path doubles back; each sample is labelled by its own |k|, so a straight cut through Γ reads 2.5, …, 0, …, 2.5.

Stacking one cut per orbital on an energy axis gives the band map:

band = calc.band_map(calc.orbitals[:20], path, hnu=30.0, fwhm=0.3)
band.energies, band.kpar, band.image     # image is (Ne, Nkpar)

Each row is cut at its own photoemission horizon — a row is an energy the analyser detects at, so the edge is a parabola, not one orbital's circle.

3D fields

psi = calc.wavefunction("HOMO", extent=8.0, points=64)
calc.write_cube("homo.cube", psi, extent=8.0)
calc.write_chgcar("CHGCAR", psi, extent=8.0)

A box fitted to the molecule costs a fraction of the cube containing it — heptacene at 0.2 Bohr is 79×211×31 fitted against 211³, 18× fewer voxels. It needs origin and spacing given explicitly, and the same ones must travel to the writer:

box = dict(points=(79, 211, 31), origin=(-4.1, -11.1, -1.6), spacing=0.106)
psi = calc.wavefunction("HOMO", **box)
calc.write_cube("homo.cube", psi, origin=box["origin"], spacing=box["spacing"])

Spectra

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

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

calc.energy_spectrum(shape="delta")                  # unbroadened sticks
calc.energy_spectrum(spin="up")                      # one channel
calc.energy_spectrum(unit="cm-1")                    # or "hartree", "nm"
calc.energy_spectrum(weights=[o.occupation for o in calc.orbitals])

shape is "gaussian", "lorentzian" or "delta"; energy_range=(low, high) restricts the window, otherwise it is sized around the peaks automatically. Each line is one unit tall, not one unit of area — turning the broadening down does not rescale the plot.

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.energy_nm, state.oscillator_strength, state.tda
state.dipole, state.dipole_magnitude   # transition dipole [a.u.]
state.dark                             # zero oscillator strength
state.holes()                          # every occupied orbital it empties

The optical spectrum those states produce:

energies, intensity = calc.absorption_spectrum(fwhm=0.1)
wavelengths, same = calc.absorption_spectrum(fwhm=0.1, unit="nm")

Each state is weighted by its oscillator strength, so a dark state lifts the curve nowhere — it is still a real state, with photoholes and momentum maps, that light simply cannot reach.

Which kind of exciton it is

character = state.character()
character.case                         # "i", "ii", "iii" or "iv"
character.holes, character.electrons, character.nto

The four structures of Kern, Windischbacher & Puschnig, Phys. Rev. B 108, 085132 (2023), Table I — and they decide what the maps look like: one map at one E_kin (i), different maps at different E_kin (ii), the same map at several (iii), or one map of a coherent sum (iv).

Photoemission from an exciton

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/")              # binaries + calculation.json + excitons.json
calc = physikmdb.load("naphthalene/")  # same object, no network at all

A saved folder carries the metadata, the excited states and the references, so it can still cite itself.

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:

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. 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:

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

(python -m physikmdb … works the same way if the script is not on your PATH.)

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.

Talking to an instance

db = physikmdb.Database()                          # the public one
db = physikmdb.Database("http://127.0.0.1:5001")   # a local dev server
db = physikmdb.Database(cache_dir="~/physikmdb", timeout=60, retries=4)

Downloads are cached between runs, one subfolder per instance. A public instance rate-limits per client, and a loop over calculations is exactly the shape that trips it — a limited request is waited out and retried rather than raised, so the loop is not your problem. retries=0 turns that off.

Units

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

Orbitals carry both: orbital.energy is eV, orbital.energy_hartree is Hartree. physikmdb.units holds the constants and the 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

plot.energy_spectrum(energies, intensity, homo=calc.homo_energy)
plot.absorption_spectrum(energies, intensity, unit="nm")
plot.line_profile(profile, corners=path.corner_distances())
plot.band_map(band)
plot.field_slice(psi, extent=8.0)
plot.save(ax, "figure.png")

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.

The modules

Three unrelated functions are called momentum_map, and they differ in what they take: Calculation.momentum_map works in eV (hnu=), kernels.momentum_map in Hartree (E_kin=), and plot.momentum_map draws one that either has already returned.

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.