The numerical representation of high-dimensional Gibbs distributions is challenging due to the curse of dimensionality manifesting through the intractable normalization constant calculations. This work addresses this challenge by performing a particle-based high-dimensional parametric density estimation subroutine, and the input to the subroutine is Gibbs samples generated by leveraging advanced sampling techniques. Specifically, to generate Gibbs samples, we employ ensemble-based annealed importance sampling, a population-based approach for sampling multimodal distributions. These samples are then processed using functional hierarchical tensor sketching, a tensor-network-based density estimation method for high-dimensional distributions, to obtain the numerical representation of the Gibbs distribution. We successfully apply the proposed approach to complex Ginzburg-Landau models with hundreds of variables. In particular, we show that the approach proposed is successful at addressing the metastability issue under difficult numerical cases.