Update docs + cleanup (#161)

* add field docs

* add references

* update introduction

* simplify nflavors and hopping_matrix

* update model docs

* add choose_field

* better type defaults

* test pad_to_nflavors

* update DQMC docs

* update crosschecks

* update HBC crosscheck

* fix missing padding

docs/make.jl

@@ -10,6 +10,7 @@ makedocs(
             "Introduction" => "DQMC/Introduction.md",
             "Model" => "DQMC/model.md",
             "Lattices" => "DQMC/lattice.md",
+            "Fields" => "DQMC/fields.md",
             "DQMC" => "DQMC/dqmc.md",
             "Configuration Recorder" => "DQMC/recorder.md",
             "Update Scheduler" => "DQMC/scheduler.md",

docs/src/DQMC/Introduction.md

@@ -1,20 +1,56 @@
 # Determinant Quantum Monte Carlo
 
-Determinant Quantum Monte Carlo is a Quantum Monte Carlo algorithm for fermionic Hamiltonians. The general idea is to use the Hubbard-Stranovich transformation to simplify the Hamiltonian to one with only quadratic fermionic degrees of freeedom. This introduces a bosonic fields for quartic terms, which can be sampled by a Monte Carlo procedure.
+Determinant Quantum Monte Carlo is a Quantum Monte Carlo algorithm for fermionic Hamiltonians. The general idea is to use the Hubbard-Stranovich transformation to simplify any term with four fermionic operators to two. This comes at the cost of introducing a bosonic field which is sampled by a Monte Carlo procedure.
 
 The minimal working example for a DQMC simulation is the following.
 
 ```julia
 using MonteCarlo
 
-model = HubbardModelAttractive(4, 2)
-dqmc = DQMC(model, beta=1.0)
+dqmc = DQMC(HubbardModel(), beta = 1.0)
 run(dqmc)
 ```
 
-This will set up and run a DQMC simulation at inverse temperature $\beta = 1.0$ using an attractive Hubbard model with a four by four square lattice. Note that by default no measurements are taken. 
+This will set up and run a DQMC simulation at inverse temperature $\beta = 1.0$ using an attractive Hubbard model with a two by two square lattice. Of course this example is not very useful. In practice you will want to modify lots of parameters and components of the simulation. We give a brief overview of the different components here. You can also check the examples section for some more involved and realistic examples.
+
+#### Model
+
+The model is a description of the Hamiltonian. Currently MonteCarlo.jl only implements the `HubbardModel` which includes nearest neighbor hoppings, a chemical potential and the Hubbard interaction. One aim of this package is to make it simple to implement different models, though currently only Hubbard interactions are supported.
+
+#### Lattice
+
+The lattice is interchangeable part of the model. We currently support a few lattices natively (Chain, Square, Cubic, Triangular) and provide a Wrapper for [LatticePhysics.jl](https://github.com/janattig/LatticePhysics.jl)
+
+#### Field
+
+The field defines how the interaction is handled, i.e. what Monte Carlo configurations look like, how the interactive term is interpreted and how updates work. We currently implement Hirsch fields and Gauß-Hermite Quadrature fields, both with a density- and magnetic-channel version. Any of these choices works with any choice of the Hubbard model, though performance and accuracy can vary.
+
+#### Updates
+
+Beyond local updates we provide a selection of global and experimental parallel updates. Like models, this part is designed to be relatively easy to extend. Most global updates just generate a new configuration and call another function to accept or reject it.
+
+The updates themselves are collected in a scheduler. We currently provide two types of them - one which simply iterates through updates and one which adaptively favors updates with high acceptance rates.
+
+#### Configuration Recorder
+
+We provide an interface to record configurations during a MonteCarlo run so that they can later be replayed. This interface also implements functionality for compressing configurations.
+
+#### Measurements
+
+We currently include measurements of the following observables:
+- Equal time Greens function (and to a less direct extend time displaced Greens functions)
+- occupations
+- charge density correlations and susceptibilities
+- magnetization in x, y and z direction
+- spin density correlations and susceptibilities in x, y and z direction
+- pairing correlation and susceptibility for generic symmetries
+- current current susceptibilities which are required for the superfluid stiffness (this can also be calculated more directly)
+- total and partial energies (quadratic terms and quartic terms separated)
+
+#### DQMC
+
+The `DQMC` struct represents the full simulation. All the components above are collected here. It also contains a bunch paramters, such as the inverse temperature `beta`, the time step discretization `delta_tau`, the number of thermazition and measurements sweeps and more.
 
-In the following pages we will discuss the various components that go into a DQMC simulation.
 
 ## Derivation
 

docs/src/DQMC/dqmc.md

@@ -23,6 +23,7 @@ Beyond this there are a couple of keyword arguments which are more involved and
 * `thermalization_measurements = Dict{Symbol, AbstractMeasurement}()`: Same as the above, but the measurements run during the thermalization stage. Might be useful to judge convergence or the number of necessary sweeps.
 * `recorder = ConfigRecorder`: A recorder for the configurations generated during the simulation.
 * `recording_rate = measure_rate`: The rate at which configurations are recorded.
+* `field::Type = choose_field(model)`: The type of the field used by the simulation
 
 Running a simulation also comes with a bunch of options via keyword arguments - most dealing with saving the simulation. The options for `run(dqmc)` include:
 
@@ -31,5 +32,7 @@ Running a simulation also comes with a bunch of options via keyword arguments -
 * `safe_before::TimeType = now() + Year(100)`: Sets a time stamp before which the simulation will cancel itself and save.
 * `safe_every::TimePeriod = Hour(10000)`: Sets a time period interval for regular saves.
 * `grace_period::TimePeriod = Minute(5)`: Sets a buffer time period for saving. By default saves at least 5 minutes before the requested time.
-* `resumable_filename = "resumable_$(datestring).jld2`: Sets the name the savefile generated from safe_before and safe_every.
-* `overwrite = true`: Enables or disables overwriting of existing files. 
+* `resumable_filename = "resumable_$datestring.jld2`: Sets the name the savefile generated from safe_before and safe_every.
+* `overwrite = true`: Enables or disables overwriting of existing files. 
+* `min_update_rate = 0.001`: Sets a minimum update acceptance rate under which a simulation is considered "failed". This will only trigger after `1 / min_update_rate` sweeps and only if all updates fall below this rate.
+* `fail_filename = "failed_$datestring.jld2"`: Sets the filename of a failed simulation.

docs/src/DQMC/fields.md

@@ -0,0 +1,62 @@
+# Fields
+
+In MonteCarlo.jl fields describe how to handle the interactive term(s) in a model/Hamiltonian. Currently we implement two transformations, the Hirsch transformation and a Gauß-Hermite Quadrature based approximation of the Hubbard Stratonovich transformation, each with a density- and spin-channel variant. 
+
+## Hirsch Transformation
+
+The Hirsch transformation applies specifically to a Hubbard interaction $U (n_{i\uparrow} - 0.5) (n_{i\downarrow} - 0.5)`. Each number operator can take the values 0 or 1, which allows us to rewrite the term as either $0.5 U (n_{i\uparrow} + n_{i\downarrow} - 1)^2 - 0.25$ or $0.5 U (n_{i\uparrow} - n_{i\downarrow})^2 + 0.25$. These two cases are the density- and spin-channel respectively. The next step is to is to introduce a bosonic field $x$ such that
+
+```math
+e^{|U|\Delta \tau \left( n_{i\uparrow} - \frac{1}{2} \right) \left(n_{i\downarrow} - \frac{1}{2} \right)} 
+= \frac{1}{2} e^{-\frac{1}{4} \Delta \tau U} \sum_{x=\pm 1} e^{\alpha x (n_{i\uparrow} + n_{i\downarrow} - \frac{1}{2})}
+```
+
+or
+
+```math
+e^{|U|\Delta \tau \left( n_{i\uparrow} - \frac{1}{2} \right) \left(n_{i\downarrow} - \frac{1}{2} \right)} 
+= \frac{1}{2} e^{-\frac{1}{4} \Delta \tau U} \sum_{x=\pm 1} e^{\alpha x (n_{i\uparrow} - n_{i\downarrow})}.
+```
+
+For both of these cases we define a positive $U$ as the attractive case. The constant $\alpha = acosh(exp(0.5 \Delta\tau U))$ in the attractive case and $\alpha = acosh(exp(-0.5 \Delta\tau U))$ in the repulsive case. Both of these can be verified by expanding the sum over $x$ and going through the different possible values of the number operator $n$.
+
+In MonteCarlo.jl the first case is implemented as `DensityHirschField` and the second as `MagneticHirschField`. Both can be used with positive and negative $U$, though $\alpha$ may become complex. This will of course influence performance. Furthermore the density channel is symmetric between spin up and down, which is used to reduced the size of matrices involved.
+
+See [Quantum Monte Carlo Methods](https://doi.org/10.1017/CBO9780511902581)
+
+## Gauß-Hermite Quadrature
+
+In this case we use Gaussian quadrature to approximate the integral generated by the Hubbard Stratonovich transformation. This is a more general approach than the Hirsch transformation, though we currently only have an implementation for Hubbard interactions.
+
+We start of the same way as with the Hirsch transformation - by rewriting the our interaction as a squared term. Generally speaking we need to arrive at $exp(-0.5\Delta\tau (c_i^\dagger V_{ij} c_j + C)^2 + C^\prime)$. Using the Hubbard Stratonovich transformation this becomes
+
+```math
+e^{-\frac{1}{2}\Delta\tau (c_i^\dagger V_{ij} c_j + C)^2 + C^\prime} = e^{-\Delta\tau C^\prime} \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{\infty} e^{-\frac{1}{2}x^2 \pm x \sqrt{\Delta\tau} (c_i^\dagger V_{ij} c_j + C)} dx
+```
+
+One method of solving this integral is Gauß-Hermite quadrature. With four weights and notes we get
+
+```math
+\frac{1}{\sqrt{2\pi}} \int_{-\infty}^{\infty} e^{-\frac{1}{2}x^2 \pm x \sqrt{\Delta\tau} (c_i^\dagger V_{ij} c_j + C)} dx = \đrac{1}{4} \sum_{x = \pm 1, \pm 2} \gamma(x) e^{\sqrt{\Delta\tau} \eta(x) (c_i^\dagger V_{ij} c_j^\dagger + C)}
+```
+
+where $\gamma$ and $\eta$ give the weights and nodes for a given $x$. For the density-channel version $V_\uparrow = V_\downarrow = \sqrt{U}$ and $C = 0.5 \sqrt{U}$, for the spin channel $V_\uparrow = - V\downarrow = \sqrt{U}$ and $C = 0$. These two cases are implemented as `DensityGHQField` and `MagneticGHQField`. Just like with the Hirsch field these two methods can create real or complex interaction matrices depending on the sign $U$ and the Density case makes use of symmetry.
+
+See [ALF](https://git.physik.uni-wuerzburg.de/ALF/ALF) Documentation (Auxiliary Field Quantum Monte Carlo section) and https://arxiv.org/pdf/2009.04491.pdf
+
+## Performance & Accuracy
+
+As a rough estimate of relative performance we can look at the runtime of a small Hubbard model (6x6 square lattice, $\beta = 1$, 10 000 sweeps total):
+
+|     Field Type      | U = +1 | U = -1 |
+| ------------------- | ------ | ------ |
+| DensityHirschField  |  6.8s  |  23s   |
+| MagneticHirschField |  51s   |  17s   |
+| DensityGHQField     |  6.2s  |  22s   |
+| MagneticGHQField    |  46s   |  15s   |
+
+As you can see the runtimes are relatively stable for the repulsive Hubbard model (negative U). This is because we either have to include spin up and down channels or work with complex matrices. In the attractive case both of these fall together. With a density-channel field we make use of spin symmetry and work with real matrices, but with a magnetic-channel field both of these slow us down dramatically.
+
+In terms of accuracy none of these should technically have a big advantage over the other other in a Hubbard model. However this doesn't generally seem to be true in practice. In one model we had issues resolving occupations accurately with `DensityHirschField` and found that spin-channel fields do a better job of resolving most observables. The exception here are spin densities, which are resolved more cleanly by density-channel fields. Overall the `MagneticGHQField` gave the best results. 
+
+In general it is probably a good idea to simulate a smaller system with all the different field options to see which yields the best combination of performance and accuracy.