# How to write a naive Dirac matrix for Lattice QCD?

I'm trying to write down the naive Dirac matrix (with fermion doubling) for a LQCD simulation with one quark, for now. I initialized the $$SU(3)$$ gauge field and the quark field. The quak field has 4 space-time indices, 1 index for spinor components (from 1 to 4) and 1 index for the color (from 1 to 3), if you agree. My question is: how to couple together all these indices in Dirac matrix? And how many indices does the Dirac matrix have? I'm thinking of something like this (in Python):

spinor=4
mu=4 #(all 4 directions in space-time)
color=3
su3=3

quark=np.zeros((Nx, Ny, Nz, Nt, spinor, color), complex)
gaugeSU3=np.zeros((Nx, Ny, Nz, Nt, mu, su3, su3), complex)


Once filled all quark and gauge fields with appropriate values (for the gauge there are 4 SU(3) matrices for each space-time point), how to treat all the indices in Dirac matrix? Simulations of SU(3) pure gauge theory I performed are coherent with the literature, so I suppose it works, but how to couple the field with quarks? I add here the code I wrote for naive fermions; It seems to be incorrect because the Dirac matrix isn't a sparse matrix but all elements are filled with values different from zero:

 def DiracMatrix(U, psi, D):
#here U is the gaugeSU3, psi is the quark field, D the Dirac
# matrix, Dirac is the spinor indices
m = 0.2

for x in range(Nx):
for y in range(Ny):
for z in range(Nz):
for t in range(Nt):
for alpha in range(Dirac):
for beta in range(Dirac):
for a in range(color):
for b in range(color):
for mu in range(4):
a_mu = [0, 0, 0, 0]
a_mu[mu] = 1
D[x, y, z, t, alpha, beta, a, b] += 0.5 * (
gamma[mu][alpha, beta]
* U[x, y, z, t, mu, a, b]
* psi[
(x + a_mu[0]) % Nx,
(y + a_mu[1]) % Ny,
(z + a_mu[2]) % Nz,
(t + a_mu[3]) % Nt,
alpha,
a,
]
- U[
(x - a_mu[0]) % Nx,
(y - a_mu[1]) % Ny,
(z - a_mu[2]) % Nz,
(t - a_mu[3]) % Nt,
mu,
a,
b,
]
.conj()
.T
* psi[
(x - a_mu[0]) % Nx,
(y - a_mu[1]) % Ny,
(z - a_mu[2]) % Nz,
(t - a_mu[3]) % Nt,
beta,
b,
]
+ m
)

return D


On a 4d cubic lattice, the naive Dirac action just has the usual 4-by-4 $$\gamma_0$$ matrix on the time-direction links and the 4-by-4 $$\gamma_1$$ on the $$x$$-direction links, etc. The gamma's talk only to the spinor indices.
You need separate $${\rm SU}(3)$$ colour matrices $$U_\mu$$ for the gauge coupling on the links in the $$\mu$$ direction. Do you not have a book that explains all this?