### Install SmoQyElPhQMC using Pkg.add

Source: https://smoqysuite.github.io/SmoQyElPhQMC.jl/stable

Alternatively, install the package programmatically by running 'using Pkg; Pkg.add("SmoQyElPhQMC")' in the Julia REPL.

```julia
julia> using Pkg; Pkg.add("SmoQyElPhQMC")
```

--------------------------------

### Install SmoQyElPhQMC via Julia package manager

Source: https://smoqysuite.github.io/SmoQyElPhQMC.jl/stable

Open the Julia REPL, enter package mode with ']', and run 'add SmoQyElPhQMC' to install the package.

```julia
julia> ]
pkg> add SmoQyElPhQMC
```

--------------------------------

### Initialize simulation parameters and data folder

Source: https://smoqysuite.github.io/SmoQyElPhQMC.jl/stable/tutorials/holstein_honeycomb_checkpoint

Initializes the simulation start time, converts time units, constructs the data folder prefix, and sets up the simulation info and directory.

```julia
    # Record when the simulation began.
    start_timestamp = time()

    # Convert runtime limit from hours to seconds.
    runtime_limit = runtime_limit * 60.0^2

    # Convert checkpoint frequency from hours to seconds.
    checkpoint_freq = checkpoint_freq * 60.0^2

    # Construct the foldername the data will be written to.
    datafolder_prefix = @sprintf "holstein_honeycomb_w%.2f_a%.2f_mu%.2f_L%d_b%.2f" Ω α μ L β

    # Get MPI process ID.
    pID = MPI.Comm_rank(comm)

    # Initialize simulation info.
    simulation_info = SimulationInfo(
        filepath = filepath,
        datafolder_prefix = datafolder_prefix,
        write_bins_concurrent = (L > 7),
        sID = sID,
        pID = pID
    )

    # Initialize the directory the data will be written to.
    initialize_datafolder(comm, simulation_info)
```

--------------------------------

### Setup QMC simulation components

Source: https://smoqysuite.github.io/SmoQyElPhQMC.jl/stable/tutorials/holstein_honeycomb_mpi

Allocates a FermionPathIntegral for both spin-up and down electrons, initializes it with electron-phonon parameters, and sets up the fermion determinant matrix, pseudofermion field calculator, KPM preconditioner, and Green's function estimator. This code is unchanged from the previous 1a) Honeycomb Holstein Model tutorial.

```julia
    # Allocate a single FermionPathIntegral for both spin-up and down electrons.
    fermion_path_integral = FermionPathIntegral(tight_binding_parameters = tight_binding_parameters, β = β, Δτ = Δτ)

    # Initialize FermionPathIntegral type to account for electron-phonon interaction.
    initialize!(fermion_path_integral, electron_phonon_parameters)

    # Initialize fermion determinant matrix. Also set the default tolerance and max iteration count
    # used in conjugate gradient (CG) solves of linear systems involving this matrix.
    fermion_det_matrix = SymFermionDetMatrix(
        fermion_path_integral,
        maxiter = maxiter, tol = tol
    )

    # Initialize pseudofermion field calculator.
    pff_calculator = PFFCalculator(electron_phonon_parameters, fermion_det_matrix)

    # Initialize KPM preconditioner.
    preconditioner = KPMPreconditioner(fermion_det_matrix, rng = rng)

    # Initialize Green's function estimator for making measurements.
    greens_estimator = GreensEstimator(fermion_det_matrix, model_geometry)
```

--------------------------------

### Run DQMC simulation command

Source: https://smoqysuite.github.io/SmoQyElPhQMC.jl/stable/tutorials/holstein_honeycomb

Example command to execute the Holstein model simulation with specific parameters for the honeycomb lattice.

```bash
> julia holstein_honeycomb.jl 1 1.0 1.5 0.0 3 4.0 5000 10000 100
```

--------------------------------

### Initialize simulation metadata

Source: https://smoqysuite.github.io/SmoQyElPhQMC.jl/stable/tutorials/holstein_honeycomb_density_tuning

Creates a metadata dictionary to store simulation parameters and tracking variables when starting a new simulation.

```julia
    # If starting a new simulation i.e. not resuming a previous simulation.
    if !simulation_info.resuming

        # Begin thermalization updates from start.
        n_therm = 1

        # Begin measurement updates from start.
        n_measurements = 1

        # Initialize random number generator
        rng = Xoshiro(seed)

        # Initialize metadata dictionary
        metadata = Dict()

        # Record simulation parameters.
        metadata["N_therm"] = N_therm  # Number of thermalization updates
        metadata["N_measurements"] = N_measurements  # Total number of measurements and measurement updates
        metadata["N_bins"] = N_bins # Number of times bin-averaged measurements are written to file
        metadata["maxiter"] = maxiter # Maximum number of conjugate gradient iterations
        metadata["tol"] = tol # Tolerance used for conjugate gradient solves
        metadata["Nt"] = Nt # Number of time-steps in HMC update
        metadata["Nrv"] = Nrv # Number of random vectors used to estimate fermionic correlation functions
        metadata["seed"] = seed  # Random seed used to initialize random number generator in simulation
        metadata["hmc_acceptance_rate"] = 0.0 # HMC acceptance rate
        metadata["reflection_acceptance_rate"] = 0.0 # Reflection update acceptance rate
        metadata["swap_acceptance_rate"] = 0.0 # Swap update acceptance rate
        metadata["hmc_iters"] = 0.0 # Avg number of CG iterations per solve in HMC update.
        metadata["reflection_iters"] = 0.0 # Avg number of CG iterations per solve in reflection update.
        metadata["swap_iters"] = 0.0 # Avg number of CG iterations per solve in swap update.
        metadata["measurement_iters"] = 0.0 # Avg number of CG iterations per solve while making measurements.
```

--------------------------------

### Import packages for density tuning simulation

Source: https://smoqysuite.github.io/SmoQyElPhQMC.jl/stable/tutorials/holstein_honeycomb_density_tuning

Imports required packages including SmoQyElPhQMC, SmoQyDQMC, and specific utilities for lattice and chemical potential tuning. This is the first step in the example, adding MuTuner.jl which is reexported by SmoQyDQMC.jl.

```julia
using SmoQyElPhQMC
using SmoQyDQMC
import SmoQyDQMC.LatticeUtilities as lu
import SmoQyDQMC.MuTuner as mt

using Random
using Printf
using MPI
```

--------------------------------

### Initialize simulation metadata for new runs

Source: https://smoqysuite.github.io/SmoQyElPhQMC.jl/stable/tutorials/holstein_honeycomb_checkpoint

Handles the logic for starting a new simulation by initializing thermalization counters, the random number generator, and the metadata dictionary.

```julia
    # If starting a new simulation i.e. not resuming a previous simulation.
    if !simulation_info.resuming

        # Begin thermalization updates from start.
        n_therm = 1

        # Begin measurement updates from start.
        n_measurements = 1

        # Initialize random number generator
        rng = Xoshiro(seed)

        # Initialize metadata dictionary
        metadata = Dict()

        # Record simulation parameters.
        metadata["N_therm"] = N_therm  # Number of thermalization updates
        metadata["N_measurements"] = N_measurements  # Total number of measurements
        metadata["N_bins"] = N_bins # Number of times bin-averaged measurements are written to file
        metadata["maxiter"] = maxiter # Maximum number of conjugate gradient iterations
        metadata["tol"] = tol # Tolerance used for conjugate gradient solves
        metadata["Nt"] = Nt # Number of time-steps in HMC update
        metadata["Nrv"] = Nrv # Number of random vectors used to estimate fermionic correlation functions
        metadata["seed"] = seed  # Random seed used to initialize random number generator in simulation
        metadata["hmc_acceptance_rate"] = 0.0 # HMC acceptance rate
        metadata["reflection_acceptance_rate"] = 0.0 # Reflection update acceptance rate
        metadata["swap_acceptance_rate"] = 0.0 # Swap update acceptance rate
        metadata["hmc_iters"] = 0.0 # Avg number of CG iterations per solve in HMC update.
        metadata["reflection_iters"] = 0.0 # Avg number of CG iterations per solve in reflection update.
        metadata["swap_iters"] = 0.0 # Avg number of CG iterations per solve in swap update.
        metadata["measurement_iters"] = 0.0 # Avg number of CG iterations per solve while making measurements.
```

--------------------------------

### Run script with mpiexecjl

Source: https://smoqysuite.github.io/SmoQyElPhQMC.jl/stable/tutorials/holstein_honeycomb_mpi

Example command to run the script with 16 MPI processes using mpiexecjl. Each process runs an independent simulation with a different random seed, and results are averaged over all walkers.

```bash
mpiexecjl -n 16 julia holstein_honeycomb_mpi.jl 1 1.0 1.5 0.0 3 4.0 5000 10000 100
```

--------------------------------

### Run simulation with mpiexecjl

Source: https://smoqysuite.github.io/SmoQyElPhQMC.jl/stable/tutorials/holstein_honeycomb_checkpoint

Command to launch the simulation with 16 MPI processes using mpiexecjl. The 10th argument (0.5) sets checkpoint_freq to 0.5 hours, so a new checkpoint is written if more than 30 minutes have passed since the last one. The same command is used to start a new simulation or resume a previous one.

```bash
mpiexecjl -n 16 julia holstein_honeycomb_checkpoint.jl 1 1.0 1.5 0.0 3 4.0 5000 10000 100 0.5
```

--------------------------------

### Initialize simulation estimators and updater

Source: https://smoqysuite.github.io/SmoQyElPhQMC.jl/stable/examples/bssh_chain

Sets up the pseudofermion calculator, KPM preconditioner, Green's function estimator, and the HMC updater.

```julia
# Initialize pseudofermion field calculator.
pff_calculator = PFFCalculator(electron_phonon_parameters, fermion_det_matrix)

# Initialize KPM preconditioner.
preconditioner = KPMPreconditioner(fermion_det_matrix, rng = rng)

# Initialize Green's function estimator for making measurements.
greens_estimator = GreensEstimator(fermion_det_matrix, model_geometry)

# Initialize Hamiltonian/Hybrid monte carlo (HMC) updater.
hmc_updater = EFAPFFHMCUpdater(
    electron_phonon_parameters = electron_phonon_parameters,
    Nt = Nt, Δt = π/(2*Nt)
)
```

--------------------------------

### Perform thermalization updates

Source: https://smoqysuite.github.io/SmoQyElPhQMC.jl/stable/examples/bssh_chain

Starts the loop to perform the specified number of thermalization updates.

```julia
# Iterate over number of thermalization updates to perform.
for update in n_therm:N_therm
```

--------------------------------

### Initialize model parameters and MuTunerLogger

Source: https://smoqysuite.github.io/SmoQyElPhQMC.jl/stable/tutorials/holstein_honeycomb_density_tuning

Initializes tight-binding and electron-phonon parameters, then sets up the MuTunerLogger using the total number of orbitals calculated via LatticeUtilities.nsites.

```julia
        # Initialize tight-binding parameters.
        tight_binding_parameters = TightBindingParameters(
            tight_binding_model = tight_binding_model,
            model_geometry = model_geometry,
            rng = rng
        )

        # Initialize electron-phonon parameters.
        electron_phonon_parameters = ElectronPhononParameters(
            β = β, Δτ = Δτ,
            electron_phonon_model = electron_phonon_model,
            tight_binding_parameters = tight_binding_parameters,
            model_geometry = model_geometry,
            rng = rng
        )

        # Initialize MuTunerLogger type that will be used to dynamically adjust the
        # chemical potential during the simulation.
        chemical_potential_tuner = mt.init_mutunerlogger(
            target_density = n,
            inverse_temperature = β,
            system_size = lu.nsites(unit_cell, lattice),
            initial_chemical_potential = μ,
            complex_sign_problem = false
        )
```

--------------------------------

### Initialize simulation environment

Source: https://smoqysuite.github.io/SmoQyElPhQMC.jl/stable/tutorials/holstein_honeycomb_density_tuning

Sets up timing, converts time units, constructs the data folder path, and initializes the simulation information object for MPI-based execution.

```julia
    # Record when the simulation began.
    start_timestamp = time()

    # Convert runtime limit from hours to seconds.
    runtime_limit = runtime_limit * 60.0^2

    # Convert checkpoint frequency from hours to seconds.
    checkpoint_freq = checkpoint_freq * 60.0^2

    # Construct the foldername the data will be written to.
    datafolder_prefix = @sprintf "holstein_honeycomb_w%.2f_a%.2f_n%.2f_L%d_b%.2f" Ω α n L β

    # Get MPI process ID.
    pID = MPI.Comm_rank(comm)

    # Initialize simulation info.
    simulation_info = SimulationInfo(
        filepath = filepath,
        datafolder_prefix = datafolder_prefix,
        write_bins_concurrent = (L > 7),
        sID = sID,
        pID = pID
    )

    # Initialize the directory the data will be written to.
    initialize_datafolder(comm, simulation_info)
```

--------------------------------

### Run OSSH simulation with SmoQyElPhQMC

Source: https://smoqysuite.github.io/SmoQyElPhQMC.jl/stable/examples/bssh_chain

This script defines the run_simulation function to initialize and execute the OSSH model simulation, including MPI communicator setup and metadata tracking.

```julia
using SmoQyElPhQMC
using SmoQyDQMC
import SmoQyDQMC.LatticeUtilities as lu

using Random
using Printf
using MPI

# Top-level function to run simulation.
function run_simulation(
    comm::MPI.Comm; # MPI communicator.
    # KEYWORD ARGUMENTS
    sID, # Simulation ID.
    Ω, # Phonon energy.
    α, # Electron-phonon coupling.
    μ, # Chemical potential.
    L, # System size.
    β, # Inverse temperature.
    N_therm, # Number of thermalization updates.
    N_measurements, # Total number of measurements and measurement updates.
    N_bins, # Number of times bin-averaged measurements are written to file.
    checkpoint_freq, # Frequency with which checkpoint files are written in hours.
    runtime_limit = Inf, # Simulation runtime limit in hours.
    Δτ = 0.05, # Discretization in imaginary time.
    Nt = 24, # Number of time-steps in HMC update.
    Nrv = 10, # Number of random vectors used to estimate fermionic correlation functions.
    tol = 1e-10, # CG iterations tolerance.
    maxiter = 10_000, # Maximum number of CG iterations.
    seed = abs(rand(Int)), # Seed for random number generator.
    filepath = "." # Filepath to where data folder will be created.
)

    # Record when the simulation began.
    start_timestamp = time()

    # Convert runtime limit from hours to seconds.
    runtime_limit = runtime_limit * 60.0^2

    # Convert checkpoint frequency from hours to seconds.
    checkpoint_freq = checkpoint_freq * 60.0^2

    # Construct the foldername the data will be written to.
    datafolder_prefix = @sprintf "bssh_chain_w%.2f_a%.2f_mu%.2f_L%d_b%.2f" Ω α μ L β

    # Get MPI process ID.
    pID = MPI.Comm_rank(comm)

    # Initialize simulation info.
    simulation_info = SimulationInfo(
        filepath = filepath,
        datafolder_prefix = datafolder_prefix,
        write_bins_concurrent = (L > 100),
        sID = sID,
        pID = pID
    )

    # Initialize the directory the data will be written to.
    initialize_datafolder(comm, simulation_info)

    # If starting a new simulation i.e. not resuming a previous simulation.
    if !simulation_info.resuming

        # Begin thermalization updates from start.
        n_therm = 1

        # Begin measurement updates from start.
        n_measurements = 1

        # Initialize random number generator
        rng = Xoshiro(seed)

        # Initialize metadata dictionary
        metadata = Dict()

        # Record simulation parameters.
        metadata["Nt"] = Nt
        metadata["N_therm"] = N_therm
        metadata["N_measurements"] = N_measurements
        metadata["N_bins"] = N_bins
        metadata["Nrv"] = Nrv
        metadata["maxiter"] = maxiter
        metadata["tol"] = tol
        metadata["seed"] = seed
        metadata["hmc_acceptance_rate"] = 0.0
        metadata["reflection_acceptance_rate"] = 0.0
        metadata["swap_acceptance_rate"] = 0.0
        metadata["hmc_iters"] = 0.0
        metadata["reflection_iters"] = 0.0
        metadata["swap_iters"] = 0.0
        metadata["measurement_iters"] = 0.0
```

--------------------------------

### Initialize Holstein Honeycomb Model (Julia)

Source: https://smoqysuite.github.io/SmoQyElPhQMC.jl/stable/tutorials/holstein_honeycomb_mpi

Defines the honeycomb lattice, tight-binding model, phonon modes, and Holstein couplings. No changes from the previous 1a) tutorial. Requires variables L, μ, Ω, α, β, Δτ, and simulation_info to be defined beforehand.

```julia
    # Define lattice vectors.
    a1 = [+3/2, +√3/2]
    a2 = [+3/2, -√3/2]

    # Define basis vectors for two orbitals in the honeycomb unit cell.
    r1 = [0.0, 0.0] # Location of first orbital in unit cell.
    r2 = [1.0, 0.0] # Location of second orbital in unit cell.

    # Define the unit cell.
    unit_cell = lu.UnitCell(
        lattice_vecs = [a1, a2],
        basis_vecs   = [r1, r2]
    )

    # Define finite lattice with periodic boundary conditions.
    lattice = lu.Lattice(
        L = [L, L],
        periodic = [true, true]
    )

    # Initialize model geometry.
    model_geometry = ModelGeometry(unit_cell, lattice)

    # Define the first nearest-neighbor bond in a honeycomb lattice.
    bond_1 = lu.Bond(orbitals = (1,2), displacement = [0,0])

    # Add the first nearest-neighbor bond in a honeycomb lattice to the model.
    bond_1_id = add_bond!(model_geometry, bond_1)

    # Define the second nearest-neighbor bond in a honeycomb lattice.
    bond_2 = lu.Bond(orbitals = (1,2), displacement = [-1,0])

    # Add the second nearest-neighbor bond in a honeycomb lattice to the model.
    bond_2_id = add_bond!(model_geometry, bond_2)

    # Define the third nearest-neighbor bond in a honeycomb lattice.
    bond_3 = lu.Bond(orbitals = (1,2), displacement = [0,-1])

    # Add the third nearest-neighbor bond in a honeycomb lattice to the model.
    bond_3_id = add_bond!(model_geometry, bond_3)

    # Set nearest-neighbor hopping amplitude to unity,
    # setting the energy scale in the model.
    t = 1.0

    # Define the honeycomb tight-binding model.
    tight_binding_model = TightBindingModel(
        model_geometry = model_geometry,
        t_bonds        = [bond_1, bond_2, bond_3], # defines hopping
        t_mean         = [t, t, t], # defines corresponding hopping amplitude
        μ              = μ, # set chemical potential
        ϵ_mean         = [0.0, 0.0] # set the (mean) on-site energy
    )

    # Initialize a null electron-phonon model.
    electron_phonon_model = ElectronPhononModel(
        model_geometry = model_geometry,
        tight_binding_model = tight_binding_model
    )

    # Define a dispersionless electron-phonon mode to live on each site in the lattice.
    phonon_1 = PhononMode(
        basis_vec = r1,
        Ω_mean = Ω
    )

    # Add the phonon mode definition to the electron-phonon model.
    phonon_1_id = add_phonon_mode!(
        electron_phonon_model = electron_phonon_model,
        phonon_mode = phonon_1
    )

    # Define a dispersionless electron-phonon mode to live on the second sublattice.
    phonon_2 = PhononMode(
        basis_vec = r2,
        Ω_mean = Ω
    )

    # Add the phonon mode definition to the electron-phonon model.
    phonon_2_id = add_phonon_mode!(
        electron_phonon_model = electron_phonon_model,
        phonon_mode = phonon_2
    )

    # Define first local Holstein coupling for first phonon mode.
    holstein_coupling_1 = HolsteinCoupling(
        model_geometry = model_geometry,
        phonon_id = phonon_1_id,
        orbital_id = 1,
        displacement = [0, 0],
        α_mean = α,
        ph_sym_form = true,
    )

    # Add the first local Holstein coupling definition to the model.
    holstein_coupling_1_id = add_holstein_coupling!(
        electron_phonon_model = electron_phonon_model,
        holstein_coupling = holstein_coupling_1,
        model_geometry = model_geometry
    )

    # Define second local Holstein coupling for second phonon mode.
    holstein_coupling_2 = HolsteinCoupling(
        model_geometry = model_geometry,
        phonon_id = phonon_2_id,
        orbital_id = 2,
        displacement = [0, 0],
        α_mean = α,
        ph_sym_form = true,
    )

    # Add the second local Holstein coupling definition to the model.
    holstein_coupling_2_id = add_holstein_coupling!(
        electron_phonon_model = electron_phonon_model,
        holstein_coupling = holstein_coupling_2,
        model_geometry = model_geometry
    )

    # Write model summary TOML file specifying Hamiltonian that will be simulated.
    model_summary(
        simulation_info = simulation_info,
        β = β, Δτ = Δτ,
        model_geometry = model_geometry,
        tight_binding_model = tight_binding_model,
        interactions = (electron_phonon_model,)
    )
```

--------------------------------

### Initialize Measurements for Honeycomb Holstein Model

Source: https://smoqysuite.github.io/SmoQyElPhQMC.jl/stable/tutorials/holstein_honeycomb_mpi

Initializes the measurement container and sets up various correlation measurements including electron and phonon Green's functions, density, pair, spin-z, and composite measurements. No changes are needed from the previous tutorial.

```julia
    # Initialize the container that measurements will be accumulated into.
    measurement_container = initialize_measurement_container(model_geometry, β, Δτ)

    # Initialize the tight-binding model related measurements, like the hopping energy.
    initialize_measurements!(measurement_container, tight_binding_model)

    # Initialize the electron-phonon interaction related measurements.
    initialize_measurements!(measurement_container, electron_phonon_model)

    # Initialize the single-particle electron Green's function measurement.
    initialize_correlation_measurements!(
        measurement_container = measurement_container,
        model_geometry = model_geometry,
        correlation = "greens",
        time_displaced = true,
        pairs = [
            # Measure green's functions for all pairs or orbitals.
            (1, 1), (2, 2), (1, 2)
        ]
    )

    # Initialize the single-particle electron Green's function measurement.
    initialize_correlation_measurements!(
        measurement_container = measurement_container,
        model_geometry = model_geometry,
        correlation = "phonon_greens",
        time_displaced = true,
        pairs = [
            # Measure green's functions for all pairs of modes.
            (1, 1), (2, 2), (1, 2)
        ]
    )

    # Initialize density correlation function measurement.
    initialize_correlation_measurements!(
        measurement_container = measurement_container,
        model_geometry = model_geometry,
        correlation = "density",
        time_displaced = false,
        integrated = true,
        pairs = [
            (1, 1), (2, 2),
        ]
    )

    # Initialize the pair correlation function measurement.
    initialize_correlation_measurements!(
        measurement_container = measurement_container,
        model_geometry = model_geometry,
        correlation = "pair",
        time_displaced = false,
        integrated = true,
        pairs = [
            # Measure local s-wave pair susceptibility associated with
            # each orbital in the unit cell.
            (1, 1), (2, 2)
        ]
    )

    # Initialize the spin-z correlation function measurement.
    initialize_correlation_measurements!(
        measurement_container = measurement_container,
        model_geometry = model_geometry,
        correlation = "spin_z",
        time_displaced = false,
        integrated = true,
        pairs = [
            (1, 1), (2, 2)
        ]
    )

    # Initialize measurement of electron Green's function traced
    # over both orbitals in the unit cell.
    initialize_composite_correlation_measurement!(
        measurement_container = measurement_container,
        model_geometry = model_geometry,
        name = "tr_greens",
        correlation = "greens",
        id_pairs = [(1,1), (2,2)],
        coefficients = [1.0, 1.0],
        time_displaced = true,
    )

    # Initialize CDW correlation measurement.
    initialize_composite_correlation_measurement!(
        measurement_container = measurement_container,
        model_geometry = model_geometry,
        name = "cdw",
        correlation = "density",
        ids = [1, 2],
        coefficients = [1.0, -1.0],
        time_displaced = false,
        integrated = true
    )
```

--------------------------------

### Initialize tight-binding and electron-phonon models

Source: https://smoqysuite.github.io/SmoQyElPhQMC.jl/stable/examples/ossh_honeycomb

Sets the hopping amplitude and constructs the tight-binding model, followed by initializing the null electron-phonon model.

```julia
# Set nearest-neighbor hopping amplitude to unity,
# setting the energy scale in the model.
t = 1.0

# Define the honeycomb tight-binding model.
tight_binding_model = TightBindingModel(
    model_geometry = model_geometry,
    t_bonds = [bond_AB_1, bond_AB_2, bond_AB_3],
    t_mean = [t, t, t],
    μ  = μ,
    ϵ_mean = [0.0, 0.0]
)

# Initialize a null electron-phonon model.
electron_phonon_model = ElectronPhononModel(
    model_geometry = model_geometry,
    tight_binding_model = tight_binding_model
)
```

--------------------------------

### Run simulation function for BSSH square model

Source: https://smoqysuite.github.io/SmoQyElPhQMC.jl/stable/examples/bssh_square

This function initializes the simulation environment, including MPI communicator setup, data folder creation, and simulation metadata for the BSSH square lattice model.

```julia
using SmoQyElPhQMC
using SmoQyDQMC
import SmoQyDQMC.LatticeUtilities as lu

using Random
using Printf
using MPI

# Top-level function to run simulation.
function run_simulation(
    comm::MPI.Comm; # MPI communicator.
    # KEYWORD ARGUMENTS
    sID, # Simulation ID.
    Ω, # Phonon energy.
    α, # Electron-phonon coupling.
    μ, # Chemical potential.
    L, # System size.
    β, # Inverse temperature.
    N_therm, # Number of thermalization updates.
    N_measurements, # Total number of measurements and measurement updates.
    N_bins, # Number of times bin-averaged measurements are written to file.
    checkpoint_freq, # Frequency with which checkpoint files are written in hours.
    runtime_limit = Inf, # Simulation runtime limit in hours.
    Δτ = 0.05, # Discretization in imaginary time.
    Nt = 24, # Number of time-steps in HMC update.
    Nrv = 10, # Number of random vectors used to estimate fermionic correlation functions.
    tol = 1e-10, # CG iterations tolerance.
    maxiter = 10_000, # Maximum number of CG iterations.
    seed = abs(rand(Int)), # Seed for random number generator.
    filepath = "." # Filepath to where data folder will be created.
)

    # Record when the simulation began.
    start_timestamp = time()

    # Convert runtime limit from hours to seconds.
    runtime_limit = runtime_limit * 60.0^2

    # Convert checkpoint frequency from hours to seconds.
    checkpoint_freq = checkpoint_freq * 60.0^2

    # Construct the foldername the data will be written to.
    datafolder_prefix = @sprintf "bssh_square_w%.2f_a%.2f_mu%.2f_L%d_b%.2f" Ω α μ L β

    # Get MPI process ID.
    pID = MPI.Comm_rank(comm)

    # Initialize simulation info.
    simulation_info = SimulationInfo(
        filepath = filepath,
        datafolder_prefix = datafolder_prefix,
        write_bins_concurrent = (L > 100),
        sID = sID,
        pID = pID
    )

    # Initialize the directory the data will be written to.
    initialize_datafolder(comm, simulation_info)

    # If starting a new simulation i.e. not resuming a previous simulation.
    if !simulation_info.resuming

        # Begin thermalization updates from start.
        n_therm = 1

        # Begin measurement updates from start.
        n_measurements = 1

        # Initialize random number generator
        rng = Xoshiro(seed)

        # Initialize metadata dictionary
        metadata = Dict()
```

--------------------------------

### Initialize KPMPreconditioner

Source: https://smoqysuite.github.io/SmoQyElPhQMC.jl/stable/tutorials/holstein_honeycomb

Initializes the KPMPreconditioner to accelerate the convergence of conjugate gradient calculations for the fermionic action.

```julia
    # Initialize KPM preconditioner.
    preconditioner = KPMPreconditioner(fermion_det_matrix, rng = rng)
```

--------------------------------

### Run script with srun on SLURM

Source: https://smoqysuite.github.io/SmoQyElPhQMC.jl/stable/tutorials/holstein_honeycomb_mpi

Example command to run the script inside a SLURM job file using srun. The command automatically detects the number of available cores requested by the job and runs the script with the appropriate number of MPI processes.

```bash
srun julia holstein_honeycomb_mpi.jl 1 1.0 1.5 0.0 3 4.0 5000 10000 100
```

--------------------------------

### Initialize Honeycomb SSH Simulation Components

Source: https://smoqysuite.github.io/SmoQyElPhQMC.jl/stable/examples/ossh_honeycomb

Sets up the necessary simulation structures including the FermionPathIntegral, determinant matrix, PFF calculator, preconditioner, and HMC updater.

```julia
# Unpack contents of checkpoint dictionary.
tight_binding_parameters = checkpoint["tight_binding_parameters"]
electron_phonon_parameters = checkpoint["electron_phonon_parameters"]
measurement_container = checkpoint["measurement_container"]
model_geometry = checkpoint["model_geometry"]
metadata = checkpoint["metadata"]
rng = checkpoint["rng"]
n_therm = checkpoint["n_therm"]
n_measurements = checkpoint["n_measurements"]
end

# Allocate a single FermionPathIntegral for both spin-up and down electrons.
fermion_path_integral = FermionPathIntegral(tight_binding_parameters = tight_binding_parameters, β = β, Δτ = Δτ)

# Initialize FermionPathIntegral type to account for electron-phonon interaction.
initialize!(fermion_path_integral, electron_phonon_parameters)

# Initialize fermion determinant matrix. Also set the default tolerance and max iteration count
# used in conjugate gradient (CG) solves of linear systems involving this matrix.
fermion_det_matrix = SymFermionDetMatrix(
    fermion_path_integral,
    maxiter = maxiter, tol = tol
)

# Initialize pseudofermion field calculator.
pff_calculator = PFFCalculator(electron_phonon_parameters, fermion_det_matrix)

# Initialize KPM preconditioner.
preconditioner = KPMPreconditioner(fermion_det_matrix, rng = rng)

# Initialize Green's function estimator for making measurements.
greens_estimator = GreensEstimator(fermion_det_matrix, model_geometry)

# Initialize Hamiltonian/Hybrid monte carlo (HMC) updater.
hmc_updater = EFAPFFHMCUpdater(
    electron_phonon_parameters = electron_phonon_parameters,
    Nt = Nt, Δt = π/(2*Nt)
)
```

--------------------------------

### Run simulation with srun

Source: https://smoqysuite.github.io/SmoQyElPhQMC.jl/stable/tutorials/holstein_honeycomb_checkpoint

Alternative command to launch the simulation using srun, typically on a cluster. The 10th argument (0.5) sets checkpoint_freq to 0.5 hours, so a new checkpoint is written if more than 30 minutes have passed since the last one. The same command is used to start a new simulation or resume a previous one.

```bash
srun julia holstein_honeycomb_checkpoint.jl 1 1.0 1.5 0.0 3 4.0 5000 10000 100 0.5
```

--------------------------------

### Initialize Parameters and Measurements

Source: https://smoqysuite.github.io/SmoQyElPhQMC.jl/stable/examples/bssh_chain

Initializes the simulation parameters and sets up the measurement container for tracking observables.

```julia
# Initialize tight-binding parameters.
tight_binding_parameters = TightBindingParameters(
    tight_binding_model = tight_binding_model,
    model_geometry = model_geometry,
    rng = rng
)

# Initialize electron-phonon parameters.
electron_phonon_parameters = ElectronPhononParameters(
    β = β, Δτ = Δτ,
    electron_phonon_model = electron_phonon_model,
    tight_binding_parameters = tight_binding_parameters,
    model_geometry = model_geometry,
    rng = rng
)

# Initialize the container that measurements will be accumulated into.
measurement_container = initialize_measurement_container(model_geometry, β, Δτ)

# Initialize the tight-binding model related measurements, like the hopping energy.
initialize_measurements!(measurement_container, tight_binding_model)

# Initialize the electron-phonon interaction related measurements.
initialize_measurements!(measurement_container, electron_phonon_model)

# Initialize the single-particle electron Green's function measurement.
initialize_correlation_measurements!(
    measurement_container = measurement_container,
    model_geometry = model_geometry,
    correlation = "greens",
    time_displaced = true,
    pairs = [
        # Measure green's functions for all pairs or orbitals.
        (1, 1),
    ]
)
```

--------------------------------

### Initialize Simulation with SimulationInfo and data folder

Source: https://smoqysuite.github.io/SmoQyElPhQMC.jl/stable/tutorials/holstein_honeycomb

Constructs the folder name based on simulation parameters and initializes the SimulationInfo object. Sets write_bins_concurrent based on lattice size L to handle small vs large system I/O.

```julia
    # Construct the foldername the data will be written to.
    datafolder_prefix = @sprintf "holstein_honeycomb_w%.2f_a%.2f_mu%.2f_L%d_b%.2f" Ω α μ L β

    # Initialize simulation info.
    simulation_info = SimulationInfo(
        filepath = filepath,
        datafolder_prefix = datafolder_prefix,
        write_bins_concurrent = (L > 7),
        sID = sID
    )

    # Initialize the directory the data will be written to.
    initialize_datafolder(simulation_info)
```

--------------------------------

### Initialize simulation components

Source: https://smoqysuite.github.io/SmoQyElPhQMC.jl/stable/examples/ossh_chain

Allocates path integrals, fermion determinant matrices, calculators, and updaters required for the simulation.

```julia
# Allocate a single FermionPathIntegral for both spin-up and down electrons.
fermion_path_integral = FermionPathIntegral(tight_binding_parameters = tight_binding_parameters, β = β, Δτ = Δτ)

# Initialize FermionPathIntegral type to account for electron-phonon interaction.
initialize!(fermion_path_integral, electron_phonon_parameters)

# Initialize fermion determinant matrix. Also set the default tolerance and max iteration count
# used in conjugate gradient (CG) solves of linear systems involving this matrix.
fermion_det_matrix = SymFermionDetMatrix(
    fermion_path_integral,
    maxiter = maxiter, tol = tol
)

# Initialize pseudofermion field calculator.
pff_calculator = PFFCalculator(electron_phonon_parameters, fermion_det_matrix)

# Initialize KPM preconditioner.
preconditioner = KPMPreconditioner(fermion_det_matrix, rng = rng)

# Initialize Green's function estimator for making measurements.
greens_estimator = GreensEstimator(fermion_det_matrix, model_geometry)

# Initialize Hamiltonian/Hybrid monte carlo (HMC) updater.
hmc_updater = EFAPFFHMCUpdater(
    electron_phonon_parameters = electron_phonon_parameters,
    Nt = Nt, Δt = π/(2*Nt)
)
```

--------------------------------

### Initialize Tight-Binding and Electron-Phonon Parameters

Source: https://smoqysuite.github.io/SmoQyElPhQMC.jl/stable/tutorials/holstein_honeycomb_mpi

Initializes tight-binding and electron-phonon parameters using the model geometry and random number generator. No changes are needed from the previous 1a) Honeycomb Holstein Model tutorial.

```julia
    # Initialize tight-binding parameters.
    tight_binding_parameters = TightBindingParameters(
        tight_binding_model = tight_binding_model,
        model_geometry = model_geometry,
        rng = rng
    )

    # Initialize electron-phonon parameters.
    electron_phonon_parameters = ElectronPhononParameters(
        β = β, Δτ = Δτ,
        electron_phonon_model = electron_phonon_model,
        tight_binding_parameters = tight_binding_parameters,
        model_geometry = model_geometry,
        rng = rng
    )
```

--------------------------------

### Initialize simulation metadata

Source: https://smoqysuite.github.io/SmoQyElPhQMC.jl/stable/tutorials/holstein_honeycomb_mpi

Sets up the random number generator and a metadata dictionary recording simulation parameters. This section is unchanged from the previous 1a) Honeycomb Holstein Model tutorial.

```julia
    # Initialize random number generator
    rng = Xoshiro(seed)

    # Initialize metadata dictionary
    metadata = Dict()

    # Record simulation parameters.
    metadata["N_therm"] = N_therm  # Number of thermalization updates
    metadata["N_measurements"] = N_measurements  # Total number of measurements and measurement updates
    metadata["N_bins"] = N_bins # Number of times bin-averaged measurements are written to file
    metadata["maxiter"] = maxiter # Maximum number of conjugate gradient iterations
    metadata["tol"] = tol # Tolerance used for conjugate gradient solves
    metadata["Nt"] = Nt # Number of time-steps in HMC update
    metadata["Nrv"] = Nrv # Number of random vectors used to estimate fermionic correlation functions
    metadata["seed"] = seed  # Random seed used to initialize random number generator in simulation
    metadata["hmc_acceptance_rate"] = 0.0 # HMC acceptance rate
    metadata["reflection_acceptance_rate"] = 0.0 # Reflection update acceptance rate
    metadata["swap_acceptance_rate"] = 0.0 # Swap update acceptance rate
    metadata["hmc_iters"] = 0.0 # Avg number of CG iterations per solve in HMC update.
    metadata["reflection_iters"] = 0.0 # Avg number of CG iterations per solve in reflection update.
    metadata["swap_iters"] = 0.0 # Avg number of CG iterations per solve in swap update.
    metadata["measurement_iters"] = 0.0 # Avg number of CG iterations per solve while making measurements.
```

--------------------------------

### Initialize measurement container and basic measurements

Source: https://smoqysuite.github.io/SmoQyElPhQMC.jl/stable/examples/ossh_honeycomb

Initializes the measurement container and the basic measurements for the tight-binding and electron-phonon models.

```julia
# Initialize the container that measurements will be accumulated into.
measurement_container = initialize_measurement_container(model_geometry, β, Δτ)

# Initialize the tight-binding model related measurements, like the hopping energy.
initialize_measurements!(measurement_container, tight_binding_model)

# Initialize the electron-phonon interaction related measurements.
initialize_measurements!(measurement_container, electron_phonon_model)
```

--------------------------------

### Configure electron-phonon model and phonon modes

Source: https://smoqysuite.github.io/SmoQyElPhQMC.jl/stable/examples/bssh_square

Initializes the electron-phonon model and adds dispersionless and frozen phonon modes.

```julia
# Initialize a null electron-phonon model.
electron_phonon_model = ElectronPhononModel(
    model_geometry = model_geometry,
    tight_binding_model = tight_binding_model
)

# Define a dispersionless phonon mode to represent vibrations in the x-direction.
phonon_x = PhononMode(
    basis_vec = [0.0,0.0],
    Ω_mean = Ω
)

# Add x-direction optical ssh phonon to electron-phonon model.
phonon_x_id = add_phonon_mode!(
    electron_phonon_model = electron_phonon_model,
    phonon_mode = phonon_x
)

# Define a dispersionless phonon mode to represent vibrations in the y-direction.
phonon_y = PhononMode(
    basis_vec = [0.0,0.0],
    Ω_mean = Ω
)

# Add y-direction optical ssh phonon to electron-phonon model.
phonon_y_id = add_phonon_mode!(
    electron_phonon_model = electron_phonon_model,
    phonon_mode = phonon_y
)

# Define frozen phonon mode with infinite mass.
fphonon = PhononMode(
    basis_vec = [0.0,0.0],
    Ω_mean = Ω,
    M = Inf # Set phonon mass to infinity.
)

# Add frozen phonon mode to model.
fphonon_id = add_phonon_mode!(
    electron_phonon_model = electron_phonon_model,
    phonon_mode = fphonon
)
```

--------------------------------

### Initialize model parameters

Source: https://smoqysuite.github.io/SmoQyElPhQMC.jl/stable/examples/ossh_honeycomb

Initializes the tight-binding and electron-phonon parameters required for the simulation.

```julia
# Initialize tight-binding parameters.
tight_binding_parameters = TightBindingParameters(
    tight_binding_model = tight_binding_model,
    model_geometry = model_geometry,
    rng = rng
)

# Initialize electron-phonon parameters.
electron_phonon_parameters = ElectronPhononParameters(
    β = β, Δτ = Δτ,
    electron_phonon_model = electron_phonon_model,
    tight_binding_parameters = tight_binding_parameters,
    model_geometry = model_geometry,
    rng = rng
)
```

--------------------------------

### Execute simulation with mpiexecjl

Source: https://smoqysuite.github.io/SmoQyElPhQMC.jl/stable/tutorials/holstein_honeycomb_density_tuning

Run the simulation using mpiexecjl with 16 processes, specifying target density 0.8 and initial chemical potential 0.0.

```bash
mpiexecjl -n 16 julia holstein_honeycomb_density_tuning.jl 1 1.0 1.5 0.8 0.0 3 4.0 5000 10000 100 0.5
```

--------------------------------

### Initialize and Execute Simulation

Source: https://smoqysuite.github.io/SmoQyElPhQMC.jl/stable/examples/bssh_chain

Initializes MPI and executes the simulation with parameters parsed from command line arguments.

```julia
# Only execute if the script is run directly from the command line.
if abspath(PROGRAM_FILE) == @__FILE__

    # Initialize MPI
    MPI.Init()

    # Initialize the MPI communicator.
    comm = MPI.COMM_WORLD

    # Run the simulation.
    run_simulation(
        comm;
        sID = parse(Int, ARGS[1]), # Simulation ID.
        Ω = parse(Float64, ARGS[2]), # Phonon energy.
        α = parse(Float64, ARGS[3]), # Electron-phonon coupling.
        μ = parse(Float64, ARGS[4]), # Chemical potential.
        L = parse(Int, ARGS[5]), # System size.
        β = parse(Float64, ARGS[6]), # Inverse temperature.
        N_therm = parse(Int, ARGS[7]), # Number of thermalization updates.
        N_measurements = parse(Int, ARGS[8]), # Total number of measurements and measurement updates.
        N_bins = parse(Int, ARGS[9]), # Number of times bin-averaged measurements are written to file.
        checkpoint_freq = parse(Float64, ARGS[10]), # Frequency with which checkpoint files are written in hours.
    )

    # Finalize MPI.
    MPI.Finalize()
end
```

--------------------------------

### Initialize Standard Measurements

Source: https://smoqysuite.github.io/SmoQyElPhQMC.jl/stable/tutorials/holstein_honeycomb

Sets up the measurement container and initializes various correlation measurements including Green's functions, density, pair, and spin-z correlations.

```julia
    # Initialize the container that measurements will be accumulated into.
    measurement_container = initialize_measurement_container(model_geometry, β, Δτ)

    # Initialize the tight-binding model related measurements, like the hopping energy.
    initialize_measurements!(measurement_container, tight_binding_model)

    # Initialize the electron-phonon interaction related measurements.
    initialize_measurements!(measurement_container, electron_phonon_model)

    # Initialize the single-particle electron Green's function measurement.
    initialize_correlation_measurements!(
        measurement_container = measurement_container,
        model_geometry = model_geometry,
        correlation = "greens",
        time_displaced = true,
        pairs = [
            # Measure green's functions for all pairs or orbitals.
            (1, 1), (2, 2), (1, 2)
        ]
    )

    # Initialize the single-particle electron Green's function measurement.
    initialize_correlation_measurements!(
        measurement_container = measurement_container,
        model_geometry = model_geometry,
        correlation = "phonon_greens",
        time_displaced = true,
        pairs = [
            # Measure green's functions for all pairs or orbitals.
            (1, 1), (2, 2), (1, 2)
        ]
    )

    # Initialize density correlation function measurement.
    initialize_correlation_measurements!(
        measurement_container = measurement_container,
        model_geometry = model_geometry,
        correlation = "density",
        time_displaced = false,
        integrated = true,
        pairs = [
            (1, 1), (2, 2),
        ]
    )

    # Initialize the pair correlation function measurement.
    initialize_correlation_measurements!(
        measurement_container = measurement_container,
        model_geometry = model_geometry,
        correlation = "pair",
        time_displaced = false,
        integrated = true,
        pairs = [
            # Measure local s-wave pair susceptibility associated with
            # each orbital in the unit cell.
            (1, 1), (2, 2)
        ]
    )

    # Initialize the spin-z correlation function measurement.
    initialize_correlation_measurements!(
        measurement_container = measurement_container,
        model_geometry = model_geometry,
        correlation = "spin_z",
        time_displaced = false,
        integrated = true,
        pairs = [
            (1, 1), (2, 2)
        ]
    )
```

--------------------------------

### Initialize simulation with MPI process ID and data folder

Source: https://smoqysuite.github.io/SmoQyElPhQMC.jl/stable/tutorials/holstein_honeycomb_mpi

Constructs the data folder name, retrieves the MPI process ID using MPI.Comm_rank, initializes SimulationInfo with the process ID, and calls initialize_datafolder with the communicator to synchronize processes.

```julia
    # Construct the foldername the data will be written to.
    datafolder_prefix = @sprintf "holstein_honeycomb_w%.2f_a%.2f_mu%.2f_L%d_b%.2f" Ω α μ L β

    # Get MPI process ID.
    pID = MPI.Comm_rank(comm)

    # Initialize simulation info.
    simulation_info = SimulationInfo(
        filepath = filepath,
        datafolder_prefix = datafolder_prefix,
        write_bins_concurrent = (L > 7),
        sID = sID,
        pID = pID
    )

    # Initialize the directory the data will be written to.
    initialize_datafolder(comm, simulation_info)
```