JuliaQCD: Portable lattice QCD package in Julia language
Yuki Nagai Akio Tomiya





nagai.yuki@mail.u-tokyo.ac.jp
Information Technology Center, The University of Tokyo, 6-2-3 Kashiwanoha, Kashiwa, Chiba 277-0882, Japan


akio@yukawa.kyoto-u.ac.jp

Department of Information and Mathematical Sciences, Tokyo Woman’s Christian University, Tokyo 167-8585, Japan


RIKEN Center for Computational Science, Kobe 650-0047, Japan




§ ABSTRACT  
We develop a new lattice gauge theory code set  JuliaQCD using the Julia language. Julia is well-suited for integrating machine learning techniques and enables rapid prototyping and execution of algorithms for four dimensional QCD and other non-Abelian gauge theories. The code leverages LLVM for high-performance execution and supports MPI for parallel computations. Julia's multiple dispatch provides a flexible and intuitive framework for development.
The code implements existing algorithms such as Hybrid Monte Carlo (HMC), many color and flavor, supports lattice fermions, smearing techniques, and full QCD simulations. 
It is designed to run efficiently across various platforms, from laptops to supercomputers, allowing for seamless scalability. The code set is currently available on GitHub <https://github.com/JuliaQCD>.


JuliaQCD: Portable lattice QCD package in Julia language
    Akio Tomiya
    September 9, 2024
========================================================





§ INTRODUCTION

Quantum Chromo-dynamics (QCD) is the fundamental theory describing the interactions of quarks and gluons, the elementary particles that constitute the subatomic world. 
Since they can be created and annihilated from the vacuum through relativistic effects, they must be handled as a quantum many-body system with uncertain particle number. 
QCD is a quantum field theory and an accurate description of our universe. 

However, perturbation theory fails to apply in QCD due to the strong coupling constant. Lattice QCD addresses this challenge by employing lattice regularization, introducing a discrete spacetime cutoff to regulate the path integral. This approach allows for the computation of quantum expectation values of physical observables.

Lattice QCD was first formulated by K. Wilson in 1974 <cit.>, and the numerical study of QCD using this method was initiated by M. Creutz <cit.>. This framework enables the detailed study of the QCD vacuum structure and provides precise insights into the Standard Model. Additionally, it facilitates investigations into phenomena such as the quark-gluon plasma and candidates for dark matter.

As we mentioned above,
to formulate lattice QCD, a spacetime lattice cutoff is introduced. This cutoff is not merely an approximation for numerical calculations; it fundamentally defines QCD itself. Traditionally, quantum field theory is formulated on a continuous spacetime, but in the presence of interactions, such theories suffer from divergences. Although these divergences can be renormalized, they are inherent in the continuum formulation.

In QCD, we aim to calculate the quantum expectation value of a gauge-invariant operator O [For certain purposes, we sometimes calculate the expectation value of a gauge variant operator. In such cases, gauge fixing is required due to Elizur's theorem <cit.>.]. In the path integral formalism, this is expressed as

    ⟨ O ⟩
    =
    1/Z∫𝒟U 𝒟q̅𝒟q e^-S_g[U] - S_f[U, q, q̅] O(U, q, q̅),

where q and q̅ are the quark and anti-quark fields, respectively, and U is the gauge field on the lattice. The term S_g[U] represents the gauge field action, and S_f[U, q, q̅] denotes the action for the quarks. The factor Z is a normalization constant ensuring that ⟨ 1 ⟩ = 1. 
One can obtain observables in QCD evaluating this integral. 

The standard algorithm used for simulating lattice QCD is the Hybrid Monte Carlo (HMC) algorithm <cit.>. HMC has become the de facto standard because of its efficiency and effectiveness in dealing with the large-dimensional integrals for lattice gauge fields with dynamical fermions. It has been successfully used in a wide range of studies, from calculating hadron spectrum <cit.> and form factors to determining the critical temperature (T_c) of QCD. The success of lattice QCD in these areas highlights the power of numerical methods in exploring non-perturbative aspects of quantum field theory.

Recent developments in machine learning (ML) have opened up new avenues for enhancing lattice QCD calculations. Many people are developing new algorithms that integrate machine learning techniques into lattice QCD simulations. The importance of prototyping in this context cannot be overstated, as we aim to create efficient algorithms for lattice QCD that leverage ML techniques. One of the central trade-offs in this process is between code simplicity and execution speed. While ease of implementation is important, a code that is too slow is unusable in practice. Therefore, we focus on optimizing the total time required for lattice QCD studies, which includes both the time spent developing the algorithm and the time taken for its execution.

Portability and parallelizability are also crucial considerations. The code must be able to run efficiently on a wide range of platforms, from laptops to supercomputers, and should be optimized to reduce the overall computational time. This allows researchers to take advantage of the available computational resources and scale their simulations as needed.

The Julia language is a new open-source scientific programming language using LLVM <cit.>. 
Unlike interpreters, Julia is a compiled language with a just-in-time (JIT) compiler. This allows it to achieve performance comparable to Fortran and C <cit.>, while maintaining the simplicity of Python. 
Julia's flexibility and speed make it an ideal choice for developing new lattice QCD algorithms, as it allows for rapid prototyping without sacrificing performance.
Thanks to these advantages, Julia is attractive language for high energy physics <cit.>.

In this paper we introduce a new lattice QCD code JuliaQCD written in Julia language.
The motivation behind the development of this code set is to minimize the combined time of algorithm implementation and execution. Existing languages such as Python are too slow for large-scale lattice QCD simulations, while C++/C/Fortran present challenges in terms of ease of use and development time. Julia strikes the right balance, offering both the speed required for intensive computations and the simplicity needed for fast prototyping. Our goal is to create a flexible and high-performance framework that can be easily adapted to new developments in both lattice QCD and machine learning.



 §.§ JuliaQCD project

JuliaQCD is a suite of first open-source lattice QCD code by Julia programming language. 
(Ultimate) High performance is not within our scope but we are focus on prototyping with good scale.
We developed this code with the following goals:

    
  * Excellent portability 
    
  * Quick and easy setup
    
  * Comprehensive suite 
    
  * Highly modifiable
    
  * Competitive performance 

The system offers excellent portability, as it runs smoothly on any machine with Julia installed, without the need for complicated setup processes. It is also quick and easy to set up, enabling users to start working in less than 5 minutes without requiring complex compilation steps. In addition, the system provides a comprehensive suite that supports both full-QCD and quenched configuration generation, along with a variety of measurement tools. Its highly modifiable structure makes it ideal for rapid prototyping and experimenting with different approaches. Finally, the system delivers competitive performance, comparable to the speed of traditional Fortran 90 codes, ensuring efficiency without sacrificing flexibility.



 §.§ Related Works

In the broader landscape, traditional codes such as the MILC code and the Columbia Physics System (CPS) have also played significant roles. Despite their tradition, these codes remain functional and are used in various research contexts. MILC, in particular, has been a cornerstone for many researchers, supporting a wide range of computations, including those utilizing GPU architectures <cit.>. CPS continues to be a relevant reference point for its C-based implementation.
FermiQCD <cit.> has contributed to the field. FermiQCD, a collection of C++ classes and parallel algorithms, facilitates fast development of lattice applications and includes optimizations for cluster computing, such as SSE2 instructions. Its focus on distributed processing has allowed it to remain an important tool in parallel lattice QCD.

OpenQCD has established itself as a flexible and advanced package for lattice QCD simulations, particularly in high-performance environments. It supports a wide range of theories, including O(a)-improved Wilson quarks, and includes specialized algorithms such as HMC and SMD. openQCD's flexibility allows for large-scale simulations and efficient use of contemporary processors, making it a widely used tool for complex lattice QCD computations.

The development of lattice QCD codes has been significantly influenced by contributions from various international communities, including notable efforts from Japan. Two key projects, the Lattice Tool Kit in Fortran 90 (LTK) <cit.> and Bridge++ <cit.>, originate from this community. LTK, primarily written in Fortran, has served as a foundational tool, though it is not optimized for contemporary high-performance computing (HPC) environments. In contrast, Bridge++ was developed focusing on providing a comprehensive and efficient toolset for lattice QCD simulations on modern supercomputers. Both projects reflect Japan's strong engagement in advancing computational techniques for lattice QCD.

In recent years, there has been a trend towards utilizing C++, like Bridge++, for its advanced features and improved performance in complex simulations. Grid <cit.> is a prime example of this trend, offering robust support for both CPU and GPU computations. The Grid Python Toolkit (GPT) extends this functionality to Python, allowing for easy integration with data analysis and machine learning workflows. Another significant development is SIMULATeQCD <cit.>, which specializes in multi-GPU implementations, enabling high-performance lattice QCD computations. Additionally, QUDA (QCD on CUDA) <cit.> serves as a sub-library that provides specialized routines for running lattice QCD computations on NVIDIA GPUs. By leveraging CUDA architectures, QUDA offers substantial speedups for specific operations within larger QCD packages, rather than functioning as a standalone QCD software suite.

The increasing adoption of Python in scientific computing has led to the development of tools like Lyncs-API and pyQCD <cit.> and also GPT. These projects provide high-level interfaces for lattice QCD calculations, combining Python's user-friendly syntax with the computational efficiency of C++ backends. Lyncs-API, as discussed by <cit.>, offers a flexible API for a variety of QCD computations, while pyQCD similarly integrates Python and C++ for efficient lattice QCD research.

This landscape illustrates the evolution of lattice QCD software from traditional Fortran and C implementations to modern, multi-language frameworks designed for the latest HPC environments. LatticeQCD.jl contributes to this ongoing evolution by leveraging Julia's capabilities, providing a new, efficient, and accessible option for researchers.



§ SOFTWARE DESCRIPTION

We support lattice gauge theory in 4-dimensional Euclidean spacetime.
We cover most of standard technologies.
Also we implement a Wizard for parameter files.
Our code LatticeQCD.jl has following functionalities.



0.44
Gauge fields


  * Optimized SU(2), SU(3)

  * General SU(N_c)

Fermion action


  * Wilson (2 flavor)

  * Staggered fermion (1-8 tastes)

  * Standard Domain-wall (2 flavor)

Measurements
    
        
  * Plaquette
        
  * Polyakov loop
        
  * Chiral condensates (Wilson, staggered)
        
  * Momentum projected pion correlator (Wilson, staggered)
        
  * R× T Wilson loop
        
  * Energy density
        
  * Topological charge (plaquette, clover, and O(a^2) improved definition)
        
  * Load & measurement mode (load and measure all configurations in a directory)
    
    Smearing
    
        
  * Stout
        
  * Gradient flow for a generic action 
    

0.44
Configuration generation algorithms


  * Cold/Hot start for SU(N_c). One instanton configuration for SU(2)

  * Heatbath for SU(N_c) & overrelaxation for a general gauge action

  * Even-odd heatbath for the plaquette action

  * Quenched HMC with SU(N_c) for a general gauge action

  * HMC (2 flavor Wilson) with SU(N_c) with a general gauge action

  * HMC (4 taste staggered fermions) with SU(N_c) with a general gauge action

  * RHMC (any flavor staggered) with SU(N_c) for a general gauge action

  * SU(N_c) stout smeared dynamical fermions 

  * Self-learning HMC with the plaquette action

I/O for gauge configurations
    
        
  * ILDG format (Binary)
        
  * JLD format (Default binary file for Julia, one of HDF5)
        
  * Text file for Bridge++ (Bridgetext)
    




If one specified other than N_f=4, 8 with the staggered fermion HMC, RHMC is automatically used. For a machine with the Apple Silicon, N_f=1-8 is available.
Parallelization is supported by MPI, which can be used in 
https://github.com/akio-tomiya/Gaugefields.jlGaugefields.jl and 
https://github.com/akio-tomiya/LatticeDiracOperators.jlLatticeDiracOperators.jl (see Appendix <ref>).
We perform consistency checks with several references. See
Appendix <ref>.



 §.§ Code Structure

JuliaQCD is a project name and it is constructed by several packages.

    
  * https://github.com/akio-tomiya/LatticeQCD.jlLatticeQCD.jl : LatticeQCD.jl is a wrapper of the following packages. Interface of this package is user friendly.
    
  * https://github.com/akio-tomiya/Wilsonloop.jlWilsonloop.jl: Wilsonloop.jl helps us to treat the Wilson loops and generic Wilson lines in any N_c and dimensions. Wilson lines can is defined symbolically.
    
  * https://github.com/akio-tomiya/Gaugefields.jlGaugefields.jl: Gaugefields.jl is a package for lattice SU(N_c) gauge fields. It handles gauge fields (links), gauge actions with MPI, and autograd. This can generate quenched configurations.
    
  * https://github.com/akio-tomiya/LatticeDiracOperators.jlLatticeDiracOperators.jl: LatticeDiracOperators.jl is a package for Dirac operators and fermions on the lattice. It handles pseudo-fermion fields with various lattice Dirac operators, fermion actions with MPI. This can generate configurations with dynamical fermions.
    
  * https://github.com/akio-tomiya/QCDMeasurements.jlQCDMeasurements.jl: QCDMeasurements.jl is a package for measuring physical quantities. It includes measurements for basic quantities like chiral condensates, plaquettes, pion correlators, and topological charge with several definitions. It also includes the gradient flow with several actions.






§ USAGE OF LATTICEQCD.JL

We show the example code to perform lattice QCD calculations using LatticeQCD.jl.
One can start a lattice QCD calculation in 5 steps.

    
  * Install Julia language from https://julialang.org/downloads/Julialang.org.  
    
  * In Julia REPL, press the  key to enter the package mode and type,
[language=JuliaLocal,style=julia]
add LatticeQCD

    and then press the “return” key. Press the “backspace” key (or “delete” key for Mac) to exit the package mode. One can get the latest version via . All dependencies will be solved automatically.

    
  * Include the package with:
[language=JuliaLocal,style=julia]
using LatticeQCD


    
  * Make a parameter file with the wizard:
[language=JuliaLocal,style=julia]
run_wizard()

    Choose parameters.

    
  * Start the simulation with the created parameter file:
[language=JuliaLocal,style=julia]
run_LQCD("my_parameters.toml")

One will get results.
One can write a parameter file.




 §.§ User interfaces

We support the following two user interfaces for LatticeQCD.jl,

    
  * Julia REPL interface (For beginners, just after the lattice QCD textbook)
    
  * General interface (Experience with another code, for batch jobs, customized purpose)

Usage 1 is explained above. 

For Usage 2, in Julia REPL, press the  key to enter the package mode and type,
[language=JuliaLocal,style=julia]
add LatticeQCD

Then, LatticeQCD.jl will be installed.

The “PARAMETER_FILE” can be created through the wizard. To use the wizard on the shell, please write the following code and save it as , 

[language=JuliaLocal,style=julia]
using LatticeQCD
run_wizard()


Then, one can run the wizard,

[language=JuliaLocal,style=julia]
julia wizard.jl


Please write the following code and save it as ,
[language=JuliaLocal,style=julia]
using LatticeQCD
run_LQCD(ARGS[1])


Then, one can execute it like this,

[language=JuliaLocal,style=julia]
julia run.jl PARAMETER_FILE




§ CONCLUSION AND OUTLOOK

In this work, we introduce JuliaQCD, a new and flexible lattice QCD code set developed in the Julia programming language. By leveraging Julia’s high-performance just-in-time compilation and ease of use, this project aims to strike a balance between rapid prototyping and computational efficiency. JuliaQCD offers a range of functionalities, including support for full-QCD simulations and various smearing techniques, while maintaining scalability across different computational environments from personal systems to high-performance computing clusters.

Our approach is motivated by the need to reduce the time required for both algorithm implementation and execution, allowing researchers to experiment with and refine techniques without sacrificing performance. The portability and ease of setup provided by JuliaQCD make it an ideal tool for researchers working at the intersection of lattice QCD and machine learning, where quick adaptation to new algorithms and methodologies is critical.

While the performance of JuliaQCD is comparable to established Fortran, its focus remains on flexibility and usability, providing a framework that can be easily extended to incorporate modern computational techniques. Future work will explore further optimizations, including deeper integration with GPU architectures and advanced machine learning techniques to further enhance the capabilities of lattice QCD simulations.

JuliaQCD represents an important step in the evolution of lattice QCD software, contributing to the broader trend of adopting more modern, user-friendly languages without compromising on performance. We anticipate that this work will help advance both the development and application of lattice QCD, facilitating new discoveries in quantum chromodynamics and related fields.



The authors thank to the https://github.com/tsuchim/Lattice-Tool-KitLattice Tool Kit written in Fortran 90.
The authors thank to Taku Izubuchi, Issaku Kanamori, Okuto Morikawa, Satoshi Terasaki and Hiromasa Watanabe for comments and contributions. 
The work of authors was partially supported by JSPS KAKENHI Grant Numbers 20K14479,
22K03539, 22H05112, and 22H05111, and MEXT as “Program for Promoting Researches on the Supercomputer
Fugaku” (Simulation for basic science: approaching the new quantum era; Grant Number JPMXP1020230411, and
Search for physics beyond the standard model using large-scale lattice QCD simulation and development of AI
technology toward next-generation lattice QCD; Grant Number JPMXP1020230409).









§ WILSONLOOP.JL

A gauge action is constructed by gauge invariant objects, Wilson loops, in discretized spacetime. Wilsonloop.jl helps us to treat with the Wilson loops and generic Wilson lines in any Nc and dimensions.
Wilsonloop.jl has the following functionalities: 

    
  * From a symbolic definition of Wilson lines, this returns SU(N_c)-valued Wilson lines as objects
    
  * Constructing all staples from given symbolic Wilson lines
    
  * Constructing derivatives of given symbolic Wilson lines (auto-grad for SU(N_c) variables)

We can easily generate a plaquette shown as 
[language=JuliaLocal,style=julia]
println("plaq")
plaq = make_plaq()
display(plaq)

The staple of the plaquette is given as 
[language=JuliaLocal,style=julia]
for μ=1:4
    println("μ = μ")
    staples = make_plaq_staple(μ)
    display(staples)
end

An arbitrary Wilson loop is constructed as 
[language=JuliaLocal,style=julia]
loop = [(1,+1),(2,+1),(1,-1),(2,-1)]
println(loop)
w = Wilsonline(loop)
println("P: ")
show(w)

Its adjoint is calculated as 
[language=JuliaLocal,style=julia]
println("P^+: ")
show(w')

The derivative of the linesdw/dU_μis calculated as 
[language=JuliaLocal,style=julia]
for μ=1:4
    dU = derive_U(w,μ)
    for i=1:length(dU)
        show(dU[i])
    end
end

Note that the derivative is a rank-4 tensor.
The derivatives are usually used for making the smearing of the gauge fields (Stout smearing can be used in Gaugefields.jl). 




§ GAUGEFIELDS.JL


This package has following functionalities

  * SU(N_c) (Nc > 1) gauge fields in 2 or 4 dimensions with arbitrary actions.
  * Z(N_c) 2-form gauge fields in 4 dimensions, which are given as 't Hooft flux.
  * U(1) gauge fields in 2 dimensions with arbitrary actions.
  * Configuration generation
  *  Heatbath
    
  *  quenched Hybrid Monte Carlo
    
  *  quenched Hybrid Monte Carlo being subject to 't Hooft twisted b.c.
with external (non-dynamical) Z(N_c) 2-form gauge fields [Thanks for O. Morikawa.]
  *  quenched Hybrid Monte Carlo for SU(N_c)/Z(N_c) gauge theory
with dynamical Z(N_c) 2-form gauge fields

  * Gradient flow via RK3
  *  Yang-Mills gradient flow
    
  *  Yang-Mills gradient flow being subject to 't Hooft twisted b.c.
    
  *  Gradient flow for SU(N_c)/Z(N_c) gauge theory

  * I/O: ILDG and Bridge++ formats are supported (c-lime will be installed implicitly with CLIME_jll )



  * MPI parallel computation(experimental)


  *  quenched HMC with MPI being subject to 't Hooft twisted b.c.


In addition, this supports followings

  *  Autograd for functions with SU(N_c) variables
    
  *  Stout smearing
    
  *  Stout force via backpropagation

We note that Autograd can be worked for general Wilson lines except for ones have overlaps. 



 §.§ File loading and saving


  §.§.§ ILDG format

We can use ILDG format, one of standard formats for configurations. 
We can read ILDG format as follows:
[language=JuliaLocal,style=julia]
using Gaugefields

NX = 4
NY = 4
NZ = 4
NT = 4
NC = 3
Nwing = 1
Dim = 4
U = Initialize_Gaugefields(NC,Nwing,NX,NY,NZ,NT,condition = "cold")
filename = "hoge.ildg"
ildg = ILDG(filename)
i = 1
L = [NX,NY,NZ,NT]
load_gaugefield!(U,i,ildg,L,NC)

With the use of the configuration, we can calculate the plaquette: 
[language=JuliaLocal,style=julia]
temp1 = similar(U[1])
temp2 = similar(U[1])

comb = 6
factor = 1/(comb*U[1].NV*U[1].NC)
@time plaq_t = calculate_Plaquette(U,temp1,temp2)*factor
println("plaq_t = plaq_t")
poly = calculate_Polyakov_loop(U,temp1,temp2) 
println("polyakov loop =(real(poly)) (imag(poly))")

We can write a configuration as the ILDG format as follows: 
[language=JuliaLocal,style=julia]
filename = "hoge.ildg"
save_binarydata(U,filename)


  §.§.§ Text format for Bridge++

Gaugefields.jl also supports a text format for Bridge++. 
A file loading and saving are expressed as follows: 
[language=JuliaLocal,style=julia]
using Gaugefields
filename = "testconf.txt"
load_BridgeText!(filename,U,L,NC)
[language=JuliaLocal,style=julia]
filename = "testconf.txt"
save_textdata(U,filename)


 §.§ Gradient flow

To smear Gauge fields is important in LatticeQCD. 
We show the codes of the Lüscher's gradient flow as follows.
[language=JuliaLocal,style=julia]
NX = 4
NY = 4
NZ = 4
NT = 4
Nwing = 0
NC = 3

U = Initialize_Gaugefields(NC, Nwing, NX, NY, NZ, NT, condition="hot")

temp1 = similar(U[1])
temp2 = similar(U[1])
comb = 6
factor = 1 / (comb * U[1].NV * U[1].NC)

g = Gradientflow(U)
for itrj = 1:100
    flow!(U, g)
    @time plaq_t = calculate_Plaquette(U, temp1, temp2) * factor
    println("itrj plaq_t =plaq_t")
    poly = calculate_Polyakov_loop(U, temp1, temp2)
    println("itrj polyakov loop =(real(poly)) (imag(poly))")
end


 §.§ Hybrid Monte Carlo

With the use of the Gaugefields.jl, we can easily make the Hybrid Monte Carlo method as following. 
[language=JuliaLocal,style=julia]
using Random
using Gaugefields
using LinearAlgebra

function calc_action(gauge_action, U, p)
    NC = U[1].NC
    Sg = -evaluate_GaugeAction(gauge_action, U) / NC 
    Sp = p * p / 2
    S = Sp + Sg
    return real(S)
end

function MDstep!(gauge_action, U, p, MDsteps, Dim, Uold, temp1, temp2)
    Δτ = 1.0 / MDsteps
    gauss_distribution!(p)
    Sold = calc_action(gauge_action, U, p)
    substitute_U!(Uold, U)

    for itrj = 1:MDsteps
        U_update!(U, p, 0.5, Δτ, Dim, gauge_action)

        P_update!(U, p, 1.0, Δτ, Dim, gauge_action, temp1, temp2)

        U_update!(U, p, 0.5, Δτ, Dim, gauge_action)
    end
    Snew = calc_action(gauge_action, U, p)
    println("Sold = Sold, Snew =Snew")
    println("Snew - Sold = (Snew-Sold)")
    ratio = min(1, exp(-Snew + Sold))
    if rand() > ratio
        substitute_U!(U, Uold)
        return false
    else
        return true
    end
end

function U_update!(U, p, ϵ, Δτ, Dim, gauge_action)
    temps = get_temporary_gaugefields(gauge_action)
    temp1 = temps[1]
    temp2 = temps[2]
    expU = temps[3]
    W = temps[4]
    for μ = 1:Dim
        exptU!(expU, ϵ * Δτ, p[μ], [temp1, temp2])
        mul!(W, expU, U[μ])
        substitute_U!(U[μ], W)
    end
end

function P_update!(U, p, ϵ, Δτ, Dim, gauge_action, temp1, temp2) # p -> p +factor*U*dSdUμ
    NC = U[1].NC
    temp = temp1
    dSdUμ = temp2
    factor = -ϵ * Δτ / (NC)
    for μ = 1:Dim
        calc_dSdUμ!(dSdUμ, gauge_action, μ, U)
        mul!(temp, U[μ], dSdUμ) # U*dSdUμ
        Traceless_antihermitian_add!(p[μ], factor, temp)
    end
end
function HMC_test_4D(NX, NY, NZ, NT, NC, β)
    Dim = 4
    Nwing = 0
    Random.seed!(123)
    U = Initialize_Gaugefields(NC, Nwing, NX, NY, NZ, NT, condition="hot", randomnumber="Reproducible")
    println(typeof(U))
    temp1 = similar(U[1])
    temp2 = similar(U[1])
    comb = 6 #4*3/2
    factor = 1 / (comb * U[1].NV * U[1].NC)
    @time plaq_t = calculate_Plaquette(U, temp1, temp2) * factor
    println("0 plaq_t =plaq_t")
    poly = calculate_Polyakov_loop(U, temp1, temp2)
    println("0 polyakov loop = (real(poly))(imag(poly))")

    gauge_action = GaugeAction(U)
    plaqloop = make_loops_fromname("plaquette")
    append!(plaqloop, plaqloop')
    β = β / 2
    push!(gauge_action, β, plaqloop)

    p = initialize_TA_Gaugefields(U) #This is a traceless-antihermitian gauge fields. This has NC^2-1 real coefficients. 
    Uold = similar(U)
    substitute_U!(Uold, U)
    MDsteps = 100
    numaccepted = 0

    numtrj = 10
    for itrj = 1:numtrj
        t = @timed begin
            accepted = MDstep!(gauge_action, U, p, MDsteps, Dim, Uold, temp1, temp2)
        end
        if get_myrank(U) == 0
            println("elapsed time for MDsteps: (t.time) [s]")
        end
        numaccepted += ifelse(accepted, 1, 0)
        if itrj 
            @time plaq_t = calculate_Plaquette(U, temp1, temp2) * factor
            println("itrj plaq_t = plaq_t")
            poly = calculate_Polyakov_loop(U, temp1, temp2)
            println("itrj polyakov loop = (real(poly))(imag(poly))")
            println("acceptance ratio ", numaccepted / itrj)
        end
    end
    return plaq_t, numaccepted / numtrj
end

function main()
    β = 5.7
    NX = 8
    NY = 8
    NZ = 8
    NT = 8
    NC = 3
    HMC_test_4D(NX, NY, NZ, NT, NC, β)
end
main()



 §.§ HMC with MPI

Here, we show the HMC with MPI. the REPL and Jupyter notebook can not be used when one wants to use MPI. At first, in Julia REPL in the package mode,
[language=JuliaLocal,style=julia]
add MPIPreferences

and
[language=JuliaLocal,style=julia]
using MPIPreferences
MPIPreferences.use_system_binary()

With the use of 
[language=JuliaLocal,style=julia]
add MPI

we can use MPI in Julia.
We show the sample code: 
[language=JuliaLocal,style=julia]

using Random
using Gaugefields
using LinearAlgebra
using MPI

if length(ARGS) < 5
    error("USAGE: ","""
    mpiexecjl -np 2 exe.jl 1 1 1 2 true
    """)
end
const pes = Tuple(parse.(Int64,ARGS[1:4]))
const mpi = parse(Bool,ARGS[5])

function calc_action(gauge_action,U,p)
    NC = U[1].NC
    Sg = -evaluate_GaugeAction(gauge_action,U)/NC #evaluate_Gauge_action(gauge_action,U) = tr(evaluate_Gaugeaction_untraced(gauge_action,U))
    Sp = p*p/2
    S = Sp + Sg
    return real(S)
end

function MDstep!(gauge_action,U,p,MDsteps,Dim,Uold,temp1,temp2)
    Δτ = 1.0/MDsteps
    gauss_distribution!(p)
    Sold = calc_action(gauge_action,U,p)
    substitute_U!(Uold,U)

    for itrj=1:MDsteps
        U_update!(U,p,0.5,Δτ,Dim,gauge_action)

        P_update!(U,p,1.0,Δτ,Dim,gauge_action,temp1,temp2)

        U_update!(U,p,0.5,Δτ,Dim,gauge_action)
    end
    Snew = calc_action(gauge_action,U,p)
    if get_myrank(U) == 0
        println("Sold = Sold, Snew =Snew")
        println("Snew - Sold = (Snew-Sold)")
    end
    ratio = min(1,exp(-Snew+Sold))
    r = rand()
    if mpi
        r = MPI.bcast(r, 0, MPI.COMM_WORLD)
    end
    #ratio = min(1,exp(Snew-Sold))
    if r > ratio
        substitute_U!(U,Uold)
        return false
    else
        return true
    end
end

function U_update!(U,p,ϵ,Δτ,Dim,gauge_action)
    temps = get_temporary_gaugefields(gauge_action)
    temp1 = temps[1]
    temp2 = temps[2]
    expU = temps[3]
    W = temps[4]

    for μ=1:Dim
        exptU!(expU,ϵ*Δτ,p[μ],[temp1,temp2])
        mul!(W,expU,U[μ])
        substitute_U!(U[μ],W)
        
    end
end

function P_update!(U,p,ϵ,Δτ,Dim,gauge_action,temp1,temp2) # p -> p +factor*U*dSdUμ
    NC = U[1].NC
    temp = temp1
    dSdUμ = temp2
    factor =  -ϵ*Δτ/(NC)

    for μ=1:Dim
        calc_dSdUμ!(dSdUμ,gauge_action,μ,U)
        mul!(temp,U[μ],dSdUμ) # U*dSdUμ
        Traceless_antihermitian_add!(p[μ],factor,temp)
    end
end


function HMC_test_4D(NX,NY,NZ,NT,NC,β)
    Dim = 4
    Nwing = 0
    Random.seed!(123)
    if mpi
        PEs = pes#(1,1,1,2)
        U = Initialize_Gaugefields(NC,Nwing,NX,NY,NZ,NT,condition = "hot",mpi=true,PEs = PEs,mpiinit = false) 
    else
        U = Initialize_Gaugefields(NC,Nwing,NX,NY,NZ,NT,condition = "hot")
    end

    if get_myrank(U) == 0
        println(typeof(U))
    end

    temp1 = similar(U[1])
    temp2 = similar(U[1])

    if Dim == 4
        comb = 6 #4*3/2
    elseif Dim == 3
        comb = 3
    elseif Dim == 2
        comb = 1
    else
        error("dimensionDim is not supported")
    end
    
    factor = 1/(comb*U[1].NV*U[1].NC)

    @time plaq_t = calculate_Plaquette(U,temp1,temp2)*factor
    if get_myrank(U) == 0
        println("0 plaq_t = plaq_t")
    end
    poly = calculate_Polyakov_loop(U,temp1,temp2) 
    if get_myrank(U) == 0
        println("0 polyakov loop =(real(poly)) (imag(poly))")
    end

    gauge_action = GaugeAction(U)
    plaqloop = make_loops_fromname("plaquette")
    append!(plaqloop,plaqloop')
    β = β/2
    push!(gauge_action,β,plaqloop)

    p = initialize_TA_Gaugefields(U) #This is a traceless-antihermitian gauge fields. This has NC^2-1 real coefficients. 
    Uold = similar(U)
    substitute_U!(Uold,U)
    MDsteps = 100
    temp1 = similar(U[1])
    temp2 = similar(U[1])
    comb = 6
    factor = 1/(comb*U[1].NV*U[1].NC)
    numaccepted = 0

    numtrj = 100
    for itrj = 1:numtrj
        t = @timed begin
            accepted = MDstep!(gauge_action,U,p,MDsteps,Dim,Uold,temp1,temp2)
        end
        if get_myrank(U) == 0
            println("elapsed time for MDsteps:(t.time) [s]")
        end
        numaccepted += ifelse(accepted,1,0)

        if itrj 
            plaq_t = calculate_Plaquette(U,temp1,temp2)*factor
            if get_myrank(U) == 0
                println("itrj plaq_t =plaq_t")
            end
            poly = calculate_Polyakov_loop(U,temp1,temp2) 
            if get_myrank(U) == 0
                println("itrj polyakov loop =(real(poly)) (imag(poly))")
                println("acceptance ratio ",numaccepted/itrj)
            end
        end
    end
    return plaq_t,numaccepted/numtrj

end

function main()
    β = 5.7
    NX = 8
    NY = 8
    NZ = 8
    NT = 8
    NC = 3
    HMC_test_4D(NX,NY,NZ,NT,NC,β)
end
main()

The command is like:
[language=JuliaLocal,style=julia]
mpiexecjl -np 2 julia mpi_sample.jl 1 1 1 2 true

We can also use MPI in LatticeDiracOperators.jl.  



§ LATTICEDIRACOPERATORS.JL

LatticeDiracOperators.jl handles fermions on a lattice. 
This package have the following functionalities:

  *  Constructing actions and its derivative for staggered fermion with 1-8 tastes with the use of the rational HMC
    
  *  Constructing actions and its derivative for Wilson fermion
    
  *  Constructing actions and its derivative for Standard Domainwall fermion
    
  *  Hybrid Monte Carlo method with fermions

With the use of the Gaugefields.jl, we can also do the HMC with stout smearing. 
This package can be regarded as the additional package of the Gaugefields.jl to treat with lattice fermions.


  §.§.§ Definition of pseudo-fermion fields

The pseudo-fermion fields can be defined as 
[language=JuliaLocal,style=julia]
using Gaugefields
using LatticeDiracOperators

NX = 4
NY = 4
NZ = 4
NT = 4
Nwing = 0
Dim = 4
NC = 3

U = Initialize_4DGaugefields(NC,Nwing,NX,NY,NZ,NT,condition = "cold")
x = Initialize_pseudofermion_fields(U[1],"Wilson")

Here, x is a pseudo fermion fields for Wilson Dirac operator. 
The element of x is x[ic,ix,iy,iz,it,ialpha]. 
ic is an index of the color. ialpha is the internal degree of the gamma matrix. 
The staggered fermions can be defined as 
[language=JuliaLocal,style=julia]
x = Initialize_pseudofermion_fields(U[1],"staggered")

If one wants to obtain the Gaussian distributed pseudo-fermions, the code is written as 
[language=JuliaLocal,style=julia]
gauss_distribution_fermion!(x)


 §.§ Definition of Dirac operators

The Dirac operators are important basic parts in lattice QCD simulations. 
The Wilson Dirac operator can be defined as
[language=JuliaLocal,style=julia]
params = Dict()
params["Dirac_operator"] = "Wilson"
params["κ"] = 0.141139
params["eps_CG"] = 1.0e-8
params["verbose_level"] = 2
D = Dirac_operator(U,x,params)

We can treat the Dirac operator as a matrix. 
Thus, we can apply the Dirac operator to the pseudo-fermion fields as follows.
[language=JuliaLocal,style=julia]
using LinearAlgebra
y = similar(x)
mul!(y,D,x)

And we can solve the equationDx = b: 
[language=JuliaLocal,style=julia]
solve_DinvX!(y,D,x)

The convergence property can be seen by setting "" flag: 
[language=JuliaLocal,style=julia]
params["verbose_level"] = 3
D = Dirac_operator(U,x,params)
gauss_distribution_fermion!(x)
solve_DinvX!(y,D,x)
println(y[1,1,1,1,1,1])

The adjoint of the Dirac operatorD^†andD^† Doperator can be defined as 
[language=JuliaLocal,style=julia]
gauss_distribution_fermion!(x)
solve_DinvX!(y,D',x)
println(y[1,1,1,1,1,1])
DdagD = DdagD_operator(U,x,params)
gauss_distribution_fermion!(x)
solve_DinvX!(y,DdagD,x) 
println(y[1,1,1,1,1,1])

We can similarly define the Dirac operator for the staggered fermion as follows:
[language=JuliaLocal,style=julia]
x = Initialize_pseudofermion_fields(U[1],"staggered")
gauss_distribution_fermion!(x)
params = Dict()
params["Dirac_operator"] = "staggered"
params["mass"] = 0.1
params["eps_CG"] = 1.0e-8
params["verbose_level"] = 2
D = Dirac_operator(U,x,params)
y = similar(x)
mul!(y,D,x)
println(y[1,1,1,1,1,1])
solve_DinvX!(y,D,x)
println(y[1,1,1,1,1,1])

The "tastes" of the staggered fermion is defined in the action.



 §.§ Definition of fermion actions

With the use of the LatticeDiracOperators.jl, we can define actions for pseudo-fermions. 
The sample codes are written as 
[language=JuliaLocal,style=julia]
NX = 4
NY = 4
NZ = 4
NT = 4
Nwing = 0
Dim = 4
NC = 3

U = Initialize_4DGaugefields(NC, Nwing, NX, NY, NZ, NT, condition="cold")
x = Initialize_pseudofermion_fields(U[1], "Wilson")
gauss_distribution_fermion!(x)

params = Dict()
params["Dirac_operator"] = "Wilson"
params["κ"] = 0.141139
params["eps_CG"] = 1.0e-8
params["verbose_level"] = 2

D = Dirac_operator(U, x, params)

parameters_action = Dict()
fermi_action = FermiAction(D, parameters_action)

Then, the fermion action with given pseudo-fermion fields is evaluated as
[language=JuliaLocal,style=julia]
Sfnew = evaluate_FermiAction(fermi_action,U,x)
println(Sfnew)

We can also calculate the derivative of the fermion actiondSf/dUas 
[language=JuliaLocal,style=julia]
UdSfdUμ = calc_UdSfdU(fermi_action,U,x)

The function  calculates theU dSf/dU.
We can also use .

In the case of the staggered fermion, we can choose "taste". 
The action is defined as 
[language=JuliaLocal,style=julia]
x = Initialize_pseudofermion_fields(U[1],"staggered")
gauss_distribution_fermion!(x)
params = Dict()
params["Dirac_operator"] = "staggered"
params["mass"] = 0.1
params["eps_CG"] = 1.0e-8
params["verbose_level"] = 2
D = Dirac_operator(U,x,params)

Nf = 2

println("Nf = Nf")
parameters_action = Dict()
parameters_action["Nf"] = Nf
fermi_action = FermiAction(D,parameters_action)
Sfnew = evaluate_FermiAction(fermi_action,U,x)
println(Sfnew)
UdSfdUμ = calc_UdSfdU(fermi_action,U,x)

This package uses the RHMC techniques to consider the tastes.



 §.§ Hybrid Monte Carlo

We show a sample code for the Hybrid Monte Carlo method with pseudo-fermion fields. 
The codes are written as 
[language=JuliaLocal,style=julia]
using Gaugefields
using LatticeDiracOperators
using LinearAlgebra
using InteractiveUtils
using Random

function MDtest!(gauge_action,U,Dim,fermi_action,η,ξ)
    p = initialize_TA_Gaugefields(U) #This is a traceless-antihermitian gauge fields. This has NC^2-1 real coefficients. 
    Uold = similar(U)
    substitute_U!(Uold,U)
    MDsteps = 10
    temp1 = similar(U[1])
    temp2 = similar(U[1])
    comb = 6
    factor = 1/(comb*U[1].NV*U[1].NC)
    numaccepted = 0
    Random.seed!(123)

    numtrj = 10
    for itrj = 1:numtrj
        @time accepted = MDstep!(gauge_action,U,p,MDsteps,Dim,Uold,fermi_action,η,ξ)
        numaccepted += ifelse(accepted,1,0)

        plaq_t = calculate_Plaquette(U,temp1,temp2)*factor
        println("itrj plaq_t =plaq_t")
        println("acceptance ratio ",numaccepted/itrj)
    end
end

function calc_action(gauge_action,U,p)
    NC = U[1].NC
    Sg = -evaluate_GaugeAction(gauge_action,U)/NC 
    Sp = p*p/2
    S = Sp + Sg
    return real(S)
end


function MDstep!(gauge_action,U,p,MDsteps,Dim,Uold,fermi_action,η,ξ)
    Δτ = 1/MDsteps
    NC,_,NN... = size(U[1])
    
    gauss_distribution!(p)
    
    substitute_U!(Uold,U)
    gauss_sampling_in_action!(ξ,U,fermi_action)
    sample_pseudofermions!(η,U,fermi_action,ξ)
    Sfold = real(dot(ξ,ξ))
    println("Sfold = Sfold")

    Sold = calc_action(gauge_action,U,p) + Sfold
    println("Sold = ",Sold)

    for itrj=1:MDsteps
        U_update!(U,p,0.5,Δτ,Dim,gauge_action)

        P_update!(U,p,1.0,Δτ,Dim,gauge_action)
        P_update_fermion!(U,p,1.0,Δτ,Dim,gauge_action,fermi_action,η)

        U_update!(U,p,0.5,Δτ,Dim,gauge_action)
    end
    Sfnew = evaluate_FermiAction(fermi_action,U,η)
    println("Sfnew =Sfnew")
    Snew = calc_action(gauge_action,U,p) + Sfnew

    println("Sold = Sold, Snew =Snew")
    println("Snew - Sold = (Snew-Sold)")
    accept = exp(Sold - Snew) >= rand()
    if accept != true #rand() > ratio
        substitute_U!(U,Uold)
        return false
    else
        return true
    end
end

function U_update!(U,p,ϵ,Δτ,Dim,gauge_action)
    temps = get_temporary_gaugefields(gauge_action)
    temp1 = temps[1]
    temp2 = temps[2]
    expU = temps[3]
    W = temps[4]

    for μ=1:Dim
        exptU!(expU,ϵ*Δτ,p[μ],[temp1,temp2])
        mul!(W,expU,U[μ])
        substitute_U!(U[μ],W)
        
    end
end

function P_update!(U,p,ϵ,Δτ,Dim,gauge_action) # p -> p +factor*U*dSdUμ
    NC = U[1].NC
    temps = get_temporary_gaugefields(gauge_action)
    dSdUμ = temps[end]
    factor =  -ϵ*Δτ/(NC)

    for μ=1:Dim
        calc_dSdUμ!(dSdUμ,gauge_action,μ,U)
        mul!(temps[1],U[μ],dSdUμ) # U*dSdUμ
        Traceless_antihermitian_add!(p[μ],factor,temps[1])
    end
end

function P_update_fermion!(U,p,ϵ,Δτ,Dim,gauge_action,fermi_action,η)
    temps = get_temporary_gaugefields(gauge_action)
    UdSfdUμ = temps[1:Dim]
    factor =  -ϵ*Δτ
    calc_UdSfdU!(UdSfdUμ,fermi_action,U,η)

    for μ=1:Dim
        Traceless_antihermitian_add!(p[μ],factor,UdSfdUμ[μ])
    end
end

function test1()
    NX = 4
    NY = 4
    NZ = 4
    NT = 4
    Nwing = 0
    Dim = 4
    NC = 3
    U = Initialize_4DGaugefields(NC,Nwing,NX,NY,NZ,NT,condition = "cold")
    gauge_action = GaugeAction(U)
    plaqloop = make_loops_fromname("plaquette")
    append!(plaqloop,plaqloop')
    β = 5.5/2
    push!(gauge_action,β,plaqloop)
    show(gauge_action)
    x = Initialize_pseudofermion_fields(U[1],"Wilson")

    params = Dict()
    params["Dirac_operator"] = "Wilson"
    params["κ"] = 0.141139
    params["eps_CG"] = 1.0e-8
    params["verbose_level"] = 2
    D = Dirac_operator(U,x,params)
    parameters_action = Dict()
    fermi_action = FermiAction(D,parameters_action)
    y = similar(x)
    MDtest!(gauge_action,U,Dim,fermi_action,x,y)
end
test1()

We can easily switch the Wilson fermion to the staggered fermions. 



§ QCDMEASUREMENTS.JL

This package has following functionalities

  *  Plaquette measurement.
    
  *  Polyakov loop measurement.

  *  Pion correlator measurement.

  *  Chiral condensate measurement.

  *  Topological charge measurement.

  *  Energy density measurement.

  *  Wilson loop measurement


 §.§ Sample


To measure variables, one have to make an instance 
[language=JuliaLocal,style=julia]
m_plaq = Plaquette_measurement(U)

By using this , one can measure and get the plaquette: 
[language=JuliaLocal,style=julia]
plaq = get_value(measure(m_plaq,U))


The example is shown as follows. 
[language=JuliaLocal,style=julia]
using QCDMeasurements
using Gaugefields
function test()
    println("SU3test")
    NX = 4
    NY = 4
    NZ = 4
    NT = 4
    Nwing = 0
    Dim = 4
    NC = 3

    U = Initialize_4DGaugefields(NC,Nwing,NX,NY,NZ,NT,condition = "cold")
    filename = "testconf.txt"
    L = [NX,NY,NZ,NT]
    load_BridgeText!(filename,U,L,NC)
    m_plaq = Plaquette_measurement(U)
    m_poly = Polyakov_measurement(U)

    plaq = get_value(measure(m_plaq,U))
    poly = get_value(measure(m_poly,U))
    println("plaq: plaq")
    println("poly:poly")

    m_energy = Energy_density_measurement(U)
    m_topo = Topological_charge_measurement(U)
    energy = get_value(measure(m_energy,U))
    topo = get_value(measure(m_topo,U))
    println("energy: energy")
    println("topo:topo")

    m_wilson = Wilson_loop_measurement(U,printvalues=true)
    wilsonloop = get_value(measure(m_wilson,U))
    println("wilson loop: ",wilsonloop)

    m_pion = Pion_correlator_measurement(U)
    m_pion_Staggered = Pion_correlator_measurement(U,fermiontype = "Staggered")
    m_pion_Wilson = Pion_correlator_measurement(U,fermiontype = "Wilson")
    pion = get_value(measure(m_pion,U))
    pion_s = get_value(measure(m_pion_Staggered,U))
    pion_w = get_value(measure(m_pion_Wilson,U))

    println("pion: pion")
    println("pion correlator with Staggered fermion:pion_s")
    println("pion correlator with  Wilson fermion: pion_w")

    m_chiral_Staggered = Chiral_condensate_measurement(U,fermiontype = "Staggered")
    m_chiral_Wilson = Chiral_condensate_measurement(U,fermiontype = "Wilson")
    chiral_s = get_value(measure(m_chiral_Staggered,U))
    chiral_w = get_value(measure(m_chiral_Wilson,U))

    println("Chiral condensate with Staggered fermion:chiral_s")
    println("Chiral condensatewith  Wilson fermion: chiral_w")


    TC_methods = ["plaquette","clover"]
    m_topo = Topological_charge_measurement(U,TC_methods = TC_methods)
    g = Gradientflow(U)
    for itrj=1:100
        flow!(U,g)
        @time plaq_t = get_value(measure(m_plaq,U))
        @time poly = get_value(measure(m_poly,U))
        println("itrj plaq_t = plaq_t")
        println("itrj polyakov loop = (real(poly))(imag(poly))")

        @time topo = get_value(measure(m_topo,U))
        print("itrj topological charge: ")
        for (key,value) in topo
            print("key value ")
        end
        println("")
    end

end
test()

One can also use the dictionary type:
[language=JuliaLocal,style=julia]
using QCDMeasurements
using Gaugefields
function SU3test()
    println("SU3test")
    NX = 4
    NY = 4
    NZ = 4
    NT = 4
    Nwing = 0
    Dim = 4
    NC = 3

    U = Initialize_4DGaugefields(NC,Nwing,NX,NY,NZ,NT,condition = "cold")
    filename = "testconf.txt"
    L = [NX,NY,NZ,NT]
    load_BridgeText!(filename,U,L,NC)

    method = Dict()
    methodname = "Eigenvalue"
    method["methodname"] = methodname
    method["fermiontype"] = "Wilson"
    κ = 0.141139
    method["hop"] =  κ
    method["nev"] = 1 #number of eigenvalues
    m = prepare_measurement_from_dict(U,method)
    value,vectors = get_value(measure(m,U)) #eigenvalues and eigenvectors
    println("methodnamevalue")
    

    method = Dict()
    methodname = "Pion_correlator"
    method["methodname"] = methodname
    method["fermiontype"] = "Staggered"
    method["mass"] = 1
    method["Nf"] = 4
    m = prepare_measurement_from_dict(U,method)
    value = get_value(measure(m,U))
    println("methodnamevalue")

    method = Dict()
    methodname = "Pion_correlator"
    method["methodname"] = methodname
    method["fermiontype"] = "Wilson"
    method["hop"] = 1
    m = prepare_measurement_from_dict(U,method)
    value = get_value(measure(m,U))
    println("methodnamevalue")


    methodsname = ["Plaquette","Polyakov_loop","Topological_charge","Chiral_condensate",
            "Pion_correlator","Energy_density","Wilson_loop","Eigenvalue"]
    method = Dict()
    for methodname in methodsname
        method["methodname"] = methodname
        m = prepare_measurement_from_dict(U,method)
        value = get_value(measure(m,U))
        if methodname == "Eigenvalue"
            println("methodname(value[1])")
        else
            println("methodname(value)")
        end
    end

end
SU3test()



§ CONSISTENCY CHECK

We compared results of LatticeQCD.jl to the following papers and codes

  * N_f=4 SU(3) staggered HMC <cit.>
  *  Quenched SU(2) improved thermodynamics <cit.>
  *  RHMC <cit.>
  *  HMC for Wilson and Clover Wilson fermions https://github.com/tsuchim/Lattice-Tool-KitLattice Tool Kit
  *  Pion correlator with the Wilson-Dirac operator <cit.>
  *  Pion correlator with the staggered Dirac operator <cit.>