We present microscopic Molecular Dynamics simulations including the efficient Ewald sum procedure to study warm and dense stellar plasmas consisting of finite-size ion charges immerse in a relativistic neutralizing electron gas. For densities typical of Supernova matter and crust in a proto-neutron star, we select a representative single ion composition and obtain the virialized equation of state (vEoS). We scrutinize the finite-size and screening corrections to the Coulomb potential appearing in the virial coefficients B2, B3 and B4 as a function of temperature. In addition, we study the thermal heat capacity at constant volume, CV, and the generalized Mayer’s relation i.e. the difference CP − CV with CP being the heat capacity at constant pressure, obtaining clear features signaling the onset of the liquid-gas phase transition. Our findings show that microscopic simulations reproduce the discontinuity in CV, whose value lies between that of idealized gas and crystallized configurations. We study the pressure isotherms marking the boundary of the metastable region before the gaseous transition takes place. The resulting vEoS displays a behaviour where effective virial coefficients include extra density dependence showing a generalized density-temperature form. As an application we parametrize pressure as a function of density and temperature under the form of an artificial neural network showing the potential of machine learning for future regression analysis in more refined multicomponent approaches. This is of interest to size the importance of these corrections in the liquid-gas phase transition in warm and dense plasma phases contributing to the cooling behaviour of early Supernova phases and proto-neutron stars.