Cubic bicontinuous phases like the double gyroid (G), double diamond (D), and plumber's nightmare (P) are of great practical interest for many emerging applications requiring highly regular nanoscale networks or porous materials. Such phases can be formed from A-B diblock copolymers by the addition of A-type homopolymer over a range of compositions and relative chain lengths. Particle-based molecular simulations were used to delineate the phase diagram in a region where self-consistent field theory predicts the presence of a G D P triple point. Since the simulation box size must be commensurate with the morphology-specific 3D unit cell size (which is not known a priori), accurate free energy estimates are required for a range of box sizes, particularly when multiple competing phases can occur at the conditions of interest. A variant of thermodynamic integration was implemented to obtain such free energies (and hence identify the stable phases and their optimal box sizes) by tracing a reversible path connecting the ordered and disordered phases. This method overcomes key limitations of free energy methods based on the evaluation of chemical potentials via molecular insertions. For the range of conditions simulated, evidence was found of D-G, D-P, and G-D-P phase coexistence, consistent with previous theoretical predictions. Our simulations also reveal key differences and similarities in the size and microstructure of the nodes and struts that make up the different types of bicontinuous networks.