Add some informative documentation

docs/src/dirac.md

@@ -1,6 +1,15 @@
 # Dirac Operators
 
-To create a Dirac operators, use the constructors below:
+Dirac operators are structs that hold a reference to the gauge background `U`, a temporary
+fermion field in case one wants to use the doubly flavoured Hermitian variant, the bare mass
+and the a Boolean indicating whether there are antiperiodic boundary conditions in the time
+direction (yes, only periodic and antiperiodic BCs are supported so far). \\
+To create a Dirac operators, use the constructors below. One thing to note is that dirac
+operators can be constructed using any `Abstractfield` and so the gauge background is always
+set to `nothing`. In order to then add a gauge background, must use the dirac operator as
+a functor on a `Gaugefield`, like `D_U = D_free(U)` (This does not overwrite the `U` in
+`D_free` but creates a new dirac operator, that references the same temporary fermion field
+as the parent).
 
 ```@docs
 WilsonDiracOperator

docs/src/gaugefields.md

@@ -1,6 +1,33 @@
-# Creating Gaugefields
+# Creating Fields on the lattice
 
-To create 4-dimensional SU(3) gauge field, use the constructor `U = Gaugefield(...)` (see below) and set the initial conditions with `identity_gauges!(U)` or `random_gauges!(U)`.
+To create 4-dimensional SU(3) gauge field, use the constructor `U = Gaugefield(...)`
+(see below) and set the initial conditions with `identity_gauges!(U)` or
+`random_gauges!(U)`. \\
+Gaugefields and the other two fields conatining group/algebra valued elements are structs
+that contain a main Array `U`, which is a 5-dimensional array of statically sized 3x3
+complex matrices, i.e., `SMatrix` objects from `StaticArrays.jl` (where arrays are stored as
+Tuples under the hood. The different backends are handles by `Kernelabstractions.jl` which
+exports a `zeros` method such that the array can be created on all supported backends with
+no extra work. \\
+The fact that the elements are statically sized immutable arrays means that, for one, there
+are no allocations when performing linear algebra operations with them and second and
+secondly, we always just override the matrices in the arrays instead of mutating them. This
+yields enormous benefits in terms of less headaches during development and lets us define
+custom linear algebra routines for SMatrices and SVectors. \\
+
+There are no plans to use more memory efficient storage schemes for SU(3) or su(3) elements
+yet. \\
+
+Fermion fields or spinors or whatever you want to call them are stored in 4-dimensional 
+arrays of `n_color * n_dirac` complex valued `SVector`s. The reason for chosing 4 instead of
+5 dimensions is that this enabled us to write routines that take care of all dirac
+components at the same time, which is more efficient (I think). \\
+When using even-odd preconditioned dirac operators, the fermion fields get wrapped in a
+struct called `EvenOdd` such that we can overload all functions on that type. Our convention
+is to define the fields on the even sites. While we haven't tested whether the following is
+actually more performant, we map all even sites to the first half of the array to have
+contiguous memory accesses. The function `eo_site` does exactly this mapping. \\
+For `Fermionfield`s we have the `ones!` and `gaussian_pseudofermions!` methods to init them.
 
 ```@docs
 Gaugefield

docs/src/index.md

@@ -2,27 +2,32 @@
 
 Inspired by the [LatticeQCD.jl](https://github.com/akio-tomiya/LatticeQCD.jl/tree/master) package by Akio Tomiya et al.
 
-## Features
+## Features:
 - Simulations of 4D-SU(3) Yang-Mills (Pure Gauge) theory
+- Simulations of full lattice QCD with arbitrary number of flavours (Staggered, Wilson-Clover)
 - [Metadynamics](https://www.researchgate.net/publication/224908601_Metadynamics_A_method_to_simulate_rare_events_and_reconstruct_the_free_energy_in_biophysics_chemistry_and_material_science)
 - [PT-MetaD](https://arxiv.org/abs/2307.04742)
 - Several update algorithms (HMC, Metropolis, Heatbath, Overrelaxation)
 - Several symplectic integrators for HMC (Leapfrog, OMF2, OMF4)
 - Gradient flow with variable integrators (Euler, RK2, RK3, RK3W7)
 - Improved Gauge actions (Symanzik tree, Iwasaki, DBW2)
 - Improved Topological charge definitions (clover, rectangle clover-improved)
+- Wilson fermions with and without clover improvement
+- Staggered fermions
+- Even-odd preconditioner
+- RHMC to simulate odd number of flavours
 - Support for CUDA and ROCm backends
 
-## Installation
-First make sure you have a Julia version at or above 1.9.0 installed. You can use [juliaup](https://github.com/JuliaLang/juliaup) for that or just install the release from the [Julia website](https://julialang.org/downloads/).
+## Installation:
+First make sure you have a Julia version 1.9.x or 1.10.x installed. You can use [juliaup](https://github.com/JuliaLang/juliaup) for that or just install the release from the [Julia website](https://julialang.org/downloads/).
 
 The package in not in the general registry. So you will have to either
-- Add the package to your Julia environment via:
+- Add the package to your Julia environment (**not recommended**) via:
 ```
 julia> ] add https://github.com/GianlucaFuwa/MetaQCD.jl
 ```
 
-or
+or (**recommended**)
 
 1. Clone this repository onto your machine.
 2. Open Julia in the directory which you cloned the repo into, with the project specific environment. This can either be done by starting Julia with the command line argument "--project" or by activating the environment within an opened Julia instance via the package manager:
@@ -44,3 +49,27 @@ pkg> instantiate
 ```
 
 If you want to use a GPU, make sure you not only have CUDA.jl or AMDGPU.jl installed, but also a fairly recent version of the CUDA Toolkit or ROCm.
+
+## Quick Start:
+1. Set parameters using one of the templates in template folder
+2. From shell, do:
+```
+julia --threads=auto run.jl parameters.toml
+```
+
+or
+
+2. Start Julia (with project):
+```
+julia --threads=auto --project
+```
+3. Import MetaQCD package:
+``` julia
+using MetaQCD
+```
+4. Begin Simulation with prepared parameter file "parameters.toml":
+``` julia
+run_sim("parameters.toml")
+```
+
+If you want to use a GPU, make sure you not only have CUDA.jl or AMDGPU.jl installed, but also a fairly recent version of the CUDA Toolkit or ROCm.

docs/src/parameters.md

@@ -1,10 +1,11 @@
 # Full Parameter list (= default):
 ```julia
-Base.@kwdef mutable struct PrintPhysicalParameters
-    L::NTuple{4, Int64} = (4, 4, 4, 4)
+Base.@kwdef mutable struct PhysicalParameters
+    # gauge parameters
+    L::NTuple{4,Int64} = (4, 4, 4, 4)
     beta::Float64 = 5.7
     NC::Int64 = 3
-    kind_of_gaction::String = "wilson"
+    gauge_action::String = "wilson"
     numtherm::Int64 = 10
     numsteps::Int64 = 100
     inital::String = "cold"
@@ -20,15 +21,35 @@ Base.@kwdef mutable struct PrintPhysicalParameters
     parity_update::Bool = false
 end
 
-Base.@kwdef mutable struct PrintBiasParameters
+Base.@kwdef mutable struct DynamicalFermionParameters
+    fermion_action::String = "none"
+    Nf::Union{Int,Vector{Int}} = 0
+    mass::Union{Float64,Vector{Float64}} = 0.0
+    wilson_r::Float64 = 1.0
+    wilson_csw::Float64 = 0.0
+    anti_periodic::Bool = true
+    cg_tol_action::Float64 = 1e-12
+    cg_tol_md::Float64 = 1e-14
+    cg_maxiters_action::Int64 = 1000
+    cg_maxiters_md::Int64 = 1000
+    rhmc_spectral_bound::NTuple{2,Float64} = (0.0, 64.0)
+    rhmc_recalc_spectral_bound::Bool = false
+    rhmc_order_action::Int64 = 15
+    rhmc_order_md::Int64 = 10
+    rhmc_prec_action::Int64 = 42
+    rhmc_prec_md::Int64 = 42
+    eo_precon::Bool = false
+end
+
+Base.@kwdef mutable struct BiasParameters
     kind_of_bias::String = "none"
     kind_of_cv::String = "clover"
     numsmears_for_cv::Int64 = 4
     rhostout_for_cv::Float64 = 0.125
-    is_static::Union{Bool, Vector{Bool}} = false
+    is_static::Union{Bool,Vector{Bool}} = false
     symmetric::Bool = false
     stride::Int64 = 1
-    cvlims::NTuple{2, Float64} = (-7, 7)
+    cvlims::NTuple{2,Float64} = (-7, 7)
     biasfactor::Float64 = Inf
     kinds_of_weights::Vector{String} = ["tiwari"]
     usebiases::Vector{String} = [""]
@@ -58,43 +79,42 @@ Base.@kwdef mutable struct PrintBiasParameters
     measure_on_all::Bool = false
 end
 
-Base.@kwdef mutable struct PrintSystemParameters
+Base.@kwdef mutable struct HMCParameters
+    hmc_trajectory::Float64 = 1
+    hmc_steps::Int64 = 10
+    hmc_friction::Float64 = 0.0
+    hmc_integrator::String = "Leapfrog"
+    hmc_numsmear_gauge::Int64 = 0
+    hmc_numsmear_fermion::Int64 = 0
+    hmc_rhostout_gauge::Float64 = 0.0
+    hmc_rhostout_fermion::Float64 = 0.0
+    hmc_logging::Bool = true
+end
+
+Base.@kwdef mutable struct GradientFlowParameters
+    flow_integrator::String = "euler"
+    flow_num::Int64 = 1
+    flow_tf::Float64 = 0.1
+    flow_steps::Int64 = 10
+    flow_measure_every::Union{Int64,Vector{Int64}} = 1
+end
+
+Base.@kwdef mutable struct SystemParameters
     backend::String = "cpu"
     float_type::String = "float64"
-    log_dir::String = ""
+    ensemble_dir::String = ""
     log_to_console::Bool = true
     verboselevel::Int64 = 1
     loadU_format::String = ""
     loadU_dir::String = ""
     loadU_fromfile::Bool = false
     loadU_filename::String = ""
-    saveU_dir::String = ""
     saveU_format::String = ""
     saveU_every::Int64 = 1
-    randomseed::Union{UInt64, Vector{UInt64}} = 0x0000000000000000
-    measurement_dir::String = ""
-    bias_dir::Union{String, Vector{String}} = ""
+    randomseed::Union{UInt64,Vector{UInt64}} = 0x0000000000000000
     overwrite::Bool = false
 end
 
-Base.@kwdef mutable struct PrintHMCParameters
-    hmc_trajectory::Float64 = 1
-    hmc_steps::Int64 = 10
-    hmc_friction::Float64 = π/2
-    hmc_integrator::String = "Leapfrog"
-    hmc_numsmear::Int64 = 0
-    hmc_rhostout::Float64 = 0.0
-end
-
-Base.@kwdef mutable struct PrintGradientFlowParameters
-    flow_integrator::String = "euler"
-    flow_num::Int64 = 1
-    flow_tf::Float64 = 0.1
-    flow_steps::Int64 = 10
-    flow_measure_every::Union{Int64, Vector{Int64}} = 1
-end
-
-Base.@kwdef mutable struct PrintMeasurementParameters
+Base.@kwdef mutable struct MeasurementParameters
     measurement_method::Vector{Dict} = Dict[]
 end
-```

docs/src/updates.md

@@ -1,16 +1,22 @@
 # Updating a Gaugefield
 
-To update a `Gaugefield` simply use the `update!` function that takes 2 arguments. 
-The first is the actual update algorithm [`update_alg`](@ref "Suppoorted Update Algorithms") and the second is the `Gaugefield` `U`.
+To update a `Gaugefield` simply use the `update!` function that takes 2 positional
+arguments and however many keyword arguments specific to the update algorithm. 
+The first positional is the actual update algorithm
+[`update_alg`](@ref "Suppoorted Update Algorithms") and the second is the `Gaugefield` `U`.
 ```julia
 U = Gaugefield(...)
 random_gauges!(U)
 
-update!(update_alg, U) 
+update!(update_alg, U; ...) 
 ```
 
 ## Supported Update Algorithms
 
+```@docs
+HMC
+```
+
 ```@docs
 Metropolis
 ```
@@ -22,7 +28,3 @@ Heatbath
 ```@docs
 Overrelaxation
 ```
-
-```@docs
-HMC
-```

docs/src/usage.md

@@ -21,4 +21,8 @@ using MetaQCD
 4. Begin Simulation with prepared parameter file "parameters.toml":
 ``` julia
 run_sim("parameters.toml")
-```
+```
+
+Logs, measurements and the lot are all written to files in the `ensembles` directory under
+the specified directory name. If no directory name is provided, one is generated based on
+time the simulation was started at.

docs/src/viz.md

@@ -1,4 +1,4 @@
-# Visualization:
+# Visualization
 We include the ability to visualize your data. For that, you have pass the the directory
 under "ensembles" that contains your measurements, creating a `MetaMeasuremnts` object
 holding all the measurements in `Dict` where the keys are symbols denoting the observable.