TLDR: When you calculate phonons, you can describe electrons at different levels of theory, typically semilocal DFT, but also hybrids or dynamical mean-field theory. Phonons do include zero-point motion, as they are essentially a set of uncoupled quantum harmonic oscillators. Enthalpy can be calculated without reference to phonons, simply adding a PV term to the Hamiltonian. Gibbs free energy is calculated by adding a PV term to the phonon calculation.
Longer answer: The starting point of any phonon calculation is the Born-Oppenheimer approximation, which allows you to separate the electron and nuclear degrees of freedom. After applying this Born-Oppenheimer approximation, you end up with two eigenvalue equations. The first one corresponds to the electrons, in which the nuclei occupy fixed positions and their coordinates only appear as parameters. This electronic eigenvalue equation is what DFT codes solve.
Your question refers to the second eigenvalue equation that results from the Born-Oppenheimer approximation, which is the nuclear equation. The Hamiltonian in this equation reads (in atomic units):
$$
\hat{H}=-\sum_i\frac{1}{2m_i}\nabla_i^2+V(\mathbf{R}),
$$
in which the sum is the kinetic energy of the nuclei, and it runs over all nuclei $i$ in the system, and the second term is the potential energy felt by the nuclei, in which $\mathbf{R}=(\mathbf{r}_1,\mathbf{r}_2,\ldots)$ is a collective variable containing all individual nuclear coordinates $\{\mathbf{r}_i\}$. The value of this potential energy, typically called the potential energy surface, at a given collective nuclear coordinate $\mathbf{R}$ is given by the electronic eigenvalue for the nuclei fixed at this coordinate. This means that, unlike the electronic equation solved in DFT for which you know the Hamiltonian, in the case of the equation for the nuclei, you don't even know the Hamiltonian, as you don't know what $V(\mathbf{R})$ is. You first need to figure out what $V(\mathbf{R})$ is, and to do this you need to solve the electronic equation many times, once at each potential value of $\mathbf{R}$. This is clearly unfeasible, as $\mathbf{R}$ spans a 3N dimensional space, where N is the number of atoms in your system.
This is where the harmonic approximation that you mention comes in. For a material, we assume that the nuclei don't move much from their equilibrium positions. This is because they are relatively heavy (compared to the electrons), so rather than exploring the entire potential $V(\mathbf{R})$, they only explore this potential in the region near their equilibrium position, which corresponds to a minimum of $V(\mathbf{R})$. To proceed, I will make a change of coordinates $\mathbf{u}_i=\mathbf{r}_i-\mathbf{r}_i^0$, to collective coordinates $\mathbf{U}$ which are relative coordinates with respect to the equilibrium coordinates $\mathbf{R}^0$. In this way, equilibrium corresponds to $V(\mathbf{U}=0)$. In the harmonic approximation, we approximate this potential by a second-order Taylor expansion about equilibrium:
$$
V(\mathbf{U})\simeq V(0)+\sum_{\alpha,\beta}\frac{\partial^2V}{\partial u_{\alpha}\partial u_{\beta} }u_{\alpha}u_{\beta}.
$$
In the sum, $\alpha$ and $\beta$ are collective indices capturing the degrees of freedom of the 3N-dimensional energy surface (cell in the crystal, atoms in the basis, and Cartesian direction). To proceed, you replace this second-order approximation to the potential into the nuclear Hamiltonian I wrote above, and you can diagonalize it in terms of phonons. This second-order approximation works very well because, in terms of phonons, it essentially allows you to replace a 3N-dimensional potential with 3N 1-dimensional potentials, and the latter are much easier to calculate.
Therefore, to go back to your question: all that goes into a phonon calculation are the second order derivatives of the potential energy surface $V(\mathbf{R})$. What information does this contain? As $V(\mathbf{R})$ corresponds to the electronic energy at $\mathbf{R}$, then the level at which you treat the electrons affects $V(\mathbf{R})$. Typically, this is calculated using semi-local DFT (LDA, GGA), but there are a few studies that calculate $V(\mathbf{R})$ using hybrid DFT, which is much more expensive but more accurate, or even other beyond-DFT methods like dynamical mean-field theory.
Once you solve the nuclear problem in terms of phonons, you can calculate the vibrational contribution at finite temperature to the Helmholtz free energy (TS term). In these calculations, the phonons do have a zero-point contribution to the energy, as they are described by a set of uncoupled quantum harmonic oscillators.
For your other specific questions:
Enthalpy: for that you simply need to add a PV term to the Hamiltonian, so no need for phonons there. The calculations aren't any more expensive than a typical DFT calculation as PV is very cheap to evaluate.
Gibbs free energy: for this you need to add both PV and TS terms, you get the TS term from the Helmholtz free energy, and the PV term from the enthalpy.
