  Accurate phase diagram calculation from molecular dynamics requires
systematic treatment and convergence of statistical averages. In this work we
propose a Gaussian process regression based framework for reconstructing the
free energy functions using data of various origin. Our framework allows for
propagating statistical uncertainty from finite molecular dynamics trajectories
to the phase diagram and automatically performing convergence with respect to
simulation parameters. Furthermore, our approach provides a way for automatic
optimal sampling in the simulation parameter space based on Bayesian
optimization approach. We validate our methodology by constructing phase
diagrams of two model systems, the Lennard-Jones and soft-core potential, and
compare the results with existing works studies and our coexistence
simulations. Finally, we construct the phase diagram of lithium at temperatures
above 300 K and pressures below 30 GPa from a machine-learning potential
trained on ab initio data. Our approach performs well when compared to
coexistence simulations and experimental results.
