Metadata-Version: 2.4
Name: physikmdb
Version: 1.8.3
Summary: Photoemission momentum maps and molecular orbitals from the PhysikMDB database.
Author-email: Dominik Brandstetter <dominik.brandstetter@uni-graz.at>
Maintainer-email: Dominik Brandstetter <dominik.brandstetter@uni-graz.at>
License-Expression: GPL-3.0-only
Project-URL: Homepage, https://physikmdb.uni-graz.at
Project-URL: Source, https://physikgitlab.uni-graz.at/brandstetter_dominik/physikmdb
Keywords: photoemission,momentum microscopy,ARPES,molecular orbitals,quantum chemistry,orbital tomography
Classifier: Development Status :: 5 - Production/Stable
Classifier: Intended Audience :: Science/Research
Classifier: Operating System :: MacOS
Classifier: Operating System :: Microsoft :: Windows
Classifier: Operating System :: POSIX :: Linux
Classifier: Programming Language :: Python :: 3
Classifier: Topic :: Scientific/Engineering :: Chemistry
Classifier: Topic :: Scientific/Engineering :: Physics
Requires-Python: >=3.10
Description-Content-Type: text/markdown
License-File: LICENSE
License-File: THIRD-PARTY-LICENSES.md
Requires-Dist: numpy>=1.24
Requires-Dist: jinja2>=3.1
Requires-Dist: certifi>=2024.2.2
Provides-Extra: plot
Requires-Dist: matplotlib>=3.7; extra == "plot"
Provides-Extra: archive
Requires-Dist: h5py>=3.10; extra == "archive"
Dynamic: license-file

# physikmdb

Photoemission momentum maps and molecular orbitals in the style of
[PhysikMDB](https://physikmdb.uni-graz.at), 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.

```bash
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

```python
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:

```python
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:

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

| | |
|---|---|
| `01_first_map.py` | one orbital, one map |
| `02_browse.py` | systems, calculations, filters, orbitals |
| `03_orbital_series.py` | a map per orbital, and plotting |
| `04_experiment.py` | tilt, substrate, polarisation, dichroism, damping |
| `05_offline_and_fields.py` | save/load, and the 3D fields |
| `06_energy_spectrum.py` | the broadened density of states |
| `07_kernels_directly.py` | the atomic-units layer, and your own basis |
| `08_your_own_calculation.py` | parse an ORCA/NWChem run of your own |
| `09_papers_and_spectrum.py` | references, transition dipoles, absorption |
| `10_band_map_and_profile.py` | k-path cuts and the band map |

## Finding a calculation

```python
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:

```python
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
```

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

```python
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:

```python
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:

```python
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:

```python
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:

```python
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

```python
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:

```python
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

```python
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:

```python
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:

```python
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:

```python
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:

```python
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

```python
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:

```python
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:

```python
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

```python
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:

```sh
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:

```sh
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

```python
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:

```bash
pip install physikmdb[plot]
```

```python
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:

```python
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

| | |
|---|---|
| `physikmdb.kernels` | the kernels as numpy functions, in atomic units |
| `physikmdb.maps` | detector blur, the horizon, k-path cuts, band maps |
| `physikmdb.spectra` | broaden a set of lines into a curve |
| `physikmdb.fields` | write a 3D field as `.cube` or CHGCAR |
| `physikmdb.binary` | read/write `basis.bin` and `coefficients.bin` |
| `physikmdb.units` | eV/Ångström ↔ Hartree/Bohr |
| `physikmdb.parsers` | read an ORCA/NWChem run into a calculation folder |
| `physikmdb.viewer` | serve one locally for a browser |
| `physikmdb.plot` | quick matplotlib views (needs `physikmdb[plot]`) |

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](https://physikmdb.uni-graz.at/license).
Data produced by ORCA remains subject to ORCA's own EULA, whoever holds it.
