# quartic_hmc

A minimal, heavily-commented, single-file Hybrid Monte Carlo (HMC) code for
the pure quartic Hermitian one-matrix model:

```
S(Phi) = N * [ (b/2) Tr(Phi^2) + c Tr(Phi^4) ]
```

`Phi` is a single `N x N` random Hermitian matrix. This is the simplest
nontrivial matrix model: for `c = 0` it reduces to the exactly-solvable
Gaussian matrix model (Wigner semicircle eigenvalue distribution), and for
`b` negative enough and `c > 0` the eigenvalue distribution splits from one
lump ("one-cut") into two ("two-cut", a double well) — a simple, cheap
example of a large-N phase transition, and a nice thing to have students
find in the output eigenvalue histogram.

## Build & run

```bash
make                         # -> ./quartic_hmc
./quartic_hmc                # run with the defaults
./quartic_hmc N=30 c=2.0 ntraj=5000 seed=42   # override any parameter
```

Every parameter listed near the top of `quartic_hmc.cpp` can be overridden
on the command line as `key=value`. Unknown keys stop the program (catches
typos).

## Files

```
quartic_hmc.cpp   everything: RNG, parameters, Action()/Force(), the
                  leapfrog HMC update, observables, output, main()
Makefile          clang++ + Armadillo build
```

The whole program is one file, organised top to bottom in seven clearly
labelled sections (random numbers, parameters, physics, HMC update,
observables, output, `main()`), so the flow can be read start to finish. It
mirrors the physics and conventions of the bigger, general-purpose
`omnicode`, but with no generic "action language" — the model is hard-coded
directly, so it's visible in one short, readable function instead of hidden
behind a parser.

## What gets measured

Only three things are recorded, on purpose:

- **the action** `S` — the standard first sanity check of any HMC run
- **`Tr(Phi^2)/N`**
- **the eigenvalues of `Phi`**

Written to `data/`:

- `obs_N50_b-4_c1.dat` — columns: trajectory, `S`, `Tr(Phi^2)/N`
- `eig_N50_b-4_c1.dat` — columns: trajectory, `N` ascending eigenvalues
- `info_N50_b-4_c1.txt` — plain-text summary of the run parameters

(the filename tag changes automatically with `N`, `b`, `c`.)

## What to look for

- **Acceptance rate** and **`<exp(-deltaH)>`** are printed periodically;
  `<exp(-deltaH)>` should be close to 1 — that is the standard correctness
  check for any HMC code (it follows from exact time-reversal symmetry of
  the leapfrog integrator, independent of the step size).
- **The force check** runs once at startup and compares the analytic force
  to a numerical derivative of the action — useful to catch a sign or
  factor-of-2 error immediately if the potential is ever changed.
- **Eigenvalue histogram**: load `eig_*.dat`, flatten all rows after
  discarding the first few (extra safety margin beyond `ntherm`), and
  histogram. For `c = 0` this should approach the Wigner semicircle; for
  `b` negative and `c` large enough it should split into two bumps.

```python
import numpy as np
import matplotlib.pyplot as plt

data = np.loadtxt("data/eig_N50_b-4_c1.dat")
eigs = data[:, 1:].flatten()          # drop the trajectory column
plt.hist(eigs, bins=100, density=True)
plt.xlabel(r"$\lambda$")
plt.ylabel(r"$\rho(\lambda)$")
plt.show()
```