more docs

docs/hamiltonian.md

@@ -77,53 +77,76 @@ assert np.allclose(np.sum(forces, axis=0), 0.0)     # forces are conservative --
 ```
 A unique feature of psiflow `Hamiltonian` instances is the ability to create a new
 hamiltonian from a linear combination of two or more existing hamiltonians.
+This is relevant for many types of free energy calculations and/or enhanced sampling
+techniques, including umbrella sampling, Hamiltonian replica exchange, or thermodynamic
+integration.
 Let us consider the particular example of [umbrella
-sampling](https://wires.onlinelibrary.wiley.com/doi/10.1002/wcms.66) using
-[PLUMED](https://www.plumed.org/):
+sampling](https://wires.onlinelibrary.wiley.com/doi/10.1002/wcms.66).
+As activated event, we consider the decay of vinyl alcohol to acetaldehyde,
+which consists of a proton jump from the oxygen to the opposite carbon:
 
-```py
-from psiflow.hamiltonians import PlumedHamiltonian
+<figure markdown="span">
+  ![Image title](hamiltonians_umbrella.svg){ width="500" }
+  <figcaption>Transformation of vinyl alcohol into acetaldehyde by means of a proton jump.
+      A reaction coordinate is constructed based on the distance of hydrogen with respect to
+      oxygen and with respect to carbon. </figcaption>
+</figure>
 
-plumed_str = ""
-bias = PlumedHamiltonian()
-
-```
+The harmonic restraint is implemented and evaluated via [PLUMED](https://www.plumed.org/).
+In psiflow, this can be done by passing a plumed input string which describes the bias
+potential into a `PlumedHamiltonian`.
 
-$$
-H = \alpha H_0 + (1 - \alpha) H_1
-$$
+```py
+from psiflow.hamiltonians import PlumedHamiltonian
 
-is very straightforward to express in code:
+plumed_str = """UNITS LENGTH=A ENERGY=kj/mol
+d_C: DISTANCE ATOMS=3,5
+d_O: DISTANCE ATOMS=1,5
+CV: COMBINE ARG=d_C,d_O COEFFICIENTS=1,-1 PERIODIC=NO
+RESTRAINT ARG=CV KAPPA=1500 AT=0.0
+"""
 
+bias = PlumedHamiltonian(plumed_str)
 
-This allows for a very natural definition of many algorithms in computational statistical physics
-(e.g. Hamiltonian replica exchange, thermodynamic integration, biased dynamics).
+```
+To add this contribution to our MACE potential, we simply sum both hamiltonians:
 
 ```py
-from psiflow.hamiltonians import PlumedHamiltonian
+potential = mace + bias
 
+# double check
+alcohol = Geometry.load('vinyl_alcohol.xyz')
+total_energy = potential.compute(alcohol, 'energy')
+mace_energy  = mace.compute(alcohol, 'energy')
+bias_energy  = bias.compute(alcohol, 'energy')
 
-plumed_input = """
+assert np.allclose(
+    total_energy.result(),
+    mace_energy.result() + bias_energy.result(),
+)
+```
 
-"""
-bias = 
+Aside from bias potentials, the combination of multiple hamiltonians is also employed in
+e.g. the calculation of anharmonic free energy corrections.
+In that case, we consider a "base" potential energy surface which is described by a
+general quadratic function (i.e. a 3Nx3N hessian matrix and a minimum-energy geometry)
+and a small perturbation which describes the difference between the quadratic
+function and the fully anharmonic potential.
+The following code snippet demonstrates the construction of mixtures of the two energy
+surfaces:
+```py
+# hessian computed via geometry optimization and finite differences
+# see sampling section
+type(hessian)  # np.ndarray
+hessian.shape  # (3n, 3n)
+type(minimum)  # Geometry
+len(minimum)   # n
 
+harmonic = Harmonic(minimum, hessian)   # create quadratic hessian; x.T @ H @ x / 2
+delta = mace - harmonic
 
-data = Dataset.load('train.xyz')
-mix = 0.5 * einstein + 0.5 * mace             # MixtureHamiltonian with E = 0.5 * E_einstein + 0.5 * E_mace
-energies_mix = mix.evaluate(data).get('energy')
+hamiltonanians = []     # linear intepolation between quadratic and MACE PES, in 10 steps
+for scale in np.linspace(0, 1, 10):
+    hamiltonians.append(hessian + scale * delta)
 
-energies_einstein = einstein.evaluate(data).get('energy')
-energies_mace     = mace.evaluate(data).get('energy')
-assert np.allclose(
-      energies_mix.result(),
-      0.5 * energies_einstein.result() + 0.5 * energies_mace.result(),
-      )
 ```
-This makes it very easy to introduce bias potentials into your simulations -- see for example the formic acid transition state [example](https://github.com/molmod/psiflow/tree/main/examples/formic_acid_transition.py).
-The following is a list of all available hamiltonians in psiflow:
-
-- `EinsteinCrystal`: A simple harmonic potential which binds atoms to a reference position.
-- `MACE`: ML potential, either pretrained as available on GitHub, or trained within psiflow (see later sections)
-- `Harmonic`: A general quadratic potential based on a Hessian matrix and an optimized geometry.
-- `PlumedHamiltonian`: a bias contribution based on a PLUMED input file.