The preparation of an equilibrium thermal state of a quantum many-body system on noisy intermediate-scale quantum (NISQ) devices is an important task in order to extend the range of applications of quantum computation. Faithful Gibbs state preparation would pave the way to investigate protocols such as thermalization and out-of-equilibrium thermodynamics, as well as providing useful resources for quantum algorithms, where sampling from Gibbs states constitutes a key subroutine. We propose a variational quantum algorithm (VQA) to prepare Gibbs states of a quantum many-body system. The novelty of our VQA consists in implementing a parameterized quantum circuit acting on two distinct, yet connected, quantum registers. The VQA evaluates the Helmholtz free energy, where the von Neumann entropy is obtained via postprocessing of computational basis measurements on one register, while the Gibbs state is prepared on the other register, via a unitary rotation in the energy basis. Finally, we benchmark our VQA by preparing Gibbs states of the transverse field Ising model and achieve remarkably high fidelities across a broad range of temperatures in statevector simulations. We also assess the performance of the VQA on IBM quantum computers, showcasing its feasibility on current NISQ devices.