We calculate per Monte Carlo evaluation on an 83 × 3 lattice the energy density e G of the gluon sector of QCD, including virtual quark loops up to the fourth power in the hopping parameter expansion. For light quarks of one flavour, we observe at T/A L ~ 95 -+ 10 a rapid variation of e G in T, accompanied by strong fluctuations from iteration to iteration, as clear signal of the deconfinement transition.The deconfinement transition, at which strongly interacting matter becomes a colour-conducting plasma, has so far been studied for pure Yang-Mills systems and for full QCD without virtual quark loops ("quenched approximation"). The concentration on Yang-Mills fields in the first studies of the phenomenon is physically meaningful: non-abelian gauge fields alone already exhibit confinement at low [1,2] and deconfinement at high [3] temperatures, at least on the lattice. The statistical mechanics of SU(N) gauge fields ,1 has therefore contributed significantly to our understanding of colour force thermodynamics. On the other hand, the neglect of virtual quark loops in the extension of full QCD is due to technical difficulties rather than to physical reasoning: the evaluation of the r X r fermion matrix, where r is some multiple of the number of lattice sites (typically 10 000 or more), quickly led to the limits of computer possibilities both in speed and in memory space. Nevertheless, it is the virtual quark loops which provide the "breaking of the string" through meson formation, so their role in deconfinement is certainly crucial. The advent of array processors has now put the inclusion of fermions within reach, and the aim of this paper is to present first results on the statistical mechanics of 1 Alexander yon Humboldt fellow, on leave from Hacettepe University, Ankara, Turkey. ,1 Fora survey see ref. A phase transition can be studied in two ways. We may consider some thermodynamic quantity -energy density, specific heat -and look for a discontinuity or singularity; or we may construct a specific order parameter to distinguish the two phases. Thus, in the Ising model, the Curie point can be characterized either by the singularity of the specific heat or through the vanishing of the spontaneous magnetization.In SU(N) Yang-Mills theory, deconfinement is related to the breaking of a global Z N symmetry [6,7]. The euclidean formulation of the partition function