even-odd wilson progress(?)

README.md

@@ -113,125 +113,3 @@ You can also create a holder of a bias potential and plot it. MetaQCD.jl creates
 bias = MetaBias(myfile, which=:mytype)
 biaspotential(bias)
 ```
-
-## Full Parameter list (= default):
-```julia
-Base.@kwdef mutable struct PhysicalParameters
-    # gauge parameters
-    L::NTuple{4,Int64} = (4, 4, 4, 4)
-    beta::Float64 = 5.7
-    NC::Int64 = 3
-    gauge_action::String = "wilson"
-    numtherm::Int64 = 10
-    numsteps::Int64 = 100
-    inital::String = "cold"
-    update_method::Vector{String} = ["HMC"]
-    hb_maxit::Int64 = 10^5
-    numheatbath::Int64 = 4
-    metro_epsilon::Float64 = 0.1
-    metro_numhits::Int64 = 1
-    metro_target_acc::Float64 = 0.5
-    eo::Bool = true
-    or_algorithm::String = "subgroups"
-    numorelax::Int64 = 0
-    parity_update::Bool = false
-end
-
-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
-    symmetric::Bool = false
-    stride::Int64 = 1
-    cvlims::NTuple{2,Float64} = (-7, 7)
-    biasfactor::Float64 = Inf
-    kinds_of_weights::Vector{String} = ["tiwari"]
-    usebiases::Vector{String} = [""]
-    write_bias_every::Int64 = 1
-    # metadynamics specific
-    bin_width::Float64 = 1e-2
-    meta_weight::Float64 = 1e-3
-    penalty_weight::Float64 = 1000.0
-    # opes specific
-    barrier::Float64 = 0.0
-    sigma0::Float64 = 0.1
-    sigma_min::Float64 = 1e-6
-    fixed_sigma::Bool = false
-    no_Z::Bool = false
-    opes_epsilon::Float64 = 0.0
-    threshold::Float64 = 1.0
-    cutoff::Float64 = 0.0
-    # for parametric
-    bias_Q::Float64 = 0.0
-    bias_A::Float64 = 0.0
-    bias_Z::Float64 = 0.0
-    # tempering specific
-    tempering_enabled::Bool = false
-    numinstances::Int64 = 1
-    swap_every::Int64 = 1
-    non_metadynamics_updates::Int64 = 1
-    measure_on_all::Bool = false
-end
-
-Base.@kwdef mutable struct HMCParameters
-    hmc_trajectory::Float64 = 1
-    hmc_steps::Int64 = 10
-    hmc_friction::Float64 = 0.0
-    hmc_integrator::String = "Leapfrog"
-    hmc_rafriction::Float64 = 1.0
-    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
-    float_type::String = "float64"
-    ensemble_dir::String = ""
-    log_to_console::Bool = true
-    verboselevel::Int64 = 1
-    save_config_format::String = ""
-    save_config_every::Int64 = 0
-    save_checkpoint_format::String = ""
-    save_checkpoint_every::Int64 = 0
-    load_config_fromfile::Bool = false
-    load_config_format::String = ""
-    load_config_path::String = ""
-    load_checkpoint_fromfile::Bool = false
-    load_checkpoint_format::String = "jld2"
-    load_checkpoint_path::String = ""
-    randomseed::Union{UInt64,Vector{UInt64}} = 0x0000000000000000
-    overwrite::Bool = false
-end
-```

docs/src/gaugefields.md

@@ -15,7 +15,7 @@ and set the initial conditions with `identity_gauges!(U)` (cold) or
 
 `Gaugefield`s, `Colorfield`s and `Expfield`s 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.
+objects from `StaticArrays.jl` (where arrays are stored as Tuples under the hood).
 
 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 secondly that we
@@ -39,9 +39,9 @@ is to define the fields on the even sites. While we haven't tested whether the f
 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.
 
-`Fermionfield`s are created in the same way as `Gaugefield`s bar the gauge action type
-parameter. For `Fermionfield`s we have the `ones!` and `gaussian_pseudofermions!` methods
-to init them.
+`Fermionfield`s are created in the same way as `Gaugefield`s with the gauge action type
+parameter being replaced by the number of Dirac indeces. For `Fermionfield`s we have the
+`ones!` and `gaussian_pseudofermions!` methods to init them.
 
 ```@docs
 Gaugefield

docs/src/usage.md

@@ -31,7 +31,7 @@ Logs, measurements and the lot are all written to files in the `ensembles` direc
 the specified directory name. If no directory name is provided, one is generated based on
 time the simulation was started at.
 
-## Build a Bias:
+## Build a Bias
 1. Set parameters using the "parameters_build.toml" example in template folder
 2. From shell, do:
 ```
@@ -53,7 +53,7 @@ using MetaQCD
 build_bias("parameters.toml")
 ```
 
-## Visualization:
+## Visualization
 We include the ability to visualize your data. For that, you just have to pass the directory where your ensemble lives:
 ```julia
 pkg> measurements = MetaMeasurements("my_ensemble")

docs/src/viz.md

@@ -18,7 +18,7 @@ For hadron correlators there is a special function `hadroncorrelator` that plots
 mean values of all time slices (without statistical uncertainties). Just specify the hadron
 whose correlator you want to see:
 ```julia
-hadroncorrelator(measurements; logscale=true, calc_meff=false, tf=0.0)
+hadroncorrelator(measurements, :pion; logscale=true, calc_meff=false, tf=0.0)
 ```
 
 You can also create a holder of a bias potential and plot it. MetaQCD.jl creates the bias

src/MetaQCD.jl

@@ -43,7 +43,7 @@ import .DiracOperators: AbstractDiracOperator, Daggered, DdaggerD, calc_fermion_
 import .DiracOperators: StaggeredDiracOperator, StaggeredEOPreDiracOperator, even_odd
 import .DiracOperators: WilsonDiracOperator, sample_pseudofermions!
 import .DiracOperators: StaggeredFermionAction, StaggeredEOPreFermionAction
-import .DiracOperators: WilsonFermionAction
+import .DiracOperators: WilsonFermionAction, WilsonEOPreFermionAction
 import .Fields: CPU, DBW2GaugeAction, IwasakiGaugeAction, SymanzikTadGaugeAction
 import .Fields: SymanzikTreeGaugeAction, WilsonGaugeAction, Plaquette, Clover
 import .Fields: Expfield, Colorfield, Gaugefield
@@ -75,7 +75,8 @@ export staple, staple_eachsite!, wilsonloop
 export Tensorfield, calc_kinetic_energy, gaussian_TA!
 export Fermionfield, StaggeredDiracOperator, Daggered, DdaggerD, WilsonDiracOperator
 export StaggeredEOPreDiracOperator, even_odd, sample_pseudofermions!
-export StaggeredFermionAction, StaggeredEOPreFermionAction, WilsonFermionAction
+export StaggeredFermionAction, StaggeredEOPreFermionAction
+export WilsonFermionAction, WilsonEOPreFermionAction
 export calc_fermion_action, gaussian_pseudofermions!
 export measure, get_value, top_charge
 export EnergyDensityMeasurement, GaugeActionMeasurement, PlaquetteMeasurement

src/diracoperators/DiracOperators.jl

@@ -14,8 +14,8 @@ using ..Utils
 import KernelAbstractions as KA
 import ..Fields: Abstractfield, EvenOdd, Fermionfield, Gaugefield, Tensorfield, clear!
 import ..Fields: check_dims, clover_square, dims, even_odd, gaussian_pseudofermions!
-import ..Fields: Checkerboard2, Sequential, @latmap, @latsum, set_source!
-import ..Fields: num_colors, num_dirac, volume
+import ..Fields: Clover, Checkerboard2, Sequential, @latmap, @latsum, set_source!
+import ..Fields: fieldstrength_eachsite!, num_colors, num_dirac, volume
 
 abstract type AbstractDiracOperator end
 abstract type AbstractFermionAction end
@@ -201,6 +201,7 @@ end
 const FERMION_ACTION = Dict{Tuple{String,Bool},Type{<:AbstractFermionAction}}(
     # boolean determines whether even-odd or not
     ("wilson", false) => WilsonFermionAction,
+    ("wilson", true) => WilsonEOPreFermionAction,
     ("staggered", false) => StaggeredFermionAction,
     ("staggered", true) => StaggeredEOPreFermionAction,
 )

src/diracoperators/staggered.jl

@@ -60,7 +60,7 @@ struct StaggeredFermionAction{Nf,TD,CT,RI1,RI2,RT} <: AbstractFermionAction
         mass;
         anti_periodic=true,
         Nf=8,
-        rhmc_spectral_bound=(mass^2, 16.0),
+        rhmc_spectral_bound=(mass^2, 6.0),
         rhmc_order_md=10,
         rhmc_prec_md=42,
         rhmc_order_action=15,

src/diracoperators/staggered_eo.jl

@@ -65,7 +65,7 @@ struct StaggeredEOPreFermionAction{Nf,TD,CT,RI1,RI2,RT} <: AbstractFermionAction
         mass;
         anti_periodic=true,
         Nf=4,
-        rhmc_spectral_bound=(mass^2, 16.0),
+        rhmc_spectral_bound=(mass^2, 6.0),
         rhmc_order_action=15,
         rhmc_prec_action=42,
         rhmc_order_md=10,
@@ -272,45 +272,14 @@ end
 function solve_dirac!(
     ψ_eo, D::T, ϕ_eo, temp1, temp2, temp3, temp4, temp5; tol=1e-14, maxiters=1000
 ) where {T<:StaggeredEOPreDiracOperator}
-    bicg_stab!(ψ_eo, D, ϕ_eo, temp1, temp2, temp3, temp4, temp5; tol=tol, maxiters=maxiters)
+    error("Not implemented yet")
+    # TODO: CGNE
     return nothing
 end
 
 # We overload LinearAlgebra.mul! instead of Gaugefields.mul! so we dont have to import
 # The Gaugefields module into CG.jl, which also allows us to use the solvers for 
 # for arbitrary arrays, not just fermion fields and dirac operators (good for testing)
-function LinearAlgebra.mul!(
-    ψ_eo::TF, D::StaggeredEOPreDiracOperator{CPU,T,TF,TG}, ϕ_eo::TF
-) where {T,TF,TG}
-    @assert TG !== Nothing "Dirac operator has no gauge background, do `D(U)`"
-    U = D.U
-    mass = T(D.mass)
-    anti = D.anti_periodic
-
-    # ψₑ = Dₒₑϕₑ
-    mul_oe!(ψ_eo, U, ϕ_eo, anti, true, false)
-    # ψₑ = DₑₒDₒₑϕₑ
-    mul_eo!(ψ_eo, U, ψ_eo, anti, false, false)
-    axpby!(mass, ϕ_eo, -1 / mass, ψ_eo) # ψₑ = Dₑₑϕₑ - DₑₒDₒₒ⁻¹Dₒₑϕₑ
-    return nothing
-end
-
-function LinearAlgebra.mul!(
-    ψ_eo::TF, D::Daggered{StaggeredEOPreDiracOperator{CPU,T,TF,TG}}, ϕ_eo::TF
-) where {T,TF,TG}
-    @assert TG !== Nothing "Dirac operator has no gauge background, do `D(U)`"
-    U = D.parent.U
-    mass = T(D.parent.mass)
-    anti = D.parent.anti_periodic
-
-    # ψₒ = Dₒₑϕₑ
-    mul_oe!(ψ_eo, U, ϕ_eo, anti, true, true)
-    # ψₑ = DₑₒDₒₑϕₑ
-    mul_eo!(ψ_eo, U, ψ_eo, anti, false, true)
-    axpby!(mass, ϕ_eo, -1 / mass, ψ_eo) # ψₑ = Dₑₑϕₑ - Dₑₒ†Dₒₒ⁻¹Dₒₑ†ϕₑ
-    return nothing
-end
-
 function LinearAlgebra.mul!(
     ψ_eo::TF, D::DdaggerD{StaggeredEOPreDiracOperator{CPU,T,TF,TG}}, ϕ_eo::TF
 ) where {T,TF,TG}

src/diracoperators/wilson.jl

@@ -152,13 +152,14 @@ function Base.show(io::IO, ::MIME"text/plain", S::WilsonFermionAction{Nf}) where
         """
         
         |  WilsonFermionAction(
-        |    Nf = $Nf
-        |    MASS = $(S.D.mass)
-        |    CSW = $(S.D.csw)
-        |    CG TOLERANCE (ACTION) = $(S.cg_tol_action)
-        |    CG TOLERANCE (MD) = $(S.cg_tol_md)
-        |    CG MAX ITERS (ACTION) = $(S.cg_maxiters_action)
-        |    CG MAX ITERS (ACTION) = $(S.cg_maxiters_md)
+        |    Nf: $Nf
+        |    MASS: $(S.D.mass)
+        |    KAPPA: $(S.D.κ)
+        |    CSW: $(S.D.csw)
+        |    CG TOLERANCE (ACTION): $(S.cg_tol_action)
+        |    CG TOLERANCE (MD): $(S.cg_tol_md)
+        |    CG MAX ITERS (ACTION): $(S.cg_maxiters_action)
+        |    CG MAX ITERS (ACTION): $(S.cg_maxiters_md)
         |    RHMC INFO (Action): $(S.rhmc_info_action)
         |    RHMC INFO (MD): $(S.rhmc_info_md))
         """
@@ -172,13 +173,14 @@ function Base.show(io::IO, S::WilsonFermionAction{Nf}) where {Nf}
         """
         
         |  WilsonFermionAction(
-        |    Nf = $Nf
-        |    MASS = $(S.D.mass)
-        |    CSW = $(S.D.csw)
-        |    CG TOLERANCE (ACTION) = $(S.cg_tol_action)
-        |    CG TOLERANCE (MD) = $(S.cg_tol_md)
-        |    CG MAX ITERS (ACTION) = $(S.cg_maxiters_action)
-        |    CG MAX ITERS (ACTION) = $(S.cg_maxiters_md)
+        |    Nf: $Nf
+        |    MASS: $(S.D.mass)
+        |    KAPPA: $(S.D.κ)
+        |    CSW: $(S.D.csw)
+        |    CG TOLERANCE (ACTION): $(S.cg_tol_action)
+        |    CG TOLERANCE (MD): $(S.cg_tol_md)
+        |    CG MAX ITERS (ACTION): $(S.cg_maxiters_action)
+        |    CG MAX ITERS (ACTION): $(S.cg_maxiters_md)
         |    RHMC INFO (Action): $(S.rhmc_info_action)
         |    RHMC INFO (MD): $(S.rhmc_info_md))
         """