  Accurate phase diagrams of multicomponent plasmas are required for the
modeling of dense stellar plasmas, such as those found in the cores of white
dwarf stars and the crusts of neutron stars. Those phase diagrams have been
computed using a variety of standard techniques, which suffer from physical and
computational limitations. Here, we present an efficient and accurate method
that overcomes the drawbacks of previously used approaches. In particular,
finite-size effects are avoided as each phase is calculated separately; the
plasma electrons and volume changes are explicitly taken into account; and
arbitrary analytic fits to simulation data are avoided. Furthermore, no
simulations at uninteresting state conditions, i.e., away from the phase
coexistence curves, are required, which improves the efficiency of the
technique. The method consists of an adaptation of the so-called Gibbs-Duhem
integration approach to electron-ion plasmas, where the coexistence curve is
determined by direct numerical integration of its underlying Clapeyron
equation. The thermodynamics properties of the coexisting phases are evaluated
separately using Monte Carlo simulations in the isobaric semi-grand canonical
ensemble. We describe this Monte Carlo-based Clapeyron integration method,
including its basic principles, our extension to electron-ion plasmas, and our
numerical implementation. We illustrate its applicability and benefits with the
calculation of the melting curve of dense C/O plasmas under conditions relevant
for white dwarf cores and provide analytic fits to implement this new melting
curve in white dwarf models. While this work focuses on the liquid-solid phase
boundary of dense two-component plasmas, a wider range of physical systems and
phase boundaries are within the scope of the Clapeyron integration method,
which had until now only been applied to simple model systems of neutral
particles.
